Stomatal and Non-Stomatal Turbulent Deposition Flux of Ozone to a Managed Peatland

: Ozone is a key trace gas in the troposphere; because it is a greenhouse gas, it is very reactive, and it is potentially toxic to humans, fauna, and vegetation. The main sink processes for ozone are chemical reactions and the turbulent deposition ﬂux to the earth’s surface. The deposition process itself is rather complex: The interactions between co-varying drivers such as the tropospheric ozone concentration, turbulence, and chemical reactions are not well understood. In the case of ozone deposition to vegetation, another aspect that must be studied is the role of stomatal regulation for a wide range of conditions. Therefore, we measured turbulent deposition ﬂuxes of ozone with the eddy covariance technique during the peak of the growing season in 2014 over a managed, rewetted peatland in NW Germany. The deposition ﬂux was large during the day (up to − 15 nmol m − 2 s − 1 ) and relatively small during the night (between − 1 and − 2 nmol m − 2 s − 1 ). Flux partitioning by applying the surface resistance analogy and further analysis showed that the stomatal uptake was smaller than non-stomatal deposition. The correction of stomatal conductance with the gross primary production (GPP) improved the estimation of day- and nighttime stomatal deposition ﬂuxes. Statistical analysis conﬁrmed that the friction velocity (u*) was the single most important driver of non-stomatal ozone deposition and that relationships with other environmental drivers are not linear and highly variable. Further research is needed to develop a better process understanding of non-stomatal ozone deposition, to quantify the role of surface deposition to the ozone budget of the atmospheric boundary layer, and to estimate uncertainties associated with the partitioning of ozone deposition into stomatal and non-stomatal ﬂuxes.


Introduction
Tropospheric ozone is an important greenhouse gas and is harmful to humans, as it impairs lung function [1]. Its average mixing ratio has more than doubled over the past 150 years [2]; the background mixing ratio today in Europe is approximately 30 ppb and is continuously increasing [3]. Ozone is not directly emitted into the atmosphere by anthropogenic emissions, but it is formed through photochemical reactions involving precursors-mainly nitrogen oxides (NO x ) and volatile organic compounds (VOC)-from both natural and anthropogenic origins [4]. The sinks of tropospheric ozone are manifold and have been studied previously in detail: Ozone reacts with various trace gases such as hydroperoxyl radicals (HO 2 ), biogenic VOC [5], and nitric oxide (NO) [6,7]. It also reacts with aerosol particles [8,9] and undergoes direct photolysis through short-wave solar radiation [6]. Another sink mechanism is its deposition on the earth's surface [10]. Once in contact with a given surface (vegetation, soil, or water), ozone will react readily and be decomposed; thus, its deposition is a terminal process [11,12]. For vegetation, two ozone deposition pathways act in parallel: First, ozone is deposited on plants' surfaces, such as on their cuticles or twig and stem surfaces. Second, ozone can enter the plants' stomata, reach the sub-stomatal cavity and thus impact photosynthesis and other physiological processes [13,14]. This may eventually cause alterations of plant morphology [15] and crop yield losses [16]. However, the stomatal ozone flux preceding these effects cannot be directly quantified on an ecosystem scale. Efforts to quantify this flux by indirect methods have only been made for a few ecosystems so far [17][18][19][20][21][22][23].
In this study, we focus on the ozone deposition process on vegetation, specifically on the partitioning between stomatal and non-stomatal fluxes [22,24,25]. We performed this study on a mid-latitude peatland ecosystem, which is managed with the goal of bringing it back to a natural or at least semi-natural state. While most studies on ozone deposition to peat-and wetlands have been conducted using mesocosm and microcosm experiments [26][27][28][29], only few studies have considered the larger picture and deal with ozone deposition to these ecosystems as a whole [30][31][32][33]. To our knowledge, there is only one study on ozone deposition to a mid-latitude peatland on the ecosystem scale [21]; thus, here we aim to add valuable information that can improve the characterization of ozone deposition to these ecosystems.
We applied the eddy covariance technique as well as the surface resistance analogy to first quantify the overall ozone deposition flux (F O 3 ) and then partition this flux into stomatal (i.e., uptake by vegetation through stomatal gas exchange, F sto ) and non-stomatal fluxes (F nsto ) during the peak of the growing season. We performed a driver analysis of non-stomatal deposition velocities and analyzed the relationship between the stomatal ozone fluxes and gross primary productivity under different environmental ozone concentrations.

