A Functional Approach to Vertical Turbulent Transport of Scalars in the Atmospheric Surface Layer

Eddy covariance has been the de facto method of analyzing scalar turbulent transport data. To refine the information available from these data, we derive a simplified version of the turbulent scalar-transport equation for the surface layer, which employs a more explicit form of signal decomposition and dispenses with Reynolds averaging in favour of an averaging operator based on the relevant scalar-flux driving variables. The resulting method, termed functional covariance, provides five areas of improvement in flux estimation: (i) Better representation of surface fluxes through closer correspondence of turbulent exchange with variations in the driving variables. (ii) An approximate 25% reduction in flux uncertainty resulting from improved independence of turbulent-flux samples. (iii) Improved data retention through less onerous quality control (stationarity) testing. (iv) Improved estimation of low-frequency flux contributions through reduced uncertainty and avoidance of driving-variable nonstationarity. (v) Potential elimination of flux-storage estimation when state driving-variables are used to define the functional-covariance flux averaging. We describe the important considerations required for application of functional covariance, apply both functional- and eddy-covariance methods to an example dataset, compare the resulting eddy- and functional-covariance estimates, and demonstrate the aforementioned benefits of functional covariance.


Introduction
The routine observation of energy and mass fluxes between the land surface and the lower atmosphere has undoubtedly been transformed by the widespread adoption of the eddycovariance technique, which has enabled the measurement of near-continuous fluxes relevant to surface areas of tens of thousands of square metres. In combination with other biophysical data, this has lead to greatly improved models of ecosystem functioning (Stockli et al. 2008;Williams et al. 2009;Groenendijk et al. 2011;Melaas et al. 2013;Stoy et al. 2013). Yet, it is also accepted that unresolved issues remain with eddy-covariance flux estimations, such as non-closure of the energy budget, inadequate accounting of low-frequency flux contributions, data losses resulting from quality-control tests, instrumental issues, and uncertainties stemming from flux-storage and advection terms. These unresolved issues complicate the interpretation of surface fluxes obtained using the eddy-covariance technique, so that minimizing their effects would further enhance our understanding of ecosystem functioning.
Of the aforementioned issues, both flux storage and advection relate to distinct terms in the scalar-transport equation, separate from those of turbulent transport. Although the storage term is normally included in the eddy-covariance technique, its measurement is often neglected or approximated (McHugh et al. 2016). Horizontal-advection terms, which are excluded from eddy covariance based on the often questionable assumption of horizontal homogeneity, are rarely measured because of the considerable effort required (Feigenwinter et al. 2004;Aubinet et al. 2005), and results have generally been inconclusive or difficult to interpret (Marcolla et al. 2014). An improved understanding of these neglected terms beyond simple parametrizations, such as the observed friction velocity u * relation for flux losses (Goulden et al. 1996;Aubinet 2008), requires the improved accuracy and precision of turbulent-flux estimations.
Some of the aforementioned issues (energy-budget closure, low-frequency fluxes, uncertainty, data losses as a result of quality control, and storage) may arise from the temporal averaging implicit in the eddy-covariance technique. For example, typical energy-budget closures of 80% have been shown to improve with increasing flux-averaging period (Mauder et al. 2006), while stationarity tests (e.g., Foken and Wichura 1996;Vickers and Mahrt 1997;Bendat and Piersol 2000) can result in the elimination of 30% of runs (Vecenaj and De Wekker 2015), and storage terms generally decrease in magnitude with increasing flux-averaging period. Other issues, such as instrumental errors (Gash and Dolman 2003;Kochendorfer et al. 2012;Horst et al. 2015;Frank et al. 2016), and the development of more precise models of flux responses to changing environmental conditions, may benefit from analyses less tied to the rigid sequential time averaging of eddy covariance.
We argue here that the temporal averaging employed by eddy covariance is a suboptimal approach for flux estimation, and, as an alternative, propose an averaging operator based on the driving variables appropriate to the scalar flux under consideration. The derivation of our methodology, referred to as functional covariance, is similar to the derivation of the eddy-covariance technique, but avoids the implicit link between the filtering and averaging operators with the standard Reynolds decomposition and averaging. Our derivation is followed by a description of the steps required to implement functional covariance, which we then apply to a sample dataset, and compare the results with fluxes determined using standard eddy-covariance methods.

