Quantifying the uncertainty of eddy covariance fluxes due to the use of different software packages and combinations of processing steps in two contrasting ecosystems

We have carried out an inter-comparison between EddyUH and EddyPro®, two public software packages for post-field processing of eddy covariance data. Datasets including carbon dioxide, methane and water vapour fluxes measured over 2 months at a wetland in southern Finland and carbon dioxide and water vapour fluxes measured over 3 months at an urban site in Helsinki were processed and analysed. The purpose was to estimate the flux uncertainty due to the use of different software packages and to evaluate the most critical processing steps, determining the largest deviations in the calculated fluxes. Turbulent fluxes calculated with a reference combination of processing steps were in good agreement, the systematic difference between the two software packages being up to 2.0 and 6.7 % for half-hour and cumulative sum values, respectively. The raw data preparation and processing steps were consistent between the software packages, and most of the deviations in the estimated fluxes were due to the flux corrections. Among the different calculation procedures analysed, the spectral correction had the biggest impact for closed-path latent heat fluxes, reaching a nocturnal median value of 15 % at the wetland site. We found up to a 43 % median value of deviation (with respect to the run with all corrections included) if the closed-path carbon dioxide flux is calculated without the dilution correction, while the methane fluxes were up to 10 % lower without both dilution and spectroscopic corrections. The Webb–Pearman– Leuning (WPL) and spectroscopic corrections were the most critical steps for open-path systems. However, we found also large spectral correction factors for the open-path methane fluxes, due to the sensor separation effect.


Introduction
The eddy covariance (EC) technique is the most direct and defensible way to measure and calculate vertical turbulent fluxes of momentum, energy and gases between the atmosphere and biosphere.During the last 3 decades, the number of long-term EC stations all over the world has increased exponentially, covering a wide range of different ecosystem types (FLUXNET, www.fluxdata.org).EC is a technique analysing high-frequency wind and scalar atmospheric data series (often called "raw data"), usually saved in hard drive devices for post-field processing and final estimations of turbulent flux values.During the past years several attempts have been made to standardise the processing methodology, at least for carbon dioxide (CO 2 ) and sensible and latent heat (LE) fluxes (Aubinet et al., 2000(Aubinet et al., , 2012;;Lee et al., 2004).However, the harmonisation of data processing is quite difficult, since most of the required steps and corrections are site-specific and instrument-specific (gas analyser and sonic anemometer).Nowadays, new and better instrumentation is available for measuring turbulent fluxes of energy and matter using the EC technique.Recent studies have compared commercially available gas analysers, focusing on precision, stability and systematic and random errors both for CO 2 fluxes (Burba et al., 2008;Ibrom et al., 2007a;Järvi et al., 2009) and methane (CH 4 ) and nitrous oxide (N 2 O) fluxes (Detto et al., 2011;Peltola et al., 2013Peltola et al., , 2014;;Rannik et al., 2015).However, only a few studies have reported an inter-comparison between EC software packages, and all focusing only on energy and CO 2 fluxes (Fratini and Mauder, 2014;Mauder et al., 2007Mauder et al., , 2008)).For example, Mauder et al. (2008) concluded that the data preparation, the coordinate rotation of Published by Copernicus Publications on behalf of the European Geosciences Union.sonic anemometer wind components and the different approaches for high-frequency spectral correction are critical processing steps, giving differences up to 10 % in their study.Fratini and Mauder (2014) compared TK3 and EddyPro ® software packages, achieving a satisfying agreement in calculated fluxes and related quality flags only after tuning the software processing steps and corrections to be similar.In fact, systematic differences in EC flux estimates strongly depend on the selection, application and order of processing steps, and the correct application, order and sometimes relevance and consequences of several processing steps are still topics under discussion (Aubinet et al., 2012;Mauder and Foken, 2006).In addition, the relevance of some processing steps and corrections depends not only on the system set-up, but also on meteorological conditions and ecosystem types (Mammarella et al., 2015;Nordbo et al., 2012).As a result, the EC processing software packages available to the community feature different implementations: some steps may be implemented using different methods (Mauder et al., 2007), while some operations and eventually further corrections suggested by recent findings are not supported by some of the software packages.This is particularly relevant for gases like CH 4 and N 2 O, for which the deployment of the EC system with easy-to-use fast response analysers have become popular only during the last decade.Therefore, neither data processing approaches have yet been standardised nor software inter-comparison studies have been published for these gas fluxes.
In this study, we have performed an inter-comparison between EddyUH and EddyPro, two public and commonly used software packages for EC data processing and calculation.The aims are to estimate the flux uncertainty due to the use of different software packages for half-hour as well as for cumulative sums, and to assess the most critical processing steps, determining the largest deviations in the calculated fluxes.We focus not only on LE and CO 2 fluxes, as it has been done in previous studies, but also on CH 4 fluxes.

Site description and measurements
The software inter-comparison was performed using datasets from two field sites in southern Finland.The first dataset was collected at the Siikaneva fen site (61 • 49.961 N, 24 • 11.567 E) during the CH 4 inter-comparison field campaign (Peltola et al. 2013).The EC data used in this study were measured during 1 May-30 June 2010 with a 3-D sonic anemometer (USA-1, Metek GmbH), two closed-path gas analysers (LI-7000, LI-COR Biosciences; G1301-f, Picarro Inc.) and one open-path gas analyser (LI-7700, LI-COR Biosciences).LI-7700 was an early prototype version of the later commercialised LI-7700.LI-7000 measured CO 2 and water vapour (H 2 O) and G1301-f measured CH 4 and H 2 O molar fractions, whereas LI-7700 measured CH 4 molar concentrations.LI-7000 and G1301-f used a shared heated sampling line that was approximately 16.8 m long (ID: 10 mm, flow rate: 24 L min −1 ).The sonic anemometer was situated 2.75 m above peat surface and the open-path LI-7700 directly below it, causing a 45 cm vertical separation between the sensors.Further details about the site and measurements can be found in Peltola et al. (2013).
The second dataset was collected between 1 July and 30 September 2010 at the Erottaja site located in Helsinki city centre (60 • 09.912 N, 24 • 56.723 E).The measurements represent densely built-up urban area, with only 5 % of the surface being covered with vegetation.The measurements are carried out in a 3.8 m high mast located on top of a 38 m high fire station tower, resulting in a total height of 41.8 m.This is a sufficient height for the EC measurements as the mean building height in the surroundings of the tower is 21.7 m.The measurement set-up consisted of an ultrasonic anemometer (USA-1, Metek GmbH) to measure the wind components and open-and closed-path infrared gas analysers (LI-7500 and LI-7200, LI-COR Biosciences) for the CO 2 and H 2 O fluctuations.Details of the measurement set-up can be found in Nordbo et al. (2013).

