Estimating local atmosphere-surface fluxes using eddy covariance and numerical ogive optimization

Regardless of study type, site topography, homogeneity and large-scale

flux estimates.Here, we present a novel comprehensive numerical scheme to identify and separate out advective contributions to exchanges in the surface layer.Comparison between the presented method and conventional methodology on observations of sensible heat, latent heat and CO 2 -fluxes from a number of sites suggests the presence of absolute flux thresholds at |Q SENS | = 30 Wm −2 , |Q LAT | = 16 Wm −2 and F CO 2 = 2.0 µmol m −2 s −1 marking clear shifts in the influence of advection.Above the thresholds, the relative difference of flux estimates δ remained fixed at δ = 5-25 % suggesting arguably negligible advection influence.Below the thresholds, however, relative difference rises to δ SENS = 51 %|88 %|225 % , δ LAT = 14 %|28 %|99 % and δ CO 2 = 41 %|83 %|521 % , where bracketed values are the 13.6th percentile, 50th percentile (the median) and the 86.4th percentile respectively, suggesting non-negligible relative influence of advection on low flux estimates.The thresholds thus serve as lower limits to local-scale flux resolvability by conventional methodology.The presented method is shown to allow for flux estimation during severe signal disruption and to yield fewer estimates for an enclosed gas analyzer during low-flux conditions suggesting the presence of a lower detection limit with this particular instrument setup, as well as a superiority of open path gas analyzers, in low-flux environments.Overall the notion of a dynamic and generally non-negligible overlap of advective and turbulent frequency-wise flux

Introduction
The eddy covariance (EC) technique allows for direct, continuous and non-invasive tower-based ecosystem-scale estimation of surface-atmosphere scalar fluxes by simultaneous sampling of atmospheric fluctuations of wind and scalars (Baldocchi, 2008).These characteristics, along with ease of operation, have promoted the widespread application of the technique in both short-term experiments and long-term monitoring network operations (e.g.FLUXNET, CarboEurope, EuroFlux, and AmeriFlux).
Reliable flux estimation in a local environment is often complicated by a number of issues relating to the large range of fluctuation-scales which drive fluxes (Stull, 1988).Fluxes driven by high-frequency fluctuations (turbulence) are inherently local in nature whereas fluxes driven by low-frequency fluctuations (advection) are associated with topographical forcing on the observed flow or large-scale meteorological phenomena, such as gravity waves, deep convection and large roll vortices (Lee et al., 2004).Traditionally the presence of a spectral gap (Stull, 1988) is assumed to exist between these contributions, allowing investigators to disentangle contributions simply by separating continuous observations into quasi-stationary intervals each yielding one flux estimate.However, the existence of a distinct spectral gap is unclear (Lee et al., 2004) and a growing body of work suggests that advection may often be non-negligible even for relatively flat sites.Furthermore studies have shown that the advection contributions are highly site-specific and characterized by significant uncertainty (Aubinet et al., 2010;Loescher et al., 2006;Yi et al., 2008).Hence, observations of atmospheric fluctuations are likely to reflect some degree of convolution between signals of local turbulent contributions and site/time-specific advection contributions.Introduction

Conclusions References
Tables Figures

