# CO2-gradient measurements using a parallel multi-analyzer setup

L. Siebicke, G. Steinfeld, T. Foken

► To cite this version:

L. Siebicke, G. Steinfeld, T. Foken. CO2-gradient measurements using a parallel multi-analyzer setup. Atmospheric Measurement Techniques, 2011, 4 (3), pp.409 - 423. 10.5194/amt-4-409-2011. hal-02652083

HAL Id: hal-02652083

https://hal.inrae.fr/hal-02652083

Submitted on 29 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CO₂-gradient measurements using a parallel multi-analyzer setup

L. Siebicke<fnref n="1" />,*, G. Steinfeld<fnref n="2" />, and T. Foken<fnref n="1" />

Received: 17 August 2010 – Published in Atmos. Meas. Tech. Discuss.: 11 October 2010

Revised: 1 February 2011 – Accepted: 2 February 2011 – Published: 1 March 2011

**Abstract.** Accurate CO₂ concentration gradient measurements are needed for the computation of advective flux terms, which are part of the full Net Ecosystem Exchange (NEE) budget equation. A typical draw back of current gradient measurement designs in advection research is the inadequate sampling of complex flow phenomena using too few observation points in space and time. To overcome this draw back, a new measurement design is presented which allows the parallel measurement of several sampling points at a high frequency. Due to the multi-analyzer nature of the design, inter-instrument bias becomes more of a concern compared to conventional setups. Therefore a statistical approach is presented which allows for accurate observations of concentration gradients, which are typically small in relation to analyzer accuracy, to be obtained. This bias correction approach applies a conditional, time dependent signal correction. The correction depends on a mixing index based on cross correlation analysis, which characterizes the degree of mixing of the atmosphere between individual sample points. The approach assumes statistical properties of probability density functions (pdf) of concentration differences between a sample point and the field average which are common to the pdf's from several sample points. The applicability of the assumptions made was tested by Large Eddy Simulation (LES) using the model PALM and could be verified for a test case of well mixed conditions. The study presents concentration time series before and after correction, measured at a 2 m height in the sub-canopy at the FLUXNET spruce forest site Waldstein-Weidenbrunnen (DE-Bay), analyzes the dependence of statistical parameters of pdf's from atmospheric parameters such as stratification, quantifies the errors and evaluates the performance of the bias correction approach.

The improvements that are achieved by applying the bias correction approach are one order of magnitude larger than possible errors associated with it, which is a strong incentive to use the correction approach. In conclusion, the presented bias correction approach is well suited for – but not limited to – horizontal gradient measurements in a multi-analyzer setup, which would not have been reliable without this approach. Finally, possible future improvements of the bias correction approach are outlined and further fields of application indicated.

## 1 Introduction

Advection is a part of Net Ecosystem Exchange (NEE) of carbon dioxide, the determination of the latter being a primary focus of a world wide network of vegetation-atmosphere exchange measuring stations, the FLUXNET (Baldocchi et al., 2001). Not only are reliable measurements of advection lacking for most FLUXNET sites, but they continue to be a challenge even for specialized advection research experiments (e.g. Aubinet et al., 2003; Staebler and Fitzjarrald, 2004; Feigenwinter et al., 2008; Aubinet et al., 2010). Advection remains further to be a major reason for the night flux problem (Finnigan, 2008). Mathematically, scalar advection is the product of the mean spatial gradient of a scalar – CO₂ in the case of the current study – and the mean wind velocity, i.e. scalar transport with the mean flow. Advection is typically addressed as vertical advection (Lee, 1998; Baldocchi et al., 2000) and horizontal advection (Baldocchi et al., 2000; Aubinet et al., 2003).

There are two main conceptually different reasons why valid and representative advection measurements are difficult to obtain. One is the instrument related accuracy, with which scalar gradients and wind vectors of the mean flow

<figure><img src="image_2.png" /></figure>

Correspondence to: L. Siebicke
(lukas.siebicke@ecofog.gf)

<footnote n="1">¹Department of Micrometeorology, University of Bayreuth, Bayreuth, Germany</footnote>

<footnote n="2">²Institute of Physics, ForWind, Center for Wind Energy Research, Carl von Ossietzky University of Oldenburg,
Oldenburg, Germany</footnote>

<footnote>*now at: Institut national de la recherche agronomique (INRA), Kourou, French Guiana</footnote>
can be measured. The other reason being undersampling of complex flow phenomena due to limited resources of real world experiments, thus yielding measurements which are not representative for a spatial (volume) and temporal (time period) mean but for a point only.

Vertical and horizontal advection pose different measurement challenges. With regards to vertical advection, reliable vertical CO₂ concentration gradients can be obtained due to vertical concentration gradients which are relatively large compared to sampling uncertainties. Measurements of vertical wind velocity are difficult to obtain, both for reasons of accuracy, precision, and resolution of sonic anemometers and particularly for reasons of the limited spatial representativeness of a point measurement. Spatially representative measurements of vertical wind speed can never be obtained from a single point measurement in complex flow, due to theoretical reasons; therefore multi-tower measurements – possibly in combination with airborne measurements – are being suggested to improve spatial representativeness of vertical wind measurements (e.g. Mahrt, 2010). Alternatively, the vertical wind velocity measurement problem is avoided by using a mass continuity approach, i.e. inferring vertical motion from horizontal divergence (e.g. Vickers and Mahrt, 2006; Montagnani et al., 2010) or a combination of the mass continuity approach and modeling (Canepa et al., 2010). Regarding horizontal advection, measurements of horizontal wind speed can be obtained with sufficiently high accuracy with sonic anemometers, even though they are often not spatially representative. In contrast, horizontal gradients are very difficult to measure with the required accuracy, because mean gradients are small in relation to instrument related uncertainty and difficult to measure at a large enough number of locations with a sufficiently high temporal resolution.

It is the main aim of this study to provide improvements for the measurement of horizontal CO₂ concentration gradients by means of a better temporal and potentially better spatial resolution. An improved resolution is needed for advection measurements in heterogeneous forests as could be shown by analyzing the effects of spatial heterogeneity and short lived phenomena on mean horizontal CO₂ concentration gradients at the site under study.

The most commonly used setup for horizontal gradient measurements is based on a switching valve system (e.g. Burns et al., 2009), which uses a single closed-path infrared gas analyzer to sample several points one after the other (“sequential approach”), returning to the same sample point once every few minutes. There is an inherent tradeoff between achievable spatial and temporal resolution. The main benefit of this setup is a common analyzer for a number of sample locations, reducing the risk of bias between those points. The current study utilizes a multi-analyzer setup, featuring an individual closed-path infrared gas analyzer for every measurement point, enabling simultaneous measurements of all points (“parallel approach”) with a high frequency. Temporal resolution is no longer parasitic to spatial resolution, the

latter depending on available resources only. With ten individual analyzers used, the spatial resolution is on the order of a sequential system. Thus the system described is capable of making measurements which are representative in the temporal domain since it can observe all relevant temporal scales of the CO₂ concentration signal.

Valid concentration measurements need to be both precise and accurate. Precision of the parallel approach used in this study is higher compared to the conventional sequential approach because there are potentially much more values available in one averaging interval, thus reducing random error. The advance in the number of values is proportional to the number of sample locations per analyzer for the sequential approach. Lower accuracy of a multi-analyzer setup compared to a single analyzer setup due to inter-instrument bias is the major drawback of the parallel approach, in addition to higher resource requirements. Bias can be reduced by careful system design and frequent calibration against accurate, known standards. Section 2.2 lists technical measures that have been taken to that end for the presented system. How to deal with the remaining bias will be the topic of the rest of the paper. The basic assumption regarding concentration differences originating from natural gradients stated in Sect. 2.4.2, which is the justification of the proposed bias correction approach, has been implicitly used by Aubinet et al. (2003) and applied for time series correction in a simple, time independent manner whereas the current study applies a conditional, time dependent signal correction. Previous studies using more than one closed path gas analyzer in a multiplexer system with multiple sampling inlets have often used co-located inlets to deal with time dependent inter-instrument bias (e.g. Sun et al., 2007), and the same procedure was applied to vertical profile measurements at the site of the current study. However, due to the characteristics of the multi-analyzer system presented in this study with only one inlet per analyzer, co-located inlets cannot be used in the same way and a new approach is needed. A number of options for inter-instrument comparison using direct measurements, which combine the setup described in the present study with the concept of co-located inlets are discussed in Siebicke (2011) in order to aid independent evaluation of the statistical calibration method presented here.