Turbulent flux calculation
The turbulent fluxes of CO 2 (F CO 2 , µmol m −2 s −1 ) and CH 4 (F CH 4 , nmol m −2 s −1 ) and sensible (H , W m −2 ) and latent (LE, W m −2 ) heat are calculated from the covariances between a respective scalar and vertical wind velocity (w) as where ρ d is the dry air density (kg m −3 ), c p the specific heat capacity of dry air (J kg −1 K −1 ), L v the latent heat of vaporisation for water (J kg −1 ), T the temperature (K) and M a and M w the molar masses of dry air and water, respectively.The terms w T , w χ CO 2 , w χ CH 4 and w χ H 2 O are the covariances between w and T, dry mole fractions of CO 2 , CH 4 and H 2 O, respectively.Overbars and primes represent temporal averaging and fluctuations, respectively.

Set-up of software runs
The results presented in this study are The eddy covariance datasets were processed using the reference combination of processing steps (Fig. 1) and available methods implemented in EddyUH and EddyPro (Table 1).The applied methods were the same for most of the steps.However, some differences between software packages were present.In EddyUH the raw data despiking was done by the difference limit method (Appendix B), while in EddyPro the Vickers and Mahrt (1997) method was used.Different experimental methods were applied for spectral correction, e.g. according to Mammarella et al. (2009) in EddyUH and Fratini et al. (2012) in EddyPro.Moreover, additional correction for water vapour cross-sensitivity (henceforth point-bypoint spectroscopic correction) of closed-path CH 4 was done in EddyUH according to Rella (2010).The same correction is not available in EddyPro.Finally, we performed four different combinations of processing steps and compared them with the reference combination in order to evaluate the effects of different calculation steps on the final flux estimates.The alternative runs were set up by modifying one step of the reference combination at a time.The first and second runs were done by excluding the spectral correction and applying the theoretical approach (instead of experimental one used in the reference combination), respectively.The Webb-Pearman-Leuning (WPL) (Webb et al., 1980, henceforth also called dilution correction in the case of closed-path sensors) and spectroscopic corrections were omitted in the third run, and finally in the final run a constant value for the time lag was used.
Flux data were quality-screened prior to analysis.CH 4 flux data were removed if the CH 4 mean mole fraction was above 5 or below 1.7 ppm.Additionally, LI-7700 fluxes were discarded if the received signal strength indicator (RSSI) was below 15.All the Siikaneva flux data were discarded if the second coordinate rotation angle (used to set w = 0) was above 10 • .Wind direction (90-180 • omitted) was used to omit periods when the measurement system at the Erottaja site was in the wake of the building.In addition, periods when there were known problems with the measurement setup at the Erottaja site were discarded.Plausibility limits were also used, since if the flux values were outside certain predefined limits they were thought to be erroneous.For Siikaneva data these limits were −50 and 160 nmol m −2 s −1 for F CH 4 , −30 and 600 W m −2 for LE and −20 and 20 µmol m −2 s −1 for F CO 2 .For Erottaja the following limits were used: −30 and 500 W m −2 for LE and −10 and 60 µmol m −2 s −1 for F CO 2 .Finally, flux data were discarded if the corresponding quality flags, as determined by EddyUH and EddyPro, were above 5 based on the Foken et al. ( 2004) flagging policy.The data coverage obtained after data screening is given in Table 2.

Flux comparison between software packages
Fluxes measured by the same system were analysed and compared as estimated by the two software packages in the reference run.In terms of regression statistics, a very good agreement between EddyUH and EddyPro was obtained for LE and F CO 2 measured by LI-7000 (Fig. 2d and e), LE and F CH 4 by G1301-f at Siikaneva (Fig. 2c and a) and F CO 2 by LI-7200 at Erottaja (Fig. 2i).For LI-7500 LE and F CO 2 no significant systematic differences between the software packages were found (Fig. 2f and h) even though the data show more scatter around the 1 : 1 line (r 2 = 0.91, RMSE = 1.29 W m −2 for LE and r 2 = 0.98, RMSE = 0.12 µmol m −2 s −1 for F CO 2 ) than the other analysers.Instead, a systematic difference was found for LE measured by the LI-7200 analyser in Erottaja, the EddyPro fluxes being 2 % higher than those calculated by EddyUH (Fig. 2g).
Finally, good agreement resulted from the LI-7700 F CH 4 , the slope being equal to unity (r 2 = 0.98 and RMSE = 0.88 nmol m −2 s −1 ; Fig. 2b).However, the sensor separation correction in EddyPro (Horst and Lenschow, 2009) caused relatively large scatter between LI-7700 fluxes, and if the correction was omitted, visually the scatter between the two software packages was reduced, although regression statistics were slightly worse (y = 0.94 x+ 0.51, RMSE = 1.21 nmol m −2 s −1 , r 2 = 0.99) (figure not shown).Applicability of the sensor separation correction for this particular dataset is discussed in Sect.4.2.
In order to further evaluate the discrepancies between the two software packages, diel patterns of flux ratio (left side in Fig. 3) and bias (right side in Fig. 3) were plotted for each flux at different processing levels.In general, the uncorrected (raw) fluxes do not show significant deviations from unity (left side of Fig. 3) or from the zero line (right side of Fig. 3), which suggests that the preparations done at the raw data level (despiking, coordinate rotation, time lag compensation) did not make a significant systematic difference to the fluxes.For LE calculated from LI-7500 data, the WPL correction tends to be slightly larger in EddyUH than in Ed-dyPro, meaning that it increases the daytime positive fluxes more (cf.WPL curve in Fig. 3r).The daytime WPL correction in EddyUH is approximately 2 W m −2 larger than in EddyPro, which corresponds to a daytime relative differ-

