Disentangling Turbulent Gas Diffusion from Non-diffusive Transport in the Boundary Layer

An analysis based on the law of linear momentum conservation demonstrates unequivocally that the mass fraction is the scalar whose gradient determines gas diffusion, both molecular and turbulent. It illustrates sizeable errors in previous micrometeorological definitions of the turbulent gas flux based on fluctuations in other scalars such as the mixing ratio or density. In deference to conservation law, we put forth a new definition for the turbulent gas flux. Net gas transport is then defined as the sum of this turbulent flux with systematic transport by the mean flow. This latter, non-diffusive flux is due to the net upward boundary-layer momentum, a Stefan flow forced by evaporation, which is the dominant surface gas exchange. A comparison with the traditional methodology shows exact agreement between the two methods regarding the net flux, but with the novelty of partitioning gas transport according to distinct physical mechanisms. The non-diffusive flux is seen to be non-negligible in general, and to dominate turbulent transport under certain conditions, with broad implications for boundary-layer meteorology.


Introduction
Surface exchange of atmospheric constituents such as greenhouse gases (GHGs) has gained in importance in the twenty first century, and is best assessed in the turbulent surface layer. Applying the intricate eddy-covariance technique, scientists from many disciplines operate "flux towers" worldwide, supplying key biogeochemical flux data at the ecosystem scale for numerous scientific purposes (Baldocchi 2020). Flux-tower data help explain and predict evolving atmospheric concentrations of CO 2 (used hereinafter as a proxy for GHGs in general), and enable land-use managers to cool climate forcing as motivated by the Kyoto and Paris climate agreements. Flux towers directly estimate the evaporation rate (E) over spatial domains that scale up readily, helping to constrain the hydrological cycle. For these and other reasons, networks have grown to include many hundreds of towers around the globe (Papale 2020), and strive to standardize eddy-covariance methodologies internationally (e.g., Franz et al. 2018;Metzger et al. 2019).
Despite the adoption of standard data processing steps, however, eddy-covariance methodologies are not yet definitive. Remarkably, even labelling an individual eddy as either rich or poor in CO 2 remains a matter of some debate. Specifically, the exact definition of the gas index-or "concentration"-to use in eddy covariance calculations is not yet resolved, as the following brief history demonstrates.
Early flux-tower researchers quantified GHG fluxes via a covariance between vertical motion and whatever gas index instruments could provide. Initially this was the constituent partial pressure (p c , Swinbank 1951;Dyer and Maher 1965;Hicks 1970), but with the emergence of infrared gas analyzers (IRGAs), the eddy covariance scalar shifted to the molar fraction (χ c , from closed-path IRGAs; Desjardins and Lemon 1974) or density (ρ c , from open-path IRGAs; Hyson and Hicks 1975;Jones and Smith 1977;Ohtaki and Matsui 1982). Then researchers realized that turbulent CO 2 fluctuations could arise due to surface exchange, not only of CO 2 , but also of water vapour and/or heat. A debate ensued over "density corrections" (Jones and Smith 1977;Bakan 1978), with consensus endorsing the Webb et al. (1980; hereafter WPL80) definition of the turbulent flux based on the covariance with the CO 2 mixing ratio (r c ). Yet the debate persisted (Fuehrer and Friehe 2002), and some micrometeorologists continued to express eddy fluxes in terms of the covariance with ρ c (Lee 1998;Finnigan et al. 2003;Finnigan 2009).
If uncertainty exists regarding the exact definition of the turbulent CO 2 flux, it is both obscured by imprecise language and entangled with the exact definition of the average vertical velocity (we use the synonyms "average" and "mean" interchangeably). Webb et al. (1980) make no clear distinction between the turbulent flux and the net flux of CO 2 . Aside from a footnote, the term "turbulent flux" appears only on their first page defining it as their paper's objective, which they specify in terms of fluctuations in r c , but as requiring corrections if ρ c fluctuations are used. Throughout the rest of the paper, they make no distinction and simply refer to "the flux". However, given that it is their purported mean vertical velocity (their Sect. 3) that defines the WPL corrections, researchers have drawn different conclusions as to the physical meaning of the corrected flux. Some interpret the WPL corrections as integral parts of the turbulent flux (Massman and Lee 2002;Lee et al. 2004;Leuning 2004;Ibrom et al. 2007) but to others (e.g., Wyngaard 1990;Liebethal and Foken 2003;Finnigan 2009) the WPL terms define transport by the mean flow that is distinct from the turbulent flux. This entanglement between the turbulent flux and the average velocity was noted by Paw U et al. (2000), who postulated the existence of two schools of thought in micrometeorology regarding turbulent transport.
From this, we can deduce with certainty that at least one of these two interpretations leads to an erroneous specification of turbulent diffusion. This is clear when considering an updraft (w′ = 0.1 m s −1 ) whose state variables combine to yield different signs for the fluctuations ρ c ′ and r c ′, as in Table 1. According to Wyngaard (1990), Liebethal and Foken (2003), and Finnigan (2009), this rising air with low CO 2 density contributes to downward turbulent CO 2 transport of w′ρ c ′ = − 11.3 µmol m −2 s −1 and any corrections for density effects would affect not turbulent transport but rather transport by the mean vertical flow. By contrast, the WPL80 specification of the turbulent flux (as we interpret it; see also Massman and Lee 2002;Ibrom et al. 2007) would have it contributing to upward turbulent transport of CO 2 with ̄ w′r c ′ = + 1.1 µmol m −2 s −1 . (The overbar denotes the average of ρ, the air density; Appendix 1 derives transport magnitudes and contextual information for the eddy described in Table 1.) Since these magnitudes differ by 12.4 µmol m −2 s −1 in terms of turbulent CO 2 transport, the meaning of the WPL corrections and their relationship to turbulent transport is disputed by a wide margin. Thus, the question remains: just what influence does this upward-moving, CO 2 -containing eddy exert on exchange of CO 2 by turbulent diffusion? In other words, how exactly should the turbulent flux of CO 2 be specified?
This paper aims to answer these key questions and disentangle diffusion from meanflow transport using an analytical framework and adhering strictly to physical conservation laws. The organization is as follows. Section 2 establishes notation and nomenclature, and the subsequent section reviews the state of knowledge regarding diffusive transport, both molecular and turbulent, including a multidisciplinary literature survey regarding Fick's first law. Section 4 presents analyses of carefully specified, theoretical cases that serve to identify unambiguously the scalar determinant of diffusion. Section 5 proposes a "new" definition of the turbulent CO 2 flux, which sums with non-diffusive transport (Kowalski 2017) to yield the net CO 2 flux. Section 6 applies both the new methodology and the WPL80 approach to representative data from a Mediterranean wetland with a broad range of gas exchange rates and compares outcomes. The penultimate section discusses these results and argues the relevance of the new methodology to different aspects of micrometeorology, followed by general conclusions in Sect. 8.