It should be noted that the term “CO₂ concentration” is used throughout this paper to describe basic principles in a consistent way. It specifically refers to “molar fraction” in units of mol·mol⁻¹ or µmol·mol⁻¹, which were used for all measured values presented herein, whereas it refers to “CO₂ density” in units of kg·m⁻³ for modelled values from the Large Eddy Simulation study (Sects. 2.5 and 3). However, further applications of the ideas about bias correction presented in this paper may prefer to describe CO₂ in terms of “mixing ratio” in units of kg·kg⁻¹ (Kowalski and Serrano-Ortiz, 2007).
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Sampling locations for sub-canopy CO₂ concentration at a 2.25 m height. M-numbers are used for reference in the text. Topography is shown by isolines with an equidistance of 0.2 m relative to 750 m a.s.l.</figcaption></figure>

# 2 Material and methods

## 2.1 Site

Measurements were carried out at the FLUXNET site Waldstein-Weidenbrunnen (DE-Bay), 50°08'31" N, 11°52'01" E, a hill site in the Fichtelgebirge Mountains in Southern Germany. The Norway spruce (Picea abies) stand is on the upper section of a forested hill, 775 m ASL, with a 3° slope facing south-west. The tree height within the footprint of the measurements is 25 m. The site is described in detail in Gerstberger et al. (2004) and a summary of background data can be found in Staudt and Foken (2007).

## 2.2 Instrumental setup

Wind vector and CO₂ concentration time series were recorded along horizontal transects at a 2.25 m height in the sub-canopy space. The spatial setup of sub-canopy sample locations is shown in Fig. 1. Ten CO₂ concentration sample points were distributed between an along slope transect from north-east to south-west (5 sample points) and an across slope transect from north-west to south-east (6 sample points), including one common point. Each point was sampled by an individual closed-path infrared gas analyzer. Instruments used were five LI-6262, one LI-6251 (LI-COR Biosciences Inc.), four BINOS (Leybold Heraeus GmbH). In addition to CO₂ concentration measurements at a 2.25 m height, sample locations M5, M6, M7, M8, M9, M10

(see Fig. 1) were equipped with sonic anemometers (USA-1, METEK GmbH) to measure wind speed, wind direction and sonic temperature at the same height. CO₂ concentration measurements are available with a frequency of 1 Hz at each sample point, sonic data were recorded at a 20 Hz frequency. To reduce the risk of systematic differences between individual closed-path gas analyzers the system was carefully designed to avoid any possible bias of the concentration measurement from differences in pressure or temperature (sample air temperature, ambient analyzer temperature, radiation). All CO₂ closed-path gas analyzers shared a common housing in a central position with controlled conditions resulting in a constant common temperature and common pressure regime. Moreover, all analyzers shared a common tailor-made automatic calibration system, using high precision reference gases (accuracy 0.1 µmol mol⁻¹). The calibration routine included an automatic calibration every 4 h using two reference concentrations, which were sampled by all ten analyzers at the same time. In addition to factory calibration, each instrument's polynomial calibration function was established on site, using multiple standards. The polynomial was checked before and during the experiment.

Individual technical measures taken to avoid systematic inter-instrument bias include the following:

- The length of tubing connecting each sample point with the corresponding gas analyzer was exactly 75 m for every point. Sample tubes used were of polyethylene-aluminum composite structure, model DEKABON 1300-M060X (Serto AG, Fuldabrück, Germany) with an inner diameter of 4 mm.

- Large diameter line intake air filters were checked regularly and replaced synchronously at all points, if necessary.

- Common ambient temperature and pressure for all gas analyzers and calibration unit, including radiation protection, active automatic temperature control by heating and cooling as well as carefully designed ambient air circulation.

- Quality control of performance of automatic temperature control system, making sure that ambient air temperatures measured at several points surrounding the gas analyzers remain within acceptable range.

- Temperature adaptation for sample lines, to allow the temperature of sample air in all sample lines to equilibrate to a common temperature prior to entering the analyzer.

- Common temperature and radiation shielding for reference gases.

- Minimization of dead volumes in calibration and valve system to ensure turbulent flow conditions and avoid contamination by previous samples.
- Flow rate of 2 L min⁻¹ (Reynolds number Re = 2520) above critical flow rate of 1.8 L min⁻¹ at critical Reynolds number (Re crit = 2320) to ensure turbulent flow conditions in all tubes, at the same time keeping the flow rate as low as possible to minimize pressure drop across the system.

- Regular flow rate check and adjustment for all sample lines.

- Bypass system to avoid back pressure effects during calibration, featuring a low pressure drop bypass flow rate control device to ensure minimum necessary bypass flow and avoid possible reverse flow and sample contamination by ambient air.

- One common pump downstream of the analyzers to reduce effects of the pump on the concentration signals and to guarantee common pressure for all analyzers, assuming equal pipe geometry of all sample lines.

- Automatic control of constant overall system flow rate by mass flow controller.

- Passive system to allow for pressure equilibration between sample cells of individual gas analyzers by connecting all analyzer outlets to a manifold with a sufficiently large diameter.

- Pre-assembly measurement and evaluation of the pressure drop caused by individual system components to ensure that associated errors of the CO₂ concentration measurements are below accepted threshold.

- Vacuum and over pressure assisted leak check for the complete system to rule out sample contamination by ambient air.

## 2.3 Data set

The data set was collected during the second intensive observation period (IOP2), 1 June to 15 July 2008 of the EGER (“Exchange processes in mountainous Regions”) experiment (Serafimovich et al., 2008). 24.6 days worth of data were used for the analysis, i.e. 1181 half hourly values taken from a window of 32.0 days (11 June to 13 July). Periods were excluded from the analysis when instruments were powered off or obviously malfunctioning.

## 2.4 Theoretical considerations regarding concentration differences

### 2.4.1 Bias

An observed concentration difference between two spatially separated sample points is the sum of a concentration difference originating from a natural atmospheric concentration

gradient and the inter instrument bias, the latter being determined by systematic (bias) and random error of the individual instruments. We will refer to this composite concentration difference also as a concentration offset, Δc. While random error of the instruments is a minor concern in the current study due to sufficiently long averaging period, instrument bias can be reduced by calibration against known standards. The calibration procedure used in this study was outlined in Sect. 2.2. The remaining bias is the sum of the error of the calibration plus the instrument drift between two consecutive calibration events. This remaining bias cannot be removed by calibration since it is intrinsic to the calibration procedure itself. However, a statistical approach detailed in Sect. 2.7 can help to distinguish between remaining bias and concentration differences originating from natural gradients based on the observed signal.

### 2.4.2 Natural concentration differences

To separate concentration differences originating from natural gradients between two spatially disjunct (i.e. up to a few tens of meters) sample points from instrument bias the following assumption is made and is the basis for bias correction used in the current study: for certain meteorological conditions the concentration time series observed simultaneously at the two locations can be statistically linked to a reference concentration which is common to both sample locations. To be more precise, under the condition of well mixed, i.e. sufficiently turbulent atmospheric conditions (hereafter “mixed” conditions) the concentration difference between the two locations which is most likely to be observed is zero. If this statement is true for the concentration difference between any two points, it can also be applied to the difference between the concentration at one sample location cᵢ, and the spatial average concentration of the sample point field ȳc(t) at a given time t. ȳc(t), which serves as a reference concentration, describes the background concentration of the sample point field at time t using the median field concentration according to Eq. (1):

$$ \tilde{c} = \begin{cases} \frac{c_{\frac{k+1}{2}}}{\frac{1}{2}(c_{\frac{k}{2}} + c_{\frac{k}{2}+1})} & k \text{ odd} \\ \frac{1}{2}(c_{\frac{k}{2}} + c_{\frac{k}{2}+1}) & k \text{ even} \end{cases} \quad (1) $$