Open-path analysers
Closed-path analysers Relative magnitude definition: . ence of 4 % (Fig. 3q).The reason was found to be a difference in one of the inputs of the WPL temperature fluctuation term (T -term), specifically in the water vapour density (ρ w ).In fact, while EddyUH uses the one from LI-7500, EddyPro calculates ρ w from meteorological relative humidity (RH) data.For closed-path analysers or other open-path fluxes there is no difference in the WPL correction, and the uncorrected and WPL corrected curves follow similar pattern.The biggest differences between the software packages were related to the spectral corrections.For closed-path LI-7000 and G1301-f the fully corrected (e.g.WPL + SC in Fig. 3) LE fluxes were approximately 7 % larger at nighttime in EddyUH than in EddyPro (Fig. 3m and k).However, during these periods the absolute difference was still below 1 W m −2 , since the night-time LE is small.This difference is due to different RH dependence of low-pass filter time constant found between the two softwares (see discus-sion in Sect.4.1).For LI-7700 F CH 4 the relative difference between the software packages was on average 12 % during the night, EddyPro fluxes being larger (Fig. 3c), which corresponds to an absolute difference of 1-4 nmol m −2 s −1 (Fig. 3d).This difference was related to the sensor separation correction, and the difference was smaller (4-10 % or 0.5-1.5 nmol m −2 s −1 ; EddyUH fluxes were larger) if the correction was not done in EddyPro (see the corresponding discussion in Sect.4.2).

Flux comparison between instruments
Fluxes measured with different gas analysers are compared as estimated using the reference combination in EddyUH and EddyPro (Fig. 4).At Erottaja, very good agreement was found between F CO 2 measured by LI-7200 and LI-7500 (Fig. 4e and f), and similar results . Median diurnal variation of flux ratio (left side) and bias (right side) for the studied instruments and variables between the two software.Light blue shows the uncorrected fluxes after despiking, coordinate rotation, detrending and time lag compensation; blue the WPL corrected fluxes (for G1301-f and LI-7700 F CH 4 this includes also spectroscopic correction); and green the WPL plus spectral corrections (WPL + SC).The WPL + SC curves represent data from the reference run, i.e. data that are fully corrected.
were obtained by using EddyUH (slope = 0.99, r 2 = 0.92, RMSE = 1.25 µmol m −2 s −1 ) and EddyPro (slope = 0.98, r 2 = 0.90, RMSE = 1.32 µmol m −2 s −1 ).LE measured by the same gas analysers show slightly weaker correspondence, the slope and RMSE being 1.02 and 8.15 W m −2 for EddyUH and 1.06 and 8.36 W m −2 for the EddyPro run, respectively (Fig. 4g and h).The small difference in LE between the two software packages is likely due to the spectral corrections (see below).Overall, the results are in agreement with Nordbo et al. (2013), who also found a better agreement between F CO 2 than LE measurements at the same site.A very good correspondence (slope = 0.97, r 2 = 1.00,RMSE = 1.94 W m −2 ) was found between LE measured by LI-7000 and G1301 systems for the EddyUH reference run at the Siikaneva site (Fig. 4c).For EddyPro similar statistics were obtained (slope = 0.96, r 2 = 1.00,RMSE = 2.56 W m −2 in Fig. 4d).By using the reference combination in the EddyUH run, a relative good agreement was also obtained between F CH 4 values measured by G1301 and LI-7700 (Fig. 4a).The regression statistics (slope = 1.09, r 2 = 0.84, RMSE = 2.64 nmol m −2 s −1 ) are similar to the ones reported in Peltola et al. (2013).A 16 % difference be-  tween the two F CH 4 values was found in the EddyPro run (Fig. 4b), because of larger spectral correction estimates for the open-path gas analyser (see Sect. 4.2).