Mathematical Operators
For some variable x, the overbar denotes its arithmetic average ( , as is common in modern notation. However, x may not necessarily be equal to the exact average ( x ), depending upon the variable x in question and how it is sampled/measured (Kowalski 2012). Following micrometeorological tradition, the prime denotes fluctuations in x about its arithmetic average (i.e., x � = x −x ). By contrast, the double prime denotes fluctuations in x about its exact average (i.e., x �� = x −x).

Gas Variables and Physical Processes
The numerous variables that quantify atmospheric component gases can confuse even experienced scientists. Perceptive meteorology students first notice this regarding humidity, whose relevant variables are described here (and in Table 2, with subscript "v" for water vapour). Since water vapour is a gas, the ideal gas law relates its partial pressure (e) to T and its density (or absolute humidity, ρ v ); density is also relevant for studying dynamics, as buoyancy's determinant. Greenhouse absorption or other radiative extinction depends on the number density (η v ), according to Beer's law. Early meteorologists (Guldberg and Hohn 1876) defined the water vapour mixing ratio (r v ) in order to track phase changes; invariant to changes in temperature (T) and pressure (p), its increases/decreases reflect net evaporation/condensation. Also, the mixing ratio for a dry-air component (r c ) is invariant to water phase changes. However, predicting phase change's direction requires knowledge of the relative humidity (U). Finally, both because it facilitates calculating the partial pressure from atmospheric state, and also for its conservation properties, the molar fraction (χ v ) is popular with plant physiologists (e.g., Farquhar and Cernusak 2012). While these associations between physical processes and variables challenge first-year students, it takes the complicated issue of diffusion to perplex scientific researchers. The challenge of specifying diffusion is underscored by both the number of gas variables that have been used to characterize turbulent fluxes (Sect. 1), and the contradictory versions of Fick's first law that are found in the literature (Sect. 3). In this regard, the appearance in Table 2 of the constituent mass fraction (f c ; or q for water vapour, termed the specific humidity) foreshadows this paper's conclusions.
Regarding nomenclature, a number of terms can describe gas transport that is due to the average fluid motion or system velocity, and not to diffusion. Irrespective of scalar gradients, the direction of such transport is in the direction of fluid momentum (i.e., downwind). Table 2 Variables quantifying gas abundance for water vapour (subscript v) and a more general constituent (c), along with their symbols, definitions, relevance, and units Comma-separated entries denote names and symbols for water vapour, and then other constituents. In the ideal gas law (p c = ρ c R c T), the constituent-specific gas constant is R c and the temperature T. The saturation mixing ratio is r v,sat . Since all molecules are equal to the ideal gas law, χ c is a fraction equivalently of pressure (µPa Pa −1 ), volume (ppmv), or moles (ppm); hence, χ c facilitates the calculation of p c from the ambient pressure p. *Some researchers, particularly in ecophysiology, define the molar fraction with reference to dry air (χ c = η c /η d ) In order to distinguish such transport from diffusion due to random motions, we will term it as due to "system motion". Other acceptable terms to describe transport caused by system motion include "gross transport", "mass flow", "convective flux" and "non-diffusive flux". However, we will avoid the term "advection" because it can be defined to depend directly on scalar gradients, as a scalar product having no direction, and thus not as transport per se (Kowalski 2017).
With regard to molecular diffusion in particular, let us recall the words of one of the great physicists of the twentieth century (Feynman et al. 1977): "We must be careful not to confuse diffusion of a gas with the gross transport that may occur due to convection currents." For some situations, this is obvious. A simple example is the case of a marker such as smoke emitted in a strong wind (large Péclet number). We should never interpret the lack of net upwind transport as implying zero diffusivity, because we must always assess diffusive transport relative to system motion. We shall recall Feynman's idea in Sect. 4 for a situation in which it is far more subtle, yet equally valid, and where it helps greatly to isolate diffusion's scalar determinant.

The State of Knowledge Regarding Gas Diffusion
This section illustrates that neither turbulent nor molecular diffusion has an unambiguous characterization to date.

Turbulent Diffusion
Careful examination of the issue of "Reynolds averaging" shows that micrometeorology has a long tradition of conflating diffusive and systematic transport. For decades, it has been the norm to define the average flux of air as (Priestley and Swinbank 1947), purporting total mass transport to be the sum of systematic and turbulent components, the former due to the mean flow with an alleged average velocity of w . Setting (1) to zero leads to the derivation of an average velocity in the direction of the heat flux (Priestley and Swinbank 1947); an equivalent approach, but using dry air only, underlies the derivation of the WPL corrections (Eqs. 11 and 12 of WPL80). Yet such an account defies Osborne Reynolds' specification of turbulence, which was rooted in the law of linear momentum conservation (LMC). A careful reading of Reynolds (1895) reveals relevant distinctions between the exact average and arithmetic average velocities, denoted here respectively as w and w (contrary to Reynolds' notation). When written using our modern notation, the consistent definition of Reynolds (1895) for the system velocity is based on average system momentum, as This can be rearranged by simple algebra to facilitate comparison with (1) as (3) w =w + 0, highlighting the key facts that total mass transport defines the average flow, and so there can be no net mass transport of air by turbulence.
These deductions are consistent with the law of LMC for a fluid or other system of particles, which defines the system velocity using a mass-weighted average of the individual component velocities, as in any introductory physics textbook (e.g., Giancoli 1984). In setting Eq. 1 to zero and thus alleging a mean velocity in the direction of the heat flux, which underlies the WPL density corrections, micrometeorology has neglected the law of LMC. By contrast, Eq. 2 respects both LMC and Feynman's idea since it exactly defines the system velocity against which we must assess transport due to random motions. In summary, while diffusion can reorder system mass, and certainly can cause transport of different constituents in different directions, any net mass flow of the system inherently defines the average velocity and therefore represents a background against which molecular or turbulent motions erratically randomize.
In addition to the correct specification of average and turbulent velocities in accordance with the law of LMC, there remains the question of which gas index from Table 2 is diffusion's determinant. Section 1 has shown that this remains unresolved at present among micrometeorologists. The basis for describing turbulent diffusion has always been an analogy with transfer by random molecular motions (Richardson 1919). Unfortunately, however, the state of knowledge regarding molecular diffusion's scalar determinant is at least as ambiguous as is that for turbulent diffusion.

