Articles | Volume 20, issue 19
Research article
06 Oct 2023
Research article |  | 06 Oct 2023

A modeling approach to investigate drivers, variability and uncertainties in O2 fluxes and O2 : CO2 exchange ratios in a temperate forest

Yuan Yan, Anne Klosterhalfen, Fernando Moyano, Matthias Cuntz, Andrew C. Manning, and Alexander Knohl

The O2 : CO2 exchange ratio (ER) between terrestrial ecosystems and the atmosphere is a key parameter for partitioning global ocean and land carbon fluxes. The long-term terrestrial ER is considered to be close to 1.10 mol of O2 consumed per mole of CO2 produced. Due to the technical challenge in measuring directly the ER of entire terrestrial ecosystems (EReco), little is known about variations in ER at hourly and seasonal scales, as well as how different components contribute to EReco. In this modeling study, we explored the variability in and drivers of EReco and evaluated the hypothetical uncertainty in determining ecosystem O2 fluxes based on current instrument precision. We adapted the one-dimensional, multilayer atmosphere–biosphere gas exchange model “CANVEG” to simulate hourly EReco from modeled O2 and CO2 fluxes in a temperate beech forest in Germany.

We found that the modeled annual mean EReco ranged from 1.06 to 1.12 mol mol−1 within the 5-year study period. Hourly EReco showed strong variations over diel and seasonal cycles and within the vertical canopy profile. The determination of ER from O2 and CO2 mole fractions in air above and within the canopy (ERconc) varied between 1.115 and 1.15 mol mol−1. CANVEG simulations also indicated that ecosystem O2 fluxes could be derived with the flux-gradient method using measured vertical gradients in scalar properties, as well as fluxes of CO2, sensible heat and latent energy derived from eddy covariance measurements. Owing to measurement uncertainties, however, the uncertainty in estimated O2 fluxes derived with the flux-gradient approach could be as high as 15 µmol m−2 s−1, which represented the 90 % quantile of the uncertainty in hourly data with a high-accuracy instrument. We also demonstrated that O2 fluxes can be used to partition net CO2 exchange fluxes into their component fluxes of photosynthesis and respiration if EReco is known. The uncertainty of the partitioned gross assimilation ranged from 1.43 to 4.88 µmol m−2 s−1 assuming a measurement uncertainty of 0.1 or 2.5 µmol m−2 s−1 for net ecosystem CO2 exchange and from 0.1 to 15 µmol m−2 s−1 for net ecosystem O2 exchange, respectively. Our analysis suggests that O2 measurements at ecosystem scale have the potential to partition net CO2 fluxes into their component fluxes, but further improvement in instrument precision is needed.

1 Introduction

Fluxes of O2 and CO2 between the terrestrial biosphere and atmosphere are inversely linked in photosynthesis, which assimilates CO2 and releases O2, and in respiration, which consumes O2 and releases CO2 (Keeling and Manning, 2014; Keeling and Shertz, 1992a; Krogh, 1919; Severinghaus, 1995). The relationship between these opposing fluxes can be described with the so-called O2 : CO2 exchange ratio (ER; see Table S1 in the Supplement for an overview of all abbreviations and variable names used here), which should be considered on various temporal and spatial scales – ranging from hourly to decadal scales temporally and from leaf to global scales spatially. Since the relationship of O2 and CO2 fluxes between the atmosphere and different carbon reservoirs (terrestrial biosphere, oceans and fossil fuels) differs on regional and global scales, these different ERs can be applied as parameters in global models in conjunction with observations of atmospheric O2 and CO2 abundances to quantify the global sinks of CO2 into the ocean and the terrestrial biosphere (Battle et al., 2000; Ishidoya et al., 2012; Keeling and Manning, 2014; Keeling and Shertz, 1992b; Tohjima et al., 2019). The global ER for the terrestrial biosphere is commonly set to 1.10 mol of O2 consumed per mole of CO2 produced (or vice versa) (Severinghaus, 1995) by assuming that this value, derived from elemental abundance data, is a representative long-term average for all land biota (Keeling and Manning, 2014; Manning and Keeling, 2006). An ER of 1.05 mol mol−1 was determined by Randerson et al. (2006) based on observed chemical compositions of plant parts for the quantification of the global carbon sink. Measurements using the oxidative ratio of organic material provided a more recent terrestrial ER estimate of 1.04 ± 0.03 mol mol−1 (Worrall et al., 2013). Using an ER of 1.05 mol mol−1 instead of 1.10 mol mol−1 in carbon budget models will contribute 0.05 Pg C yr−1 more to the global land carbon sink and an equivalent amount less to the ocean carbon sink (Keeling and Manning, 2014), indicating that the ER needs to be well constrained when parameterized in global ocean and land carbon cycle models.

On an ecosystem scale, a mole-fraction-based and a flux-based O2 : CO2 ratio can be considered (Ishidoya et al., 2013; Seibt et al., 2004). The former is defined as the fluctuations in the mole fraction of O2 per mole fraction of CO2 in the atmosphere (ERconc). Thus, ERconc is usually derived from the slopes of linear regressions between observed atmospheric O2 and CO2 mole fractions (Battle et al., 2019; Ishidoya et al., 2013; Seibt et al., 2004). Battle et al. (2019) observed an average ERconc=1.08± 0.007 mol mol−1 in a mixed deciduous forest over a 6-year period with temporal variations on a 6 h basis ranging between 0.85 and 1.15 mol mol−1. Measurements of canopy air O2 and CO2 mole fractions at two different forest sites yielded ERconc estimations between 1.01 and 1.03 mol mol−1 averaged over 24 h periods and between 1.14 and 1.19 mol mol−1 during daytime only (Seibt et al., 2004). Ishidoya et al. (2013) obtained differing ERconc at two heights within a cool temperate deciduous forest, reflecting variations in ERconc with canopy height. Furthermore, they observed different ERconc during daytime (0.87 mol mol−1) and nighttime (1.03 mol mol−1) in summer, indicating a significant variation in ERconc over the diel period (Ishidoya et al., 2013). Faassen et al. (2023) found much higher ERconc over 24 h (2.05 ± 0.03 mol mol−1) than for daytime only (1.10 ± 0.12 mol mol−1) and nighttime only (1.22 ± 0.02 mol mol−1) due to variations in the boundary layer height during the measurement period.

The flux-based O2 : CO2 ratio is defined as the O2 flux per CO2 flux between an ecosystem and the atmosphere (EReco). Flux estimates can be described as the net turbulent exchange or the overall net exchange (turbulent plus storage flux), and we focused on the latter in this study. Very few studies have attempted to quantify EReco because measuring O2 fluxes at ecosystem scale is still a major challenge. Since O2 and CO2 are strongly anti-correlated in the processes of photosynthesis and respiration, changes in both scalars are very similar in absolute numbers, typically on the order of a few parts per million. However, the relative changes in O2 are much smaller than in CO2 owing to the much higher atmospheric abundance (around 210 000 ppm for O2 and around 400 ppm for CO2), making O2 measurements at sufficient precision and accuracy technically challenging. Thus, previous studies resorted to, for instance, the flux-gradient method, chamber measurements and modeling approaches. Ishidoya et al. (2015) determined a daily mean net turbulent ER = 0.86 mol mol−1 based on O2 and CO2 gradient measurements. Faassen et al. (2023) reported daytime and nighttime EReco as 0.92 ± 0.17 and 1.03 ± 0.05 mol mol−1, respectively. In general, EReco depends on the elemental composition and reduction state of organic material and on the temporal variation and spatial distribution of sinks and sources of ecosystem flux components (Seibt et al., 2004). As described by Battle et al. (2019), the dynamics and interrelations of the various sinks and sources within the ecosystem, each with their own EReco, result in the mixed signal ERconc.

Current micrometeorological approaches to measure gas exchange between ecosystems and the atmosphere include eddy covariance, flux-gradient and eddy accumulation methods, which could all theoretically be used to determine ecosystem O2 fluxes. The applicability of the eddy covariance technique for O2 flux estimation, however, requires high precision at a high measurement frequency (10–20 Hz). Except for a homemade, non-commercial vacuum ultraviolet (VUV) absorption analyzer (Stephens et al., 2003), no suitable instrument exists so far. The application of the eddy accumulation method is also technically challenging and has not yet been applied to O2 (Emad and Siebicke, 2023a, b).

With the flux-gradient method, O2 fluxes can be inferred from an O2 gradient above a canopy and from an eddy diffusivity (K), which can be derived based on additional CO2, sensible heat or latent heat flux measurements (Baldocchi et al., 1988). This method assumes that heat and mass are transported in a similar manner between two adjacent levels above the canopy (Baldocchi et al., 1988). The method's applicability is again particularly challenging for O2 estimates owing to the typically large measurement uncertainty in relation to the small O2 gradient. One approach to increase the measurement-to-noise ratio is to move the lower inlet of the gradient measurement closer to or even inside the canopy. This approach, however, violates the assumption of the flux-gradient method owing to infrequent but predominantly large eddies within the canopy, counter-gradient fluxes and possible non-differentiable gradients (Raupach, 1989; Wilson, 1989). The flux-gradient method has already been used for O2 flux estimation above a cool temperate forest (Ishidoya et al., 2015), an urban canopy (Ishidoya et al., 2020) and a boreal forest (Faassen et al., 2023). The theoretical limits of the flux-gradient method for O2 fluxes given current instrument precision and accuracy, however, have not yet been fully explored.

Chamber-level gas exchange measurements provide an alternative approach to measure the ER of individual components such as leaf, stem and soil, which could be scaled up to ecosystem level. Branch and soil chamber measurements in a German temperate forest showed an average ER of leaf net assimilation (ERAn; net assimilation defined as carboxylation minus photorespiration and dark respiration) between 1.08 ± 0.16 and 1.19 ± 0.12 mol mol−1, as well as an ER of soil respiration (ERsoil) of 0.94 ± 0.04 mol mol−1 (Seibt et al., 2004). In a cool temperate deciduous forest in Japan, chamber measurements indicated an ERAn=1.02± 0.03 and ERsoil=1.11± 0.01 mol mol−1 (Ishidoya et al., 2013). Hilman et al. (2019) measured an average ER of stem respiration (ERstem) between 0.97 and 1.95 mol mol−1 for tropical, temperate and Mediterranean trees with a closed-flow chamber system with two continuous flow analyzers.

The ER variability in assimilation and respiration fluxes found in these studies provides a potential approach to partition net CO2 fluxes into their components following similar approaches based on stable isotopes in CO2 (Knohl and Buchmann, 2005; Ogee et al., 2004; Wehr and Saleska, 2015; Zobitz et al., 2007). Using simultaneous measurements of net ecosystem O2 and CO2 fluxes and considering the ER for the photosynthetic and respiratory processes in a canopy and at the soil surface, two mass balance equations can be written for O2 and CO2 (see Eq. 1 below). Hourly or half-hourly ER would need to agree with the typical time step of flux estimates derived with the eddy covariance technique, which is the standard method of measuring gas exchange between land surfaces and the atmosphere (Baldocchi et al., 2001; Goulden et al., 1996). Theoretically, such an O2-based partitioning method only works for periods when the ER of gross assimilation (ERA) and ER of gross ecosystem respiration (ERR) differ because a second independent mass balance equation is needed to yield CO2 fluxes of assimilation (FA) and respiration (FR). According to Ogee et al. (2004), the difference in ER has to be large enough to obtain a reasonable accuracy in the partitioned net CO2 fluxes. Consequently, an analysis of temporal dynamics in ERA and ERR is necessary in order to evaluate the possibility of applying O2 observations in a CO2 flux-partitioning approach.

The contribution of flux components to the temporal and spatial variability in overall ecosystem O2 fluxes can also be explored by modeling approaches. For example, net turbulent ER was simulated with a simple one-box model with daily time steps by assuming that O2 and CO2 mole fractions are spatially constant and temporally variable within the canopy (Ishidoya et al., 2015; Seibt et al., 2004). These simulations indicated that variations in net turbulent ER are influenced not only by leaf and soil fluxes but also by turbulence inside and outside the canopy (Seibt et al., 2004). To explore the drivers of ER variations at the ecosystem scale, more precise turbulence effects need to be considered. However, simple one-box models assume uniform and well-mixed air columns throughout the canopy, with the result that modeled ER lacks variations for different layers within the canopy.