with k = 1...n observations (c₁, c₂, ..., cₖ) being the concentration measurements (c₁(t), c₂(t), ..., cₙ(t)) at n locations sorted in ascending order. The statistical measure describing the concentration difference most likely to be observed is the mode of the probability density distribution (pdf) of the concentration differences cᵢ(t) – ȳc(t), which is assumed to be close to zero under the condition of well mixed i.e. sufficiently turbulent atmospheric conditions.

This is illustrated in Fig. 2b for two hypothetical time series c₁(t) = 7, 6, 5, 5, 8, 5, 4, 6, 5, 6 and c₂(t) = 7, 6, 7, 5, 3, 5, 4, 5, 6, 5, displayed in Fig. 2a. The characteristics of turbulence justify the assumed mode of the pdf to be
<figure><img src="image_4.png" /></figure>

**Fig. 2.** Hypothetical concentration time series $c_1(t)$ and $c_2(t)$ with time $t \in [1, 10]$ (**a**, **c**), and corresponding frequency and density distributions of concentration differences $c_i(t) - \tilde{c}(t)$ (**b**, **d**) for mixed conditions (**a**, **b**) and for non mixed conditions (**c**, **d**). Regarding the density distributions in Fig. b and d, the histogram bars indicate the frequency for binwidths of 1.0, the solid line is a kernel density estimation generated with the same tools which were used for density estimation of measured concentration data as described in Sect. 2.7.

close to zero, i.e. turbulence consists of temporal perturbations of a mean state which are stochastic and relatively short in duration compared to the observation period. The mode is zero even though the time series $c_1(t)$ and $c_2(t)$ given in Fig. 2a have a different mean (temporal mean indicated by overline): $\overline{c_1(t)} = 5.7$ and $\overline{c_2(t)} = 5.3$, and even though the mean of the concentration difference $c_i(t) - \tilde{c}(t)$ is different from zero: $c_1(t) - \tilde{c}(t) = 0.2$ and $c_2(t) - \tilde{c}(t) = -0.2$.

For atmospheric conditions without turbulent mixing (hereafter “non mixed” conditions) above stated assumption does not need to be fullfilled. Since there is no effective mechanism of mixing, two sample locations can be continuously exposed to air masses with different concentrations – see concentration time series $c_1(t) = 4, 3, 2, 2, 5, 2, 1, 3, 2, 3$ and $c_2(t) = 8, 7, 8, 6, 4, 6, 5, 6, 7, 6$ in Fig. 2c – i.e. there is a persistent natural gradient and no common background concentration is observed at both sample points. Thus, the two points will most frequently sample a concentration difference which represents this gradient, and the mode of the probability density distribution is non zero, Fig. 2d.

All combinations of the well mixed and non mixed case are possible. It depends on turbulence statistics and the

length of the time series incorporated in the probability density distribution whether mixing is sufficient to produce a mode of the pdf close to zero or not. A method to quantify the degree of mixing will be presented in Sect. 2.6.

## 2.5 Large Eddy Simulation

An idealized Large Eddy Simulation (LES) case study was performed in order to check whether the assumption made in Sect. 2.4.2 is true, i.e. whether the mode of the probability density function of the difference $c_i(t) - \tilde{c}(t)$ between the scalar concentration at one sample location $c_i(t)$ and the scalar concentration averaged over all sample points $\tilde{c}(t)$ is close to zero for well mixed conditions, even in the case that the distribution of sources and sinks of the scalar is not homogeneous and a mean spatial concentration gradient $\frac{\partial \tilde{c}}{\partial y}$ with concentration $c$ and horizontal distance $y$ exists. Large Eddy Simulation is a tool that is used to study turbulence related processes in the atmospheric boundary layer. It can therefore be used to extract statistical properties of turbulence for the well mixed case. The simulation does not intend to perfectly mimic subcanopy conditions but to test general statistical properties of turbulent mixing, i.e. whether strong turbulent mixing is able to allow the average field background concentration $\tilde{c}(t)$ to emerge as the dominant mode of the probability density function rather than local sources or sinks producing the dominant mode.

The LES model used in this study is the Parallelised LES Model (PALM) that has been developed at the Institute of Meteorology and Climatology at the Leibniz University in Hannover, Germany. Detailed information on the LES approach, model equations and numerical schemes applied in PALM are given in Raasch and Schröter (2001) or – continuously updated – on-line on the homepage of the PALM group (Raasch, 2010). In our applications of PALM presented here an additional prognostic equation for a scalar quantity was solved so that the temporal development of a scalar concentration field with distributed sources and sinks could be simulated.

Two simulations with different setups were carried out for this study. In our first simulation (“case A”) a horizontally homogeneous distribution of scalar *sources* and *sinks* was prescribed. However, the scalar concentration *field* was initialized with a horizontal concentration gradient. This setup resulted in a temporally decaying horizontal concentration gradient due to turbulent mixing.

In the second simulation (“case B”) the initial *field* of scalar concentration and the *sinks* of scalar concentration were prescribed to be horizontally homogeneous. However, the *sources* of scalar concentration were horizontally heterogeneously distributed. This setup resulted in a temporally evolving concentration gradient.

In the following paragraphs the setup of the two simulation runs will be described in detail. First of all those settings
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Setup of Large Eddy Simulation study. Virtual sensor locations (a), Source-sink distribution (b) and background concentration gradient (c). Grid spacing: 5 m.</figcaption></figure>

common to both simulations will be reported before pointing out details regarding the scalar concentration field which are specific to each setup.

In both simulations the model domain consisted of 640 × 640 × 256 grid points and a basic grid spacing of 5 m was used. Above a height of 1000 m the grid was stretched vertically until a maximum grid size of 20 m was reached. The total extension of the model domain was 3.2 × 3.2 × 2 km.

Both LES simulations were initialized with wind profiles that were obtained from a one-dimensional prerun in order to accelerate the transition to a stationary state in the three-dimensional main run. The geostrophic wind ($u_g, v_g$) was prescribed as (3 m s⁻¹, 0 m s⁻¹) while $u$ and $v$ correspond to the x- and y-direction, respectively. The roughness length was 0.1 m. At the bottom boundary of the model domain a near-surface heat flux of 0.01 K m s⁻¹ was prescribed, so that a convective boundary layer with an Obukhov-length in the range between −40 and −50 m developed with time. The Coriolis parameter corresponded to a geographical latitude of 55°.

Sources and sinks of the scalar were switched on as soon as the simulation had reached a quasi-stationary state, i.e. after a spin-up time of 2 h. The sources of the scalar were situated at a height of 2.5 m and distributed homogeneously over the total horizontal extension of the model domain. The sinks of the scalar were also distributed over the total horizontal extension of the model domain but situated at a height of 27.5 m.

In both simulations time series of scalar concentration were recorded at 16 observation points within the x y cross section of the model domain at a height of 17.5 m beginning from the first release of scalar quantity until the end of the LES 7200 s later. Data from these time series could be used in order to calculate the differences between the concentration at a single sample point and the concentration averaged over all sample points as required in order to check the va-

lidity of the assumption made in Sect. 2.4.2. Figure 3 shows the locations of these observation points of the two LES. The coordinates of the 16 observation points were composed out of the x-coordinates (760 m, 785 m, 810 m, 835 m) and y-coordinates (760 m, 785 m, 810 m, 835 m). Thus, the distance between two adjacent observation points along the x- or y-direction was 25 m.

In case A the initial scalar concentration field showed a gradient along the y-direction. The initial concentration increased by $3.038 \times 10^{-7}$ kg m⁻⁴ from $y=0$ to $y = \frac{L_y}{2}$, while it decreased by $3.038 \times 10^{-7}$ kg m⁻⁴ from $y = \frac{L_y}{2}$ to $y = L_y$ ($L_y$ is the length of the model domain along the y-direction). It is worth mentioning that the prescribed gradients are equivalent to $\pm0.16\mu mol \ mol^{-1} \ m^{-1}$ which deliberately has been chosen to represent the maximum of gradients observed in the field at the site under study and published for other sites (Aubinet et al., 2003; Heinesch et al., 2007) during stable stratification, even though gradients are smaller during neutral and unstable stratification, i.e. the stratification regime present in the LES. In that sense, the LES with strong gradients tests a worst case scenario.

