Impact of coordinate rotation on eddy covariance fluxes at complex sites

The choice of coordinate system to calculate eddy covariance fluxes becomes particularly relevant at complex measurement sites. The traditional way is to perform double rotation (DR) of the coordinate system i


Introduction
The eddy covariance (EC) method is the most direct and defensible way to measure vertical turbulent fluxes of momentum, energy and gases between the atmosphere and biosphere.In order to avoid methodological errors in flux calculations, the coordinate rotation and flux averaging should be performed in a consistent way (Finnigan et al., 2003(Finnigan et al., , 2004)).Coordinate rotation performed according to velocity statistics over a defined time period T (typically 30 min) is a common practice to define the coordinate system for flux calculation.The rotation could be done as three consecutive rotations (triple-rotation, TR) where the first two rotations are done such that the x-axis is oriented along the mean wind and the average vertical wind speed becomes zero.The third rotation was proposed by McMillen (1988)

Table 1
Description of sites and main topographic features.

Site Type
Location Refernece 1 Kumpula, Station for Measuring Ecosystem -Atmosphere Relationships, SMEAR III About four kilometers northeast from downtown of Helsinki, Finland.The heterogeneous surrounding area consists of buildings, patchy forest and low vegetation, parking lots and roads.
The surface can be classified to three sectors.The building sector in the northern direction (320 -40°N) is dominated by university campus and the Finnish Meteorological Institute buildings (mean height 20 m) in the vicinity of the tower.One of the main roads leading to Helsinki city centre passes the road sector (40 -180N°) with the closest distance between the road and the tower being 150 m.The area between is covered with deciduous forest.In the vegetation sector (180 -320N°), most of the surface is covered with vegetation in the University Botanical Garden and in the City Allotment Garden.
60°12.17′N, 24°57.671′E,26 m a.s.l.Its topography is characterized by an alpine slope with a mean slope of about 11°falling in a N-S direction.About 60m upslope to the north of the main tower, a pasture disturbs the fetch in the main nighttime wind direction.In addition to the S-SE slope (toward Isarco Valley) there is another slope west (270°N) toward Sarentino valley.
The meteorological conditions are dominated by a very persistent local slope wind system with upslope S-SW winds during the day and downslope NNW winds during the night.