Theory
In the scalar-transport equation, the source/sink term S c of some scalar ρ b c is expressed as the sum of the time rate of change of the scalar, the three-dimensional advective transport, and three-dimensional molecular diffusion (Foken et al. 2012a) where c is a scalar, which may be a mixing ratio or the enthalpy of air, and u is the threedimensional velocity vector. The K term represents the scalar molecular diffusivity, and the ρ b term is the density of the static, background, atmospheric composition. Justification for the use of the density ρ b instead of the dry air density is given in Appendix 3. The standard approach in converting Eq. 1 into a form more tractable for experimental measurement employs Reynolds decomposition, which separates variables into their true mean and deviations from that mean (i.e. the fluctuating component). Nominally, this separation is based on an infinite ensemble of measurements (Kaimal and Finnigan 1994) such that the true mean is the ensemble mean a , and the prime term a † represents deviations from that ensemble mean, a = a + a † . (2) As ensemble measurements are restricted to theoretical or laboratory environments, we replace the infinite ensemble of measurements with an infinite time series of measurements, and employ the ergodic hypothesis to assert that turbulence advected past the measurement location can be treated as an ensemble of measurements.
In practice, infinite time series are impossible, with the available data limited by the experiment length or duration. A further temporal restriction of experimental data involves breaking the data into, typically, 30-min periods, upon which Reynolds decompositon is performed under the assumption of stationarity. Such treatment of the time series results in a more complex representation than that suggested by Eq. 2. It is common to have, in addition to the implicit infinite average value a , a mean associated with a measured dataset of finite lengthȃ (i.e. the experiment mean), a low-frequency fluctuating componentã (i.e. the flux-averaging-period mean), and a high-frequency fluctuating component a , a = a +ȃ + a + a . ( Equation 3 represents a typical time-series decomposition implicit in eddy-covariance studies, though a more complex decomposition would be required if additional filtering were to be applied. For the standard eddy-covariance approach, the low-frequency fluctuating component a often takes the form of a block average (i.e. a). It is worth noting that only the a term on the right-hand side (r.h.s.) of Eq. 3 is a constant function of time. The other terms are time-varying components and thus do not, in the strictest sense, qualify as mean terms under the Reynolds rules of averaging. For a single experiment, we may simplify Eq. 3 by combining the a andȃ terms on the basis that we do not have an infinite dataset with which to evaluate the first two terms on the r.h.s. of Eq 3, viz a = a +ȃ. (4) Assuming our experimental mean a is a small deviation from the infinite mean a , (i.e. a ≈ 0), we employ a as our Reynolds mean value, which is an assumption that obviously improves with increasing experiment length. The resulting decomposition implicit in turbulent exchange field studies then becomes

123
For convenience in our derivation, we temporarily replace the two temporally fluctuating quantities in Eq. 5 with a single fluctuating quantity, allowing Eq. 5 to be expressed as We next define our averaging operator. Under Reynolds averaging, the fluctuating terms are, by definition, deviations from the Reynolds mean. For standard eddy-covariance approaches, this results in the relation a = 0, with the equivalent relation for Eqs. 5 and 6 being a * = 0. Reynolds averaging thus provides the convenient relation a b = 0, which, along with the assumption of stationarity that negates low-frequency flux contributions in the product of means ab = 0, enables the great simplification seen in the eddy-covariance equation.
Given the temporal averaging imposed by standard eddy-covariance techniques and the fact that true stationarity is rare, we have chosen to abandon Reynolds averaging, and instead employ an averaging operator that better isolates the various states of the typically nonstationary environmental drivers controlling the scalar flux. While the specifics of our averaging operator, designated by a dashed overbar, are purposefully left undefined, it can be considered as representing an average over some non-overlapping data grouping, of which the sequential time grouping employed by eddy covariance is one option.
To begin the derivation, we use the simplified version of Eq. 7 to expand the time rate-ofchange term of the continuity equation which, upon application of our averaging operator, enables us to remove the time rate-ofchange term of the experiment-mean background static density, giving We next apply Eq. 7 to the scalar-transport equation (Eq. 1) for which the diffusion term has been assumed negligible (Leuning 2004), and apply the decomposition only to the time rateof-change term and the gradient term's scalar variable. Upon removal of the rate-of-change of constant values and application of our averaging operator, then Using Eq. 9, we remove the third and fourth terms on the r.h.s. of Eq. 10. Following Leuning (2004), we also remove the first term on the r.h.s. of Eq. 10 under the assumption that the temperature-and pressure-induced density fluctuations do not correlate with fluctuations in the scalar mixing ratio, giving If we expand the density and velocity in the advection term in Eq. 11, and make the assumption of horizontal homogeneity, the resulting equation becomes We then remove components containing the experimental mean vertical velocity from the second term on the r.h.s. of Eq. 12, giving which assumes that the experimental mean vertical velocity component has been forced to zero through coordinate rotation.
In standard eddy-covariance derivations, the second term in the vertical partial derivative of Eq. 13 is assumed negligible and is removed. However, for the signal decomposition employed in our derivation, the fluctuating variables contain both high-and low-frequency fluctuating components, so that for the remainder of this derivation we assume this term is negligible.
Upon replacement of our simplification of the fluctuating components (Eq. 6), we express Eq. 13 as which we integrate vertically from the surface to the measurement height h. For clarity, the nominal horizontal integration commonly presented in eddy-covariance derivations is not applied. Under the assumption of horizontal homogeneity, horizontal integration represents an inconsequential operation. Assuming no sources or sinks within the atmospheric column above the surface (S c,z>0 = 0), the equation becomes Equation 15 defines the surface flux as a function of the vertical turbulent transport of scalars with no assumption about the form of filtering used to determine the low-frequency fluctuations ( ) or the form of the averaging operator ( ). We refer to this computational approach as "functional covariance" to distinguish it from the eddy-covariance method. The rationale for this description becomes apparent via its application described below. We note that the form of Eq. 15 is similar to that derived via eddy-covariance approaches despite the non-specific nature of the filtering and averaging operators. This is a characteristic previously noted by Germano (1992) when filtering the Navier-Stokes equations in the context of largeeddy simulation.
An eddy-covariance form of scalar vertical turbulent transport can be obtained from the functional-covariance equation by applying a sequential block-average filter to determine the low-frequency fluctuations ( = ) with corresponding sequential block averaging ( = ) of the fluctuating components. With this filter and averaging, Reynolds averaging rules may be applied, allowing elimination of the interaction covariances (wc and w c), so that Eq. 15 becomes This version of the eddy-covariance equation differs from standard presentations on three points. First, the density term represents the experimental mean instead of the flux-averagingperiod mean. This difference stems from assumptions on the true mean of the Reynolds decomposition, and is likely to have, at most, a few percent effect on the resultant fluxes. Second, the 'storage' term is the averaging-period mean of the time rate-of-change of the highfrequency fluctuating component and not the time rate-of-change of the averaging-period mean (Foken et al. 2012a). This difference appears to stem from an incorrect application, in some eddy-covariance derivations, of Reynolds averaging rules to the time rate-of-change term under the ergodic hypothesis (see, Finnigan 2006). The third difference is the inclusion of a term representing the product of means w c, which is often excluded from eddy-covariance equations because the low-frequency fluxes it represents are not of interest, or because it is (often incorrectly) assumed zero or erroneous.
One final comment is required regarding the use of fluxes determined by the functionalcovariance method (Eq. 15). Because this method does not predefine the averaging operator, the possibility exists for an unequal number of samples to contribute to different averaging realizations. As a result, any subsequent averaging of the functional-covariance fluxes must account for this disparity in sampling representation. As an example, the average flux for the entire experimental period is determined from individual functional covariances using the relationship where S c is the average flux for the entire dataset, S c is a single functional-covariance flux estimate, K is the number of functional-covariance flux estimates, k i is the number of samples in the i th flux estimate, and k is the average number of samples in all flux estimates. For eddy covariance, in which the number of samples in each flux estimate is identical, a straightforward averaging of the fluxes is possible.

Methods
Here we describe the data and methods needed to implement the functional covariance, with details of the computational methods presented in Appendix 1.

Data Requirements
The data required for calculating functional-covariance fluxes are the same as those required for eddy covariance (e.g. the parameters u, v, w, T , and CO 2 concentration and H 2 O concentration, all measured at 10 Hz), with the additional requirement that variables that control the rate of scalar fluxes of interest be measured at a rate that captures their temporal rate-ofchange. For the fluxes of energy and CO 2 , this can be largely achieved by recording radiation parameters at 10 Hz, while other scalar fluxes may depend more strongly on other controlling variables. The selection of variables is discussed in Sect. 3.5.

Dataset
Data are taken from the Griffin experiment site (Clement et al. 2012), and were taken from one of two flux systems mounted at 38 m above a Sitka spruce plantation with a canopy height of 25 m. The instrumentation consisted of a sonic anemometer/thermometer (Gill,Solent R3), an enclosed open-path infrared gas analyzer (Licor, LI7500) (Clement et al. 2009) and a photosynthetically active radiation (PAR) sensor (Macam, SD101QV). Data were collected at 10 Hz and recorded using a Campbell Scientific CR3000 datalogger. The response time of this and similar silicon radiation sensors have been shown to be sufficiently fast for the purpose of this analysis (Clement 2004;Suehrcke et al. 1990).

Separation of High-and Low-Frequency Fluctuations
The primary analysis makes use of the sequence of 30-min mean values a to separate the low-and high-frequency fluctuating quantities ( a and a ) from the time series a, This filter (a) was chosen purely for reasons of convenience, and other filters with different time constants or shapes may be considered.
To assess the effect of increased low-frequency contributions to the turbulent exchange, we replaced a with values of the 30-min means that are low-pass filtered using a convolved Hann window with the window sizes of 3, 6, 12, 24, 48, 72 or 96 h. These low-pass-filtered means were used in lieu of the 30-min mean in Eq. 18 to separate the high-and low-frequency fluctuating time series.

Averaging Operator
In contrast to eddy covariance, which inherently employs sequential time averaging, functional covariance must specify an averaging operator. We suggest an averaging operator based on environmental drivers relevant to the flux being examined, and expressed as a conditional average of some variable a dependent on the driving variables (b, c, . . .), Using this notation, the variable a represents one of the covariance terms from Eq. 15, while b, c, etc. represent the environmental variables controlling the flux associated with the variable a. For each functional-covariance average, we must also specify the corresponding subrange of each environmental variable (e.g. B i to B i+1 ). The details of selecting the appropriate dependent variables, and their ranges, are discussed further in the following sections. It should be noted that eddy-covariance time averaging is also a conditional average using a form of Eq. 19 in which the only dependent variable is time, The difference between functional and eddy covariances can be visualized graphically in Fig. 1, where the eddy covariance can be considered as a one-dimensional array of cells, while functional covariance can be seen as a multi-dimensional array of cells, with each cell representing a flux estimate based on the subrange of the corresponding dependent-variable dimension.

Variable Selection for Defining the Averaging Operator
The variables defining the averaging operator (Eq. 19) comprises the incoming radiation flux R (PAR flux was employed in this analysis), air temperature T , vapour pressure deficit Because Griffin is a monoculture coniferous forest in a maritime environment that has reached canopy closure, we have not included the leaf area and soil moisture, as both are relatively constant and, thus, likely minor determinants of variations in flux values. However, as has been suggested previously (Clement et al. 2012), horizontal advection associated with katabatic and anabatic flows, or pressure-driven advection (Finnigan and Belcher 2004), are believed to significantly affect the flux rates at this site. The drivers of these advective effects are probably a complex combination of topography, canopy structure, and atmospheric flow conditions. Because the definitive controlling variables of advection are not well understood, we have selected the simpler proxy variables of wind speed and wind direction in an attempt to account for the various possible advective effects on the fluxes.
In general, proxy variables should be avoided when using functional covariance because the flux response to the proxy variable may differ from its response to the corresponding true driver. However, proxy drivers may be required until knowledge of the true drivers and associated measurements are available. By comparison, eddy covariance uses time as a proxy for all driving variables, and subsequently attempts to sort out the cause-and-effect relationship between fluxes and drivers using both temporally-averaged fluxes and drivers under the assumption of stationarity for both.
Note that if a significant driving variable is not included in the functional-covariance averaging operator, the neglected variable may obfuscate the responses of known variables. It is, therefore, necessary for users of this approach to have an understanding of, or at least a good hypothesis of, the variables driving the flux. Similar misattribution of flux effects may also occur in eddy-covariance post-analysis if the correct driving variables are not examined.
While knowledge of the driving variables may seem a weakness in functional covariance, an ex post facto comparison of functional-and eddy-covariance fluxes may be used to iden-tify unknown driving variables, as shown below. Obviously, such data exploration is labour intensive, making it preferable to establish the driving variables a priori.

Averaging-Operator Granularity
In the standard eddy-covariance method, a temporal-averaging period of 30 min has been established as capturing most scales of turbulent transport under the assumption of stationarity (Foken and Wichura 1996). Lenschow et al. (1994) show (their Eqs. 50 and 59) that eddycovariance measurements using sampling intervals much greater than the integral time scale (i.e. disjunct sampling) have a random uncertainty that decreases with the inverse of the square root of the number of independent samples (1/ √ N ). For contiguous, high-frequency samples (e.g. eddy covariance), the decrease in uncertainty with increasing sample size is slower because of temporal autocorrelation.
Data sampling by the functional-covariance method for a particular flux estimate combines contiguous and widely-separated temporal samples. For an approximate comparison of the error estimates from eddy and functional covariances, we assume that independent samples must be separated by at least 10 integral time scales, and that the integral time scale is approximately 1 s. Under these assumptions, a typical 30-min eddy-covariance averaging period sampled at 10 Hz contains 180 independent samples, with an associated random uncertainty of 7.5%. To estimate the number of independent samples in a functional-covariance flux, we examined the data contributing to a turbulent flux for a single set of averaging-operator parameter ranges (PAR flux between 300 and 350 µmol m −2 s −1 , air temperature between 8 and 12 • C, wind speed between 1 and 3 m s −1 , wind direction between 180 • and 240 • , and vapour pressure deficit between 0.5 and 1.0 kPa). For the Griffin experiment dataset, 162,730 samples (i.e. of 10-Hz data) meet these parameter criteria, and are obtained from 358 different 30-min periods spanning the two-year dataset; an examination indicated there are 3227 independent samples. As this subset of data is approximately 10 times larger than a standard 30-min eddy-covariance dataset (i.e. 18,000 samples), we estimate a functional-covariance dataset of corresponding size having 323 independent samples. This independent sample size corresponds to a random-error estimate of 5.5%, suggesting that the random-flux uncertainty of the functional-covariance method is 25% more efficient than that of the eddy-covariance method.
Importantly, more data may be added to a functional-covariance flux estimate to reduce its random uncertainty, while doing the same for an eddy-covariance flux entails increasing the averaging period, with the concomitant potential for violating the assumption of stationarity. Note also that a finer granularity in the averaging operator applied to the driving variables would likely increase the number of widely temporally-separated samples, thereby reducing uncertainty-though this improvement comes at the expense of fewer samples in each average flux. This suggests there may be an optimum in the subrange granularity of averaging-operator driving variables.
The preceding analysis provides some guidance in defining the averaging-operator bin sizes required by the functional-covariance method. We know that 400 independent samples are required to provide a flux estimate with a random uncertainty of 5%. The preceding calculations show that approximately 50 samples of 10-Hz data are required to obtain one independent sample for functional covariance. Under these assumptions, averaging bins of > 20,000 data samples are required to estimate fluxes to within 5% random uncertainty.
Unfortunately, driving variables used to define the functional-covariance averaging operator are unlikely to be uniformly distributed over their ranges. Unless elaborate binning schemes are adopted, the number of samples contributing to each functional-covariance flux may be considerably larger or smaller than the 20,000 samples desired for a 5% random uncertainty. However, if we optimistically assume uniform sample distributions over the various driving-variable dimensions, we can use the desired 20,000 samples per bin to identify binning strategies.
To begin with, we note that for one year of 10-Hz data, there are sufficient instantaneous samples for approximately 16,000 averaging bins, with each containing the desired 20,000 raw data samples required to achieve the desired 5% random uncertainty. If we assign seven binning ranges for each of our five averaging-operator variables (PAR flux, air temperature, vapour pressure deficit, wind direction, and wind speed) we would have 16,807 averaging bins. Obviously, some variables have a greater influence on fluxes, while regions of some of the variable ranges must be more closely sampled because of the more rapid response of the flux to changes in that variable.
The variable binning boundaries employed in our analysis are listed in Table 1, which produces 17,010 bins, exceeding the recommended number of bins for one year of uniformlydistributed data, but sufficient for the two years of data being analyzed. Although the binning strategy (Table 1) employed for our example dataset appears appropriate, it highlights the problem of adding more averaging-operator variables or having a finer binning granularity. Both increasing the number of driving variables or employing finer granularity result in large increases in the number of averaging bins, and a corresponding increase in flux uncertainty. For eddy covariance, the equivalent problem corresponds to employing shorter averaging periods.
Complex systems, such as towers located in regions with great spatial or temporal heterogeneity or airborne turbulent-flux measurements of multiple ecosystems, likely require many more driving variables (Metzger et al. 2013), and possibly a finer granularity than that used here. Such measurement systems requiring a large number of averaging bins require exceedingly large datasets to provide reasonably small flux uncertainties over the range of all variables employed. Note that the use of eddy covariance instead of functional covariance under such circumstances does not address this problem, as an even greater amount of data would be required to obtain the equivalent uncertainties, as described above.

Data Alignment
When using eddy covariance, one important data treatment is the temporal alignment of the vertical velocity component and scalar time series, which is commonly referred to as lag removal (Moncrieff et al. 1997;Rebmann et al. 2012). Because sonic-anemometer wind speeds and temperatures are nominally sampled coincidentally, temporal alignment is usually applied only to scalars sampled by other sensors, particularly trace gases sampled with closed-path analyzers. The functional-covariance approach additionally requires the temporal alignment of the driving variables used by the averaging operator. Firstly, it bears stating a conceptual example of what we believe is happening.
For a flux tower above a horizontal uniform surface, a step change of radiation flux in the flux footprint at time t 0 induces rapid changes in some driving variables (e.g. surface and air temperatures and vapour pressure deficit) and minimal changes in other driving variables (e.g. soil moisture, leaf area), with a corresponding rapid change in the scalar fluxes (Pearcy 1990;Sassenrath-Cole and Pearcy 1994;O'Sullivan et al. 2013;Sun et al. 2016). Assuming the wind direction at cloud height is roughly aligned with that near the surface, then both the radiation field and the near-surface changes in atmospheric properties would be advected towards the flux tower. However, the radiation-flux field advects more rapidly than the near-surface atmospheric properties because of differences in wind speed at cloud level and at the surface, and because of the rate of turbulent diffusion from the surface to the measurement height. Therefore, the near-surface atmospheric fields would lag the radiation fields for measurements made on the tower. For example, for our dataset, flux measurements are made 20 m above the zero-plane displacement, and we assume an average flux-footprint distance of 100 m, a friction velocity of 0.3 m s −1 , a cloud speed of 15 m s −1 , and an average wind speed of 2 m s −1 between the surface and the height of measurement. Therefore, the radiation-flux field would take 7 s to reach the tower, but the turbulence field would take a longer 83 s to reach the tower, giving an approximate 76-s lag of the turbulence field behind the radiation-flux field as measured at the tower.
While analysis can identify temporal lags between radiation and its influence on nearsurface atmospheric properties in the flux footprint as observed at the flux tower, there are many drivers not advected with the flow field, such as leaf and soil temperature, stomatal conductance, soil moisture, leaf area, canopy structure and topography. Some of these variables (soil temperature, soil moisture, leaf area, canopy structure and topography) typically change so slowly that 30-min average values, intermittent measurements, or even one-off characterizations may be sufficient for their representation. Other variables (leaf and soilsurface temperature, stomatal conductance) may change rapidly, and we must either correlate their effect with atmospheric properties, or devise cunning plans to measure their real-time, in situ, responses. There may also be intermittent events that produce rapid changes in normally slowly-changing environmental drivers, such as rain events affecting soil moisture and temperature. Such events would require more frequent sampling of variables that would nominally be sampled slowly.
In our analysis, we assume driving variables, other than radiation, are either well represented by atmospheric properties, or that spatially-variable, but temporally-static drivers can be characterized using the wind direction as a proxy variable. We realize that the topic of data alignment of driving variables and fluxes is only briefly discussed, and suggest more in-depth research of its effect on functional-and eddy-covariance fluxes.

Method Application
The calculations of functional and eddy covariances and the resampling of functional covariances as time-series fluxes, as described in Appendix 2, were performed using EdiRe version 1.5.0.75. Subsequent analysis of the results used Python version 3.5.2.

Results and Discussion
We first examine the magnitude of nonstationarity of radiation flux and then describe the effects of data alignment, aggregation and low-frequency contributions on functionalcovariance fluxes. We then present a comparison between eddy-covariance and functionalcovariance fluxes as both time series and functional responses, followed by comparisons of uncertainty estimates and low-frequency flux components.

Stationarity of Drivers
Much of the rationale for developing functional covariance arises from stationarity issues. Typically, nonstationarity is considered to cause a misrepresentation of fluxes because of neglected low-frequency turbulence contributions. While this problem may be solved by increasing eddy-covariance averaging periods to include these low-frequency contributions, the more common approach is to eliminate fluxes from nonstationary averaging periods.
Similarly, nonstationarity of driving variables can result in the misinterpretation of eddycovariance fluxes, which are normally interpreted using the mean environmental drivers from the corresponding averaging period. If these drivers are nonstationary, there will be a misrepresentation of the environmental responses of the eddy-covariance fluxes if the fluxdriver relationship is nonlinear. This problem is especially true for radiation flux, which has strong influences on most scalar fluxes, and can vary greatly over short time periods. Nonstationarity of drivers cannot be resolved by extending the temporal averaging period, as can be done when addressing large-scale stochastic variability in eddy-covariance fluxes, because such an approach would further degrade the relation of eddy-covariance fluxes to the driving variable.
To demonstrate this issue, in Fig. 2 we present radiation-flux nonstationarity plotted as the mean probability distribution of (10-Hz) incoming PAR flux in relation to the corresponding 30-min-mean incoming PAR flux. For moderate to high levels of 30-min-mean radiation flux, Enhancement of fluxes in response to lag removal of high-frequency data relative to high-frequencyradiation measurements. The curves show the percentage change of fluxes at different lag times for moderate (left) and high (right) radiation levels the actual distribution of PAR flux varies significantly from the mean value. For mean PAR flux 1200 µmol m −2 s −1 , the actual distribution of PAR flux is strongly skewed to lower levels while at higher mean PAR flux it is positively skewed. It is notable that, on average, ≈ 10% of radiation-flux data at moderate to high levels of mean radiation flux differ from that mean by more than ≈ 500 µmol m −2 s −1 .
Despite the fast responses of photosynthesis to radiation changes (Pearcy 1990;Sassenrath-Cole and Pearcy 1994;Sun et al. 2016), the variation in high-frequency radiation is generally not considered when monitoring ecosystem CO 2 fluxes, but is more common in solar-energy studies, where a similar behaviour to that in Fig. 2 has been observed (Soubdhan et al. 2009). The effect of radiation-flux nonstationarity on the CO 2 flux, F c , is discussed in Sect. 4.6.
By employing functional covariance in which instantaneous data are grouped by similar environmental conditions, problems caused by nonstationarity of environmental drivers are minimized, and the nonstationarity of turbulence in the usual sense ceases to be relevant. However, the general concept of stationarity is still valid in that turbulent statistics, under a specific set of environmental conditions, should remain time independent, unless there is some environmental factor forcing a change in turbulence characteristics not captured by the driving variables employed in the functional-covariance calculation. In this situation, the variable causing that temporal variation in turbulence characteristics should be incorporated into the functional-covariance calculation.

Data Alignment
The lag in transit time of atmospheric properties relative to the radiation field during advection between the footprint and tower appears to not have been quantified previously, though the phenomenon is intuitively straightforward. By removing this lag, functional-covariance fluxes become more closely associated with their driving variables (e.g. radiation) as observed in Fig. 3, which shows the relative changes in functional-covariance fluxes for moderate and high radiation-flux conditions, as a function of the lag time. Figure 3 shows that a lag time of about 90 s gives the maximum enhancement in the response of the functional-covariance flux to varying radiation; this lag corresponds closely to the 'back-of-the-envelope' calculations presented in Sect. 3.7. Closer examination reveals that the optimum lag is different for the four fluxes examined, with the optimum lags for the parameters F c , H , u w and λE being 80, 98, 140, and 165 s, respectively. The observed differences in the lag times likely reflect differences in the vertical source distribution through the canopy (Ross and Harman 2015), with CO 2 uptake greatest at the top of the canopy, and water-vapour sources strongest deeper in the canopy. The longer lag times for fluxes originating lower in the canopy (e.g. for latent heat λE) will be a combination of the increased time taken for radiation-flux changes to affect changes deep in the canopy plus the increased time taken for the affected scalar to traverse the greater distance to the flux tower. For short canopies with little vertical separation of flux sources, it is predicted that this lag time is similar for all fluxes. Subsequent analysis of our example dataset assumes a constant 90-s lag of all atmospheric fields behind the radiation-flux field. A more thorough analysis would impose a flux-specific lag for all variables upon which a flux depends.
While the observed lag of fluxes relative to radiation-flux is important for the functionalcovariance method, it may be less so for standard eddy-covariance approaches-unless there are anomalous radiation-flux changes in the first/last 1-2 min of a run. For short eddycovariance averaging periods (e.g. 5 to 15 min), any discrepancy caused by not accounting for this lag may result in improper relationships between radiation flux and fluxes. If longer lag times exist for other driving variables, these would likely affect both eddy and functional covariances. Although not investigated, we suggest that the observed differences in lag times may have important ramifications for the assumption of similarity of the transport of scalars over tall vegetation.
We also emphasize that our treatment of data alignment deals only with the radiation flux, and that there may be lag effects resulting from events associated with other drivers not captured by our analysis. We suggest this topic requires further detailed research.

Averaging Contributions
One of the attractions of functional covariance is the possibility of accumulating more data under a particular averaging condition than would normally be included in a standard eddy-covariance averaging period. However, in defining our averaging-operator bins, we optimistically assumed the data would be uniformly distributed over the range of each driving variable. Examination of the actual data accumulated in the averaging bins shows that only 22% of the bins accumulated more 10-Hz samples than the 18,000 samples used to evaluate a standard eddy-covariance flux. Indeed, approximately 50% of the bins contained < 1800 raw data samples (i.e. < 3 min of data).
There are several ways of interpreting this observation. A pessimistic view dismisses functional covariance as inadequate for obtaining precise flux estimates for driving-variable combinations having low sample counts. This view can be countered by the fact that these poorly-represented driving-variable combinations are completely unrepresented by standard eddy-covariance approaches. Although the flux uncertainty may be high, at least the functional-covariance approach provides an estimate of fluxes for these niche conditions, which would likely be obfuscated by temporal averaging under the standard eddy-covariance procedure.
A more optimistic view seizes upon this information to identify the driving-variable conditions lacking turbulent-transport data, and either uses that knowledge to refine the analysis,  Table 1) better design experiments, or predict the duration of a field experiment required to obtain sufficient samples for a more robust flux estimate. To this end, we have examined the fraction of driving-variable averaging bins containing < 5000 data samples (corresponding to ≈ 10% flux uncertainty) and > 20,000 data samples (5% uncertainty) as a function of the parameter ranges specified in Table 1. Apparent from the analysis shown in Fig. 4 is that the majority of flux estimates have fewer than the 5000 samples needed to achieve a 10% level of uncertainty. It is interesting to note that, for this example dataset, the poorest representation of the fluxes occurs at a high radiation flux, low air temperature, high vapour pressure deficit, and low wind speed. While the distribution of representation by wind-direction is quite uniform, it should be noted that the wind-direction averaging bins were selected to optimize representation, while other bins were more uniformly distributed to well capture the nature of the associated environmental responses.
An important conclusion from this analysis is that, even for this relatively simple ecosystem, more than eight years of data may be required to characterize the majority of functional-covariance flux estimates to within 5%. Coarsening the granularity of some of the parameter ranges (e.g, the PAR flux) may reduce data requirements, though finer granularity in other parameters (i.e. vapour pressure deficit) may offset those gains. Additionally, experiments lasting as long as eight years are likely to require characterization of other flux dependencies (e.g. soil moisture, leaf area). A further option for reducing the uncertainty without extending the length of experiments is to provide parallel sampling (Hill et al. 2016) on the assumption that driving variables and ecosystem characteristics are the same when It seems that the characterization of the environmental responses of functional-covariance fluxes inevitably requires compromises on accuracy or quality, or increased effort or cost if a low uncertainty is desired for the majority of environmental driving-variable conditions. While this conclusion seems pessimistic, it is noted that eddy-covariance fluxes likely require significantly more data to achieve the same accuracy obtained by functional covariance when assessing fluxes for the same range of environmental conditions.

Influence of High-and Low-Frequency Fluctuation-Separation Filters
Increasing the low-frequency content in the turbulent exchange by modifying the filter separating high-and low-frequency fluctuations (Eq. 18) results in substantial changes to flux estimates. The percentage change in diurnal and nocturnal functional-covariance fluxes in relation to the filter length is shown in Fig. 5, illustrating the substantial increases in diurnal values of H and F c with increasingly lower-frequency contributions. (Note that the filterwindow length corresponds to flux contributions out to half that period, e.g. a 48-h window length roughly corresponds to low-frequency flux contributions up to 24 h). The diurnal values of the fluxes λE and u w show slight increases (2%) for low-frequency contributions up to a 12-h filter-window length, but exhibit reductions for lower-frequency contributions.
The results of Foken (2008) show a qualitatively similar, near-linear increase in H and non-linear increase/decrease in λE with increasingly lower-frequency flux contributions. Quantitatively, however, the effect observed by Foken (2008) is substantially larger; a 90% increase in the value of H compared with our 25% increase. Earlier studies (Sakai et al. 2001;Finnigan et al. 2003;Mauder et al. 2006) have also noted the potential contributions of low-frequency turbulent transport.
Current thought is that the low-frequency flux contributions are associated with circulations tied to landscape heterogeneity (Foken 2008;Stoy et al. 2013). While such spatial Fig. 6 Effect of increasing the low-frequency content of sensible heat flux H in relation to the driving variables of the functional-covariance flux. Each curve represents the average percentage increase in H calculated using the specified high-pass-filter window size relative to H calculated using a 30-min high-pass-filter window size variability may plausibly be associated with low-frequency fluxes, a more refined understanding of the behaviour of low-frequency fluxes is required to identify the controlling variables that predict when and where such low-frequency contributions occur. A proper analysis of low-frequency contributions likely requires an understanding of, and accounting for, the potential driving variables at distances currently considered far from the tower location (Malhi et al. 2004).
Functional-covariance fluxes also enable us to examine the low-frequency flux contributions in relation to the driving variables (Table 1). While unlikely to be optimal for predicting low-frequency fluxes, these variables enable insight into the behaviour of these low-frequency fluxes. Figure 6 shows the average percentage changes in diurnal functional-covariance sensible heat fluxes, for a range of low-frequency filters, plotted against the driving variables used in their determination. The low-frequency flux contributions do not change substantially with radiation flux, and although no relationship of low-frequency flux with air temperature or vapour pressure deficit would be expected, the functional-covariance results show nonlinear dependencies. Further analysis is required to determine if these are real or artifacts of the available data, its presentation, or unknown correlations with a neglected driving variable. If low-frequency flux contributions are tied to landscape characteristics, it is more likely that they may be related to the wind direction and wind speed. Figure 6 suggests that lowfrequency contributions increase the value of H at low to moderate wind speeds, but reduce fluxes at high and very low wind speeds. The directional dependence of low-frequency contributions is complex, but reflects topographically-induced flow effects for this site (Clement et al. 2012). Although the land cover at the Griffin site is not ideally uniform, it compares favourably to some of the agricultural examples shown in Foken (2008), suggesting that L L Fig. 7 Comparison of a sample of functional-covariance (FC) and eddy-covariance (EC) time series of turbulent fluxes. High-frequency turbulent fluxes determined using 30-min and 24-h filter windows are presented in the left and right columns, respectively low-frequency contributions at the Griffin site may be more dependent upon topography than land-cover variability.

Comparison of Functional-and Eddy-Covariance Fluxes as Time Series
Functional-covariance fluxes calculated using 30-min and 24-h filter windows are converted to time series (see Eq. 29) to compare with the time series of eddy-covariance fluxes calculated using the same high-frequency turbulent fluctuations (Eq. 5), with a sample presented in Fig. 7.
While the temporal patterns are similar for functional-and eddy-covariance fluxes, closer examination shows discrepancies for individual 30-min-average fluxes, and also more persistent discrepancies lasting several hours. Therefore, the immediate question 'which flux is correct?' arises, which, as usual, may be answered with 'it depends'. If all relevant driving variables are well represented in the formulation of the functional-covariance fluxes, then any differences between functional-and eddy-covariance fluxes will be associated with stochastic-atmospheric or random instrumental variability, and the functional-covariance fluxes will be more appropriate. However, if a significant driving-variable is not incorporated into the functional-covariance calculations, eddy-covariance fluxes will be more accurate during periods when the neglected variable has a significant effect on the flux.
We believe that the discrepancies in the fluxes of H and λE observed in Fig. 7 arise because of poor representation of canopy water-vapour conductance in the functional-covariance flux estimates. Examination of the flux time series suggests that discrepancies tend to be greater on days following measurable precipitation. A wetness sensor deployed on-site does not show a strong correspondence with discrepancies between eddy-and functional-covariance fluxes, indicating the need for a better method of independently and continuously measuring canopy wetness or conductance.
An important deduction from Fig. 7 is that discrepancies between the functional-and eddycovariance fluxes provide an opportunity to identify hitherto neglected driving variables, and L L Fig. 8 Comparisons of average diel patterns of the fluxes H , λE, F c and u w from the eddy-covariance (EC) and functional-covariance (FC) methods their period of importance. While this may be a fruitful avenue of investigation, variables with small to moderate effects on fluxes may be difficult to distinguish from background stochastic noise. Returning to Fig. 7, we note that functional-covariance fluxes with a 30-min filter are a smoother version of the time series of eddy-covariance fluxes, which is expected, since more random variability is removed from function-covariance fluxes than from eddy-covariance fluxes. The functional-covariance fluxes with 24-h low-frequency information show a general increase in the flux magnitude, and greater variability relative to 30-min-filtered fluxes. However, this variability is much less than that of the corresponding eddy-covariance fluxes containing corresponding low-frequency information. The variability in functionalcovariance fluxes may represent driving-variable temporal variations, uncertainty caused by insufficient sampling, or incorrect driving-variable specification.
The net effect of using functional-covariance fluxes is revealed in the average diel pattern of fluxes shown in Fig. 8. Functional-covariance fluxes are larger in the morning relative to eddy-covariance fluxes, with a tendency for a relative reduction during midday. This pattern exists for both 30-min and 24-h low-frequency filtering, though the value of H shows a midday enhancement for 24-h filtering.
The differences between the average diel patterns of functional-and eddy-covariance fluxes shown in Fig. 9 enable a better understanding of the observed differences. Also plotted is a model of the flux change-in-storage term developed for this site (Clement 2004). The air-L L Fig. 9 The differences between average diel patterns of functional-covariance (FC) and eddy-covariance (EC) fluxes. Black and grey lines represent functional-covariance fluxes obtained using 30-min and 24-h lowfrequency filters, respectively. The dashed line represents the change-in-storage flux term from Clement et al. (2012) column storage terms included in the model have been increased by a factor of 2.7 to account for the change in canopy height (9 m for the storage models and currently 25 m). The flux differences in Fig. 9 suggest that the functional-covariance method produces fluxes close to the eddy-covariance fluxes adjusted for the change-in-storage term. While the difference in functional-and eddy-covariance fluxes presented in Fig. 9 appears to capture the behaviour of the morning change in storage in the model, the afternoon peak precedes that of the model storage term by several hours and is of a greater magnitude. Note that there appears to be a change-in-storage for momentum fluxes, though this effect is small, being on the order of a few percent. Based on these observations, it is suggested that use of the functional-covariance method inherently accounts for the change-in-storage flux term, so that storage measurements are not required when using this method. Further testing should be done to verify this observation.

Comparison of Functional-and Eddy-Covariance Fluxes as Environmental Responses
Another method of comparing functional-and eddy-covariance fluxes is along one or more of the driving-variable dimensions used in defining the functional-covariance fluxes. For our example dataset, we compare F c values because of the well-known environmental responses for this flux. We begin by comparing the functional-covariance fluxes (0.5-and 24-h filtered) with eddy-covariance fluxes as a function of the PAR flux in Fig. 10, where the light-response curves are shown in the top two graphs. Curves are shown for functional-covariance fluxes with 30-min and 24-h filtering, as well as eddy-covariance fluxes and 30-min filter-window functional-covariance fluxes resampled as run-period fluxes (i.e. Eq. 29). The light-response curves are similar for the four versions of the CO 2 flux, though the magnitude is larger for Although the light-response curves are generally similar, closer inspection shows that the functional-covariance curve is smoother under high radiation-flux, and has a sharper transition at a PAR flux ≈ 750 µmol m −2 s −1 , which more clearly separates the light-limited and substrate-limited CO 2 flux regimes. The bottom left graph shows the percentage difference of the functional-covariance light-response curves relative to the eddy-covariance curve, illustrating that the quantum-use efficiency of the functional-covariance flux exceeds that of the eddy-covariance flux by ≈ 6% at all light levels except for a region between about 900 and 1300 µmol m −2 s −1 .
These observations suggest that the light-response curves of the eddy-covariance flux may be underestimating the light-use efficiency of the CO 2 flux. The clearer separation of the light-and substrate-limited regions of the light-response curves for functional covariance (Fig. 10) suggests that fitting of the light-response models to these curves may provide better estimates of quantum yield and the value of A max . Obviously, further testing of this observation is required to ensure it is not a site-specific characteristic.
The observed differences in functional-and eddy-covariance light-response curves enables interpretation of the contribution of different light levels to the total flux. The bottom right graph of Fig. 10 presents the relative importance of the flux F c as a function of the radiation flux (the relative importance is calculated as the normalized product of radiation-flux probability and the light-response curves), showing a larger carbon uptake at lower and very high radiation flux by the functional-covariance model, which provides a contrasting interpretation of light use in this forest.
The method of calculation of the light-response curves in Fig. 10 should be considered when predicting fluxes using the corresponding models. For both models, the radiation-flux conditions for the prediction period should be the same as those under which the model was developed. The eddy-covariance light-response model was formed using nonstationary radiation-flux conditions, and should only be applied to the radiation flux with similar nonstationarity. Likewise, the functional-covariance light-response curve is formed using (quasi-) stationary radiation-flux data, and, thus, should be applied to stationary radiation-flux values or to instantaneous radiation-flux samples. From this consideration, it is clear that CO 2 fluxes predicted using a functional-covariance light-response curve applied to instantaneous radiation-flux values produce a more accurate estimate than those from typical application of eddy-covariance light-response curves. Figures 11 and 12 provide details of the wind-direction dependency of the F c lightresponse curve. In Fig. 11, the top two graphs show the light-response-curve anomalies as a function of wind direction for functional-covariance (left graph) and eddy-covariance (right graph) fluxes. The anomaly is calculated as the difference between the light-response curve as a function of wind direction and the curve for all wind directions ( Fig. 10 top left). The curves are calculated using 20 radiation-flux values and 5 • wind-direction bins. The bottom two graphs of Fig. 11 show the logarithmic representation of the number of instantaneous samples in each radiation-flux/wind-direction bin. It is noteworthy that the data coverage of the functional-covariance flux is available for all direction and radiation fluxes, while the eddy-covariance fluxes are not present in some of the averaging bins because of the inherent within-run averaging of radiation and wind direction.
Both the functional-and eddy-covariance light-response-curve anomalies (Fig. 11) show a distinct wind-direction pattern with a positive anomaly for wind directions of ≈ 280 • (indicating a minimum in the value of A max ) and a negative anomaly for wind directions of about 200 • . The variability in the anomalies for easterly flow may be associated with the lower amount of data available for these wind directions, implying that more data are needed to characterize these regions.
The effects of the other functional-covariance averaging parameters (air temperature T , wind speed U , and vapour pressure deficit D) on the wind-direction anomalies are shown in Fig. 12, where the anomaly curves are further divided by low and high values of the associated parameters (i.e. T >< 6 • C, U >< 2 m s −1 , D >< 0.5 kPa). Figure 12 shows that the anomalies at wind directions ≈ 200 • and 280 • persist at both high and low values of the other functional-covariance flux-driving variables. A strong air temperature dependence of the parameter A max can be inferred from the general differences in the top two plots, though the wind-direction anomaly can still be observed.
While an in-depth analysis of the cause of a wind-direction dependency in the value of F c is beyond the scope, we simply note that the use of functional covariance appears to provide a clearer picture of the underlying phenomena because of its ability to more completely cover the parameter space being investigated.

Uncertainty Estimates
Although we have discussed the improved flux-uncertainty estimates using functional covariance and have described the method of calculation, the previous results have been presented without error estimates, and so are presented here. For each functional-covariance averaging bin, uncertainty values are calculated as the standard deviation of the instantaneous products of primes. As described in Appendix 2, corresponding eddy-covariance uncertainties are calculated using fluxes grouped by the same driving-variable ranges as used in calculating the functional-covariance fluxes. For the eddycovariance flux uncertainties, it is necessary to adjust the uncertainties so that they represent sample, instead of averaging-period, uncertainties.
Functional-and eddy-covariance-flux uncertainties are only compared when at least three eddy-covariance fluxes are present for the specific combinations of driving-variable parameters. Figure 13 (top left) shows a comparison of the number of qualifying samples in functional-and eddy-covariance-flux averaging bins expressed as the equivalent number of 30-min runs. The deviation from a 1:1 relationship exists because of the averaging of a range of environmental conditions spanning multiple parameter bins within eddy-covariance run periods.
The differences in flux uncertainties (eddy versus functional) normalized by the median functional-covariance uncertainty are plotted as a histogram in Fig. 13 (top right). If the uncertainties from both methods are similar, then the resulting distribution would be centered on zero. Figure 13 shows that the distribution is positively skewed with a median value of 0.29, Fig. 12 CO 2 flux light-response-curve anomalies from eddy-covariance and functional-covariance fluxes for data separated by low (left column) and high (right column) values of air temperature, wind speed and vapour pressure deficit. The colour scale for these plots is the same as for the corresponding plots in Fig. 11 suggesting that eddy-flux uncertainties are approximately 30% larger than those obtained by the functional-covariance method, as predicted in Sect. 3.6.
The bottom two graphs of Fig. 13 show examples of the F c uncertainty determined by the functional-covariance calculations. The bottom left graph shows a time series of the flux F c with an uncertainty expressed as one standard error. Note that the standard error has been increased by a factor of √ 50 to account for the reduction in sample size associated with sample independence, as described in Sect. 3.6.
The bottom right graph shows F c light-response curves for high and low wind speeds for air temperatures between 8 and 12 • C, with the uncertainties representing two standard errors (again adjusted for independent samples). The curves indicate significantly higher uptake under low-wind-speed conditions at moderate PAR-flux levels. This reduction in CO 2 flux at high wind speeds may reflect an increase in stomatal resistance as a result of improved canopy-atmosphere coupling, hinting at the need for a humidity driving variable within the canopy in the functional-covariance calculation.

Low-Frequency Vertical-Flux Terms
Because the functional-covariance approach does not depend upon Reynolds averaging rules, the resulting vertical turbulent-transport equation (Eq. 15) retains three additional terms for low-frequency transport (e.g. the terms w c, wc , and w c). The functional-covariance term w c is identical in nature to the eddy-covariance 'vertical-advection' term w c, with both being descriptors of the contribution of low-frequency fluctuations and not the true mean-advection terms. The interaction terms (w c and wc ) are a consequence of the flexibility of using mismatched averaging and filtering operators. The same terms would be required by eddy covariance if block-average filtering and matched-temporal averaging were not employed.
We have plotted these functional-covariance terms for the fluxes H , F c , λE and uw in Figs. 14 and 15, with each resampled as 30-min averages before forming the diel-average curves presented. For all fluxes (Figs. 14, 15), the primary contribution to the vertical flux is through the high-frequency turbulent-transport term (e.g. w c ). Also, for all fluxes, the lowfrequency vertical-velocity-component/high-frequency scalar-interaction term wc contains almost no mean contribution when averaged by the time of day. In contrast, the low-frequency scalar/high-frequency vertical-velocity component and the low-frequency covariance terms ). Each of the terms is determined as a functional-covariance flux matrix and resampled as sequential 30-min-average values before averaging by hour of the day. The left and right columns contain terms for the sensible heat and CO 2 fluxes, respectively. Each graph contains curves for functional covariances determined using 30-min, 24-h and 48-h filter windows for the determination of high-frequency fluctuations (w c and w c) indicate varying levels of flux contribution over the day. Initial impressions are that the interaction term w c tends to vary about the zero-contribution line, while the low-frequency covariance term w c appears to provide a bias with the same sign as the high-frequency turbulent flux. As seen in the three filter-window curves in each plot, this low-frequency information is mostly transferred from the low-frequency covariance w c to the high-frequency covariance w c when longer filter windows are employed.
Because the filter used to separate the high-and low-frequency fluctuations should simply redistribute the fluxes between the four transport terms of Eq. 15, we expect these terms to sum to the same result regardless of the filter employed. To examine this effect, we plot the sum of the four transport terms for each of the four fluxes as a diel curve in Fig. 16, illustrating that for the fluxes F c , λE, and largely for uw, the total curves are identical or nearly identical, regardless of the filter used to isolate the high-frequency fluctuations. In contrast, the curves for the flux H show distinctly different diel patterns depending upon the filter used to obtain the high-frequency turbulence.
For the flux H , the use of longer filter windows results in an increase in the magnitude of daytime and nocturnal fluxes. The lack of this effect in the fluxes F c and λE, and similar but minor characteristics in the flux uw, suggests that this is not determined by ecosystem energetics, but is likely instead an instrument artifact (Richiardone et al. 2012). The known influence of temperature upon sonic-transducer behaviour leads us to propose that this effect due to correlations of temperature-induced low-frequency contributions in the sonic-anemometer Fig. 15 As Fig. 14, but the left and right columns present terms for latent heat and momentum fluxes, respectively signals w, T and U . What is not clear is whether this low-frequency redistribution results in a more erroneous flux or more correct flux estimate. If the observed H values with lowfrequency contributions are more correct, this behaviour suggests that underestimations of the energy budget may be more attributable to sensible-heat-flux measurements, and not equally or proportionally attributable to both sensible and latent heat fluxes (Twine et al. 2000). This hypothesized effect needs to be examined for other sonic-anemometer models to determine its persistence.

