Estimating surface fluxes using eddy covariance and numerical ogive optimization

. Estimating representative surface ﬂuxes using eddy covariance leads invariably to questions concerning inclusion or exclusion of low-frequency ﬂux contributions. For studies where ﬂuxes are linked to local physical parameters and up-scaled through numerical modelling efforts, low-frequency contributions interfere with our ability to isolate local biogeochemical processes of interest, as represented by turbulent ﬂuxes. No method currently exists to disentangle low-frequency contributions on ﬂux estimates. Here, we present a novel comprehensive numerical scheme to identify and separate


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 (e.g., Baldocchi, 2008). These characteristics, along with ease of operation, have promoted the widespread application of the technique in both short-term experiments and longterm monitoring network operations (e.g., FLUXNET, Car-boEurope, 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 are associated with e.g., topographical forcing on the observed flow, or large-scale meteorological phenomena, including 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
low-frequency contributions may often be non-negligible, even for relatively flat sites. Furthermore studies have shown that the low-frequency contributions are highly site-specific and characterized by significant uncertainty 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 low-frequency contributions.
The importance of including vertical low-frequency 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), while other studies suggest otherwise . Kanda et al. (2004) demonstrated that, although the systematic bias decreased when including low-frequency contributions, the variance between flux estimates increases greatly. In other words, any single flux estimate becomes vulnerable to random low-frequency contributions, and thus increasingly difficult to interpret in terms of local surface fluxes. Moreover, it has been commented that horizontal low-frequency contributions, which are typically assumed negligible, may become significant during conditions of low turbulence intensity and gravitational flows (Yi et al., 2008) as well as during flow disturbance associated with complex topography (Zeri et al., 2010).
Accordingly, we can distinguish between two principal applications of the EC technique: (1) process-oriented studies in which fluxes are being linked to local biogeochemical processes for parametric insight into universal causal flux relationships and up-scaled through numerical modelling, and (2) long-term net ecosystem-exchange studies in which the flux estimates are understood to be site-specific, applying only for the unique conditions of a particular ecosystem. 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) adjusting the flux averaging time to strike an appropriate balance between adequate sampling of the turbulent flux contribution while avoiding excessive inclusion of low-frequency contributions (Sun et al., 2006); (2) ensuring horizontal homogeneous conditions within the foot print of the flux; (3) estimating vertical low-frequency contribution by performing profile measurements of fluxes on a single tower (Lee, 1998;Leuning et al., 2008) and filtering out observations reflecting excessive low-frequency influence (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 (Hojstrup, 1981(Hojstrup, , 1982Hunt 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. In the absence of a distinct spectral gap between contributions, separating flux contributions by adjusting the flux 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. In the limit of low absolute covariance (i.e., small fluxes) a relative large variance of the co-spectra estimates complicate the comparison between observed and theoretical co-spectra. While such cases could be treated as reflecting observations approaching the detection limit of the system, and discarded accordingly, they are important for exchange studies over low-flux surfaces such as sea ice, creating a demand for a new approach here.
Ensuring co-spectral similarity requires a number of site/system-specific empirical co-spectral corrections to account for high-frequency non-white noise/dampening produced by the presence of the EC system in the observed flow as well as signal dampening in closed path systems, greatly complicating the approach (Aubinet et al., 2000;De Ligne et al., 2010;Kaimal, 1968;Massman and Ibrom, 2008;Moncrieff et al., 1997;Moore, 1986;Silverman, 1968). Matching the co-spectral peak solves the issue of excessive scaling offset mentioned above, but increases the risk of subjective analysis.
Here, we present a novel method for estimating locally meaningful atmosphere-surface fluxes despite lowfrequency influences, using a single eddy covariance system and a numerical modelling scheme for ogive optimization. Accordingly we call this method ogive optimization. Ogive optimization makes no assumptions regarding optimal flux averaging time or the presence of a spectral gap and improves the flux estimates by also considering contributions in the very high/low frequency ranges. To evaluate the method, we applied it on eddy covariance observations of sensible heat, latent heat and CO 2 flux at five sites covering different ranges of fluxes, ecosystem types and topographical conditions. Results were compared with the conventional eddy covariance method both in terms of flux estimate yield and flux difference relative to flux strength.

Eddy covariance and spectral analysis
The theory of eddy covariance is well established (e.g., 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) 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 dry and r c = mol c mol −1 dry are the dry mixing ratio of humidity and CO 2 concentration scalars, respectively. The terms w θ v , w r q and w r c are the covariance between turbulent fluctuations of vertical wind and turbulent fluctuations of potential virtual temperature and dry mixing ratio of humidity and CO 2 concentration, respectively. For Eq. (1) to truly represent the vertical fluxes a number of assumptions should be met during field operation. Principal among these are stationarity of the observation ∂ξ ∂t = 0 , horizontal homogeneity ∂ξ ∂x 1 = 0 & ∂ξ ∂x 2 = 0 , mass conservation w = 0 & ∂ ς j ∂x j = 0 , negligible density flux ρ ρ 1 and a vertical constant flux layer (e.g., dF CO 2 dz = 0) (Foken and Wichura, 1996). 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 frequencydependent 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 eddy covariance system on the flow, oscillations of the tower (or ship), 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, (2) The principal use of ogives is to estimate the optimal flux averaging time as the point of convergence of cumulative cospectral energy to an asymptote (Berger et al., 2001;Foken et al., 2006). However, low-frequency influences may result in ogives which instead either converge to an extremum or diverge, depending on the direction of low-frequency fluxes. Such conditions may arise in the absence of a distinct spectral gap during significant overlap of high-frequency (turbulent) and low-frequency flux contributions. Depending on the severity of the deviation from asymptotic behaviour, an optimal averaging time can be impossible to determine. Such cases are conventionally considered to be in-stationary and no flux estimation is possible. showing examples of how low-frequency influence could 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 low-frequency flux contributions are separated by a spectral gap, allowing investigators to isolate the former simply by choosing an appropriate flux 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). Both turbulent fluxes and the low-frequency contribution, shown in Fig. 1b as a blue region of ogive divergence relative to asymptotic convergence, may be positive or negative, though the former has been illustrated as positive here.

Why eddy covariance often fails to capture local fluxes
Given the unclear existence of a spectral gap (Lee et al., 2004), however, another more general situation is the case (Fig. 1c) of overlapping contributions from low-frequency motions, turbulence and site and instrument-specific nonwhite noise/dampening. One way to strike a balance between adequate inclusion of the turbulent contribution and exclusion of excessive low-frequency influence is by adjusting T 1 . Typically a fixed averaging time is set for an entire experiment (here 30 min is shown) and the flux errors are assumed, or tested (e.g., Novick et al., 2014), to be negligible. In the high-frequency end of the spectrum, instrument response limitations may prevent observation of the smallest scales of turbulence contribution (Here 10 Hz is shown). Furthermore, instrument-specific non-white noise and/or 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;Massman and Ibrom, 2008;Kaimal, 1968;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 influence of dampening and non-white noise on the ogive distribution occurs in the reverse (Fig. 1d) relative to co-spectral space (Fig. 1c).
Observations reflecting excessive low-frequency influence, relative to the turbulent contribution, (Fig. 1e) are typically discarded. This is because strong relative lowfrequency influence results in non-negligible flux contribution to the overall estimate and further obstructs any efforts to separate contributions by adjusting the flux 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 low-frequency influence 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. This prohibits 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 fulfil 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 require-ments that T 1 should be long enough to reduce random error (Berger et al., 2001;Lenschow and Stankov, 1986) and short enough to avoid low-frequency influence associated with non-stationarity (Foken and Wichura, 1996;Vickers and Mahrt, 1997). However, as noted, adjustment of T 1 will not generally allow for separation of turbulent and low-frequency 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 low-frequency contributions.
The following is an iterative scheme for developing averaging intervals based on basic data quality requirements. Data collected during a field-experiment is considered 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 (i) within the range T min + T 1 2 to T max − T 1 2 in intervals of T 1 , and data set lengths T w (i) within the range 10 to 60 min in 5 min intervals. Here T 1 is set according to desired temporal resolution of flux estimates. For this study we use sitespecific settings of T 1 = 5 min, T 1 = 15 min and the conventional T 1 = 30 min to strike a balance between desired temporal resolution and computational cost of running the ogive optimization method. The minimum data set 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 data set length is chosen to be 60 min to ensure approximate stationarity of the local turbulencedriven 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 data set around T w (j ), which passes the assessment criteria. These include absence of instrument diagnostics errors, absence of long data gaps, favourable mean wind direction and reasonably narrow range of wind directions (≤ 60 • ). Minor spiking (≤ 1 %) is corrected based on the median of surrounding data points, and the data set is discarded otherwise (for spiking > 1 %) according to Vickers and Mahrt (1997). Other quality assessment includes the requirement that momentum fluxes should be negative. 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 also be affected by low-frequency contributions. As direction, 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 horizontal along-wind component, and simply verify 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 instrument-specific non-white noise, dampening and low-frequency influences. Typically around four es-Atmos. Chem. Phys., 15, 2081-2103, 2015 www.atmos-chem-phys.net/15/2081/2015/ timates of T w (j ) are evaluated by the iterative bisectional algorithm for each T m (i) before the optimal T w (j ) is determined. Finally, signals with very rapid evolutions such as transient signals in dynamic systems like eddy covariance observations may undergo abrupt changes associated with observational interference 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 leading to semi-permanent changes (Lee et al., 2004;Mahrt, 1991;Vickers and Mahrt, 1997). In this study we perform the Haar analysis for the data set T w (j ), given by the bisectional algorithm, and observations are sorted in three categories (good, soft flag, hard flag) according to the presence and severity of signal discontinuities (Vickers and Mahrt, 1997). Conventionally these events are thought to preclude flux estimation on the basis of stationarity violations. Accordingly we discard such data sets (hard flags) when applying the conventional eddy covariance method. 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. Therefore no flux estimates derived using the ogive optimization are discarded, unless visually inspected.

Mass ogive calculation
As noted, adjusting the flux averaging time will not generally allow for separation of turbulent and low-frequency flux contributions. Subtracting a running mean from observed signals, as opposed to the conventional linear detrending, allows for enhanced filtering of low-frequency contributions alone (Sakai et al., 2001;Mcmillen, 1988). Consequently some combination of data-set length (averaging time) and running mean window size might allow for filtering out of low-frequency contributions while retaining turbulent contributions. Note that both adjusting the flux averaging time and subtracting a running mean from the observed signal may, in many cases, provide sufficient separation of turbulent fluxes and low-frequency contributions. Here we apply both to arrive at a more generally applicable approach. We visualize this concept by calculating co-spectra, and corresponding ogives, for a very large number of data permutations and derive 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 10 min and the maximum time available (60 min in this example) and 200 linear increments for the running mean window in the range of 1 min to half the length of the data set in question (30 min in this example). 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 asymptote. A classic ogive shape. The flux estimate for a regular 30 min linear detrended data set (red) appears representative of the overall ogive pattern as well. The presence of haze on the graph below −60 Wm −2 suggests that a small part of the permutations states are affected by low-frequency motions. The presence of these large-scale motions is in part supported by the Haar analysis, which has soft-flagged the temperature signal (Fig. 2b). Additionally, inspection of a running covariance with a 5 min window (Fig. 2c) indicates the onset of increased flux variability in the range 20-60 min.

The model and the optimization method
Unfortunately not all ogive density maps indicate as well defined fluxes as shown in Fig. 2. In such cases, answering the overall question of most likely flux requires the fitting, or optimization, of an ogive model to the ogive density map. With the introduction of an optimization aspect the advantage of performing the analysis for ogives, as opposed to co-spectra, becomes clear. In the limit of low absolute covariance (i.e., small fluxes), co-spectra typically become increasingly characterized by both positive and negative frequency-wise flux contributions. The co-spectral model, however, can only account for fluxes in one direction. Observed and modelled ogives, in contrast, are able to describe and account for this bidirectionality.
The basic premise in our model solution is that a region exists in the mid-to-high frequency range of the ogive representation, which is least impacted by instrument-specific non-white noise, dampening and low-frequency influences. This was illustrated in Fig. 1b, d, f. While such a region is clearly evident in Fig. 2a for f > 5 × 10 −2 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 the lower observational bound (f = 2.5 × 10 −4 , or 60 min). While the cospectral 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) f 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 4/3 power law in the inertial subrange. Subsequently an ogive distribution is calculated using Eq. (2). We set A  the low-frequency end-point f 1 (equivalent to the averaging time T 1 ) of the ogive distribution to a variable parameter F 0 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/under estimation of covariance due to inclusion of low-frequency contributions or the use of inadequate averaging times. Similarly, in the high-frequency range the problem may arise in the form of under-estimation of covariance due 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 flux, 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 F cor = F + F LF1 − F HF1 , where F is the uncorrected flux. Note that F HF1 is subtracted as it is of opposite sign relative to F . 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 a 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 as F cor = F + F LF2 − F HF2 .

Combined the corrections amount to
where O 2 is adopted as shorthand mathematical notation for ogive optimization. In practice the corrections are accounted for by combining F and all low-frequency contributions into one parameter: F 0 = F + F LF1 + F LF2 , which controls the shape of the ogive model, and by adding a high-frequency vertical offset F HF = − (F HF1 + F HF2 ). In total we are thus left with four tunable parameters: F 0 , F HF , µ and f x , for which the final model-estimated flux is 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 parameter combination, for which the model ogive follows optimally the strongest densities in the density map (Fig. 2a). 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 include optimization of a random parameter guess within reasonable parameter bounds using a fast local optimization algorithm and a slower, but robust, Darwinian evolution-style global optimization algorithm called Differential Evolution (Storn and Price, 1997). The optimization is performed for a number of frequency intervals and the final solution is chosen by subjective visual inspection. An in-depth explanation of these steps can be found in Appendix A. One intriguing consequence of including a modelling 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 ecosystem, topography and flux strengths (Q SENS , Q LAT and F CO 2 ) are investigated (Fig. 4).

High-flux environment
The Abisko field-site (Fig. 4a) is located in Stordalen (68 • 21.248 N, 19 • 3.02 E), a mixed mire 10 km east of Abisko in subarctic Sweden, the site of a number of past and ongoing studies on mire carbon fluxes (e.g., Christensen et al., 2012;Jackowicz-Korczynski et al., 2010). The measuring mast is situated on the coastal edge of a minerotrophic fen dominated by sedges and a lake environment. Wind patterns consistently alternate between the upwind fen-environment signal towards the west and the upwind lake-environment signal towards the east. Hence we treat the observations separately as fen and lake sites, respectively. Continuous eddy covariance observations were conducted from 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.9 m height and an LI-7500 open path gas analyzer (LI-COR ® , Lincoln, NE, USA) on a boom extending towards the southwest at 2.5 m height, with a 0.43 m horizontal off-set along the boom and a slight tilt of the instrument relative to the vertical plane to allow water dripping. 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 T ABI 1 = 30 min.

Low-flux environment
Young Sound (Fig. 4c) 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 polynya at the mouth of the fjord (Rysgaard et al., 2003). Continuous eddy covariance observations were conducted at three sites within the fjord system in the period 20 March to 27 April 2012. Two separate fieldstations, one static and one mobile, were used at three different locations (ICEI, POLYI and DNB Raw data were logged at 20 Hz. The sonic anemometer and gas analyzer were displaced horizontally by 0.4 m in orthogonal alignment to the prevailing 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 shore-adjacent 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.7 m relative to the snow surface 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. Raw data were logged at 10 Hz. Sea-ice and snow-cover thickness was approximately 110 and 75 cm for the ICEI and DNB sites, respectively, and approximately 25 and 20 cm for the POLYI site. Average air temperature increased from −35 to −15 • C during 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 low-frequency motions due to their locations in a fjord surrounded by mountains. Flux estimates were evaluated at the three sites for the following intervals: T POLYI 1 = T ICEI 1 = 5 min and T DNB 1 = 15 min. The higher resolutions of flux estimates, relative to the Abisko and RIMI sites, were chosen for the purpose of another study concerning CO 2 fluxes on sea ice.

Instrument corrections and post-processing
During post-processing, 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 the following: an empirical angle of attack correction (Nakai and Shimoyama, 2012), and humidity and crosswind corrections (Liu et al., 2001;Schotanus et al., 1983). We convert all observations to mixing ratios (Burba et al., 2012) using the Webb-Pearman-Leuning 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, POLYI and ICEI) 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. Coordinate rotation, linear de-trending and iterative de-spiking of raw data is performed according to Vickers and Mahrt (1997). 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. Here, we locate the optimal cross-correlation automatically based on the incident horizontal wind flow and the specific geometry of each covariance system.

Examples of ogive optimization performance
In the following, we describe several typical cases observed and the associated performance of the ogive optimization method.
1. Near-absence of low-frequency influence is observed leading to a strong similarity between the ogive density pattern, the 30 min linear detrended ogive and the modelled ogive. This is illustrated in Fig. 5a for a case of sensible heat flux at the Abisko site. Disregarding the high-frequency component associated with extrapolation of model results, seen here to contribute F HF ≈ −5 Wm −2 to the overall modelled sensible heat flux (Eq. 4), the standard 30 min linear detrending approach will suffice to provide the turbulent flux estimate. The observational period in question was characterized by neutral atmospheric stability, high wind speed U = 8.25 ms −1 and winds originating from the fen area (Fig. 5a), along with fairly constant temperature 8.2 ± 0.4 • (Fig. 5b) and a slight gradual increase in uptake (Fig. 5c). 2. Cases where non-negligible low-frequency influence on the flux estimate is observed for CO 2 flux (Fig. 6). The low-frequency contribution is seen to be positive just like the turbulent flux (Fig. 6a). The ogive optimization method is seen to separate the turbulent and the low-frequency contributions completely, yielding only the locally meaningful turbulent flux. The observational period in question was characterized by a slightly unstable atmospheric stability zL −1 = −0.19, moderate wind speed U = 3.65 ms −1 and wind originating from the lake area (Fig. 6a), along with a slightly increasing atmospheric temperature and a varying CO 2 concentration (Fig. 6b). A consequent marked increase in flux covariance around 35-45 min is evident in Fig. 6c. 3. An example of ambivalence caused by bimodality in the ogive density pattern is illustrated in Fig. 7, for a case of sensible heat flux at the Abisko site. Such cases indicate that fluxes are changing within the sampling period. The ogive optimization method is seen to capture the turbulent flux contribution with the strongest data density. Had both modes been of equal ogive density, the choice of mode during subjective evaluation would be based on the quality of the model ogive optimization, and the length of the time-series responsible for the modes. If both ogive models were equally good, the choice would fall on the mode produced by the ogives which consist of shorter time-series as they represent a more instantaneous flux estimate relative to the mode produced by longer time-series. The sampling period in question was characterized by a slightly stable atmo-sphere zL −1 = 0.12, low wind speed U = 2 ms −1 and wind originating from the fen area (Fig. 7a), along with a steady decline in atmospheric temperature (Fig. 7b) and a strong variation in flux covariance (Fig. 7c).

The inadequacy of applying a fixed averaging interval
for flux estimation becomes apparent in Fig. 8 seen to be too long and also to increasingly reflect low-frequency interference (Fig. 8a). This is consistent with an abrupt increase in atmospheric temperature (Fig. 8b) and decrease in covariance (Fig. 8c) around 30-40 min. The ogive optimization method identifies the appropriate flux estimate (Fig. 8a), whereas the standard 30 min linear detrending method fails on account of in-stationarity. In addition, the case is a perfect example of how co-spectral evaluation of frequency-wise contributions can be misleading (Fig. 8a, inner plot). The observational period in question was characterized by a stable atmosphere zL −1 = 0.29, low wind speed U = 2.6 ms −1 and wind originating from the fen area (Fig. 8a).

The inadequacy of applying a fixed averaging interval
for flux estimation becomes apparent in Fig. 8, for a case of sensible heat flux at the Daneborg site. Here, the ogive density pattern is seen to reflect a gradual evolution in the ogive flux pattern with increasing averaging time. The standard 30 min averaging time is seen to be too long and also to increasingly reflect low-frequency interference (Fig. 8a). This is consistent with an abrupt increase in atmospheric temperature (Fig. 8b) and decrease in covariance (Fig. 8c) around 30-40 min. The ogive optimization method identifies the appropriate flux estimate (Fig. 8a), whereas the standard 30 min linear detrending method fails on account of in-stationarity. In addition, the case is a perfect ex- ample of how co-spectral evaluation of frequency-wise contributions can be misleading (Fig. 8a, inner plot). The observational period in question was characterized by a stable atmosphere zL −1 = 0.29, low wind speed U = 2.6 ms −1 and wind originating from the fen area (Fig. 8a).

Atmos
6. Signals may be degraded for a number of reasons such as instrument failure, electronic interference etc. Such a case is illustrated in Fig. 9, for a case of CO 2 flux at the Abisko site. Here a brief drop in atmospheric CO 2 concentration, hard flagged by the Haar analysis (Fig. 9b), gives rise to an intermittent 3-fold increase in flux covariance (Fig. 9c), ultimately resulting in the lowfrequency influences illustrated in Fig. 9a. Nonetheless, the ogive optimization method is seen to identify the actual prevalent flux during this period. The observational period in question was characterized by nearneutral atmospheric stability zL −1 = −0.07, moderate wind speeds U = 5 ms −1 and wind originating from the lake area (Fig. 9a), along with a near-constant air temperature T air = 6.9 • (Fig. 9b).
7. During conditions of strong high-frequency dampening caused by the use of a closed path instrument, the ogive optimization method automatically shifts the high-frequency bound on optimization towards lower frequencies to avoid influence of the dampened frequencies during optimization. This is illustrated in Fig. 10 for a case of latent heat flux at the ICEI site. Here the upper optimization bound is shifted back to 1 Hz thus allowing for an accurate description of the high-frequency fluxes as well (Fig. 10a, inner plot). The observational period in question was characterized by a slightly unstable atmosphere zL −1 = −0.08 and moderate wind speeds U = 5 ms −1 (Fig. 10a), along with gradual increases in both atmospheric H 2 O content and atmospheric temper-ature (Fig. 10b) and a gradual increase in flux covariance (Fig. 10c).

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 lowfrequency contributions, the inadequacy of the fixed averaging interval and the extrapolation of modelled ogives into unobservable high/low frequencies. δ is evaluated within i = 10 intervals of absolute flux estimates F 30 min i as the standard deviation of difference in flux estimate relative to the mean absolute ogive optimization flux estimate within respective intervals: where square brackets [ ] i signify flux estimates native to interval i of the equivalent absolute flux estimates by the standard eddy covariance method F 30 min i . Estimates of relative flux difference δ are shown logarithmically in Fig. 11 for all three scalar flux types, at all five 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 three 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).
As hypothesized in Sect. 2.2, the average relative flux difference is seen to be very high for small absolute flux estimates, peaking at δ SENS = 80 %, δ LAT = Figure 11. Relative difference in percent (see Eq. 5) is shown logarithmically as a function of absolute flux estimate for all investigated sites. Also shown are the median (red line), standard deviation (light gray area) and 25-75 % percentile (dark gray area) of the relative differences. In the bottom of the figure, histograms of absolute ogive optimization flux estimate ranges are shown for each site. Numbers indicated to the left of the histograms are the respective maximum values. 23 % and δ CO 2 = 98 % for the lowest absolute flux estimates. The variation in δ is quite high for low absolute flux estimates, with the 13.6 and 86.4 percentiles of δ reaching as much as δ SENS = 40 − 208 %, δ LAT = 10 − 98 % and δ CO 2 = 52 − 538 %. For larger absolute flux estimates (|Q SENS | > 40 Wm −2 , |Q LAT | > 20 Wm −2 and F CO 2 > 100 mmol m −2 d −1 ) the relative difference is seen for all three flux types to drop and level off to a near-stable range of δ = 5-20 %. These absolute flux thresholds thus mark clear shifts between non-negligible low-frequency contributions on one side and plausibly negligible low-frequency contributions on the other.
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 local-scale flux resolvability by the standard 30 min linear detrending ap-proach, or as an argument for the application of enhanced flux estimation techniques such as the presented method. For the presented observations the consequences are illustrated by the histograms of the different sites (Fig. 11). Although the location of the flux threshold is a bit unclear for latent heat flux, estimation of locally meaningful fluxes at the three sea-ice sites Daneborg, POLYI and ICEI is essentially impossible without accounting for low-frequency 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 CO 2 = 400 mmol m −2 d −1 and that most flux estimates from RIMI are from the morning or late evening/night, which explains the range of relatively small fluxes.
The relative flux difference was furthermore investigated in terms of atmospheric stability (Fig. 12). Though variation in δ is significant for all flux types, average δ appears to be lowest between slightly unstable zL −1 ≈ −0.2 and neutral conditions (δ ≈ 10-20 %). In contrast δ is significantly larger for zL −1 > 0.2 (20 < δ < 1000 %). This is consistent with the current consensus on turbulence spectra (Kaimal et al., 1972;Olesen et al., 1984): for strongly unstable conditions zL −1 < −0.2 all spectra have increased low frequency components (Hojstrup, 1982), which would have been filtered out using the ogive optimization method. For strongly stable conditions zL −1 > 0.2 the turbulence spectral intensity is often small relative to low frequency variation associated with meso-scale variability (Larsen et al., 1980;Vickers and Mahrt, 2003). Exactly for neutral and slightly unstable conditions boundary layer turbulence structure is at its simplest being dominated by shear produced turbulence that is best described by the standard spectral expressions, being the background for the standard eddy-correlation flux determination methods (Kaimal et al., 1972;Olesen et al., 1984).
For many flux estimates the vertical wind speed signal or the scalar signal are non-stationary to the point of prohibiting a flux estimation using traditional methodology. Hence the ogive optimization method may also provide a greater number of flux estimates. This is shown in Table 1 to generally be true for the Abisko and Daneborg sites, both of which characterized by degraded signal quality at times. Sites RIMI and POLYI are inconclusive in this respect and the conventional method appears superior in the case of ICEI. The latter may be related to the very low fluxes observed for this site (Fig. 11) suggesting the presence of a detection limit for the ogive optimization method when using the particular instrument setup at ICEI (LI-7200 enclosed gas analyzer) within the respective ranges |Q SENS | < 25 Wm −2 , |Q LAT | < 3 Wm −2 and F CO 2 < 10 mmol m −2 d −1 . No similar characteristics indicate the presence of a detection limit of the ogive optimization method at any of the other sites (open path gas analyzers), suggesting superiority of open path instruments in very low flux environments when using the ogive optimization method. 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, indicating a puzzling caveat to current theory. The occurrences challenge the notion that fluxes should be of same sign regardless of incident eddy scales. One explanation might be that vertical low-frequency contributions represent only one part of a net low-frequency contribution and hence is balanced by a horizontal component. Indeed the horizontal lowfrequency 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 lowfrequency contributions, and their influence on turbulent flux estimates, is necessary.