Back Close
Full The importance of including vertical advection contributions in studies is debated.For instance, some studies suggest that inclusion may improve closure in energy and carbon-balance studies (Finnigan et al., 2003;Mahrt, 1998;Sakai et al., 2001;von Randow et al., 2002), which are otherwise often characterized by lack thereof (Wilson et al., 2002), while other studies suggest otherwise (Aubinet et al., 2010) and that inclusion may also enhance the variance of flux estimates, thus complicating interpretation in terms of local instantaneous surface fluxes (Kanda et al., 2004).Moreover, it has been commented that horizontal advection, which is typically assumed negligible, may become significant during certain conditions (Yi et al., 2008;Zeri et al., 2010), suggesting that proper inclusion of advection contributions require a more sophisticated array of EC systems forming a mass balance control volume in which all components of the mass balance may be adequately observed.Indeed the very nature of observations reflecting vertical advection of opposite sign relative to the turbulent contribution, as will be shown in this study, complicates the notion that fluxes should be of same sign regardless of incident eddy scales, unless vertical advection is taken to be only one part of a net advection contribution and hence perhaps balanced by a horizontal advective component.
Accordingly, we can distinguish between two principal applications of the EC technique: process-oriented studies in which fluxes are being linked to local biochemical processes for parametric insight into universal causal flux-relationships and up-scaled through numerical modeling efforts, and long-term net ecosystem-exchange studies in which the flux estimates are understood to be site-specific, applying only for the unique conditions of ecosystem heterogeneity, topography and large-scale meteorological flows experienced during the study.That is: studies in which investigators seek to disentangle signals to obtain the turbulence contribution only, and studies in which ideally both the vertical and horizontal advection contribution is included, and flux averaging times extended accordingly (Sakai et al., 2001), as a site-specific flux.This study will focus on the former and we will refer to the turbulence driven fluxes as locally meaningful fluxes, following Lee et al. (2004) For process-oriented studies, a number of typical approaches exist to estimate locally meaningful fluxes.These include: (1) tuning the flux averaging time to strike an appropriate balance between adequate sampling of the turbulent flux contribution while avoiding excessive inclusion of large-scale advection contributions (Sun et al., 2006); (2) optimizing site-location and discarding data reflecting critical wind-directions, to limit topographical forcing on the observed flow; (3) estimating vertical advection by performing profile measurements of fluxes on a single tower (Lee, 1998;Leuning et al., 2008) and filtering out observations reflecting excessive advection (Novick et al., 2014); (4) filtering observations based on co-spectral similarity with theoretical co-spectra assumed to represent local flux distributions for ideal site-conditions during a number of different surface-layer conditions (Hojstrup, 1981(Hojstrup, , 1982;;Hunt et al., 1985;Kaimal, 1978;Kaimal et al., 1972;Moore, 1986;Moraes, 1988;Moraes and Epstein, 1987;Olesen et al., 1984); (5) estimating the ideal turbulent contribution by matching the observed co-spectral peak with that of a theoretical distribution (Sorensen and Larsen, 2010).
While each method has its merits none is universally applicable and without its caveats.(1) In the absence of a distinct spectral gap between contributions, separating flux contributions by tuning the averaging time will inevitably fail.Moreover given the evolving nature of the natural flow a proposed spectral gap is likely to change in character over time, indicating that setting a fixed averaging time for an entire experiment inevitably causes some misrepresentation of fluxes.(2) As stated above, optimizing site location and acceptable wind-directions may in fact be a futile effort in removing the influence of advection.Furthermore it will serve to further support an existing bias in study-sites towards flat homogeneous sites.(3) Estimating the advection using observations on a single tower can produce uncertain estimates (Aubinet et al., 2010;Loescher et al., 2006;Yi et al., 2008) and requires a minimum of two EC systems, which due to instrument cost might not be a feasible option for many researchers.(4) Theoretical co-spectra depend on atmospheric stability, which in turn depend on the flux of sensible heat (Kaimal et al., 1972); a term we know to be af-Introduction

Conclusions References
Tables Figures

Back Close
Full fected by advection contributions.A similar circular dependency arises when attempting to scale observed co-spectra to a scalar flux so as to verify similarity with the theoretical co-spectra.Consequently, similarity of observed co-spectra with theoretical co-spectra may be impossible to evaluate in the presence of a strong advective contribution.Furthermore, in the limit of low absolute covariance (i.e.: small fluxes) large relative variance of the co-spectrum results in frequency-wise contributions which vary in sign, despite frequency-wise averaging, thus complicating an unambiguous comparison between observed co-spectra and theoretical co-spectra.While such cases could be treated as reflecting observations approaching the lower detection limit of the system, and discarded accordingly, they are ubiquitous for exchange studies in environments characterized by very low fluxes, such as during night-time or over sea-ice creating a demand for a new approach here.Finally, ensuring co-spectral similarity requires a number of site/system-specific empirical co-spectral corrections to account for high-frequency noise/dampening produced by the presence of the EC system in the observed flow as well as signal dampening in closed path systems (Aubinet et al., 2000;De Ligne et al., 2010;Kaimal, 1968;Massman and Ibrom, 2008;Moncrieff et al., 1997;Moore, 1986;Silverman, 1968), greatly complicating the approach.(5) Matching the co-spectral peak solves the issue of excessive scaling-offset mentioned above, but places an increased responsibility on the subjective evaluation of the investigator as well as retaining the circular dependency on stability and becoming increasingly ambiguous in the limit of low absolute covariance.
Here, we present a novel method for estimating locally meaningful atmospheresurface fluxes despite advection influences, using a single eddy covariance system and a numerical modeling scheme for Ogive optimization.Accordingly the method is called Ogive optimization.Ogive optimization makes no assumptions regarding optimal averaging time or the presence of a spectral gap, improves the flux estimates by also considering contributions in the high/low frequency ranges which due to instrument limitations and advective influences cannot be observed directly and allows for very high temporal resolution of flux evolution at less expense in terms of flux independence due Introduction

Conclusions References
Tables Figures

Back Close
Full to a modelling aspect.Furthermore, it requires no application of site/system-specific co-spectral corrections, is not affected by any circular dependencies for scaling of cospectra and decreases significantly the spectral flux ambiguity in the limit of low absolute fluxes.To evaluate the method, we applied it on EC observations of sensible heat, latent heat and CO 2 -flux at five sites covering different ranges of fluxes, ecosystemtypes and levels of topographical interference.Results were compared with conventional method both in terms of flux estimate yield and flux difference (i.e.advection influence) relative to flux strength.
2 Method and theory

Eddy covariance and spectral analysis
The theory of eddy covariance is well established (Baldocchi, 2008).Average surface fluxes of sensible heat, latent heat and CO 2 may be estimated over a large upwind area (Kljun et al., 2004;Kormann and Meixner, 2001) using tower mounted fast response instruments by: where Q SENS is the sensible heat flux, c p is the specific heat of dry air, ρ d is the mass density of dry air, e.g.w = w +w is the Reynolds decomposition of vertical wind-speed into its average ( w) and turbulent w components θ v is potential virtual temperature, Q LAT is the latent heat flux, L h is the latent heat of vaporization, c ds is the molar concentration of dry air mol dry m −3 , F CO 2 is the CO 2 -flux and r q = mol q mol 1 and negligible flux changes with height (e.g.dF CO 2 dz = 0).Here ρ is the air density, x j represents the three axes of observed flow, ς = {u, v, w} is the wind vectors, s = θ v , r q , r c is the scalars of interest and ξ = (ς, s) is the latter two combined.
Flux estimates (Eq. 1) may be decomposed into frequency-dependent contributions, called co-spectra Co ws (f ), between vertical wind-velocity w and the scalars of interest s for frequencies f .Deviations of observed co-spectra from theoretical co-spectra (Kaimal et al., 1972;Moore, 1986) can be linked to a number of issues including influence of the EC system on the flow, oscillations of the tower (or ship), advection, topographical forcing on the flow, etc., and is often used to filter out observations characterized by excessive non-turbulent influence (Novick et al., 2014).
Subsequently we may perform an Ogive analysis (Desjardins et al., 1989;Foken et al., 2006;Lee et al., 2004).The analysis requires the same basic assumptions and involves the cumulative summation of co-spectral energy, starting from the highest frequencies, The principal use of Ogives is to estimate the optimal observation time as the point of convergence of cumulative co-spectral energy to an asymptote (Berger et al., 2001).A feature, shown in this study to be common for observations in environments dominated by advection, yet scarcely described in the literature (Foken et al., 2006), is convergence to an extremum, reflecting a change in the direction of fluxes with lowerfrequency contributions.Such conditions may be understood as reflecting an equilibrium point between two conflicting flux contributions, e.g.turbulence-and advection-Figures

Back Close
Full In the ideal case (Fig. 1a), turbulent and advective flux contributions are separated by a spectral gap, allowing investigators to isolate the former simply by choosing an appropriate averaging time T 1 and using fast response instruments recording at frequency T 2 .Accordingly, the corresponding Ogive distribution is seen to converge to a stable flux estimate within T 1 (Fig. 1b).Note that this case reflects positive turbulent fluxes while advection may drive both positive and negative fluxes, shown here as a blue region of Ogive-divergence in Fig. 1b.
More typical, however, is the case (Fig. 1c) of overlapping contributions from advection, turbulence and site and instrument-specific noise/dampening.One way to strike a balance between adequate inclusion of the turbulent contribution and exclusion of excessive influence from the advective contribution is by tuning T 1 .Typically a fixed averaging time is set for an entire experiment (here 30 min is shown) and the inevitable flux errors are assumed, or tested (e.g.Novick et al., 2014), to be negligible.In the highfrequency end of the spectrum, instrument response limitations may prevent observation of the smallest scales of turbulence contribution (Here 10 Hz is shown, reflecting the Nyquist frequency for an instrument operating at 20 Hz).Furthermore noise may either be assumed, or visually inspected, to be negligible and site-specific dampening may at times be reduced by application of site-specific co-spectral corrections, called transfer functions (Aubinet et al., 2000;De Ligne et al., 2010;Kaimal, 1968;Massman and Ibrom, 2008;Moncrieff et al., 1997;Moore, 1986;Silverman, 1968).Accordingly the Ogive distribution indicates negligible influence on the flux estimates for a 30 min averaging time and an instrument response time of 20 Hz (Fig. 1d).Note that the influ-Introduction

Conclusions References
Tables Figures

Back Close
Full ence of dampening and noise on the Ogive distribution occurs in the reverse (Fig. 1d) relative to co-spectral space (Fig. 1c).
Observations reflecting excessive advection, relative to the turbulent contribution, (Fig. 1e) are typically discarded.This is because strong relative advection results in non-negligible flux contribution to the overall estimate and further obstructs any efforts to separate contributions by tuning the averaging time (Fig. 1f).The use of the term relative in this context refers to the fact that an identical problem can arise despite modest advection when estimating fluxes in a low-flux environment.Flux estimation in such environments are often further complicated by a high ratio of co-spectral variance to actual turbulent flux contribution, prohibiting unambiguous evaluation of similarity between observed co-spectra and theoretical co-spectrum distributions, as well as proper estimation of the co-spectral peak (Sorensen and Larsen, 2010).

Formation of averaging intervals
In order to fulfill the stationarity requirement described in Sect.2.1, continuous observations are typically subdivided into averaging intervals.Averaging interval time T 1 has conventionally been assumed constant, based on the requirements that T 1 should be long enough to reduce random error (Berger et al., 2001;Lenschow and Stankov, 1986) and short enough to avoid non-stationarity associated with advection (Foken and Wichura, 1996;Vickers and Mahrt, 1997).However, as noted, tuning of T 1 will not generally allow for separation of turbulent and advective flux contributions.Here, we propose a method that makes no assumption regarding the presence of a spectral gap.Instead we require averaging time to be as long as necessary while ensuring stationarity of the local processes, irrespective of the temporal evolution of advection contributions.Hence the question of averaging time becomes less a question of separating contributions, and more a question of identifying, in advance, the dynamic characteristics of the flux evolution in the local environment.
The following is an iterative scheme for the formation of averaging intervals based on basic data quality requirements.Data collected during a field-experiment is consid-21397 Introduction

Conclusions References
Tables Figures

Back Close
Full Here α is set according to desired temporal resolution of flux-estimates.For this study we use site-specific settings of α = 5min, α = 15 min and the conventional α = 30 min to strike a balance between desired temporal resolution and computational cost of running the full Ogive optimization method.The minimum dataset length is chosen to be 10 min for the Ogive function to yield statistically representative estimates of the scales of turbulence-driven fluxes and the maximum dataset length is chosen to be 60 min to ensure approximate stationarity of the local turbulence-driven fluxes.
Using an iterative bisectional algorithm for enhanced computational speed, combinations of T m (i ) and T w (j ) are evaluated with regard to a number of basic quality assessment criteria to obtain the longest dataset around T w (j ), which passes the assessment criteria.These include absence of instrument diagnostics errors, absence of long data gaps, favorable mean wind-direction and reasonably narrow range of wind-directions (≤ 60 • ).Minor spiking (≤ 1 %) is corrected based on the median of surrounding datapoints and the dataset is discarded otherwise (for spiking > 1 %) according to Vickers and Mahrt (1997).Other quality assessment includes the requirement that momentum fluxes be negative for non-static surfaces.Evidence to the contrary would imply a disconnection between the upwind surface processes and the point of observation on the tower, a condition typically associated with low wind conditions.However, momentum fluxes may be affected by advection, implying a circular dependency.As sign, but not strength, of the turbulent momentum flux is relevant we may simplify by calculating an Ogive distribution Og wU based on the momentum flux co-spectrum Co wU , where U is the wind velocity, and verifying that Og wU < 0 at the mid-range natural frequency of f m = 0.1 Hz.The choice of f m reflects a spectral region least impacted by noise, damp-Introduction

Conclusions References
Tables Figures

Back Close
Full ening and advection.Typically around 4 estimates of T w (j ) are evaluated for each T m (i ) before the optimal T w (j ) is determined.Lastly, signals with very rapid evolutions such as transient signals in dynamic systems like eddy covariance observations may undergo abrupt changes such as a jump, or a sharp change in the first or second derivative, associated with observational interference by e.g.electrical interference or instrument error.These are referred to as dropouts and discontinuities in Vickers and Mahrt (1997).Global transforms, like the Fourier transform, are usually not able to detect these events.In contrast, Wavelet transforms such as the Haar transform, permit a localized evolutionary spectral study of signals, thus allowing for detection of subtle signal discontinuities which lead to semi-permanent changes (Lee et al., 2004;Mahrt, 1991;Vickers and Mahrt, 1997).Observations are sorted in three categories (good, soft flag, hard flag) according to the presence and severity of such events (Vickers and Mahrt, 1997).Conventionally these events are thought to preclude flux estimation on the basis of stationarity violations, and such datasets (hard flags) are discarded when applying the conventional methods.However, as will be shown, we find in this study that the Ogive optimization method allows for convincing flux estimation in many cases of soft and hard flags, suggesting superiority in applicability relative to existing methods.In this study we perform the Haar analysis for the optimal window-size T w (j ) and filter out any conventional method flux estimates accordingly while the Ogive optimization flux estimates are allowed, pending visual inspection.

Mass Ogive calculation
As noted, simply tuning the averaging time will be insufficient in achieving an accurate estimate of turbulent fluxes.What is needed is a second parameter to control the low-frequency contribution alone.Subtracting a running mean from observed signals, as opposed to the conventional linear detrending, allows for enhanced filtering of low-Introduction

Conclusions References
Tables Figures

Back Close
Full frequency contributions alone (Sakai et al., 2001;Mcmillen, 1988).Consequently some combination of averaging time and running mean resolution might allow for filtering out of advection contributions while retaining turbulent contributions.Here we visualize this concept by calculating co-spectra, and corresponding Ogives, for a very large number of data pertubations and deriving a map of the resulting Ogive density pattern (Fig. 2a).
The figure illustrates the density of 10 000 individual Sensible heat-flux Ogives based on the following data perturbations: 50 linear increments on the averaging time axis between 5 min and the maximum time available (60 min in this example) and 200 increments for the running mean resolution in the range of one second to half the length of the dataset in question.Increment density for the latter is shifted towards the lowresolution range (large running mean interval) of the running mean, which is of greater interest to us.The standard 30 min linear detrended Ogive is marked in red.
What is clear immediately in this particular example is the strong consistency between individual Ogive representations.This suggests that the fluxes are very well defined for this particular period with an actual flux around −55 Wm −2 following the convergence to a horizontal assymptote.A classic Ogive shape.The flux estimate for a regular 30 min linear detrended dataset (red) appears representative of the overall Ogive pattern as well.The presence of haze below −60 Wm −2 suggests that a small part of the pertubed states are affected by advection.This makes sense as the linearly detrended 60 min dataset illustrated in Fig. 2c appears affected by some large-scale oscillations, which for low-resolution running mean permutations are likely to have a significant impact.The influence of large-scale oscillations is supported by the Haar analysis, which has soft-flagged the temperature signal accordingly.

The model and the optimization method
It is clear from dampening and advection influences.This was illustrated in the Fig. 1b, d, f.While such a region is clearly evident in Fig. 2a for f > 10 −1.3 the example also illustrates the real advantage of performing a large number of data perturbations and deriving a density map of possible solutions: the most likely Ogive distribution of the observation in question stands out as a very well defined pattern which for this particular observation allows us to extend our understanding of the ongoing fluxes all the way to f = 10 −3.6 , or 60 min.While the co-spectral peak method (Sorensen and Larsen, 2010) bases its flux estimation on one point within this representation (i.e. the peak) we base our flux estimate on the entire range to enhance the certainty.
To describe the most likely flux resulting from a given Ogive density pattern we apply the generalized co-spectral distribution model (Lee et al., 2004) where a number of parameters are tuned to change the appearance of the co-spectral distribution: A 0 is a normalization parameter, µ is a broadness parameter controlling the shape of the spectrum, f is the natural frequency, f x is a horizontal offset of the distribution and m = 3 4 is a constant describing co-spectra characterized by a −7 3 power law in the inertial subrange.Subsequently an Ogive distribution is calculated using (Eq.2).We set A 0 = 1 and instead scale the low-frequency end-point f 1 (equivalent to the averaging time T 1 ) of the Ogive distribution to a variable flux parameter F 0 so as to allow for more direct flux control.This is particularly convenient when formulating reasonable limits on fluxes for the optimization algorithm described below.
One important aspect considered is the concept of local fluxes that cannot be observed directly.The problem may arise in the low-frequency range as over/underestimation of covariance due to inclusion of advection or the use of inadequate averaging times, and in the high-frequency range as under-estimation of covariance due Figures

Back Close
Full to inadequate sensor frequency, attenuation and distortion by both the spatial averaging of the sensors and the sampling and filtering of the sensor electronics.This is illustrated in Fig. 3a.Actual fluxes, represented by an ideal theoretical co-spectrum (red line) is shown alongside a corresponding Ogive (black line).In the case of insufficient observation time (here 20 min) and observational frequency (here f nyquist = 5 hz) the missing range of observed fluxes can be illustrated as an orange area below the co-spectrum and as equivalent blue ranges in the negative and positive Ogive axes respectively.The corrected flux is shown as a dashed black line and is derived as , where F 0 is the uncorrected flux.Note that F HF1 is subtracted as it is of opposite sign relative to F 0 .Secondly, because theoretical models are empirical representations, verified only within a certain frequency range, it becomes tempting to investigate how much, if any, of the fluxes are being left out by such restriction in observational range, assuming that the model is valid outside this frequency range.
The consequent extrapolation of model results beyond actual observed frequencies is illustrated in Fig. 3b and the corrected flux (dashed black line) may similarly be derived In practice we include both these considerations by combining F 0 and all lowfrequency contributions into one parameter: F LF = F 0 + F LF1 + F LF2 and adding a highfrequency offset F HF = − (F HF1 + F HF2 ).In total we are thus left with four tunable parameters: F LF , F HF , µ and f x , for which the final model-estimated flux is: where O 2 is adopted as shorthand mathematical notation for Ogive optimization.Our goal is to tune the parameters of the Ogive model to achieve an optimal fit to the density map.That is, to find the solution, which follows optimally the strongest densities in the density map (Fig. 2a) using a fast local optimization algorithm and a slower, but robust, darwinian evolutionstyle global optimization algorithm called Differential Evolution (Storn and Price, 1997).
The optimization is performed for 18 frequency intervals, seen as intervals of optimization weights in Fig. 4a, and the final solution is chosen by subjective visual inspection (Fig. 4b).A number of objective approaches to visual inspection were investigated but none were found, at this stage, to best subjective visual inspection.Further development of this aspect is of continued interest as subjective visual inspection, aside from being a very time consuming process, may result in personal bias on final flux estimates.
One intriguing consequence of including a modeling and optimization aspect is that the inevitable occurrence of overlapping data-intervals does not relate linearly to interdependency of successive flux estimates, suggesting that the Ogive optimization approach allows for very high temporal resolution of flux evolution at less expense in terms of flux independence.

Field sites
To evaluate the Ogive optimization method, five sites, reflecting different environments in terms of both ecosystem, topography and flux strengths, are investigated.

High-flux environment
The Abisko field-site is located in Stordalen ( 68 an LI-7500 open path gas analyzer (LI-COR ® , Lincoln, NE, USA) on a boom extending towards the southwest at 2.50 m height, with a 0.43 m horizontal offset along the boom and a slight tilt of the instrument relative to the vertical plane to allow water dripping.Raw data were logged at 10 Hz.Observations reflecting wind origins along the boom axis were filtered out to limit flow distortion.Flux estimates were evaluated in intervals of: α RIMI = 30 min.

Intermediate-flux environment
The RIMI site is an active FLUXNET site located in a large flat homogeneous grassland area (55 LI-7500 open path gas analyzer (LI-COR ® , Lincoln, NE, USA) mounted on the same boom at heights of 2.2 m and 2.1 m respectively extending from the side of a 10 m mast, with a horizontal offset along the boom of 0.40 m.Observations reflecting wind origins along the boom axis were filtered out to limit flow distortion.Flux estimates were evaluated in intervals of: α ABI = 30 min.

Low-flux environment
Young Sound is the entrance of a 7 km wide fjord in NE Greenland characterized by thick fast sea-ice within the fjord and an ice-free polynia at the mouth of the fjord (Rysgaard et al., 2003).Continuous EC observations were conducted at three sites within • relative to the horizontal plane at a height of 2.7 m.The sonic anemometer and gas analyzer were displaced horizontally by 0.4 m in orthogonal alignment to the prevailing along-sound wind direction, so as to limit the instrument flow distortion and temporal offset between simultaneous signals.
In addition to filtering for tower based flow distortion, observations from the shoreadjacent DNB site reflecting wind-directions associated with the shoreline were likewise filtered out due to anthropogenic interference.For the static tower a Gill Windmaster pro sonic anemometer (Gill Instruments ® , Lymington UK) was mounted at a height of 3.66 m and an enclosed LI-7200 gas analyzer (LI-COR ® , Lincoln, NE, USA) was mounted with a 65 cm inlet tube terminating directly under the sonic anemometer.Seaice thickness and snow-cover was approximately 110 cm / 75 cm respectively for the ICE1 and DNB sites, and approximately 25 cm / 20 cm respectively for the POLY1 site.Average air temperature climbed steadily from −35 • to −15 • throughout the period.
As such all sites were expected to be characterized by significantly smaller turbulent fluxes relative to the Abisko and the RIMI sites while simultaneously being subjected to varying degrees of large-scale topographical advection due to their locations in a fjord surrounded by mountains.Flux estimates were evaluated at the three sites for the following intervals: α POLY1 = α ICE1 = 5 min and α DNB = 15 min.Introduction

Conclusions References
Tables Figures

Back Close
Full

Instrument corrections
A number of instrument-specific corrections are needed to adjust for instrument-bias.
For the sonic anemometers (Gill R2, Gill R3, Gill Windmaster Pro and METEK USA-1) these include: An empirical angle of attack correction (Nakai and Shimoyama, 2012) known to enhance fluxes over aerodynamically smooth surfaces.Humidity and crosswind corrections (Liu et al., 2001;Schotanus et al., 1983) accounting for the effect of humidity fluctuation and wind components normal to the measurement path, on estimation of the air temperature by sonic anemometer.Coordinate rotation, linear de-trending and iterative de-spiking of raw data is performed according to Vickers and Mahrt (1997).
We convert all observations to mixing ratios (Burba et al., 2012) using the WPL correction when necessary (Sahlee et al., 2008;Webb et al., 1980) as recommended by Ibrom et al. (2007).The need for instrument heating corrections (Burba et al., 2008) associated with operation of the open path LI-7500 in a cold environment (Daneborg, POLY1 and ICE1) is alleviated by using the newer LI-7500A with a "cold"-setting correcting observations down to −25 • (Burba et al., 2011).Sites featuring the LI-7500 (Abisko and RIMI) never reached sufficiently cold temperatures to warrant instrument heating corrections during this experiment.Temporal offset between sensor signals is typically corrected based on a maximum cross-covariance analysis (Berger et al., 2001;Fan et al., 1990).In the limit of low absolute covariance, however, actual temporal offset may be obscured by secondary cross-covariance optima.Furthermore, the use of open path instruments precludes the assumption of a fixed offset value.Consequently we narrow the search for optimal cross-correlation to a range described by incident horizontal flow and the specific geometry of each EC system, a methodological correction, which proved to be robust during automatic operation.Introduction

Conclusions References
Tables Figures

Back Close
Full 3 Results and discussion

Examples of Ogive optimization performance
In the following we describe several typical cases observed and the associated performance of the Ogive optimization method.(1) Commonly observed is the near-absence of advection and strong similarity between the Ogive density pattern, the 30 min linear detrended Ogive and the modelled Ogive.This is illustrated in Fig. 5 for a case of sensible heat flux.Not counting the high-frequency component associated with extrapolation of model results, seen here to contribute ≈ 5 Wm −2 to the overall modelled sensible heat flux, the standard 30 min linear detrending approach will suffice to provide the turbulent flux estimate.( 2) An example of non-negligible advective contribution on the low-frequency part of the Ogive spectrum is shown in Fig. 6.The Ogive optimization method is seen to separate the turbulent and the advective contributions completely, yielding the locally meaningful turbulent flux, regardless of the sign of the advective contribution relative to that of the turbulent contribution.( 3) An example of ambivalence caused by a bimodality in the Ogive density pattern is shown in Fig. 7.
Such cases indicate that fluxes are changing within the data-window in question.The Ogive optimization method is seen to appropriately yield the turbulent flux contribution with the strongest density influence.In the case of Fig. 7a the method is seen to also filter out any advection contributions.(4) The inadequacy of applying a fixed averaging interval for flux estimation becomes readily apparent in Fig. 8. Here, both Ogive density patterns are seen to reflect a gradual evolution in the Ogive flux pattern with increasing averaging time.In Fig. 8a the standard 30 min averaging time does not produce the well-defined Ogive flux pattern that is seen to emerge around 60 min averaging time.
In contrast, the standard 30 min averaging time is seen to be too long and also to increasingly reflect advective interference in Fig. 8b.In both cases the Ogive optimization method identifies the appropriate flux estimate, whereas the standard 30 min linear detrending method fails on account of instationarity.(5) The Ogive optimization method is seen in Fig. 9  quality reflecting e.g.instrument failure or electronic interference.( 6) During conditions of strong high-frequency dampening caused by the use of a closed path instrument, the high-frequency bound on optimization is automatically shifted towards lower frequencies to avoid influence of the dampened frequencies during optimization.As illustrated in Fig. 10 this allows for accurate description of the missing high-frequency flux signal.

Comparison of Ogive optimization and the conventional method
The difference in flux estimates of the standard 30 min linear detrending approach and the Ogive optimization method is associated with both the inclusion/exclusion of advective contributions, the inadequacy of the fixed averaging interval in capturing a representative estimate of the fluxes and the extrapolation of modelled Ogives into un-observable high/low-frequencies.The relative flux difference δ is evaluated within i = 10 intervals of absolute Ogive optimization flux estimates F O 2 i as the standard deviation of difference in flux estimate relative to the mean absolute Ogive optimization flux estimate within respective intervals: As hypothesized in Sect.2.2, the overall relative flux difference is seen to be very high for small absolute flux estimates, peaking at δ SENS = 51 %|88 %|225 % , δ LAT = 14 %|28 %|99 % and δ CO 2 = 41 %|83 %|521 % for the lowest absolute flux estimates, where bracketed values are the 13.6th percentile, 50th percentile (the median) and the 86.4th percentile respectively.Moving towards higher absolute flux estimates the relative difference is seen for all three flux types to drop and level off to a near-stable range of δ = 5-25 %.The onset of near-stable relative difference occurs approximately following the absolute flux thresholds marking clear shifts between non-negligible advection contributions on one side and plausibly negligible advection contributions on the other.For latent heat flux, the location of a universal threshold is less clear with the overall pattern suggesting a shift in advection influence at Q O 2 LAT = 16 Wm −2 but all low-flux sites (RIMI, Daneborg, POLY1 and ICE1) suggesting a shift in advection influence as low as Depending on perspective and the character of observed fluxes at a particular site the described thresholds may either serve as an indicator of a lower limit to localscale flux resolvability by the standard 30 min linear detrending approach, or as an argument for the application of enhanced flux estimation techniques such as the presented method.In the presented cases the consequences are illustrated well by the histograms of the different sites (Fig. 11).Although the location of the flux threshold is a bit unclear for latent heatflux, estimation of locally meaningful fluxes at the three sea-ice sites Daneborg, POLY1 and ICE1 is essentially impossible without accounting for advection contributions.The same applies for sensible heat flux at the Abisko lake site, latent heat flux at the grassland site RIMI and CO 2 -flux at both the Abisko lake site and the RIMI site.Note that only the Abisko Fen environment showed a dynamic range in excess of F O 2 CO 2 = 3.5 µmol m −2 s −1 and that most flux estimates from RIMI are from the morning or late evening/night, which explains the range of relatively small fluxes.Introduction

Conclusions References
Tables Figures

Back Close
Full For many flux estimates the vertical wind speed signal, or the scalar signal in particular, are non-stationary to the point of prohibiting a flux estimation using traditional methodology.Hence the presented method may also provide a greater number of flux estimates.This is shown in Table 1 to generally be true for the Abisko site and the Daneborg site, both of which characterized by degraded signal quality at times.Sites RIMI and POLY1 are inconclusive in this respect and the conventional method appears superior in the case of ICE1.The latter may be related to the very low fluxes observed for this site (Fig. 11) suggesting perhaps the presence of a lower detection limit when using the Ogive optimization method and an LI-7200 enclosed gas analyzer within the Shifts in flux direction towards low-frequency advective contributions (e.g.Fig. 6b) were found to be common in this study.To our knowledge such occurrences are not described by any existing theoretical framework, indicating a puzzling caveat to current theory.The occurrences greatly complicates the notion that fluxes should be of same sign regardless of incident eddy scales, unless vertical advection is taken to be only one part of a net advection contribution and hence perhaps balanced by a horizontal advective component.Indeed the horizontal advective component has been shown to be significant during certain conditions (Yi et al., 2008;Zeri et al., 2010), despite typically being assumed negligible.The finding indicates that further investigation of the interplay between advective contributions, and their influence on turbulent flux estimates, is necessary.

Conclusions
The presented Ogive optimization method has been shown to successfully separate local/non-local flux contributions as well as enhance flux estimation by both investiga-Introduction

Conclusions References
Tables Figures

Back Close
Full Lastly, low-frequency shifts in flux direction were found to be common in this study.To our knowledge such occurrences are not described by any existing theoretical framework.Based on studies indicating non-negligible horizontal advection during certain conditions (Yi et al., 2008;Zeri et al., 2010) we hypothesize a more intricate balancing interplay between vertical and advective flux contributions which, if confirmed, suggests the need for more sophisticated EC system arrays if advective contributions are to be accurately included (i.e.: site-specific studies).If exclusion of advection is desired (i.e.: universal process oriented studies), the presented method should be unaffected by these questions.
Future work might involve addressing the computational demands of the method (currently 2-3 min processing time for each flux estimate in a MATLAB-native parallel processing scheme), implementing automatic selection of the optimal frequencyinterval for Ogive optimization as opposed to subjective evaluation by visual inspection, implementing model capability towards improved handling of tower/ship oscillations and the presence of gaps in the data (Barnhart et al., 2012), as well as making the Ogive optimization approach, and all presented aspects thereof, available in open source MATLAB code form.Future studies might include an extensive parametrical data-analysis of well-defined (co)spectra (e.g.Fig. 5) for a wide range of observational conditions to expand and consolidate our current understanding and formulation of theoretical universally applicable (co)spectra (Hojstrup, 1981(Hojstrup, , 1982;;Hunt et al., 1985;Kaimal, 1978;Kaimal et al., 1972;Moore, 1986;Moraes, 1988;Moraes and Epstein, 1987;Olesen et al., 1984), as well as application of the method in a sophisticated multi-tower setup to investigate more closely the influence of both vertical and horizontal advection on the results presented.Introduction

Conclusions References
Tables Figures

Back Close
Full As described in Sect.2.4.2 our goal is to tune the four final model parameters F LF , F HF , µ and f x to achieve the optimal fit between a modelled Ogive and the Ogive density map (e.g.Fig. 2).The process is called optimization and involves the following steps: 1.A random guess of parameters is made within a set of reasonable bounds.The speed and accuracy of any optimization method involving preset bounds, depend greatly on the reasonable choice of these bounds.Here we set 0.05 < µ < 1 and −2 < log (f x ) < 1. 2. Based on the initial guess, a local solution is determined using the MATLAB function fminsearchbnd which is a nelder-mead polytope direct search optimization algorithm.The algorithm is fast for problems of low dimensionality such as ours, but not certain to converge to a global solution.The goal is to perform a rough, but fast, improvement of the random guess to limit processing time for the next step, which is far more labor intensive.
3. Based on the local optimization produced by fminsearchbnd, a global solution is determined using the Differential Evolution (DE) algorithm (Storn and Price, 1997).Differential evolution is a simple and reliable evolutionary population-based search technique, which has been successfully applied on a wide range of problems in a variety of scientific fields (Mallipeddi et al., 2011) Full Clearly the problem, as posed to the optimization algorithm, lacks an element mirroring our basic sense of intuition.Different schemes to address this issue were investigated, though none proved robust enough at this time to compete with basic subjective evaluation during visual inspection.Note the gradual decrease in optimization weighting of high-frequency Ogive density (Fig. 4a), which has been added to limit any influence of high-frequency noise and dampening during the optimization.The high-frequency limit of the fitting interval is furthermore allowed to move to lower frequencies for closed-path instruments to account for excessive dampening of the high-frequency end of the spectrum often observed with this type of instrument.Full and below which advection influence is non-negligible with average relative differences climbing as high as δ SENS = 51 %|88 %|225 % , δ LAT = 14 %|28 %|99 % and δ CO 2 = 41 %|83 %|521 % where bracketed values are the 13.6th percentile, 50th percentile (the median) and the 86.4th percentile respectively.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dry mixing ratio of humidity and CO 2 concentration scalars respectively.Equation (1) is valid only during adherence to a number of assumptions which ideally should be met during field operation.Principal among these are Figure 1 illustrates a number of observational situations showing examples of how advection could adversely affect our ability to capture local fluxes.In the figure, situations are shown using both co-spectral and Ogive plots.In the ideal case (Fig.1a), turbulent and advective flux contributions are separated by a spectral gap, allowing investigators to isolate the former simply by choosing an appropriate averaging time T 1 and using fast response instruments recording at frequency T 2 .Accordingly, the corresponding Ogive distribution is seen to converge to a Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ered continuous with end-points T min and T max , despite the presence of gaps.A range of possible subsets are preliminarily determined based on a linear series of interval midpoints T m Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 2 that answering the overall question of most likely flux, requires the application of an Ogive model.Particularly because not all Ogive density maps are as well defined.The basic premise in our model solution is that a region exists in the midto-high frequency range of the Ogive representation which is least impacted by noise, Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | . A number of local and global optimization techniques were investigated in terms of accuracy and speed.The final steps taken in optimizing the model with respect to the Ogive density map are explained in detail in Appendix A. The steps include optimization of a random guess within reasonable parameter bounds 21402 Discussion Paper | Discussion Paper | Discussion Paper | • 21.248 N, 19 • 3.02 E), a mixed mire 10 km East of Abisko in subarctic Sweden.The measuring mast is situated at the edge of a minerotrophic fen dominated by sedges.Wind patterns consistently alternate between an upwind fen-environment signal towards the west and an upwind lakeenvironment signal towards the east.Hence we treat the observations separately as fen-and lake-sites respectively.Continuous EC observations were conducted in the period 2 July to 1 August 2012.Site instruments include an R3 sonic anemometer (Gill Instruments ® , Lymington UK) mounted on top of the mast at 2.92 m height and Discussion Paper | Discussion Paper | Discussion Paper | the fjord system in the period 20 March 2012 to 27 April 2012.Two separate field-stations, one static and one mobile, were used at three different locations (ICE1, Discussion Paper | Discussion Paper | Discussion Paper | POLY1 and DNB).ICE1 (74 • 18.576 N, 20 • 18.275 W) was located 2 km from the coastline during the period 20-27 March and DNB (74 • 18.566 N, 20 • 13.998 W) was located approximately 200 m from the coastline during the period 30 March to 27 April.POLY1 (74 • 13.883 N, 20 • 07.758 W) was located at the mouth of the sound close to the ice-free polyna region during the period 24-27 March.For the mobile tower (POLY1 and DNB) a METEK USA-1 sonic anemometer (METEK ® , Elmshorn, Germany) was mounted at a height of 3.1 m and an LI-7500A open path gas analyzer (LI-COR ® , Lincoln, NE, USA) was mounted at an angle of 70 Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | to allow for trustworthy flux estimates despite severely degraded signal 21407 Discussion Paper | Discussion Paper | Discussion Paper | [ ] i signify flux estimates native to interval i of the equivalent absolute Ogive optimization flux estimates F O 2 i .Estimates of relative flux difference δ is shown logarithmically in Fig. 11 for all three scalar flux types, at all 6 observation sites and for all 10 intervals of the respective Ogive optimization flux ranges.Outliers have been excluded from the flux ranges shown in the bottom of the figure to ensure a minimum of 3 flux estimates within the largest absolute flux-estimate bin and resolution of the resulting δ estimates have been doubled by spline interpolation.The median relative difference is shown (red line) along with standard deviation (light gray area) and 25-75 % percentile range (dark gray area).Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and F CO 2 < 0.1 µmol m −2 s −1 ranges.No similar features indicate the presence of a lower detection limit for the LI-7500 (Abisko & RIMI) and the LI-7500A (Daneborg & POLY1) gas analyzers based on these numbers, suggesting superiority of open path instruments in very low flux environments.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |tion of a large range of averaging times and running means, and extrapolation of model results into otherwise unobservable frequency ranges.The presented method makes no assumptions concerning appropriate averaging time or the presence of a spectral gap, does not require the application of site/system-specific spectral corrections and allows for very high temporal resolution of flux evolution at less expense in terms of flux independence due to a modelling aspect.Comparison between conventional methodology and the presented method on observations of sensible heat, latent heat and CO 2 -fluxes from a number of sites suggests the presence of an overall absolute flux relative difference δ remains fixed at an arguably negligible δ = 5-25 % and below which relative difference rises to δ SENS = 51 %|88 %|225 % , δ LAT = 14 %|28 %|99 % and δ CO 2 = 41 %|83 %|521 % for the lowest absolute flux estimates, where bracketed values are the 13.6th percentile, 50th percentile (the median) and the 86.4th percentile respectively.This suggests a non-negligible relative influence of advection on low flux estimates, essentially preventing local-scale estimation of small fluxes using traditional methods.The method has furthermore been shown to allow for flux estimation despite severe signal disruption, resulting in increased flux estimate yield during such conditions, while simultaneously yielding less flux estimates for an LI-7200 enclosed gas analyzer during very low-flux conditions suggesting the possible presence of a lower detection limit in the|Q SENS | < 25 Wm −2 , |Q LAT | < 3 Wm −2 and F CO 2 < 0.1 µmol m −2 s −1ranges with this particular instrument setup, as well as a superiority of open path instruments in low-flux environments.The study suggests favourable application of the Ogive optimization method in most environments, particularly in environments characterized by small fluxes such as during night-time or over sea-ice.Overall the notion of a dynamic and generally non-negligible overlap of advective and turbulent flux contributions is confirmed, suggesting the inevitable indiscriminate inclusion/exclusion of both when setting a fixed averaging time according to conventional methodology.Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |and survivability within the probability space of a multivariate problem.Mutations are governed by predefined mathematical relations, called strategies, which depend on crossover probability CR = [0, 1] and differential weight F ∈ [0, 2], and survivability relates to the change in probability (i.e. the sum of Ogive density below a given Ogive model solution) between two generations.The performance of the optimization algorithm varies with each problem and depends greatly on the choice of strategy and algorithm parameters.For the purpose of optimizing the algorithm performance a number of cases were investigated using various strategies and a large number of parameter variations resulting in the application of the strategy called DE/best/1/exp and parameters NP = 40, CR = 1 and F = 0.8 for a maximum of 100 iterations.4.Often optimizing a smaller subset of the problem is an advantage, particularly during advective conditions which persist despite data-perturbation.One such case is shown in Fig.4b.The Ogive density map reflects an advective flux contribution of opposite sign relative to the turbulent flux contribution.Solutions derived for a range of 18 frequency intervals, shown in terms of optimization weighting in the lower pane, can be seen in the upper pane as green lines.All solutions based on frequency intervals with lower bound before or after the Ogive density peak are seen to underestimate the actual undisturbed turbulent flux (marked in blue).
Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .Figure 8 .Figure 11 .
Figure 1.The figure illustrates a number of typical observational situations in terms of cospectra (top row) and Ogives (bottom row).Shown are the turbulent fluxes (red) and advective/noise/dampening components (blue).(a, b) In the ideal case a spectral gap separates advection and turbulent fluxes and a fixed averaging interval allows for disentanglement of contributions.(c, d) In the typical case advection and turbulent fluxes overlap resulting in some flux misrepresentation often assumed to be negligible.(e, f) In the case of strong relative advection, the advection contribution cannot be neglected and the observation is typically discarded.
The bounds for F LF and F HF are a bit more complicated.If > 80 % of the summed density map is located on, say, the positive side (suggestingF LF > 0), bounds on F LF and F HF are set as 0 < F LF < |2R + | and − |R + | < F HF < 0,where R + is the 95th percentile range of the positive side of the density map.Reverse bounds are applied if > 80 % of the summed density map is located on the negative side (− |2R − | < F LF < 0 and 0 < F HF < |R − |), and both cases are run if neither side contains > 80 % of the summed density map.

Table 1 .
Number of flux estimates from the conventional method (N 30 min ), the Ogive optimization method (N O 2 ) and the number of combined pairs of estimates (N Both ) used to determine the relative flux-estimate differences illustrated in Fig.11.