Technical note: New insights into stomatal oxygen transport viewed as a multicomponent diffusion process

We investigate oxygen (O₂) transport through stomata, focusing on its interaction with water vapour (H₂O) flux. The dominant H₂O flux exerts a drag force on other gases, a well-studied effect in the ternary air–water vapor–carbon dioxide (CO₂) system but unexplored for O₂ transport. This study aims to: (1) apply the Stefan-Maxwell equations to a quaternary system of H₂O, O₂, CO₂, and N₂; (2) identify conditions where O₂ transport from stomata to the atmosphere occurs against its mole fraction gradient (“uphill”); and (3) derive an expression linking the O₂ mole fraction in sub-stomatal air spaces (x_oi) to that in the atmosphere (x_oa) based on atmospheric relative humidity.

Our theoretical results, constrained by typical values from previous flux observations of this quaternary system, reveal distinct transport regimes defined by the mole flux ratio of H₂O and O₂ ($F_{\mathrm{w}}/F_{\mathrm{o}}$). Uphill O₂ diffusion occurs in the common regime where $F_{\mathrm{w}}/F_{\mathrm{o}}\gg \mathrm{1}$, and internal O₂ mole fraction increases towards its atmospheric value as relative humidity tends to 100 %. These theoretical results offer a framework for interpreting laboratory and field experiments on stomatal O₂ exchange under stagnant atmospheric or low Reynolds number conditions and can support the development of more physically accurate models of leaf–atmosphere oxygen exchange.

1 Introduction

The terrestrial biosphere exchange of water (H₂O), carbon dioxide (CO₂) and oxygen (O₂), plays a crucial role in the global water and carbon cycles. These cycles integrate processes across scales, from the stomatal level to the global atmosphere. Transport between stomata and the atmosphere represents the smallest-scale process governing the exchange of greenhouse gases. When scaled up to leaf and canopy levels, this stomata-atmosphere interaction modifies surface fluxes, which directly influence the diurnal variability of carbon tracers (Faassen et al., 2023) and impacts weather (Jacobs and de Bruin, 1997; Boussetta et al., 2013) and climate-carbon patterns (Bonan et al., 2024). Advancing our understanding of these interactions requires a more integrated approach that bridges ecophysiology and transport dynamics across multiple scales (Vilà-Guerau de Arellano et al., 2023; Miralles et al., 2025). In doing so, we improve the representation of carbon tracers across a wide range of scales.

While the exchange of CO₂ (Goudriaan, 1977; von Caemmerer and Farquhar, 1981) and H₂O (Jarvis and McNaughton, 1986), typically quantified through molar flux density, has been extensively studied and well-characterized, the transport of O₂ between stomata and the atmosphere has received comparatively little attention. Recent advances now allow biosphere-atmosphere exchange of O₂ to be measured with increasing accuracy (Stephens et al., 2007; Ishidoya et al., 2013; Battle et al., 2019; Faassen et al., 2023), supported by conceptual land-atmosphere modelling studies (Yan et al., 2023; Faassen et al., 2024). This provides a unique opportunity to study the H₂O-CO₂-O₂ system through a complete and integrated approach, advancing our understanding of processes related to photosynthesis and the transport of carbon tracers. Notably, O₂ transport exhibits distinct characteristics: O₂ exists at a much higher mole fraction in the atmosphere than either CO₂ or water vapour, yet its flux is of the same order of magnitude as that of CO₂ and typically two to three orders of magnitude smaller than the water vapour flux. As a result, the O₂ flux may be influenced by fluxes of other gases. The significance of this issue can strongly affect the interpretation of O₂ exchange on stomatal level, as was recently highlighted by Vesala (2024), whereas on the scale of canopy exchange the O₂-transport is dominated by turbulent fluxes (Ishidoya et al., 2013; Faassen et al., 2023).

The stomatal–atmosphere exchange is affected by molecular interaction between gases (also referred to as species). This exchange is primarily driven by the dominant effect of water (i.e. stomatal transpiration) relative to the individual fluxes of O₂ and CO₂. Specifically, transpiration at the stomatal level generates drag forces (also referred to as frictional forces) and a bulk displacement of air known as Stefan (or molar) flow. This process leads to gradients between the sub-stomatal cavity and the atmosphere mole fractions, thereby modifying their diffusive transport. Traditionally, the relationship between the flux of a given species and its mole fraction gradient is described by Fick’s law, but the additional drag forces associated with multicomponent interactions can lead to the emergence of more complex (non-Fickian) flux–gradient relationships. Accurately interpreting measurements of O₂, CO₂, and H₂O, as well as their representation in models, requires these species–species interactions to be taken into account, e.g. Jarman (1974). In this framework, the Stefan–Maxwell equations provide a more complete description than Fick’s law, unifying the effects of molar flow, interspecies drag (friction), and diffusive flux into a single formulation. Fick’s law can be deduced from the Stefan–Maxwell equations as an approximation.

To address this in the case of O₂, we investigate the transport of O₂, CO₂, H₂O, and the inert carrier nitrogen (N₂) between the leaf intercellular space within the stomata and the atmosphere. Our approach is based on the well-established theory of multicomponent diffusion (Jarman, 1974). Mathematically, the theory is encapsulated in the Stefan–Maxwell equations, which describe the net balance between the thermodynamic force (i.e. mole fraction gradient) on a given species that drives its diffusion and the momentum exchange (frictional drag) arising from binary collisions with all the other species, which depends on their relative molar velocities. Previous studies have shown that Stefan flow can play a significant role in CO₂ transport at the stomatal level (Jarman, 1974), and may also be relevant at larger atmospheric scales – from the canopy to the atmospheric boundary layer – as well as for key compounds involved in the photosynthesis cycle, such as O₂ (Kowalski et al., 2021). However, a comprehensive mathematical description of oxygen exchange in relation to water vapour, carbon dioxide, and nitrogen at the stomatal level remains absent.

We adopt a two-step approach. First, we analyze a binary mixture following the method described in Lushnikov et al. (1994) where one species diffuses in a stagnant carrier species. Secondly, we examine a multicomponent (quaternary) system consisting of H₂O, O₂, CO₂ diffusing in stagnant N₂, building on previous analyses of ternary systems, such as H₂O-CO₂-air (Jarman, 1974; von Caemmerer and Farquhar, 1981). Our key objectives are:

1. To derive interrelated molar flux density and gradient expressions for O₂, H₂O, and CO₂ that account for the net balance between drag and thermodynamic forces.

2. To identify general conditions under which the transport of O₂ occurs against its mole fraction gradient, referred here as “uphill”, due to the drag effect of transpiration.

3. To derive an expression for the internal stomatal O₂ mole fraction as a function of atmospheric relative humidity.

To achieve these objectives, we apply the Stefan-Maxwell equations to the system composed of N₂-O₂-H₂O-CO₂. We focus on the transport processes between stomata and the atmosphere, without explicitly accounting for the biochemistry of O₂ sources and sinks within the leaf. While our study is based on theoretical calculations, these are constrained by representative measurements of water vapour, carbon dioxide, and oxygen fluxes at the leaf level (Vilà-Guerau de Arellano et al., 2020), which are subsequently scaled down to the stomatal level at which the Stefan-Maxwell equations apply.

The rest of the paper is structured as follows. In Sect. 2 we outline the three key elements of multicomponent diffusion theory that were used in this study: molar balance, decomposition of species transport into molar flow and molecular diffusive flux components, and the Stefan-Maxwell equations relating species mole fraction gradients to species flux densities. In Sect. 3, as a simple reference case, we use this theory to derive expressions for gas transport in a binary system consisting of one tracer species diffusing in stagnant air (e.g. Jarman, 1974). In Sect. 4 we extend this theory to a more realistic quaternary system consisting of O₂, H₂O and CO₂ diffusing in stagnant N₂. In particular, we derive a simple expression linking the O₂ mole fraction gradient to relative humidity. We then identify the general regime conditions of O₂ transport with respect to the water vapour flux. We discuss the implications of these results in Sect. 5.

2 Multicomponent diffusion theory

To study the exchange of O₂ between a single stoma and a stagnant atmosphere relating to H₂O, CO₂ and N₂, we derive a theory that is applied to gas transport through a single stomatal pore. Our derivations are at the stomatal level, but can be scaled up to the leaf level. When discussing our results numerically using representative leaf-scale flux observations (Vilà-Guerau de Arellano et al., 2020), we convert from stomatal-scale flux densities (mol m⁻² pore s⁻¹) to leaf-scale flux densities (mol m⁻² leaf s⁻¹). In short, by multiplying the stomatal-scale flux density by the cross-sectional area of a pore (m² pore), we obtain the flux through a single stoma. Multiplying this value by the stomatal density (number of stomata per m² of leaf surface) then yields the flux density at the leaf scale.

