Interactions between nocturnal turbulent flux, storage and advection at an "ideal" eucalypt woodland site

. While the eddy covariance technique has become an important technique for estimating long-term ecosystem carbon balance, under certain conditions the measured turbulent ﬂux of CO 2 at a given height above an ecosystem does not represent the true surface ﬂux. Proﬁle systems have been deployed to measure periodic storage of CO 2 below the measurement height, but have not been widely adopted. This is most likely due to the additional expense and complexity and possibly also the perception, given that net storage over intervals exceeding 24 h is generally negligible, that these measurements are not particularly important. In this study, we used a 3-year record of net ecosystem exchange of CO 2 and simultaneous measurements of CO 2 storage to ascertain the relative contributions of turbulent CO 2 ﬂux, storage, and advection (calculated as a residual quantity) to the nocturnal CO 2 balance and to quantify the effect of neglecting storage. The conditions at the site are in relative terms highly favourable for eddy covariance measurements, yet we found a substantial contribution ( ∼ 40 %) of advection to nocturnal turbulent ﬂux underestimation. The most likely mechanism for advection is cooling-induced drainage ﬂows, the effects of which were observed in the storage measurements. The remaining ∼ 60 % of ﬂux underestimation was due to storage of CO 2 . We also showed that substantial underestimation of carbon uptake (approximately 80 gC m − 2 a − 1 , or 25 % of annual carbon uptake) arose when standard methods ( u ∗ ﬁltering) of nocturnal ﬂux correction were implemented in the absence of storage estimates. These biases were reduced to approximately 40–45 gC m − 2 a − 1 when the ﬁlter was applied over the entire diel period, but they were nonetheless large relative to quantiﬁable uncertainties in the data. Neglect of storage also distorted the relationships between the CO 2 exchange processes (respiration and photosynthesis) and their key controls (light and temperature respectively). We conclude that the addition of storage measurements to eddy covariance sites with all but the lowest measurement heights should be a high priority for the ﬂux measurement community.