State Forest
Located in south-central Indiana, USA.
Trees in the tower footprint 60-80 years old.The forest is a secondary successional broadleaf forest within the maple-beech to oak hickory transition zone of the eastern deciduous forest.
The area surrounding the tower is dominated by a ridge-ravine topography with relative relief less than 60 m over length scales of several hundred meters.The tower is located on a ridge which extends in NW-SE direction.The terrain slopes downward to the SW and NE, with two gullies cut into these slopes.Each gully is oriented approximately radially to the tower and has an elevation drop of approximately 80 m over a distance of 2-3 km (slopes in the range of 2% to 4%).Within the immediate vicinity of the tower (i.e.250 meter radius), local slopes range between 3 -6% along the gullies, and <2% elsewhere.The principal wind direction is from the W-SW.
39°19′22.8"N,86°24′48.0"W,275 m a.s.l.Froelich et al. (2005) 5 Hyytiälä, SMEAR II Located in Hyytiälä, Southern Finland.Scots pine with mixture of Norway spruce and birch, 55-years old.The tower is located on the southern part of a very low spur ridge which extends for about 1 km along the northsouth direction.It degrades to W-SW direction, where at a distance of about 700 m is an oblong lake, situated at 150 m above sea level.The slope becomes more gentle in E-SE direction, where the elevation difference is up to 10-15 m over a distance of 1.5 km.On the northern side of the measuring tower, the terrain is almost flat reaching the ridge top (193 m a.s.l.) over a distance of 500 m.
During the nights decoupling of the flow over and below the forest has been observed leading to night-time drainage of the carbon dioxide (Alekseychik et al., 2013).
179 m a.s.l.Hari and Kulmala (2005) to minimize the vertical momentum flux in cross-wind direction (e.g., Kaimal and Finnigan, 1994).Finnigan et al. (2004) however have shown that the third rotation utilising the stress tensor cannot be reliably defined for the complex real flows.Therefore, only the DR of coordinate frame over the averaging period is recommended in practice.
An alternative approach is to define the plane to which the velocity vectors position the closest in the least square sense over a longer period of time (Wilczak et al., 2001).Since the method is applied by aggregating the velocity vectors over a long period of time, it can be interpreted as finding the ensemble-averaged plane and is therefore frequently referred to as the planar-fit (PF) method.Wilczak et al. (2001) propose a linear fitting procedure to define the alignment of the plane, interpreting the intercept of the fit as the anemometer offset.Mammarella et al. (2007) have performed the fitting procedure according to the wind direction, obtaining a plane for each narrow sector around a measurement tower.This variant of the coordinate rotation method is usually referred as sector-wise PF (Sabbatini et al., 2018).Finnigan (2004) has analysed in detail the choices of coordinate systems and concluded that the PF method is likely to result in more stable flux results over steep sloping topographies than the DR method.However, application of the PF method to each particular site requires decision-making about the sector sizes, thresholds to exclude abnormal conditions and possible other factors affecting the flow conditions at the site.
The flow can be tilted with respect the anemometer due to several reasons.The most obvious reason is the tilted positioning of the anemometer with respect the local terrain or the flow field.The instrumental tower can also deflect the flow to the extent that can lead to tilted coordinates (e.g., Dabbert, 1968, Le andHu, 2002).Further reasons for the measured wind field tilt can arise from the three-dimensional nature of air motion over horizontally inhomogeneous canopy or terrain (e.g., ejections and sweeps, Bohrer et al., 2009).Such tilting of the measured wind in anemometer's co-ordinate system leads to a non-zero measured vertical wind speed that is attributable to the mean wind but is detected as the vertical component due to anemometer's tilt with respect to the flow field.
It is possible that a mean vertical wind component that is not related to horizontal heterogeneity can exist.This component is a residual mean vertical velocity in the long-term average coordinates.In addition to the previously mentioned anemometer offset, there can be several meteorological reasons for the residual mean vertical velocity in the long-term coordinates (e.g.Lee, 1998;Lee et al., 2004).A nonzero vertical wind may exist due to local thermal circulation and convective conditions.Vertical velocities as large as 10 cm s −1 are possible at 40 m above the ground within the descending air column according to the simulation study of Moeng (1984).Mahrt and Ek (1993) found that a persistent ascending motion existed in the daytime over the heterogeneous landscape dominated by forest and cropland, with a mean value of 2 cm s −1 at a height of 125 m above the ground.In case of horizontal flow convergence/divergence, the (30-min) mean velocity vector will no longer be parallel to the terrain surface.The divergence along the slope due to the gravitational acceleration is compensated by a downward air motion and has been shown to yield a mean vertical velocity of a few cm s −1 at a measurement height of tens of meters (e.g., Lee, 1998).
Different interpretations and use of vertical velocities in defining the coordinate systems for flux calculations can result in systematic differences in fluxes.Lee et al. (2004) demonstrated that the tilt error of a few degrees can result in up to 5% flux bias.Such errors can depend on the time of day and thus have a potential to end up with significant systematic bias in net ecosystem exchange that is typically the sum of daytime uptake and night-time emission of CO 2 .Therefore, the choice of the appropriate coordinate system is of high importance for complex sites.
Attempts have been made to improve the surface-atmosphere exchange estimates by using the residual vertical velocities with respect to the long-term average to derive the missing flux components.Lee (1998), Baldocchi et al. (2000) and Paw U et al. (2000) combine the residual with the continuity equation to estimate the contribution of vertical advection to the surface layer mass balance.Mammarella et al. (2007) analysed the role of vertical advection in night-time CO 2 exchange at the forested site in Southern Finland and found that accounting for the term improved the CO 2 exchange to biologically feasible levels under low turbulence conditions.Rannik et al. (2009) observed the same for night-time ozone deposition at the same site.Several other studies have considered vertical as well as horizontal advection terms in carbon dioxide budget.Aubinet et al. (2003) found that the advective fluxes were frequently opposite in direction and of similar order of magnitude.Due to large variability and uncertainty the advective fluxes did not enable to achieve feasible exchange estimates.More extensive measurements found that the advective terms were site specific and the size of the advective fluxes was of the same order of magnitude as turbulent fluxes (Feigenwinter et al., 2008;Leuning et al., 2008).Finnigan et al. (2003) instead attributed the variation in mean vertical wind speed to lowfrequency component of air motions and proposed to calculate the contribution of the respective cospectral densities to the total flux.The issue is still unresolved and, because the residual vertical wind speed is highly variable, it is a very challenging task to attribute its value to the flow mechanism and evaluate the relevance of the vertical motion to the surface-atmosphere exchange at each particular observation period.
The current study investigates the alternative coordinate systems and residual vertical velocities defined by the flow fields by using the measurement data from a variety of sites ranging from complex urban and forest sites to nearly ideal forest and peatland sites.In particular, the aims were i) to evaluate the goodness of fit of PF coordinate systems and the influences of atmospheric stability and diurnal variation on the residual mean vertical winds, ii) to provide guidelines for choosing the method of defining the coordinate systems for flux calculations, iii) to assess the variability of night-time vertical advection estimated from the vertical wind relative to long-term coordinate plane, and iv) to analyse the impact of the choice of coordinate system on cumulative flux estimates.

Double rotation of fluxes
The commonly used coordinate system is the stream-lined coordinate system as determined by flow statistics over the averaging period.Relevant coordinate rotation is performed preferably as two consecutive rotations (Finnigan et al., 2004).The tilt rotation angle is defined by the horizontal and vertical wind components in the anemometer's coordinate system by Kaimal and Finnigan (1994 where 2 denotes the horizontal wind speed at the wind direction θ and overbar denotes time averaging over period T.

Ensemble-averaged coordinate systems
The ensemble-averaged methods attempt to define the ensemble averaged coordinate systems that correspond in general to flow streamlines.Wilczak et al. (2001) The PF method approximates the streamlined coordinate system by a linear relationship