As in case B the initial mean scalar concentration prior to imposing the additional spatial gradients in case A was $6.997 \times 10^{-4}$ kg m⁻³. Note that this is equivalent to $378 \mu mol \ mol^{-1} \ CO_2$, which was the background concentration observed at the experimental test site described in Sect. 2.1. The resulting initial concentration field is shown in Fig. 3c. The source strength was set to $8.8 \times 10^{-8}$ kg m⁻³ s⁻¹, while the sink had a strength of $-8.8 \times 10^{-8}$ kg m⁻³ s⁻¹. It was chosen to correspond to a typical maximum daytime Net Ecosystem Exchange of $-20 \mu mol \ m^{-2} \ s^{-1}$ observed at the measurement site.

In case B (horizontally homogeneous initial concentration field) a basic source strength of $4.4 \times 10^{-8}$ kg m⁻³ s⁻¹ was prescribed at $y=0$ and $y = L_y$. Between $y=0$ and $y = \frac{L_y}{2}$ the gradient of the source strength, $\frac{\partial s}{\partial y}$,
was $5.5 \times 10^{-11} \text{ kg m}^{-4} \text{ s}^{-1}$ for $y < \frac{L_y}{2}$, while it was $-5.5 \times 10^{-11} \text{ kg m}^{-4} \text{ s}^{-1}$ between $y = \frac{L_y}{2}$ and $y = L_y$. Thus, the mean horizontal source strength was exactly the same as in case A. The sink strength was prescribed as in case A (and thus again approx. equivalent to a typical maximum daytime Net Ecosystem Exchange of $-20 \text{ µmol m}^{-2} \text{ s}^{-1}$ observed at the site).

It is obvious that the Large Eddy Simulations presented here are an idealization and do not account for the complexity of the given forest site, particularly because they do not fully account for the forest canopy. However, we would like to stress that the purpose of the simulation is to test the idealized case of turbulent mixing given realistic physical values of scalar concentration gradients and a vertical source and sink distribution that does mimic sources at the forest floor and sinks in the forest canopy with respect to their vertical distribution and their intensity. Verifying and accepting the assumption made in Sect. 2.4.2 first for an idealized case is necessary before addressing measurements from the more complex forest setting. Whether conditions in the forest at any given time show sufficient mixing is not evaluated by LES but by the application of an empirical mixing index (see Sect. 2.6) which is based on measured data.

## 2.6 Mixing index

A “mixing index” MI was formulated to quantify the degree of mixing between the real world sample points given in Fig. 1. A threshold value $MI_c$ was then used to separate well mixed conditions from not sufficiently mixed conditions. The mixing index MI was based on cross correlation $R(\tau)$, e.g. $R_{c_1c_2}(\tau)$ of the simultaneous concentration time series $c_1(t)$ and $c_2(t)$ at spatially separated sample locations, normalized by their mean variance $\sigma^2$. The cross correlation function is given as:

$$R_{c_1c_2}(\tau) = \frac{1}{T_F} \int_{-T_F/2}^{T_F/2} c_1(t) \cdot c_2(t+\tau) dt \quad (2)$$

with time lag $\tau$ between concentration time series $c_1(t)$ and $c_2(t)$, $T_F$ being the length of the time window of $c_1(t)$ and $c_2(t)$ and $\tau \in [-T_F, T_F]$. MI then writes:

$$\text{MI} = \max(|R_{c_1c_2}(\tau)|) \cdot \left( \frac{\sigma_{c_1}^2 + \sigma_{c_2}^2}{2} \right)^{-1} \quad (3)$$

More specifically, MI was calculated using the mean cross correlation of CO₂ concentration time series $c_5$ and $c_6$ recorded at a sample point pair oriented along the terrain slope (locations M5, M6) and $c_5$ and $c_8$ recorded at a sample point pair oriented across the slope (M5, M8) divided by the mean field variance of all concentration time series $c_5, c_6, \dots, c_{14}$ at sample locations M5, M6, \dots, M14 using a window length of $T_F = 60$ min.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Density distribution of mixing index MI (solid line). Dashed lines at $MI=0.06$ and $MI=0.12$ enclose range for sensible choices of a critical mixing index $MI_c$, (a). Diurnal course of mixing index on 29 June 2008 (solid line) and $MI_c$ (dashed line), (b). MI is representative for the whole sample point field (see Sect. 2.6 for details of the calculation).</figcaption></figure>

The critical mixing index $MI_c$ was empirically inferred from the density distribution of MI given in Fig. 4a. Sensible values were found to be in the range $MI_c \in [0.06, 0.12]$, corresponding to a sharp bend in the density distribution separating MI's representative of well mixed daytime conditions (distribution tail to the right of $MI_c$ in Fig. 4a) from low MI's representative of night time conditions with little mixing (distribution peak to the left of $MI_c$ in Fig. 4a). Figure 4b presents a typical diurnal cycle of the mixing index which is clearly separated into mixed and non mixed conditions by $MI_c$.

## 2.7 Bias correction

A statistical correction was applied to the CO₂ concentration time series $c_i(t)$ of every individual analyzer (which had previously been calibrated against known reference gas standards) to correct for remaining instrument related bias $\Delta c_i$. This yields the statistically bias corrected time series $c_{i,corr}(t)$ according to Eq. (4):

$$c_{i,corr}(t) = c_i(t) - \Delta c_i \quad \text{for } MI \ge MI_c \quad (4)$$

Instrument related bias of the CO₂ concentration signal was observed to vary over time. It is therefore appropriate to apply a bias correction that is time dependent, too. Analyzer specific values of instrument bias $\Delta c_i$ were computed for every 60-minute interval $T_F$ of the concentration time series $c_i(t)$ by finding the mode (max(density)) of the probability density distribution (pdf) of the instantaneous concentration differences of an individual analyzer $c_i(t)$ relative to the field average concentration $\tilde{c}(t)$ according to:

$$\Delta c_i = \max(\text{pdf}(c_i(t) - \tilde{c}(t))) \quad (5)$$

with $\tilde{c}(t)$ defined in Eq. (1) and the total number of analyzers $n=10$.
Identifying the mode of the pdf requires a robust estimate of the distribution. A comparison of histogram based and kernel-density-estimator based approaches showed that the latter are superior in terms of robustness relative to scatter in the distribution, which is a valuable feature particularly for limited sample sizes. Density estimates were generated using a moving window Gaussian kernel for smoothing (Wand and Jones, 1995). The optimal width of the window was adaptively and automatically found using pilot-density-estimates according to Sheather and Jones (1991), implemented in the *dpik* function of the *KernSmooth* library (Ripley, 2009) provided with R (R Development Core Team, 2009), also providing the *bkde* function which was used to estimate the density.

Having found an individual bias value for every analyzer, the mixing index was checked to decide whether concentration time series correction was applicable. For well mixed conditions, i.e. $MI \ge MI_c$, the observed 60-minute concentration time series $c_i(t)$ of every analyzer was shifted by the analyzer specific bias value $\Delta c_i$ found for the given 60-minute interval, yielding the bias corrected concentration time series $c_{i,corr}(t)$ according to Eq. (4). For $MI < MI_c$ the correction was applied using the last valid bias value satisfying $MI \ge MI_c$.

In order to verify that concentration offsets $\Delta c_i$ found are related to slow drift of the analyzers (instrument bias) rather than driven by meteorological forcing of natural concentration gradients, a regression analysis was performed studying the correlation of $\Delta c_i$ versus ambient air temperature, pressure and atmospheric stability $\zeta$, respectively. The stability parameter $\zeta$ is defined as $\zeta = (z - d)L^{-1}$ with measurement height $z$, displacement height $d$ and Obukhov-length $L$. No significant correlation was found between the concentration offset and the three meteorological parameters, which is an indication that the calculated offset $\Delta c_i$ is dominated by instrument bias and should therefore be removed with the proposed conditional bias correction approach, respecting $MI \ge MI_c$.