Introduction
Over the past 2 decades, eddy covariance measurements have been widely adopted as a tool for aggregate flux measurement (Baldocchi, 2003), and there are now over 650 operational monitoring sites registered with the international flux network (Fluxnet, http://fluxnet.ornl.gov).Within the Australian regional network (OzFlux, http://www.ozflux.org.au),there are 29 active sites (Beringer et al., 2016).The use of the eddy covariance technique allows continuous automated monitoring of mass and energy fluxes, and long-term multi-Published by Copernicus Publications on behalf of the European Geosciences Union.
site datasets have yielded valuable ecological insights in recent years (Baldocchi, 2008).
It has long been documented that eddy covariance measurements are prone to underestimation of the true surface efflux of CO 2 at night (e.g.Goulden et al., 1996).The key processes associated with this underestimation are storage and advection (Aubinet et al., 2012).In the former process, CO 2 may be stored below the measurement height under calm conditions.However, since the CO 2 mole fraction must be approximately preserved over longer timescales (e.g.24 h +), nocturnally respired and accumulated CO 2 is generally released with the initiation of buoyancy-generated turbulence following sunrise, such that net storage is zero.
Advection involves mean transfer of CO 2 due to the development of horizontal and vertical gradients in scalar fields; the primary process by which this occurs is the initiation of gravity-driven drainage currents on sloping terrain (Aubinet, 2008).In contrast to storage, this mechanism generally results in a net loss of CO 2 from the observing system.Since substantial respiration occurs nocturnally, this loss causes a selective systematic error towards CO 2 uptake for cumulative carbon budgets at timescales > 24 h.These drainage flows have been observed to occur on slopes of < 1 • (Aubinet et al., 2003;Staebler and Fitzjarrald, 2004).Most non-agricultural ecosystem measurement sites are on sloped terrain because historically, flat, arable land has been cleared for agriculture.Drainage-induced advection is therefore thought to occur to varying extents at most sites.
The storage term can be calculated using relatively simple instrumentation that measures CO 2 concentrations along a vertical profile between the eddy covariance instrumentation and the ground (Yang et al., 2007).In contrast, attempts to measure advection have involved deployment of complex instrumentation and delivered highly uncertain results (Aubinet et al., 2010;Leuning et al., 2008).Thus, indirect approaches to data correction have been devised, the most common of which is the identification of a threshold below which the nocturnal turbulent CO 2 flux declines with turbulent activity (as expressed by friction velocity (u * )) (Goulden et al., 1996).Since there should be no relationship between u * and ecosystem respiration, this decline is interpreted as an increase in the storage and advection terms at the expense of the turbulent flux.Data below this threshold (herein u * th ) are discarded and replaced using functional relationships between known physical respiratory drivers (primarily temperature) and nocturnal CO 2 fluxes.This u * filtering has been criticised on theoretical and practical grounds (Aubinet et al., 2012;Aubinet and Feigenwinter, 2010;Van Gorsel et al., 2007), but remains the most widely adopted approach to nocturnal data correction.It should be applied to the sum of the measured turbulent flux and storage terms (Papale, 2006), but it is quite commonly applied to measurements of turbulent flux in isolation, rather than the sum of the turbulent flux and storage terms.This is because only 10-30 % of sites globally have deployed profile systems to measure storage (D.Papale, personal communication, 23 November 2015).In Australia, only four of the 29 active sites have profile measurements, whereas ≥ 15 sites have canopies of sufficient height to warrant them (using measurement height > 3 m as a threshold for requirement).
As noted, while the net storage over time is approximately zero, it is generally positive nocturnally since CO 2 accumulates below the measurement height.In the morning, the sign reverses due to two processes: (1) the turbulent transfer of the accumulated CO 2 upwards through the measurement plane and (2) photosynthetic CO 2 uptake by the canopy.Thus, neglect of storage means that nocturnal respiratory CO 2 release is underestimated, but this is balanced by underestimation of morning ecosystem photosynthetic CO 2 uptake.
However, when nocturnal u * filtering is applied, it implicitly accounts for both storage and advection.Since there is no corresponding correction for the reversal of the storage term after sunrise, the requirement for the storage term to be approximately zero is violated, and the nocturnally respired CO 2 is effectively counted twice (Aubinet et al., 2002).This unavoidably biases measurements towards net CO 2 efflux and also affects the apparent relationship between ecosystem CO 2 fluxes and climatic controls.
Given the number of sites that do not have profile systems, it is thus important to quantify the effects of failing to measure storage.In this study, we use a 3-year record of CO 2 exchange (including storage measurements) for an Australian eucalypt woodland to investigate the interaction between nocturnal turbulent flux, storage and advection.We devise a simple method to infer the magnitude of advection and in turn quantify the apportionment of the nocturnal ecosystem CO 2 source between turbulent flux, storage and advection.We quantify the biases in annual carbon exchange that arise from neglecting the storage term and discuss their effects on interpretation of CO 2 fluxes in the context of climatic drivers.Given the significant additional investment and complexity associated with the construction and deployment of profile systems alongside eddy covariance systems, it might be argued that the incurred bias of neglecting storage could be ignored if it is small relative to other measurement uncertainties.We therefore also propagate the errors associated with determination of u * th , random measurement error and imputation (gap-filling) error to annual estimates of net carbon exchange and assess their magnitude relative to biases due to neglect of storage.

Site description
The site became operational in December 2011 (see Table 1 for site characteristics).It was established as part of a project investigating the concurrent effects of catchment reafforestation on biodiversity, carbon sequestration and stream wa- 15.9 Mean annual precipitation (mm): long term  g 560 a By biomass (although also by number of individuals in the case of the overstorey).b Determined from tree surveys in combination with allometric relationships developed for the relevant species in the immediate vicinity (see Paul et al., 2013, for further details).c Determined from allometric equations relating aboveground to belowground biomass (see Paul et al., 2014, for further details).d Estimated from direct field survey and conversion of biomass to carbon using carbon : biomass ratio of 0.45 (Chapin III et al., 2002).e Determined from direct survey of soils (to 40 cm depth) and subsequent laboratory analysis (see Cunningham et al., 2015, for further details).f Determined from site observational data for the years 2012-2014.g From nearest long-term rainfall measurement site (Mangalore Airport, Bureau of Meteorology station ID 088109).
ter yields.From an eddy covariance perspective, the site is considered "ideal" in that it is relatively flat and homogeneous within the footprint area, and the canopy is very open (leaf area index ≈ 1).While the tower is situated on a slope of approximate south-easterly aspect (Fig. 1), the slope is generally less than 1 • .Vegetation is relatively homogeneous, consisting of a sparse eucalypt overstorey (dominant species is E. microcarpa, stand leaf area index (LAI) ≈ 1; see Table 1) and a sparse shrub understorey.The open canopy reduces the potential for strong decoupling of the subcanopy space from the overlying air, and underestimation of respiration by above-canopy eddy covariance systems relative to direct chamber measurements correspondingly decreases with declining leaf area (Speckman et al., 2015).Mean annual temperature for 2012-2014 was 15.9 • C (minimum and maximum temperatures of −3.1 and 45.0 • C).Average annual rainfall for the nearest long-term rainfall measurement site (Mangalore Airport, Bureau of Meteorology station ID 088109) is 560 mm (1971-2000 average).

Instrumentation
The eddy covariance (herein EC) method was used to measure CO 2 fluxes at 36 m (mean (±SD) canopy height was 15.3 (±6.4) m, but emergent individual trees were up to 30 m).The EC method requires fast-response instrumentation to measure simultaneous variations in scalar (here CO 2 ) and vector (three-dimensional wind velocities) quantities.A Campbell Scientific (Logan, USA) CSAT3 sonic anemometer (Table 2) was used to measure wind velocities and a LI-COR (Lincoln, USA) LI7500 infrared gas analyser (herein IRGA) to measure CO 2 and H 2 O vapour mole fractions.EC data were logged at 10 Hz (post-processing is described below), and 30 min averages for radiant and subsurface energy fluxes and standard meteorology (temperature, humidity, wind speed and direction, rainfall, barometric pressure) were also logged.All data were transferred telemetrically to a central server at Monash University, Clayton, Australia.
A custom-built profile system using a LI-COR LI840 IRGA measured changes in CO 2 storage below the EC measurement height.The system consisted of an array of six gas intakes (configured logarithmically in the vertical -0.5, 2, 4, 8, 16, 36 m) connected to a sampling system (Fig. 2) via tubing (of equal length, 18 m, for all heights; enclosure www.biogeosciences.net/14/3027/2017/Biogeosciences, 14, 3027-3050, 2017  Raw data retention was generally high over 2012-2014, with the station and all critical instruments running continuously except for 23 August-25 October 13, during which the profile system was damaged and had to be partially rebuilt.Data for this period are excluded from the analysis.

Data processing and analysis
Post processing of EC data (including quality assurance and quality control) was undertaken using OzFluxQC, a software package developed by the OzFlux community (primarily P. Isaac, the data manager of OzFlux) in the Python programming language (Isaac et al., 2016).The Python programming language was used for all subsequent data analysis.Using block averaging, 30 min fluxes were calculated from the 10 Hz data (Moncrieff et al., 2004).Corrections applied to the raw data included two-dimensional coordinate rotation (Lee et al., 2004), in which the coordinate frame is rotated to force first the mean cross-wind, then the mean vertical wind, to zero over the measurement period; frequency attenuation corrections, which account for loss of covariance associated with high-frequency cutoff, sensor separation, and path averaging of then instrumentation (Massman and Clement, 2004); and density corrections associated with the effects of surface sensible and latent heat fluxes (Webb et al., 1980).Data were processed to level 4. In OzFluxQC, this represents the point at which all quality assurance-quality control (QA-QC), except u * filtering, has been applied to flux data, and meteorological data have undergone QA-QC and gap filling (Isaac et al., 2016).
We identified the u * threshold (u * th ) using change-point detection (CPD) following Barr et al. (2013).The time series was divided into multiple temperature-stratified samples, and a two-phase linear regression model was fitted to all possible change points (change in CO 2 exchange as a function of u * ) for each sample.If the change point that minimises the sum of squares error shows statistically significant improvement over a null model (no change point), the change point (i.e.u * th ) is retained (subject to additional quality control criteria, as described by Barr et al., 2013).A bootstrapping procedure (in which the data for each year were randomly sampled with replacement 1000 times; see Papale, 2006) was used to generate a probability distribution for u * th , of which the mean and 95 % confidence interval provide a best estimate and uncertainty interval for u * th .Since we are assessing the effects of neglecting profile measurements in this study, we identified u * thresholds for the turbulent CO 2 flux alone and for the sum of turbulent CO 2 flux and CO 2 storage.
The rate of change of CO 2 storage was calculated from the difference between quasi-instantaneous (2 min) vertical concentration profiles at the tower at the beginning and end of the flux averaging period (Finnigan, 2006).We adopted the approach for storage calculation of Yang et al. (2007): Here, C/ t is the time rate of change of CO 2 molar density (µmol CO 2 m −3 ), k is the profile level, z is height above the surface (m) and n is the number of profile levels.Mole fraction reported by the IRGA was converted to CO 2 molar density using the ideal gas law (temperature measurements were drawn from instruments co-located with the air inlet for each profile level, whereas pressure measurements were drawn from ground level).The 2 min period preceding the beginning and the end of the 30 min period (e.g.1128-1130 and 1158-1200 for the 1130-1200 period) were used to calculate C. Given that the pump draw was simultaneously divided across six lines, there was a lag > 1 min, so that the sampling was approximately temporally centred at each half hour.The average of the time derivative for two levels (k and k − 1) is the best estimate of the time derivative for the layer that has k and k − 1 as its upper and lower boundaries, except for the lower layer, for which it is assumed that C/ t for k = 1 is representative for the layer.Layers were scaled according to the layer thickness and the storage term represented the sum over all layers.

Analytical framework
We assess the CO 2 -carbon balance in the familiar context of a notional control volume with an orthogonal coordinate system (Finnigan et al., 2003), the mass balance of which (neglecting horizontal turbulent flux divergence) is Here, NEE (net ecosystem exchange of CO 2 ) is the true source term, term 1 is the turbulent flux across the upper horizontal plane of the control volume at instrument height h m (w is the vertical velocity, and overbar and prime denote mean and quasi-instantaneous fluctuation from mean, respectively), term 2 is the storage term integrated over a finite time period (t) and control volume depth (z), term 3 is the sum of the advection components in the horizontal dimensions (x and y, with corresponding vectors u and v), and term 4 is the vertical advection.We adopted the standard micrometeorological convention in which NEE is positive (negative) when the net transfer of carbon is from ecosystem to atmosphere (atmosphere to ecosystem).For brevity, throughout we use the term "carbon balance" to refer to long-term (annual) ecosystem-atmosphere exchange of CO 2 , and we do not quantify the much smaller exchanges of non-CO 2 gas species (CH 4 and volatile organic compounds) and potential net lateral aquatic transfer of dissolved organic and inorganic carbon.In the following text, Eq. ( 2) is simplified to Here, F c and S c are the turbulent CO 2 flux and storage terms, and the vertical and horizontal advection terms have been collapsed to a single advection term (A c ).During the day, when turbulence is well-developed, the turbulent flux (F c ) is generally the dominant term, but at night the other terms may become dominant under weak mixing.Following the identification of u * th , data were rejected where u * < u * th .There is some debate as to whether data filtering should be confined to the nocturnal period or whether it should also be applied across the full diel cycle (Papale, 2006).While the main emphasis in this study is on the nocturnal case, we assess the effects of inclusion of additional diurnal filtering on www.biogeosciences.net/14/3027/2017/Biogeosciences, 14, 3027-3050, 2017 annual carbon balances.While nocturnal advection was not measured, we took a similar approach to that of Clement et al. (2012) and inferred advection as the residual of the terms in Eq. (3).Nocturnally, NEE is equivalent to ecosystem respiration (herein ER).While ER is unknown when u * < u * th , it can be estimated ( ER) using an empirical model, the parameters of which are optimised for periods in which the sum of turbulent flux and storage approximates ER (i.e. when u * > u * th ).Equation (3) thus becomes

Model selection and imputation
Because our approach implicitly assumes that the model is a reliable estimator of the true respiratory source term, a robust model selection procedure is paramount to a reliable analysis.We used a simple set of empirical functions describing the response of respiration to temperature and soil moisture.Optimisation of model parameters was undertaken using the robust non-linear least squares implementation in the Python SciPy package.The optimised functions were used to estimate ER for u * < u * th (and for subsequent gap filling).Akaike's information criterion (AIC) was used to determine the most appropriate model.AIC is a likelihood criterion that penalises the addition of model parameters, thereby discouraging the use of models of spurious complexity, and is formulated as where L is the likelihood of the fitted model.We tested two approaches to model fitting.In the first, we simply fitted parameters to the entire dataset, and in the second we fitted some parameters at annual time steps and additional parameters at arbitrary time steps (see below).The second approach arguably allows respiratory variation associated with neglected independent variables to be implicitly accounted for.We used AIC to objectively rank the performance of models using both approaches, as well as to determine the ideal time step in the second approach.We tested linear-, simple exponential-and Arrhenius-type temperature response functions, in combination with sigmoid soil moisture response functions as described in Table 3. Linear functions have little theoretical justification as models for respiratory temperature response, and it has also been argued by Lloyd and Taylor (1994) that Arrhenius-type models have more theoretical basis than simple exponential models.However, it is possible that underlying ecosystem respiratory dynamics may be sufficiently obscured by the random error inherent in eddy covariance data that more theoretically justifiable models may nonetheless not be justifiable in practice.Moreover, since respiration is the sum of aboveground (autotrophic) and belowground (autotrophic + heterotrophic) components subject to differing climatic conditions, we also tested two-compartment models, similar to the approach of Clement et al. (2012) and Swanson and Flanagan (2001).
The model with the best performance overall was a combined temperature (Lloyd and Taylor, 1994;herein L&T) and soil moisture response model of the form For the temperature response function (term 1), T is the measured air temperature (see discussion below), R 10 is the respiration rate at the reference temperature (T ref ), E o is an activation energy parameter that controls temperature sensitivity, and T 0 is the temperature at which metabolic activity approaches zero.T 0 and T ref are fixed at −46.02 and 10 • C respectively (Lloyd and Taylor 1994) since the unconstrained version of the function is over-parameterised (Reichstein et al., 2005;Richardson and Hollinger, 2005).The soil moisture response function (term 2, adapted from Richardson et al., 2007) is a sigmoid scalar response function that is constrained to the interval [0, 1] and effectively modifies R 10 for the effects of low soil moisture.The function only accounts for the effects of low soil moisture (θ 1 and θ 2 parameters control the x intercept and gradient respectively).We tested an alternative soil moisture response function (see Table 3) that accounted for effects of both low and high soil moisture, but performance was poor.This is most likely due to several factors: (1) the site is primarily water-limited, such that there is limited data available to meaningfully model effects of high soil moisture; (2) the wettest period coincides with low soil temperature (winter) such that nocturnal signal : noise ratio for NEE is expected to be small, making it difficult to identify saturation effects; and (3) the alternate two-parameter soil moisture response function sacrifices flexibility at low soil moisture, which is likely to be the primary moisture-related constraint on respiration in this ecosystem.
The two-compartment L&T model also performed poorly, with and without soil moisture.This may reflect some combination of the low signal : noise ratio at this relatively low productivity site (see below) and/or the low soil carbon storage at this site (Table 1).Low soil carbon does not necessarily imply that the belowground respiratory contribution must be small since it may be driven by a combination of autotrophic (root) and heterotrophic (supported by exudates from the roots) respiration.However, we also tested numerous weighted combinations of soil and air temperature as inputs for Eq. ( 6) and found that (i) a weighting of 0 for soil temperature and 1 for air temperature produced the lowest RMSE (which will also translate to the lowest AIC since there is no change to the number of parameters) for all air temperature heights tested (0.5, 2, 4, 8, 16, 36 m) and (ii) the worst-performing temperature was soil temperature, and the best-performing was 36 m air temperature.Lloyd and Taylor (1994 This is well above the average canopy height (15.3 m, see Table 1).This is unexpected since the temperature that yields the best model prediction should be the temperature associated with the ecosystem stratum where the maximum respiratory production is occurring.However, the 36 m temperature measurement shows very similar dynamics to the surface radiant temperature, both as a function of time of night and of u * (Fig. 3).The surface radiant temperature in this case aggregates ground, trunk space and canopy level temperatures (since the canopy is very open, LAI ≈ 1.1), and when used as the temperature driver for model optimisation, it produced an almost identical RMSE to 36 m temperature (0.953 vs. 0.951 µmolCO 2 m −2 s −1 respectively), lower than for all other temperature measures.We hypothesise that a substantial proportion of CO 2 may be sourced from the litter layer, which contains over 4 times more carbon than the upper 40 cm of the soil (Table 1).Litter contributions to soil respiration exceeding 50 % have been reported in the literature (e.g.Xiao et al., 2014).
Other authors (e.g.Lasslop et al., 2010;Reichstein et al., 2005) have used the temperature response function in isolation, fitting E o at an annual timescale and allowing R 10 to vary on a smaller time step.Thus, variations in R 10 may implicitly capture the effects of environmental controls (such as soil moisture) not explicitly included in the model.This was considered inappropriate for this ecosystem since there may be rapid responses of ER to episodic rainfall following dry conditions at timescales shorter than the chosen step.Moreover, because soil moisture constraints coincide with high summer temperatures, this can force the E o parameter to extremely small values (i.e.such that the respiratory response is close to zero) when the complete seasonal cycle is included in the parameterisation.We nonetheless found that the best performing version of the model above (i.e.lowest AIC) used an annual time step for fitting of the E o , θ 1 , and θ 2 parameters; simultaneous fitting of soil moisture response ensures that the E o parameter is not artificially reduced by effects of soil moisture.The remaining R 10 parameter was then fitted on a 7-day time step (with R 10 for intervening days linearly interpolated), which was chosen because it yielded the lowest AIC (Fig. 4a; the only other candidate model with a likelihood even close to -although nonetheless inferior to -this model was the exponential model combined with the same sigmoid soil moisture response function, fitted with the same time step).This suggests that there is real subannual variation remaining in the ER signal at frequencies of ≤ 7 days, which may include synoptically or seasonally driven effects.While additional real respiratory variance may also be present at frequencies > 7 days, the results suggest that they cannot be separated from the noise.
While the coefficient of determination (r 2 ) for the chosen model was low (0.18), this appears to be primarily a function of the relatively low respiratory CO 2 production associated with this low-productivity site.Since there is irreducible random error in eddy covariance and profile measurements (Hollinger and Richardson, 2005), even a "perfect" model cannot account for 100 % of variance in observations.Since random error as a function of flux magnitude has a non-zero intercept and modest gradient (see Sect. 3.5 and Fig. 13), the relatively low nocturnal respiration at the site results in low signal : noise ratio, nocturnal EC measurements at the site.This is compounded when NEE represents the sum of the turbulent and storage fluxes because the latter substantially additionally increases random error (see Finnigan, 2006, for further discussion).Thus, higher r 2 values observed in the literature (e.g.> 0.6; Clement et al., 2012) most likely pertain to more productive ecosystems with a strong seasonal cycle.To test the effect of random error on model parameterisation, we used a procedure to characterise the random error in our measurements (see Sect. 2.6 below), then superimposed this error on the pure model signal produced using the parameters derived from the optimisation in combination with the observational temperature and soil moisture data (Fig. 4b).At higher temperatures (> 22 • C), the model tends to some-what overestimate ER.We suspect that this reflects the fact that these temperatures occur during a period of prolonged dry conditions when the effects of low soil moisture could have sustained effects (e.g.episodes of drought-induced leaf senescence) that are not captured through an instantaneous soil moisture response.These instances represent < 3 % of all data.
The variance in the synthetic data explained by the model (13 %) was actually smaller than for the observational data (18 %).This apparent paradox is explained by the fact that the methodology for characterising random error is known to overestimate random error (by a factor of 1.5 to 2; Billesbach, 2011; Dragoni et al., 2007).Reducing the random error by these amounts results in r 2 of 0.25 and 0.37 respectively (see Fig. 4b inset).Thus, a perfect model degraded with a conservative estimate of random error would still only explain 37 % of the variance in NEE at maximum.Correspondingly, up to 20 % (38-18 %) of the unexplained variance in our observations may be associated with model error (due to missing drivers or mischaracterisation of the relationship between drivers and response).
We nonetheless ascertained that the R 10 parameter from the model could be reliably extracted at the chosen time step even with the inflated random error estimate.We ran 10 4 trials, in each of which ER estimated by the model was degraded with a realisation of random error and the parameters were then estimated again (Fig. 4c).The 95 % confidence interval shows that the dynamics of R 10 were reproduced even when the data were degraded (again suggesting that the variability in R 10 is real, rather than an artefact of random error).Estimated cumulative annual ER ± resulting uncertainty (95 % confidence interval for the mean) was 995 ± 32 (3.2 %), 905 ± 57 (6.3 %) and 1028 ± 34 gC m −2 a −1 (3.3 %) for 2012, 2013 and 2014, respectively (the higher error in 2013 reflects missing data during late winter-mid-spring, as previously noted).This approach implicitly assumes that the initial model is a true descriptor of reality.Therefore, we can conclude little about the potential role of the above-noted model errors, but it does allow us to conclude that random error alone is not a barrier to reasonably robust parameter estimation, and that associated parameter and corresponding ER estimation errors are relatively small.
While the emphasis in this study is on the effects of nocturnal data treatment, gap filling was also required for daytime to assess the effects of nocturnal data treatment on annual NEE.We used a Michaelis-Menten-type rectangular hyperbolic model (Ruimy et al., 1995) of modified form (Falge et al., 2001) to estimate NEE, where ER was calculated from Eq. ( 5) using daytime temperatures in conjunction with nocturnally derived parameter estimates: Here, α is the initial slope of the photosynthetic light response, Q is photosynthetic photon flux density, and β is photosynthetic capacity at 2000 µmol photons m −2 s −1 .The same step size and interpolation procedure as for the nocturnal fitting of the respiration model were used.We adopted the additional light-response model criterion in which A opt is modified to include a non-linear scaling factor to account for the effects of vapour pressure deficit (VPD) on stomatal conductance (Lasslop et al., 2010): Here, VPD 0 is a threshold value above which stomatal conductance becomes sensitive to VPD, and k is a fitted parameter defining the β response to VPD.

Uncertainty estimation
We quantified sources of uncertainty in the data arising from random measurement and model error.We calculated random error from a daily differencing procedure (Hollinger and Richardson, 2005).When differences in critical drivers are sufficiently small (< 35 W m −2 for insolation, < 3 • C for air temperature, and < 1 m s −1 for wind speed), differences between NEE data pairs separated by 24 h were considered to represent random error.Since random error in EC data is heteroscedastic and follows the Laplace distribution, we calculated the standard deviation of the error (σ [δ]) for j samples in iu * bins as www.biogeosciences.net/14/3027/2017/Biogeosciences, 14, 3027-3050, 2017 σ [δ] was regressed on the flux magnitude to derive a linear relationship that was in turn used to estimate σ [δ] for each datum.Monte-Carlo-style simulation was used to translate these estimates into annual uncertainty, whereby for each of the 10 4 simulations, estimates of random error for all observational data were aggregated over 1 year.This yielded a normal distribution of uncertainties, the 2σ bounds of which were taken as the annual uncertainty.
With respect to model error, we followed Keith et al. (2009) in which, for day and night conditions, a subsample of 10 3 records with observational data was randomly selected from the annual dataset.A proportion of the observational data in the subsample equal to the observed proportion of data missing in the annual time series was then replaced with model estimates, and the summed difference between the complete observational and gap-filled subsamples was calculated and expressed as a proportion of the observational sum.This was repeated 10 4 times, and the 2σ bounds of the proportional error were calculated and then applied to the annual sum to produce an absolute annual model error.
To combine uncertainties, we assumed independence of the random and model estimates and sum in quadrature: Here, ε tot , ε r and ε m are combined total, random and model uncertainty respectively.
Other unquantified sources of error may also be present.Barr et al. (2013) argued that one of the key sources of uncertainty in annual NEE is the estimation of u * th .We used their bootstrapping approach to derive robust confidence intervals for u * th , and we included estimates of the effect on annual NEE of setting u * th to the upper and lower bounds of the 95 % confidence interval for u * th .

Contribution of mass balance components to nocturnal carbon dynamics
A nocturnal relationship was observed between friction velocity (u * ) and F c (Fig. 5), with F c declining quasilinearly below u * th = 0.42 m s −1 (identified using change-point analysis of F c ) and approaching zero at zero turbulence.There was little seasonal variation in u * th (data not shown), and the estimates and uncertainty bounds for all years were similar (Table 4).Given that the primary abiotic respiration controls, temperature and soil moisture, showed no relationship with u * except at u * < 0.08 m s −1 (which is linked to declines in temperature), we interpreted this as an increase in the nonturbulent terms of the mass balance.
The decline in F c with declining u * was accompanied by a corresponding increase in S c because as turbulence is progressively suppressed, CO 2 is expected (excluding advective losses) to be increasingly stored below the measurement height.The strong sensitivity of CO 2 accumulation to u * th below the measurement height was observed in raw time series data under varying u * (Fig. 6); the CO 2 mole fraction responded very sensitively to variations in u * , and the effect of u * crossing u * th is evident.
However, S c was not the only important additional term in the nocturnal mass balance.A u * -dependent decline in F c + S c was also evident (Fig. 7) below a threshold of u * = 0.32 m s −1 (identified using change-point analysis of F c ), al-  though this was inherently more uncertain due to the much higher random error in storage relative to the EC measurements (Table 4).There is no plausible explanation for u *dependent declines in the respiratory CO 2 source (again, excluding changes in relevant controls), and so we infer that this represents CO 2 losses associated with the remaining mass balance terms (i.e.advection).This can be quantitatively estimated as a residual following Eq.( 4).Note that parameter optimisation of the temperature response function used F c + S c as the target variable because S c is observed to be non-zero when u * u * th (see Fig. 5, and subsequent discussion).
The terms in Eq. ( 4) are plotted as a function of u * in Fig. 7.The inferred advection estimate increased rapidly below u * threshold = 0.32 m s −1 and was comparable to S c at the lowest u * values.This indicates that under the calmest conditions F c approached zero, and approximately half of the CO 2 respired by the ecosystem was stored below the measurement height, while the remainder was advected away.Integrated over the interval 0 < u * < u * th , S c accounted for 61 % of the difference between F c and ER, with the other 39 % attributed to the advection components (A c ).This indicates that even on relatively flat terrain, such as that observed at this site, the nocturnal advection term is significant.
Since advection processes are not necessarily confined to the night, we also conducted the same analysis for the daytime.Daytime data were filtered using the nocturnally derived u * threshold derived for the sum of the turbulent flux and storage terms (0.32 m s −1 ).Over most of the range of u * values, on average there was little departure between the observations and model (Fig. 8), except at very low u * (and to a lesser degree at the highest values of u * ).Assuming that advection is responsible for model-data departures as a function of u * , this indicates some possible residual early morn-  ing advection since this is when the overwhelming majority of low u * daytime conditions occur.

Inferred advection mechanisms
S c for the individual layers is presented in Fig. 9.There was a clear decline in S c for all layers below 8 m when u * was less than approximately 0.25 m s −1 .In contrast, storage in the higher 8-16 and 16-36 m layers continued to increase near linearly.We hypothesise that these results indicate the onset of drainage flows at low levels under stable conditions, causing horizontal advective losses of CO 2 preferentially from the lower layers of the control volume.Since the presence of advection is inferred, it is not possible to directly assess the relative contributions of horizontal and vertical advection, but several factors indicate that horizontal advection is most likely the dominant mechanism.Both the presence of a long and relatively consistent upward slope to the northwest of the tower (3-4 km to the ridgeline; see Fig. 1) and often clear and calm nocturnal conditions associated with the site's low altitude and continentality are conducive to the development of terrain-induced flows.This may be offset to some degree by the open, sparse canopy structure, which, while aiding surface radiant heat loss, also inevitably reduces decoupling of subcanopy meteorological conditions from the overlying atmosphere (Speckman et al., 2015).
The presence of katabatic flows alone does not imply advective CO 2 loss from the control volume -there must also be horizontal inhomogeneity in the CO 2 scalar field (see Eq. 2).This may arise when the upwind flux source area (or footprint) extends beyond the limits of homogeneous vegetation cover.If the footprint at Whroo exceeds the distance to the local ridgeline under katabatic flow conditions, relatively CO 2 -depleted air will be mixed downward and entrained into the surface flow.The mid-slope location of the site is also key to horizontal advective CO 2 loss; by contrast, for towers in valley bottoms, CO 2 would instead be expected to accumulate, likely resulting in advective CO 2 gain and corresponding increases in S c .
Moreover, drainage flows are generally confined to the trunk space below the canopy (Aubinet et al., 2003).At the study site, mean canopy height was 15.3 ± 6.4 m (SD, Table 1).Conservatively assuming that the canopy comprises the upper 30 % of tree height, drainage flows may be confined to depths of 10 m, comparable to commonly reported values (Goulden et al., 2006;Mahrt et al., 2001).The ongoing increases in storage in the 8-16 and 16-36 m layers suggest that the CO 2 source in these layers originated primarily from the vegetation rather than vertical transfer from lower layers.
In the interval between the u * thresholds for F c alone and F c + S c (i.e.0.32 ≤ u * ≤ 0.42 m s −1 ), A c was not significantly different to zero (see Fig. 7).The linear relationship between each of the lower layers (0-0.5, 0.5-2, 2-4 and 4-8 m) and the mean 8-36 m layer in the interval 0.32 ≤ u * ≤ 0.42 m s −1 (i.e. the change-point-derived u * threshold for F c + S c and F c ) can be used to approximate the expected rate of change for those layers in the absence of advection when u * < 0.32 m s −1 .Extrapolation of this linear relationship to conditions where u * < 0.32 m s −1 provides an estimate of the expected magnitude of S c in the absence of advection.If drainage flows are the primary advective mechanism, then the sum of F c and the linearly adjusted storage term should approximate ER.The correction to the storage in the 0-8 m layers when u * < 0.32 m s −1 increased the 0-36 m storage term such that for the interval 0 < u * < 0.32 m s −1 the mass balance was approximately closed because ER − F c ≈ S c to within the uncertainty (95 % confidence interval) of the bin means over this interval (Fig. 10).This indicates that the decline in S c at lower layers was of approximately the same magnitude as the inferred advection, consistent with the presence of low-level drainage flows removing CO 2 from the control volume.
The presence of katabatic flow may be observed as vertical shear in wind profiles, with terrain-aligned flow at low levels (below canopy) and synoptically aligned flow at higher levels (above canopy).Unfortunately, the installed wind instrumentation (RM Young Wind Sentry set; see Table 2) lacks the requisite resolution (minimum detectable wind speed = 0.5 m s −1 ; minimum speed required to effect 10 • deflection in wind direction = 0.8 m s −1 ) to detect weak drainage flows characteristic of moderate vegetated terrain (typically less than 0.5 m s −1 ; Aubinet et al., 2003;Goulden et al., 2006;Mahrt et al., 2001).While temperatures at lower levels also declined (slightly) more rapidly than at higher levels (Fig. 3), this cannot plausibly explain low-level storage declines.While cooler lowerlevel temperatures and lower u * are causally interrelated, they also have opposing effects on S c , whereas decreasing temperatures reduce ER and decreasing u * increases the diversion of respired CO 2 into S c relative to F c (absent A c ).Thus, the steady temperature decline relative to u * cannot explain the abrupt change in sign of the relationship between storage and u * at low levels.

Effects of correction methods on diel and annual NEE
The importance of the mass balance components to the diel carbon balance is presented in Fig. 11 (note the data pre-sented in the main figure -and discussed below, except where indicated -are based on the use of a diel u * filter).
Given that the contribution of S c must average approximately zero over the diurnal cycle, annual NEE sums for both F c and F c + S c were comparable: approximately −450, −400 and −560 gC m −2 a −1 for 2012, 2013 and 2014 respectively (Table 5; small differences (< 20 gC m −2 a −1 ) were observed for F c versus F c + S c , reflecting small differences in parameters of functions used for gap filling).However, the amplitude of the diel cycle also increased with the addition of S c (as well as a small daytime phase shift in peak CO 2 uptake: from midday -synchronous with the solar radiative peakon average, to about 1100).The application of the u * correction increased the nocturnal respiration estimate substantially, indicating the presence of advection under weak turbulence, as previously discussed.As previously noted, there was also evidence of continued advection following sunrise, which is expected given that substantial surface heating is re- The inset plot in Fig. 11 shows the difference in NEE dynamics over the diel cycle arising from application of the 24 h versus nocturnal-only u * filter.For F c + S c , the difference was small and mostly confined to the early morning, when, as noted above, the additional filtering most likely accounts for a small amount of residual advection between 06:00 and 08:00.
The corresponding net decrease in annual C uptake was small (< ±10 gC m −2 a −1 ); however, it was also inconsistent between years, indicating the probable presence of more complex dynamics during this transition period.
Given that the majority of sites both internationally and in Australia do not have profile measurement systems, below we discuss the effects of neglect of S c because this is the de facto approach taken for sites without profile systems.As a secondary option, a single-point storage term (herein S c_pt ) can be derived from the EC gas analyser.This will increasingly underestimate storage for taller towers, where much CO 2 accumulates within the control volume and is subject to substantial error (Gu et al., 2012;Yang et al., 2007), but may nonetheless potentially reduce bias.
On average, S c_pt performed very poorly nocturnally (relative to S c ), strongly underestimating CO 2 accumulation within the control volume (Fig. 12).This is expected since the profile system is required because the accumulation of CO 2 principally occurs below the measurement height in the absence of well-developed turbulence.Performance was also poor during the early morning transition when the change in sign of S c_pt lagged S c by up to 2 h.However, when summed with F c and filtered for low u * conditions (u * th for F c + S c_pt as determined by change-point analysis was 0.33, 0.30 and 0.35 m s −1 for 2012, 2013 and 2014 respectively), performance was substantially improved, as long as the filtering was applied across the diel cycle.This corrected the nocturnal data and excluded a substantial proportion of the problematic early morning data such that the diel dynamics were relatively similar to u * -corrected F c +S c (Fig. 11).Compared with nocturnal-only u * filtering, diel u * filtering of F c +S c_pt slightly increased estimated CO 2 release before sunset, but substantially increased CO 2 uptake just after sunrise (inset plot in Fig. 11).In fact, the peak in the increased uptake after sunrise (−0.6 µmol CO 2 m −2 s −1 ) was similar in magnitude to the shortfall in storage from S c_pt relative to S c , thus bringing the early morning NEE dynamics for diel u * -corrected F c +S c_pt into approximate agreement with the NEE estimate derived using S c from the profile system.
This correspondingly yielded estimates of relatively low bias (± < 20 gC m −2 a −1 ) in the diel u * -corrected estimates relative to F c + S c for the years 2012 and 2014, but caused more substantial overestimation of the net carbon sink for 2013.For both 2012 and 2014, the reported annual C sink was increased because much of the data during the early morning when S c_pt tended to be too high was removed and replaced with data modelled using higher u * conditions (under which S c_pt was a smaller component of the mass balance relative to F c , thus reducing bias).The inconsistency in 2013 arose due to the combined effects of error on model parameterisation and subsequent imputation, in combination with the effects of missing data.Random error estimates (derived from Eq. 8) are shown in Fig. 13; at low flux magnitudes the error associated with F c + S c_pt was large relative to that for www.biogeosciences.net/14/3027/2017/Biogeosciences, 14, 3027-3050, 2017 F c alone or F c + S c .This may be due to a mismatch in the timing of peak S c_pt relative to S c ; small day-to-day variations in this timing may -in conjunction with noise in S c_pt -result in large variations in NEE estimates for a given set of environmental conditions (as essentially demonstrated in Fig. 13).Large errors in point-based storage estimates were also reported by Gu et al. (2012) for a forest ecosystem.This in turn affects the parameter estimates for R 10 in the respiration function (Eq.5) and α in the light-response function (Eq.6) in particular.There was a prolonged period of missing data during 2013 (late August to late October), and the estimate of NEE during that period is strongly dependent on the parameter estimates immediately preceding and succeeding the data gap.When these estimates are subject to greater error, the interpolated values in the gap are similarly subject to greater error.This effect can also be seen in the larger uncertainty bounds for R 10 estimates derived from u * -filtered F c + S c presented in Fig. 4c.Correction of F c alone (using diel u * filtering) yielded comparable nocturnal NEE to that estimated from u *corrected F c +S c because it implicitly accounted for the nocturnal effects of both storage and advection by eliminating data where those terms were large (Fig. 11).There was still a small shortfall nocturnally because, particularly in the early evening, S c was on average 0 when u * > u * th , such that the respiration model used to fill the gaps was unavoidably optimised with a biased estimate of F c .The NEE estimate during the early morning was improved relative to uncorrected F c , again because this removed data during the period when S c was relatively large and F c was a correspondingly biased estimator of NEE.The effect of filtering low u * conditions during the day can also be seen in the inset panel of Fig. 11 and had a similar effect to daytime filtering of S c + F c_pt .How-ever, as noted with the nocturnal case, in the early morning the magnitude of S c was initially large even after u * > u * th .As a result, the F c data used for optimisation were still biased.
Why was S c non-zero when u * > u * th ?In theory, in the absence of a net CO 2 source-sink or advective CO 2 loss, the presence of a vertical gradient in CO 2 density at the morning resumption of turbulence must result in turbulent efflux of CO 2 from the control volume, which in turn must -as dictated by mass conservation -be balanced by a change in storage of equal and opposite magnitude.Under such conditions this would be expected to manifest as a spike in F c (for example, see Aubinet et al., 2012, and Fig. 5.2 therein).However, the net source-sink term is rarely zero.In the morning, photosynthetic activity began at approximately 06:00 on average, as evidenced by the sudden decline in S c at this time (Fig. 12), crossing zero at approximately the time that u * reached its diel minimum.The ecosystemlevel light compensation point (when gross primary production (GPP) ≈ ER, as indicated by F c + S c = 0) occurred between 07:00 and 08:00, with net photosynthetic CO 2 uptake thereafter (see Fig. 11).By consuming CO 2 within the control volume, photosynthetic activity necessarily reduces the vertical CO 2 gradient and thus commensurately reduces the turbulent flux across the upper plane of the control volume (the situation is more complex than this: when turbulence is weak, storage changes within the canopy may not have a substantial effect on the CO 2 gradient across the upper plane of the control volume, such that F c and S c are not immediately coupled).As such, there was no identifiable spike in F c .However, whether the change in storage occurs as a result of turbulent ventilation of CO 2 from the control volume or photosynthetic uptake, this change is missed by the EC measurement system and NEE estimates based on F c alone are consequently biased high.Over time, turbulent transfer and photosynthetic consumption steadily eliminate the accumulated CO 2 surplus in the control volume, but the process is not instantaneous, thus the presence of large initial early morning S c even when u * > u * th .
Early morning photosynthesis may be particularly pronounced in eucalypt-dominated ecosystems because the characteristically pendulous (in some species up to 75 % of mature leaves typically hang at angles > 80 • from horizontal; Pereira et al., 1987), amphistomatous leaves evolved to maximise incident radiation at low sun angles.This shifts photosynthetic activity towards periods with lower vapour pressure deficit, reducing water losses (James and Bell, 1996).Mutual shading would partially counteract this effect at low sun angles, but this may have less effect in systems with sparse canopies such as the woodland in this study.An Australian temperate eucalypt forest site with long-term turbulent flux and CO 2 profile measurements also showed no morning spike in CO 2 efflux (Van Gorsel et al., 2007, and Fig. 4 therein).
Effectively, the reverse process began before sundown.The net surface CO 2 sink transitioned to a source between 17:00 and 18:00, as GPP declined relative to ER.This is because temperature (the key driver of respiration) lags insolation (the key driver of photosynthesis), while higher vapour pressure deficit may additionally induce partial stomatal closure.This was accompanied by a rise in the quantity of CO 2 partitioned into S c despite u * > u * th .This is most likely again related to the effects of the prior conditions on the vertical CO 2 gradient.The photosynthetic drawdown of CO 2 during the day resulted in a CO 2 deficit within the control vol-ume (relative to air at higher levels) that persisted until the evening.The minimum CO 2 mole fraction occurred at the afternoon zero crossing of S c , which is simply the time rate of change of CO 2 mole fraction, scaled appropriately.Thus, the change in sign of storage may result partially from the reversed concentration gradient (driving mixing of CO 2 downward from aloft), but by the late afternoon the system was an increasing net CO 2 source, such that a large proportion of the CO 2 entering the control volume from the ecosystem acted to increase the CO 2 mole fraction within the control volume.Again, this is inevitably missed by the EC system, and thus NEE estimates based on F c alone are necessarily biased low.Over time, turbulent transfer and photosynthetic consumption steadily eliminate the accumulated CO 2 deficit in the control volume, but the process is not instantaneous, thus the presence of large initial early evening S c even when u * > u * th .
Note however that this early evening bias does not affect annual NEE estimates because it is balanced against equivalent daytime biases.The diel dynamics of storage can be divided into two 12 h periods, over each of which the sum of S c ≈ 0: 21:00-09:00 and 09:00-21:00 (these times are chosen because for 21:00-09:00, CO 2 mole fraction on average was greater than its mean value, whereas for 09:00-21:00, it was less).In the first period, if only F c is used to estimate NEE, a potential double-counting problem arises in the morning because of the fact that the neglected change in storage partially reflects the removal of nocturnally accumulated respired carbon, which has already been accounted for by the nocturnal correction.In the second period, the underestimation of true NEE due to neglect of S c is balanced across the period (the sign and magnitude of S c between 09:00 and www.biogeosciences.net/14/3027/2017/Biogeosciences, 14, 3027-3050, 2017 approximately 16:00 balances the period between 16:00 and 21:00).Thus, the corresponding annual NEE estimates derived from diel u * -filtered F c were in all years biased towards net carbon efflux (by 47, 39 and 45 gC m −2 a −1 for 2012, 2013 and 2014 respectively) relative to the corresponding best estimate of NEE (i.e.diel u * corrected F c + S c ).In comparison to uncorrected estimates of annual NEE (i.e.gap-filled raw data), this reduced the absolute bias but reversed the sign.The performance was more consistent but on average more biased than for F c + S c_pt .The higher bias is expected because although using the diel u * filter resulted in comparable improvement in NEE across the diel cycle (see inset of Fig. 11) for both F c and F c + S c_pt , the improvement in F c was not added to an existing (albeit imperfect) estimate of S c as was the case with S c_pt .Thus, some proportion of the contribution of the morning negative storage term to the mass balance was still missed, resulting in underestimation of CO 2 uptake on average.
Far more problematic for estimation of annual NEE was the effect of applying a nocturnal u * filter only.This unmasked the entire morning period (when u * < u * th ) during which the magnitude of S c was relatively large and F c correspondingly diminished.The effect was such that applying the correction changed the sign of the bias (as in the above case), but also (on average) increased its magnitude.Conversely, the uncorrected NEE estimates were biased towards uptake by 66, 61 and 65 gC m −2 a −1 and nocturnal u * -corrected F c was biased towards efflux by 79, 60 and 92 gC m −2 a −1 .In other words, in the absence of storage estimates, a more accurate estimate of annual NEE would be obtained at this site if the nocturnal u * correction were not applied.This is expected given the fact that S c contributed more to the nocturnal underestimation of F c than advection.Whereas advected CO 2 is (in this case) lost to the sensing system, the stored CO 2 causes an equivalent offset in the turbulent flux measurements when it is evacuated from the control volume in the morning (via the mechanisms described above), thereby biasing the system towards efflux.The mechanism is as described above, but is more severe when u * filtering is not extended to the day.If the proportional contributions of S c and A c are approximately equivalent, the application of nocturnal u * filtering to F c alone would be expected to preserve the magnitude but reverse the sign of the induced bias.In theory, where S c is dominant (in this study we ascertained a contribution of approximately 61 %), the bias will necessarily be worsened by this treatment (since the quantity that needs to be corrected -A c -is smaller than the quantity that should not be corrected but implicitly is -S c ).By the same logic, where A c is dominant, the bias should be reduced.In practice, sites with severe advection problems may also be subject to flow decoupling that renders u * filtering ineffective (e.g.Van Gorsel et al., 2007), but the issue is not confined to the u * filtering methodology per se.Any method of correcting nocturnal data that does not take account of the subsequent daytime effects imposed by the nocturnally stored CO 2 will be subject to the same issues.

Effects of neglecting CO 2 storage on physiological interpretation of data
From the perspective of deriving annual NEE sums, it is daytime rather than nocturnal measurements that are more critical; applying u * filtering to either F c or F c + S c resulted in similar estimates of nocturnal NEE on average.However, the effects of neglecting S c depend on time of night.F c underestimated NEE following sunset, even where u * > u * th for F c .
A secondary nocturnal problem is also recognised in the literature: when u * increases following extended calm periods, stored CO 2 is vented from the control volume, which artificially inflates F c relative to the true source term, the extent of which effect will depend on the importance of advection (Aubinet et al., 2012).When a u * threshold is imposed, such periods are likely to be included in the retained data.This has the opposite effect to the early evening effect and is more likely to be problematic later in the evening when stable stratification and substantial storage of respired CO 2 is more likely.Both effects were observed in the nocturnal progression of S c for periods when u * > u * th (Fig. 14): S c was on average > 0 for the first 4-6 h after sunset and < 0 afterwards.On balance, the effect in this study was to increase slightly the estimation of ER.This explains why S c was slightly positive when u * > u * th in Fig. 5, and why F c was slightly lower than F c +S c at night (primarily in the early evening) in Fig. 11.
However, given that temperature decreases during the evening, this suggests that the slope of temperature response functions will be slightly increased for F c + S c versus F c alone.Given that the optimisation procedure minimises the prediction error, this may not have a large quantitative effect averaged over the evening, but interpretation of system response to temperature is distorted.Moreover, extrapolation beyond the parameterisation domain (e.g.estimation of daytime ER) may result in substantial error because distortion of function parameters (e.g.E o and R 10 in Eq. 5) will potentially result in systematic error (because the function optimised using F c will underestimate NEE at high temperatures).Any systematic error in estimated daytime ER will then necessarily propagate to estimation of GPP (commonly calculated as NEE -ER).Because these errors are offsetting, this is not likely to have a large effect on annual NEE estimates.
Similar distortion of response to insolation occurs during the day.The addition of S c substantially affects diurnal NEE dynamics, particularly during the morning, which affects the interpretation of the controls on NEE.For example, Fig. 15 shows the difference in radiation use efficiency (RUE, here simply defined as the ratio of mean NEE to mean insolation) during daylight hours for F c alone versus F c + S c .RUE was higher in the early morning and declined more sharply when NEE = F c + S c .Such declines in RUE are often associated with stomatal response to increasing VPD, and so the importance of this driver may be missed or minimised when S c is not measured.Filtering low u * conditions during the day substantially improved the light response of F c , again because it removed a large proportion of the data during periods when S c was large.It also moderately increased the RUE of F c +S c in the early morning, again most likely reflecting the removal of the small previously identified effect of residual advection in the early morning when u * < u * th .
Application of light-response function analysis to daytime data to extract either photosynthetic or respiratory parameters is problematic in the absence of storage measurements because F c = NEE during most of the day.The estimation of ER (and quantum efficiency) derived from light-response function analysis (e.g.Gilmanov et al., 2003;Lasslop et al., 2010) is strongly dependent on the magnitude of observed NEE when insolation is low (sunrise and sunset), and thus the effect of neglecting the storage term may be particularly distorting to these parameters (Aubinet et al., 2012).

Sources of uncertainty
While bias is to be avoided where possible, if the known effect of bias is small relative to the uncertainty in NEE, then it may not be of particular concern.Here we analyse the magnitude of several uncertainty sources relative to the magnitude of bias (given the inconsistency problems with the point-based storage estimate, here we only compare F c with F c + S c ).One of the largest sources of uncertainty in annual NEE estimates is expected to derive from uncertainty in u * th (Barr et al., 2013;Papale, 2006).We propagated this uncertainty to annual NEE (for both F c and F c + S c ) by filtering and gap-filling the data using the lower and upper bounds of the 95 % confidence interval (CI) of the normally distributed population (N = 10 3) of u * th derived from CPD (Table 4).Much larger effects were evident for the lower uncertainty bound (µ − 2σ , where µ is the best estimate for u * th , and σ is the standard deviation), which is to be expected because systematic errors in nocturnal flux measurement occur at low u * .However, there should be no systematic variation in NEE when u * > u * th .The direct effect of the upper uncertainty bound (µ + 2σ ) is expected to be minimal.While the reduction of nocturnal data availability for higher u * th is expected to increase parameter imputation uncertainty, the effect here was minor, with annual NEE for u * th = µ and u * th = µ + 2σ differing by < 10 gC m −2 a −1 in all years.
The uncertainty in u * th was greater for F c + S c than for F c alone due to the additional random error inherent in S c (see Finnigan, 2006, for further discussion), which feeds into the change-point detection process.This propagated to larger uncertainty in the lower bound for annual NEE (50-75 gC m −2 a −1 compared to 20-40 gC m −2 a −1 for F c ).This is because for F c +S c , the lower bound of the u * th uncertainty is below the first percentile of the nocturnal data, such that the full effect of advection was propagated to annual NEE, but for F c alone, it was closer to the 40th percentile, such that only a small proportion of storage and advection occurred in the interval between u * th = µ and u * th = µ − 2σ .The best estimates of NEE for F c alone versus F c + S c are sufficiently different that their respective uncertainty ranges do not overlap, indicating that biases are large relative to u * -induced uncertainty.
The large lower-bound uncertainty in annual NEE derived from F c + S c is very likely overestimated.If the lack of a relationship between u * and F c above u * th indicates that additional terms in the mass balance are negligible, then S c should on average approach zero at u * th .Given that this is what we observed (Fig. 5), and the measurement system (and associated measurement errors) for S c is independent of that for F c , this is an independent validation of u * th .It is not clear how such information might be used in the context of frequentist statistical analysis, but it strongly suggests that the uncertainty bounds for NEE that include the effects of u * th uncertainty presented here are unrealistically large.
The NEE uncertainty contribution of combined random and model error was small (≤ 30 gC m −2 a −1 ) by comparison (Table 6), < 10 % of annual NEE for each year.Modelinduced uncertainty was generally larger than random uncertainty, but the difference was more pronounced at night.This is because far more nocturnal data were removed by u * filtering than daytime data.The removal of observational data correspondingly increases the model error (because more model values are used) and reduces random error (because random error is only compounded across the cases for which there are observations).In contrast, during the day, there are less model data, and random error on individual data points is generally larger because of larger fluxes in combination with the heteroscedasticity of random error (see Fig. 13).
Annual NEE uncertainty due to random and model error was slightly greater for F c + S c than for F c alone due to additional random error to that for F c arising from the addition of S c (Fig. 13).This was largely nocturnally determined because the storage term was smaller and less variable during the daytime when fluxes were largest.The increased annual uncertainty of F c +S c was largely due to higher model uncertainty, which most likely reflects the propagation of random error to model uncertainty through its effects on non-linear parameter estimation.This was most pronounced nocturnally because of the lower signal-to-noise ratio of F c + S c relative to the daytime.
The method of calculating model-induced uncertainty may overestimate model error.It compounds observation-model data differences, but the observational data already contain random error so this necessarily inflates the propagated model uncertainty above errors associated with missing driver information or systematic measurement error.This problem is then partially propagated to the daytime because random error contributes to uncertainty in the nocturnally derived parameters of Eq. ( 5), which are then used to calculate the ER component of daytime NEE (Eq.6).Conversely, since the model-data comparison inevitably only occurs for periods where observation are available, it cannot be ascertained whether model and observational data would depart more substantially under conditions where the observational data are missing (which are often more extreme than the conditions for which observational data are available).The total error may thus be somewhat constrained.
The interdependence of model and random error also technically renders invalid the assumption of independence in Eq. ( 9).However, the effect is to increase rather than to decrease uncertainty, which as noted above was nonetheless small.While there are methods for separating model and random components (see Dragoni et al., 2007), this generally requires co-location of two instrument arrays.However, the daily differencing procedure we used is known to overestimate error by up to a factor of 2 (Billesbach, 2011;Dragoni et al., 2007) due both to potential wind-dependent temporal variations in source-sink strength (which may materially af- fect annual NEE estimates; Griebel et al., 2016) and because some signal is captured in the differencing procedure.
The random error in S c may potentially be reduced by redesign of the intake ports of the profile system.At present, air is drawn through a single port (via a filter) into the sampling tubes.This is similar to the design of other commercially available systems.Multiple spatially distributed intake ports may act to smooth the effects of random eddies, although it is imperative that this does not result in a large increase in the effective averaging time for the system, which, as noted by Finnigan (2006), may reduce the magnitude of the calculated storage term.
It should be emphasised that there are numerous sources of uncertainty that have not been quantified here.Perhaps the most important of these are systematic errors in the measurements themselves, which may be an extremely important source of true uncertainty (Lasslop et al., 2008).Thus, the uncertainties reported here for F c + S c should also not be formally interpreted as total uncertainty in the true sourcesink term, but as the uncertainty contributed by a subset of quantifiable errors.

Conclusions
We used a simple method to infer nocturnal advection from measurements combined with a simple and widely used empirical temperature response respiration model (modified for the effects of low soil moisture).While the nocturnal signalto-noise ratio was low in this relatively unproductive ecosystem, our results suggest that relatively accurate parameter estimation -and corresponding estimation of ER -is still possible.Even at our very flat site, approximately 40 % of nocturnal CO 2 flux underestimation was attributable to advection.Observation of reductions in storage at lower levels (within 8 m of the surface) in response to declining u * indicate that the most likely advective mechanisms are ter-rain drainage flows.High-resolution measurements of wind directions within the control volume would be invaluable for directly detecting the presence of these flows and are planned for this site.
For the ideal case of observational NEE measurements consisting of the sum of turbulent flux and profile-measured storage, we found that correcting for nocturnal respiratory underestimation (due to advection) using u * filtering reduced the cumulative annual CO 2 -carbon sink estimate by 10-20 %.Applying the u * filter across the diel cycle rather than nocturnally for this series had a small effect on cumulative annual NEE, very slightly reducing carbon uptake on average (though not consistently for all years).The bias occurred in the early morning and likely indicates residual drainageinduced horizontal advection after sunrise.
In the absence of storage estimates from a profile system, single-point estimated storage at this site generally biased cumulative annual NEE estimates towards efflux.Nocturnal u * filtering removed many of the erroneous nocturnal storage data, such that on average nocturnal estimates were not substantially biased.However, many large biases (towards efflux) occurred annually when the u * filtering was only applied nocturnally.Application of the diel u * filter reduced bias substantially.However, the single-point storage estimate greatly increased the error in measurements and is likely to result in large uncertainties where significant data gaps are present.
The use of F c alone to estimate NEE also resulted in consistent bias towards efflux, but the bias was on average larger, although more consistent, relative to the use of the singlepoint estimate.As above, the bias was far worse when u * filtering was only applied nocturnally.In fact, in this case the bias was larger than if the correction was not applied at all.This is likely to be the case for any site where the quantity of CO 2 diverted into storage overnight exceeds the quantity diverted into advection.

Figure 2 .
Figure 2. (a) Schematic and (b) photographic layout of profile gas analysis system.

Figure 3 .
Figure 3. Dependence of temperature on (a) u * and (b) time of day (temperature measurement for −0.05 m represents soil temperature, and SRT is surface radiant temperature).

Figure 4 .
Figure 4. (a) Akaike's information criterion (AIC, left-hand axis) and number of parameters (right-hand axis) for model with time-stepvarying R 10 component.(b) Comparison of temperature response of observational data with modelled data degraded with data-derived random error (dashed lines represent ±1 SD, observational points are shown for reference).(c) Interpolated R 10 estimates derived from observational data optimisation using 7-day time step (black line), and 95 % confidence interval for the mean (grey shading) for R 10 values extracted from 10 4 trials in which ER predicted from the empirical model was degraded with data-derived random error.

Figure 5 .
Figure 5. (LH axis) dependence of mean measured nocturnal carbon mass balance components (turbulent flux (F c ), storage (S c ) and F c +S c ) on friction velocity (u * ) ; (RH axis) air temperature at EC instrumentation height (36 m) and volumetric soil moisture content at 10 cm.Vertical lines denote u * th for both F c and F c + S c , as labelled.

Figure 6 .
Figure 6.(LH axis) profile system time series (date labels correspond to midnight) of CO 2 mole fraction (solid coloured lines); (RH axis) time series of u * measured at 36 m (dotted grey line).Note that horizontal dashed line marks u * th for F c .

Figure 7 .
Figure 7. Dependence of measured (F c , S c ), model-estimated ( ER) and inferred (A c ) mass balance components on friction velocity (u * , the grey shaded area represents the 95 % confidence interval for the sample bin mean of the inferred advection components).Vertical lines denote u * th for both F c and F c + S c , as per Fig. 5.

Figure 8 .
Figure 8. Dependence of measured (F c , S c ), model-estimated ( NEE) and inferred (A c ) mass balance components on friction velocity (u * , the grey shaded area represents the 95 % confidence interval for the sample bin mean of the inferred advection components).Vertical lines denote u * th for both F c and F c + S c , as per Fig. 5.

Figure 9 .
Figure 9. Dependence of measured storage components on friction velocity for individual layers.

Figure 10 .
Figure 10.Dependence of corrected storage components and Re −F c on friction velocity (dashed lines represent corrected storage estimates, stippled region represents difference between measured (blue) and corrected (dashed blue) 0-36 m storage, shaded regions represent 95 % CI for the u * bin mean).

Figure 11 .
Figure 11.Effects of storage addition and diel u * filtering on diel mean NEE dynamics (inset shows differences in NEE arising from application of diel versus nocturnal-only u * filter; calculated as nocturnal-filtered series subtracted from diel-filtered series).

Figure 12 .
Figure 12.Mean diurnal cycle of aggregate and component S c (LH axis, solid horizontal line represents mean S c ) and u * (RH axis, dashed horizontal line represents change-point-derived nocturnal u * threshold for F c ); note that day length as indicated by insolation is > 12 h due to missing data for several months during winter-spring 2013, slightly biasing the data towards a longer photoperiod.

Figure 13 .
Figure 13.Standard deviation of estimated random error (σ [δ]) as a function of flux magnitude for turbulent flux (F c ), turbulent flux plus profile-based storage estimate (F c + S c ) and turbulent flux plus point-based storage estimate (F c + S c_pt ).

Figure 14 .
Figure 14.Dependence of S c (including only data where u * > u * th ) on time after sunset (upper panel, dotted line is air temperature); cumulative percentage of total nocturnal S c observations (lower panel).

Figure 15 .
Figure 15.Effects of addition of storage term and application of daytime u * filtering on radiation use efficiency (RUE).

Table 3 .
Statistical ranking of candidate models for ER estimation (mode refers to the temporal fitting method, k is the number of parameters, RMSE is root mean square error, r 2 is coefficient of determination, AIC is Akaike's information criterion, AIC is the AIC difference relative to the best candidate, and w is Akaike weight.Note that additional candidate models were tested but were excluded if parameter estimation failed).

Table 4 .
Lower 95 % CI bound (µ − 2σ ), mean (µ), and upper 95 % CI bound (µ + 2σ ) of Gaussian PDF of u * th (derived from change-point detection of bootstrapped samples; see Methods), data percentile (i.e. percentage data excluded for each u * th ) and resulting imputed annual estimate of NEE.Note that (i) µ − 2σ set to zero if < 0 (e.g.F c + S c in 2013), (ii) respiration and light-response function analysis could not find a solution for F c + S c in 2013 when u * = 0.73 (insufficient data for robust statistical fit).

Table 6 .
Number of data available (after quality control and u * filtering), annual uncertainty due to random and model error for daytime and nighttime, and summed annual uncertainty (gC m −2 a −1 ) for diel u * -filtered F c and F c + S c .