Molecular Diffusion
Molecular diffusion is a "rather complex phenomenon" (Batchelor 1967) that many scientific disciplines describe incorrectly and furthermore as if it were straightforward and intuitive. This assertion is proven by the numerous, conflicting definitions of Fick's first law that can be found in scientific textbooks. All define the diffusive flux as proportional to (the negative of) a concentration gradient, but with "concentration" defined variously as the number density (η c ), density (ρ c ), mass fraction (f c ), or molar fraction (χ c ).
Generally, the meaning of "concentration" in Fick's first law varies according to scientific discipline. Many texts in physics (Feynman et al. 1977;Giancoli 1984) and physical chemistry (Mortimer 2000; Atkins and Paula 2006; Hofmann 2018) specify the concentration gradient directing diffusion of constituent (c) as ∇η c , while scientists studying the gas exchanges of biological systems, such as plant physiologists, are likely to write Fick's law using either ∇ρ c (Monteith 1973;Jones 1983) or ∇η c (Nobel 2005). Engineers seem to be largely consistent in specifying Fick's law using ∇f c (Geankoplis 1993;Kreith et al. 1999;Lienhard and Lienhard 2000;Bird et al. 2002), whereas fluid dynamics textbooks exhibit the gamut of variability, specifying the concentration gradient as ∇ρ c (Smits 2000 . Clearly, a definitive specification of diffusion's scalar determinant is needed, and this is particularly so in the case of the atmospheric boundary layer, which is both compressible and variable in composition (hence molecular mass) such that any two of the four gradients (∇η c , ∇ρ c , ∇f c , or ∇χ c ) may have opposite signs.