Planar-fit method of
where the coefficients b 0, 1, 2 can be obtained by linear regression over an ensemble of observations.The method is usually applied sector-wise i.e. the regression coefficients are determined for a set of observations within wind direction interval Δθ.Further in this paper we refer to sector-wise PF method as SPF.
The PF method defines a different set of coordinate rotations and the fluxes obtained in the PF coordinate system can be obtained by the procedures described in Wilczak et al. (2001).The coefficient b 0 can be interpreted as the instrumental offset, and in such case the comparable tilt angle to a rotation ϕ DR can be obtained from (3) 2015) Ross and Grant (2015) propose a continuous fit method that approximates the tilt angle ϕ DR over the ensemble of observations as the Fourier series

Method by Ross and Grant (
. The coefficients α n and β n are obtained by fitting the periodic function ϕ CF of wind direction θ (in radians) to the tilt rotation angles obtained by the DR method.
The authors argue that the potential disadvantage of the method is missing the possibility of quantifying an offset term in vertical wind speed.As we see later, neglecting the offset term can lead to systematic difference in calculated fluxes.An alternative continuous planar fit (CPF) is constructed here by defining the relationship in terms of wind components as Note that this representation is similar to (2) when writing . This approach approximates in principle the PF method as the continuous function of wind direction and accounts for the wind-direction dependent intercept term that is not correlated with horizontal wind.The coordinate rotation tilt angle defined by this approach is given by (5)

Sites and measurements
Measurements from several EC sites reported in Table 1 were used.The sites were categorized in the current study into complex urban (Kumpula), complex forests (Coweeta, Renon, Morgan-Monroe), medium complexity forests (Hyytiälä, Svartberget), low complexity forests (Sodankylä, Sorø) and ideal sites (Norunda, Abisko-Stordalen) based on the topographic features and flow patterns earlier observed at the sites.Table 2 summarises the measurements setups at all sites.

Data processing
Turbulent fluxes and other statistics reported in the study were calculated over 30 min averaging period by block averaging approach.The analysis of coordinate systems and rotations was based on the 30 min average turbulent wind speed records available in anemometer's coordinate system.For some sites and variables (summarised in Table 5 for the systematic differences in fluxes when rotations to different coordinate systems were performed), the covariances between all wind speed components and scalar concentrations were available.Those covariances, including after rotation to a final coordinate system, were not corrected for underestimation due to low and high pass filtering.
Prior to flux calculation, raw data despiking was performed.The fluxes used in the analysis of the effect of vertical advection were corrected for frequency response underestimation at low and high frequencies according to commonly accepted procedures (Mammarella et al, 2009;Sabbatini et al., 2018).For flux calculations the software routinely used for flux processing were applied (Table 2).
Fittings of the long-term average coordinate systems were done using data with the stability parameter confined by −1.5< ζ <1.5, the friction velocity u * > 0.05 m s −1 and 2.5% of observations with lowest and highest wind speed conditions (in total 5%) were omitted.
From measured profiles (Table 2) the storage change term was estimated by calculating the central three-point difference in time of the storage of CO 2 below the flux measurement level.The concentration was assumed to change linearly with height in between the observation points and extrapolated to the surface by using the slope of the lowest two levels.The vertical advection of CO 2 was estimated according to the method proposed by Lee (1998) by using the vertical wind speeds with respect to long-term average flow streamlines.

Results
Throughout the results section, the Sorø, Hyytiälä and Renon sites were selected to represent the sites of low, medium and high complexity and are discussed in more detail.The other sites are discussed when relevant or in conjunction of impact of different coordinate rotations on turbulent fluxes, which requires availability of flux components in all three directions.Such data were available for Kumpula, Hyytiälä, MM and Sorø.

Coordinate rotations and vertical wind speeds
The SPF method is usually applied to selected wind direction intervals such that the coefficients are independent for sectors.Fig. 1, left panels, show the SPF method coefficients for Sorø site; the coefficients were fitted to 22.5°sectors.At Sorø the intercept coefficient b 0 was mostly within ± 0.05 m s −1 and not different from zero if the 95% confidence intervals were considered.Only at around 280°a larger value of b 0 occurs.The coefficients b 1 and b 2 were relatively small.At Hyytiälä the coefficient b 0 showed some variation with wind direction, being higher for the northern direction.Apart from that direction, the intercept was close to zero within the uncertainty intervals.The coefficients b 1 and b 2 were larger in magnitude compared to Sorø.We will discuss the relation of this variation to topography in the following section.According to regression statistic, the coefficient of determination, the multilinear fit showed moderate-to-good fit and the linear regression models were significant (p-value < 0.05) for all direction intervals at Hyytiälä.For Sorø the coefficient of determination was low for almost all directions and the fitted model turned to be insignificant for a range of wind directions.This can be associated to lower number of observations.Also, if the terrain is flat then no coordinate rotation is essentially needed and the model turns to be insignificant.Another site with low R 2 and insignificant planar fit regressions in wide range of directions was MM (Figs.A3).We tested sector wise fitting of the re- (reducing the number of coefficients) to MM and Sorø data, which however did not help to improve significance of the fit.
At Renon site the variation range of all regression coefficients was much larger (Fig 1, note different vertical scale compared to Sorø and Hyytiälä plots).In particular, the coefficient b 1 experienced strong transition at 45-90 o direction.Note that in this direction interval the flow field can be affected by the tower structure (Table 2).The correlation statistics indicated also weak relationship between the x-y wind components and the vertical wind at that direction, and the regression  LGR, GGA-911-0011-0004: 0.35, 0.69, 1.22, 1.88, 2.89 m Boom orientation to the north appeared insignificant.This was partly due to lower number of observations.
Note that, due to the site complexity, the fraction of stability classes depended strongly on wind direction at the Renon site.Therefore, also the regressions established for fixed wind direction intervals were dominated by certain stability classes.If the fitted planes were dependent on the stability, then the obtained regressions were biased towards the dominating stabilities.Dependence of PF on stability will be further studied in Section 4.2.
Fig. 2 plots the residuals W W ¯¯PF of the SPF across wind direction and time of day or stability planes.For Sorø and Hyytiälä the residuals did not show clear dependence on stability in any of the directions indicating that the plane approximated by the SPF was independent on stability.However, some pattern of dominant positive and negative residuals emerged on the diurnal variation (vs.wind direction) plot.For the Renon site the residuals tended to be negative at the stable conditions.In diurnal variation patterns clear positive "domains" existed during the day-time and at nights negative values dominated.The pattern was however wind-direction dependent and indicated dominating upward and downward vertical motions due to topography.The prevailing slope at the site is towards S-SE and we see dominating diurnal variation in the respective wind direction range.The positive upward wind speeds were dominant for day-time hours.Subsidence dominated at nights but also during the day-time negative vertical winds prevailed in the direction from 45 to 90°(can be due to disturbance of tower on wind field).Note that the residuals seemed to exhibit stronger patterns on wind direction vs. hour plots compared to stability.This suggested that the vertical motion pattern is not directly related to local stability but rather to development of diurnal flow patterns, which can be attributed to boundary layer state.Also, the weaker stability dependence appeared mainly due to unstable condition.This suggests that upslope winds in unstable conditions were much less important than downslope winds, driven by catabatic flows.During the study period, the meteorological conditions at Renon were dominated by a very persistent local slope wind system with upslope S-SW winds during the day and downslope NNW winds during the night.This could be observed also on the patterns shown in Fig. 2  We performed also CPF fitting according to Eq. ( 4).We used least squares liner fitting with n n , sin( ), cos( ), , sin ., ) as the independent variables, enabling to obtain also the confidence intervals of the estimated coefficients via OLS regression statistics.For all sites we performed first fitting with number of harmonics N = 8, then excluded the higher harmonics which did not have any significant values and performed fitting again.Table 3 presents the final number N used for CPF at each site.
Additionally, we performed fitting of the CPF method assuming no wind direction dependence of the intercept term (taking in Eq. ( 4)), denoting this approach as CPF0.The coefficient would represent the anemometer drift term if no other mechanisms of mean vertical motions independent of the mean wind would exist.In reality we know that flow disturbances and meteorology such as day-time upward motions, large-scale subsidence or night-time gravity-induced down-slope flows can induce the non-zero intercept.The coefficient α 0 appeared significant for all sites except Hyytiälä and Sorø (Table 3).
Table 3 presents also the root mean square deviations (RMSD) of the residual vertical wind obtained for the CPF method.We tested and found that that the RMSD values were nearly the same for all methods including CPF0 (not shown).However, significant variation among sites existed.