Our key conclusions do not depend on this conversion factor, since this it cancels out when scaling (or downscaling) from the stomata to the atmosphere. For simplicity we consider one-dimensional gas transport through a stomatal pore connecting the stomatal interior to the external atmosphere. Appendix B describes the geometrical assumptions underlying this simplification. The theory presented in the next section can readily be generalized to three spatial dimensions. Finally, for clarity, Table A1 in Appendix A, provides symbol definitions, units, and representative parameter values, including typical values for mole fractions, prescribed molar flux densities based on observations, and molecular diffusivity coefficients.

2.1 Molar balance

Molar balance of a tracer provides the overall context for the description of multicomponent diffusion. The molar balance equation for c_α, the concentration of chemical species α (units mol m⁻³), reads

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial c_{\mathit{\alpha }}}{\partial t}}=-{\displaystyle \frac{\partial F_{\mathit{\alpha }}}{\partial l}}=R_{\mathit{\alpha }},\end{array}$$

where F_α (mol m⁻² s⁻¹) and R_α are, respectively, the stomatal-scale molar flux density and the net production rate (i.e.. sources – sinks) of species α, t is time and l is the distance along the pore (see Fig. B1). For n species, the total mole fraction of the mixture is

$$\begin{array}{}\text{(2)}& c_{\mathrm{t}}=\sum _{\mathit{\alpha }=\mathrm{1}}^n{c}_{\mathit{\alpha }}.\end{array}$$

and Eq. (1) can then be re-expressed in terms of the mole fraction $x_{\mathit{\alpha }}=c_{\mathit{\alpha }}/c_{\mathrm{t}}$ (t is here the total amount of species), as

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}{\displaystyle \frac{\partial F_{\mathit{\alpha }}}{\partial l}}+r_{\mathit{\alpha }},\end{array}$$

where $r_{\mathit{\alpha }}=R_{\mathit{\alpha }}/c_{\mathrm{t}}$.

2.2 Decomposition of fluxes into molar flow and diffusive flux

For each chemical species α, we can decompose its molar flux density F_α into a component related to the average flow of the mixture as a whole (molar flow, also called Stefan flow) and a diffusive component related to its transport relative to this average (diffusive flux). We first note that F_α is related to the species velocity u_α (relative to the stationary laboratory frame) by

$$\begin{array}{}\text{(4)}& F_{\mathit{\alpha }}=c_{\mathit{\alpha }}{u}_{\mathit{\alpha }}\end{array}$$

The mole-average velocity of the mixture (u) is defined as the sum of the species velocities weighted by their respective mole fractions:

$$\begin{array}{}\text{(5)}& u=\sum _{\mathit{\alpha }=\mathrm{1}}^n{x}_{\mathit{\alpha }}{u}_{\mathit{\alpha }},\end{array}$$

By substituting $x_{\mathit{\alpha }}=c_{\mathit{\alpha }}/c_{\mathrm{t}}$ into Eq. (5) and using Eq. (4), u can be expressed in terms of the total molar flux density of the mixture (F_t):

$$\begin{array}{}\text{(6)}& u={\displaystyle \frac{F_{\mathrm{t}}}{c_{\mathrm{t}}}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& F_{\mathrm{t}}=\sum _{\mathit{\alpha }=\mathrm{1}}^n{F}_{\mathit{\alpha }},\end{array}$$

We can then decompose F_α (Eq. 4) into the sum of a Stefan flow defined by

$$\begin{array}{}\text{(8)}& s_{\mathit{\alpha }}=c_{\mathit{\alpha }}u,\end{array}$$

and a diffusive flux defined by

$$\begin{array}{}\text{(9)}& f_{\mathit{\alpha }}=c_{\mathit{\alpha }}(u_{\mathit{\alpha }}-u).\end{array}$$

Using Eqs. (4) and (9), we can find an expression of the flux F_α related to the molar flow and the diffusive flux:

$$\begin{array}{}\text{(10)}& F_{\mathit{\alpha }}=c_{\mathit{\alpha }}{u}_{\mathit{\alpha }}=s_{\mathit{\alpha }}+f_{\mathit{\alpha }}.\end{array}$$

From these relationships we note that the total diffusive flux of the mixture is zero, i.e.. ${\sum }_{\mathit{\alpha }=\mathrm{1}}^n{f}_{\mathit{\alpha }}=\mathrm{0}$.

Whereas Eqs. (8) and (9) define Stefan flow and diffusion in terms of species velocities, we can use Eqs. (6) and (10) to re-express these components in terms of the species molar flux densities, giving

$$\begin{array}{}\text{(11)}& s_{\mathit{\alpha }}=x_{\mathit{\alpha }}{F}_{\mathrm{t}}\end{array}$$

and

$$\begin{array}{}\text{(12)}& f_{\mathit{\alpha }}=F_{\mathit{\alpha }}-x_{\mathit{\alpha }}{F}_{\mathrm{t}}.\end{array}$$

From Eq. (11) we note that the Stefan flow of a given species is proportional to its mole fraction so that, for a given total flux density F_t, we expect Stefan flow to be much larger for O₂ (x_o = 0.21) than for H₂O ($x_{\mathrm{w}}={\mathrm{10}}^{-\mathrm{2}}$) or CO₂ ($x_{\mathrm{c}}=\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-\mathrm{4}}$) (see Table A for additional values). Using the decomposition of Eqs. (8)–(10), we can then re-express molar balance (Eq. 3) as

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial t}}+{\displaystyle \frac{\partial \left(x_{\mathit{\alpha }}u\right)}{\partial l}}=-{\displaystyle \frac{\mathrm{1}}{c}}{\displaystyle \frac{\partial f_{\mathit{\alpha }}}{\partial l}}+r_{\mathit{\alpha }}.\end{array}$$

where the first and second terms on the right-hand side represent the contributions to molar transport from Stefan flow and diffusion, respectively.

To complete our analysis of the molar O₂-equation, and focusing on the molar balance for oxygen, Appendix D shows Eq. (13) in a dimensionless form. The analysis shows that the diffusive term can be ignored at atmospheric scales – due to the dominance of turbulence with respect the molecular diffusion transport. However, at the stomatal scale, characterized by low Reynolds numbers (rougly speaking maximum values of 10⁴, the Stefan flow and the diffusive terms needs to be included. It is at this stomatal scale, and in order to solve Eq. (13) for the mole fractions x_α, we need to relate the Stefan flow and diffusion terms to the species mole fractions. This linkage is supplied by the Stefan-Maxwell equations, and it is at this stage that the need for a multicomponent approach becomes apparent.

2.3 Stefan-Maxwell equations

Stefan-Maxwell equations provide a framework to study the dependence of the flux as a function of the gradients of other of all other species. For a multicomponent system consisting of n chemical species, the relationship between the one-dimensional mole fraction gradient of species α and the species flux densities is given by the Stefan-Maxwell equation (Curtiss and Hirschfelder, 1949; Kalkkinen et al., 1991; Bird et al., 2007):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial l}}& =\sum _{\mathit{\beta }=\mathrm{1},\mathit{\beta }\ne \mathit{\alpha }}^n{\displaystyle \frac{x_{\mathit{\alpha }}{F}_{\mathit{\beta }}-x_{\mathit{\beta }}{F}_{\mathit{\alpha }}}{c{D}_{\mathit{\alpha }\mathit{\beta }}}}=\sum _{\mathit{\beta }=\mathrm{1},\mathit{\beta }\ne \mathit{\alpha }}^n{\displaystyle \frac{x_{\mathit{\alpha }}{f}_{\mathit{\beta }}-x_{\mathit{\beta }}{f}_{\mathit{\alpha }}}{c{D}_{\mathit{\alpha }\mathit{\beta }}}}\\ & =\sum _{\mathit{\beta }=\mathrm{1},\mathit{\beta }\ne \mathit{\alpha }}^n{\displaystyle \frac{x_{\mathit{\alpha }}{x}_{\mathit{\beta }}(u_{\mathit{\beta }}-u_{\mathit{\alpha }})}{D_{\mathit{\alpha }\mathit{\beta }}}},\end{array}\end{array}$$

where D_αβ is the binary diffusivity of species α and β, and the summations are over all species β other than species α. The second equality follows (i) substituting F_α and F_β in the numerator using Eq. (12) and (ii) cancellation of the Stefan flow components (x_αF_t and x_βF_t). The last and third equality follows substituting F_α and F_β using Eq. (4).

Physically, Eq. (14) describes the net balance of forces on species α. The left-hand side is the thermodynamic driving force for diffusion of species α while the right-hand side is the sum of the drag (frictional) forces due to molecular collisions between species α and all the other species. The drag force on species α due to collisions with species β is proportional to the product of their mole fractions (x_αx_β) and their velocity difference (u_β−u_α). The inverse of the binary diffusivity ($\mathrm{1}/D_{\mathit{\alpha }\mathit{\beta }}$) then represents an effective drag coefficient. A useful feature of the Stefan-Maxwell equations is that the binary diffusivities are virtually independent of species mole fractions (Bird et al., 2007). As shown below (Sects. 3 and 4), Eq. (14) can be inverted to express the molar flux density of species α (F_α) in terms of the spatial gradients of all other species in the system. This facilitates the comparison with Fick's law. It also shows the implications of solving the molar balance equation (Eq. 13) using a multicomponent approach that accounts for the interactions between the different species.

