InnFLUX – an open-source code for conventional and disjunct eddy covariance analysis of trace gas measurements: an urban test case

We describe and test a new versatile software tool for processing eddy covariance and disjunct eddy covariance flux data. We present an evaluation based on urban non-methane volatile organic compound (NMVOC) measurements using a proton transfer reaction quadrupole interface time-of-flight mass spectrometer (PTR-QiTOF-MS) at the Innsbruck Atmospheric Observatory. The code is based on MATLAB® and can be easily configured to process highfrequency, low-frequency and disjunct data. It can be applied to a wide range of analytical setups for NMVOC and other trace gas measurements, and is tailored towards the application of noisy data, where lag time corrections become challenging. Several corrections and quality control routines are implemented to obtain the most reliable results. The software is open source, so it can be extended and adjusted to specific purposes. We demonstrate the capabilities of the code based on a large urban dataset collected in Innsbruck, Austria, where three-dimensional winds and ambient concentrations of NMVOCs and auxiliary trace gases were sampled with high temporal resolution above an urban canopy. Concomitant measurements of 12C and 13C isotopic NMVOC fluxes allow testing algorithms used for determination of flux limits of detection (LOD) and lag time analysis. We use the high-frequency NMVOC dataset to generate a set of disjunct data and compare these results with the true eddy covariance method. The presented analysis allows testing the theory of disjunct eddy covariance (DEC) in an urban environment. Our findings confirm that the disjunct eddy covariance method can be a reliable tool, even in complex urban environments when fast sensors are not available, but that the increase in random error impedes the ability to detect small fluxes due to higher flux LODs.


S.1 Flow Chart
The workflow of the innFLUX code is illustrated in Figure S1, which depicts a flow chart of the individual processing steps.
The sequence of steps needing user input or feedback starts with the first step (step 0), which is only necessary when the directional planar fit option is used for anemometer tilt correction.
Step 0 is followed by step 1 and 2. Steps 0, 1 and 2 correspond to the MATLAB scripts innFLUX_step0, innFLUX_step1 and innFLUX_step2, respectively. Input, intermediary 5 and result output files are displayed in orange, and manual steps to be performed by the user are displayed in blue.

S.2 Comparison with EddyPro®
We performed a comparison between flux results obtained by the innFLUX code and EddyPro. Performing a rigorous 10 comparison of various eddy covariance software packages can be a complex task, as pointed out by Fratini and Mauder (2014) in their in-depth comparison of EddyPro and TK3. While such a detailed comparison would be out of the scope of the present publication, we show some basic evaluations between innFLUX and EddyPro. To perform a meaningful comparison we choose available parameters in EddyPro in a way to match the parameters used for innFLUX as closely as possible. For anemometer tilt correction, we used the double rotation method, because it should be implemented the same way in EddyPro as in innFLUX. We did not apply any automatic spectral corrections. Results were filtered for high quality fluxes, using only fluxes of quality class 3 or below on the 1-9 scale introduced by Foken et al., 2004. Figure S2 shows a scatter plot of the CO2 flux calculated by innFLUX against the flux calculated by EddyPro. Figure S3 shows a corresponding plot for the 5 sensible heat flux. Overall innFLUX and EddyPro compare will with a slope of 1.02 (0.95) and an r 2 of 0.97 (0.99) for CO2 (sensible heat) fluxes.

S.3 Calibration setup used for the PTR-QiTOF-MS
The PTR-QiTOF-MS was calibrated using a calibration gas standard CC447665 from Apel-Riemer Environmental, Inc. The standard was used for sensitivity calibrations of the PTR-QiTOF-MS. The standard contained mixing ratios of individual compounds between 1 and 5 ppmv. The uncertainty of mixing ratios of the pure compounds was taken as 5%. A continuous, controlled flow (Bronkhorst Mass Flowcontroller F-200CV-002, 5 sccm, Kalrez plunger) was dynamically diluted into VOC 15 free air (catalytically scrubbed, Parker Zero Air Generator, model 1001) to produce calibration mixtures in the low double digit ppbv range. Table S1 lists the compounds of the calibration gas used in this study.

S.4 Spectral Analysis and Corrections
Co-spectra of scalars and vertical wind speed are a standard output of innFLUX and are calculated by means of fast Fourier transform (FFT). Co-spectra normalized by the integral over the entire frequency range (equals the covariance), frequency 5 scaled and bin-averaged into 60 log-spaced intervals over 6 decades of non-dimensional frequency are stored for each 30-minute period. Only co-spectra of half-hour periods were considered for that spectral analysis that a) pass at a quality class (Foken et al. 2012b) of three or less, b) pass the u* filtering and c) show a flux signal-to-noise ratio (SNR) of three or more to exclude spectra from periods where the covariance is very noisy. These scaled, non-dimensional 10 co-spectra can be averaged for each scalar (e.g. for toluene measured by PTR-Qi-TOF at 93.06988 m/z). An analytical , ≡ 3 4 ⁄ (S4.2) (see Lee et al. 2010 and references therein) was fitted to the mean co-spectrum to define a model spectrum for each scalar. Figure S4 shows the mean cospectrum of toluene (open circles) and the fitted model cospectrum (red). Fitting the model 15 cospectrum was weighted strongly towards a two decade wide window around the cospectral maximum (black line in Figure   S4).

Figure S4. Model spectrum (red line) in Equ. S4.2 was fitted to the average of normalized, frequency scaled toluene cospectra that passed quality checks a-c (N=800) (open circles). The fit was weighted around the cospectral maximum (black line).
For each half hour, a model spectrum ( , ̅ ) is derived from the non-dimensional model spectrum in dimensional 5 frequency based on the mean wind speed and the measurement height above the displacement plane (see red line in bottom panel of Figure 7 in the main text).
Transfer functions (see top panel of Figure 7 in the main text) for sensor separation, path averaging of the sonic anemometer, path averaging in PTR-MS drift tube, and gas inlet tube attenuation were calculated for each half-hour as described by with being the length of the drift tube and the flow velocity of the sample gas (not the ions!) along the drift tube.
Note that depends on the gas flow rate through the drift tube, drift pressure and temperature and the cross section area and therefore needs to be determined for the specific geometry and settings parameters of the respective PTR-MS instrument. Essentially, over is the characteristic exchange time of the drift tube gas volume. 5 Tube attenuation: If the tube flow is turbulent, the tube attenuation is described by with the geometry parameters and being inner tube radius and tube length, being the volumetric tube flow, and being the Reynold's number, which is a function of , and kinetic viscosity of the sample air. Note that one needs to make sure that is constant and/or measured. Inlet sections with different geometries or flows are treated as separate transfer 10 function. For the VOC flux measurement in the presented setup the inlet manifold transporting sample gas from near the anemometer to the laboratory and the subsampling of the PTR-Qi-TOF through a PEEK capillary were treated separately, but Figure 7 shows that the attenuation in the capillary is negligible. The multiplication of the model co-spectra with the product of the transfer functions results in attenuated co-spectra. The relative error on the flux introduced by high frequency attenuation (low-pass filtering) is given by  Foken et al. (2012a). Note that the high 20 frequency loss depends on operational settings, geometry parameters as well as the mean wind speed and direction of each flux averaging period. The measured co-variance of each scalar of each flux averaging period is corrected for the considered high frequency losses by dividing it by the ratio of the integrals in equation S4.8. The mean high frequency loss can be determined from the slope of a linear regression of the uncorrected versus corrected flux ( Figure S5). The geometry and operational parameters for the PTR-Qi-TOF flux measurement setup is listed in Table S2. 5