Dependence on stability and time of day
To further evaluate dependence of the PF plane on stability we classified data into stability classes and performed fitting of CPF on each class separately, by using the constant intercept coefficient with = 0 n 0, and = 0 n 0, in Eq. (4).For Sorø, the intercept terms α 0 were identified for the stability classes as following: 0.042, 0.033, 0.029, −0.005, 0.009 for Strongly Unstable (SU), Unstable (US), Near Neutral (NN), Stable (ST), and Strongly Stable (SS), respectively, and for Hyytiälä 0.037, 0.015, 0.010, −0.001, −0.012.For Renon, the following average intercept terms were determined with wider range of variation from unstable to stable classes: 0.057, −0.025, −0.028, −0.054, −0.067 for SU, US, NN, ST, SS, respectively.The coefficients indicate variation of the average vertical speed from negative under stable conditions to positive under unstable.
The wind direction variation of the slope term results from topographically and land-cover heterogeneity induced flow patterns.For Sorø the variation of b S with wind direction is small and might result from slightly tilted positioning of the sonic anemometer with respect to flow streamlines.For Hyytiälä the pattern looks also fairly sinusoidal with minor modulation by the terrain effects.The average pattern of b S for Renon is more complex.The measurement tower at Renon is located on a slope in north-south direction, which implies positive b S values with southern winds and negative with northern.The observed pattern is however more complex, being influenced by local factors in the vicinity of the tower such as pasture located 60 m upslope to the north of the main tower.
In case of Hyytiälä, the slope term was systematically lower for the Fig. 2. Residuals (W W ¯¯PF , in m s −1 ) after fitting the SPF method.The upper plots represent dependence on wind direction and diurnal variation, the lower plots on stability on vertical axis.Left panels are for Sorø, middle for Hyytiälä and the right panels for Renon.