3 Binary mixture: one species diffusing in stagnant air

We first apply Eq. (14) to a binary system (n=2) consisting of one tracer species (e.g. α= H₂O, O₂ or CO₂) diffusing in stagnant air. Here we define air composed solely by nitrogen and indicated by the subscript a. This provides a simplified but useful reference case against which the more realistic quaternary system will be compared (Sect. 4). In this case there is only one term on the right-hand side of Eq. (14 (β=α). Substituting F_a=0 (stagnant air, i.e.. u_a=0, where a means air) and $x_a=\mathrm{1}-x_{\mathit{\alpha }}$, then solving for F_α gives:

$$\begin{array}{}\text{(15)}& F_{\mathit{\alpha }}=-{\displaystyle \frac{c_{\mathrm{t}}{D}_{\mathit{\alpha }a}}{\mathrm{1}-x_{\mathit{\alpha }}}}\left({\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial l}}\right),\end{array}$$

where D_αa is the binary diffusivity of species α in air. The total molar flux density of the binary mixture is simply $F_{\mathrm{t}}=F_{\mathit{\alpha }}+F_{\mathrm{a}}=F_{\mathit{\alpha }}$ (since F_a=0), so from Eqs. (11), (12) and (15) the Stefan flow and diffusive components of F_α are given by

$$\begin{array}{}\text{(16)}& s_{\mathit{\alpha }}=x_{\mathit{\alpha }}{F}_{\mathit{\alpha }}=-{\displaystyle \frac{x_{\mathit{\alpha }}}{\mathrm{1}-x_{\mathit{\alpha }}}}{c}_{\mathrm{t}}{D}_{\mathit{\alpha }a}\left({\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial l}}\right)\end{array}$$

and

$$\begin{array}{}\text{(17)}& f_{\mathit{\alpha }}=(\mathrm{1}-x_{\mathit{\alpha }})F_{\mathit{\alpha }}=-c_{\mathrm{t}}{D}_{\mathit{\alpha }a}\left({\displaystyle \frac{\partial x_{\mathit{\alpha }}}{\partial l}}\right)\end{array}$$

Eq. (15) shows that for this simple binary system, the molar flux density of the diffusing species (F_α) is proportional to its mole fraction gradient $\left(\frac{\partial x_{\mathit{\alpha }}}{\partial l}\right)$, i.e. a Fickian diffusion law, with an effective diffusion coefficient $D_{\mathit{\alpha }a}/(\mathrm{1}-x_{\mathit{\alpha }})$ that depends on the mole fraction x_α. In Sect. 4, we contrast this result with the four-component system (N₂-O₂-H₂O-CO₂) for which, in general, the molar flux densities of water and oxygen (F_w and F_o) each depend on both mole fraction gradients $\left(\left(\frac{\partial x_{\mathrm{w}}}{\partial l}\right)\mathrm{and}\left(\frac{\partial x_{\mathrm{o}}}{\partial l}\right)\right)$, leading to non-Fickian diffusion. In Sect. 4 we will also show that, for water vapour (α=w), the Fickian diffusion result of Eq. (15) is a special case of the expression for F_w in the four-component system, in the limiting case that the water vapour flux density is much greater than that for oxygen or carbon dioxide (i.e.. F_w≫F_o or F_c), which is a realistic assumption. For oxygen (α=o), however, this is not the case (Sect. 4).

As anticipated from Eq. (11), when x_α is much smaller than 1 (i.e. for a dilute tracer), the contribution of Stefan flow (Eq. 16) becomes negligible, and F_α is dominated by the diffusive component f_α (Eq. 17). Therefore, the dependence of the effective diffusion coefficient on x_α can be interpreted as a correction to Eq. (17) arising from the influence of Stefan flow. This correction is significant for O₂. For CO₂, the correction is smaller but still non-negligible, as demonstrated in the next section, where we apply the full derivation of the Stefan–Maxwell equations (Eq. 14) to the four-component system under study.

4 Quaternary mixture: water vapour, oxygen and carbon dioxide diffusing in stagnant nitrogen

We now apply the Stefan–Maxwell equations (Eq. 14) to a four-component mixture consisting of water vapour, oxygen, and carbon dioxide diffusing in stagnant nitrogen (labelled by species indices $\mathit{\alpha }=w,o,c$, and n, respectively). Figure 1 shows the flux directions of the four components (Fig. 1a), their orders of magnitude, and the decomposition of the O₂-flux into a Stefan flow and a diffusive molar flux component (Fig. 1b). If we consider the flux components in Fig. 1b from the perspective of the physiological constraints, f_o depends on the photosynthesis rate and stomatal conductance. s_o depends on the transpiration rate, that is on water vapour pressure deficit and on stomatal conductance. If the conductance decreases, both f_o and s_o decrease, too.

Figure 1 Multicomponent fluxes in a system composed of water vapour (w), oxygen (o), and carbon dioxide (c) diffusing in stagnant nitrogen (n). The mole fraction at the sub-stomatal cavity and the atmosphere are expressed by the subscripts i and a, respectively. (a) Direction of the fluxes of the four compounds as investigated in our study, with the dominant role by the water vapor molar flux. (b) Decomposition of O₂-flux into a Stefan flow component and a diffusive molar flux component. The dashed lines indicate that the Stefan flow and diffusive flux are much larger than the O₂ flux as shown by the theory. L is the length and r the radius of the cylindrical symmetric stomatal pore, respectively. The O₂ flux is prescribed based on observations of the CO₂ flux at leaf level. More details on the geometry can be found at Fig. B1.

In Eq. (14), for each species α there are now three terms in the sum over all other species β≠α on the right-hand side. For H₂O, O₂ and CO₂ ($\mathit{\alpha }=w,o,c$) we obtain the following expressions for their mole fraction gradients in terms of the stomatal-scale molar flux densities:

$$\begin{array}{}\text{(18a)}& {\displaystyle \frac{\partial x_{\mathrm{w}}}{\partial l}}={\displaystyle \frac{x_{\mathrm{w}}{F}_{\mathrm{o}}-x_{\mathrm{o}}{F}_{\mathrm{w}}}{c_t{D}_{\mathrm{ow}}}}+{\displaystyle \frac{x_{\mathrm{w}}{F}_{\mathrm{n}}-x_n{F}_{\mathrm{w}}}{c_t{D}_{\mathrm{wn}}}}+{\displaystyle \frac{x_{\mathrm{w}}{F}_{\mathrm{c}}-x_{\mathrm{c}}{F}_{\mathrm{w}}}{c_t{D}_{\mathrm{wc}}}}\text{(18b)}& {\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial l}}={\displaystyle \frac{x_{\mathrm{o}}{F}_{\mathrm{w}}-x_{\mathrm{w}}{F}_{\mathrm{o}}}{c_t{D}_{\mathrm{ow}}}}+{\displaystyle \frac{x_{\mathrm{o}}{F}_{\mathrm{n}}-x_n{F}_{\mathrm{o}}}{c_t{D}_{\mathrm{on}}}}+{\displaystyle \frac{x_{\mathrm{o}}{F}_{\mathrm{c}}-x_{\mathrm{c}}{F}_{\mathrm{o}}}{c_t{D}_{\mathrm{oc}}}}\text{(18c)}& {\displaystyle \frac{\partial x_{\mathrm{c}}}{\partial l}}={\displaystyle \frac{x_{\mathrm{c}}{F}_{\mathrm{o}}-x_{\mathrm{o}}{F}_{\mathrm{c}}}{c_t{D}_{\mathrm{oc}}}}+{\displaystyle \frac{x_{\mathrm{c}}{F}_{\mathrm{n}}-x_n{F}_{\mathrm{c}}}{c_t{D}_{\mathrm{cn}}}}+{\displaystyle \frac{x_{\mathrm{c}}{F}_{\mathrm{w}}-x_{\mathrm{w}}{F}_{\mathrm{c}}}{c_t{D}_{\mathrm{wc}}}}.\end{array}$$

4.1 Non-Fickian relationship between H₂O and O₂ gradients and fluxes

The derived equations relating gradients to multiple fluxes (Eqs. 18a–18c) simplify if we assume that (1) on stoichiometric grounds the CO₂ molar flux density is equal and opposite to the O₂ molar flux density ($F_{\mathrm{c}}=-F_{\mathrm{o}}$), and (2) N₂ is stagnant (F_n=0). Focusing now on H₂O and O₂, and making the substitutions $F_{\mathrm{c}}=-F_{\mathrm{o}}$, F_n=0 and $x_{\mathrm{n}}=\mathrm{1}-x_{\mathrm{o}}-x_{\mathrm{w}}-x_{\mathrm{c}}$, Eqs. (18a)–(18c) can be written in the form:

$$\begin{array}{}\text{(19a)}& {\displaystyle \frac{\partial x_{\mathrm{w}}}{\partial l}}={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}(A{F}_{\mathrm{w}}+B{F}_{\mathrm{o}})\text{(19b)}& {\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial l}}={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}(C{F}_{\mathrm{w}}+D{F}_{\mathrm{o}}),\end{array}$$

where

$$\begin{array}{}\text{(20a)}& {\displaystyle }A=-{\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{ow}}}}-{\displaystyle \frac{\mathrm{1}-x_{\mathrm{o}}-x_{\mathrm{w}}-x_{\mathrm{c}}}{D_{\mathrm{wn}}}}-{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{wc}}}}\text{(20b)}& {\displaystyle }B=x_{\mathrm{w}}\left({\displaystyle \frac{\mathrm{1}}{D_{\mathrm{ow}}}}-{\displaystyle \frac{\mathrm{1}}{D_{\mathrm{wc}}}}\right)\text{(20c)}& {\displaystyle }C={\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{ow}}}}\text{(20d)}& {\displaystyle }D=-{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{ow}}}}-{\displaystyle \frac{\mathrm{1}-x_{\mathrm{o}}-x_{\mathrm{w}}-x_{\mathrm{c}}}{D_{\mathrm{on}}}}-{\displaystyle \frac{x_{\mathrm{o}}+x_{\mathrm{c}}}{D_{\mathrm{oc}}}}.\end{array}$$

The right-hand sides of Eqs. (19a)–(19b) only depend on F_w and F_o. Eqs. (19a)–(19b) can be inverted (e.g. using the standard rules of matrix algebra) to express the fluxes in terms of the gradients:

$$\begin{array}{}\text{(21a)}& {\displaystyle }{F}_{\mathrm{w}}={\displaystyle \frac{c_t}{AD-BC}}\left\{D{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{w}}}{\mathrm{\Delta }l}}-B{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{o}}}{\mathrm{\Delta }l}}\right\}\text{(21b)}& {\displaystyle }{F}_{\mathrm{o}}={\displaystyle \frac{c_t}{AD-BC}}\left\{-C{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{w}}}{\mathrm{\Delta }l}}+A{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{o}}}{\mathrm{\Delta }l}}\right\}.\end{array}$$

Therefore, in general, each of the H₂O and O₂ fluxes F_w and F_o depends on both gradients – and vice-versa, as seen from Eqs. (19a)–(19b) – giving a non-Fickian diffusion law between fluxes and gradients.

This is a general characteristics of the Stefan-Maxwell equations (Eq. 14), whether expressed in terms of F_α or f_α. The physical interpretation is that while $\partial x_{\mathit{\alpha }}/\partial l$ represents the thermodynamic force driving the diffusion of species α in the absence of other species, within a mixture this force is counterbalanced by the total drag force exerted on species α by all other species. This drag force depends on the relative velocities – and consequently on the flux densities – of all species in the mixture. The Stefan-Maxwell equations, whether expressed in terms of F_α or f_α, describe this balance of forces. As a result, when solving for F_α or f_α by inverting the Stefan-Maxwell equations, F_α or f_α depend on the mole fraction gradients of all species, not just on $\partial x_{\mathit{\alpha }}/\partial l$.

However, in the present case, Eqs. (21a)–(21b) simplify further if we note that typically F_w is much larger than F_o, so that in the right-hand sides of Eqs. (19a)–(19b) the terms involving F_o are much smaller than the terms involving F_w, and to a first approximation can therefore be neglected, leading to

$$\begin{array}{}\text{(22a)}& {\displaystyle \frac{\partial x_{\mathrm{w}}}{\partial l}}\approx {\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}A{F}_{\mathrm{w}}\approx {\displaystyle \frac{\mathrm{1}-x_{\mathrm{w}}}{c_{\mathrm{t}}{D}_{\mathrm{wa}}}}{F}_{\mathrm{w}}\text{(22b)}& {\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial l}}\approx {\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}C{F}_{\mathrm{w}}\approx {\displaystyle \frac{x_{\mathrm{o}}}{c_{\mathrm{t}}{D}_{\mathrm{wa}}}}{F}_{\mathrm{w}},\end{array}$$

where we have substituted the terms A and C using Eqs. (20a) and (20b). For the second equalities in Eqs. (22a)–(22b) we have made the further approximations in Eq. (20a): (1) x_c is typically much smaller than x_o and x_w, and (2) $D_{\mathrm{wn}}\approx D_{\mathrm{ow}}\approx D_{\mathrm{wa}}$, where D_wa is the molecular diffusivity of air (Table A1).

Equation (22a) reproduces the same Fickian relationship between the H₂O flux and H₂O gradient that was obtained previously for the binary case of H₂O diffusing in stagnant air (Eq. 15 with α=w). In contrast, Eq. (22b) implies that the O₂ gradient is proportional to the H₂O flux, and not the O₂ flux.

4.2 Diffusive and molar flows

As shown by Fig. 1b is interesting to calculate the contribution of the Stefan flow and diffusive molar flux to the flux. Using Eqs. (7), (11) and (12) we can quantify the Stefan flow and diffusive components of the molar flux densities for a given species α. In the quaternary system, since $F_{\mathrm{c}}+F_{\mathrm{o}}=F_{\mathrm{n}}=\mathrm{0}$ – recall that we assume $F_{\mathrm{c}}=-F_{\mathrm{o}}$ –, the total molar flux density of the mixture (F_t) (Eq. 7) is therefore equal to the water flux density, i.e. F_t=F_w. From Eqs. (11) and (12), therefore:

$$\begin{array}{}\text{(23)}& s_{\mathit{\alpha }}=x_{\mathit{\alpha }}{F}_{\mathrm{w}}\end{array}$$

and

$$\begin{array}{}\text{(24)}& f_{\mathit{\alpha }}=F_{\mathit{\alpha }}-x_{\mathit{\alpha }}{F}_{\mathrm{w}}\end{array}$$

These relationships between the stomatal-scale flux densities s_α, f_α and F_α (mol m⁻² pore s⁻¹) also apply to the leaf-scale flux densities. Here we use the following conversions as follows simply by multiplying both sides of Eqs. (23)–(24) by a_sρ_s: $s_{\mathit{\alpha }}^{\mathrm{leaf}}=a_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}{s}_{\mathit{\alpha }}$, $f_{\mathit{\alpha }}^{\mathrm{leaf}}=a_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}{f}_{\mathit{\alpha }}$ and $F_{\mathit{\alpha }}^{\mathrm{leaf}}=a_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}{F}_{\mathit{\alpha }}$ (mol m⁻² leaf s⁻¹), where a_s (m² pore) is the cross-sectional area of the stomatal aperture and ρ_s (number of stomata m⁻² leaf) is stomatal leaf density.

This allows us to gauge the relative importance of the Stefan and diffusive components at leaf scales using representative values of $F_{\mathit{\alpha }}^{\mathrm{leaf}}$. Table 1 gives values of $s_{\mathit{\alpha }}^{\mathrm{leaf}}$ and $f_{\mathit{\alpha }}^{\mathrm{leaf}}$ for H₂O, O₂ and CO₂. Also shown are the corresponding scaled gradients $a_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}\partial x_{\mathit{\alpha }}/\partial l$ calculated by multiplying both sides of Eqs. (18a)–(18c) by a_sρ_s and substituting $F_{\mathit{\alpha }}^{\mathrm{leaf}}$ into the right-hand sides of Eqs. (18a)–(18c).

Table 1 Stefan flow ($s_{\mathit{\alpha }}^{\mathrm{leaf}}$) and diffusive molar flux density ($f_{\mathit{\alpha }}^{\mathrm{leaf}}$) components of the leaf-scale H₂O, O₂ and CO₂ molar flux densities ($F_{\mathit{\alpha }}^{\mathrm{leaf}}$), calculated by multiplying both sides of Eqs. (23) and (24) by a_sρ_s and substituting typical values of $F_{\mathit{\alpha }}^{\mathrm{leaf}}=a_s{\mathit{\rho }}_s{F}_{\mathit{\alpha }}$, where a_s (m² pore) is stomatal aperture, ρ_s (number of stomata m⁻² leaf) is stomatal leaf density and F_α is the stomatal-scale flux density. Also shown are the scaled mole fraction gradients $a_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{s}}\partial x_{\mathit{\alpha }}/\partial l$ calculated from Eqs. (18a)–(18c) (multiplied through by a_sρ_s) using typical values of $F_{\mathit{\alpha }}^{\mathrm{leaf}}$ and the parameter values in Table A1. We use the following mole fractions for H₂O = 1 $\times {\mathrm{10}}^{-\mathrm{2}}$, O₂ = 0.21 and CO₂ = 4$\times {\mathrm{10}}^{-\mathrm{4}}$. For this value of H₂O, it corresponds to a specific humidity equal to 16.1 g_w kg${}_{\mathrm{air}}^{-\mathrm{1}}$ (representative high atmospheric humid conditions like the tropical forest).