Study Site and Study Period
In NW Germany, natural peat has been widely cut through the past centuries and is still cut today. Once drained, a large portion of the soils have been used for intensive agricultural activities. Over the past decades, some of the cut areas have started to be restored in order to regain the natural properties of these ecosystems, such as reestablishing the groundwater level to be near the soil surface.
Our study site is located within the Mittleres Wietingsmoor (MWM) and north of the village Freistatt within the Diepholzer Moorniederung, a peatland consisting mainly of bogs. This Diepholz peatland extends over 1575 ha, of which 1176 ha are under high-priority nature protection. Until the 1970s, peat was cut manually here; after that and up until 1995, it was cut using industrial machinery. Overall, up to 1.3 m of peat was cut. In 1982, restoration started in some areas and was extended and intensified after 1999 [34,35].
The actual site (52 • is within an open, rewetted grassland with Juncus effuses being the dominant species. In wetter areas, we found Eriophorum vaginatum and Sphagnum cuspidatum, while Holcus lanatus appeared in dryer patches. The average vegetation height increased from 81 cm to 85 cm from the beginning to the end of our experimental period. The area is kept free of woody vegetation through regular sheep grazing.
The experimental period was from 30 June through 4 August 2014. Due to rather frequent and sometimes extensive power losses before 15 July, we limit our data analysis to the period 15 July through 4 August 2014.

Instrumentation
The eddy covariance setup consisted of three main components. First, a CSAT3 (Campbell Scientific Logan, UT, USA) ultrasonic anemometer was employed to measure the 3D wind vector and the sonic temperature. It was mounted at a height of z = 2.42 m above ground level and oriented to the southwest. Second, a LI-7200 CO 2 /H 2 O enclosed-path gas analyzer (LI-COR, Inc., Lincoln, NE, USA) measured the CO 2 and the water vapor concentration. Its intake tube was 0.5 m long and had an inner diameter of 5.3 mm. The inlet was mounted 4 cm to the north and 22 cm above the center of the ultrasonic anemometer. Third, a CLD 88 O3 (ECO PHYSICS AG, Dürnten, Switzerland) ozone analyzer was used to measure the mixing ratio of O 3 (χ O 3 ) with the gas phase chemiluminescence method after O 3 reacted with NO, which was offered in surplus. The instrument has a measurement range from 0.5 ppb to 5000 ppb and the lag time within the analyzer is below 1 second. The inlet of a 16 mm inner diameter Teflon tube was located 13 cm to the northeast of the center of the ultrasonic anemometer. The data acquisition frequency for all instruments described so far was 5 Hz. Whenever needed, ozone mixing ratios (χ O 3 , in units ppb) were converted to concentrations (c O 3 , in units of µg m −3 or nmol m −3 ) by using actual air temperature (T air ) and pressure data.
Further, low-frequency measurements were made for air temperature (Tair) and relative air humidity (rH) with an HMP45A (Vaisala, Helsinki, Finnland) within a radiation shield at 2 m above ground. A four-channel CNR 4 radiometer (Kipp & Zonen B.V., Delft, Netherlands) was employed for the measurement of longwave and shortwave radiation (upwelling and downwelling, respectively), and a Young Typ 61302 V sensor measured the air pressure. These measurements were collected at 1 Hz and aggregated to 10-min averages. The photograph in Figure 1 (left) shows the instrumented field site.
Atmosphere 2017, 8,175 3 of 17 the ultrasonic anemometer. Third, a CLD 88 O3 (ECO PHYSICS AG, Dürnten, Switzerland) ozone analyzer was used to measure the mixing ratio of O3 ( ) with the gas phase chemiluminescence method after O3 reacted with NO, which was offered in surplus. The instrument has a measurement range from 0.5 ppb to 5000 ppb and the lag time within the analyzer is below 1 second. The inlet of a 16 mm inner diameter Teflon tube was located 13 cm to the northeast of the center of the ultrasonic anemometer. The data acquisition frequency for all instruments described so far was 5 Hz. Whenever needed, ozone mixing ratios ( , in units ppb) were converted to concentrations ( , in units of µg m −3 or nmol m −3 ) by using actual air temperature (Tair) and pressure data. Further, low-frequency measurements were made for air temperature (Tair) and relative air humidity (rH) with an HMP45A (Vaisala, Helsinki, Finnland) within a radiation shield at 2 m above ground. A four-channel CNR 4 radiometer (Kipp & Zonen B.V., Delft, Netherlands) was employed for the measurement of longwave and shortwave radiation (upwelling and downwelling, respectively), and a Young Typ 61302 V sensor measured the air pressure. These measurements were collected at 1 Hz and aggregated to 10-minute averages. The photograph in Figure 1 (left) shows the instrumented field site.

Eddy Covariance
The sensible heat flux (H), latent heat flux (LE), flux of CO2 ( ) and were calculated with the software EddyPro 6.2.0 (Licor Bioscience) [37] (see Equation (1), section 2.3.2, for the calculation of the ozone flux). Raw data processing and flux calculation included de-spiking [38], double rotation of the 3D wind components, block averaging, lag correction with covariance maximization within defined bounds, corrections for sensor separation according to Horst and Lenschow [39], highfrequency corrections according to Ibrom et al. [40], and high-pass filtering effects according to Moncrieff et al. [41]. Additionally, the ozone fluxes were corrected for the dilution of the ozone concentration by water vapor, which was not measured within the instrument [42]. The calculated fluxes were quality controlled for stationarity and integral turbulence characteristics [43]. The CO2 storage term was estimated using the 1-point storage correction and then added to the CO2 fluxes to derive the net ecosystem exchange (NEE). Nocturnal data were filtered for low friction velocity (u*) values. The u* threshold (0.1 m s −1 ) was estimated following Papale et al. [44].

Eddy Covariance
The sensible heat flux (H), latent heat flux (LE), flux of CO 2 (F CO 2 ) and F O 3 were calculated with the software EddyPro 6.2.0 (Licor Bioscience) [37] (see Equation (1), Section 2.3.2, for the calculation of the ozone flux). Raw data processing and flux calculation included de-spiking [38], double rotation of the 3D wind components, block averaging, lag correction with covariance maximization within defined bounds, corrections for sensor separation according to Horst and Lenschow [39], high-frequency corrections according to Ibrom et al. [40], and high-pass filtering effects according to Moncrieff et al. [41]. Additionally, the ozone fluxes were corrected for the dilution of the ozone concentration by water vapor, which was not measured within the instrument [42]. The calculated fluxes were quality controlled for stationarity and integral turbulence characteristics [43]. The CO 2 storage term was estimated using the 1-point storage correction and then added to the CO 2 fluxes to derive the net ecosystem exchange (NEE). Nocturnal data were filtered for low friction velocity (u*) values. The u* threshold (0.1 m s −1 ) was estimated following Papale et al. [44].
Data quality control was performed according to Mauder and Foken [45] and resulted in three classes: QC_0 data are considered to be of very good quality, QC_1 data are of acceptable quality, and QC_2 data should not be considered for scientific data analysis. Further, data sets that are associated with low u* values (u* < 0.1 m s −1 ) are considered problematic because the respective turbulence is not sufficiently established. After removal of QC_2 data and those with low u* values (see Section 3.1 for details), the remaining data set was gap-filled with marginal distribution sampling as described in Reichstein et al. [46]. Statistical analyses (Sections 3.2 and 3.3) were performed with non-gap-filled data.
Subsequently, NEE was partitioned into gross primary productivity (GPP) and ecosystem respiration (R eco ) by using daytime flux partitioning with light response curves [47] as implemented in REddyProc 1.0.0 for R [48].
The flux footprints were calculated as described in Kljun et al. [36] for half-hour fluxes when the u* threshold was exceeded. Nighttime fluxes were calculated whenever the global radiation (i.e., the shortwave, downwelling radiation R g ) was below or equal to 5 W m −2 , while daytime conditions were defined for global radiation values above 200 W m −2 . With this classification, the transition periods around sunrise and sunset were excluded from the analysis.
Co-spectra were computed for the vertical wind component in conjunction with temperature and O 3 , respectively, for each half hour. Frequency-class binned averages were then used to calculate ensemble co-spectra from good-quality data and well-developed turbulence conditions for individual stability classes. Reference co-spectra for temperature and vertical wind were calculated according to Kaimal and Finnigan [49].
The stability of the boundary layer was determined from the ultrasonic anemometer data as (z-d)/L with L being the Monin-Obukhov length and d the displacement height (d = 0.67 · z). Neutral stability conditions were defined as −0.03 ≤ z/L ≤ 0.03, stable conditions as z/L > 0.03, and unstable conditions as z/L < −0.03, respectively.

Resistance Terms and Deposition Velocity
F O 3 was measured with the eddy covariance method as described above on the basis of high-frequency measurement of c O 3 and the vertical wind component w (Equation (1) in [24]), where the primes describe the deviations of individual data points from the half-hour means and the overbar indicates the half-hour averages. Further, F O 3 can be expressed in terms of concentration gradients and resistance terms (Equation (4) in [24]): Here, c O 3 is the ambient O 3 concentration at measurement height and c 0 is the concentration within the sub-stomatal cavities, which can be assumed to be zero [50]. The total resistance against the O 3 deposition R tot is equal to the sum R a + R b + R c (Equation (2)), where the aerodynamic resistance (R a ) and the bulk leaf boundary layer resistance for ozone (R b ) are computed from ultrasonic anemometer data as (Equations (8) and (9) in [50]) and (Equation (6) Here, σ v is the standard deviation of the lateral wind component, u* is the friction velocity, u is the horizontal wind velocity, Sc = 1.07 is the Schmidt number for ozone and Pr = 0.72 is the Prandtl number [51]. The surface resistance R c can now be calculated from Equation (2) with c 0 = 0 and by using the measured ozone flux F O 3 (Equation (1)) and the ambient ozone concentration c O 3 : Following the reasoning of Gerosa et al. [25] and Fares et al. [17], F O 3 can be divided into F sto and F nsto (Equation (9) in [24]), F sto happens together with the exchange flux of CO 2 and H 2 O vapor through the plant stomata. F nsto includes all other O 3 sinks in the vegetation including plant, water, soil surfaces, and chemical reactions. It is further presumed that the flux-resistance concept can be applied for both terms on the right-hand side of Equation (6) separately while using an ozone concentration near the surfaces c O 3 sur f (Equation (9) in [24]): Note that c O 3 sur f is not measured but will be computed as a residual once R sto and R nsto are known (see below, Equations (8)- (12). The surface resistance R c , which is known from Equation (5), can be partitioned into two resistances acting in parallel, the stomatal (R sto ) and the non-stomatal (R nsto ) resistances (Equation (9) in [21]): To calculate the stomatal resistance of ozone, R sto , we use the estimate of the stomatal conductance for water vapor from the Penman-Monteith approach (g sto,PM ) and the molecular diffusivities for water vapor D H 2 O and ozone D O 3 (Equation (4) in [11]), where β is the Bowen ratio H / ET, s is the slope of the water vapor saturation curve (Pa K −1 ), and γ is the psychrometric constant (kPa K −1 ). R b,H 2 O is the bulk leaf boundary layer resistance for water vapor, which is computed from Equation (4) with the Schmidt number for water vapor (0.68; [51]). We calculated the water vapor pressure deficit (VPD) in Equation (9) from the difference between es T sur f , which is saturation vapor pressure (kPa) at the temperature at canopy surface level T surf , and e, which is the vapor pressure (kPa) at measurement level z. This procedure was proposed by Gerosa et al.
(Equation (12) in [24]), who suggest calculating es T sur f as and this calculation is valid as long as no direct evaporation takes place. To account for "computational contamination" of the stomatal water vapor conductance in Equation (9) from direct evaporation, we applied a correction after Lamaud et al. [11]. For this purpose, the slope of the regression line α between g sto,PM and GPP is computed for conditions during daytime and rH < 60% only. Under these conditions, the direct evaporation from wet surfaces is at its minimum (zero), and the stomatal conductance for water vapor and CO 2 are related by the constant α. This ratio is used to compute the stomatal conductance of ozone for all conditions (Equation (5) in [11]): Now, R sto can be directly calculated as and R nsto can be calculated from Equation (8), c O 3 sur f from Equation (7), and the flux partitioning can be derived from Equation (6). This routine was applied for each 30-min interval of the measurement period. The deposition velocity v d is typically considered as the inverse of R tot , while, by convention, a minus sign is introduced into Equation (11) (Equation (7) in [50]): Additionally, the deposition velocity can be calculated from measured data directly as (Equation (9) in [11]) In analogy to eqs. 8 and 11, ozone deposition velocities towards the stomata are defined as v d,sto = 1 R sto , and towards all non-stomatal surfaces as v d,nsto = 1 R nsto . Note that v d,sto equals g sto .

Statistical Driver Analysis
Statistical analyses were performed with R 3.2.2 [48]. A relative importance analysis was carried out using the relaimpo package [52] with the Lindeman, Merenda and Gold method (LMG). The contribution of each predictor variable (see below) to the total R 2 between the modelled and the measured non-stomatal deposition velocity was calculated by first averaging the contributions of each variable according to different ordering of predictor variables within one model, and then averaging the results of step 1 for each variable across multiple linear models of different sizes (different numbers of predictor variables included). This procedure needs to be applied when the predictor variables are correlated to each other, as their order influences the individual contributions to the total R 2 within the multiple linear model [52]. The following predictor variables were used for the relative importance analysis of F nsto , χ O 3 , VPD, T air , u*, LE, GPP, R g and rH. The selection was based on variables which were identified as potential predictors for ozone fluxes in the literature. To better understand the relationships between the drivers and v d,nsto , the latter was binned according to the various drivers for graphical display.

Regional Ozone Analysis
To understand the spatial representativeness of the ozone concentrations in the study region MWM, regression and correlation analyses of hourly c O 3 were performed between MWM and two air quality measurement stations of the State of Lower-Saxony, which are operated by the Lufthygienisches Überwachungssystem Niedersachsen (LÜN). The LÜN stations Südoldenburg (SO) and Allertal (AT) were 61 and 68 km in northwestern and east-northeastern directions from MWM, respectively. Data from the LÜN stations were provided as hourly concentrations. For this analysis, the MWM data were averaged to hourly c O 3 .

Flux Data Quality Control
The eddy covariance method is widely employed to analyze fluxes of energy (H, LE) and F CO 2 , and respective data quality assurance protocols are well established. We applied these protocols in our data analysis as described in Section 2.3 for the computation of all fluxes including those of O 3 . The resulting data gaps are reported in Table 1. Altogether, 12% of the time series showed data gaps from power outages at the site. Additionally, 5-13% of the data were flagged because they did not fulfill the requirements of stationarity and/or integral turbulence characteristics (QC_2). Further, 18-23% of the good-quality data were flagged because of low u* values. The analysis of co-spectra is an additional diagnostic procedure to evaluate the quality of flux estimates in terms of spectral characteristics of the turbulent transport and the instrumental setup. The dashed blue line in Figure 2 shows a theoretical co-spectrum model [49] of the vertical wind component with temperature, which represents H. Deviations from the modelled co-spectrum in the high-frequency range result from dampening of ambient fluctuations within the intake tube and the time lag resulting from sensor separation [53], while low-and mid-frequency deviations are rather caused by non-stationarities, trends, or by low-frequency fluctuations within the analyzed time series (e.g., [54]). High similarity of measured data with the modelled co-spectra thus indicates a good coverage of the ambient turbulence by the employed measurement setup. Figure 2 shows that the ensemble co-spectra w O 3 are in good agreement with the low-to mid-frequency range of the modelled co-spectrum, while the high-frequency range is dampened and becomes noise dominated at normalized frequencies above 2.3. On the other hand, the co-spectra w'T' were noise dominated only at normalized frequencies above 3.0, but not dampened ( Figure 2). Dampening effects of the ozone flux, as shown in Figure 2, are accounted for in each half-hour flux computation through spectral corrections as described in Section 2.3. The spectral correction factors had a median of 1.23, an inner quartile range of 0.13, and the 10th and 90th percentiles were 1.12 and 1.40.
coverage of the ambient turbulence by the employed measurement setup. Figure 2 shows that the ensemble co-spectra are in good agreement with the low-to mid-frequency range of the modelled co-spectrum, while the high-frequency range is dampened and becomes noise dominated at normalized frequencies above 2.3. On the other hand, the co-spectra w'T' were noise dominated only at normalized frequencies above 3.0, but not dampened (Figure 2). Dampening effects of the ozone flux, as shown in Figure 2, are accounted for in each half-hour flux computation through spectral corrections as described in section 2.3. The spectral correction factors had a median of 1.23, an inner quartile range of 0.13, and the 10 th and 90 th percentiles were 1.12 and 1.40.

General Meteorology and Flux Pattern
An eleven-day time series of selected meteorological, ecosystem flux, and ozone data is shown in Figure 3. Clear diurnal cycles are evident for all depicted parameters. The solar radiation (top panel) indicates that most days were partly cloudy. Tair reached maxima between 20.7 °C and 28.4 °C,

General Meteorology and Flux Pattern
An eleven-day time series of selected meteorological, ecosystem flux, and ozone data is shown in  Figure S1) for details). The slopes of the correlation lines revealed that the concentration measurements at MWM were about 25% to 30% lower than those at the LÜN stations. The nighttime data fit generally better between the sites than the peak concentrations during the days.
During the nights, the near-surface ozone decreased due to non-stomatal deposition (see Section 3.3 for more details) from a shallow, stable boundary layer to the surface, and potentially due to chemical reactions with gases such as NO. During the days, χ O 3 increased rapidly in the morning due to break-up of the stable nocturnal boundary layer (Figure 3 drop of z/L from positive to negative values) and down-mixing of ozone from the residual layer. Furthermore, rather slow production of additional ozone occurred during the course of the days [55].  The stomatal and non-stomatal deposition velocities are shown in the second panel from below ( Figure 3). The stomatal deposition velocity is driven by GPP (Equation (11)) and is thus zero at night. The non-stomatal deposition velocity is generally larger than v d,sto . On 1 August during the daytime, v d,sto is larger than v d,nsto , and v d,nsto even becomes slightly negative for individual half hours. Even though χ O 3 is rather high, with values around 50 ppb, F O 3 is low with daytime values of 4-6 nmol m −2 s −1 . The mismatch between GPP and F O 3 (low F O 3 at high GPP) caused v d,nsto to become negative. We could not find any evidence for instrumental malfunctioning that could have caused the low F O 3 . The wind direction was south-southwest, which was not common during our measurement period. Therefore, we can only assume that spatial heterogeneity in surface properties (i.e., much lower LAI or open water surfaces) likely caused the low F O 3 on this day.
The measured ozone fluxes are shown in the bottom panel of Figure 3. They were negative throughout with large deposition fluxes up to about −15 nmol m −2 s −1 occurring during the midday hours.

O 3 -Flux Partitioning
The turbulent deposition flux of ozone, F O 3 was partitioned into F sto and F nsto (Equation (6)) for each half hour. Figure 4 shows the median diurnal courses of these individual fluxes. F nsto is generally larger than F sto . Both F sto and F nsto follow a diurnal pattern. F sto is driven by the conductance of the stomata (and thus related to GPP) and by c O 3 . It is zero during the nights, and the total F O 3 is established by F nsto . F nsto is larger than F sto for most of the daytime as well, except during the afternoon hours, when F nsto and F sto are about equal to each other.

O3-Flux Partitioning
The turbulent deposition flux of ozone, was partitioned into and (Equation (6)) for each half hour. Figure 4 shows the median diurnal courses of these individual fluxes. is generally larger than . Both and follow a diurnal pattern. is driven by the conductance of the stomata (and thus related to GPP) and by . It is zero during the nights, and the total is established by . is larger than for most of the daytime as well, except during the afternoon hours, when and are about equal to each other. For situations when no evaporation is present, , (Equation (9)) can be used to calculate instead of (Equation (11)). Using this simplification leads to an overestimation of the amount of the nocturnal by about 1 nmol m −2 s −1 , (nocturnal ≈ −1 nmol m −2 s −1 ) and a root mean square error (RMSE) between the simplified and the correct approach of 1.3 nmol m −2 s −1 (details not shown in detail). This overestimate stems from nocturnal evaporation, of which a considerable portion was interpreted to be a stomatal flux when applying the simplification. For daytime conditions, the RMSE between the two methods was 1.7 nmol m −2 s −1 and as calculated with the simplified approach is mostly higher, which suggests that also during daytime, a certain amount of evaporation falsely contributed to the computed transpiration flux and thus to . Therefore, taking the GPP into account when calculating the stomatal conductance (Equation (11)) and thus leads to much better results.  For situations when no evaporation is present, g sto,PM (Equation (9)) can be used to calculate F sto instead of g sto (Equation (11)). Using this simplification leads to an overestimation of the amount of the nocturnal F sto by about 1 nmol m −2 s −1 , (nocturnal F sto ≈ −1 nmol m −2 s −1 ) and a root mean square error (RMSE) between the simplified and the correct approach of 1.3 nmol m −2 s −1 (details not shown in detail). This overestimate stems from nocturnal evaporation, of which a considerable portion was interpreted to be a stomatal flux when applying the simplification. For daytime conditions, the RMSE between the two methods was 1.7 nmol m −2 s −1 and F sto as calculated with the simplified approach is mostly higher, which suggests that also during daytime, a certain amount of evaporation falsely contributed to the computed transpiration flux and thus to F sto . Therefore, taking the GPP into account when calculating the stomatal conductance (Equation (11)) and thus F sto leads to much better results.
The partitioning of F O 3 into F sto and F nsto allows us to disentangle different processes driving F O 3 . Figure 5 shows a plot of the stomatal flux F sto versus GPP with χ O 3 used to color-code individual data points. As GPP and χ O 3 are the main drivers of F sto , a dependence of the slopes of regressions for various χ O 3 classes is to be expected in this plot. The larger χ O 3 is, the larger (with larger negative numbers) F sto is. However, there remains a scatter in Figure 5 due to the fact that it is not the ambient ozone concentration c O 3 that is used to compute F sto but the surface ozone concentration c O 3sur f (Equation (7)). This surface concentration is computed for each half hour and it introduces the variability in the computation of F sto from measured data. The partitioning of into and allows us to disentangle different processes driving . Figure 5 shows a plot of the stomatal flux versus GPP with used to color-code individual data points. As GPP and are the main drivers of , a dependence of the slopes of regressions for various classes is to be expected in this plot. The larger is, the larger (with larger negative numbers) is. However, there remains a scatter in Figure 5 due to the fact that it is not the ambient ozone concentration that is used to compute but the surface ozone concentration (Equation (7)). This surface concentration is computed for each half hour and it introduces the variability in the computation of from measured data. The more important non-stomatal flux is the residual of flux computations (Equation (6)). According to general understanding, it should be associated with the surface area available for deposition and its properties (e.g., LAI and surface wetness) and atmospheric turbulence (expressed as u*). Figure 6 presents the results of the driver analysis for , . Overall (lowest panel in Figure  6), an R 2 of about 26% was achieved, of which u* and Rg together account for 50%. The low R² value shows that a large portion of the variability of cannot be explained with the data and analytical tools employed. A graphical presentation of the relationships between , and environmental drivers ( Figure 6, top 6 panels) shows that most relationships are not linear and contain much scatter, which explains the low value of R² for the multiple linear model. The more important non-stomatal flux F nsto is the residual of flux computations (Equation (6)). According to general understanding, it should be associated with the surface area available for deposition and its properties (e.g., LAI and surface wetness) and atmospheric turbulence (expressed as u*). Figure 6 presents the results of the driver analysis for v d,nsto . Overall (lowest panel in Figure 6), an R 2 of about 26% was achieved, of which u* and R g together account for 50%. The low R 2 value shows that a large portion of the variability of F nsto cannot be explained with the data and analytical tools employed. A graphical presentation of the relationships between v d,nsto and environmental drivers ( Figure 6, top 6 panels) shows that most relationships are not linear and contain much scatter, which explains the low value of R 2 for the multiple linear model. Nevertheless, the observed relationships between , and environmental drivers at our site are in agreement with the work of Altimir et al. [56], Lamaud et al. [11], Fares et al. [24], Coyle et al. [23], and Hogg et al. [57]. Even though all these measurements were performed in different ecosystems, the respective results agree with ours in that the relationships between fluxes and drivers are not linear and that the scatter in these relationships is large.
Following a suggestion given by Fowler et al. and Cape et al. [21,58], we calculated the dependency of to thermal decomposition at the surface, and we found that our site had 25% higher activation energy as compared to the Auchencorth Moss in Fowler et al. [21].

O3-Flux Partitioning Uncertainties
Generally speaking, the partitioning of into and works well. Nevertheless, there are several sources of uncertainty that have not yet been reported in the literature. When using the Penman-Monteith approach to calculate the stomatal conductance (or stomatal resistance) for H2O (ep. 9), the molecular diffusivities of H2O and O3, of and , are employed. The values of this Nevertheless, the observed relationships between v d,nsto and environmental drivers at our site are in agreement with the work of Altimir et al. [56], Lamaud et al. [11], Fares et al. [24], Coyle et al. [23], and Hogg et al. [57]. Even though all these measurements were performed in different ecosystems, the respective results agree with ours in that the relationships between fluxes and drivers are not linear and that the scatter in these relationships is large.
Following a suggestion given by Fowler et al. and Cape et al. [21,58], we calculated the dependency of F nsto to thermal decomposition at the surface, and we found that our site had 25% higher activation energy as compared to the Auchencorth Moss in Fowler et al. [21].

O 3 -Flux Partitioning Uncertainties
Generally speaking, the partitioning of F O 3 into F sto and F nsto works well. Nevertheless, there are several sources of uncertainty that have not yet been reported in the literature. When using the Penman-Monteith approach to calculate the stomatal conductance (or stomatal resistance) for H 2 O (ep. 9), the molecular diffusivities of H 2 O and O 3 , of D H 2 O and D O 3 , are employed. The values of this ratio as employed by various authors [50,51,59] vary between 1.67 and 1.51, which translates to a variation of F sto by about 10%.
The exclusion of direct evaporation from plant surfaces and open water surfaces is stated, and Lamaud et al. [11] proposed a method to account for this issue by using the relationship between GPP and stomatal conductance under conditions with relative humidity below 60% (Equation (11)). However, the effect of soil evaporation can not be corrected for with this method. Most likely, soil evaporation can only be excluded if additional soil evaporation measurements are performed.
Finally, additional uncertainty originates from the GPP data which were used when applying the correction according to Equation (11). In our case, we separated the net ecosystem flux (NEE) of CO 2 into GPP and ecosystem respiration (R eco ) by estimating the daytime R eco from light response curves [47]. This yields a GPP of 0 µmol m −2 s −1 during the night. When using nighttime F CO 2 data instead and extrapolating these values to the daytime through a temperature relationship [46], the nocturnal GPP is not necessarily zero. Rather, it fluctuates around 0 and could eventually cause negative stomatal conductance during the night. Additionally, the daytime GPP values as calculated from these two methods would differ from each other as well.

Conclusions
Ozone deposition fluxes (F O 3 ) were measured during the peak of the growing season at a peatland in NW Germany. The employed eddy covariance setup for the O 3 fluxes and employing the gas phase chemiluminescence method after reaction of O 3 with NO showed good co-spectral characteristics even though the signal was damped in the mid-to high-frequency range, most likely caused by the intake tube.
The ozone fluxes were partitioned into stomatal (F sto ) and non-stomatal (F nsto ) fluxes by using resistance analogy. The use of the Penman-Monteith approach to estimate F sto led to an overestimate of the stomatal conductance and thus of v d,sto and F sto . The inclusion of GPP into the computational routine led to a clear improvement of the estimate of F sto . The stomatal flux is thus rather well understood and can thus be parameterized and modeled with the stomatal conductance and the ambient ozone concentration. For further measurements of F O 3 to vegetated surfaces, we recommend measuring F CO 2 and applying the corrections employing GPP to determine F sto . However, further research needs to be performed in order to quantify the effect of various routines to separate NEE into GPP and R eco on the flux partitioning of ozone.
The non-stomatal flux, which is computed as the residual F nsto = F O 3 − F sto (Equation (6)), was larger than F sto , is less well defined, and more difficult to parameterize or to model. Atmospheric turbulence plays the single most important role in non-stomatal ozone flux; however, the overall coefficient of determination between the flux and drivers was rather poor, i.e., only 26%. Further research on ozone deposition should focus on longer time series covering all seasons, on various types of ecosystems, and other types of surfaces in order to better understand the non-stomatal ozone deposition flux.