Effects of different flux processing combinations
The impact of calculation steps on the final flux estimates in EddyUH and EddyPro is presented as deviation (in %) from the reference run for Siikaneva (Fig. 5) and Erottaja (Fig. 6) datasets.For closed-path systems in Siikaneva, the run without spectral correction had the largest effect on LE, being (in the EddyUH run) 14 and 9 % lower for G1301-f, and 16 and 12 % lower for the LI-7000 system, during night-time and daytime, respectively (Fig. 5c and d).When compared to the EddyUH run, slightly smaller deviations (10 and 8 % for G1301-f and 12 and 11 % for the LI-7000) were found in the EddyPro run, but the trend in the deviations was consistent between the software packages.A similar range of deviations was found in Erottaja LE measured by LI-7200 (Fig. 6b).However, performing the theoretical spectral correction on the Siikaneva LE had a minimal effect (Fig. 5c and d), while the deviations found in LE measured by LI-7200 at Erottaja ranged between 6 and 3 % (Fig. 6b).The use of a constant time lag produces a small effect on LE, except during nighttime when the higher RH increases the sorption of H 2 O in the sampling line walls (Nordbo et al., 2014), and H 2 O time lag becomes larger than its nominal value.At Siikaneva the theoretical spectral correction produces up to 7 % higher F and F CH 4 measured by LI-7000 and G1301-f systems, respectively, with respect to the experimental spectral correction (used in the reference run).If the LI-7000 F CO 2 is calculated without performing the dilution correction (e.g. using the wet mole fraction), we obtained 43 % higher daytime CO 2 uptake with both software packages, and about 3 % lower positive fluxes during night-time.The same correction also has a relevant impact on daytime F CH 4 measured by G1301-f, resulting in 10 and 6 % deviations in EddyUH and EddyPro, respectively.
For open-path systems the critical step is represented by the WPL (and spectroscopic) correction.In Erottaja, the net CO 2 emission calculated from LI-7500 data without WPL correction is underestimated by 38 and 37 % during daytime in EddyUH and EddyPro, respectively (Fig. 6c).Instead, the nocturnal F CO 2 is 8 % smaller than from the reference run.
Although the effect of no WPL is lower on LE, the calculated deviations are still relevant, being 15 and 12 % during daytime (in EddyUH and EddyPro, respectively), and 4 and 3 % during night-time.Finally, in Siikaneva the LI-7700 F CH 4 calculated in EddyUH without the combined WPL and spectroscopic correction shows −71 and 20 % deviations from the reference run during daytime and night-time.The same median values estimated in EddyPro are −67 and 18 % (Fig. 5b).In addition, it can be seen from the same figure that the deviation of nocturnal F CH 4 calculated without spectral correction (no spec run) in EddyPro is larger (−21 %) than the one obtained in EddyUH (−8 %).This is consistent with the fact that EddyPro gives larger nocturnal spectral correction factors with respect to EddyUH (see Sect. 4.1 and Fig. 3).

Differences between cumulative fluxes
Cumulative sums of non-gap-filled flux time series are mostly within ±2 %, which suggests that there is no significant systematic bias between the two software packages (Table 2).Biggest relative differences were in LI-7700 F CH 4 and LI-7500 H 2 O flux cumulative sums, −6.7 and −5.3 %, respectively, meaning that the cumulative values estimated with EddyPro were somewhat larger than with EddyUH.For LI-7700 F CH 4 , if the sensor separation correction was omitted, the relative difference was 1.0 % (EddyUH fluxes larger).The smallest relative difference was obtained for G1301-f cumulative F CH 4 , 0.03 %.The absolute differences between the cumulative CH 4 fluxes at Siikaneva fen during the period May-June 2010 were −0.02 g (CH 4 ) m −2 (LI-7700, EddyPro larger) and less than 0.01 g (CH 4 ) m −2 (G1301-f) (Table 2).The difference between cumulative CO 2 fluxes was −4 g (CO 2 ) m −2 ; Ed-dyUH was showing slightly higher CO 2 uptake.This originated from the fact that EddyPro estimated approximately 1-3 % higher respiration at night and EddyUH calculated < 1 % higher uptake during daytime (cf.Fig. 3e).These differences were caused by the spectral corrections and they inflicted the observed deviation between the cumulative LI-7000 CO 2 fluxes.At the Erottaja urban site during July-September 2010 EddyUH showed slightly lower cumulative CO 2 emission (−2 g (CO 2 ) m −2 ) for LI-7200, whereas for LI-7500 the difference was larger (−27 g (CO 2 ) m −2 ).Similarly, as in the case of LI-7000 CO 2 fluxes, here the difference observed between LI-7500 cumulative F CO 2 was also likely caused by the spectral corrections (cf.Fig. 3i).The cumulative H 2 O fluxes were within 2 mm.However, the data coverage should be considered when evaluating the significance of these absolute differences.The data coverage of the Siikaneva measurements was between 73 % (CH 4 , G1301-f) and 21 % (CH 4 , LI-7700).At the Erottaja site lower data coverage was obtained (between 37 % (CO 2 , LI-7500) and 29 % (H 2 O, LI-7500 and LI-7200)).

Discussion
We have performed an inter-comparison between EddyUH and EddyPro, two public software packages for EC flux calculation.Both software packages feature up-to-date methods for EC raw data processing steps and corrections.Flux data as estimated by the reference combinations (Table 1) were in good agreement.In general, there were not significant systematic differences in the uncorrected fluxes calculated by EddyUH and EddyPro (Fig. 3).This suggests that an optimal choice for the raw data preparation and processing schemes leads systematic biases to be avoided in the fluxes between the two software packages.The most significant differences between the software packages only occurred after the flux corrections, and the impacts of these steps are further discussed below for fluxes measured by closed-and open-path systems.
I. Mammarella et