Because, even under mixed conditions, natural concentration differences could account for a (very small) part of the observed concentration offsets $\Delta c_i$, an error analysis was performed. The aim was to quantify the benefit of the application of the bias correction approach in a hypothetical “worst case” scenario, i.e. assuming that observed concentration offsets $\Delta c_i$ are solely determined by natural concentration differences rather than instrument bias. A relative error is defined in Eq. (6), describing the ratio of the error possibly attributed to the bias correction approach to the improvement achieved by the correction, which can be expressed as the span of the range of instrument bias (“drift span”). This relative error writes:

$$ \text{error}_{\text{rel}} = \frac{Q_1(\Delta \text{off}_i) - Q_4(\Delta \text{off}_i)}{\max(\text{off}_i) - \min(\text{off}_i)} \quad (6) $$

with the change of the concentration offset $\Delta c_i$ between two

consecutive 60-minute intervals $\Delta off_i = \Delta c_i(t) - \Delta c_i(t - 60 \text{ min})$ and with $Q_1$ and $Q_4$ being the 25% and 75% quartiles of the density distribution, respectively, which reflect a typical range of $\Delta off_i$.

## 2.8 Net ecosystem exchange and horizontal advection

This section indicates the relevance of measurements of CO₂ concentration gradients for the quantification of the exchange of CO₂ across the vegetation-atmosphere interface, i.e. the Net Ecosystem Exchange of CO₂ (NEE). NEE can be calculated according to the following formula (Aubinet et al., 2003; Feigenwinter et al., 2004, and others):

$$ \text{NEE} = \frac{1}{V_m} \int_{0}^{h} \left( \frac{\partial \bar{c}}{\partial t} \right) dz + \frac{1}{V_m} \left( \bar{w}' c' \right)_h + \frac{1}{V_m} \int_{0}^{h} \left( \bar{w}(z) \frac{\partial \bar{c}}{\partial z} + \bar{c}(z) \frac{\partial \bar{w}}{\partial z} \right) dz + \frac{1}{V_m} \int_{0}^{h} \left( \bar{u}(z) \frac{\partial \bar{c}}{\partial x} + \bar{v}(z) \frac{\partial \bar{c}}{\partial y} \right) dz \quad (7) $$

with the molar volume of dry air $V_m$, CO₂ concentration $c$, time $t$, horizontal distances $x$ and $y$, vertical distance above ground $z$, height of the control volume $h$, horizontal wind velocity $u$ along the x-direction, horizontal wind velocity $v$ along the y-direction and vertical wind velocity $w$ along the z-direction. Over-bars denote temporal means and primes denote the temporal fluctuations relative to the temporal mean. The terms on the right hand side of Eq. (7) are the change of storage (term I), the vertical turbulent flux (term II), vertical advection (term IIIa), vertical mass flow from the surface e.g. due to evaporation (term IIIb) according to Webb et al. (1980), and horizontal advection (term IV). The form of NEE presented in Eq. (7) excludes the horizontal variation of the vertical turbulent flux and the horizontal variation of vertical advection. Eq. (7) further neglects the flux divergence term: $\frac{1}{V_m} \int_{0}^{h} \left( \frac{\partial (\bar{u}'c')}{\partial x} + \frac{\partial (\bar{v}'c')}{\partial y} \right) dz$. Term II and sometimes terms I and III on the right hand side of Eq. (7) are central components of routine flux measurements at many sites and will not be discussed here. In contrast, term IV, the observation of which is challenging and has only been realized in a limited number of experiments, shall be addressed here.

Accurate observations of horizontal concentration gradients of CO₂ are important for the determination of horizontal advection $F_{HA}$, because $F_{HA}$ is the product of the horizontal wind velocity and the horizontal concentration gradient of the scalar CO₂ according to Eq. (8):

$$ F_{HA} = \frac{1}{V_m} \int_{0}^{h} \left( \bar{u}(z) \frac{\partial \bar{c}}{\partial x} + \bar{v}(z) \frac{\partial \bar{c}}{\partial y} \right) dz. \quad (8) $$
The fact that density distributions of concentration differ-
ences can have a mode of zero and a non zero mean, as seen
in Fig. 5b, is a crucial feature for the computation of horizon-
tal advection, because only a non zero mean gradient $\frac{\partial \tilde{c}}{\partial x} \neq 0$
and/or $\frac{\partial \tilde{c}}{\partial y} \neq 0$ can generate a non zero horizontal advection
term $F_{HA}$.

3 Results

After presenting results of the LES study, which contribute to
the acceptance of the assumptions stated in Sect. 2.4.2, this
section subsequently presents results of measured CO₂ con-
centration time series and gradients before and after applying
the conditional bias correction as well as statistics about the
improvement which can be achieved by the correction. Fur-
thermore, observed concentration differences are put in the
context of atmospheric stratification.

The results of the LES study demonstrate that for the given
simulation the assumption stated in Sect. 2.4.2 is valid, i.e.
the mode of the density distribution of the concentration dif-
ference between any sample point and the sample point field
average is essentially zero, Fig. 5a. Since both case A and
case B lead to the same conclusion, only data of case B are
shown in Fig. 5. Observed deviations of the density distri-
bution mode from zero are insignificant, with the maximum
deviation, considering all instrument's distributions, divided
by the mean distance of the sample point from the sample
point field center, accounting for a 2.0% fraction only of the
prescribed concentration gradient in the LES (case A). For
case B the maximum deviations of the mode from zero were
+0.015 and -0.025 µmol mol⁻¹. Dividing this range of dis-
tribution modes by the range of the distributions means yields
a fraction of 0.15. Considering the small gradients under well
mixed conditions, this is a very small error. Conditions with
large gradients are not an issue because they are excluded by
the mixing index filter and are not used to determine concen-
tration offsets when applying the bias correction approach.

The given deviations of the pdf's modes translate to an
error attributed to estimates of the horizontal advective flux
component, if estimates are based on concentration measure-
ments corrected using the bias correction approach and thus
removing the small deviation of the mode from zero. This
potential error in the advection estimate is small compared to
other uncertainties typically associated with advection esti-
mates, e.g. due to an insufficient number of sampling points
in space such as the often limited number of observation
height levels of horizontal gradients.

An important feature of the density distributions shown is
their skewness, separating mode and mean of a given distri-
bution as illustrated in Fig. 5b for two selected sample points.
The difference in the mean values of the density distributions
is due to the concentration gradient and source-sink distri-
bution prescribed in the LES. It thus demonstrates that it is
possible to compute advective flux terms even from distri-

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Density distribution of LES modelled concentration differences $c_i(t) - \tilde{c}(t)$ of a point measurement $c_i(t)$ relative to the field average concentration $\tilde{c}(t)$ for concentration time series $c_1(t), c_2(t), \dots, c_{16}(t)$ and $n=16$ sensor locations 1, 2, \dots, 16 (**a**), and for $c_1(t)$ and $c_{12}(t)$ at sensor locations 1 and 12 (**b**). Note that the density distributions of $c_1(t) - \tilde{c}(t)$ and $c_{12}(t) - \tilde{c}(t)$ have a common mode but different mean.</figcaption></figure>

butions with mode equal to zero, since the mean gradient,
which is necessary to compute F<sub>HA</sub> according to Eq. (8), is
expressed in the mean which does not need to be zero even
though the mode is essentially zero.

In order to evaluate the performance of the bias correc-
tion, Fig. 6a shows the CO₂ concentration evolution dur-
ing one day measured at ten locations in the sub-canopy on
29 June 2008 without bias correction but including calibra-
tion using known reference gas standards. Figure 6b presents
the same data after applying the bias correction. The com-
parison of the two figures clearly demonstrates that the bias
correction is able to remove systematic concentration offsets
between different analyzers in the uncorrected measurements
(Fig. 6a). The offsets are most obvious during well mixed
daytime conditions – when natural concentration differences
are relatively small – and could be eliminated successfully in
the bias corrected time series at all times of the day (Fig. 6b).