Conclusions
The use of meteorological methods to estimate mass and energy exchange between the surface and the atmosphere has endowed environmental researchers with a powerful set of tools for understanding both natural and anthropogenic ecosystems. While all ecosystem flux-measurement methods have their place, eddy covariance has gained particular favour because of its directness and scales of representation-both temporal and spatial. However, the widespread, intensive employment of eddy covariance has revealed limitations in this method. Here, we have started from the same theoretical foundation and derived an alternate representation of scalar vertical turbulent transport in the surface layer to address some of the weaknesses found in the eddy-covariance method, and refer to this new approach as functional covariance.
Our derivation of the functional-covariance method differs from that for eddy covariance on several points by employing, we believe, a more correct interpretation of the Reynolds Fig. 16 The sum of all four functional-covariance-flux terms presented in Figs. 14 and 15 for the fluxes of sensible heat, latent heat, CO 2 and momentum. Each panel contains curves for functional covariances determined using 30-min, 24-h and 48-h filter windows for the determination of high-frequency fluctuations decomposition of typical time-series signals. In our derivation, we decouple the averaging operator from the filter used for separating high-and low-frequency fluctuating signals, allowing us to apply averaging operators other than the temporal averaging employed by eddy covariance. Freed to choose other averaging methods, we suggest averaging of turbulent fluctuations by variables representing the primary drivers of the flux being calculated, which results in an n-dimensional array of average fluxes and their uncertainties, with each dimension associated with a driving variable. This array of functional-response fluxes may be used directly as an empirical model, or may be resampled to create a time series of flux estimates. When viewed in the context of functional covariance, eddy covariance represents a single configuration of the functional-covariance approach, in which time is the driving variable and the filter determining the high-frequency fluctuations and the averaging operator are explicitly-linked temporal averages.
The functional-covariance method has requirements beyond those of standard eddycovariance procedures. Perhaps the most important requirement is the need for an understanding of the processes and driving variables controlling the flux being calculated. The user must measure the most relevant of these driving variables at a rate that captures their temporal variability, and ensure proper temporal alignment (i.e. lag removal) of all signals. Researchers must also make decisions on the granularity of flux averaging along the driving-variable dimensions of the averaging array, in a way that balances the identification of environmental-response characteristics with data availability, in order to identify fluxresponse characteristics with sufficient certainty. Given these additional data requirements, functional-covariance fluxes have several advantages over standard eddy-covariance fluxes, with the five key areas of improvement being representation, uncertainty, data retention, low-frequency fluxes, and storage.
Functional-covariance fluxes are found to better represent the response of scalar fluxes to environmental conditions, which is an improvement arising from the closer correspondence of the turbulent flux to its driving variables because the temporal smoothing of the eddycovariance averaging operator is no longer employed. For our example analysis using CO 2 flux data, this improvement is shown to better separate the light and substrate-limited regions of the light-response curve, resulting in increases in quantum-use efficiency of ≈ 5% over important regions of the light-response curve. We also show that the more complete coverage of parameter space afforded by the functional-covariance method may assist in interpretation of the environmental responses of fluxes to drivers.
Estimates of the functional-covariance-flux uncertainty are shown to be approximately 25% smaller than those obtained from eddy-covariance methods, which is the result of the greater independence of samples contributing to average fluxes. Again, this improvement is founded on the move from averaging of temporally-contiguous data to averaging based on the driving variables.
For eddy covariance, stationarity testing can remove ≈ 30% of raw data, while the functional-covariance method groups data by similar environmental conditions, resulting in implicit stationarity conditions. This retention of data, previously lost to nonstationarity quality-control tests, results in a reduction in flux uncertainty, particularly for the niche environmental conditions that these data may represent.
As a result of the reduced uncertainties and the ability to add more data to an average flux estimate without the risk of nonstationarity, low-frequency fluxes obtained using the functional-covariance method are found to be more precise than those obtained using eddycovariance approaches. The low-frequency fluxes for our example dataset exhibit substantial low-frequency flux contributions for variations extending out to 24 h. The low-frequency contributions amount to a ≈ 25% increase in fluxes of sensible heat and CO 2 , which is an increase considerably less than the 90% increase in the flux H observed by Foken (2008). Examination of the low-frequency components of the sensible heat flux raises the possibility that instrument effects may be influencing these low-frequency flux contributions.
The residuals between functional-and eddy-covariance fluxes suggest that the flux-storage component is largely accounted for by the functional-covariance turbulent-transport term, which suggests that measurements of the the flux-storage component may not be required when using the functional-covariance method. It is hypothesized that the use of state variables when formulating the functional-covariance averaging operator results in both positive and negative rates-of-change being incorporated into the same flux estimate. If the assumption that the rates-of-change of the scalar of interest are not preferentially linked to any particular state of a driving variable, the flux-storage term should tend to zero with sufficient sampling, with the corresponding flux being attributed to the turbulent-transport term. This potential significant improvement in experimental measurement methodology requires further research and verification.
The preceding benefits depend on a competent understanding of the causes of the scalar flux and measurement of the relevant driving variables at appropriate rates, as well as the temporal alignment of the driving variables with the turbulent fluctuations. However, we suggest that a comparison of fluxes from functional-and eddy-covariance methods may provide a tool for identifying periods or conditions under which other, previously unidentified, drivers affect fluxes. These may be in the form of variables that affect the emission of the relevant scalar from the surface or may be in the form of the previously-ignored component of the scalar-transport equation (e.g. advection). With further development, the functional-covariance method will provide a useful tool for improving scalar-flux modelling.
In summary, we feel that the functional-covariance approach to calculating scalar fluxes is both a generalization and an advancement on the already-popular eddy-covariance method of flux estimation. We suggest that the additional efforts required to implement functional covariance are outweighed by the benefits it provides in the way of improved comprehension and certainty of flux estimates.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Turbulent-Flux and Uncertainty Computations
Our evaluation of turbulent fluxes using the functional-covariance method (Eq. 15) required several passes through the dataset, with the number of passes depending on pre-existing information, the choice of high-frequency-fluctuation filtering, and the desired form of the final fluxes (response functions versus time series).
As with eddy covariance, a first pass of the data determines the coordinate-rotation coefficients. A second pass, with coordinate rotation applied, determines both experimental mean a, and 30-min mean a values of scalars and velocity components. The 30-min-mean values are subsequently used as a filter to separate lowa and high frequency a fluctuations as described in Sect. 3.3.
Determination of the temporal lag of scalars and vectors relative to radiation flux is performed as part of the second pass, but could be performed in multiple passes depending on the configuration of the analysis software. Alternately, an estimated radiation-flux lag value can be employed similar to our 'back-of-the-envelope' calculation in Sect. 3.7, thus requiring no additional passes through the dataset. Current understanding of the lags associated with driving variables is poor, such that efficient means of determination requires considerable further research.
In the next pass, we set experimental mean values, apply coordinate rotation, remove the lag of scalars and vectors relative to the radiation flux, and separate the time-series into highand low-frequency fluctuations. We then create binning indices based on each of the variables used in defining the averaging operator. Using temperature as an example, we can define the corresponding averaging-bin index, All 10-Hz samples of the components of Eq. 15 (e.g. w c ) are then assigned to averaging bins using the driving-variable indices. The average flux components are each stored in a The turbulent fluxes and their uncertainties defined by Eqs. 25 through 27 characterize the flux responses to the environmental variables used in their definition. However, these fluxes still require corrections for high-frequency attenuation (Moore 1986;Horst 2000;Massman 2000) and density effects if the scalar is measured as a molar density instead of mixing ratio (Webb et al. 1980;Leuning 2004).
The frequency-response correction varies slightly from standard approaches. First, because the flux is not descended from a time period, it is not possible to employ frequencyresponse corrections based on temporally-corresponding empirical cospectra (Foken et al. 2012b), thus requiring a spectral model based on the environmental conditions. The second issue is that the low-frequency flux loss caused by a limited averaging period is now defined by the low-frequency cut-off of the filter used to define the high-frequency fluctuations (a in Eq. 18). Density correction can be applied as designed if corresponding fluxes of sensible and latent heat are calculated using the same functional-covariance approach. Here we do not apply flux corrections, and only make comparisons between methods using uncorrected turbulent-exchange estimates, while assuming that any difference in corrections of the two fluxes is small relative to other effects.