Multilayer atmosphere–biosphere models such as CANVEG (Baldocchi, 1997; Baldocchi and Wilson, 2001) differ from one-box models in that they are designed to represent the temporal and (vertical) spatial scale of an eddy covariance tower. Therefore, they are a good simulator to test and examine new types of observations (Oikawa et al., 2017). CANVEG includes within-canopy transport of CO2, water vapor and energy (Baldocchi, 1997; Baldocchi and Wilson, 2001), with the result that if it were adapted to O2 processes, one could evaluate the accuracy of different flux measurement techniques such as eddy covariance or flux-gradient approaches. Published ER values of gross and/or net assimilation, stem respiration, and soil respiration can be employed as parameters to derive component-specific O2 fluxes from existing modeled CO2 fluxes. Thus, concurrent O2 and CO2 fluxes and ER can be plausibly simulated for multiple canopy layers and for the whole ecosystem, with which we can analyze the main drivers of modeled ER values, their diel and seasonal variability, and vertical variations. In addition, concurrently simulated mole fraction profiles – a function of turbulent dispersion and the strength and location of scalar sources and sinks – enable us to test the precision of the flux-gradient method for O2 flux estimation while choosing various measurement heights inside and above the canopy. Furthermore, the performance of an O2-based source-partitioning method can be evaluated based on model simulations.

Given these considerations, we defined the following objectives for this study: (1) to implement atmosphere–biosphere O2 : CO2 exchange ratios for various ecosystem components in the multilayer CANVEG model; (2) to explore temporal and spatial variations in O2 : CO2 exchange ratios at ecosystem scale, as well as the underlying main drivers at ecosystem scale; (3) to evaluate the potential precision of the flux-gradient approach to obtain O2 fluxes; and (4) to evaluate the feasibility of O2 flux measurements for CO2 flux partitioning.

2 Methods

2.1 Site description