Conclusions
The presented ogive optimization method has been shown to successfully separate local from non-local flux contributions. In addition, it enhances flux estimation by both investigation of a large range of averaging times and running mean detrending, and extrapolation of optimized ogive model results. The method makes no assumptions concerning ap-propriate averaging time or the presence of a spectral gap, does not require the application of transfer functions and allows for very high temporal resolution of flux evolution. For high flux rates (|Q SENS | > 40 Wm −2 , |Q LAT | > 20 Wm −2 and F CO 2 > 100 mmol m −2 d −1 ) we found that the average relative difference between fluxes estimated by ogive optimization and the conventional method was low (5-20 %) suggesting negligible low-frequency influence and that both methods capture the turbulent fluxes equally well. For flux rates below these thresholds, however, the average relative difference between flux estimates was found to be very high (23-98 %) suggesting non-negligible low-frequency influence and that the conventional method fails in separating low-frequency influences from the turbulent fluxes. The average relative flux difference was found to be lowest (10-20 %) for slightly unstable and neutral atmospheric stabilities zL −1 = −0.2 to zL −1 = 0. In contrast δ was significantly larger (20-1000 %) for zL −1 > 0.2. This is consistent with current consensus on turbulence spectra. Furthermore, the ogive optimization model has been shown to allow for flux estimation despite signal disruption. Fewer flux estimates could be derived relative to the conventional method for an LI-7200 enclosed gas analyzer in very low flux conditions, suggesting the possible presence of a detection limit in the |Q SENS | < 25 Wm −2 , |Q LAT | < 3 Wm −2 and F CO 2 < 10 mmol m −2 d −1 ranges 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 over sea ice. Overall, the notion of a dynamic and generally nonnegligible overlap of low-frequency and turbulent flux contributions is confirmed.
Finally, 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 low-frequency contributions during certain conditions (Yi et al., 2008;Zeri et al., 2010) we hypothesize a more intricate balancing interplay between vertical and horizontal low-frequency flux contributions which, if confirmed, suggests the need for more sophisticated eddy covariance system arrays if low-frequency contributions are to be accurately included (i.e., for site-specific studies). If exclusion of low-frequency contributions is desired (i.e., for universalprocess-oriented studies), the presented method should be unaffected by these questions. As described in Sect. 2.4.2 our goal is to tune the four final model parameters F 0 , 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:

Atmos
1. A random guess of parameters is made within a set of reasonable bounds. The speed and accuracy of any optimization method involving pre-set bounds depend greatly on the reasonable choice of these bounds. Here we set 0.05 < µ < 1 and −2 < log (f x ) < 1. The bounds for F 0 and F HF are a bit more complicated. If > 80 % of the summed density map is located on, say, the positive side (suggesting F 0 > 0), bounds on F 0 and F HF are set as 0 < F 0 < |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 0 < 0 and 0 < F HF < |R − |), and optimization is performed twice, using the two different sets of bounds, if neither side contains > 80 % of the summed density map.
2. Think of optimizing a model ogive to an ogive density map as choosing a path between two points in the Pyrenees for which you travel at the highest possible average altitude, all the while being constrained to a certain type of path (the ogive form and the associated parameters). In more technical terms optimizing the four parameters of the model ogive may be thought of as locating the point in a parameter-wise four-dimensional probability space, for which the net ogive density reached along the path is the highest. In this context we seek a global, as opposed to local, solution within the probability space formed by the four parameters. Based on the initial random guess, a local solution is determined using the MATLAB function fminsearchbnd (available through the Mathworks ® file exchange) 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 computationally expensive.
3. Based on the local optimization of parameters 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). Inspired by Darwinian evolutionary theory it optimizes a problem by iteratively improving a population of NP candidate solutions (agents) based on random candidate mutation (motion) 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 observational 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. If enough agents (NP) are initiated and allowed to evolve throughout the probability space sufficiently long (iterations) the DE algorithm is certain to locate a global solution (optimal ogive parameters).
4. Often optimizing a smaller subset of the problem is an advantage, particularly during low-frequency interference which persists despite data perturbation in the mass ogive phase. One such case is shown in Fig. 13. Optimizing in subsets is achieved by subdividing the problem into 18 frequency interval weights in the range 0 to 1, signifying 0 to 100 % influence of a given part of the density map on the optimization output (Fig. 13a, black lines). Corresponding solutions for the 18 frequency interval weights are shown in Fig. 13b (green lines). All solutions based on frequency intervals with lower bound before or after the ogive density peak f ≈ 7 × 10 −4 are seen to underestimate the actual undisturbed turbulent flux. Accordingly an appropriate solution (blue line) may be estimated within the subset solution for which the frequency interval features the ogive density peak as its lower bound. Essentially, the optimization 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. 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. Note the gradual decrease in optimization weighting of high-frequency ogive density (Fig. 13a), which has been added to limit any influence of high-frequency instrument-specific non-white noise and dampening during the optimization. The highfrequency limit of the fitting interval is furthermore al- lowed 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. The latter is illustrated in the results and discussion section of this study.

Code availability
The executable code of our procedure, ogive optimization, will be made available and can be acquired by e-mailing the corresponding author (jasi@envs.au.dk or lls@bios.au.dk). The program is coded in MATLAB and is optimized for use with the parallel computing toolbox.