Critical steps and recommendations for closed-path systems
EC measurements from three closed-path systems (LI-7000, LI-7200 and G1301-f) were processed, and the impact of different processing step combinations was analysed using runs performed with EddyUH and EddyPro.Among the different calculation procedures analysed, the spectral correction was the most relevant for the closed-path LE measurements at the two sites, the median values being between 6 and 16 %.On average, the use of theoretical spectral correction gave up to 6 % lower LE in Erottaja, while in Siikaneva the deviation respect to the reference run was generally below 3 %.We determined a stronger RH dependence of low-pass filter time constant in Erottaja than in Siikaneva (Fig. 7) caused by the non-heated sampling line there (Nordbo et al., 2013).Moreover, because of different approaches in the spectral correction methods, the low-pass filter time constants estimated by EddyPro in Erottaja were larger than those estimated by EddyUH for RH values lower than 80 % (Fig. 7c).This may explain the 2 % difference between LI-7200 LE as estimated by EddyPro and EddyUH (Fig. 2g).Although the relative magnitudes of spectral corrections are not directly comparable with other studies, they are in the same range as previously reported (Fig. 1).Surprisingly, at Siikaneva the theoretical spectral correction gave on average 7 % higher F CO 2 and F CH 4 than the experimental spectral correction, and the result was consistent between the software packages.A possible explanation could be that in the site-specific co-spectral model used in the reference run, the co-spectral peak frequency n m = 0.056 estimated in unstable conditions is shifted to lower frequencies respect to the Kaimal et al. (1972)-based atmospheric surface layer (ASL) co-spectral model used in the theoretical spectral correction (Moncrieff et al., 1997).In addition, in stable conditions the stability dependence of the estimated n m is less pronounced, and it does not follow the ASL parameterisation strictly (data not shown).This would result in smaller spectral corrections when using the site-specific co-spectral model in Eq. (C2).Generally it is thought that the theoretical approach easily underestimates the high-frequency spectral attenuation for closed-path systems, because of potential variations in the mass flow rates and relative uncertainty in the flow regime (Aubinet et al., 2000).Recent studies, both in the lab and in the field, have also demonstrated the important effects that different dust filters and rain caps (typically mounted at the tube inlet) have on the EC system cut-off frequency (Aubinet et al., 2016;Metzger et al., 2016).Such effects are not included in the theoretical approach.However, we have shown in this study how the use of a site-specific co-spectral model can reverse the relative magnitude of spectral correction calculated with the experimental approach compared to the one obtained via the theoretical approach.
In addition, we found that in Siikaneva the LI-7000 F CO 2 was greatly affected by the dilution effect due to large H 2 O fluxes (average daytime value of LE equals 170 W m −2 ) caused by the wet surface conditions and the presence of vegetation at the site.The same correction affected less F CH 4 Table 2. Cumulative sums estimated with EddyUH and EddyPro.Relative differences and data coverage are also shown.Data were not gap-filled prior to calculation of the cumulative sums.

Siikaneva fen (May-June 2010)
Erottaja urban site (July-September 2010) measured by G1301-f, because the flux to concentration ratio for CH 4 was larger than the one for CO 2 .In Erottaja the same effect on F CO 2 measured by LI-7200 was much lower (on average 2 %), because of the 3 times smaller daytime values of LE at this urban site.Besides this, the sampling line was not heated in Erottaja, and different magnitude of the dilution correction between the two sites was expected.In fact, although the heating (and insulation) of the sampling line decreased the H 2 O low-pass filter effect (especially at increasing values of RH), at the same time it increased the H 2 O fluctuations in the sampled air, leading to larger dilution correction.It is common to think that the WPL and spectroscopic corrections have small importance for closed-path analysers, since the temperature fluctuations are dampened in the sampling tube (Leuning and Judd, 1996;Rannik et al., 1997).However, this depends on the magnitude of H 2 O fluctuations, which, as we have demonstrated here, is related to the ecosystem type and system set-up.Fortunately, current closed-path gas analysers also report H 2 O turbulent signals, and the measured gas mole fractions can be readily converted into dry mole fractions, either in the analyser internal software or in the post-field data processing by using the pointby-point dilution and spectroscopic corrections.Finally, for our sites and datasets, the use of nominal constant time lag was only an issue for nocturnal LE, when the absorption effect on H 2 O in the sampling line became more relevant, determining an increase of H 2 O time lag, and thus flux underestimations of 2 and 7 % in Siikaneva and 3 % in Erottaja, respectively (see Figs. 5c and d, 6b).Daytime deviations were very small because of the strategy adopted for searching the H 2 O time lag (see Appendix B).

Critical steps and recommendations for open-path systems
For open-path gas analysers the WPL correction is the most critical step, and although it depends on ecosystem type, season and target gas, the WPL correction terms can often surpass the magnitude of the flux itself and also change the sign of the measured target gas flux (e.g.Peltola et al., 2013).Thus it is critical to perform this correction accurately, especially when small target gas fluxes co-occur with large H and LE fluxes.Here, we have shown that LI-7500 F CO 2 values, measured at Erottaja urban site, are underestimated by 38 % on average if the correction is omitted, and the impact of the combined WPL and spectroscopic correction is even larger for LI-7700 F CH 4 measured at the Siikaneva fen site (Sect.3.3).A 2 % difference of daytime values of H , as calculated by the two software packages, was found in Siikaneva (data not shown).The difference stems from the fact that in EddyUH a site-specific co-spectral model is used in the spectral correction, whereas EddyPro uses the co-spectral model from Moncrieff et al. (1997) in doing the spectral correction to H . Nevertheless, the effect on WPL and spectroscopic correction terms for LI-7700 F CH 4 was negligible, as seen in Fig. 3.If the response time which characterises the measurement system ability to measure the flux contribution of small eddies, i.e. high frequencies, is determined using power spectra, as done in EddyPro (Fratini et al., 2012), then the highfrequency dampening caused by spatial sensor separation needs to be estimated separately.In EddyPro it was done using the method proposed by Horst and Lenschow (2009), while EddyUH uses cospectra to estimate the measurement system's high-frequency response, and thus no additional correction for sensor separation is needed.The Horst and Lenschow (2009) method is based on co-spectral peak frequency (n m ) parameterisations against the stability parameter ζ , in addition to the ASL co-spectral model (as presented in Horst, 1997).Using these assumptions, they derived a dependence between the signal dampening due to sensor separation and the co-spectral peak wavenumber and sensor separation in crosswind, along-wind and vertical directions.
For LI-7700 at the Siikaneva site, it was shown that this correction resulted in systematic differences between the software packages (Sect. 3.4) and between the two co-located CH 4 instruments (Sect.3.2).The correction method seemed to overcorrect LI-7700 CH 4 fluxes, which resulted in CH 4 fluxes that are too high.As mentioned above, the correction method relies on n m vs. ζ parameterisations and Horst (1997) co-spectral model and, if these do not comply with the spectral characteristics of turbulence observed at the site, then the correction will be biased.Furthermore, LI-7700 was situated significantly below the sonic anemometer (0.45 m) when compared with the sonic measurement height (2.75 m), and possibly, in such cases, the correction method does not perform well.Nevertheless, the difference observed in this study emphasises the need for accurate spectral corrections and the importance of minimising the sensor separation when constructing an EC measurement set-up.The sensor separation corrections are especially important for open-path analysers, since they cannot be mounted very close to the sonic anemometer due to their size and the flow distortion they may create.