The meteorological and plant-specific ecophysiological measurements used in our model simulation were derived from the Leinefelde FLUXNET tower site (DE-Lnf; located in central Germany (511942′′ N, 102204′′ E; 450 m a.s.l.; Anthoni et al., 2004). The vegetation at the site is an even-aged managed beech stand (Fagus sylvatica L.) with an age of approximately 130 years (Tamrakar et al., 2018). Between 2002 and 2016, the mean annual temperature was 8.3 ± 0.7 C and the average cumulative annual precipitation was 600 ± 150 mm (Braden-Behrens et al., 2019). The canopy height (ht) was 37.5 m, and the effective leaf area index (LAI) was at maximum 4.8 m2 m−2 in the growing season in 2015 (Braden-Behrens et al., 2017).

Meteorological variables are continuously measured including air temperature, air humidity, direct and diffuse global radiation, photosynthetic photon flux density, wind velocity, air pressure, vapor pressure deficit, precipitation, atmospheric CO2 mole fraction (CO2 atm), soil temperature, and soil moisture. Also, fluxes of net ecosystem CO2 exchange (FCO2), sensible heat (H) and latent heat (LE) are obtained with the eddy covariance technique at 44 m above the ground level (Anthoni et al., 2004). The meteorological variables were used as input data for our model simulations, while the flux estimates were storage-term-corrected and then used for model calibration and validation (see below). In this paper, upward fluxes (release to the atmosphere) are presented as positive quantities and downward fluxes (uptake by the ecosystem) as negative quantities. Thus, O2 fluxes always have opposite signs to their corresponding CO2 fluxes, which is in line with micrometeorological conventions.

2.2 Model description and model setup

We used the one-dimensional, multilayer atmosphere–biosphere gas exchange model “CANVEG”, described by Baldocchi (1997) and Baldocchi and Wilson (2001). The model domain included 120 model layers above the ground, in which the lower 40 aboveground layers covered the entire canopy, while the bottom layer represented the soil surface for the description of soil carbon and energy fluxes. The domain also included 10 belowground soil layers; however, this study did not consider processes within the soil column in any detail. CANVEG used hourly meteorological variables as drivers, as well as site-specific parameters (see Table 1), to simulate atmosphere–biosphere water vapor, CO2 and energy fluxes within and above the forest canopy.

Table 1Model parameters adjusted to the study site Leinefelde, Germany.

 Parameters were calibrated with eddy covariance measurements of hourly FCO2, FR, H and LE in 2012 and 2013 via the Markov chain Monte Carlo (MCMC) method.

Download Print Version | Download XLSX

The carbon, water and energy modules in CANVEG have been validated for various environmental conditions and forest types (Baldocchi, 1997; Baldocchi et al., 2002, 1999). Moreover, CANVEG has previously been applied to an unmanaged beech-dominated forest site only 30 km away from the site of this study (Knohl and Baldocchi, 2008) and has recently been used to simulate the isotopic composition of carbon assimilates at Leinefelde (Braden-Behrens et al., 2019). We translated the original C code (Baldocchi, 1997) to Fortran 90, which was then used for further implementations.

Atmospheric O2 mole fraction (O2 atm) as an input for the model was deduced from a fixed O2 : CO2 mole ratio of 1.15 mol mol−1 and continuous CO2 mole fraction measurements at the site (Table 1). The fixed O2 : CO2 mole ratio was derived from measurements at the University of Göttingen from November 2017 to January 2018 using a high-precision O2 measurement system developed by Penelope Pickers (University of East Anglia, UK), which is very similar to the system described in Pickers et al. (2017). For these measurements, the correlation between O2 and CO2 mole fractions had R2=0.99.

Some model parameters regarding leaf photosynthesis, stomatal conductance and soil respiration were fitted to the actual site conditions via the Markov chain Monte Carlo (MCMC) method (Van Oijen et al., 2005). Eddy covariance measurements of hourly FCO2, H and LE, as well as the estimated ecosystem respiration (FR) in 2012 and 2013, were used to calibrate the model parameters (Table 1). The years 2014–2016 were used for model validation. The leaf phenology parameters, including day of year (DOY) for the start of leaf growth, end of leaf growth, start of leaf fall and end of leaf fall (leafout, leaffull, leaffall and leaffall_complete), were derived from daily camera images in 2015 above the canopy. LAI during the course of a year was simulated based on these four parameters: the DOY range before leafout and after leaffall_complete was defined as winter when LAI = zero, and the DOY range between leaffull and leaffall was defined as summer when LAI = 4.8 m2 m−2. During spring (leafout< DOY  < leaffull) and during autumn (leaffall< DOY < leaffall_complete) LAI increased or decreased linearly, respectively. The maximum LAI of 4.8 m2 m−2 and the LAI fraction (fLAI) at five different heights in the canopy were measured using a LI-2000 plant canopy analyzer (LI-COR Biosciences GmbH, Germany) in 2015 (Braden-Behrens et al., 2017). The vertical LAI profile was assumed to follow a beta distribution, which was fitted to the observed fLAI (Table 1). This relationship between LAI and height (z) allocates leaves mainly in the upper canopy (z/ht 0.45) with almost no leaves in the bottom canopy (Fig. 1a). The wood area index (WAI) consisted of the branches (80 % of total WAI) and the stems (20 % of total WAI). The branches were situated in the upper canopy (z/ht 0.45) following the same distribution algorithm as LAI, while in the lower canopy (z/ht< 0.45), the fraction of stem WAI per layer to total stem area was deduced from the fraction of stem diameter per layer to the diameter at breast height (fDBH) as a function of height (z): fDBH= 102  2.6z+ 0.08z2 0.0023z3 (Grundner et al., 1952). This setup of the forest canopy including leaf phenology and the vertical LAI and WAI profiles was used for all years of the model run. All site-specific parameters used in this study are listed in Table 1.

Figure 1(a) Distribution of vertical leaf and wood area indices (LAI and WAI in m2 m−2 per canopy layer) used in the CANVEG model, derived from measurements at the Leinefelde study site (Braden-Behrens et al., 2017). The y axis is the ratio of the height in the canopy (z) to the top of the canopy (ht). (b) O2 : CO2 exchange ratio of leaf dark respiration (ERrd in mol mol−1) as a function of leaf temperature (Tleaf in C) after Tcherkez et al. (2003). The dashed red lines indicate the range of Tleaf and corresponding ERrd in this study.


For the simulation of net ecosystem O2 fluxes (FO2), values of ER had to be chosen: the input parameter of ERA was set to 1.00 mol mol−1 (Table 1) by assuming that photosynthesis produces glucose (C6H12O6), resulting in equal O2 and CO2 fluxes. The ER of canopy respiration was attributed to the ER of leaf dark respiration (ERrd) and stem respiration (ERstem). ERstem was fixed to 1.04 mol mol−1 (Randerson et al., 2006), while the ERrd was set to increase with leaf temperature (Tleaf; Fig. 1b) according to Tcherkez et al. (2003). ERsoil was set to 1.10 mol mol−1 (Randerson et al., 2006; Severinghaus, 1995). To quantify the dependency of the CANVEG model regarding these fixed ER parameters, we also conducted a sensitivity analysis, where we changed ERA, ERstem and ERsoil each by ± 10 % and estimated the resulting relative changes in simulated O2 fluxes. Furthermore, the impact of changed ER parameters was also investigated in the following parts of this study (see Sect. 2.3 and 2.5 below).

To validate the model, we used eddy covariance measurements of FCO2, H and LE from 2014 to 2016. To quantify the model performance, we calculated the slope, intercept and the coefficient of determination (R2) of a linear regression between modeled and observed FCO2, H and LE, as well as the root mean square error (RMSE).

2.3 Model simulations of flux- and mole-fraction-based exchange ratios

In CANVEG, CO2 fluxes are simulated for the leaf, stem and soil components. The O2 fluxes of each component are estimated by scaling each corresponding CO2 flux by its ER. Respiratory CO2 fluxes are defined to be positive, while assimilation CO2 fluxes are negative. O2 fluxes always have the opposite sign from the corresponding CO2 fluxes, which would result in negative ER values. However, we have defined all ER parameters to be positive by including the factor (1) in all relevant equations (see below) to be consistent with most published literature concerning O2 : CO2 exchange ratios (Ishidoya et al., 2013; Seibt et al., 2004). Another way of considering this is that the ERs are the ratios of moles of O2 consumed per mole of CO2 produced (or moles of O2 produced per mole of CO2 consumed).

The O2 and CO2 ecosystem fluxes are the balance of the simulated fluxes of gross assimilation (FA, carboxylation minus photorespiration) and gross ecosystem respiration (FR). The latter consists of leaf dark respiration (Frd), stem respiration (Fstem) and soil respiration (Fsoil, consisting of 50 % respiration by heterotrophs and 50 % by autotrophs):

(1) F CO 2 = F A + F rd + F stem + F soil = F A + F R F O 2 = - F A ER A - F rd ER rd - F stem ER stem - F soil ER soil = - F A ER A - F R ER R ,

where ERA, ERrd, ERstem and ERsoil are given as model parameters (see Sect. 2.2.). The simulated FO2 and FCO2 include the storage fluxes associated with changes in O2 and CO2 mole fractions in the canopy air space because they were inferred by integrating fluxes for all canopy layers. In general, the CANVEG model only considered dry mole fractions of O2 and CO2. Usually, O2 measurements are reported in per meg, which describes the change in the O2 to N2 ratio relative to a reference. To convert from parts per million to per meg, the factor 1/0.2095= 4.8 per meg ppm−1 can be used, where 0.2095 represents the O2 mole fraction of air (in mol mol−1). In this study, we chose mole fraction as the unit for O2 to be consistent with regard to the calculation of O2 : CO2 exchange ratios, which are usually presented in moles per mole (mol mol−1).

For the model simulations, ER could be obtained for the entire ecosystem, for the net assimilation at the leaf level or for only respiratory processes by considering the simulations of the corresponding flux components. The ER of the overall ecosystem (EReco) in hourly time steps was calculated as the ratio of the hourly FO2 and FCO2 (including storage terms) summed up over the entire canopy height; that is,

(2) ER eco = - F O 2 F CO 2 .

EReco for specific canopy heights (ERecoz) was derived as the slope of linear regressions fitted to O2 and CO2 fluxes of multiple simulated time steps for each canopy layer.

Furthermore, the simulated ERs of net O2 and CO2 assimilation (ERAn) and of all respiratory fluxes (ERR) were derived as


Moreover, we assessed the impact of the model parameters ERA, ERstem and ERsoil by changing each by ± 10 % on estimates for EReco and ERAn within the sensitivity analysis.

The atmospheric O2 mole fraction at each canopy layer was also computed by CANVEG, analogous to that done for CO2 mole fraction (Baldocchi, 1997). CANVEG estimated atmospheric mole fraction per layer as a function of multilayer gas flux diffusion determined by a Lagrangian dispersion matrix (Baldocchi, 1992) and the atmospheric background gas mole fraction. The mole-fraction-based ER (ERconc) and ERconc at specific canopy heights (ERconcz) were defined as the ratio between the fluctuations in O2 and CO2 mole fractions, and both were calculated as the slopes of linear regressions fitted to hourly atmospheric O2 versus CO2 mole fractions for the growing seasons (the days of year with leaves in the canopy, between leafout and leaffall_complete) of all simulation years (Battle et al., 2019; Ishidoya et al., 2013; Seibt et al., 2004). Thus, we obtained ERecoz and ERconcz with the same approach by deriving the slopes of hourly data to allow for a comparison.

2.4 Evaluation of the flux-gradient method to obtain O2 fluxes

The CANVEG simulations of ecosystem O2 fluxes and O2 mole fraction gradients provided the opportunity to test the applicability of the flux-gradient approach to estimate FO2. We assumed the flux-gradient measurement system could be installed both above the canopy and close to the forest floor. We especially aimed at testing the performance of the flux-gradient method based on current typical instrument performance for O2 measurements. The turbulent O2 (FO2), CO2 (FCO2), sensible heat (H) and latent heat (LE) fluxes are related to vertical scalar gradients as follows (Meredith et al., 2014):

(5) F O 2 = - K o Δ o Δ z ρ n F CO 2 = - K c Δ c Δ z ρ n H = - K T Δ T Δ z ρ m c p LE = - K v Δ v Δ z λ ,

where Δz (m) is the vertical height difference between the two measurement heights; ΔT, Δv, Δc and Δo denote the difference in air temperature (K), water vapor (kg m−3), CO2 dry-air mole fraction (ppm) and O2 dry-air mole fraction (ppm) between measurement heights, respectively; ρn and ρm are the molar density (mol m−3) and mass density of the air (kg m−3), respectively; cp is the specific heat capacity of air (J kg−1 K−1); and λ is the latent heat of evaporation (J kg−1). The superscript tilde in the flux nomenclatures denotes turbulent fluxes (without storage fluxes). Ko, Kc, KT and Kv (m2 s−1) are the eddy diffusivities of the relevant scalars. Assuming that heat and mass are transported in a similar way between two adjacent levels above the canopy and so assuming that Ko=Kc=KT=Kv (Baldocchi et al., 1988), then O2 fluxes can be estimated with each of the following equations:

(6) F O 2 , c = F CO 2 Δ o Δ c F O 2 , T = H Δ o ρ n Δ T ρ m c p F O 2 , v = LE Δ o ρ n Δ v λ .

From simulations of FCO2, H and LE and vertical scalar profiles, we derived FO2 from FO2 plus the storage term based on the flux-gradient method and compared these to the directly modeled FO2 (Eq. 1). Here, the subscripts c, T and v denote which flux and scalar are used (CO2 mole fraction, air temperature or water vapor, respectively).

There are usually three main sources of error in the flux-gradient method: (1) the uncertainty in the vertical gradient (that is, the gradient of O2 mole fraction, ΔoΔz) resulting from the precision and accuracy of the measurement instruments; (2) the magnitude of the mole fraction difference (Δc, ΔT or Δv) between the two measurement heights, which is usually small when the measurement heights are too close to each other or when the atmosphere is well mixed; and (3) the measurement uncertainty in the turbulent fluxes (FCO2, H or LE), which we assumed to be zero because we applied here only our simulated turbulent fluxes. So here, we quantified the extent of the first two sources of uncertainty and defined conditions when the flux-gradient method could perform satisfactorily to obtain FO2. The influence of the first uncertainty was evaluated by adding a “measurement error” to Δo, where the uncertainty was assumed to be normally distributed with a mean of zero and a standard deviation of ± 0.7 ppm (3.36 per meg) based on typical measurement uncertainty of the O2 mole fraction instrument used to derive the fixed atmospheric O2 : CO2 ratio (Pickers et al., 2017). Then the difference between the FO2 derived via the flux-gradient method with and without the measurement uncertainty (σFO2) was evaluated.

The second uncertainty due to the magnitude in the gradient as a function of Δz was analyzed by estimating FO2 based on the flux gradient between a top measurement height at 2 times the canopy height in our model setup and each layer below until the soil surface (z/ht=0). The top measurement height was set to z/ht=2 following customary recommendations for the setup of eddy covariance towers following Rebmann et al. (2018). We also included measurement heights inside the canopy, where the vertical profiles are mostly nonlinear due to scalar sources and sinks, to illustrate the effect of violating the assumptions of the flux-gradient method. For comparison, the difference between the FO2 estimations derived by the flux-gradient method (FO2,(c,T,v), based on FCO2, H or LE and their respective vertical scalar profile) and by model simulations (FO2,CANVEG) was calculated:

(7) Δ F O 2 , ( c , T , v ) = F O 2 , ( c , T , v ) - F O 2 , CANVEG ,

where ΔFO2 is the difference for the application between the top measurement height (z/ht=2) and each layer below.

Finally, we also tested a three-height flux-gradient method based on the recent study of Faassen et al. (2023). They derived scalar concentrations at three heights (z/ht=0.9, 3.7 and 6.9 with ht = 18 m), fitted a quadratic scalar–height relationship and expressed the vertical gradient as the first derivative of z (see Eqs. 10 and 11 in Faassen et al., 2023). In our study, we selected the three heights at z/ht=1.05, 1.45 and 2 with ht =37.5 m to be with all heights above the canopy.

2.5 Uncertainties in partitioning net ecosystem CO2 fluxes based on O2 fluxes

The net ecosystem CO2 exchange (FCO2) consists of two different components: gross assimilation (FA) and gross ecosystem respiration (FR). Similar to the stable-isotope flux-partitioning approach (Bowling et al., 2001; Knohl and Buchmann, 2005; Ogee et al., 2004; Oikawa et al., 2017; Yakir and Wang, 1996), O2 and CO2 flux mass balance equations can be written as shown in Eq. (1), where FCO2 is the observed ecosystem flux from eddy covariance measurements and FO2 is obtained by multiplying FCO2 by the modeled EReco in CANVEG following Eq. (2) (owing to the lack of actual FO2 measurements). We treated these mass balance equations as a probabilistic process assuming terms on the right-hand side are uncertainty quantities with a priori values (FAb,FRb,ERAb,ERRb) and uncertainties (σFAb,σFRb,σERAb,σERRb). Fluxes and exchange ratios, i.e., FA, FR, ERA and ERR, can then be calculated that minimize the differences between the left-hand side observations and the right-hand side “model” under consideration of their uncertainties, leading to a posteriori quantities (FA,FR,ERA,ERR) with corresponding uncertainties (σFA,σFR,σERA,σERR). A cost function (J) was then written as a linear system with all the differences weighted by the corresponding a priori uncertainties:

(8) J = 1 2 [ F A + F R - F CO 2 σ F CO 2 2 + F A ER A + F R ER R - F O 2 σ F O 2 2 + F A - F A b σ F A b 2 + F R - F R b σ F R b 2 + ER A - ER A b σ ER A b 2 + ER R - ER R b σ ER R b 2 ] .

The last four terms allow a solution to be defined with fewer equations than unknowns. The a posteriori values and uncertainties were returned at minimum J with predefined a priori values and uncertainties (Table 2). For the J function with multiple variables as in our case, the a posteriori means of any parameter, x, were found along the gradient of each variable where its Jacobian equaled zero (Jx=0; Tarantola, 2004), while the corresponding a posteriori uncertainties were expressed as the square root of the inverse Hessian at the minimum (2Jx2; Tarantola, 2004):

(9) σ F A σ F R σ ER A σ ER R = 2 J F A 2 2 J F A F R 2 J F A ER A 2 J F A ER R 2 J F R F A 2 J F R 2 2 J F R ER A 2 J F R ER R 2 J ER A F A 2 J ER A F R 2 J ER A 2 2 J ER A ER R 2 J ER R F A 2 J ER R F R 2 J ER R ER A 2 J ER R 2 - 1 .

By assuming no correlations among the variables, only the diagonal elements of the Hessian were used in a posteriori uncertainty calculations.

Table 2Assigned a priori values and uncertainties to build the cost function, J, for the uncertainty estimation of using O2 fluxes to partition net CO2 fluxes.

Download Print Version | Download XLSX

We evaluated the a posteriori uncertainties of partitioned photosynthetic fluxes on a typical day during summer (4 July 2012) with assigned a priori uncertainties. The a priori uncertainty of gross assimilation (σFAb) was set to 10 µmol m−2 s−1 and of ecosystem respiration (σFRb) to 5 µmol m−2 s−1, following Ogee et al. (2004) assuming less constraint on a posteriori results (Table 2). The uncertainty of the net CO2 fluxes (σFCO2) was derived from Mann and Lenschow's model (Lenschow et al., 1994) and calculated for our site to be 2.5 µmol m−2 s−1 (Braden-Behrens et al., 2019). We also examined if σFA could be reduced if more accurate net CO2 fluxes were measured (σFCO2=0.5µmol m−2 s−1).

The uncertainty of measured ecosystem O2 fluxes (σFO2) is unknown to us. Consequently, we used the results from the flux-gradient method evaluation (Sect. 2.4.). In order to clearly quantify the effect of σFO2 and σERA on flux-partitioning precision, we defined a σFO2 series ranging from 0.1 to 15 µmol m−2 s−1, representing the 90 % quantile of random Δo measurement uncertainty (see Sect. 2.4.), and a series of σERAb ranging from 0.001 to 0.1 mol mol−1. σERRb was fixed to either 0.05 or 0.001 mol mol−1.

Moreover, we assessed the impact of the model parameters ERA, ERstem and ERsoil on the source-partitioning results by changing each by ± 10 % and by estimating the absolute change in the a posteriori σFA.

3 Results

3.1 Model performance

The model generally showed similar performance for FCO2, H and LE during both calibration and validation (Fig. 2), indicating robust model behavior as a multilayer canopy flux simulator. The model validation for FCO2 (R2=0.82, slope of 1.016) was generally better than for H (R2=0.7, slope of 0.879) and LE (R2=0.77, slope of 1.02) (Fig. 2b, d and f). The disagreement between modeled and measured FCO2 indicated some uncertainties in the parameters for soil and stem respiration, as well as phenology, in the model equations. The similar scale but opposite sign of y intercepts for H and LE calibration simulations (Fig. 2c and e) indicated underestimation in H and the same amount of overestimation in LE. The slopes deviating from 1 for H and LE could come from a non-closure of the energy balance in the eddy covariance observations.

Figure 2Comparison of (a, b) net ecosystem CO2 flux (FCO2), (c, d) sensible heat flux (H), and (e, f) latent heat flux (LE) from 2012–2016 between model simulations (y axes) and eddy covariance observations (x axes). The left column shows all hourly data points for the calibration period (2012–2013), and the right column shows all hourly data points for the validation period (2014–2016). The linear regression line function, coefficient of determination (R2) and the root mean square error (RMSE) are included in each panel. The dashed lines are the 1:1 lines.


Due to potential variations in the ER model parameters (which were here taken from literature), we conducted a sensitivity analysis to show how these parameters affected the modeled FO2. If ERA was increased or decreased by 10 %, the modeled FO2 sum of the entire study period increased or decreased on average by 20.3 % correspondingly. Similarly, a change by plus or minus 10 % increments in ERsoil and ERstem caused the FO2 sum to decrease or increase by 8.6 % and 1.7 %, respectively. These results directly follow Eq. (1), where the derivative with respect to a specific ER gives the corresponding flux in percent. Oxygen fluxes were hence most sensitive to the ER of the largest carbon fluxes.

3.2 Temporal dynamics of O2 : CO2 exchange ratios

The median of the hourly ecosystem O2 : CO2 exchange ratio (EReco) throughout the simulation period (2012–2016) was 1.08 mol mol−1, where the annual medians did not differ between years. The annual mean EReco ranged from 1.06 to 1.12 mol mol−1 across the 5 years. Hourly EReco also varied seasonally and within the diel course, as shown as an example for the year 2012 in Fig. 3a. During the non-growing season, EReco values were constrained between 1.04 and 1.10 mol mol−1, representing a mixture of the prevailing stem and soil respiration processes. During the growing season, EReco was close to 1.00 mol mol−1 during daylight hours due to the dominance of photosynthetic processes and sometimes even smaller than 1.00 mol mol−1 when daytime FO2 was smaller than daytime FCO2. This could occur with ERA= 1.00 mol mol−1 and ERstem, ERsoil and ERrd >1.00 mol mol−1 (following Eq. 1) when more O2 was consumed than CO2 released for the respiratory fluxes, and thus the magnitude of net FO2 decreased. During nighttime in the growing season, EReco was > 1.00 mol mol−1, representing a mixture of stem, soil and leaf dark respiration. For transition periods (sunrise and sunset), with flux magnitudes close to zero, EReco values were very high, owing to very small FCO2. Because EReco is a ratio, values could get extremely large and approach infinity as FCO2 approaches zero. However, since corresponding FO2 values were also very low, these EReco values had very little effect on median and mean EReco of the overall ecosystem over a longer time period.

Within the sensitivity analysis, the initial annual median EReco of 1.08 mol mol−1 changed only by up to 0.02 mol mol−1 due to the change in ERA or ERstem by ± 10 %. Increasing or decreasing ERsoil had the largest impact, where median EReco increased or decreased to 1.00 or 1.17 mol mol−1, respectively. Also here, the interannual difference between years was very small. A similar pattern could be found for the annual mean EReco, which varied between 1.04 and 1.15 mol mol−1 depending on ERA and ERstem and varied even between 1.00 and 1.24 mol mol−1 due to ERsoil.

The median and mean of hourly O2 : CO2 net assimilation ratio (ERAn) were 0.99 and 0.97 mol mol−1, respectively, for all growing seasons during the simulation period and did not vary between years. In the sensitivity analysis, ERAn was only slightly impacted by changes in the model parameter of ERA (ERstem and ERsoil had no impact). Again, the seasonal and diel variations in ERAn in the year 2012 of the original simulation are shown in Fig. 3b as an example. During nighttime, ERAn was equivalent to ERrd and thus also dependent on Tleaf (Fig. 1b). With low Tleaf at the beginning or end of the growing season, ERAn was often smaller than 0.90 mol mol−1. During daytime, when the magnitude of FA was usually much larger than the magnitude of the opposing flux Frd, ERAn was negatively correlated to Tleaf. Note that Frd and ERrd responded differently to Tleaf; that is, Frd was a fraction of Vcmax, which had an optimal temperature at 27 C (Table 1), while ERrd was positively correlated with Tleaf (Fig. 1b). Consequently, during periods with high Tleaf and low irradiation, Frd was small, but ERrd was large, and the magnitude of the O2 flux of leaf respiration was larger than the magnitude of the CO2 flux with |-Frd ERrd|>|Frd|. Moreover, |-FA ERA| and |FA| were small with ERA=1.00 mol mol−1. It followed that under these conditions and given model implementation, ERAn described the ratio of O2 uptake and CO2 uptake (both fluxes with the same sign), when more O2 was consumed due to leaf dark respiration than released by assimilation (|-Frd ERrd|>|-FA ERA|). In addition, because values of FA were below zero and values of Frd were greater than zero, values of ERAn (Eq. 3) laid mostly not between ERA and ERrd. Similar to EReco, high variations in ERAn were usually found during transition periods with low flux magnitudes.

Figure 3Temporal variations in (a) the exchange ratio of net ecosystem fluxes (EReco, mol mol−1) and (b) the exchange ratio of net assimilation (ERAn, mol mol−1) by hour of day and day of year in 2012. The exchange ratios were calculated as the ratio of the hourly FO2 and FCO2 (including storage terms) summed up over the entire canopy height. As a guide, 1 July is day 183.


3.3 Vertical profiles of O2 : CO2 flux and mole fraction ratios

The vertical profiles of EReco and ERconc differed temporally and spatially. Figure 4 shows the diel vertical profiles of ERecoz and ERconcz averaged over all growing seasons from 2012–2016 (between leafout and leaffall_complete). The mean diel ERecoz ranged from 0.985 to 1.10 mol mol−1 (Fig. 4a). ERecoz at the ground and bottom layers (z/ht0.35) showed very little variability across the day, reflecting the dominance of stem and soil respiration with fixed values of ERsoil and ERstem (Fig. 4a). The upper levels of the canopy showed ERecoz between 0.99 and 1.04 mol mol−1 during the daylight period (06:00 to 20:00; all times are UTC + 1) due to the dominating fluxes of assimilation and stem respiration. The leaf dark respiration did not have a large impact on averaged daytime ERecoz. Moreover, the defined LAI and WAI distributions (Fig. 1a) were represented in the vertical profile of ERecoz, whereas the top canopy contained a larger proportion of sunlit leaves (z/ht>0.75) than the middle part (0.35 <z/ht< 0.75). Hence, ERecoz in the top canopy was influenced more by fluxes of assimilation in daytime hours and was close to 1.00 mol mol−1. Between z/ht=0.3 and z/ht=0.5, ERecoz was larger than 1.06 mol mol−1 during daytime due to higher respiratory processes than assimilation affected by low radiation and relatively high temperatures. The ERecoz during nighttime (approximately before 06:00 and after 20:00) of the upper and middle canopy was usually larger than 1.04 mol mol−1 due to respiratory fluxes.

Figure 4Comparison of the diel dynamics of the height-dependent O2 : CO2 flux and mole fraction ratios averaged over all growing seasons (day of year 110 to 320) from 2012–2016. (a) Vertical profile of the O2 : CO2 flux ratio inside the canopy (ERecoz, mol mol−1), including the whole canopy domain and the soil component (z/ht=0). (b) Vertical profile of the O2 : CO2 mole fraction ratio inside the canopy (ERconcz, mol mol−1), including the whole canopy domain. The exchange ratios for specific canopy heights were derived as the slope of linear regressions fitted to O2 and CO2 fluxes or dry-air mole fractions of multiple simulated time steps for each canopy layer.


The mean diel ERconcz showed relatively small variations ranging from 1.115 to 1.15 mol mol−1 (Fig. 4b) and thus closely matched the prescribed atmospheric O2 : CO2 mole fraction slope of 1.15 (Table 1). Especially during nighttime (before 06:00 and after 20:00), ERconcz was mainly driven by the atmospheric O2 and CO2 background levels. However, bottom layers showed slightly lower values of ERconcz, down to 1.12 mol mol−1, owing to an accumulation of CO2 close to the soil surface produced by soil respiration and low turbulence. During daytime, the canopy air column was well mixed due to stronger turbulence. Nevertheless, ERconcz values were slightly lower in the top canopy layers towards late afternoon and sunset, caused by prevailing canopy respiration.

3.4 Evaluation of the flux-gradient method to obtain O2 fluxes

The vertical profiles of air temperature, water vapor, CO2 and O2 mole fractions were modeled for the entire CANVEG domain including 40 canopy layers and 80 atmosphere layers above the canopy. Figure 5 shows examples of vertical profiles for 12:00 to 13:00 (daytime) and 23:00 to 00:00 (nighttime) on 4 July 2012, an arbitrarily chosen sunny day. Generally, during daytime the vertical profiles within the canopy (Fig. 5a and c) were mostly induced by radiative transfer, leaf photosynthesis, transpiration and autotrophic respiration, which were influenced by the vertical LAI and WAI distributions (Fig. 1a). Furthermore, soil evaporation and respiration resulted in higher water vapor and CO2 mole fractions close to the soil surface. For the layers above the canopy (z/ht>1), the profiles changed monotonically. Daytime O2 and CO2 profiles (Fig. 5c) showed a mirrored shape because the O2 and CO2 fluxes were contributing inversely to the atmospheric mole fractions. Nighttime water vapor and CO2 profiles (Fig. 5d and d) showed a continuous decrease with height and the O2 profile a continuous increase due to the dominance of soil evaporation and soil, with stem and leaf respiration in the lower layers being a sink for O2. During nighttime, air temperature (Fig. 5b) was slightly lower at the canopy top than inside the canopy due to higher energy loss by the emission of longwave radiation.

Figure 5Vertical profiles of (a, b) air temperature (T) and water vapor (v) and (c, d) CO2 and O2 mole fractions of the entire model domain, where O2 mole fractions are shown as the difference from 209 750 ppm (EquO2, 209 750 ppm was derived as the intercept of the relationship between measured atmospheric O2 and CO2 mole fractions; see Table 1). (e, f) ΔFO2,(c,T,v) that resulted from Eq. (7) (Sect. 2.4). The left panels (a)(c) and (e) show mean profiles for 12:00 to 13:00 (daytime) and the right panels (c), (d) and (f) for 23:00 to 00:00 (nighttime), all for 4 July 2012. The flux-gradient method was applied for the gradients between a top measurement height at z/ht= 2 and each layer below and is based on profiles and fluxes of CO2, H and LE (ΔFO2,c, ΔFO2,T and ΔFO2,v).


Based on these modeled vertical profiles and the corresponding flux (FCO2, H or LE, respectively), O2 fluxes were calculated with the flux-gradient method and compared to the modeled O2 fluxes from CANVEG, both corrected for the storage term. So in the following we always describe the ecosystem fluxes (turbulent fluxes plus storage terms). Figure 5e and f show the difference between the various flux-gradient methods derived and modeled FO2 (ΔFO2,(c,T,v); Eq. 7) for the respective simulation hours, when the scalar gradients were derived from two heights (Sect. 2.4). An FO2 estimate and a ΔFO2,(c,T,v) value were obtained for each layer. Generally, ΔFO2 derived with the flux-gradient method based on the CO2 profile (ΔFO2,c) was lower than ΔFO2 derived from the temperature and water vapor profile (ΔFO2,T, ΔFO2,v; Fig. 5e and f). For daytime conditions (Fig. 5e), the mean ΔFO2,c, ΔFO2,T and ΔFO2,v above the canopy were 0.030 ± 0.09, 1.55 ± 0.54 and 4.26 ± 0.63 µmol m−2 s−1, respectively (Table 3). There was little vertical variation in ΔFO2,(c,T,v) above the canopy for nighttime (Fig. 5f). Here, the mean ΔFO2,c, ΔFO2,T and ΔFO2,v were 0.53 ± 0.04, 1.98 ± 0.20 and 0.47 ± 0.24 µmol m−2 s−1, respectively. By applying the three-height flux-gradient method based on Faassen et al. (2023), ΔFO2 for the daytime hours had a similar magnitude for ΔFO2,c at 0.13 µmol m−2 s−1 and for ΔFO2,v at 4.31 µmol m−2 s−1 and was larger for ΔFO2,T at 4.72 µmol m−2 s−1. The corresponding nighttime ΔFO2,c, ΔFO2,T and ΔFO2,v derived from the three-height flux-gradient method were 0.50, 2.41 and 0.66 µmol m−2 s−1, indicating the similar ΔFO2,c and ΔFO2,v but a larger magnitude of ΔFO2,T than with the two-height flux-gradient method.

Table 3Difference between the FO2 estimations derived by the flux-gradient method (FO2,(c,T,v), based on FCO2, H or LE and their respective vertical scalar profile) and by model simulations (FO2,CANVEG) for above-canopy fluxes and for day- and nighttime individually. Results of the two-height approach are shown as the mean and standard deviation of flux gradients derived between z/ht= 2 and each layer below or above the canopy. Also results of the three-height approach are shown, where the flux gradient was derived between three fixed heights (z/ht= 1.05, 1.45 and 2 with ht = 37.5 m).

Download Print Version | Download XLSX

The ΔFO2,(c,T,v) within the canopy during daytime increased and was highly variable for all three methods due to the presence of sources and sinks, as well as nonlinearity of the gradients (Fig. 5e). ΔFO2,c and ΔFO2,T showed hyperbolic shapes with very low (<50 µmol m−2 s−1) and high values (> 50 µmol m−2 s−1), where the CO2 dry-air mole fractions or the temperatures, respectively, were very close to the conditions at the top measurement height, and so the gradients were very small. The sudden jumps from large positive to large negative values were caused by the change in signs of Δc and ΔT.

To guarantee a large gradient, the heights with z/ht=2 and z/ht=1.05 were used in inferring FO2 from vertical CO2, temperature and water vapor gradients for the following analysis. Figure 6a, b and c show the median diel courses of ΔFO2,c, ΔFO2,T and ΔFO2,v for all growing seasons from 2012–2016. Assuming that with these heights the gradients were large enough, the inferred FO2 agreed well with modeled FO2 for ΔFO2,c throughout the median diel course, ranging from 0.45 to 0.15 µmol m−2 s−1 (Fig. 6a). The medians of ΔFO2,T and ΔFO2,v indicated that FO2,T was overestimated by up to 1.59 µmol m−2 s−1 and FO2,v underestimated by up to 5.43 µmol m−2 s−1 during daytime hours (Fig. 6b and c). The standard deviations of ΔFO2,(c,T,v) reflected the diel variation in turbulent conditions and vertical gradients, which were also dependent on the eddy diffusivity. The nighttime standard deviation of ΔFO2,v was relatively large but smaller for ΔFO2,T. The latter produced more outliers during daytime, especially during times of sunrise and sunset. The standard deviation of ΔFO2,c was relatively low and usually < 10 µmol m−2 s−1 across all times of the day except at 08:00, 12:00 and 19:00 (Fig. 6a).

Figure 6(a, b, c) Median diel cycles of the differences between O2 fluxes derived by the flux-gradient method and by CANVEG simulation (ΔFO2,(c,T,v)) for all growing seasons from 2012–2016. The flux-gradient method was applied for the gradients between z/ht= 2 and z/ht= 1.05 and is based on profiles and fluxes of (a) CO2, (b) H and (c) LE (ΔFO2,c, ΔFO2,T and ΔFO2,v). The error bars indicate the standard deviation of ΔFO2,(c,T,v) by hour. (d) Histogram of uncertainties in FO2 (σFO2) derived by the flux-gradient method based on CO2 profiles and fluxes, in which a random uncertainty in O2 mole fractions (Δo) was included. The uncertainty in Δo followed a normal distribution with a mean of 0 and a standard deviation of 0.7 ppm (Pickers et al., 2017). In order to include daytime hours with an active canopy for the estimation of σFO2, Δo 1 ppm was used as a filter, assuming higher-oxygen dry-air mole fractions close to the canopy than in the top domain layers.


The above analysis evaluated the flux-gradient method solely regarding the characteristics and dynamics of various scalar gradients. Moreover, accurate and precise measurements of the scalars are also necessary for the satisfactory performance of this method. We added a random uncertainty to our modeled O2 mole fractions to simulate gradient measurements with the current instrument uncertainty (Δo in Eq. 6). Figure 6d shows the distribution of the differences (σFO2) between the FO2 estimates based on the flux-gradient method including a random measurement uncertainty in Δo or not. For this analysis, only hourly time steps within all growing seasons from 2012–2016 were chosen with Δo 1 ppm, when O2 mole fractions increased with decreasing height above the canopy due to prevailing gross assimilation over respirations during daytime. The median of resulting σFO2 was 0.20 µmol m−2 s−1 and thus very close to zero. Here, we extracted the 10 % and 90 % quantile of σFO2=-14.2 and 14.5 µmol m−2 s−1. Thus, we used 15 µmol m−2 s−1 as the upper limit of σFO2 in the evaluation of the flux-partitioning approach (Sect. 3.5).

3.5 Uncertainties in partitioning net ecosystem CO2 fluxes based on O2 fluxes

For the test day from 07:00 to 19:00 on 4 July 2012, model output of hourly FO2 was used to derive the main CO2 flux components. The a posteriori uncertainties of the partitioned fluxes of gross assimilation (σFA) decreased significantly with decreasing uncertainties of σERA and σFO2, indicating the importance of reducing errors in ER and O2 flux measurements (Fig. 7). The a priori uncertainties had strong effects on a posteriori uncertainties because a large σFAb allowed large |FA-FAb| to reach a minimum J value and vice versa (Eq. 8). Without the constraints of a priori uncertainties (Fig. 7a), σFA reached 193 µmol m−2 s−1 at its maximum and then reduced with smaller σFO2 and σERA to 28 µmol m−2 s−1, which was still larger than the a priori value (Table 2). If a priori uncertainties (σFAb, σFRb, σERAb, σERRb) were included (Fig. 7b and c), σFA was much lower. When assuming an uncertainty for the net CO2 fluxes (σFCO2) of 2.5 µmol m−2 s−1, σFA showed very little variation and ranged between 4.74 and 4.88 µmol m−2 s−1, remaining close to the minimum of the chosen a priori uncertainty in FA and FR (Fig. 7b). When assuming more accurate FCO2 and ERR measurements with σFCO2=0.5µmol m−2 s−1 and σERR=0.001 mol mol−1, σFA was reduced to a minimum of 1.43 µmol m−2 s−1 (Fig. 7c). A moderate level of a priori uncertainties in O2 fluxes and ERA (bold black lines in Fig. 7c) resulted in σFA=4.48µmol m−2 s−1 for our test day. In this case, the partitioned FA was 28.3 µmol m−2 s−1, which was about 6 % lower than the estimated gross assimilation obtained with the eddy covariance technique (FAb=30.2µmol m−2 s−1). With regard to the sensitivity analysis, σFA was only slightly impacted by ERA. σFA ranged from 1.42 to 4.83 µmol m−2 s−1 for the case of the lower a priori uncertainty (with σFCO2=0.5µmol m−2 s−1 and σERR=0.001 mol mol−1).

Figure 7Uncertainty in partitioned gross assimilation CO2 flux (FA) determined from eddy covariance net ecosystem CO2 flux (FCO2) with net ecosystem O2 flux (FO2), O2 : CO2 ratio of gross assimilation (ERA) and ecosystem respiration (ERR) on 4 July 2012. (a) Optimized a posteriori uncertainty of FA (σFA) without a priori FA values and uncertainties. (b) Optimized σFA including all of the a priori terms in the J function as written in Eq. (8), with a priori uncertainty of FCO2 (σFCO2) = 2.5 µmol m−2 s−1 and a priori uncertainty of ERR (σERRb) = 0.05 mol mol−1. (c) Same cost function as for (b) but with σFCO2=0.5µmol m−2 s−1 and σERRb= 0.001 mol mol−1. The bold black lines show the practical optimization test with σERAb and σFO2 around 0.01 mol mol−1 and 15 µmol m−2 s−1, respectively (see Fig. 6d).


4 Discussions

4.1 Model setup and model performance

We added O2 : CO2 exchange ratios and O2 flux processes into the one-dimensional, multilayer atmosphere–biosphere gas exchange model CANVEG. To represent natural atmosphere–ecosystem exchange satisfactorily, we first calibrated and validated the model based on eddy covariance CO2 and energy flux observations from a temperate deciduous forest in Leinefelde, Germany, from 2012–2016. In a previous study, model performance was evaluated based on hourly CO2, water vapor and energy fluxes in temperate oak forests (Baldocchi and Wilson, 2001). That evaluation, for hourly FCO2, yielded a slope of 1.09 of the regression between observations and simulation with R2= 0.82, which is comparable to our results (slope of 1.02 and R2=0.82; Fig. 2b). The model application in a deciduous temperate forest in central Germany (Knohl and Baldocchi, 2008) also showed a high match between hourly modeled and measured FCO2 (slope of 0.997, R2= 0.857). In addition, Hanson et al. (2004) compared the CANVEG model with seven other stand-level models where CANVEG performed very well (slope of 0.93, R2=0.82) based on simulated FCO2. In our study, the comparison between hourly LE simulation and observations obtained a regressed slope of 1.02 and R2=0.77 (Fig. 2f), indicating a better model performance than for daily evapotranspiration by Hanson et al. (2004) (slope of 1.17, R2=0.73). Knohl and Baldocchi (2008) found a slope of 0.926 and R2=0.825 for hourly LE simulation and a slope of 1.021 and R2=0.869 for hourly H simulation, indicating an underestimation of LE and a small overestimation of H. In our study, we observed an overestimation of LE and underestimation of H. The model performance (with regard to the slope, R2 and RMSE) in the energy fluxes was generally lower than for CO2 flux simulations because fitted parameters mainly affected the CO2 fluxes and leaf assimilation (Table 1). By adjusting the assimilation rate, only transpiration was changed, which then had an impact on LE and H. The non-unity slope of H and LE could also point to the non-closure of the energy balance in the eddy covariance observations.

Furthermore, the modeling error could be caused by the implemented soil respiration algorithm, which did not consider the influence of soil water changes. Moreover, parameters for soil respiration were only calibrated based on eddy covariance observations (FCO2 and FR) on the ecosystem scale, where independent chamber measurements would be beneficial. Moreover, an error in the seasonality of carbon and energy fluxes could be introduced by the uncertainty in leaf growth phenology and annual LAI. Although we simulated fluxes from 2012–2016, the total full-leaf LAI and leaf growth phenology parameters (Table 1) were only measured in the year 2015 and kept constant across the modeling period (Table 1). Adjusting LAI annually would only affect the timing of the fluxes but not the overall O2 : CO2 exchange ratio (ER) pattern.

This study used fixed ER parameter values owing to the lack of direct chamber O2 and CO2 flux measurements for leaf, stem and soil flux components at our study site. The O2 : CO2 exchange ratio of gross assimilation (ERA) was set to 1.00 mol mol−1 (Table 1), describing the production of carbohydrates by gross assimilation. Busch et al. (2018) described how plants use nitrogen while assimilating CO2, resulting in carbon loss from the photorespiratory pathway in the form of glycine and serine. Since nitrogen assimilation increases O2 emissions but has smaller effects on CO2 uptake, incorporating nitrogen assimilation in the Farquhar et al. (1980) photosynthesis model would help to represent photosynthetic O2 emissions more mechanistically in models. In this case, environmental conditions such as nitrogen fertilization and utilization would cause different ERA values.

Studies obtaining exchange ratios of O2 and CO2 via chamber measurements at the soil or stem scale often state the so-called apparent respiratory quotient (ARQ), which is defined as the ratio of CO2 efflux to O2 uptake (Angert et al., 2012; Helm et al., 2021; Hilman and Angert, 2016; Hilman et al., 2022). Thus, ARQ could be compared to our ERsoil or ERstem by taking the inverse of ARQ, which is the CO2 : O2 conductance ratio, following Hilman and Angert (2016). However, ARQ is also influenced by biotic and abiotic non-respiratory processes such as dissolution and refixation of respired CO2 in the xylem sap (Angert et al., 2012; Hilman and Angert, 2016; Hilman et al., 2022), so we expect differences between the various quantities. Furthermore, studies discuss the so-called oxidative ratio (OR) based on the elemental analysis of organic material. OR is based on the stoichiometry of the respiratory product or net synthesized biomass, which represents the oxidation state of respiratory substances (Hilman et al., 2022; Juergensen et al., 2021).

All ARQ values from the cited references were converted to ERstem or ERsoil for easier comparison. The ERstem parameter of 1.04 mol mol−1 used in this study was derived by Randerson et al. (2006) based on the OR of chemical compositions (lipid, lignin, protein, soluble phenolic, etc.) assigned to woody stems. Hilman and Angert (2016) measured a mean ERstem=1.47 mol mol−1 (ARQ = 0.68 ± 0.04 mol mol−1) with direct continuous measurements for an apple tree. In addition, ERstem also showed variations between 1.22 and 1.61 mol mol−1 (ARQ = 0.62 to 0.82 mol mol−1) during the measurement period (Hilman and Angert, 2016). The ERstem varied between 1.28 and 2.56 mol mol−1 (ARQ = 0.39 to 0.78 mol mol−1) with the mean of 1.69 mol mol−1 (ARQ = 0.59 mol mol−1) among tropical, temperate and Mediterranean forests (Hilman et al., 2019). In addition, dry or wet environmental conditions lead to a seasonal variation in ERstem (Angert et al., 2012).

The global OR of soils is suggested to be equal 1.10 ± 0.05 (Severinghaus, 1995). According to Hockaday et al. (2015), the soil OR is 1.006 at ambient CO2 level and increases to 1.054 with elevated CO2 level. Worrall et al. (2013) also derived a global soil OR = 1.04. Seibt et al. (2004) obtained an ERsoil=0.94 mol mol−1 with field chamber measurements, while Ishidoya et al. (2013) obtained ERsoil=1.11 mol mol−1. ERsoil also showed seasonal variations from about 1.11 mol mol−1 (ARQ = 0.9 mol mol−1) during late spring and summer to about 1.43 mol mol−1 (ARQ = 0.7 mol mol−1) during winter in a Mediterranean mixed conifer forest (Hicks Pries et al., 2020). Depending on ecosystem type, such as alpine areas, temperate, Mediterranean or tropical forests, and on sampling strategies, such as sampling of soil air or bulk soil, measured ERsoil varied between 0.88 to 4.35 mol mol−1 (ARQ = 0.23 to 1.14 mol mol−1) (Angert et al., 2012, 2015; Hilman et al., 2022). These variabilities related to seasons, forest types and ecosystem processes strongly suggest that site specific ERstem and ERsoil should be used in O2 flux simulations. A logarithmic relationship between soil ARQ and soil temperature, as found by Hilman et al. (2022), could also be introduced to future soil O2 flux models. Due to this high variance between derived ERs of these different studies, we conducted a sensitivity analysis by changing ERA, ERstem or ERsoil by ± 10 % to show how these parameters affected the modeled FO2, EReco and ERAn. Furthermore, we assessed the impact of these model parameters on the source-partitioning results. In summary, the model simulations showed a small sensitivity towards the model parameter settings. The modeled FO2 sum was mostly sensitive to ERA, which corresponded to the largest flux component. EReco and ERAn changed by less than 10 % in each case. The uncertainty in the source-partitioning results were mostly driven by the uncertainty of O2 flux estimates (σFO2) and much less by the ER parameters. Generally, all model simulations yielded the same tendency and pattern of exchange ratios.

4.2 Temporal and vertical dynamics of O2 : CO2 exchange ratios

The O2 : CO2 flux exchange ratio (EReco) quantifies the simultaneous canopy–atmosphere gas exchange of the whole ecosystem. We obtained EReco by aggregating simulated O2 and CO2 fluxes of all canopy layers and taking the ratio or by deriving the slopes of linear regressions fitted to O2 and CO2 fluxes of multiple simulated time steps for each canopy layer (ERecoz). The temporal variations in EReco arose from diel and seasonal variations in the flux contributions of gross assimilation and respiration to net ecosystem O2 and CO2 exchange. Since assimilation and respiration are two individual processes which are influenced by two differing main drivers – photosynthetic photon flux density and temperature – they usually show shifted diel cycles. Furthermore, fluxes from respiration consist of various components originating from various sources (e.g., respiration by heterotrophs, leaves or roots), which can also differ in their diel cycles, in their ER and in their proportions of total O2 and CO2 ecosystem fluxes. Further studies should obtain ER independently with respective chamber measurements in order to separate environmental effects (e.g., radiation, temperature, humidity) on each componential O2 and CO2 flux.

The EReco contains information about the turbulent flux exchange, as well as the O2 and CO2 storage terms between soil surface and measurement height. Our study focused on the whole ecosystem O2 and CO2 exchange ratio including storage terms. Annual mean EReco ranged from 1.06 to 1.12 mol mol−1 within the 5 years, and estimates of ERecoz varied between 0.99 and 1.10 mol mol−1 with height in the canopy (Fig. 4a). Seibt et al. (2004) reported daytime net turbulent ER (considering turbulent fluxes and not including storage terms) between 1.26 and 1.38 mol mol−1, which they derived with a one-box model. Next to the inclusion or exclusion of storage terms and the usage of different models, differences between Seibt et al.'s (2004) work and ours could also be caused by the difference in considered time periods: our simulations covered 5-year growing seasons of O2 and CO2 fluxes between the canopy and the atmosphere, and Seibt et al. (2004) focused on July and August between 1999 and 2001. Moreover, we used different componential ER parameters (Table 1) in our simulations.

Diel ERAn variations reflected separate responses of gross assimilation and leaf dark respiration to temperature. The median and mean of hourly ERAn were 0.99 and 0.96 mol mol−1, respectively, for all growing seasons during the study period. However, ERAn showed extreme values during transition hours with low flux magnitudes (Fig. 3b). Ishidoya et al. (2013) found ERAn values close to 1.02 mol mol−1 via leaf chamber measurements. According to Seibt et al. (2004), ERAn ranged between 1.04 and 1.20 mol mol−1, observed also via chamber measurements when flux rates were between 2 and 5 µmol m−2 s−1. A lower flux rate (1.7 µmol m−2 s−1) led to a higher variability in ERAn (Seibt et al., 2004). The divergence between our ERAn estimates (which were close to 1.00 mol mol−1) and to the chamber measurements could be caused by the utilization of varying nitrogen sources that would increase ERAn (Seibt et al., 2004).

The mole-based O2 : CO2 exchange ratio (ERconc) is determined by the atmospheric background mole fractions of O2 and CO2, by the distributions and dynamics of sources and sinks, and by the turbulence inside the canopy. ERconc is usually derived based on the slopes of Deming regressions of observed O2 and CO2 mole fractions, accounting for uncertainty in both variables (Battle et al., 2019; Ishidoya et al., 2020). Our results of ERconc and EReco confirmed that ERconc cannot represent simultaneous O2 and CO2 exchange as EReco, which was also recently found by Faassen et al. (2023). We also estimated ERconcz for each canopy layer representing O2 and CO2 mole fractions of air at certain canopy heights. The mean diel ERconcz showed only very small variations ranging from 1.12 to 1.15 mol mol−1 within the diel course. Battle et al. (2019) observed an average ERconc=1.081± 0.007 mol mol−1 in a mixed deciduous forest over a 6-year period and ERconc=1.03± 0.01 mol mol−1 on two summer days in July 2007. Their ERconc measurements also showed temporal variations on a 6 h basis between 0.85 and 1.15 mol mol−1. Seibt et al. (2004) measured and modeled ERconc during day- and nighttime at several sites and obtained values varying between 1.04 and 1.19 mol mol−1. Ishidoya et al. (2013) observed daily average ERconc=0.94± 0.01 mol mol−1, with daytime ERconc=0.87± 0.02 mol mol−1 and nighttime ERconc=1.03± 0.02 mol mol−1. Ishidoya et al. (2013) also built a one-box canopy O2 : CO2 budget model applying the same parameter values ERA=1.00 mol mol−1 and ERR=1.10 mol mol−1 as our study. Their observed daytime ERconc=0.87 mol mol−1 agrees with their modeled net turbulent ER = 0.89 mol mol−1. Our modeled ERconcz estimates showed a lower temporal variability within the mean diel course than in the cited studies. This is to a large part due to background O2 that was fixed to 1.15 of atmospheric CO2 mole fractions (Table 2). One would expect, though, that this ratio might be lower during summer and most probably has also a diel cycle. Future work could include continuous measurements at the site resulting in a varying background value and potentially larger diel and seasonal variability. It is also possible that mixing in CANVEG was too strong, with the result that modeled ERconcz was excessively influenced by the background value. This could be improved in future by comparing modeled temperature, water vapor and CO2 mole fractions with measured mole fractions at different canopy heights, which have become standard measurements at eddy covariance sites in forests now.

4.3 Estimation of ecosystem O2 fluxes and applications

Eddy covariance measurements, as typically conducted for CO2 fluxes, are currently not possible for O2 fluxes because no sufficiently fast and precise O2 analyzer is commercially available yet (except for a self-made, non-commercial vacuum ultraviolet (VUV) absorption analyzer developed by Stephens et al., 2003). Requirements would be a precision of below 1 ppm against a background concentration of 210 000 ppm on a high, turbulence-resolving measurement frequency (Keeling and Manning, 2014). However, vertical profiles of air temperature, water vapor, CO2 and O2 mole fractions can already be obtained with high precision. With our modeled vertical profiles, we determined O2 fluxes based on the flux-gradient approach, testing various profile setups and the necessary instrument precision for O2 mole fraction measurements (Figs. 5 and 6). By choosing various heights to derive the mole fraction gradients, we confirmed that the selected heights should both be above the canopy. This guarantees that the profiles are differentiable as there are no sources or sinks between sampling heights and that the eddy diffusivity of O2 is the same as that of the other corresponding scalars (Baldocchi et al., 1988). In addition, the mole fraction difference between the two heights should be as large as possible to decrease the uncertainty in O2 flux estimates. Here, we selected, amongst others, heights at z/ht= 1.05 and 2 to obtain large gradients. Faassen et al. (2023) applied the flux-gradient method to estimate O2 fluxes in a boreal forest with a canopy height of 18 m. Their measurements were conducted between 23 and 125 m for the vertical scalar gradient, reaching about 7 times the canopy height. Such a large distance between measurement heights in a profile system is usually only feasible for cropland, grassland or peatland study sites with low vegetation. For high vegetation, such as forest sites, a tall tower is needed (as in Faassen et al., 2023). However, by choosing two measurement heights with a large distance (e.g., multiple tens of meters), the difference between the footprint extensions of each height also becomes large, potentially resulting in erroneous flux estimates. If, for instance, the vertical CO2 gradient could be doubled, the uncertainty in FO2 fluxes caused by the measurement uncertainty of O2 gradients would be reduced by 50 % according to Eq. (6).

The median differences between FO2 derived with the flux-gradient method and modeled FO2 (ΔFO2,(c,T,v)) were generally < 5.5 µmol m−2 s−1, independent of which scalar concentrations and fluxes were used for the latter. However, ΔFO2,v and ΔFO2,T deviated more from zero during daytime, indicating that FO2 estimates based on LE and the water vapor profile and on H and the temperature profile would lead to underestimation or overestimation, respectively, during daytime by the flux-gradient method (Fig. 6). The FO2 estimates during nighttime were more uncertain based on temperature and water vapor, as indicated by large standard deviations. These “outliers” occurred due to vertical gradients that were too small, caused by a low activity of sources and sinks and/or insufficient turbulence. The flux-gradient method based on CO2 mole fractions and fluxes yielded FO2 estimates in better agreement with modeled FO2. But this was probably because the O2 sources and sinks were highly correlated to CO2 processes due to the O2 modeling setup and constant ER (Eq. 1). Consequently, it is still recommended to use all the available gas or energy gradients to derive O2 fluxes with the flux-gradient methods, and then choose the most appropriate method (if this is possible) for various times during the day or year depending on the magnitude of the gradients, the quality of flux measurements and the turbulence. The magnitude of the gradients could additionally be increased for each scalar by choosing scalar-specific measurement heights.

The flux-gradient method has already been used for O2 flux estimation above a cool temperate forest (Ishidoya et al., 2015), an urban canopy (Ishidoya et al., 2020) and a boreal forest (Faassen et al., 2023). The latter study applied a three-height flux-gradient approach, where they estimated the eddy diffusivity K based on CO2 and temperature measurements at three heights and applied a vertical O2 gradient between two heights. We also tested this three-height flux-gradient approach based on our model simulations, but we assumed that all scalars including O2 were measured at three heights. Based on our simulations, we could not observe an improvement in the flux estimation due to the inclusion of three measurement heights in the flux-gradient method instead of two heights.

Uncertainty of O2 mole fraction estimates resulted in a median close to zero for the uncertainty σFO2. The uncertainty in O2 mole fraction estimates was selected randomly following a normal distribution in the model simulations. Our analysis showed that the flux-gradient method has the potential for FO2 estimation, but we also found that estimated FO2 could be over- or underestimated by up to ± 5.5 µmol m−2 s−1. To make the flux-gradient method more precise, the vertical scalar gradient should be as large as possible and flux and profile measurements as precise as possible. To achieve this, on the one hand, a larger distance between measurement heights is needed (not possible over large forest stands but applicable for crop-, grass- and peatland), and on the other hand, a higher measurement precision is necessary to reduce the uncertainty in scalar gradient measurements.

In general, mass is transported in air due to diffusive and non-diffusive processes. Diffusive transport can be induced due to random turbulent or molecular motions acting against a gradient. As shown in Fig. 5, an exemplary vertical profile or gradient of CO2 mole fraction regarding dry air shows a higher mole fraction close to the soil surface due to respiratory processes and a lower mole fraction within the forest canopy due to net assimilation during daytime. Above the canopy the CO2 dry-air mole fraction increases slightly again within the boundary layer. The vertical O2 profile is mirrored to this CO2 profile (when dry-air mole fractions are considered). Because of the processes of evaporation and transpiration from the soil surface and canopy, water vapor is also added to the air column, and the vertical water vapor profile usually shows a decreasing water vapor mole fraction with increasing height. The addition of water vapor molecules to an air package dilutes the other molecules in that air package such as N2, O2 and CO2 by replacing some of them. Thus, the ratio between the number of O2 or CO2 molecules and total number of air molecules (equal to mole fraction regarding moist air) decreases, and therefore the vertical O2 and CO2 gradients change. Furthermore, due to the addition of water vapor molecules, other air molecules are being displaced and moved away from the evaporating surface. This displacement effect yields a non-diffusive transport (also known as Stefan flow) that does not necessarily follow a gradient (Kowalski, 2017; Kowalski et al., 2021). The magnitudes of the dilution and displacement effects depend on the mass fraction of each gas (number and weight of molecules per mass of air), where O2 is more affected than CO2 due to its high abundance (Kowalski et al., 2021). Considering the above-described vertical profile, O2 diffuses downwards towards the evaporating surface following the increased gradient due to the dilution effect. However, this downward motion can be offset by the displacement effect.

To analyze the transport of and the relationship between O2 and CO2 molecules, the dilution and displacement effects have to be considered – also in relation to the turbulent transport. The magnitudes and directions of diffusive (turbulence and molecular diffusion) and non-diffusive transport are variable and need to be quantified experimentally for various atmospheric conditions, ecosystems and heights above the ecosystems. Thus, the significance and impacts of the various transport types are unknown and currently under discussion. Considering the many open questions regarding non-diffusive transport, we have not implemented the Stefan flow within CANVEG until now.

The CANVEG model considers mole fractions regarding dry air (removing all the water vapor) for O2 and CO2, and therefore the dilution effect is excluded from the model simulations, and vertical gradients do not change due to the process of evapotranspiration. This allows for the comparison with O2 measurements, where it is common practice to cryogenically dry the air before analysis for O2 (Pickers et al., 2017). The non-diffusive transport (Stefan flow) would play a role in our study within the application of the flux-gradient method and the estimation of ERconc. By the modification of the vertical gradients due to the non-diffusive transport, flux estimates based on the flux-gradient method would differ (personal communication with Andrew Kowalski). However, our study considered mostly net ecosystem fluxes in this application. Further, Kowalski et al. (2021) determined that the Webb, Pearman and Leuning (WPL) methodology, based on perturbations in the dry-air mass fraction, correctly estimated biogeochemical fluxes (for both water vapor and CO2) despite incorrectly describing transport mechanisms. Therefore, the WPL methodology predicts that artificially eliminating the effects of water vapor (dilution and displacement) and expressing each gas with reference to dry air will yield the equivalent flux-gradient relationships. Furthermore, by assuming all scalars (temperature, water vapor, CO2 and O2) are transported similarly (and thus assuming the eddy diffusivities Ko, Kc, KT and Kv are the same), we have added an additional uncertainty. Also due to the change in the vertical gradients, the estimation of ERconc will be affected because the displacement by evapotranspiration has a different impact on CO2 and O2. However, again for the mole fractions regarding dry air, the effect should be small. Also, the estimated ERconc (and also EReco) were reasonable and in line with the current understanding of the process.

Lastly, in our study, FCO2 obtained with the eddy covariance technique was source partitioned based on simulated FO2, and the uncertainty in gross assimilation (σFA) was evaluated. By estimating CO2 flux components following the same approach based on stable isotopes in CO2, Knohl and Buchmann (2005) derived a σFA for instantaneous half-hourly data of 6 µmol m−2 s−1 assuming CO2 and 13CO2 flux uncertainties of 0.5 µmol m−2 s−1 and 25 ‰ µmol m−2 s−1, respectively. A σFA of around 4 µmol m−2 s−1 was found with a higher uncertainty in CO2 fluxes (σFCO2) of 2 µmol m−2 s−1 by Ogee et al. (2004). Our study obtained comparable results under similar levels of σFCO2=2.5µmol m−2 s−1 (Fig. 7b). However, Ogee et al. (2004) found that σFA can be reduced to 2 µmol m−2 s−1 when an isotopic disequilibrium is larger than 0.004 (Fig. 6 in Ogee et al., 2004, p. 11). We obtained ERA and ERR disequilibrium (|ERA ERR|) of around 0.086 mol mol−1 but still could not improve our σFA under normal σFCO2 levels (Fig. 7b). This was probably because our uncertainty in O2 fluxes (σFO2) was much larger (up to 15 µmol m−2 s−1) relative to the ERA and ERR disequilibrium (0.086 mol mol−1) (Ogee et al., 2004). Here, we derived σFO2 based on the analysis of the flux-gradient application (Fig. 6d). Thus, a higher precision in FO2 estimates and/or a larger ERA and ERR disequilibrium  0.086 mol mol−1 is needed. In our simulations, the disequilibrium had low variation due to the fixed ER parameters. A small variation was only introduced by the variable ERrd due to leaf temperature. Implementing variable ER parameters (depending on, for example, environmental conditions) in the model or obtaining real ER values by measurements could be beneficial. Figure 7c shows that by improving the precision in FCO2 and ERR estimates and not so much in FO2 also yields a lower uncertainty in estimates of gross assimilation. Faassen et al. (2023) changed EReco by ± 0.20 mol mol−1, which resulted in a change in partitioned FA of 6.7 %. However, they compared their partitioned flux components with FA and FR derived using the eddy covariance method by assuming that the latter describes the “true” values. They emphasize the importance of a correct estimate for EReco. Our evaluation of σFA was based on assigned a priori uncertainties to all elements which were independent of the flux values (Ogee et al., 2004). Compared with eddy covariance data, our partitioned FA also differed by about 6 %, which is comparable to flux-partitioning results of Faassen et al. (2023).

5 Conclusions

We implemented O2 : CO2 exchange ratios in the CANVEG multilayer ecosystem–atmosphere gas exchange model to enable hourly ecosystem O2 flux simulations. The simulated ecosystem O2 : CO2 exchange ratio (EReco) showed strong diel and seasonal variations. The annual mean EReco ranged from 1.06 to 1.12 mol mol−1 during the 5-year study period and depended significantly on our assumptions about the fixed model parameters describing the exchange ratios of the ecosystem components: leaves, stem and soil (ERA, ERstem, ERsoil). In particular, changes in ERsoil by ± 10 % yielded annual mean EReco from 1.00 up to 1.24 mol mol−1. We also found that hourly EReco and exchange ratios of net assimilation (ERAn) exhibited high variability during transition periods (e.g., during sunrise and sunset) with low flux magnitudes.

According to our simulations, it is feasible to derive ecosystem O2 fluxes with the flux-gradient approach based on sensible heat, latent heat and CO2 turbulent flux measurements under field conditions, when the vertical gradients are measured between 1.05 to 2 times the canopy height. In particular, the vertical O2 difference should be larger than 1 ppm. However, including uncertainty in O2 mole fraction measurements of 0.7 ppm would increase the uncertainty in O2 flux estimates by up to 15 µmol m−2 s−1. The precision of the source-partitioning application was driven by a priori uncertainties of O2 and CO2 flux and ERA and ERR measurements. With an ERA and ERR disequilibrium (|ERA ERR|) of about 0.086 mol mol−1, the uncertainty of partitioned gross assimilation can be constrained to < 5 µmol m−2 s−1 by narrowing the uncertainty of CO2 measurements and ERR estimates to 2.5 µmol m−2 s−1 and 0.05 mol mol−1. O2 flux measurements and additional information on the exchange ratios of gross assimilation and ecosystem respiration (ERA, ERR), for example obtained by chamber measurements, can thus be used as a source-partitioning approach for net CO2 fluxes.

Our model study highlighted the potential temporal and spatially variability in O2 : CO2 exchange ratios of various ecosystem components and the drivers of O2 fluxes at a forest study site. Furthermore, we provided guidance to micrometeorological approaches, such as the flux-gradient method, to obtain sufficient O2 flux estimates depending on measurement setup and on current instrument precision. We further tested the usage of O2 flux estimates to source partition net CO2 fluxes. Further understanding of the relationship between environmental drivers and O2 fluxes and O2 : CO2 exchange ratios, as well as continuous and long-term observations based on, for example, long-term chamber measurements, will greatly help to improve our ecosystem model and our understanding of the carbon cycle in terrestrial ecosystems.

Data availability

The data sets and model code generated and/or analyzed during the current study are available from the corresponding author.


The supplement related to this article is available online at:

Author contributions

AlK and ACM acquired project funding. YY and AlK designed the work. YY conducted the main analysis, and AnK and FM aided in its interpretation. MC provided the model code. YY primarily wrote the paper with input from all coauthors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This research was funded by the European Research Council under the European Union's Horizon 2020 research and innovation program (grant agreement no. 682512 – OXYFLUX). We thank the staff from the Bioclimatology Group of the University of Göttingen, especially Dietmar Fellert, Frank Tiedemann, Edgar Tunsch and Marek Peksa, for their continuous support in data acquisition and instrument maintenance. We thank Penelope A. Pickers, Emmanuel Blei, Julian Deventer and Mattia Bonazza for building O2 measurement systems and providing atmospheric O2 and CO2 data. We also thank Jelka Braden-Behrens for obtaining leaf area measurements, Ashehad Ali for suggestions on model parameter calibrations, Jan Muhr for interpretations of ER results, and Rijan Tamrakar and Christian Markwitz for preparing the meteorological and eddy covariance data. Also, we thank the forest manager Ulrich Breitenstein for allowing the experimental setup at the Leinefelde site. We also thank Andrew S. Kowalski for his important input and suggestions regarding the process of non-diffusive transport during the open discussion. Finally, we thank the reviewers for their helpful and constructive comments that improved the quality of this paper.

Financial support

This research has been supported by the European Research Council, H2020 European Research Council (grant no. 682512).

Review statement

This paper was edited by Christopher Still and reviewed by two anonymous referees.


Angert, A., Muhr, J., Negron Juarez, R., Alegria Muñoz, W., Kraemer, G., Ramirez Santillan, J., Barkan, E., Mazeh, S., Chambers, J. Q., and Trumbore, S. E.: Internal respiration of Amazon tree stems greatly exceeds external CO2 efflux, Biogeosciences, 9, 4979–4991,, 2012. 

Angert, A., Yakir, D., Rodeghiero, M., Preisler, Y., Davidson, E., and Weiner, T.: Using O2 to study the relationships between soil CO2 efflux and soil respiration, Biogeosciences, 12, 2089–2099,, 2015. 

Anthoni, P. M., Knohl, A., Rebmann, C., Freibauer, A., Mund, M., Ziegler, W., Kolle, O., and Schulze, E. D.: Forest and agricultural land-use-dependent CO2 exchange in Thuringia, Germany, Glob. Change Biol., 10, 2005–2019,, 2004. 

Baldocchi, D.: A Lagrangian random-walk model for simulating water vapor, CO2 and sensible heat flux densities and scalar profiles over and within a soybean canopy, Bound.-Lay. Meteorol., 61, 113–144,, 1992. 

Baldocchi, D.: Measuring and modelling carbon dioxide and water vapour exchange over a temperate broad-leaved forest during the 1995 summer drought, Plant Cell Environ., 20, 1108–1122,, 1997. 

Baldocchi, D., Falge, E., Gu, L. H., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X. H., Malhi, Y., Meyers, T., Munger, W., Oechel, W., U, K. T. P., Pilegaard, K., Schmid, H. P., Valentini, R., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities, B. Am Meteorol. Soc., 82, 2415–2434,<2415:FANTTS>2.3.CO;2, 2001. 

Baldocchi, D. D. and Wilson, K. B.: Modeling CO2 and water vapor exchange of a temperate broadleaved forest across hourly to decadal time scales, Ecol. Model., 142, 155–184,, 2001. 

Baldocchi, D. D., Hincks, B. B., and Meyers, T. P.: Measuring biosphere-atmosphere exchanges of biologically related gases with micrometeorological methods, Ecology, 69, 1331–1340,, 1988. 

Baldocchi, D. D., Fuentes, J. D., Bowling, D. R., Turnipseed, A. A., and Monson, R. K.: Scaling isoprene fluxes from leaves to canopies: Test cases over a boreal aspen and a mixed species temperate forest, J. Appl. Meteorol., 38, 885–898,<0885:SIFFLT>2.0.CO;2, 1999. 

Baldocchi, D. D., Wilson, K. B., and Gu, L. H.: How the environment, canopy structure and canopy physiological functioning influence carbon, water and energy fluxes of a temperate broad-leaved deciduous forest-an assessment with the biophysical model CANOAK, Tree Physiol., 22, 1065–1077,, 2002. 

Battle, M., Bender, M. L., Tans, P. P., White, J. W. C., Ellis, J. T., Conway, T., and Francey, R. J.: Global carbon sinks and their variability inferred from atmospheric O2 and δ13C, Science, 287, 2467–2470, 2000. 

Battle, M. O., Munger, J. W., Conley, M., Sofen, E., Perry, R., Hart, R., Davis, Z., Scheckman, J., Woogerd, J., Graeter, K., Seekins, S., David, S., and Carpenter, J.: Atmospheric measurements of the terrestrial O2 : CO2 exchange ratio of a midlatitude forest, Atmos. Chem. Phys., 19, 8687–8701,, 2019. 

Bowling, D. R., Tans, P. P., and Monson, R. K.: Partitioning net ecosystem carbon exchange with isotopic fluxes of CO2, Glob. Change Biol., 7, 127–145,, 2001. 

Braden-Behrens, J., Yan, Y., and Knohl, A.: A new instrument for stable isotope measurements of 13C and 18O in CO2 – instrument performance and ecological application of the Delta Ray IRIS analyzer, Atmos. Meas. Tech., 10, 4537–4560,, 2017. 

Braden-Behrens, J., Markwitz, C., and Knohl, A.: Eddy covariance measurements of the dual-isotope composition of evapotranspiration, Agr. Forest Meteorol., 269, 203–219,, 2019. 

Busch, F. A., Sage, R. F., and Farquhar, G. D.: Plants increase CO2 uptake by assimilating nitrogen via the photorespiratory pathway, Nat. Plant., 4, 46–54,, 2018. 

Collatz, G. J., Ball, J. T., Grivet, C., and Berry, J. A.: Physiological and environmental-regulation of stomatal conductance, photosynthesis and transpiration – a model that includes a laminar boundary-layer, Agr. Forest Meteorol., 54, 107–136,, 1991. 

Emad, A. and Siebicke, L.: True eddy accumulation – Part 1: Solutions to the problem of non-vanishing mean vertical wind velocity, Atmos. Meas. Tech., 16, 29–40,, 2023a. 

Emad, A. and Siebicke, L.: True eddy accumulation – Part 2: Theory and experiment of the short-time eddy accumulation method, Atmos. Meas. Tech., 16, 41–55,, 2023b. 

Faassen, K. A. P., Nguyen, L. N. T., Broekema, E. R., Kers, B. A. M., Mammarella, I., Vesala, T., Pickers, P. A., Manning, A. C., Vilà-Guerau de Arellano, J., Meijer, H. A. J., Peters, W., and Luijkx, I. T.: Diurnal variability of atmospheric O2, CO2, and their exchange ratio above a boreal forest in southern Finland, Atmos. Chem. Phys., 23, 851–876,, 2023. 

Farquhar, G. D., Caemmerer, S. V., and Berry, J. A.: A biochemical-model of photosynthetic CO2 assimilation in leaves of C3 species, Planta, 149, 78–90,, 1980. 

Goulden, M. L., Munger, J. W., Fan, S. M., Daube, B. C., and Wofsy, S. C.: Measurements of carbon sequestration by long-term eddy covariance: Methods and a critical evaluation of accuracy, Glob. Change Biol., 2, 169–182,, 1996. 

Grundner, F., Schober, R., and Schwappach, A. F.: Massentafeln zur Bestimmung des Holzgehaltes stehender Waldbäume und Waldbestände: nach den Arbeiten der deutschen und österreichischen forstlichen Versuchsanstalten, Parey, 1952. 

Hanson, P. J., Amthor, J. S., Wullschleger, S. D., Wilson, K., Grant, R. F., Hartley, A., Hui, D., Hunt, J., E Raymond, Johnson, D. W., and Kimball, J. S.: Oak forest carbon and water simulations: model intercomparisons and evaluations against independent data, Ecol. Monogr., 74, 443–489,, 2004. 

Helm, J., Hartmann, H., Göbel, M., Hilman, B., Herrera Ramírez, D., and Muhr, J.: Low-cost chamber design for simultaneous CO2 and O2 flux measurements between tree stems and the atmosphere, Tree Physiol., 41, 1767–1780,, 2021. 

Hicks Pries, C., Angert, A., Castanha, C., Hilman, B., and Torn, M. S.: Using respiration quotients to track changing sources of soil respiration seasonally and with experimental warming, Biogeosciences, 17, 3045–3055,, 2020. 

Hilman, B. and Angert, A.: Measuring the ratio of CO2 efflux to O2 influx in tree stem respiration, Tree Physiol., 36, 1422–1431,, 2016. 

Hilman, B., Muhr, J., Trumbore, S. E., Kunert, N., Carbone, M. S., Yuval, P., Wright, S. J., Moreno, G., Pérez-Priego, O., Migliavacca, M., Carrara, A., Grünzweig, J. M., Osem, Y., Weiner, T., and Angert, A.: Comparison of CO2 and O2 fluxes demonstrate retention of respired CO2 in tree stems from a range of tree species, Biogeosciences, 16, 177–191,, 2019. 

Hilman, B., Weiner, T., Haran, T., Masiello, C. A., Gao, X., and Angert, A.: The apparent respiratory quotient of soils and tree stems and the processes that control It, J. Geophys. Res.-Biogeo., 127, e2021JG006676,, 2022. 

Hockaday, W. C., Gallagher, M. E., Masiello, C. A., Baldock, J. A., Iversen, C. M., and Norby, R. J.: Forest soil carbon oxidation state and oxidative ratio responses to elevated CO2, J. Geophys. Res.-Biogeo., 120, 1797–1811,, 2015. 

Ishidoya, S., Morimoto, S., Aoki, S., Taguchi, S., Goto, D., Murayama, S., and Nakazawa, T.: Oceanic and terrestrial biospheric CO2 uptake estimated from atmospheric potential oxygen observed at Ny-Alesund, Svalbard, and Syowa, Antarctica, Tellus B, 64, 18924,, 2012. 

Ishidoya, S., Murayama, S., Takamura, C., Kondo, H., Saigusa, N., Goto, D., Morimoto, S., Aoki, N., Aoki, S., and Nakazawa, T.: O2 : CO2 exchange ratios observed in a cool temperate deciduous forest ecosystem of central Japan, Tellus B, 65, 21120,, 2013. 

Ishidoya, S., Murayama, S., Kondo, H., Saigusa, N., Kishimoto-Mo, A. W., and Yamamoto, S.: Observation of O2 : CO2 exchange ratio for net turbulent fluxes and its application to forest carbon cycles, Ecol. Res., 30, 225–234,, 2015. 

Ishidoya, S., Sugawara, H., Terao, Y., Kaneyasu, N., Aoki, N., Tsuboi, K., and Kondo, H.: O2 : CO2 exchange ratio for net turbulent flux observed in an urban area of Tokyo, Japan, and its application to an evaluation of anthropogenic CO2 emissions, Atmos. Chem. Phys., 20, 5293–5308,, 2020. 

Juergensen, J., Muhr, J., and Knohl, A.: Variations of the oxidative ratio across ecosystem components and seasons in a managed temperate beech forest (Leinefelde, Germany), Forests, 12, 1693,, 2021. 

Keeling, R. F. and Manning, A. C.: 5.15 – Studies of Recent Changes in Atmospheric O2 Content, in: Treatise on Geochemistry, 2nd Edn., edited by: Holland, H. D. and Turekian, K. K., Elsevier, Oxford, 385–404,, 2014. 

Keeling, R. F. and Shertz, S. R.: Seasonal and interannual variations in atmospheric oxygen and implications for the global carbon-cycle, Nature, 358, 723–727, 1992a. 

Keeling, R. F. and Shertz, S. R.: Seasonal and interannual variations in atmospheric oxygen and implications for the global carbon cycle, Nature, 358, 723–727,, 1992b. 

Knohl, A. and Baldocchi, D. D.: Effects of diffuse radiation on canopy gas exchange processes in a forest ecosystem, J. Geophys. Res.-Biogeo., 113, G02023,, 2008. 

Knohl, A. and Buchmann, N.: Partitioning the net CO2 flux of a deciduous forest into respiration and assimilation using stable carbon isotopes, Global Biogeochem. Cy., 19, GB4008,, 2005. 

Kowalski, A. S.: The boundary condition for vertical velocity and its interdependence with surface gas exchange, Atmos. Chem. Phys., 17, 8177–8187,, 2017. 

Kowalski, A. S., Serrano-Ortiz, P., Miranda-García, G., and Fratini, G.: Disentangling turbulent gas diffusion from non-diffusive transport in the boundary layer, Bound.-Lay. Meteorol., 179, 347–367,, 2021. 

Krogh, A.: The Composition of the Atmosphere: An Account of Preliminary Investigations and a Programme, Kongelige Danske Videnskabernes Selskab: Matematisk – Fysiske Meddelelser, Vol. 1, 12, A.F. Høst, 1919. 

Lenschow, D. H., Mann, J., and Kristensen, L.: How long is long enough when measuring fluxes and other turbulence statistics, J. Atmos. Ocean Tech., 11, 661–673,<0661:HLILEW>2.0.CO;2, 1994. 

Manning, A. C. and Keeling, R. F.: Global oceanic and land biotic carbon sinks from the Scripps atmospheric oxygen flask sampling network, Tellus B, 58, 95–116,, 2006. 

Meredith, L., Commane, R., Munger, J., Dunn, A., Tang, J., Wofsy, S., and Prinn, R.: Ecosystem fluxes of hydrogen: a comparison of flux-gradient methods, Atmos. Meas. Tech., 7, 2787–2805,, 2014. 

Ogee, J., Peylin, P., Cuntz, M., Bariac, T., Brunet, Y., Berbigier, P., Richard, P., and Ciais, P.: Partitioning net ecosystem carbon exchange into net assimilation and respiration with canopy-scale isotopic measurements: An error propagation analysis with 13CO2 and CO18O data, Global Biogeochem. Cy., 18, GB2019,, 2004. 

Oikawa, P. Y., Sturtevant, C., Knox, S. H., Verfaillie, J., Huang, Y. W., and Baldocchi, D. D.: Revisiting the partitioning of net ecosystem exchange of CO2 into photosynthesis and respiration with simultaneous flux measurements of 13CO2 and CO2, soil respiration and a biophysical model, CANVEG, Agr. Forest Meteorol., 234/235, 149–163,, 2017. 

Pickers, P. A., Manning, A. C., Sturges, W. T., Le Quéré, C., Mikaloff Fletcher, S. E., Wilson, P. A., and Etchells, A. J.: In situ measurements of atmospheric O2 and CO2 reveal an unexpected O2 signal over the tropical Atlantic Ocean, Global Biogeochem. Cy., 31, 1289–1305,, 2017. 

Randerson, J., Masiello, C., Still, C., Rahn, T., Poorter, H., and Field, C.: Is carbon within the global terrestrial biosphere becoming more oxidized? Implications for trends in atmospheric O2, Glob. Change Biol., 12, 260–271,, 2006. 

Raupach, M. R.: Applying Lagrangian fluid mechanics to infer scalar source distributions from concentration profiles in plant canopies, Agr. Forest Meteorol., 47, 85–108,, 1989. 

Rebmann, C., Aubinet, M., Schmid, H., Arriga, N., Aurela, M., Burba, G., Clement, R., De Ligne, A., Fratini, G., Gielen, B., Grace, J., Graf, A., Gross, P., Haapanala, S., Herbst, M., Hortnagl, L., Ibrom, A., Joly, L., Kljun, N., Kolle, O., Kowalski, A., Lindroth, A., Loustau, D., Mammarella, I., Mauder, M., Merbold, L., Metzger, S., Molder, M., Montagnani, L., Papale, D., Pavelka, M., Peichl, M., Roland, M., Serrano-Ortiz, P., Siebicke, L., Steinbrecher, R., Tuovinen, J. P., Vesala, T., Wohlfahrt, G., and Franz, D.: ICOS eddy covariance flux-station site setup: a review, Int. Agrophys., 32, 471–494,, 2018.  

Seibt, U., Brand, W. A., Heimann, M., Lloyd, J., Severinghaus, J. P., and Wingate, L.: Observations of O2 : CO2 exchange ratios during ecosystem gas exchange, Global Biogeochem. Cy., 18, GB4024,, 2004. 

Severinghaus, J.: Studies of the Terrestrial O2 and Carbon Cycles in Sand Dune Gases and in Biosphere 2, PhD thesis,, 1995. 

Stephens, B. B., Keeling, R. F., and Paplawsky, W. J.: Shipboard measurements of atmospheric oxygen using a vacuum-ultraviolet absorption technique, Tellus B, 55, 857–878,, 2003. 

Tamrakar, R., Rayment, M. B., Moyano, F., Mund, M., and Knohl, A.: Implications of structural diversity for seasonal and annual carbon dioxide fluxes in two temperate deciduous forests, Agr. Forest Meteorol., 263, 465–476,, 2018. 

Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, 57–80,, 2004. 

Tcherkez, G., Nogues, S., Bleton, J., Cornic, G., Badeck, F., and Ghashghaie, J.: Metabolic origin of carbon isotope composition of leaf dark-respired CO2 in French bean, Plant Physiol., 131, 237–244,, 2003. 

Tohjima, Y., Mukai, H., Machida, T., Hoshina, Y., and Nakaoka, S.-I.: Global carbon budgets estimated from atmospheric O2/ N2 and CO2 observations in the western Pacific region over a 15-year period, Atmos. Chem. Phys., 19, 9269–9285,, 2019. 

Van Oijen, M., Rougier, J., and Smith, R.: Bayesian calibration of process-based forest models: bridging the gap between models and data, Tree Physiol., 25, 915–927,, 2005. 

Wehr, R. and Saleska, S. R.: An improved isotopic method for partitioning net ecosystem-atmosphere CO2 exchange, Agr. Forest Meteorol., 214, 515–531,, 2015. 

Wilson, J.: Turbulent transport within the plant canopy, in: The Forest-Atmosphere Interaction, Springer, 443–480,, 1989. 

Worrall, F., Clay, G. D., Masiello, C. A., and Mynheer, G.: Estimating the oxidative ratio of the global terrestrial biosphere carbon, Biogeochemistry, 115, 23–32,, 2013. 

Yakir, D. and Wang, X. F.: Fluxes of CO2 and water between terrestrial vegetation and the atmosphere estimated from isotope measurements, Nature, 380, 515–517,, 1996. 

Zobitz, J. M., Burns, S. P., Ogee, J., Reichstein, M., and Bowling, R.: Partitioning net ecosystem exchange of CO2: A comparison of a Bayesian/isotope approach to environmental regression methods, J. Geophys. Res.-Biogeo., 112, G03013,, 2007. 

Short summary
A better understanding of O2 fluxes, their exchange ratios with CO2 and their interrelations with environmental conditions would provide further insights into biogeochemical ecosystem processes. We, therefore, used the multilayer canopy model CANVEG to simulate and analyze the flux exchange for our forest study site for 2012–2016. Based on these simulations, we further successfully tested the application of various micrometeorological methods and the prospects of real O2 flux measurements.
Final-revised paper