Table 3
The number of harmonics N chosen for the CPF methods and the root-meansquare-deviation (RMSD) of residuals.The intercept α 0 applies to CPF0 (the method with only one intercept coefficient invariant of wind direction).SU class than for the other classes (Fig. 3).We investigated this further in detail and identified that the regression coefficients were not significant in SU case.This implied that the larger value of b I compensated for lower b S and therefore the SU case could not be considered significantly different than other stability classes.
The observed dependence of the intercept coefficient on stability can bias the estimated coordinate frame in case of imbalanced data distribution in stability classes.In practice it is difficult to judge what is the most relevant coordinate frame but one could argue that the flow under neutral stratification conditions would set the best frame for flux calculations.We performed also fitting of the CPF method by choosing the weights of observations such that the distribution of observations in stability classes for each 22.5°wind direction interval corresponded to that of overall dataset at each site.The results for CPF with balanced stability classes are presented for Renon in Fig. 4. Consistently with previous results the slope coefficient b S did not differ for two fitting approaches.However, the intercept term b I appeared systematically different for wind directions from approx.135 to 225°.This was the interval with largest imbalance in stability distribution (Fig. 1).We did not observe such a systematic difference for any other site.

Mean vertical velocity and advective term
Different fitting methods can result in differences in vertical tilt angles (as derived from Eqs. ( 1), ( 3) and ( 5)).The biggest difference can exist between DR and other methods if the intercept terms are significant.This is due to the fact that the PF methods apply the residual term as the offset and remove from the vertical wind.The DR method does not account for the offset and can result in systematically different tilt angles for flux calculation (Fig. 5).The differences for Sorø and Hyytiälä are relatively small, but for Renon it can be up to 20°in unfavourable direction (around 90°).
The previously observed diurnal variation of the intercept term of coordinate fitting results in diurnal pattern of the residual term calculated as W W ¯¯PF (Fig. 5).The systematic average diurnal pattern in residual vertical wind exists at all three sites, being the smallest in  The residual vertical wind has been interpreted as the advective motion in vertical direction and related to vertical advection in CO 2 exchange (e.g., Lee, 1998).Fig. 6 shows the average diurnal variations of flux components for three sites, and dependence of the components on u * at night.For Hyytiälä, the vertical advection compensated for the net ecosystem exchange (NEE) underestimation at low turbulence conditions.This has been previously shown by Mammarella et al. (2007).For Sorø the turbulent flux is approximately constant with u * down to 0.3 m s −1 .The vertical advection term seems not to improve the NEE estimate neither at Sorø nor at Renon sites.At the Renon site the vertical advection term derived from the SPF method seemed to over-compensate at low turbulence conditions (u * < 0.2 m s −1 ).We tested also if the CPF0 method, which does not allow for wind direction dependent variation in intercept, results in different flux estimates compared to the SPF method (not shown).The observed difference was marginal, even though elevated vertical advection estimates were observed at Sorø and Renon sites under low u * conditions.
When looking at the wind direction dependence of the covariance and vertical advection terms, the advection term exhibited large variability with wind direction at both sites (Fig. 8), in particular at Renon.It was difficult to evaluate whether this variability in VA was truly dependent of wind direction as indicated by Fig. 8 or partly induced by large uncertainties in estimation of the term.But it was quite evident that accounting for the vertical advection in NEE estimation without considering the horizontal advection term does not improve the NEE estimates at Sorø and Renon.But also at Hyytiälä, where on the average the correction seemed to work (Fig. 6), the wind direction variability did not allow to make conclusion that the method properly accounted for underestimation by the EC method or was artificially obtained by cancellation of errors over wind directions.