Conclusions
We have estimated and analysed the flux uncertainty due to the use of two software packages, using datasets 2 and 3 months long including CO 2 , CH 4 and LE fluxes measured over a wetland and an urban site in Finland.Outputs from Ed-dyUH and EddyPro, two popular software packages for postfield processing of eddy covariance data, were compared.We evaluated the most critical processing steps, determining the largest deviations in the calculated fluxes.We found that the raw data preparation and processing steps were consistent between the software packages, and most of the deviations in the estimated fluxes were due to the flux corrections.Among the different calculation procedures analysed, the spectral correction was the most relevant for closed-path LE fluxes, reaching a night-time median value of 15 % at the wetland site.We found up to 43 % deviation (with respect the reference run) if the closed-path CO 2 flux is calculated without the dilution correction, while the CH 4 fluxes were up to 10 % lower without dilution and spectroscopic corrections.The WPL (and spectroscopic) correction was the most critical step for open-path systems.However, we also found large spectral correction factors for the open-path CH 4 fluxes, due to the sensor separation effect.Turbulent fluxes calculated with a reference combination of processing steps were in good agreement, the systematic difference between the two software packages being up to 2 and 6.7 % for half-hour and 2-month cumulative sum values, respectively.This result is an improvement with respect to earlier software intercomparison studies (e.g.Mauder et al., 2008), and it suggests that a consistent choice of implemented methods for the postfield processing steps can minimise the systematic flux uncertainty due to the usage of different software packages.Finally, it is recommended for future studies that the impacts of processing steps on fluxes are investigated in more detail, including validation of new methods and corrections across different types of compounds/instruments and ecosystems.

Data availability
Data are freely available upon request from the authors.