Analyses
The aim of this section is to definitively describe diffusion, both turbulent and molecular, specifically in terms of the appropriate constituent-gas abundance index (i.e., the "concentration" that figures in flux-gradient theories such as Fick's first law). Several theoretical examples, which we term cases, are carefully defined to eliminate incorrect scalar indices from consideration. Let us keep in mind Feynman's requirement to assess diffusive effects relative to system motion, even when the stream velocity is quite small, as in the following subtle examples.

General Framework
The following basis will be used commonly in a set of four cases to be studied. Let us consider an isothermal system (I) that is a cubic chamber completely isolated from its environment, which is the rest of the universe (II). The components of I are its six walls (I.A) and the fluid therein contained (I.B), which is an ideal gas composed of an equal number of moles of heavier gas (I.B.1; left) and lighter gas (I.B.2; right), initially unmingled on either side of a halfway division (dashed line) as depicted in Fig. 1. If we were to suppose the division to be a barrier, with equal volume and number of molecules on either side, then the ideal gas law would imply that each component gas exert the same pressure on the barrier. Instead, however, let us exclude a barrier and examine how the two fluid components mix in the absence of any forces or fields external to I.
The following analyses, valid whether diffusion is purely molecular or predominantly turbulent (whatever the Péclet number), examine the mechanisms that mix I.B. We know from both experience and the second law of thermodynamics that I.B will evolve towards a well-mixed state of equilibrium (no gradients), with uniform values for every gas species of η c , ρ c , f c , and χ c (as well as the mixing ratio, r c ) throughout the container, and one might suppose this to be an instance of pure diffusion. However, three cases of particular interest will show that this is not generally so and will furthermore help to identify the scalar concentration whose gradient determines diffusion. Subsequently, a fourth case will highlight the relevance of these analyses to the atmospheric boundary layer.

Theoretical Cases
The first case shows that systematic transport and diffusion generally operate in tandem to achieve mixing, even for the kinematically equilibrated framework described above. In CASE 1, both fluid components are hydrogen (H 2 ) but with an isotopic distinction, with deuterium (4 g mol −1 ) on the left and protium (2 g mol −1 ) on the right. These two gases are chosen merely for illustrative purposes, with the difference in mass being key to illustrating the transport mechanisms involved in mixing.
Paying heed to Feynman, we should characterise system motion before examining the effects of diffusion. The asterisk in Fig. 1 marks the initial centre of mass of I.B, at two-thirds of the distance from the centre of mass of protium to that of deuterium. Mixing shifts the centre of mass of fluid I.B to the centre of the container, undeniably defining systematic motion to the right. In other words, a mass flow develops within the container that tends to drag along any objects embedded in the flow, including all gas molecules. This mass flow is due, not to a net flow of molecules, but rather to a net flow of H 2 neutrons: two per deuterium molecule, versus zero per protium molecule.
Newtonian physics aptly describes why the fluid I.B, initially at rest, accelerates to the right, and then decelerates to come to a rest at a new position. Let a barrier (at the dashed line in Fig. 1) that is not part of system I initially separate the two gas components I.B.1 and I.B.2, but vanish such that mixing can begin at the instant t 0 . (Prior to t 0 , the pressure is uniform throughout system I.B; afterwards, the pressure cannot be defined until a new equilibrium is reached.) Upon the barrier's disappearance, the force per unit area on the left and right container walls (I.A) remains constant until some later instant t 1 marking the first collision of a (faster) protium molecule with the left container wall (rather than the right wall, where it would have hit had the barrier not been removed). This and other novel protium collisions with the left wall push the container I.A to the left (action), and by Newton's third law the fluid I.B receives a force to the right (reaction) and thus accelerates to the right, according to Newton's second law. At some later instant t 2 , the first collision of a (slower) deuterium molecule with the right wall pushes I.A to the right (action), causing a force to act on I.B to the left (reaction) and therefore its deceleration as it begins to come to a rest at its final position.
Taking a step back to examine the larger picture, we note that the system (I) has no interaction with the rest of the universe (II), which is its environment. Therefore, it undergoes no accelerations, although the same is not true about its components. As noted above, the fluid (I.B) moves to the right and so the container (I.A) must move to the left, such that the system (I) remains at rest, according to Newton's first law. This illustrates why a mass-centred framework is essential for defining the fluid velocity; a molecule-centred (kinematic) framework would fail to explain the movement of the container in the context of Newton's laws. If the mass of the gas is negligible in comparison with that of the container (or if the container is fixed to a gigantic mass such as the Earth), then the container can be approximated as at rest while its contents shift to the right.
This case shows how dynamics outperforms kinematics in precisely describing motion, much as Newton improved upon Galileo's description of falling bodies. A kinematic description of CASE 1, based on the ideal gas law, would erroneously suggest no net fluid motion during the swap of protium and deuterium molecules, neglecting the importance of mass in dynamics. Both the ideal gas law and the state variables p and T are kinematic in nature, being independent of molecular mass, as the following summary of kinetic theory shows. The two isothermal isotopes have identical molecular kinetic energies (i.e., Ts) when the root-mean-square velocity of the protium molecules is a factor √ 2 (about 41%) greater than that of the deuterium molecules, which have twice the mass. Similarly, when colliding with a container wall during equilibrium, the average change in momentum of a deuterium molecule is 41% greater than that of a protium molecule, but the latter collisions occur 41% more often such that the average momentum transfer per unit time and per unit wall-surface area (i.e., p) is the same for each isotope. In other words, the ideal gas law describes molecular kinematics, but is an inapt basis for appreciating dynamics.
With systematic motion understood, we can begin to address the issue of random transport (diffusion). Relative to the systematic motion or mass flow described above, protium diffuses upstream and deuterium downstream. The final equilibrium is achieved for each component at the same moment, because protium diffusion is more rapid, in accordance with Graham's law. More to the point, adding traces of helium (He; 4 g mol −1 ) to the mixture and contrasting its scalar gradients will reveal the correct scalar determinant of diffusion.
The second case clearly eliminates several gas indices from consideration as the scalar whose gradient determines diffusion. In CASE 2, fluid I.B has a uniform helium molar fraction χ He = 10 ppm, with the remaining 99.999% of molecules being H 2 , again with deuterium on the left and protium on the right. Regarding scalar gradients, this implies ∇χ He = 0, ∇η He = 0 and ∇ρ He = 0, but ∇f He > 0. In terms of dynamics, CASE 2 is practically identical to CASE 1; with the substitution of mere traces of He, the molecular mass of gas I.B.1 is unchanged, and that of gas I.B.2 increases only marginally. Figure 1 still serves perfectly, although the initial position of the asterisk must now shift imperceptibly to the right, versus CASE 1. Here it is very illustrative to compare transport types for He as the fluid I.B is mixed. As noted above, the mass flow (net flux of H 2 neutrons) drags He to the right. However, the overall position of He does not change, since both the initial conditions and the well-mixed equilibrium require uniform distributions of χ He , ρ He and η He throughout I.B. Therefore, the mass flow of He to the right must be offset by He diffusion to the left (i.e., the motion of He relative to fluid I.B) in order for net He transport to be zero. In short, the He remains in place by diffusing upstream. Thus, CASE 2 clearly demonstrates a situation with a negative diffusive flux due to a positive ∇f He , but null values of ∇χ He , ∇η He and ∇ρ He , demonstrating that only ∇f He can correctly determine diffusion.
The third case reinforces the relevance of mass in a situation with no He diffusion. In CASE 3, fluid I.B has a uniform f He = 10 mg kg −1 , with the other 999,990 mg kg −1 being H 2 that is once again deuterium on the left and protium on the right. Now the gradients ∇χ He , ∇ρ He and ∇η He are all negative (with double values of these scalars on the left versus the right) and yet no He diffusion occurs, as is consistent with the null value of ∇f He .
In fact, f He is spatially constant and remains so during the intermingling, even as deuterium diffuses downstream and protium upstream, because f He is conserved through fluid motions (conservation of a mass fraction is demonstrated mathematically in the appendix of Kowalski and Argüeso 2011), excluding the possibility of any He diffusion. Rather, the net He transport to the right is purely systematic in nature, and due simply to the drag on He caused by the net flow of H 2 neutrons.
These cases highlight the need to examine diffusion against a background of mass flow for any constituent (c), and thus from the perspective of dynamics. In this context, the gradients of kinematic variables χ c (µmol mol −1 ) and η c (mol m −3 )-which arise naturally out of ideal gas law applications-fail to indicate diffusion's direction because they neglect the key role played by mass in dynamics. On the other hand, the gradients of absolute scalar amounts ρ c (kg m −3 ) and η c (again) fail to describe diffusion because they focus exclusively on the constituent of interest and neglect the relevant context of system mass. In conclusion, as the lone "concentration" variable that satisfies the requirement of quantifying the scalar relative to the overall fluid mass, the mass fraction (f c ) is diffusion's definitive determinant.
A final case shows the relevance of the preceding analyses to the atmospheric boundary layer. In CASE 4, fluid I.B is composed on the left (I.B.1) of dry air (29 g mol −1 ), and on the right (I.B.2) of moist air (28.8 g mol −1 ; with a water vapour mixing ratio of r v = 10 g kg −1 ), each with a CO 2 mass fraction of f c = 600 mg kg −1 . We can readily calculate the left/right values of the CO 2 mixing ratio (r c ; mg kg −1 ) of 600/606. Similarly, the isothermal left/right values of the CO 2 molar fraction χ c (µmol mol −1 ) of 395.5/392.7 imply a gradient in the CO 2 density ρ c (mg m −3 ), which is more than 0.7% greater on the left. (Appendix 2 calculates these values.) However, from analysis of the foregoing cases, it is clear that this is a situation with no CO 2 diffusion. Thus, the answer to the question posed in Sect. 1 is that, of the two differing interpretations of the relationship between the WPL corrections and turbulent fluxes, neither is correct. In fact, the "data" in Table 1 were generated by taking f c ′ = 0 as a starting point, and thereby ensuring an eddy with no influence at all on turbulent CO 2 diffusion (hence the artificial precision).

The Flaws in Previous Turbulent-Flux Specifications
The above analyses highlight errors in the prevailing micrometeorological definitions of transport by turbulence, which are entangled with systematic transport by the mean flow. Fundamental principles relevant to quantifying surface exchange include conservation of both mass and momentum. While the former has been used to justify defining turbulent CO 2 diffusion based on fluctuations in ρ c ′ (Lee 1998;Paw U et al. 2000;Finnigan et al. 2003), it is the latter that must be applied to distinguish between physical transport processes. For the conditions defined in Table 1, specification based on ρ c ′ (Wyngaard 1990;Liebethal and Foken 2003;Finnigan 2009) underestimates the turbulent CO 2 flux by 11.3 µmol m −2 s −1 because it focuses exclusively on CO 2 and neglects the context of local system mass, demonstrated relevant in the previous section. By contrast, the WPL80 specification based on r c ′ overestimates the turbulent CO 2 flux by 1.1 µmol m −2 s −1 due to several errors.
The first error in the WPL80 flux-partitioning scheme is the use of arithmetic averaging to infer a spurious mean velocity in the direction of the heat flux as a consequence of the erroneous mass flux partitioning in their Eq. 12, which is analogous to our Eq. 1 but for dry air only. In comparison with the most basic definition of average velocity-as displacement divided by time elapsed-this velocity has been shown to be an artefact of inexact averaging (Kowalski 2012). Alternatively, since Eq. 12 of WPL80 applies to dry air (of constant composition), the error of using arithmetic averaging to calculate the average velocity is clear when viewed in the context of statistical sampling theory. Sonic anemometers have fixed sensing volumes, and therefore sample variable population sizes: generally, the colder the eddy, the greater the number of molecules sampled. Giving all samples equal weight as in arithmetic averaging therefore admits sample bias in favour of warm eddies. Since warm eddies tend to rise and cool eddies to descend when the heat flux is positive, this bias causes over-estimation of the average vertical velocity, and is the explanation for the alleged mean vertical motion in the direction of the heat flux.
A second error regards the evaporation-based velocity component defined by WPL80 in their Eq. 14. This term due to the water vapour flux is in error because it represents an attempt to describe dynamics based on the ideal gas law, expressed in Eq. 8c of WPL80. The previous section shows that such a kinematic description, neglecting the lower molecular mass of water vapour (versus dry air), leads to a bias in system momentum. In the absence of temperature fluctuations, Eq. 8c of WPL80 states that an increase in water vapour results in an equal decrease in dry air, on a molar basis. By contrast, we define turbulent perturbations on a (fractional) mass basis. The difference between the two approaches thus includes a factor of proportionality that is the ratio of molecular masses that WPL80 define as µ. For isothermal conditions (no heat flux), this µ is essentially the difference between the Stefan flow (Kowalski 2017) and Eq. 14 of WPL80.
Thirdly, the exclusion of water vapour from the denominator of the CO 2 mixing ratio, which WPL80 use to define the turbulent flux, neglects the intrinsically bilateral nature of diffusion. As the dominant surface exchange process, evaporation promotes not only upward diffusion of water vapour but also downward diffusion of dry air, including CO 2 . This is aptly illustrated by considering the case of steady-state evaporation with null CO 2 exchange, such that the net CO 2 flux (F c ) is exactly zero. Because of space occupied by newly introduced water vapour, air near the evaporating surface is diluted with regard to CO 2 (or any dry air species), whose turbulent flux is therefore downward. However, evaporated water vapour not only dilutes but also displaces other gas species, and so the downward turbulent flux is accompanied by upward systematic transport of CO 2 by the Stefan flow. The two components cancel each other to yield F c = 0. The WPL80 mixing-ratiobased definition gets the net flux correct (unaffected by evaporation, r c is constant in this example), but falsely characterises the net flux as wholly turbulent.

A "New" Definition of the Turbulent CO 2 Flux
The following formulation presents the total CO 2 flux F c (kg m −2 s −1 ) as the sum of two independent (disentangled) fluxes, the non-diffusive flux due to transport by systematic motion (F c,ndiff ) and the turbulent or diffusive flux (F c,diff ) where w s = E − is the Stefan flow velocity (m s −1 ; Kowalski 2017), E the evaporation rate (kg m −2 s −1 ) and ̄c the average CO 2 density (kg m −3 ). Fluctuations (double primes) of both vertical wind velocity and CO 2 mass fractions are calculated about mean values obtained via density-weighted averaging (Sect. 2.1 and Eqs. 8-9 of Kowalski 2012). The final overbar on the right-hand side of (4), representing arithmetic averaging, is justified because it overlies conserved CO 2 (eddy) momentum.
Since the moist air density ( ) is not directly measurable at 10 Hz, it is determined by adding the average air density ( ̄ , computed from thermohygrometer and barometer data) to fluctuations that we estimate from sonic anemometer and gas analyser data. This we do, neglecting turbulent pressure fluctuations (i.e., p = p ), via a procedure that initializes the temperature (T) using the sonic temperature (T s ) and then iterates as follows These steps are the ideal gas law for water vapour (Eq. 5a), where R v is the specific gas constant for water vapour, and the definitions of the specific humidity (Eq. 5b) and sonic temperature (Eq. 5c). For the dataset analysed here, convergence to five or more significant digits for all variables required a maximum of three iterations. From these, the virtual temperature (T v ) is defined as and the ideal gas law for moist air defines Density fluctuations (ρ′) are then determined by subtracting from ρ its arithmetic average, which is not otherwise used.
An analogous formulation is proposed for estimating water vapour fluxes. However, because the computation of w s requires knowledge of the water vapour flux itself, another iterative process is necessary to estimate both E and w s simultaneously, with E initialized to the traditional WPL80 water vapour flux, using the formulation where E is equal to the net water vapour flux (F v ) that is also disentangled into its non-diffusive (F v,ndiff ) and diffusive (F v,diff ) transport components. The iteration is completed when values of both E and w s converge within a specified tolerance. Our tests show that three to five iterations are usually sufficient for convergence to five significant digits.
Finally, and in a similar fashion, the sensible heat flux H (W m −2 ) can be estimated as where c p (J kg −1 K −1 ) is specific heat of moist air at constant pressure.

Experimental Site and Materials
The measurements used to compare the novel flux-calculation methodology with traditional (WPL80) procedures come from a Mediterranean wetland with very large flux (5b) q = 0.622 e p (5c) T = T s (1 + 0.51q) .
The site combines a number of characteristics that make it very suitable for flux-tower measurements. The ground is very level (usually flooded), with an expansive patch of monospecific vegetation (the common reed; phragmites australis) providing hundreds of metres of fetch in every wind direction. The low tower height (6 m; over reeds that grow to 3 m by mid-summer) thus ensures a homogeneous tower footprint. The combination of abundant water and nutrients and southeast Spain's sunny, dry climate make for very large H 2 O and CO 2 flux magnitudes during the growing season. Finally, the site is quite windy, near the head of the Lecrín valley with its numerous wind turbines.
Tower-top measurements at 10 Hz include the wind vector components and sonic temperature from a sonic anemometer (CSAT-3, Campbell Scientific, Logan, UT, USA), and densities of H 2 O and CO 2 from an open-path IRGA (LI-7500; LI-COR Inc., Lincoln, Nebraska, USA) that also measures the atmospheric pressure. Just below at 5-m, a thermohygrometer (HMP 45C, Campbell Scientific, Logan, Utah, USA) reports the air temperature and relative humidity to a data logger (CR3000, Campbell Scientific) that stores half-hour averages.
We compare the traditional (WPL80) versus disentangling methodologies in a consistent coordinate system. This was determined by two sequential rotations to correct for anemometer tilt (Finnigan et al. 2003) and force arithmetic mean velocities (prior to density corrections) to zero: first w (attack angle) and then v (wind direction). Annual integrations of net ecosystem exchange estimated for both methodologies were calculated filling gaps using the marginal distribution sampling technique (Reichstein et al. 2005). We examine experimental data for 2015, with a percentage of rejected CO 2 fluxes of 45% (74% of the rejected data occurred during night-time periods).

Results
For representative fluxes across multiple seasons, we compare net CO 2 fluxes calculated using the new methodology with traditional calculations from the WPL80 methodology. Emphasis here is on fluxes of CO 2 , because both water vapour and sensible heat fluxes from the two methodologies were virtually indistinguishable, and because non-diffusive transport of water vapour is only a tiny fraction of total water vapour transport (Kowalski 2017). In every case, F c represents the flux calculated with the new methodology, and F WPL the traditional calculation.
The net CO 2 flux calculated by the two methodologies agrees quite exactly, as shown in Fig. 2 for a full year (2015) of data. Such agreement extends to annual integration of net ecosystem exchange for the Padul dataset, which for 2015 was a net uptake of 209 gC m −2 for the disentangled methodology, versus 207 gC m −2 for the traditional WPL80 methodology.
Since the two methodologies for calculating the net CO 2 flux, disentangling (F c ) and traditional (F WPL ), agree so perfectly, in the following results only F c is presented. However, in Figs. 3 and 4 the black curve can be interpreted as representing either the net flux from our methodology (F c ) or the turbulent flux defined by WPL80 (F WPL ), which are indistinguishable when plotted. By contrast, the red and blue curves are unique to the disentangling methodology, and represent transport by Stefan flow and by turbulence, respectively.
The two methodologies differ regarding the physical mechanisms of transport (Fig. 3). Compared with the new methodology that is grounded in conservation law, the WPL80  methodology (black) underestimates the magnitude of the daytime turbulent flux by about 10% during early summer. At midday, the net flux of about − 20 µmol m −2 s −1 is disentangled into about − 22 µmol m −2 s −1 of turbulent transport (blue), and + 2 µmol m −2 s −1 of transport by the Stefan flow (red). Other seasons show even greater relative discrepancies between the two methods.
As flux magnitudes evolve over the seasons from their summer evaporative and photosynthetic maxima, the relative decrease in CO 2 uptake overrides the percent reduction in evaporation, such that the relevance of non-diffusive transport escalates (Fig. 4). Comparing Figs. 3a and 4a, from early summer to the transition-to-senescence period that dominates the ecosystem's behaviour as an annual carbon sink (Serrano-Ortiz et al. 2020), daytime CO 2 uptake (negative F c ; black) decreases by nearly an order of magnitude, while the non-diffusive flux (approximately proportional to evaporation; red) falls by only 50%, with the result that the magnitude of non-diffusive CO 2 transport is roughly one-third that of the total CO 2 flux. In this case, regarding downward turbulent CO 2 transport, the WPL80 estimate (peaking near 3 µmol m −2 s −1 ; black) underestimates the turbulent flux (reaching 4 µmol m −2 s −1 ; blue) by about 25% (Fig. 4a). The discrepancies are even greater in late April, when the reeds begin to bud; CO 2 fluxes are then dominated by respiratory decomposition of the previous year's organic material and photosynthesis is only strong enough to reduce midday CO 2 emissions to near zero (Fig. 4b). Non-diffusive CO 2 transport (proportional to evaporation; red) is modest, but still enough to dominate net CO 2 transport (black) during the morning hours, when the turbulent flux (blue) is near zero and is extremely overestimated by the WPL80 methodology. During afternoon hours, the diffusive (blue) and non-diffusive (red) fluxes are of similar magnitude, implying that the WPL80 methodology overestimates turbulent transport by 100%.

A Critique of the Webb, Pearman, and Leuning Methodology
Traditional flux calculations produce accurate determinations of the net CO 2 flux, but leave the mechanisms of physical transport entangled. The net fluxes computed by the two methodologies agree because they use the same equation, however expressed differently, whether as w c in Eq. 19 of WPL80 or as wf c in our Eq. 4. The key difference lies in how that net flux is decomposed into transport by the mean flow, versus transport by turbulence. The WPL80 decomposition into turbulent and non-turbulent ("mean") contributions has neglected the law of linear momentum conservation (LMC) in two important ways. First, the average velocity should be defined consistent with LMC as by Osborne Reynolds (1895). By contrast, w is an arbitrary, near-zero value of w that has no physical relevance. Second, the scalar variable whose fluctuations define diffusive transport is neither the density (ρ c ) nor the mixing ratio (r c ), but rather the mass fraction (f c ). Nonetheless, based as it is on the ideal gas law, the WPL80 methodology provides an accurate accounting of CO 2 molecules exchanged at the surface, such that many fluxes that have been published over the last few decades will not require revision.
It is interesting to note that WPL80 derived "alternative expressions" for fluxes of CO 2 and water vapour that coincide neatly with our decomposition, as long as we suppose turbulent fluctuations to be defined in terms of the mass fraction. Then, Eq. 23 of WPL80 states that the turbulent water vapour flux must be divided by a factor (1 − q) to calculate net evaporation, in agreement with the derivation by Kowalski (2017) that q is the fraction of the vapour flux that is non-diffusive. For the CO 2 flux, again following this supposition, Eq. 22 of WPL80 specifies net CO 2 exchange as the sum of the turbulent flux and the product of the average CO 2 density with a velocity that is equal to the evaporation rate divided by the average air density, mirroring our Eq. 4. Again, the key difference regards exactly how the turbulent flux is defined.
Generalizing from the specific case of CO 2 to other gases, we should note that Eq. 4 specifies the non-diffusive flux as a function of the evaporation rate and the density of the gas in question. It does not depend on the exchange rate of that gas. Thus, the relative importance of non-diffusive transport depends on both gas quantity and reactivity. An extreme example here is that of O 2 , whose atmospheric abundance exceeds that of CO 2 by two orders of magnitude, yet has a comparable exchange rate. As a result, its upward transport by systematic boundary-layer motion often dwarfs its net flux (i.e., |F c,ndiff | ≫ |F c | ), making its diffusion sizeable and downward despite photosynthetic O 2 production at the surface. (Emission of water vapour is so much greater than that of O 2 , that the net effect of surface exchange that of O 2 dilution.) The examination of the relative magnitudes of F c,ndiff and F c,diff in composing F c for GHGs other than CO 2 will be addressed in a future paper.

Implications for Micrometeorology
If some net fluxes published using WPL80 methodologies were accurate, as the above agreement indicates, others that made use of Monin-Obukhov similarity theory (MOST) may include non-negligible errors. The scope of MOST is to relate turbulent-not netfluxes to other boundary-layer variables. Thus, flux-gradient relationships apply to the diffusive (F c,diff ) flux component only, and not strictly to the net flux. The same can be said for numerous applications of second-order moments for scalars fluctuating in turbulence (Detto and Katul 2007), which should be characterized in terms of the statistics of the mass fraction. These include spectra and cospectra, for which similarity should be expected between the temperature, e.g., and the mass fraction-but not necessarily the mixing ratio or density-of the gas of interest.
Spectral corrections are a particular concern in this regard. To illustrate this, let us consider a surface with 90 W m −2 of latent heat flux (E = 2 mmol m −2 s −1 ) and 2 µmol m −2 s −1 of O 2 emission (modest photosynthesis), into air with 70% relative humidity, 20.95% O 2 , and a density of 1.22 kg m −3 . For simplicity, we will specify no sensible heat flux. Table 3 presents Table 3 Total gas fluxes and flux components according to the two methodologies for 2 μmol m −2 s −1 O 2 emissions over an evaporating surface with 90 W m −2 of latent heat flux and no sensible heat flux Dark grey rows (total fluxes) are specified from the text. Light grey rows are calculated according to each methodology. Unshaded rows are calculated simply as the difference of the previous two flux values according to the two methodologies compared in this paper, which disagree substantially regarding O 2 flux decomposition. Now let us hypothesise an open-path instrument that measures O 2 density fluctuations with a response time that is insufficient to track the smallest eddies, unlike the fast-response IRGA measuring water vapour alongside, and so would be incapable of producing the data in Table 3 without correction for high-frequency loss. Let us also suppose that empirical determinations have shown a 10% flux loss, e.g., by applying the instrument's cospectral transfer function by low-pass filtering to fluxes that are measured without losses. Frequency response corrections based on MOST cospectral similarity (e.g., Massman 2000) could be applied directly to turbulent fluxes measured with such a slow sensor, but not to any other flux or component in Table 3, including the net flux and WPL corrections.
As a final implication, recent attempts to partition net gas exchanges into contributing biogeochemical processes-such as breaking net CO 2 exchange into photosynthetic and respiratory components, or separating the net water vapour flux into soil evaporation and transpiration (Thomas et al. 2008;Scanlon and Kustas 2010;Skaggs et al. 2018;Stoy et al. 2019)-are based on MOST but truly applicable only to turbulent scalar fluctuations (i.e., those of the mass fraction). The relevance of non-diffusive fluxes to such applications remains to be explored.

Conclusions
Surface-normal transport of CO 2 can be turbulent and/or non-diffusive, with the net flux being the sum of these two components representing distinct physical transport mechanisms. The law of linear momentum conservation must be respected when distinguishing between these two types of transport. Non-diffusive transport depends on the momentum of air due to water vapour exchange (Stefan flow), and the CO 2 density. Diffusive transport, whether turbulent or molecular, is proportional to the concentration gradient, where the "concentration" must be defined as the constituent mass fraction (f c ). The turbulent flux should be defined as the (density-weighted) covariance between f c and the vertical velocity.
where R d is the specific gas constant for dry air. Finally, we can calculate the CO 2 mixing ratio (r c ) as the ratio of ρ c (measured directly) to ρ d (calculated as above) Then the eddy flux defined by perturbations in r c is straightforward as d w′r c ′, with d from (11), using data from column 2 in Table 1, w′ as specified (0.1 m s −1 ), and r c ′ determined via the difference between mixing ratios calculated for the updraft and the boundary-layer average.
For context, we note the energy fluxes associated with the eddy in Table 1. Calculations for water vapour, similar to those above for CO 2 , reveal a latent heat flux of 210 W m −2 while the sensible heat flux amounts to 238 W m −2 .

Appendix 2: CO 2 Concentrations for Case 4
The following calculations lead to the values of the CO 2 mixing ratio (r c ) and CO 2 molar fraction (χ c ) for CASE 4, referring to Fig. 1.
For the dry air on the left side, with no water vapour present, the mass fraction and mixing ratio are equal (r c = 600 mg kg −1 ), and the molar fraction can be determined simply by a conversion of units For the moist air on the right side of Fig. 1, it is necessary to recall the relationship between the specific gas constants for moist air (R m ) and dry air (R d ), as a function of the mixing ratio (r v ), as (Rogers and Yau 1989), which derives from their molecular masses (M m and M d ) For r v = 10 g kg −1 , and using M d = 29 g mol −1 , this yields M m = 28.8 g mol −1 . The mass of moist air is the sum of the mass of dry air and water vapour and when using the definition of the water vapour mixing ratio this becomes Solving (17) for the dry air mass and substituting this into the definition of the CO 2 mixing ratio yields