As shown in Eq. (22a) the gradient of H₂O is also dependent of the mole fraction of water. Under very dry conditions, for instance air masses characterized by mole fractions H₂O = 10⁻³ equal to the specific humidity of 1.6 g_w kg${}_{\mathrm{air}}^{-\mathrm{1}}$, the scalar gradient is −9.46 m⁻¹ compared to the control value of −9.37⁻¹ m⁻¹ as shown in Table 1. These results stress the already-noted importance of Stefan flow for O₂ transport relative to Stefan flow for H₂O and CO₂. For H₂O, characterized by a positive value (i.e. stomata to atmosphere) Stefan flow (100 mol m⁻² leaf s⁻¹) is dominated by an even larger diffusive flow in the same direction (9900 mol m⁻² leaf s⁻¹). For CO₂, the positive Stefan flow (4 mol m⁻² leaf s⁻¹) is offset by diffusion in the opposite direction (−14 mol m⁻² leaf s⁻¹), but yet relevant. For O₂ the diffusive and Stefan components are individually two orders of magnitude larger than the prescribed F_o but are almost equal in magnitude and opposite in direction to each other.

Another consequence of the relative importance of Stefan flow in O₂ transport is that, unlike H₂O and CO₂, the flux of O₂ does not necessarily follow the opposite direction of its mole fraction gradient. Specifically, both F_w and F_c are opposite in sign to their respective mole fraction gradients, i.e. following qualitatively the Fickian law relationship between the flux and the gradient. This means that water vapour and CO₂ are transported “downhill”, from regions of high to low mole fractions, although deviating from Fick's law due to the influence of the fluxes of other species.

This is not the case for O₂. As shown in Eq. (22b), the O₂ gradient is directly proportional to, and in the same direction as, the O₂ Stefan flow rate, s_o = x_oF_w. This represents a clear departure from classical Fickian diffusion, where flux and gradient are expected to have opposite direction. We explore the conditions for uphill O₂ transport more generally in Sect. 4.4.

4.3 Relationship between internal O₂ mole fraction and relative humidity

From the gradient-flux relationships given by Eqs. (22a)–(22b) (which are based on the realistic approximations F_w≫F_o, x_c≪x_o and x_w, and $D_{\mathrm{wn}}\approx D_{\mathrm{ow}}\approx D_{\mathrm{wa}})$, it follows that the gradients in O₂ and H₂O satisfy the relationship

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathrm{1}}{x_{\mathrm{o}}}}{\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial l}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-x_{\mathrm{w}}}}{\displaystyle \frac{\partial x_{\mathrm{w}}}{\partial l}}.\end{array}$$

By integrating Eq. (25 from l=0 (internal stomata) to l=L (atmosphere, see Fig. B1), it follows that

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{x_{\mathrm{oi}}}{x_{\mathrm{oa}}}}\approx {\displaystyle \frac{\mathrm{1}-x_{\mathrm{wi}}}{\mathrm{1}-x_{\mathrm{wa}}}}\approx {\displaystyle \frac{\mathrm{1}-\mathit{\eta }{q}_{\mathrm{s}}}{\mathrm{1}-\mathit{\eta }{q}_{\mathrm{s}}\mathrm{RH}}},\end{array}$$

where x_oi and x_oa are, respectively, the internal and atmospheric O₂ mole fractions, x_wi and x_wa being the corresponding H₂O mole fractions. In the second equality we have expressed the latter in terms of the atmospheric saturated specific humidity (q_s), the atmospheric relative humidity RH (assuming RH = 1 internally), and the air-to-water molecular weight ratio (η=1.61 kg_a kg${}_{\mathrm{w}}^{-\mathrm{1}}$) to convert from mole fractions to relative humidity in expression (26).

From the expression (26) we infer that from normal atmospheric values of relative humidity (RH < 1), it follows from Eq. (26) that the internal O₂ mole fraction is smaller than its atmospheric value. Figure 2 shows the predicted ratio of internal to atmospheric O₂ mole fractions as a function of relative humidity. For all the values of RH, the mole fraction in the sub-stomatal cavity is lower than the one of the atmosphere. The lowest internal O₂ mole fraction occurs when RH = 0 (largest water demand by the atmosphere), when the predicted internal O₂ fraction is approximately 3 % lower than its atmospheric value. The predicted internal and atmospheric O₂ mole fraction equalise at RH = 1. However, when RH = 1 (lowest water demand) the assumption F_w≫F_o that underlies Eq. (26) is no longer valid. As a result, the simple relationship of Eq. (26) breaks down, and Eq. (22b) must be replaced by the exact expression given by Eq. (19b), i.e. including the term DF_o; since the coefficient D is negative (Eq. 20d). The effect of including this term is to decrease the O₂ gradient (i.e. increase the ratio $x_{\mathrm{oi}}/x_{\mathrm{oa}}$) so that the internal and atmospheric O₂ mole fractions equalize at a value of RH slightly less than 1 (data not shown).

Figure 2 The ratio of the internal oxygen to the atmospheric oxygen mole fraction as a function of the relative humidity calculated using Eq. (26), and a saturation specific humidity of $q_{\mathrm{s}}=\mathrm{19.6}\cdot {\mathrm{10}}^{-\mathrm{3}}$ kg_w kg${}_{\mathrm{a}}^{-\mathrm{1}}$ calculates at the temperature 298.15 K.

4.4 Uphill vs downhill O₂ transport: three regimes

In the previous sections, we have discussed the dominant role of the H₂O flux with respect the O₂transport. Here, we examine more generally the conditions under which O₂ transport (F_o) and/or its diffusive component (f_o) occurs “uphill”, i.e. against the O₂ mole fraction gradient, due to the dominant drag-force effect of water-driven Stefan flow. For F_o>0 (i.e. O₂ transport in the direction stoma to atmosphere), uphill O₂ transport occurs when $\partial x_{\mathrm{o}}/\partial l>\mathrm{0}$. From the exact expression for the O₂ gradient given by Eq. (19b), this occurs when $C{F}_{\mathrm{w}}+D{F}_{\mathrm{o}}>\mathrm{0}$. The expression reads:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\displaystyle \frac{F_{\mathrm{w}}}{F_{\mathrm{o}}}}& >{\displaystyle \frac{-D}{C}}={\displaystyle \frac{x_{\mathrm{w}}}{x_{\mathrm{o}}}}+{\displaystyle \frac{D_{\mathrm{ow}}}{D_{\mathrm{on}}}}{\displaystyle \frac{\mathrm{1}-x_{\mathrm{o}}-x_{\mathrm{w}}-x_{\mathrm{c}}}{x_{\mathrm{o}}}}\\ & +{\displaystyle \frac{D_{\mathrm{ow}}}{D_{\mathrm{oc}}}}{\displaystyle \frac{x_{\mathrm{o}}+x_{\mathrm{c}}}{x_{\mathrm{o}}}}={\displaystyle \frac{\mathrm{1}}{x_{\mathrm{o}}}}(\mathrm{1}+R)\end{array}\end{array}$$

where in the last expression R is defined by

$$\begin{array}{}\text{(28)}& \begin{array}{rl}R& =\left({\displaystyle \frac{D_{\mathrm{ow}}}{D_{\mathrm{on}}}}-\mathrm{1}\right)\left(\mathrm{1}-x_{\mathrm{o}}-x_{\mathrm{w}}-x_{\mathrm{c}}\right)\\ & +\left({\displaystyle \frac{D_{\mathrm{ow}}}{D_{\mathrm{oc}}}}-\mathrm{1}\right)\left(x_{\mathrm{o}}+x_{\mathrm{c}}\right).\end{array}\end{array}$$

Since D_ow>D_on and D_ow>D_oc (Table 1) it follows that R>0. With the assumed values of F_w, F_o and parameter values in Table 1, the condition Eq. (27) is comfortably satisfied, the left- and right-hand sides being 100 and ≈ 5.5 respectively. This condition will be violated when transpiration is sufficiently small compared to O₂ transport (e.g. when RH is close to 1).

From Eq. (24) with α=0, the diffusive O₂ flux density is negative (f_o<0) whenever

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{F_{\mathrm{w}}}{F_{\mathrm{o}}}}>{\displaystyle \frac{\mathrm{1}}{x_{\mathrm{o}}}}.\end{array}$$

Thus, as the ratio of water to oxygen flux densities ($F_{\mathrm{w}}/F_{\mathrm{o}}$) increases, we can identify three oxygen transport regimes (Table 2), distinguished by the relative orientation of the O₂ gradient ($\partial x_{\mathrm{o}}/\partial l$), O₂ flux density (F_o) and diffusive O₂ flux density (f_o).