Appendix 2: Comparisons of Functional Covariance with Eddy Covariance
There are two options if we wish to make comparisons between functional-covariance and eddy-covariance fluxes. The two options arise because eddy-covariance fluxes exist in the time domain, while functional covariances are created in a multi-dimensional space defined by the averaging operator.
The first option is to average eddy-covariance fluxes by some or all of the parameters used in defining the functional covariance. In this option, we average all available eddy-covariance run-mean fluxes based on run-mean driving variables using the same procedure described in Sect. 5. For this version of the averaging operator, Q represents the total number of eddycovariance-flux values, and the j indices are the same as those described by Eq. 21 but based on run-mean variables instead of instantaneous values, resulting in P values being averaged in each cell of the matrix M, This approach should result in dimensionally-identical and comparable response functions for both methods, but may be problematic because the of limited availability of mean driving parameters corresponding to niche driving-variable combinations realized by the functional-covariance approach. It may be necessary to average over only some of the dimensions of the functional-covariance responses to make comparisons between the two approaches. The second option for comparing functional-and eddy-covariance fluxes is to use the instantaneous values of the driving variables to sample the multi-dimensional array of the functional-covariance fluxes, and average those values over the 30-min period corresponding to an eddy-covariance flux, which creates a run-period functional covariance, This approach has the benefit that the functional-covariance uncertainties may be sampled in the same manner and propagated over the eddy-covariance run period to produce a runspecific uncertainty on the run-averaged functional-covariance flux. It must also be noted that reconstructing a time series of fluxes from the functionalcovariance multi-dimensional flux array using Eq. 29 does not limit the user to producing average fluxes for consecutive 30-min periods. Indeed, estimates of time-interval-averaged fluxes can be determined for periods of any length, limited only by the corresponding availability of driving variables.
For comparison, we also calculated an eddy-covariance run-based uncertainty similar to the method of Hollinger and Richardson (2005) using fluxes from similar environmental conditions (i.e. the same conditions used for the functional-covariance binning). However, we did not limit these eddy-covariance flux errors to daily separations as it was believed that the effect of other temporally-changing controlling variables would have been small relative to those employed. Because the uncertainty obtained from the functional-covariance method provides a standard-deviation for instantaneous values, we adjusted the eddy-covariance rungrouped uncertainties by a factor of √ 18000/K to provide an equivalent sample uncertainty (where K is the number of samples in the functional-covariance flux estimate).