Inter-instrument bias leads to relatively constant offsets
between individual concentration measurements $c_i(t)$ dur-
ing daytime conditions (Fig. 6a), exactly matching the pe-
riod of a high mixing index (Fig. 4b). The minor importance
of concentration differences due to natural gradients during
well mixed conditions is the reason why inter-instrument
bias becomes the prominent component of observed inter-
instrument concentration differences (compare also Fig. 9
and Fig. 10). Well mixed conditions with $\text{MI} \ge \text{MI}_\text{c}$ and
$\text{MI}_\text{c} = 0.13$ were observed every day during the experiment,
accounting for 30% of the entire data set. There are a few
cases where mixed conditions are present for short isolated
periods (e.g. one or two 60-minute MI values) at transition
times in the early morning or sometimes in the early evening.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Calibrated CO₂ concentration time series, *before* bias correction (a) and *after* bias correction with MIc = 0.12 (b), measured at ten sub-canopy locations M5, M6, ..., M14 at a 2.25 m height with a 1 Hz resolution on 29 June 2008.</figcaption></figure>

While Fig. 6a and 6b presented CO₂ concentration time series before and after bias correction on 29 June 2008, Fig. 7 displays an example of corresponding density distributions of concentration differences during a well mixed 60-minute period at midday of the same day, which were used during bias correction. Probability density distributions with analyzer specific non-zero distribution modes in the uncorrected data of Fig. 7a have been shifted by their mode so that the new mode of the distributions is equal to zero after bias correction (Fig. 7b). Figure 7b also emphasizes sample location specific differences of the distribution shape, such as different skewness and kurtosis, which is an effect of natural concentration gradients being unique for every sample location.

Having discussed probability density distributions above for an ideal case with mixed conditions, Fig. 8 demonstrates the effect of atmospheric stratification (ζ) and the degree of mixing (MI) on the shape of selected 60-minute probability density distributions of concentration differences, which mark typical conditions during the course of a fair weather

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Density distribution of concentration differences $c_i(t) - \tilde{c}(t)$ using measured 60-minute concentration timeseries $c_i(t)$ *before* bias correction (a) and *after* bias correction (b). Number of sample locations $n = 10$ (M5, M6, ..., M14), sampling resolution 1 Hz, on 29 June 2008, 12:00–13:00. Legend indicates measurement locations according to Fig. 1.</figcaption></figure>

day, 29 June 2008. Distributions of the well mixed case in Fig. 8c are unimodal and show high kurtosis. This is beneficial for the reliable estimation of the mode, which is necessary for bias correction. High kurtosis is a consequence of small natural horizontal and vertical gradients during well mixed conditions in the middle of the day. Figure 8b and d represent transition periods between night and day and between day and night, respectively, while Fig. 8a and e are examples of night time conditions, with Fig. 8e being a representative example for conditions with katabatic sub-canopy drainage flow under very stable conditions. The kurtosis of the distributions correlates with ζ (indicator for atmospheric stratification) as well as with MI (indicator for turbulent mixing), the result being that kurtosis decreases and skewness often increases with increasing stability parameter ζ and decreasing mixing index MI. This is due to large horizontal and vertical scalar concentration gradients during
<figure><img src="image_10.png" /></figure>

**Fig. 8.** Density distribution of concentration differences $c_i(t) - \tilde{c}(t)$ using bias corrected measured 60-minute concentration timeseries $c_i(t)$, number of sample locations $n = 10$, for five typical cases over the course of the day on 29 June 2008 with varying stability parameter $\zeta$ (measured at a 36 m height) and mixing index MI (according to Eq. 3), night time, 01:00–02:00, $\zeta = -0.16$, MI=0.015 (a), night-day transition, 07:00–08:00, $\zeta = 0.65$, MI=0.020 (b), daytime, 12:00–13:00, $\zeta = -0.27$, MI=0.218 (c), day-night transition, 19:00–20:00, $\zeta = 0.06$, MI=0.010 (d) and nighttime with katabatic drainage flow, 22:00–23:00, $\zeta = 19.50$, MI=0.016 (e). Legend indicates measurement locations according to Fig. 1.

such conditions, also potentially causing multimodal distributions (Fig. 8b and d), which can lead to disambiguities concerning the relevant mode if they were to be used for bias correction, which they are not due to the mixing index condition. However, the effect of atmospheric stability $\zeta$ is not uniform, meaning that multiple modes and skewed distributions (Fig. 8b) and low kurtosis (Fig. 8d) are more pronounced during transition periods with moderate vertical exchange, whereas the night time cases such as Fig. 8e with the highest stability parameter $\zeta$ and least vertical exchange are less skewed and more homogeneous with respect to different sample locations. The absence of vertical exchange results in horizontally relatively homogeneous sub-canopy scalar concentrations even though there are large vertical gradients.

Figures 9 and 10 demonstrate that observed concentration offsets $\Delta c_i$ can be separated into offsets which are mainly determined by instrument bias alone and into offsets which are determined by instrument bias as well as significant natural concentration differences. Figure 9a shows offset time series over two days with a succession of mixed daytime conditions (approx. 8 h to 16 h) with little scatter in the offset time series when natural gradients are small and offsets are mainly con-

<figure><img src="image_11.png" /></figure>

**Fig. 9.** Time series of the modes of density distributions of concentration differences $c_i(t) - \tilde{c}(t)$ (see Fig. 7a for example distributions for one 60-minute time step) for 10 sampling locations on 29 June and 30 June 2008, *before* filtering with mixing index (a), and *after* filtering with mixing index $MI_c = 0.13$ (b). Modes from periods which satisfy $MI < MI_c$ are not used during bias correction (grey mask). The last mode at a time with $MI \ge MI_c$ is used instead (solid lines).

trolled by low frequency instrument bias and night time conditions with high scatter and large absolute values in the offset time series when natural gradients are the predominant cause. After applying the mixing index to filter the offset time series, those periods with predominant natural gradients were effectively excluded (Fig. 9b). The remaining offsets are controlled by instrument bias and can therefore be used in the bias correction approach.

The different offset characteristics during daytime and nighttime described above are due to the dependence of natural concentration differences on the mixing index and atmospheric stability, both of which have a distinct daily cycle. Figure 10 illustrates the dependence of concentration offsets on the mixing index MI. For low values of MI, natural horizontal gradients are large, as a result of horizontal source heterogeneities and potential mixing of a vertical concentration
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Modes of 60-minute density distributions of concentration differences $c_i(t) - \bar{c}$ minus analyzer drift (mode of 24 h pdf of 60-minute modes subtracted daily) versus mixing index MI. Grey points indicate measurements, the solid line marks the 25% and 75% quantiles for mixing index binwidths of 0.025 with the circles centered at each bin. Dashed line at MI=0.12 indicates a sensible choice for the critical mixing index $MI_c$.</figcaption></figure>

profile with large vertical gradients. For larger values of MI, particularly for $MI > MI_c$, offsets are relatively small. Figure 10 displays the dependence of only natural concentration differences on MI. The graph shows offsets with instrument bias removed by subtracting a low frequency component (see figure caption for details). The fact that this technique does not perfectly separate natural concentration differences from instrument bias explains the scatter and outliers in Fig. 10 which are present even at higher values of MI. The majority of data points (indicated by solid lines for the 25% and 75% quantiles) in Fig. 10 is quite close to zero concentration difference for higher values of MI ($MI > MI_c$). That indicates that there are no major natural concentration differences under those conditions which could unintentionally be removed by the bias correction approach.

To compare the benefits of the bias correction approach with potential errors, Table 1 displays results of an error analysis, listing the potential for improvement by using the bias correction approach (drift span), an estimate of the potential absolute error ($Q_4(\Delta off_i) - Q_1(\Delta off_i)$) and the relative error ($error_{rel}$) for ten sampling locations. Values of the relative error are on the order of 10%, which is a satisfying result, keeping in mind that those are “worst case” values pretending that offsets during mixed conditions, i.e. when the bias correction is applied, were purely caused by natural gradients, which they are not in reality. Therefore the true relative error will be even smaller than values given in Table 1 for $error_{rel}$.