Table 2 Three O₂ transport regimes distinguished by the relative orientations of the O₂ gradient ($\partial x_{\mathrm{o}}/\partial l$), O₂ flux density (F_o, assumed positive), and diffusive O₂ flux density (f_o), as determined by the value of the flux ratio ($F_{\mathrm{w}}/F_{\mathrm{o}}$) according to Eqs. (27)–(28) and (29). Positive values of $\partial x_{\mathrm{o}}/\partial l$ indicate that the internal O₂ mole fraction is less than atmospheric O₂ mole fraction. Positive f_o or F_o show that flux direction is from stoma to the atmosphere.

Table 2 outlines three regimes for water vapour and oxygen transport, defined by the direction of their respective mole fraction gradients. These regimes are characterised as uphill – where O₂ is transported from the stomata to the atmosphere despite an opposing mole fraction gradient – and downhill, where the flux and gradient are in opposite directions, consistent with classical Fickian diffusion. As noted above, due to the typically large value of the flux ratio ($F_{\mathrm{w}}/F_{\mathrm{o}}$), and since $(\mathrm{1}+R)/x_{\mathrm{o}}\approx \mathrm{5.5}$, the O₂ transport between the stoma and the atmosphere typically falls within regime 3, resulting in uphill transport of O₂ (F_o) in the direction from low to high O₂ mole fraction, as already seen in Table 2. In this regime, the molar flux density of H₂O drags O₂ molecules from the stoma to the atmosphere, driven by the O₂ Stefan flow s_o=x_oF_w. This drag effect leads to lower O₂ mole fraction in the stoma than in the atmosphere (x_oi<x_oa), and a negative diffusive O₂ molar flux density (f_o) (i.e. from the atmosphere to the stomata) that counterbalances but does not completely offset positive Stefan flow (Table 2) since we prescribe that the O₂ flux is from the stoma into the atmosphere. As $F_{\mathrm{w}}/F_{\mathrm{o}}$ decreases into regime 2, the O₂ gradient changes sign but the diffusive flux f_o remains negative (i.e. Stefan flow is still important), resulting in uphill diffusion from low to high O₂ mole fractions; this reflects the non-Fickian character of the flux-gradient relationships in general. At very low values of $F_{\mathrm{w}}/F_{\mathrm{o}}$ (regime 1), Stefan flow becomes very small, f_o changes sign, and O₂ transport is dominated by positive diffusion from the stoma to the atmosphere. It is important to note that the characterization of these three regimes does not account for the roles of photorespiration or dark respiration, and potential sources of O₂ within the sub-stomatal cavity are omitted from the analysis. Similarly, transport processes driven by thermal differences between the stomata and the atmosphere (thermodiffusion) are omitted.

5 Discussion

Our research extends the studies of Jarman (1974) and von Caemmerer and Farquhar (1981), which examined the transport of a ternary gas mixture comprising air, water vapour, and carbon dioxide. Here, air is composed by combining O₂ and N₂. Jarman (1974) also briefly considered the CO₂-H₂O-O₂-N₂ quaternary system considered in our study, but only with the aim of verifying the insignificant effect of O₂ transport on internal CO₂ mole fraction. For completeness, in Appendix D we demonstrate that our equations for the quaternary system (Eqs. 18a–18c) reduce to these ternary formulations.

Following Jarman (1974), in Appendix E we evaluated the percentage corrections to the O₂ gradient due to Stefan molar flow, using a ternary system (CO₂, H₂ and air, the latter compost by O₂ and N₂) as the reference system (cf. Sect. 3). As shown in Table E1, our results are consistent with those in Table 1 and confirm the key role of the H₂O molar flux density in determining the O₂ gradient.

Although previous studies have recognized that transport of O₂ is distinguished by its higher mole fraction and that Stefan flow could be relevant across scales from the canopy to the atmospheric boundary layer (Kowalski et al., 2021), our analysis of the fundamental molar balance equation demonstrates that the primary effects of Stefan flow emerge at the stomatal scale. To this end Appendix C shows a non-dimensional analysis of the governing equation of O₂ calculating the O₂-transport under representative Reynolds numbers measured in a greenhouse and at the top of a rainforest canopy. Connecting to previous studies of O₂ transport, Kowalski (2025) has proposed a simple model for the effects of water vapour flow on O₂ transport which assumes, a priori and without rigorous physical foundation, that the sub-stomatal O₂ pressure deficit (relative to the atmosphere) is equal to the atmospheric O₂ mole fraction multiplied by the water vapour pressure surplus in the sub-stomatal cavity. Since partial pressure is proportional to mole fraction, this assumption implies the relationship.

$$\begin{array}{}\text{(30)}& x_{\mathrm{oa}}-x_{\mathrm{oi}}=x_{\mathrm{oa}}(x_{\mathrm{wi}}-x_{\mathrm{wa}}),\end{array}$$

Interestingly, Eq. (26) of our study is equivalent to the relationship

$$\begin{array}{}\text{(31)}& x_{\mathrm{oa}}-x_{\mathrm{oi}}={\displaystyle \frac{x_{\mathrm{oa}}}{\mathrm{1}-x_{\mathrm{wa}}}}(x_{\mathrm{wi}}-x_{\mathrm{wa}}),\end{array}$$

which, because x_wa≪1, is practically identical to Eq. (30). However, in contrast to Eq. (30), our result (Eq. 31) has been derived from the fundamental physics of multicomponent diffusion (Stefan-Maxwell equations), and is itself an approximation (which assumes in particular that F_w≫F_o). Our underlying theory not only provides a physical justification for Eq. (29), but also a rigorous framework for calculating corrections to it under a wider range of conditions.

Our study does not yet address the physiological connection between sources and sinks of O₂ and CO₂ driven by photosynthesis, photorespiration and dark respiration. A logical next step would be to couple our theoretical expressions for the non-Fickian relationship between gradients and fluxes to biochemical models of photosynthetic CO₂ fixation (e.g. Farquhar, 1982) and O₂ transport dynamics during photosynthesis and photorespiration. Traditional models often over-simplify photorespiration and dark respiration by neglecting the dynamic interactions between different gases, key processes in the transport stomata-atmosphere.

Our theoretical findings – particularly the influence of H₂O molar flux density on internal O₂ mole fraction – suggest that incorporating a more physically realistic representation of O₂ dynamics could improve the coupling with models of photosynthetic processes, leading to more accurate estimates of O₂ transport's impact on photosynthetic efficiency. A better understanding of these processes could refine the predictions of photosynthesis under varying environmental conditions (Zhang et al., 2024), where photorespiration critically limits carbon assimilation. In a similar vein, it will be convenient to study the relationship between the O₂ flux density with respect to the CO₂ flux relationship. In our study, we assume on stoichiometric grounds that the O₂ flux density is equal and opposite to the CO₂ molar flux density. However, we do not have observational evidence that this occurs at the leaf level. Besides the Stefan (molar) flow, additional processes such as thermodiffusion (Curtiss and Hirschfelder, 1949; Griffani et al., 2024) may also contribute to gas transport, particularly in the presence of leaf–air temperature gradients. In the present work, we did not account for these thermal effects, but they may become relevant for leaves exposed to strong radiation loads or under fluctuating microclimatic conditions, such as cloud passages or sunflecks (van Diepen et al., 2025). An open question therefore remains regarding the relative importance of thermodiffusion in comparison to Stefan flow and molecular diffusion. To this end, a future aim is to represent O₂ biophysical processes independently from CO₂ at all relevant atmospheric scales.

One key prediction in our study – that the internal stomatal O₂ mole fraction is typically lower than the atmospheric O₂ mole fraction (Eq. 26, Fig. 2) – offers an opportunity to experimentally test the underlying multicomponent diffusion theory on which it is based. Such a test would require the precise measurement of internal O₂ mole fraction under stagnant or low Reynolds number atmospheric conditions and without photosynthesis. Figure 2 suggests that such measurements would be best carried out under a wide range of relative humidity, and expecting the larger gradient values under low relative humidity values. Experimental deviations from the predictions shown in Fig. 2 would also be of value, because they would signal missing physics not included in our current approach – which assumes the ideal gas approximation and isothermal conditions –; crucially, the underlying theory of the Stefan-Maxwell equations provides a rigorous thermodynamic framework for including such extensions (e.g. to non-ideal gases and thermal forces).

6 Conclusions

Oxygen transport between the sub-stomatal cavity and the atmosphere exhibits two key characteristics that require particular consideration: (i) the relatively high mole fraction of oxygen compared to other primary components of the photosynthesis system, and (ii) the significant influence of the water molar flux density (transpiration) on oxygen transport. The Stefan-Maxwell equations provided a rigorous physical framework for examining those characteristics quantitatively. Based on these equations, our findings emphasize the interplay between fluxes, gradients, and diffusivity coefficients within the quaternary N₂-O₂-H₂O-CO₂ system, and provide new insights into the dynamics of oxygen stomatal gas exchange. The main conclusions are:

1. The water vapour molar flux density generates a drag force that drives the movement of oxygen from the stomatal interior to the atmosphere. The large water-driven Stefan or molar flow of oxygen is almost (but not exactly) counterbalanced by a diffusive O₂ flux from the atmosphere to the stomatal interior.