Random and systematic differences
The linear fit statistics were used to evaluate the relationship of fluxes rotated to different coordinate systems (Fig. 9).Here MM was selected as the third site because for Renon the three flux components in unrotated coordinate system were not available.In addition, Kumpula data was used for the analysis in this section.We performed the orthogonal linear least square fitting (called also as the total least squares) which makes no assumption on the dependent and independent variables and the coefficient of the inverse relationship is given just by the inverse of the slope (P1 in Table 4).
At Sorø and Hyytiälä, the correlation between the flux estimates was much higher than at MM site (Fig. 9).In Table 4 we present also the RMSE of the fit residuals, which are also largest for MM.Note that larger systematic difference (slope P1 different from one) exists only between the DR and other methods, being largest for Kumpula.
The linear regression statistics were used also to identify outliers that deviate more than expected at pre-defined confidence level (set to 95%) (Fig. 9).Note that better correlation did not imply necessarily smaller amount of outliers because also narrower confidence intervals were obtained in such case.Nevertheless, the approach would serve useful to identify outliers that originate from the choice of the coordinate system.We also performed flux quality calculations according to flux stationarity criterion (Foken and Wichura, 1996) for Hyytiälä site.This was based on the relative flux error estimate and labelled mostly small-in-value fluxes as non-stationary, in general not overlapping with the outliers as identified by linear fit statistics.
The intercepts of the planar fitting are treated usually as the offsets and subtracted from the vertical wind component to obtain the vertical wind speed in anemometer's coordinate system.Therefore, this assumption changes the average plane into which the coordinates are rotated to calculate the vertical EC flux.The vertical tilt angles obtained by using the SPF and CPF methods were systematically different from the DR approach for several sites, in particular the more complex sites Coweeta, MM and Kumpula (Figs.A2, A3, A1, respectively), but also Sodankylä and Abisko-Stordalen (Figs.A5 and A7).At Abisko-Stordalen the deviation at around 180°can be the result of tower disturbance on the flow.
The substantial difference in vertical coordinate rotation angles is expected to affect the EC flux estimates in a systematic way.The vertical component of the turbulent flux, after second rotation, can be expressed as where the terms on the right-hand side are in the sonic anemometer coordinate frame and the term + u c v c cos sin represents the hor- izontal flux after rotation to a system aligned with the main wind direction.It is a composite of vertical and horizontal flux components and systematic difference in the tilt angle ϕ determines how much the vertical and horizontal turbulent flux components contribute to the final flux after rotation.If the coordinate system is incorrectly defined, the vertical turbulent flux after rotation is contaminated by the horizontal turbulent flux, since the two flux components are correlated.Furthermore, the horizontal turbulent flux is larger in magnitude by a factor of 2 or more, depending on the scalar in question.
Due to sensitivity of flux estimates on the tilt angle (see Fig. B1 as an example for Kumpula) and systematic differences in tilt angles as derived from different methods (Fig. 5), also systematic differences were expected in the average or cumulative fluxes.Table 5 shows the cumulative fluxes for four sites.At Kumpula, the DR method resulted in fluxes that were 5 to 10% higher compared to other methods.At Hyytiälä and other sites the respective systematic differences were smaller.The variants of the PF methods resulted in cumulative fluxes that were as a rule all very close.Deviations from this were the CO 2 fluxes at Kumpula urban site (CPF results deviating from SPF), the CO 2 fluxes at MM site where the difference was up to 2.2%, and the H 2 O fluxes at Sorø site with differences up to 2.7%.

Uncertainty of the coordinate plane
Fig. 10 compares the aggregated over all wind directions distribution statistics of the difference between the tilt angles as obtained from DR and SPF methods.In addition to the systematic differences of the median values from the zero line, which are consistent with the coefficients α 0 reported in Table 3, the significant spread of the difference at Coweeta and Renon sites indicates large variability of deviations with wind direction as well as time of day.For other sites the variation range between the 10th and 90th percentiles remained below 5°.
The OLS regression approach enables to calculate the confidence intervals for the dependent variable i.e. for W ¯in case of PF approach.The confidence intervals can be interpreted as the error estimates of the long-term average vertical wind speed in anemometer's coordinates system (and could be used in determination of the uncertainty of fitted coordinate plane).The error in W ¯further translates into uncertainty in vertical advection if the wind speed relative to PF plane is used in calculation.Fig. 11 presents the obtained average errors (as presented by 95% confidence intervals) in vertical wind speed.Note that the errors presented here are the confidence bounds for the fitted functions only (i.e.Eqs. ( 2) and ( 4)) and do not contain the estimates of the random uncertainty for each observation.The CPF method enables to obtain lower uncertainty in predicted vertical winds (and thus also in coordinate system defined by the fitted plane).This is in part due to less parameters in the model (if compared with the SPF over all wind directions) but also because the method is not as sensitive as the SPF method to the limited observations (or other causes of fitting uncertainty) in particular wind direction intervals.