<table><caption>**Table 1.** Offset statistics and error analysis for ten sample points, i.e. ten analyzers, demonstrating the correction potential of the bias correction approach (“drift span”), typical values for the maximum error possibly attributed to the bias correction approach for a “worst case” scenario (from quartile $Q_1(\Delta off_i)$ to quartile $Q_4(\Delta off_i)$) and their ratio, i.e. the relative error $error_{rel}$ according to Eq. (6). See Sect. 2.7 for definition of the terms.</caption><thead><tr><th>Sample point</th><th>drift span</th><th>Q<sub>1</sub>(Δoff<sub>i</sub>) [μmol mol<sup>-1</sup>]</th><th>Q<sub>4</sub>(Δoff<sub>i</sub>) [μmol mol<sup>-1</sup>]</th><th>error<sub>rel</sub> [] [μmol mol<sup>-1</sup>]</th></tr></thead><tbody><tr><td>M5</td><td>9.4</td><td>-0.46</td><td>0.41</td><td>0.09</td></tr><tr><td>M6</td><td>10.4</td><td>-0.86</td><td>0.72</td><td>0.15</td></tr><tr><td>M7</td><td>8.1</td><td>-0.53</td><td>0.47</td><td>0.12</td></tr><tr><td>M8</td><td>7.2</td><td>-0.7</td><td>0.49</td><td>0.16</td></tr><tr><td>M9</td><td>6.1</td><td>-0.36</td><td>0.41</td><td>0.13</td></tr><tr><td>M10</td><td>23.1</td><td>-1.11</td><td>0.93</td><td>0.09</td></tr><tr><td>M11</td><td>9.2</td><td>-0.89</td><td>0.66</td><td>0.17</td></tr><tr><td>M12</td><td>23</td><td>-0.43</td><td>0.45</td><td>0.04</td></tr><tr><td>M13</td><td>14.3</td><td>-0.55</td><td>0.54</td><td>0.08</td></tr><tr><td>M14</td><td>12.4</td><td>-0.58</td><td>0.56</td><td>0.09</td></tr></tbody></table>

## 4 Discussion

There are three issues connected to the quality of the bias correction approach which shall be discussed in this section: potential underestimation of natural concentration differences (signal loss), tradeoff between limiting instrument drift and limiting signal loss, and finding the appropriate window length when applying the bias correction.

The previous section has shown that the improvements that were achieved by applying the bias correction approach are one order of magnitude larger than possible errors associated with it, which is a strong incentive to use the correction approach. However, there is potential for losing part of the natural concentration gradients when applying the correction, due to possibly imperfect separation of instrument bias and concentration differences originating from a natural gradient, even during mixed conditions. A quantification of this phenomenon was given in Table 1. The acceptance of this relatively small potential error when applying the bias correction approach needs to be compared to errors which are likely to be attributed to the gradient measurements with no correction applied. It is known from various advection experiments that instrument related bias between sampling points can be on the order of the natural horizontal concentration gradients, particularly at relatively homogeneous sites. This in turn can lead to considerable overestimation of the absolute value of horizontal advection, which is one of the reasons why including the horizontal advection flux term in the Net Ecosystem Exchange (NEE) budget equation often leads to increased scatter of NEE and does not necessarily produce reliable NEE estimates. As a consequence, NEE is often computed using the turbulent and storage fluxes only. We suggest that rather than including a noisy and potentially too large advection estimate in the NEE equation, it is better
to include a bias corrected estimate of horizontal advection. Doing so and at the same time accounting for vertical advection – the same arguments apply here as to avoid overestimation and noise – should give more realistic NEE estimates than those obtained from turbulent and storage flux alone.

When applying the bias correction, a balance should be found between limiting the effect of instrument drift on the gradient measurements and signal loss by potential underestimation of natural gradients. This balance can be tuned by the choice of the value for the critical mixing index $MI_c$. A high value of $MI_c$ better preserves natural gradients because bias correction values are only determined from data during well mixed conditions and therefore can not eliminate natural gradients during other conditions, particularly at night when natural gradients are typically large. A low value of $MI_c$ removes instrument drift more thoroughly since bias correction values can be found more often, i.e. from well mixed as well as partly mixed conditions. Therefore, we recommend to choose a higher $MI_c$ the more stable the analyzer is and just low enough to allow the instrument to “survive” periods during which no bias correction values can be found (i.e. nighttime) using previously established correction values (inherited from daytime) without facing prohibitive instrument drift during those periods.

The third issue is finding the appropriate window length $T_F$ when applying the bias correction. This is the length of the time series used to compute density distributions of concentration differences (pdf) to find their mode as outlined in Sect. 2.7. For this study the window length was chosen to be $T_F = 60$ min. The higher the instrument drift is, the shorter this window has to be in order to find a mode which is representative for the instrument bias during that time window and not affected by a significant trend of the bias. On the other hand, choosing the window as long as possible helps to preserve natural gradients which are persistent for longer periods, since persistent natural gradients with periods longer than the window length and present during non mixed conditions, and therefore affecting the mode of the pdf, are removed by the bias correction for $MI > MI_c$. However, we can conclude from the data that it is not satisfactory to choose an infinite window length (such as the time constant bias correction applied by Aubinet et al., 2003) in order to preserve natural gradients because observed instrument bias is subject to drift over time. Given the window length of 60-minutes used in the current study, the concentration difference error due to signal loss of natural concentration differences during the day has been shown to be smaller than the error of the concentration offset which would be caused by the drift of the instrument bias if the latter was corrected by a time constant bias value. Future studies should test window lengths larger than 60-minutes, particularly when using more stable analyzers.

Future work on the improvement of the bias correction approach should include a refined condition to test which data should be used when determining the pdf and the bias.

Rather than using fixed 60-minute intervals to determine $MI$ and accepting all data in a 60-minute interval satisfying $MI \ge MI_c$, a more fine grained selection of data entering the pdf can be used to select only those parts of the time series which have common properties at more than one sample point for a time period on the order of the duration of coherent structures, i.e. seconds to minutes. Among the tools which can be used to find common properties within the time series are cross correlation analysis and pattern recognition. Thereby only data with similar concentration at several sample points will enter the pdf. This helps to exclude the influence of natural gradients on the mode of the pdf, which will then be determined by instrument bias alone. Such short term correlation of time series at several sample points by tracking individual structures in the time series should be done for sample point pairs rather than using properties of the complete sample point field. These pair wise correlations then need to be linked together by choosing different configurations of sample point pairs and combining their information.

Future work should also test the applicability of the bias correction approach to sensor networks with a possibly large number of sampling points. The approach could be used when working with sensors which have a relatively high resolution but suffer from low accuracy. Such sensors can deliver the fine structure (high frequency part) of the time series. The bias correction approach corrects constant and drifting instrument bias (low frequency part) and thus ensures the relative accuracy of the measurements.

## 5 Conclusions

This paper has presented a measurement design capable of addressing the issue of inadequate sampling of natural concentration gradients in the temporal domain – a common characteristic of many advection measurement setups – by increasing the temporal resolution of the gradient measurements. Observing gradients with a sufficiently high temporal resolution and therefore capturing as much information as possible over a large range of temporal scales is crucial for reliable advection estimates computed from concentration gradients. In order to produce accurate gradient measurements in a multi-analyzer setup, an approach was presented which adequately addresses the problem of interinstrument bias. It was shown that the uncertainties associated with this approach are one order of magnitude smaller compared to the benefit achieved for the given setup. For completeness it should be stated that it is always advisable to avoid instrument bias as far as possible by appropriate technical measures, e.g. sampling system design and calibration against known reference gas standards (conventional calibration). However, presented statistical calibration method is independent of conventional calibration in the sense that it solely deals with the remaining bias after conventional calibration and that the functioning of the statistical calibration
is to the largest extend independent of the magnitude of this remaining bias. The proposed bias correction approach is therefore a suitable tool at least for multi-analyzer setups measuring horizontal gradients at one height, given a certain proximity of individual sampling locations. There might also be benefits from applying the bias correction approach to sequentially measured data from switching valve systems in a single-analyzer setup. It should be tested in the future whether the bias correction approach can be transferred to measurements of vertical gradients, although care has to be taken due to strong systematic vertical gradients particularly at night in the case of CO₂ concentration. The concept outlined in the current paper should not be limited to measurements of CO₂ concentration but be useful for the accurate observation of gradients of other scalars, too. Furthermore, it need not be limited to gradient measurements for the computation of advective flux components but is worth considering for any gradient based flux measurement application. Finally, the bias correction approach might be useful for the relative adjustment of signal levels between individual sensors in any kind of sensor network that samples phenomena which – at least part of the time – lead to common characteristics of the observed signal at several locations in the network. We therefore propose to test presented method at other experimental sites measuring CO₂ concentration gradients and to explore above mentioned additional applications.