2. We identify three distinct regimes of O₂ transport depending on the ratio of water to oxygen molar flux densities. These regimes are distinguished by the relative orientations of the O₂ molar flux density, the diffusive component, and gradient in O₂ mole fraction between the sub-stomatal cavity and the atmosphere. Under typical conditions, oxygen moves uphill from stoma to the atmosphere from low to high O₂ mole fraction, whereas the diffusive component of O₂ transport moves from the atmosphere to the stomata (downhill).

3. For this typical regime, we derived a simple analytical expression for the ratio of the internal to atmospheric O₂ mole fractions as a function of relative humidity, which explicitly accounts for the drag effect of water-driven Stefan flow. This expression offers the possibility to interpret and explain experiments measuring simultaneous O₂ mole fractions in the stomata and the atmosphere, grounded in underlying multicomponent diffusion theory. If future experiments indeed, measure a lower stomatal internal oxygen mole fraction compared to the atmosphere, our results would offer a robust physical explanation.

Our theoretical results establish a baseline for developing a general multicomponent theory of stomatal gas exchange, which can also be extended to other carbon tracers, such as stable isotopologues. By incorporating isotopic variants of CO₂ (e.g., ¹³CO₂, C¹⁸OO), the framework can help disentangle the relative contributions of different biochemical and physical processes – such as diffusion, carboxylation, and photorespiration – to overall gas exchange. This extension opens new avenues for interpreting isotopic measurements in both laboratory and field settings, providing a physically sounded representation of the transport of carbon tracers across scales and more comprehensive understanding of leaf-level carbon and oxygen dynamics under varying environmental conditions.

Appendix A: List of variables, symbols, values and units

Table A1 presents the definitions of the variables used in this study, including their symbols, units, and characteristic values applied in the calculations. The calculations are performed at the stomatal level, with units of [mol s⁻¹], while the input variables are based on representative observations taken at the leaf level, where the units are [mol m⁻² s⁻¹] (Jarman, 1974; Vilà-Guerau de Arellano et al., 2020). Leaf-level values are downscaled to the stomatal level by dividing them by the leaf stomatal density, which represents the number of stomata per unit area of leaf. “water” is referring to water vapour.

Table A1 Definitions of the variables, symbols, characteristic values, and units used in the study. Molecular diffusivities are calculated at 293 K using the expression proposed by Reid et al. (1987).

¹ All the molecular diffusivities are calculated Reid et al. (1987).

Appendix B: Geometry and scaling of one-dimensional gas exchange between stomata and the atmosphere

Figure B1 shows the one-dimensional representation of the stomata and the atmosphere in which the Stefan-Maxwell equations are applied. The transport of all the species is through the stomatal pore that connects between the sub-stomatal cavity and the atmosphere. This geometry and assumptions are similar to the binary and quaternary mixtures (see Fig. 1).

Figure B1 Sketch of the geometry of the stomata and atmosphere systems where the oxygen exchange takes place. The spatial one-dimensionality is represented by the coordinate l and the characteristic length scale L. Based on observations, we assume that this length scale is larger than the radius of the stomatal pore r, but it is not a theoretical requirement.

Appendix C: Non-dimensional O₂ governing equation including Stefan flow effects

A key question in our analysis is to determine when a multicomponent diffusion representation is necessary to accurately describe oxygen transport. To address this, we consider the one-dimensional molar balance equation for O₂ and non-dimensionalize it using characteristic length and velocity scales. For the binary system of O₂ diffusing in stagnant N₂ (Sect. 3), Eq. (13) becomes:

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial t}}=-{\displaystyle \frac{\partial \left(x_{\mathrm{o}}u\right)}{\partial l}}+{\displaystyle \frac{\partial }{\partial l}}\left(D_{\mathrm{no}}\left({\displaystyle \frac{\mathrm{d}{x}_{\mathrm{o}}}{\mathrm{d}l}}\right)\right)+r_{\mathrm{o}}.\end{array}$$

in which we use Eq. (17) to represent the diffusive molar flux density (f_o) as a function of the gradient of oxygen $\left(\frac{\mathrm{d}{x}_{\mathrm{o}}}{\mathrm{d}l}\right)$. This equation has units of (mol_o mol${}_{\mathrm{air}}^{-\mathrm{1}})$ s⁻¹. Note that in the multicomponent diffusion approach, this equation must be solved simultaneously with the governing equations for u (Navier-Stokes equation) and for the mole fraction of water vapour x_w.

If Eq. (C1) is made non-dimensional using the characteristic length ℒ and velocity 𝒰 inertial scales, the molecular diffusion term will be divided by the Reynolds number, defined as the ratio of inertial acceleration to the molecular diffusion (ν_a) and it reads $Re=\mathcal{L}\mathcal{U}/{\mathit{\nu }}_{\mathrm{a}}$; and the Schmidt number, defined as the ratio of momentum diffusivity (kinematic viscosity) to mass diffusivity ($Sc={\mathit{\nu }}_{\mathrm{a}}/D_{\mathrm{no}}$) (Vilà-Guerau de Arellano et al., 2004). To derive this non-dimensional equation, we introduce the following non-dimensional variables: $\stackrel{\mathrm{̃}}{t}=(\mathcal{U}/\mathcal{L})t$, $\stackrel{\mathrm{̃}}{u}=u/\mathcal{U}$, and $\stackrel{\mathrm{̃}}{l}=l/\mathcal{L}$. Substituting these non-dimensional variables into Eq. (C1) and dividing the equation by the factor $(\mathcal{L}/\mathcal{U})$, we obtain the non-dimensional equation for oxygen. The final non-dimensional equation for oxygen is:

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial \stackrel{\mathrm{̃}}{t}}}=-{\displaystyle \frac{\partial \left(x_{\mathrm{o}}\stackrel{\mathrm{̃}}{u}\right)}{\partial \stackrel{\mathrm{̃}}{l}}}+{\displaystyle \frac{\partial }{\partial \stackrel{\mathrm{̃}}{l}}}\left((ReSc{)}^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{d}{x}_{\mathrm{o}}}{\mathrm{d}\stackrel{\mathrm{̃}}{l}}}\right)\right)+{\stackrel{\mathrm{̃}}{r}}_{\mathrm{o}},\end{array}$$

where the source term is ${\stackrel{\mathrm{̃}}{r}}_{\mathrm{o}}=r_{\mathrm{o}}(\mathcal{L}/\mathcal{U})$. The diffusive term (second term on the right-hand side) is now expressed as a function of the Reynolds number, defined as the ratio of inertial acceleration ($Re=\mathcal{L}\mathcal{U}/{\mathit{\nu }}_{\mathrm{a}}$), and the Schmidt number, defined as the ratio of momentum diffusivity (kinematic viscosity) to mass diffusivity ($Sc={\mathit{\nu }}_{\mathrm{a}}/D_{\mathrm{no}}$). Additionally, the source/sink term is made dimensionless by introducing a representative source/sink scale.

Using typical values for the flow velocity near the leaf surface (ranging from 0.05 to 0.15 m s⁻¹ taken by under greenhouse conditions (Kimura et al., 2020), the length scale (order millimeters) and a kinematic viscosity of air $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹, we obtain Reynolds numbers on the order of 1 to 100. For larger leaves (0.1 m) and more realistic wind conditions – see, for instance, values of 1 m s⁻¹ taken at the canopy top in the Amazonian rain forest (González-Armas et al., 2025) – the Reynolds number is 6666. These greenhouse and forest Reynolds number values indicate that the flow remains laminar around the leaf surface. As a result, we need to retain the molecular diffusion term (2nd right-hand side term) in the governing equation for O₂, and describe this O₂ transport process using a multidiffusion component as described here.

Appendix D: Equivalence ternary and quaternary systems

In this Appendix, we study the equivalence between the equations derived for the quaternary system (18a), (18b) and (18c), and the ternary-system equations derived by Jarman (1974) and further applied by von Caemmerer and Farquhar (1981). We start by rewriting the gradients of water and carbon dioxide as a function of the molar flux densities. The general expression for water vapour is as follows:

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}{\displaystyle \frac{\partial x_{\mathrm{w}}}{\partial l}}& ={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}[-({\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{ow}}}}+{\displaystyle \frac{x_n}{D_{\mathrm{wn}}}}+{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{wc}}}})F_{\mathrm{w}}+{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{ow}}}}{F}_{\mathrm{o}}\\ & +{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{wn}}}}{F}_{\mathrm{n}}+{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{wc}}}}{F}_{\mathrm{c}}]\end{array}\end{array}$$

and the gradient for carbon dioxide reads:

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}{\displaystyle \frac{\partial x_{\mathrm{c}}}{\partial l}}& ={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}[-({\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{ow}}}}+{\displaystyle \frac{x_n}{D_{\mathrm{wn}}}}+{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{wc}}}})F_{\mathrm{c}}+{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{oc}}}}{F}_{\mathrm{o}}\\ & +{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{cn}}}}{F}_{\mathrm{n}}+{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{wc}}}}{F}_{\mathrm{w}}]\end{array}\end{array}$$