Discussion and conclusions
Commonly used methods for coordinate rotation in EC calculations are the double rotation (DR) and planar-fitting (PF) methods.We evaluated the commonly used sector-wise PF method as well as two variants of the continuous PF method where the wind direction variation of the intercept and wind-speed dependent coefficients were represented as the continuous functions in the form of Fourier series.The fitting was performed by using OLS approach enabling to obtain the uncertainties of the fitted coefficients and choose the appropriate number of harmonics for the CPF for each site.In general, larger number of harmonics was relevant to more complex sites.
The residual variation of vertical winds was very different at the studied sites (Table 3).The random variability of average vertical wind speed should be in the order of 0.01-0.02m s −1 depending on the conditions.Over rough surfaces such as forests the random uncertainty is expected to be larger and over smooth surfaces and low measurement heights smaller.Indeed at Abisko-Stordalen, which was a peatland site with 2.2 m measurement height, the RMSD of the residual vertical wind was within the expected limits.For Renon the RMSD was the highest but for other sites there seemed to be no clear pattern of the residual variance according to site complexity.The vertical wind velocity relative to the PF coordinates measured above canopy was observed to be predominantly negative at night and positive during the day at several sites (see Froelich et al., 2005 for MM andMammarella et al., 2007 for Hyytiälä).This is also seen in Fig. A3 for MM and in Fig. 6 for Sorø, Hyytiälä and Renon.For other sites the average diurnal variation in W ¯was less evident.
We also analysed the vertical advection terms derived from the residual vertical wind speed relative to PF coordinate systems.At Hyytiälä, where on the average the vertical advection seemed to improve the night-time NEE estimates (Mammarella et al., 2007), high variability with wind direction was observed.In the presence of steep vertical gradients in CO 2 concentration, which are frequently realized at night, positive vertical CO 2 advection fluxes were observed at MM when average vertical wind was negative (Su et al., 2008).This was in agreement with Fig. A3, which additionally indicated extreme variability of the advection term with friction velocity in the wind direction interval from approx.160 to 225°.For other forested complex sites (see Fig. 9. Comparison of CO 2 fluxes at three sites after rotating into different coordinate systems.Left panels are for Sorø, middle panels for Hyytiälä and right panels for MM site (units in µmol m −2 s −1 ).Linear fits are performed without the assumption of independent and dependent variables.Percentage of outliers (red circles) according to linear regression statistics for Sorø are 2.83, 3.51 and 3.04% (from top to down), for Hyytiälä 1.33, 5.47 and 0.63% and for MM 3.35, 1.74 and 1.64%, respectively.Additional quality screening was performed for flux stationarity by using the threshold value 1 (Foken and Wichura, 1996).Data points with flux stationarity exceeding 1 are labelled with magenta circles (Hyytiälä only).

Table 4
Orthogonal least squares linear regression slope coefficients (P1) and rootmean-square-errors (RMSE) for CO 2 fluxes at Kumpula, Hyytiälä, MM and Sorø sites.Fitting was performed with no constant term.The coefficients P1 are presented with respect to SPF method and with as many decimal digits as relevant considering 95% confidence intervals of the coefficient (1.00 denotes the coefficients which did not differ significantly from one).CPF0 denotes the variant of the CPF method with intercept coefficient independent of wind direction.and A6).Thus, vertical advection estimates derived from residual vertical velocity cannot be attributed to NEE with confidence at any of the sites without considering also the horizontal advection.Similar conclusion was made by Aubinet et al. (2010) who studied in detail the vertical as well as horizontal advection at three sites with different topographies, including Renon and Norunda.They also concluded that accurate measurement of full mass balance is difficult due to large uncertainties related to sampling.The uncertainty in determination of coordinate frame obtained by the PF approach contributes to such uncertainties and should be minimized as presented in the current study.
Our results further suggest that cumulative fluxes obtained by DR method can be systematically biased up to 5% in most cases, but could be even 10% depending on site and flux in question.This is related to the intercept term, and when different from zero, the DR method relies on a coordinate system that systematically deviates from the PF methods.The difference and its direction depend on the sign of the intercept, but the average vertical tilts defined by the DR and PF coordinate systems have been observed to be different from a few to up to 10°at very complex sites.The sites showing several degrees differences between the DR and PF coordinate systems at some direction intervals were Renon (Fig. 6), Coweeta (Fig. A2), MM (Fig. A3) and Kumpula (Fig. A1).For less complex sites, the difference was in general not more than a few degrees.However, at the ideal Abisko-Stordalen site, a large difference could be noticed at the direction interval 135-225°(Fig.A7).Fig. 10 summarizes the statistics for the deviations between the tilt angles defined by the DR and SPF methods.For complex sites the median deviates the most from zero.Two sites with large variation range of the differences are Coweeta and Renon.For MM the range between the first and last deciles remains below 5°, as for the sites classified into categories other than high complexity.The large variation range of the rotation tilt angles obtained from DR method at the complex sites is not feasible.Therefore, we conclude that the PF method is more relevant for defining the coordinate system for flux calculation at complex sites.
Systematic differences between the fluxes obtained from the variants of the PF coordinate systems are in general small.We also noticed that the linear relationships in Eq. ( 2), when applied sector-wise and defining the SPF coordinate system for flux calculation, can be statistically weak and/or with low coefficients of determination.Insignificant relationships can be the result of a limited number of observations at certain wind directions or indicate weak relationship between the mean vertical and horizontal winds.Therefore, we suggest to use the linear fit statistics for evaluation of significance of the PF.
We also studied how the unbalanced distribution of observations in different stability classes within wind direction intervals can affect the PF results.The difference was evaluated between the fittings performed to original dataset and the one with stability-balanced weights (for details see Section 4.2).The only site where we observed a systematic difference in the obtained intercept term was Renon (Fig. 4) where the flow patterns and directions are strongly affected by topography.Even though strong stability dependence of the coordinate systems seems not common among sites, we argue that similar test serves useful at complex flux sites.If not done, the fitted coordinate frame might be biased towards prevailing stability conditions, leading also potentially to biased coordinate frame for flux calculations.
The regression statistics between the flux estimates rotated to different coordinate systems can be used also to identify outliers resulting from the choice of the coordinate system.The outliers can be diagnosed by using the residual intervals of the flux observations calculated in different coordinate planes.Such outliers appeared different than those flagged by flux stationarity criteria (e.g., Foken and Wichura, 1996).Flux stationarity is essentially a relative measure and therefore the fluxes small in absolute value tend to be flagged as non-stationary by this criterion.The outlier diagnosis by regression statistics helps to identify deviating observations and can be used in evaluation of respective flux uncertainties.The uncertainty arising from the choice of the coordinate system is not typically considered in evaluation of uncertainties of EC fluxes (e.g., Mauder et al., 2013).
In summary, we recommend that eddy covariance practitioners choose the simplest PF coordinate system.The PF coordinates should be used instead of DR because of the potential systematic difference in coordinate frames, and as a consequence in fluxes, if non-zero vertical wind component exists that is not correlating with horizontal motion.The PF frame defining the coordinate system for flux calculation should  be based on significant statistical relationships.We recommend to use the CPF method in cases when the number of observations in some wind direction interval is low or the obtained fit is insignificant due to large variance in measurements.The method enabled to obtain the PF with lower uncertainty as compared to the SPF method and therefore also the coordinate frame defined by the method is less uncertain.The high residual variation in vertical wind can be indication of low-frequency motions.Such motions can contribute to vertical exchange as proposed by Finnigan (2004) and the significance of such transport mechanism should be evaluated for each site separately.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Dependence of SPF coefficients on wind direction.The upper subplots were obtained when fitting all 3 coefficients to data in 22.5 o sectors.Middle subplots: The coefficient of determination (R 2 ) and significance (p-value) of the linear fit.Lower subplots: Number of observations for unstable (US, −1.5< ς<−0.3),neutral (N, −0.3< ς<0.3) and stable (S, −0.3< ς<1.5) stratification classes.Top left panels are for Sorø, top right for Hyytiälä and the lower panels for Renon.
in direction intervals from 135 to 225 o .Day-time positive values of residuals in W were observed mostly at around 180-315 o with parallel negative regions at nights.