*Acknowledgements.* The authors wish to acknowledge the help and technical support performed by the staff of the Bayreuth Center for Ecology and Environmental Research (BayCEER) of the University of Bayreuth. The experiment was funded by the German Science Foundation (FO 226/16-1, ME 2100/4-1, ZE 792/4-1). We would further like to acknowledge S. Raasch and his group for developing and providing the LES tool PALM. The LES runs were carried out on the HLRN-II-system of the HLRN (North German Alliance for Supercomputing).

Edited by: T. von Clarmann

## References

Aubinet, M., Heinesch, B., and Yernaux, M.: Horizontal and vertical CO₂ advection in a sloping forest, Bound.-Lay. Meteorol., 108, 397–417, 2003.

Aubinet, M., Feigenwinter, C., Heinesch, B., Bernhofer, C., Canepa, E., Lindroth, A., Montagnani, L., Rebmann, C., Sedlak, P., and Gorsel, E. V.: Direct advection measurements do not help to solve the night-time CO₂ closure problem: Evidence from three different forests, Agric. For. Meteorol., 150, 655–664, 2010.

Baldocchi, D., Finnigan, J. J., Wilson, K., Paw U, K. T., and Falge, E.: On measuring net ecosystem carbon exchange over tall vegetation on complex terrain, Bound.-Lay. Meteorol., 96, 257–291, 2000.

Baldocchi, D. D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y.,

Meyers, T., Munger, W., Oechel, W., Paw, K. T., Pilegaard, K., Schmid, H. P., Valentini, R., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: FLUXNET: A New Tool to Study the Temporal and Spatial Variability of Ecosystem-Scale Carbon Dioxide, Water Vapor, and Energy Flux Densities, B. Am. Meteorol. Soc., 82, 2415–2434, 2001.

Burns, S. P., Delany, A. C., Sun, J., Stephens, B. B., Oncley, S. P., Maclean, G. D., Semmer, S. R., Schröter, J., and Ruppert, J.: A Programmable Portable Trace-Gas Measuring System and an Evaluation of Calibration Techniques for In-Situ Carbon Dioxide Measurements, J. Atmos. Oceanic Tech., 26, 291–316, 2009.

Canepa, E., Georgieva, E., Manca, G., and Feigenwinter, C.: Application of a mass consistent flow model to study the CO₂ mass balance of forests, Agric. For. Meteorol., 150, 712–723, 2010.

Feigenwinter, C., Bernhofer, C., and Vogt, R.: The influence of advection on the short term CO₂-budget in and above a forest canopy, Bound.-Lay. Meteorol., 113, 201–224, 2004.

Feigenwinter, C., Bernhofer, C., Eichelmann, U., Heinesch, B., Hertel, M., Janous, D., Kolle, O., Lagergren, F., Lindroth, A., Minerbi, S., Moderow, U., Molder, M., Montagnani, L., Queck, R., Rebmann, C., Vestin, P., Yernaux, M., Zeri, M., Ziegler, W., and Aubinet, M.: Comparison of horizontal and vertical advective CO₂ fluxes at three forest sites, Agric. For. Meteorol., 148, 12–24, 2008.

Finnigan, J.: An introduction to flux measurements in difficult conditions, Ecol. Appl., 18, 1340–1350, 2008.

Gerstberger, P., Foken, T., and Kalbitz, K.: The Lehstenbach and Steinkreuz Catchments in NE Bavaria, Germany, in: Biogeochemistry of Forested Catchments in a Changing Environment: A German Case Study, edited by: Matzner, E., 172, 15–41, Springer, Heidelberg, 2004.

Heinesch, B., Yernaux, M., and Aubinet, M.: Some methodological questions concerning advection measurements: a case study, Bound.-Lay. Meteorol., 122, 457–478, 2007.

Kowalski, A. and Serrano-Ortiz, P.: On the relationship between the eddy covariance, the turbulent flux, and surface exchange for a trace gas such as CO₂, Bound.-Lay. Meteorol., 124, 129–141, 2007.

Lee, X.: On micrometeorological observations of surface-air exchange over tall vegetation, Agric. For. Meteorol., 91, 39–49, 1998.

Mahrt, L.: Computing turbulent fluxes near the surface: Needed improvements, Agric. For. Meteorol., 150, 501–509, 2010.

Montagnani, L., Manca, G., Canepa, E., and Georgieva, E.: Assessing the method-specific differences in quantification of CO₂ advection at three forest sites during the ADVEX campaign, Agric. For. Meteorol., 150, 702–711, 2010.

R Development Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org, ISBN 3-900051-07-0, 2009.

Raasch, S.: PALM group available online under, http://palm.muk.uni-hannover.de, 2010.

Raasch, S. and Schröter, M.: PALM – A large-eddy simulation model performing on massively parallel computers, Meteorologische Zeitschrift, 10, 363–372, 2001.

Ripley, B.: KernSmooth: Functions for kernel smoothing for Wand & Jones (1995), http://CRAN.R-project.org/package=KernSmooth, S original by: Matt Wand. R port by: Brian Ripley.
R package version 2.23-3, 2009.

Serafimovich, A., Siebicke, L., Staudt, K., Lüers, J., Hunner, M., Gerken, T., Schier, S., Biermann, T., Rütz, F., von Buttlar, J., Riederer, M., Falge, E., Mayer, J.-C., and Foken, T.: ExchanGE processes in mountainous Regions (EGER) – Documentation of the Intensive Observation Period (IOP2) June, 1 to July, 15 2008, Arbeitsergebnisse, 37, University of Bayreuth, Department of Micrometeorology, ISSN 1614-8916, 180 pp., 2008.

Sheather, S. J. and Jones, M. C.: A reliable data-based bandwidth selection method for kernel density estimation, J. R. Stat. Soc., Series B, 53, 683–690, 1991.

Siebicke, L.: Interactive comment on “CO₂ gradient measurements using a parallel multi-analyzer setup” by L. Siebicke et al., Atmos. Meas. Tech. Discuss., 3, C2390–2406, 2011.

Staebler, R. and Fitzjarrald, D.: Observing subcanopy CO₂ advection, Agric. For. Meteorol., 122, 139–156, 2004.

Staudt, K. and Foken, T.: Documentation of reference data for the experimental areas of the Bayreuth Centre for Ecology and Environmental Research (BayCEER) at the Waldstein site, Arbeitsergebnisse 35, University of Bayreuth, Department of Micrometeorology, ISSN 1614-8916, 37 pp., 2007.

Sun, J., Burns, S. P., Delany, A. C., Oncley, S. P., Turnipseed, A. A., Stephens, B. B., Lenschow, D. H., LeMone, M. A., Monson, R. K., and Anderson, D. E.: CO₂ transport over complex terrain, Agric. For. Meteorol., 145, 1–21, 2007.

Vickers, D. and Mahrt, L.: Contrasting mean vertical motion from tilt correction methods and mass continuity, Agric. For. Meteorol., 138, 93–103, 2006.

Wand, M. P. and Jones, M. C.: Kernel Smoothing, Chapman and Hall, London, 1995.

Webb, E. K., Pearman, G. I., and Leuning, R.: Correction of flux measurements for density effects due to heat and water vapour transfer, Q. J. R. Meteorol. Soc., 106, 85–100, 1980.