EddyUH
(https://www.atm.helsinki.fi/Eddy_Covariance/EddyUHsoftware.php) is written in MATLAB and includes a graphical user interface (GUI).In order to advance methodological issues (concerning especially CH 4 and N 2 O fluxes), besides standardised procedures, the most recent corrections and methods have also been implemented in EddyUH (see Table A1).
Post-field EC data processing with EddyUH is done through user-defined projects.A project in this context means a certain time period of data from a certain site that are processed with certain user-defined processing methods.These methods are determined by the user using the GUI and are saved in a set-up file where the site specifics and measurement system characteristics among other things are also defined.Therefore, all the processed data are always related to the saved project.The same project may include up to five different gas analysers combined with the same ultra-sonic anemometer, giving the possibility for the user to process several raw datasets at the same time.The software includes  a number of modules, which operate at different levels of post-processing (Fig. A1).Preliminary fluxes are calculated in the preprocessor, where the first level of processing is done to the raw dataset.Then several corrections are applied in the flux-calculation module, and the final fluxes are calculated (Table A1).In order to optimise the processing and properly apply all needed corrections, several software tools are available (Fig. A1).Co-spectral data are used in the high-frequency spectral transfer function estimator, where the low-pass filter time constant is experimentally estimated for each gas according to Mammarella et al. (2009).This approach is particularly relevant in the case of closed-path systems.Further, the time lag optimiser is a useful tool to verify the correctness of the chosen time lag window (and eventually refine it) for each gas, as well as to determine the varying window boundaries for H 2 O as explained in Appendix B. Finally, other modules are available in EddyUH, for example for analysing the estimated spectra and cospectra, for calculating flux uncertainties and for determining the flux footprint statistics according to the Kormann and Meixner (2001) model (Fig. A1).
Table A1.List of implemented methods for data processing in EddyUH.

Quality control and spike detection
Raw data quality tests (Vickers and Mahrt, 1997), diagnostic flags; several despiking methods (Rebmann et al., 2012) Conversion to dry mole fraction for closed-path gas analysers Dilution correction point by point Spectroscopic correction Rella et al. (2010) Calculation of turbulent fluctuations Block averaging, linear detrending and autoregressive running mean filter (Rebmann et al., 2012) Coordinate rotation Planar fit (Wilczak et al., 2001), sector-wise planar fit, 2-D rotation (Rebmann et al., 2012) Crosswind correction of sonic temperature Liu et al. (2001) Time lag determination Constant time lag, cross-covariance maximisation, time lag optimisation Quality statistics Flux steady-state and integral turbulence characteristics (Foken and Wichura, 1996), instrumental noise (Lenschow et al., 2000) and random flux error (Finkelstein and Sims, 2001) Corrections to the covariances Implemented methods

Stability parameter
Eq. (C1) High-frequency loss Theoretical method (Moncrieff et al., 1997); experimental method, e.g.empirical estimation of TF H and cospectra model (Mammarella et al., 2009) Low-frequency loss Rannik and Vesala (1999) Humidity The raw data are quality-flagged according to physical plausible ranges of high-frequency values of each variable, diagnostic parameters (if available) and several tests, as described in Vickers and Mahrt (1997).Further spikes are then detected, and commonly, this is done applying the Vickers and Mahrt (1997) method.However, other methods also exist in EddyUH, e.g. the difference limits method, which compare the difference between consecutive data points to a given threshold for each raw data time series (see Rebmann et al., 2012 for more details).If the time series contains too many spikes, the data might be useless and a flux should not be calculated for the averaging period of interest (commonly 30 min).Foken ( 2008) suggests excluding time series with more than 1 % spikes from further analysis.
B2 Point-by-point conversion to dry mole fraction for closed-path gas analysers (dilution correction) Current closed-path gas analysers measure H 2 O inside the sampling cell, making the conversion of gas mole fraction relative to dry air possible through a point-by-point dilution correction.
B3 Point-by-point spectroscopic correction for closed-path gas analysers In addition, H 2 O affects the shape and width of an absorption line used to estimate gas concentration via pressure broadening.This cross-interference can be corrected with the socalled spectroscopic correction (e.g.Rella et al., 2010).Many of the new laser-based gas analysers output dry mole fraction and thus no dilution or spectroscopic corrections are needed during data post-processing.However, for older gas analysers these corrections are needed.

B4 Coordinate rotation
A coordinate rotation is applied to the wind velocity components, in order to align the x axis parallel to the mean wind direction and to set the mean vertical velocity equal to zero.This is done according to common practice with two alternative approaches, the so-called 2-D rotation (Rebmann et al., 2012) or the planar-fit method (Wilczak et al., 2001).

B5 Calculation of turbulent fluctuations
In order to extract the turbulent fluctuations from the measured time series, the mean values are subtracted from the time series.There are three methods available in EddyUH, i.e. block averaging, linear detrending and autoregressive filtering (Rebmann et al., 2012).Of these three methods, only block averaging fulfils the Reynolds averaging rules.
All methods attenuate the low-frequency part of the cospectra.Block averaging has the smallest effect on the cospectra, whereas autoregressive filtering attenuates the cospectra the most (Rannik and Vesala, 1999).Linear detrending and autoregressive filtering are methods used to remove unwanted low-frequency variation (trend) in the signal (e.g.Mammarella et al., 2010).However, often block averaging is recommended.
B6 Crosswind correction of sonic temperature Sonic anemometers calculate sonic temperature T s based on three paths, and thus crosswind should be taken into account.
The correction can be applied point by point to the temperature fluctuations (Liu et al., 2001;Eq. 10) or to the temperature covariance (Liu et al., 2001;Eq. 12).Note that some sonic anemometers might include this correction in their internal firmware.

B7 Time lag determination and adjustment
The gas signal measured by closed-path analysers usually lags behind the wind speed measurement made with the sonic anemometer.The time lag can be estimated theoretically if sampling tube length and diameter are known, in addition to the flow rate in the tube.However for H 2 O the time lag depends on relative humidity due to adsorption and desorption of water on the tube walls (Ibrom et al., 2007a;Mammarella et al., 2009;Massman and Ibrom, 2008).The lag between open-path gas analyser measurement and sonic anemometer measurement is caused by the sensor separation: the further away the gas analyser is from the anemometer, the longer the time lag between the two measurements is.It also depends on wind speed and direction.(Foken and Wichura, 1996), instrumental noise (Lenschow et al., 2000) and random flux error (Finkelstein and Sims, 2001).All these data are saved in monthly binary files, and then used by other modules (Fig. A1).

Appendix C: Corrections to the covariances in EddyUH
In the second level of processing, several corrections must be applied to the 30 min covariances, and the set of corrections are different for closed-and open-path systems (see Fig. 1).At this stage the estimated covariances are used to calculate the stability parameter defined as

C1 Spectral correction
Flux loss at high frequency is due to the incapability of the measurement system to detect small-scale variation.The inadequate frequency response, sensor separation and line averaging, and, in closed-path systems, the air sampling trough tubes and filters are the main reasons causing co-spectral attenuation.On the other hand, flux loss at low frequency is due to limited averaging period (30 min) and trend removal.The frequency response correction is usually performed based on a priori knowledge of the system transfer function and the unattenuated cospectrum, e.g.
Here CF is the estimated spectral correction factor, C ws the normalised unattenuated cospectrum, f the frequency and TF = TF H • TF L the total transfer function.The correction, performed by multiplying the covariance by the factor CF, always increases the flux magnitude.The low-frequency correction depends on the method used for calculating the turbulent fluctuations (see Appendix B), and is performed using theoretically derived formulations for TF L (Rannik and Vesala, 1999).
The high-frequency transfer function TF H can be derived either theoretically or experimentally (Foken et al., 2012).The correction is different for momentum flux, sensible and latent heat fluxes and other gas fluxes, and it differs between open-and closed-path EC systems.In the theoretical approach, the ASL co-spectral models (Moncrieff et al., 1997) are used, and TF H is calculated as superposition of specific transfer functions representing different causes of flux loss, whose formulas can be found in Leuning and Judd (1996), Moncrieff et al. (1997) and Moore (1986).This approach works fine for correcting the momentum and sensible heat flux, as well as for gas fluxes measured by open-path systems.Alternatively, the experimental approach can be used, where the model cospectra and TF H are estimated using in situ measurements.Different methods have been proposed for retrieving the TF H from the measured power spectra or cospectra of T s and the target gas dry mole fraction (Fratini et al., 2012;Ibrom et al., 2007a;Mammarella et al., 2009;Nordbo et al., 2014).In EddyUH the method by Mammarella et al. (2009) is included.Many studies have used the theoretical approach because it is simpler to apply, while the experimental approach requires site-and sensor-specific investigations.

C2 WPL correction
For measurements done with an open-path gas analyser, fluctuations in air density cause apparent variations in measured scalar concentration, and this needs to be corrected according to Webb et al. (1980).The correction is performed to 30 min fluxes of the target gas, and the H 2 O and temperature covariances used in the correction should correspond to situations in ambient air, e.g.those that are fully corrected.For closed-path gas analyser the correction can be done to the covariances as well as an alternative of the point-by-point dilution correction applied to the raw data.The correction is less critical than for an open-path gas analyser due to the fact that the temperature fluctuations are usually dampened in the sampling line (Rannik et al., 1997).Then, only H 2 O fluctuations are relevant, and it is preferable that H 2 O is measured in the same cell as the gas whose flux is corrected.If this is not the case, then external H 2 O measurements can be used, provided that the H 2 O covariance is modified to correspond to circumstances in the measurement cell (Ibrom et al., 2007b;Peltola et al., 2014).

C3 Spectroscopic correction
Gas molar concentration measurements carried out with instruments based on laser spectroscopy (like LI-7700 and G1301-f analysers) also require corrections for spectroscopic effects that affect measured values, in addition to the abovementioned WPL or dilution corrections.As these spectroscopic effects are related to the changes in shape of the absorption line, due to the changes in gas temperature, H 2 O and pressure, one can incorporate the spectroscopic effects into WPL terms (modified equation) (McDermitt et al., 2011).
When estimating the CH 4 fluxes using the LI-7700, one should always use the temperature coming from the sonic anemometer and not those recorded by the in-path thermocouple of the LI-7700 instrument.Furthermore, one should always use the uncorrected CH 4 molar concentration (mmol m −3 ) for flux measurements, while H 2 O fluxes necessary to compute CH 4 fluxes should be acquired with an H 2 O analyser (in our study LI-7000).
For closed-path systems measuring H 2 O in the same sampling cell, it is simple to do the correction to the raw data point by point (see above).In case H 2 O is not measured internally, it is preferably to dry the air samples or to have external H 2 O measurements.In the latter case the correction can be applied to the half-hourly averaged fluxes according to the method proposed by Peltola et al. (2014).

C4 Humidity correction of sonic temperature
The correction is based on the transformation of sonic temperature (T s ) to actual air temperature (T).In EddyUH, the updated version (van Dijk et al., 2004) of the original Schotanus et al. (1983) correction is implemented.Following the derivation in van Dijk et al. (2004) the temperature covariance is calculated as w T = (1 − 0.51q) w T s − 0.51T w q , (C3) where w T s and w q are the final sonic temperature and H 2 O covariances (e.g. after spectral correction), and q is specific humidity (kg (H 2 O) kg (moist air) −1 ).The covariance w T is then used in Eq. ( 3) to calculate H, while the covariance w T s is used to recalculate the stability parameter in Eq. (C1).

C5 Iteration of corrections
Corrections to the covariances are repeated in an iteration loop until the flux change is smaller than 0.01 % (see Fig. 1).
In EddyUH these steps are performed in the "flux calculation" module, including the estimates of flux density according to Eqs. ( 1)-( 4).

Figure 1 .Figure 2 .
Figure1.EC data processing scheme for open-and closed-path gas analyser data.Relative magnitude of each processing step is also reported, according to this and other studies.Note the different ways with different sign conventions (see footnotes) in which the relative magnitudes are calculated in different studies, where F raw is the uncorrected flux, F run the flux without the corresponding correction, F corr the flux after the corresponding correction and F ref the fully corrected flux.The reported percentages demonstrate the typical impact that the correction methods have on the fluxes and should be considered only as indicative, since they depend on measurement site and set-up characteristics.

Figure 4 .
Figure 4. Scatter plots of F CH 4 measured by G1301-f and LI-7700 (a, b), LE measured by G1301-f and LI-7000 (c, d), F CO 2 measured by LI-7500 and LI-7200 (e, f) and LE measured by LI-7500 and LI-7200 (g, h).Subplots in the left column show fluxes calculated with EddyUH, and those in the right column fluxes calculated with EddyPro.Each dot represents a 30 min data value.Dashed lines indicate the 1 : 1 line, and red solid lines the linear regression to the data.

Figure 5 .Figure 6 .
Figure 5.Effect of different calculation procedure of the estimated flux at the Siikaneva site, presented as deviation (in %) from the reference run (ref).Deviation is defined as (run-ref)/ref, where run refers to the run performed with no spectral correction (no spec), theoretical spectral correction (theor spec), no WPL (or dilution) and spectroscopic correction (no WPL) and using a constant time lag (const lag).Bars indicate the median values, and error bars denote 25th and 75th percentiles.Note the different scale on y axis in the subplots (b, e) compared to (a, c, d).

Figure 7 .
Figure 7. Relative humidity (RH) dependence of low-pass filter time constant as estimated with the two software packages.Note the different scale on y axis in the subplot (c).
et al. (1980) for open-path gas analysers and Ibrom et al. (2007b) for closedpath gas analysers Spectroscopic correction Based on McDermitt et al. (2011) for open-path gas analysers and Peltola et al. (2014) for closed-path gas analysers B1 Quality control and despiking )where z is the measurement height (m), d the displacement height, (m), L the Obukhov length (m), T p the potential temperature (K), u * = −1 ), g the acceleration due to the gravity (m s −2 ) and κ the von Karman constant.
EddyUHsoftware.php, and it was originally developed in order to harmonise data processing among several EC sites operated by the group.Later, the software package was made publicly available.The EddyUH processing flow chart (Fig.A1) is introduced in Appendix A, while a brief description of post-field data processing operations and methods, as presented in TableA1, is given in Appendices B and C.The EddyUH software was compared against EddyPro, perhaps the most used software in the EC flux community, developed by LI-COR Biosciences Inc. (Lincoln, NE, USA).It is freely available and well documented (www.licor.com/eddypro).
based on Ed-dyUH version 1.7 and EddyPro version 5.2.1.EddyUH is a software package for EC raw data processing, developed by the Micrometeorology Research Group at the Department of Physics, University of Helsinki (Finland).It is freely downloadable from https://www.atm.helsinki.fi/Eddy_Covariance/

Table 1 .
Software set-ups for the reference combination.
al.: Quantifying the uncertainty of eddy covariance fluxes the preprocessor, statistics of the determined time lags can be evaluated with the time lag optimiser, in addition to the H 2 O time lag RH dependence.Later these results can be used in the final flux calculation step.For H 2 O a RH-dependent lag window can be used.For other gases if the determined lag deviates more than 3σ from the mean time lag (where the averaging period is defined by the user), then the mean time lag is used.This approach limits the possibility for erroneous time lags, which may occur if no clear maximum in the crosscovariance can be found.