Appendix 3: Justification for the Use of Static Background Air Instead of Dry Air in the Scalar-Transport Equation
To begin with, we assert that the scalar-transport equation must be nonspecific in that it should be applicable to the transport of any scalar component in an atmosphere of any composition. This is necessary to avoid the use of different equations, or derivations therefrom, to define fluxes.
With regard to the density term in the scalar-transport equation, we assert that the common use of dry-air density is employed to achieve an atmosphere of constant composition within which the variations of the trace gas of interest (e.g. CO 2 ) can be defined. Specification of the dry-air density thus avoids flux errors introduced by coincident water-vapour fluxes, as outlined in Webb et al. (1980). We contend that the use of dry-air density is incorrect in that it violates our first assertion of generality of the scalar-transport equation. A simple thought experiment is sufficient to demonstrate this argument.
First, we consider a hypothetical simplified atmosphere (nitrogen, argon, water vapour, and carbon dioxide) in which we measure the fluxes of CO 2 . In the presence of large fluxes of water vapour, use of dry-air density in the transport equation would indeed provide a static background atmospheric composition in which the CO 2 fluxes could be evaluated. However, if we additionally stipulate that there are also large fluxes of argon, specifying a dry atmosphere would not render a density term of constant composition because of the large argon fluxes. In this simple situation, it is clearly appropriate to specify an argon-and water-vapour-free atmosphere when defining our scalar-transport equation.
Using the preceding argumentation, it is straightforward to assert that the flux of any atmospheric constituent, other than the scalar of interest, could contribute to errors in the determination of the scalar of interest. Considering all atmospheric components experience some level of exchange, we are left with the conundrum of specifying an "atmosphere free" atmosphere for which we can specify a density to use in the scalar-transport equation.
To address this conundrum, we hypothesize that atmospheric constituents can be divided into an active fraction that fluctuates as a result of exchange with the surface ρ a , and a static, background fraction ρ b , which is a constant component of the atmosphere. The constant component of each constituent can be identified as the minimum observable constituent concentration at the experiment location. From this hypothesis, our background static density ρ b can be defined as where ρ i represents each of the n atmospheric constituents, and the time period T is that observable by the experiment being conducted.
Considering the simple situation in which water vapour is the only contaminating flux in the atmosphere being measured, we note that ρ b > ρ dr y because ρ b = ρ − ρ v a , whereas ρ dr y = ρ − (ρ v a + ρ v b ), where ρ v a and ρ v b are the active and background-static fractions of the water vapour, respectively. For most natural ecosystems, the effect of this change is likely very small (< 2%) and the dry-air density can be used without significant loss of precision. We have included this modification for theoretical reasons. However, the corollary of this argument is that the 'active' fraction of water vapour ρ v a may be required in the density corrections described by Webb et al. (1980); this possibility requires further investigation.