Fig. 4 .
Fig. 4. CPF method applied to Renon data by using the observed distribution of observations and with balancing weights in 22.5 o sectors such that the distribution would correspond to that over all directions in stability classes US (−1.5< ς<-0.3),N (N, −0.3< ς<0.3) and S (−0.3< ς<1.5).

Fig. 6 .
Fig. 6.Turbulent CO 2 flux F, the storage change flux ST, the vertical advection VA, and their sum (Tot = F + ST + VA) for Sorø, Hyytiälä and Renon sites.Diurnal variation with VA term was calculated using the SPF method (upper plots) and variation of night-time fluxes with friction velocity (lower panels).All fluxes are in µmol m −2 s −1 .

Fig. 7 .
Fig. 7. Turbulent flux F, vertical advection VA (using SPF residuals), and total (= F + ST + VA) for three sites.Variation with wind direction and time of day.All fluxes are in µmol m −2 s −1 .

Fig. 8 for
Fig.8for Renon and A2 for Coweeta) the night-time vertical advection term showed also very high variability.The same applied to the sites with medium and low complexity (for example Sorø and Norunda, see Figs.8 and A6).Thus, vertical advection estimates derived from residual vertical velocity cannot be attributed to NEE with confidence at any of the sites without considering also the horizontal advection.Similar conclusion was made byAubinet et al. (2010) who studied in detail the vertical as well as horizontal advection at three sites with different topographies, including Renon and Norunda.They also concluded that accurate measurement of full mass balance is difficult due to large uncertainties related to sampling.The uncertainty in determination of coordinate frame obtained by the PF approach contributes to such uncertainties and should be minimized as presented in the current study.Our results further suggest that cumulative fluxes obtained by DR method can be systematically biased up to 5% in most cases, but could be even 10% depending on site and flux in question.This is related to the intercept term, and when different from zero, the DR method relies on a coordinate system that systematically deviates from the PF methods.The difference and its direction depend on the sign of the intercept, but the average vertical tilts defined by the DR and PF coordinate systems have been observed to be different from a few to up to 10°at very complex sites.The sites showing several degrees differences between the DR and PF coordinate systems at some direction intervals were Renon (Fig.6), Coweeta (Fig.A2), MM (Fig.A3) and Kumpula (Fig.A1).For less complex sites, the difference was in general not more than a few degrees.However, at the ideal Abisko-Stordalen site, a large difference could be noticed at the direction interval 135-225°(Fig.A7).Fig.10summarizes the statistics for the deviations between the tilt angles defined by the DR and SPF methods.For complex sites the median deviates the most from zero.Two sites with large variation

Fig. B1 .
Fig. B1.Average diurnal variation of turbulent fluxes (left panels) and sensitivity to vertical tilt angle (right panels) for Kumpula site.

Table 2
Measurement sites and instrumental setups.

Table 5
Cumulative fluxes for four sites over the periods indicated in Table2, obtained by using different methods of coordinate rotation.The percent difference is calculated with respect the SPF method.No gap-filling of missing data was performed.