To connect and compare with expressions B9 and B10 at (von Caemmerer and Farquhar, 1981) we need to introduce the following assumptions: (1) in the ternary system used by Jarman (1974); von Caemmerer and Farquhar (1981), the fluxes F_o and F_n are grouped in a flux of air F_a. There is a mole fraction gradient of air across the stomata. However, they assume there net flux of air is zero, i.e. F_a=0 at the interface by Jarman (1974); von Caemmerer and Farquhar (1981). (2) the mole fraction of oxygen and nitrogen are added $x_{\mathrm{o}}+x_{\mathrm{n}}=x_a$, (3) the diffusion coefficients are related as follow: D_ow=D_nw and D_oc=D_nc, and (4) we approximate the partial derivatives of x_o, x_w and l with increments Δ.

Using these previous assumptions, we define the conductance of water vapour in air, carbon dioxide in air, and carbon dioxide in water vapour as follows:

$$\begin{array}{}\text{(D3)}& g_{\mathrm{aw}}={\displaystyle \frac{c_{\mathrm{t}}{D}_{\mathrm{wa}}}{\mathrm{\Delta }l}}{g}_{\mathrm{ac}}={\displaystyle \frac{c_{\mathrm{t}}{D}_{\mathrm{ca}}}{\mathrm{\Delta }l}}{g}_{\mathrm{ac}}={\displaystyle \frac{c_{\mathrm{t}}{D}_{\mathrm{wc}}}{\mathrm{\Delta }l}}.\end{array}$$

In applying these assumptions, the equations become the following:

$$\begin{array}{}\text{(D4)}& \mathrm{\Delta }{x}_{\mathrm{w}}=-\left({\displaystyle \frac{x_a}{g_{\mathrm{aw}}}}+{\displaystyle \frac{x_{\mathrm{c}}}{g_{\mathrm{wc}}}}\right)F_{\mathrm{w}}+{\displaystyle \frac{x_{\mathrm{w}}}{g_{\mathrm{wc}}}}{F}_{\mathrm{c}},\end{array}$$

and

$$\begin{array}{}\text{(D5)}& \mathrm{\Delta }{x}_{\mathrm{c}}=-\left({\displaystyle \frac{x_a}{g_{\mathrm{ac}}}}+{\displaystyle \frac{x_{\mathrm{w}}}{g_{\mathrm{wc}}}}\right)F_{\mathrm{c}}+{\displaystyle \frac{x_{\mathrm{c}}}{g_{\mathrm{wc}}}}{F}_{\mathrm{w}}.\end{array}$$

These expressions are equivalent to the ones found by Jarman (1974) (see expressions of the CO₂-gradient at page 932). However, to obtain the final equivalence with Eqs. (B9) and (B10) at von Caemmerer and Farquhar (1981), we need to modify the notation. von Caemmerer and Farquhar (1981) defined water vapour and carbon dioxide gradients and fluxes in the following way: $F_{\mathrm{w}}=\mathrm{\Delta }{x}_{\mathrm{w}}=x_{\mathrm{wi}}-x_{\mathrm{wa}}$ and $F_{\mathrm{c}}\approx \mathrm{\Delta }{x}_{\mathrm{c}}=-x_{\mathrm{ca}}+x_{\mathrm{ci}}$. Following this sign convention, we obtain:

$$\begin{array}{}\text{(D6)}& x_{\mathrm{wi}}-x_{\mathrm{wa}}=\left({\displaystyle \frac{x_a}{g_{\mathrm{aw}}}}+{\displaystyle \frac{x_{\mathrm{c}}}{g_{\mathrm{wc}}}}\right)F_{\mathrm{w}}+{\displaystyle \frac{x_{\mathrm{w}}}{g_{\mathrm{wc}}}}{F}_{\mathrm{c}},\end{array}$$

and

$$\begin{array}{}\text{(D7)}& x_{\mathrm{ci}}-x_{\mathrm{ca}}=-\left({\displaystyle \frac{x_a}{g_{\mathrm{ac}}}}+{\displaystyle \frac{x_{\mathrm{w}}}{g_{\mathrm{wc}}}}\right)F_{\mathrm{c}}-{\displaystyle \frac{x_{\mathrm{c}}}{g_{\mathrm{wc}}}}{F}_{\mathrm{w}}.\end{array}$$

The late expressions are identical to Eqs. (B9) and (B10) at von Caemmerer and Farquhar (1981) probing that the ternary solutions derived by Jarman (1974) can be retrieved from the more general fourth-equation systems here derived.

Appendix E: Quaternary versus binary mixture

From Eq. (18b), and after simplifying by assuming the net flux of nitrogen F_n=0, we group the various terms as functions of the mole diffusive fluxes F_o, F_w and F_o. This allows us to derive the correction contributions to the single mole diffusive expression that follows Fick's law. Rearranging Eq. (18b), we obtain:

$$\begin{array}{}\text{(E1)}& \begin{array}{rl}{\displaystyle \frac{\partial x_{\mathrm{o}}}{\partial l}}& ={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}[-(\stackrel{\mathit{\text{a}}}{\overbrace{{\displaystyle \frac{x_{\mathrm{n}}}{D_{\mathrm{on}}}}}}+\stackrel{\mathit{\text{b}}}{\overbrace{{\displaystyle \frac{x_{\mathrm{w}}}{D_{\mathrm{ow}}}}}}+\stackrel{\mathit{\text{c}}}{\overbrace{{\displaystyle \frac{x_{\mathrm{c}}}{D_{\mathrm{oc}}}}}})F_{\mathrm{o}}\\ & +\left(\stackrel{\mathit{\text{d}}}{\overbrace{{\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{ow}}}}}}\right)F_{\mathrm{w}}+\left(\stackrel{\mathit{\text{e}}}{\overbrace{{\displaystyle \frac{x_{\mathrm{o}}}{D_{\mathrm{o}}}}}}\right)F_{\mathrm{c}}]\\ & ={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{t}}}}[-(a+b+c)F_{\mathrm{o}}+d{F}_{\mathrm{w}}+e{F}_{\mathrm{c}}]\end{array}\end{array}$$

Following the procedure outlined by Jarman (1974) to calculate the correction terms, we determine the percentage of the correction factor for the O₂-gradient. We take term (a) as the reference, representing the interaction between oxygen molecules relative to the compound with the larger mole fraction – in this case, nitrogen. For the remaining terms, we assume that the correction factors are defined as follows:

$$\begin{array}{}\text{(E2a)}& {\displaystyle }\begin{array}{rl}{C}_{\mathrm{1}}& =\mathrm{100}{\displaystyle \frac{b}{a}}\mathit{\%}=\mathrm{100}{\displaystyle \frac{x_{\mathrm{w}}{D}_{\mathrm{on}}}{x_{\mathrm{n}}{D}_{\mathrm{ow}}}}\mathit{\%}{C}_{\mathrm{2}}=\mathrm{100}{\displaystyle \frac{c}{a}}\mathit{\%}\\ & =\mathrm{100}{\displaystyle \frac{x_{\mathrm{c}}{D}_{\mathrm{on}}}{x_{\mathrm{n}}{D}_{\mathrm{oc}}}}\mathit{\%}\end{array}\text{(E2b)}& {\displaystyle }\begin{array}{rl}{C}_{\mathrm{3}}& =-\mathrm{100}{\displaystyle \frac{d}{a}}\mathit{\%}=-\mathrm{100}{\displaystyle \frac{x_{\mathrm{o}}{D}_{\mathrm{on}}}{x_{\mathrm{n}}{D}_{\mathrm{ow}}}}{\displaystyle \frac{F_{\mathrm{w}}}{F_{\mathrm{o}}}}\mathit{\%}{C}_{\mathrm{4}}\\ & =-\mathrm{100}{\displaystyle \frac{e}{a}}\mathit{\%}=-\mathrm{100}{\displaystyle \frac{x_{\mathrm{c}}{D}_{\mathrm{on}}}{x_{\mathrm{n}}{D}_{\mathrm{oc}}}}{\displaystyle \frac{F_{\mathrm{c}}}{F_{\mathrm{o}}}}\mathit{\%}\end{array}\end{array}$$

Table E1 shows the correction values using the values of the molar flux densities and molecular diffusivity coefficients shown in Table A1. The very large value of the coefficient C₃ corroborates the dominance of the Stefan flow.

As shown in Table E1 only C₃ is the largest correction factor that strongly depends on the water vapour molar flux density. These results corroborate the key role played by the Stefan flow as shown in Table 1.

Table E1 Correction factors to the O₂ molar flux densities based on Eqs. (E2a) and (E2b) as a function of three different water vapour molar flux densities. The correction factors are dimensionless and in %. In bold the larger correction factors.