No impact of tropospheric ozone on the gross primary productivity of a Belgian pine forest

High stomatal ozone (O3) uptake has been shown to negatively affect crop yields and the growth of tree seedlings. However, little is known about the effect of O3 on the carbon uptake by mature forest trees. This study investigated the effect of high O3 events on gross primary productivity (GPP) for a Scots pine stand near Antwerp, Belgium over the period 1998–2013. Stomatal O3 fluxes were modelled using in situ O3 mixing ratio measurements and a multiplicative stomatal model, which was parameterised and validated for this Scots pine stand. Ozone-induced GPP reduction is most likely to occur during or shortly after days with high stomatal O3 uptake. Therefore, a GPP model within an artificial neural network was parameterised for days with low stomatal O3 uptake rates and used to simulate GPP during periods of high stomatal O3 uptake. Possible negative effects of high stomatal O3 uptake on GPP would then result in an overestimation of GPP by the model during or after high stomatal O3 uptake events. The O3 effects on GPP were linked to AOT40 and POD1. Although the critical levels for both indices were exceeded in every single year, no significant negative effects of O3 on GPP were found, and no correlations between GPP residuals and AOT40 and POD1 were found. Overall, we conclude that no O3 effects were detected on the carbon uptake by this Scots pine stand.


Introduction
Tropospheric ozone (O 3 ) is a secondary air pollutant that has the potential to negatively affect vegetation, leading to reduced growth and carbon sequestration potential (ICP Vegetation, 2011; Subramanian et al., 2015).Background concen-trations of tropospheric O 3 have increased 30 % since preindustrial times (Young et al., 2013) and are projected to further increase considerably until about 2050 (IPCC, 2007).Depending on the scenarios, background O 3 levels might either increase or decrease after 2050 (IPCC, 2007).
In recent years, many studies have been conducted to investigate the mechanisms underlying the O 3 impacts on vegetation.Ozone reduces plant growth by altering photosynthetic rates, carbohydrate production, carbon sequestration, carbon allocation, and carbon translocation (Beedlow et al., 2004;Ashmore, 2005;Wittig et al., 2009).Once O 3 enters the leaves through the stomata, it can affect plant growth through direct cellular damage (Mauzerall and Wang, 2001), leading to visible leaf injury and reduced leaf longevity (Li et al., 2016).In response to O 3 , respiratory processes increase, which will also affect the tree's carbon balance (Ainsworth et al., 2012).Skärby et al. (1987) proved that the dark respiration of Scots pine shoots increased after long-term exposure to a low level of O 3 .Protective responses, such as compensation (e. g. repair of injured tissue), avoidance (e. g. stomatal closure), and tolerance (e. g. alteration of metabolic pathways), all consume carbon; hence, resistance to O 3 damage costs energy.The size of this cost affects the amount of carbon remaining to support growth (Skärby et al., 1998).
To assess the impact of O 3 , several indices have been created: AOT40 (ppb h), the cumulative O 3 mixing ratio in excess of a threshold of 40 ppb, and POD y , the accumulated O 3 flux above a flux threshold y (nmol m −2 s −1 ).Critical levels are quantitative estimates of exposure to O 3 above which direct adverse effects may occur (CLRTAP, 2015); these have been determined for the indices based on O 3 dose-response relationships from fumigation experiments with enhanced O 3 mixing ratios (Karlsson et al., 2004).The magnitude of the O 3 impact on plants depends on the intensity of O 3 exposure, environmental factors influencing both plant photosynthesis and the O 3 flux to plant surfaces, and plant speciesspecific defensive mechanisms (Musselman and Massman, 1999).Because of the variable plant responses to similar O 3 mixing ratios, the question arises as to whether widely applicable tolerable limits of O 3 mixing ratios exist (Skärby et al., 1998).
While high stomatal O 3 fluxes have been shown to affect the yield of crops and the growth of tree seedlings and saplings (e.g.Büker et al., 2015), little is known about the effect on mature forest trees.When scaling up the results from seedlings to mature trees, the resulting data should be viewed with caution due to differences in energy budgets, canopy : root balances, and architecture and carbon allocation patterns (McLaughlin et al., 2007;Huttunen and Manninen, 2013).In addition to the uncertainties related to the up-scaling from seedlings to mature trees, data from controlled experiments should also be used with caution because trees can react differently in field conditions (Skärby et al., 1998).The effect of O 3 uptake on carbon uptake under ambient O 3 mixing ratios by trees has hardly been studied in situ.Some studies showed reductions in plant growth due to stomatal O 3 uptake (Zapletal et al., 2011;Fares et al., 2013;Yue and Unger, 2014), while other studies did not show any effect (Samuelson, 1994;Zona et al., 2014).The presence of a stomatal O 3 uptake effect was found to be species-and site-specific, and there is a clear need for more studies investigating the effect of O 3 on carbon uptake by mature trees in the field (Huttunen and Manninen, 2013).
In this study, we investigated the effect of O 3 at ambient levels on the gross primary productivity (GPP) of a mature Scots pine stand in Flanders, Belgium over a period of 14 growing seasons between 1998 and 2013.The investigation of O 3 effects on GPP is relevant because GPP represents the first step in the process of C assimilation and quantifies the rate at which C substrate is provided for growth, wood production, et cetera.Critical levels of AOT40 and POD 1 are being exceeded for this stand (Neirynck et al., 2012), indicating a potential effect of O 3 on tree productivity already at current ambient levels.To detect O 3 effects on GPP, we adopted a modelling approach that involved simulating GPP with a model with an O 3 -damage-free parameterisation and evaluating model overestimations of GPP.We used an artificial neural network (ANN) to model GPP.ANNs are a power tool to process multidimensional data in which complex nonlinear interrelationships between the parameters can be expected.ANNs are successfully used in remote sensing, evolutionary ecology, et cetera, and have previously been used to model GPP (Lek and Guegan, 1999;Rochelle-Newall et al., 2007;Akhand et al., 2016;Liu et al., 2016).In this study, we used ANNs since they do not employ predefined model conditions compared to conventional statistical models.

Study area
The study area consisted of a 2 ha Scots pine stand in a 150 ha coniferous and deciduous forest named "De Inslag", situated in Brasschaat (+51 • 18 33 N, +04 • 31 14 E) northeast of the Antwerp agglomeration and east-northeast of the Antwerp harbour (Neirynck et al., 2008).The site has a temperate maritime climate with a mean annual temperature of 11 • C and a mean annual precipitation of 830 mm (Neirynck et al., 2008).The soil has been classified as Albic Hypoluvic Arenosol (Gielen et al., 2011), a moderately wet sandy soil with a distinct humus and/or iron B horizon (Janssens et al., 1999).The sandy layer overlays a clay layer, which is situated at a depth of 0.7-2 m.As a result of the poor drainage, groundwater depth is typically high, fluctuating between 0.5 and 2 m (Carrara et al., 2003).
The pine stand was planted in 1929 (Neirynck et al., 2008).Until the autumn of 1999, when the forest was thinned, the tree density amounted to 542 trees ha −1 .The thinning decreased the tree density to 376 trees ha −1 .The average canopy height is 21.4 m (Op de Beeck et al., 2010).With a peak in leaf area index (LAI) of 1.3 ± 0.5 m 2 m −2 in 2007 (Op de Beeck et al., 2010) and an average LAI of 1.2 ± 0.5 m 2 m −2 in the period 1998-2007, the stand canopy is very sparse.Only two needle-age classes are present: current-year needles and 1-year-old needles (Op de Beeck et al., 2010).
The stand is part of the ICP Forests Level II and the Fluxnet CarboEurope-IP networks and is equipped with a 41 m tall instrumentation tower.Meteorological measurements and measurements of ecosystem CO 2 exchange with the eddy covariance technique have been conducted at the site on a continuous basis since 1996 (Gielen et al., 2013).

Measurements
The study covered the period 1998-2013, with the years 1999 and 2003 excluded due to poor data quality or coverage.

Meteorology
Air temperature (T air ; • C) and humidity (RH; %) were measured with a PT100 and a HMP230 dew point transmitter (Vaisala, Helsinki, Finland) in aspirated radiation shields mounted on the tower at 2, 24, and 40 m of height.Wind speed (WS, m s −1 ) was measured with a cup anemometer (LISA; Siggelkow GmbH, Hamburg, Germany) at 24, 32, and 40 m of height.Ingoing and outgoing shortwave and long-wave radiation were measured at the top of the tower with a CNR1 radiometer and a CMP6 pyranometer (Kipp and Zonen, Delft, The Netherlands).Rainfall was registered by a tipping bucket rain gauge (NINA precipitation pulse transmitter; Siggelkow GmbH, Hamburg, Germany).Both T air and RH were used to calculate the vapour pressure deficit (VPD; kPa).The soil temperature (T soil ; • C) was measured at 9 cm below the soil surface with temperature probes (DPS-404; Oxsensis, Didcot, UK).The soil water content (SWC; m 3 m −3 ) was measured at 25 cm below the soil surface with time domain reflectometers (CS616; Campbell Scientific, Logan, Utah, USA).Instant SWC was read manually from the reflectometers every 3 to 14 days and values were interpolated to obtain daily estimates taking into account water inputs via precipitation (Gielen et al., 2010).The soil water potential (SWP; MPa) was derived from the SWC measurements with the model of van Genuchten (van Genuchten, 1980).All meteorological variables (except SWC and rainfall) were measured every 10 s and half-hourly means were calculated.Data gaps were filled with data from nearby weather stations.

Ozone mixing ratio
The O 3 mixing ratio ([O 3 ]; ppb) was measured at a 10 s resolution above the canopy at 24 m of height with a UV photometric analyzer (TEI 49i; Thermo Fisher Scientific, Waltham, MA, USA) and converted to half-hourly averages.Data gaps were filled with [O 3 ] measurements done at 40 m of height.If these were not available, gaps were filled with [O 3 ] measurements from a nearby weather station from the Flemish Environmental Agency (VMM) at Luchtbal, which is less than 10 km from the site.

Leaf area index
A continuous time series with daily LAI values was reconstructed for the pine stand based on the historical data.The general approach was to keep the seasonal pattern measured in 2009 by Op de Beeck et al. (2010) fixed for each year and to scale it year per year to the seasonal maximum LAI (LAI max ).LAI max was measured with the LAI-2050 (LI-COR, Lincoln, Nebraska, USA) in 1997 and 2003 by Gond et al. (1999) and Konôpka et al. (2005), respectively, and with digital hemispherical photography in 2007 by Op de Beeck et al. (2010).To ensure consistency across the time series, the measurements were corrected for clumping using a factor of 0.83 (Jonckheere et al., 2005).The three measurements of LAI max were interpolated linearly to derive LAI max values for the missing years.The thinning event in 1999 was accounted for by subtracting the removed leaf biomass, determined with allometric relations from Yuste et al. (2005) and specific leaf area measurements from Op de Beeck et al. (2010).

Gross primary productivity
Gross primary productivity (µmol C m −2 s −1 ) was derived from net ecosystem exchange (NEE) measured with the eddy covariance technique and following the standard data quality procedures as explained in Appendix A. Half-hourly aver-aged values of GPP were derived for the 14 growing seasons of the study period and integrated into daily and growing season totals.

Stomatal conductance
The measurements of stomatal conductance to H 2 O (g st,H 2 O ) were done at needle level during the summers of 2007 (Op de Beeck et al., 2010) and 2013 to obtain data for the parameterisation of the multiplicative stomatal model used in the calculation of stomatal O 3 fluxes (see Sects.2.3 and 2.4).The two summers were marked by quite different environmental conditions: cold and wet in 2007 and warm and dry in 2013.Measurements were carried out with the LI-6400 gas exchange system (LI-COR, Lincoln, Nebraska, USA) and included diurnal stomatal courses as well as stomatal responses to PAR, T air , and VPD.Measurements were carried out on sets of three or four live fascicles, i.e. six to eight needles, which were enclosed in the LI-6400's leaf chamber while attached to the tree.Twenty-six needle sets were measured in total, equally divided between current-year and 1-year-old needles.Each needle set was harvested after being measured, and the hemi-surface needle area was determined in order to express g st,H 2 O on the correct needle area basis.Needle area was derived from needle dimensions (length and width at the top, middle, and base), assuming a hemi-circular crosssectional needle area.The measurements of g st,H 2 O were converted to stomatal conductance to O 3 (g st ) by multiplying g st,H 2 O with the ratio of the molecular diffusivities of water vapour and O 3 in the air (0.61).

Calculation of stomatal O 3 fluxes
The stomatal O 3 fluxes were calculated at a half-hourly resolution from a continuous series of half-hourly [O 3 ] meteorology and daily LAI with an electric analogue model built from three resistances in series: where R tot is the total resistance to O 3 , R aero is the aerodynamic resistance to O 3 , R bl is the quasi-laminar boundary layer resistance to O 3 , and R can is the canopy resistance to O 3 (all expressed in s m −1 ).
The aerodynamic resistance was calculated following Grünhage (2002) with where κ is the von Karman constant (0.43), u * (m s −1 ) is the friction velocity, L is the Obukhov length, z is the [O 3 ] measurement height (24 m), d is the zero plane displacement (= 0.1 h), z 0 is the momentum roughness parameter (= 0.65 h), h is the canopy height, and h is the atmospheric stability function.This function is calculated using the set of coefficients published by Dyer (1974): L. T. Verryckt et al.: No impact of tropospheric ozone on the GPP of a Belgian pine forest for unstable atmospheric stratification (L < 0 m) for stable atmospheric stratification (L > 0 m) and for neutral atmospheric stratification (|L| → ∞).
The quasi-laminar boundary layer resistance was calculated following Baldocchi et al. (1987) with where κ is the von Karman constant (0.43), u * (m s −1 ) is the friction velocity, which is derived from the measured momentum fluxes, Sc is the Schmidt number (1.07 for O 3 ), and Pr is the Prandtl number (0.72 for O 3 ).
The canopy resistance was calculated from a stomatal resistance (R st ) and a non-stomatal resistance (R nst ) mounted in parallel: The stomatal resistance R st was calculated with an algorithm that divides the pine canopy into eight horizontal leaf layers, with the LAI being divided equally between the layers, that simulates the transfer of radiation through the layered canopy.The algorithm then calculates the stomatal resistance for the sunlit and shaded area fraction of each leaf layer with the multiplicative stomatal model described by Jarvis (1976) and reformulated by Emberson et al. (2000).
Resistance values are then integrated over all layers to obtain canopy level R st .The algorithm is explained in more detail in Op de Beeck et al. (2010).The version of the multiplicative stomatal model used in this study is described in detail in Appendix B. This model was given a site-specific parameterisation as explained in Sect.2.4.The non-stomatal resistance R nst was assumed to be constant in time and set to 279 s m −1 .This value was derived from the long-term O 3 flux measurements in Brasschaat (Neirynck et al., 2012).
Total and stomatal O 3 fluxes (F tot and F st ; nmol m −2 s −1 ) were calculated on a half-hourly basis with where 44.64 is the molar density of air in mol m −3 at an air pressure of 101.3 kPa and an air temperature of 0 • C, used here to convert flux units from m s −1 to mol m −2 s −1 .Halfhourly fluxes were aggregated to daily and yearly values.

Parameterisation and validation of the multiplicative stomatal model
The multiplicative stomatal model was parameterised and validated against the dataset of g st measurements collected at the site.This dataset included, in addition to the measured g st , also PAR, T air , VPD, and SWP and was split into a parameterisation set and a validation set by grouping the odd and even rows of data after being ranked by PAR.The parameterisation was done by optimising the model parameters with the function "lsqcurvefit" in Matlab (Matlab, Natick, MA, USA; Statistics Toolbox Release, 2013a), which finds the best parameter values starting from an initial value and which can be used to fit non-linear functions with more than two independent variables.The parameters of the boundary functions f PAR , f T air , f VPD , and f SWP were optimised separately, starting from initial values that were estimated visually from plots of g st versus each of the input variables (PAR, T air , VPD, and SWP).The phenology function f phen was set to 1 for the parameterisation of f PAR , f T air , f VPD , and f SWP since g st had been measured on mature needles only.We included f phen in the final model to estimate the stomatal O 3 fluxes over the growing season (Appendix B).The parameterised model was then tested against the validation dataset.The model performance was evaluated with the linear regression y = ax + b fitted to the plot of the measured versus the modelled g st and with the following set of performance statistics: the coefficient of determination (R 2 ), mean bias (MB), relative mean error (RME), Willmott's index of agreement (d), model efficiency (ME), and root mean square error (RMSE) and its systematic (RMSE s ) and unsystematic (RMSE u ) components.These statistics are explained briefly in Appendix C. To visually evaluate the goodnessof-fit of each boundary function, the modelled g st was plotted versus each of the input variables and the corresponding boundary function added to the scatter plot.

Detecting O 3 effects on GPP
We adopted a modelling approach to detect possible O 3 effects on GPP.Under the assumption that O 3 -induced GPP reduction is most likely to occur during and shortly after days of high stomatal O 3 fluxes, we parameterised a GPP model against a dataset from which such days where removed and then we simulated the daily and growing season GPP with this supposedly O 3 -damage-free model.A reduction in GPP due to O 3 would become apparent as a model overestimation of daily GPP for the days on which an O 3 effect was assumed, and possibly also as an overestimation of growing season GPP.The physiological mechanism beyond the assumption made above is that at high stomatal O 3 fluxes, the trees' defensive mechanisms cannot detoxify all O 3 entering the needles and damage is caused to the photosynthetic apparatus (Dizengremel, 2001;Matyssek and Sandermann, 2003).This leads to decreased gross photosynthetic rates and GPP.The damage is repaired afterwards when the stomatal O 3 load decreases.
We used a feed-forward back propagation artificial neural network (ANN) as a GPP model in Matlab (Matlab, Natick, MA, USA).The ANN contained 10 nodes organised in 1 layer, which came out as the best performing network after comparing networks containing different numbers of nodes and/or layers (data not shown).The default settings of the Matlab Neural Network Toolbox were used.A normalisation process was applied for training and testing the data; data were scaled to [−1 1] based on the lowest and highest value in the dataset.We used the Levenberg-Marquardt algorithm to train the ANN for 1000 iterations (Marquardt, 1963).The progress of the training procedure was monitored using the mean square error (MSE) of the network.The daily GPP data were used as dependent target variables in the ANN.The input variables were year, day of the year, T min , T max , T mean , average VPD, SWC, R g , average T soil , and average WS.The daily totals of the variables were used, with the exception of VPD, T soil , and WS, for which the daily averaged values were used.The individual weights of these parameters on our model were estimated by replacing each input variable with a random permutation of its values.This was done for the GPP model, as described above, and a GPP model containing O 3 as an input variable to test whether O 3 had any explanatory power on GPP.
To obtain an O 3 -damage-free GPP model, the days for which an O 3 effect on GPP was expected were removed from the dataset.We assumed that if an O 3 effect occurs, it would occur on the days with the highest stomatal O 3 fluxes.Because the defensive capacity of the pine trees was not quantified, and hence the O 3 load above which O 3 would affect GPP was not known, we repeated the analysis trice by removing the days with the 2, 5, and 10 % highest stomatal O 3 fluxes.Because the results for a 2 and 10 % cut-off were equal to those for a 5 % cut-off, we report only the results for a 5 % cut-off.The model was trained with two-thirds of the remaining dataset, while the other one-third was used to test the model.This O 3 -damage-free model was then run with the full dataset.
The model overestimation of daily GPP was evaluated (1) from the linear regression on the data of the measured versus the modelled GPP for the days on which an O 3 effect was assumed, testing whether the regression slope and intercept were different from 1 and 0, respectively, and (2) by comparing the measured and the modelled daily GPP for these days by means of a paired-samples t test or a Wilcoxon signed-rank test if differences were not normally distributed (Shapiro-Wilk test).A significant outcome of this test in combination with a regression slope significantly lower than 1 (and an intercept not different from 0) would together point to a significant overestimation of GPP.Furthermore, (3) the regression slope and intercept were compared with the slope and intercept of the regression fitted to the dataset used to train and test the GPP model.This was done to evaluate whether GPP estimations for the days on which we assumed an O 3 effect were, in relative terms, significantly higher than GPP estimations for the days used for model training and testing.This would become apparent as a significantly lower slope (with an intercept no different from 0).The model overestimation of the growing season GPP was evaluated with the first two tests above on the growing season data.Additionally, the residuals of growing season GPP (model -measurement) were plotted against AOT40, POD 1 , and total growing season stomatal O 3 uptake and linear regression lines were fitted.It was tested whether the regression slope and intercept were significantly different from 0 to assess the presence of a statistically significant O 3 dose-response relationship.
Since it may take some time for trees to repair the damage to the photosynthetic apparatus induced by O 3 , O 3 effects might last several days after a peak of O 3 exposure.They might thus not be detected with the model parameterised as explained above.To account for such a sustained O 3 effect, the modelling was repeated, now not only excluding the days with the highest stomatal O 3 fluxes from the dataset for model training but also the following days.The modelling was repeated with three different such delay periods: the first, the first 2, and the first 6 days following each flux peak.The results were evaluated with the same statistical tests as mentioned above.Because the results were similar for the three delay periods, only the results for the 2-day period are shown.
High O 3 events are often coupled with specific meteorological conditions, i.e. high radiation and air temperatures.Since the dataset for the model training had been compiled by removing the days with the highest stomatal O 3 fluxes, it was not unlikely that these conditions were underrepresented in the training dataset.If so, this could induce a bias in the model response to radiation and temperature and possibly result in overestimations of GPP for the days on which an O 3 effect was expected, which we then might wrongly attribute to O 3 .To evaluate the risk for such model bias, we compared the frequency distribution, range of radiation, T min , T max , T mean , and VPD between the training dataset and the dataset with the days on which we expected an O 3 effect.
One of the assumptions in our approach is that the O 3 effects on GPP are short term, i.e. they last just a few days, and are hence not carried over.The presence of a carry-over effect would compromise the validity of our approach.We  can rule out a carry-over effect by testing whether trees exposed to low stomatal O 3 fluxes late in the growing season behave in the same way as when exposed to similarly low O 3 fluxes early in the growing season.To test this, we compiled a dataset that contained only the days after the first major peak of stomatal O 3 flux in the growing season.From this period, we further selected only the days with low stomatal O 3 fluxes for which no short-term O 3 effect was expected.In other words, we excluded the days with a peak of stomatal O 3 flux plus the 6 following days.We trained the GPP model with these data and then predicted GPP for the days before the first major O 3 peak in each growing season.If a carryover effect was present, or at least an effect induced during the first major O 3 flux peak, it would be somehow included in the trained model.This would then underestimate GPP for the days before each first major O 3 peak, where a carry-over effect has assumptively not yet occurred.The model underestimation of GPP was evaluated from a linear regression on the data of the measured versus the modelled GPP, testing whether the regression slope and intercept were different from 1 and 0, respectively.This slope and this intercept were also compared with the slope and intercept of the regression line fitted to the training data.Also, the measured and modelled GPP were compared with a paired-samples t test or a Wilcoxon signed-rank test whether differences were not normally distributed (Shapiro-Wilk test).
All statistics were performed with R version 3.2.3(The R Project for Statistical Computing Core Team, Vienna, Austria 2015) at a significance level of p = 0.05.

Measurements: meteorology, GPP, and LAI
Figure 1 shows a fingerprint of the multi-annual average diel and seasonal patterns of the main meteorological variables T air , incoming global radiation (R g ), and VPD as well as measured GPP.This figure gives a good overview of how meteorology and GPP typically changed over time in this forest; inter-annual anomalies from the average patterns can be found in Fig. S1 in the Supplement.Distinct daily and seasonal patterns can be observed; for example, reaching the highest values in summer in the afternoon.Similar patterns can also be observed in GPP, which basically follows the pattern of R g .As seen in Fig. 1, the photosynthetic period extends, on average, from day of the year 115 (the end of April) till day of the year 300 (the end of October).The time series of precipitation and SWP are provided in Fig. 2, while the seasonal LAI courses are shown for each year in Fig. 3.The yearly maximum LAI ranged from 1.4 to 1.9 m 2 m −2 .The thinning of the forest in 1999 can clearly be observed in the LAI pattern.After the thinning, the canopy never fully closed.

Multiplicative stomatal model and simulated O 3 fluxes
The optimised parameter values of the model are presented in Table 1.The different statistics to evaluate the model performance are presented in Table 2, and this is for both the parameterisation and the validation dataset.For the parameterisation dataset, the measured data were plotted against modelled g st and plotted in Fig. 4a.The slope of the linear fit was not significantly different from 1 (p = 0.87) and the intercept was not significantly different from 0 (p = 0.81).
The model evaluation for the validation dataset was equally good as for the parameterisation dataset (Table 2).Also, in the linear fit for the validation set (Fig. 4, b), the slope was  not significantly different from 1 (p = 0.98) and the intercept was not significantly different from 0 (p = 0.70).
Figure 5 shows the scatter plots of the measured g st versus each of the model input variables PAR, T air , VPD, and SWP and the fitted boundary function for each plot.
The average daily O 3 fluxes for the different years are presented in Fig. S2.Daily F st ranges from 1.12 to 1.52 nmol O 3 m −2 day −1 .In 2011, the daily F st was the lowest, while the highest values were observed in 2002.The annual average ratio F st /F tot varied between 24 and 28 % (Fig. S2).We observed the lowest ratios at the beginning and at the end of the growing season.Above-average ratios were observed at the peak of the growing season.

Ozone effects on GPP
Figure 6 shows the frequency distributions of R g , T min , T max , T mean , and VPD for the training dataset and the dataset with days on which we assumed an O 3 effect.Days in the latter dataset are generally more concentrated in the upper half of each variable's range.The training dataset includes more days in the lower half, but conditions of high radiation, tem- perature, or VPD do not seem to be underrepresented as the dataset also included a substantial number of days in the higher part.For all variables, the variable range of the dataset with days for which we assumed an O 3 effect is a fully contained range of the training dataset.
All parameters in the GPP model were ranked according to their contribution to GPP prediction (Table 3).Global radiation is the most important parameter in defining GPP, with a mean square error (MSE) of 37 500.81mol m −2 s −1 , followed by day of the year (30 240.61 mol m −2 s −1 ) and year (27 486.63 mol m −2 s −1 ).The maximum air temperature and VPD contribute equally to the model with an MSE of about 15 300 mol m −2 s −1 .Wind velocity, T min , and SWC contribute the least to GPP.Ozone as an input variable had an MSE of 11 885.73 mol m −2 s −1 (Table 3b) and contributed the least with an MSE similar to the overall model (10 019.30mol m −2 s −1 ).
To test for carry-over O 3 effects, we evaluated and compared the linear regressions of the measured versus the modelled GPP of a dataset with low O 3 fluxes after the first ma-  jor O 3 flux peak in the growing season and a dataset before this peak (Fig. 7).For both regressions, the intercept and slope were not significantly different from 0 and 1 respectively (training: p slope = 1, p intercept = 1; testing: p slope = 0.83, p intercept = 0.44).The slopes were also not significantly different from each other (p = 0.86), and neither were the intercepts (p = 0.53).
Figure 8 shows the measured versus the modelled daily GPP for the model trained without the days with the highest stomatal O 3 fluxes (GPP model 1) and the model trained to also test for lag effects (GPP model 2).Both models reproduced daily GPP well for the dataset against which they were trained and tested, as indicated by the high R 2 values and the fitted regression lines falling on the 1 : 1 line (Fig. 8a, b).For both models, the regression slope for the dataset with the days on which we assumed an O 3 effect was significantly lower than 1 and the intercept significantly higher than 0 (Fig. 8c, d).For GPP model 1, the regression slopes were not significantly different between the two datasets (p = 0.46), but the intercepts were (p < 0.05).For GPP model 2, both the regression slopes and intercepts differed significantly (p < 0.001 and p < 0.001).However, a Wilcoxon signed-rank test showed for both models that the modelled daily GPP was not significantly higher than the measured daily GPP for the days on which an O 3 effect was assumed (p = 0.83 and p = 0.64, respectively).Also, a paired-samples t test showed for both models that the modelled growing season GPP was not significantly higher than the measured growing season GPP (p = 0.93 and p = 0.55, respectively).The slope and intercept of the linear regression line were not significantly different from 1 and 0 (Fig. 8e, f).
No statistically significant correlations were found between the model residuals of growing season GPP and total stomatal O 3 uptake (F st ), AOT40, and POD 1 (Fig. 9). Figure S4 shows the relation between the measured growing season GPP, F st , AOT40, and POD 1 , as well as the time series of these parameters.No statistically significant correlations were found between the measured growing season GPP and AOT40 and POD 1 .A significantly negative correlation was found between the measured growing season GPP and F st (p = 0.006) due to a decline in GPP at a low F st (F st < 70 mmol O 3 m −2 ).  4 Discussion

Multiplicative stomatal model
All statistics shown in Table 2 clearly indicated that the fitted multiplicative stomatal model performed well.For both the parameterisation and the validation datasets, the model explained 72 % of the variance in g st .For both datasets, the slope and intercept of the linear regression lines of the measured versus the modelled g st were not significantly different from 1 and 0, respectively (Fig. 4).Moreover, the model efficiency (ME in Table 2) of 0.72 and the Wilmott's index (d) close to 1 both indicate that the modelled values matched the measured values well.A good model provides low root mean square error (RMSE), while the systematic component (RMSE s ) should approach zero and the unsystematic component (RMSE u ) should approach the RMSE (Willmott et al., 1985), which was the case for this model.A low mean bias (MB) and a low mean relative error (MRE) further indicated a very good performance.The good performance of the model can also be observed in Fig. 5, in which the boundary lines represent the response of g st to the independent variables when other variables were not limiting.The boundary lines fitted close to the data points, which is an indication of a good model because the multiplicative stomatal model is based on the assumption that the variables act more or less multiplicatively and independently from each other (Grüters et al., 1995).
As explained in the mapping manual of the Convention on Long-Range Transboundary Air Pollution (CLRTAP), Scots pine is the representative species to assess the risk of O 3 damage to coniferous forests in Atlantic central Europe (CLRTAP, 2015).This risk is assessed on the basis of O 3 doses calculated with the DO 3 SE algorithm, which employs a Jarvis-type stomatal model that has been parameterised for Scots pine based on a compilation of primary and secondary data (Emberson et al., 2007;Büker et al., 2015;CLRTAP, 2015).The parameterisation for our Scots pine stand differs in some numbers from the one used in the DO 3 SE algorithm.The most remarkable difference is that the g max of the Scots pines in Brasschaat is much lower (0.14 vs. 0.18 mol O 3 m −2 s −1 ).This low g max may imply that, during episodes of high O 3 mixing ratio, the Brasschaat site is unlikely to take up very high amounts of O 3 (Altimir et al., 2004;Emberson et al., 2007).This may have contributed to the absence of a clear O 3 response at our site.A second difference is that the stomata of the pine trees re- main opened at night (g min = 0.02 mol O 3 m −2 s −1 ), while the DO 3 SE model simulates full stomatal closure.Furthermore, the response to temperature for our Scots pine stand is shifted slightly higher (T opt = 25 vs. 20 • C) and the response to soil drought is much stronger (SWC max = −0.19 vs. −0.7 MPa and SWC min = −1.18 vs. −1.5 MPa).From these differences, it can be inferred that stomatal O 3 uptake rates at the Brasschaat site are considerably lower than would be simulated with the DO 3 SE model for generic Scots pine.This highlights the importance of a site-specific parameterisation when aiming to assess stomatal O 3 loads at site level.

Stomatal O 3 fluxes
The stomatal O 3 flux contributed on average 26 % to the total O 3 flux over the study period (Fig. S2).This fraction is similar to the 21 % stomatal O 3 flux in a Danish Norway spruce stand (Mikkelsen et al., 2004) and the 30 % stomatal O 3 flux in Quercus ilex in Italy (Vitale et al., 2005;Gerosa et al., 2005).Cieslik (2004) showed that in southern Europe, the stomatal O 3 flux of different vegetation types, such as pine forests and Mediterranean shrubs, is typically less than 50 % of the total O 3 flux.A 5-year study on a Mediterranean Pinus ponderosa stand showed a stomatal O 3 flux contribution of 57 % (Fares et al., 2010).Clearly, species-and site-specific differences, such as tree age or micro-climate, are introducing large variability in stomatal O 3 uptake (Neirynck et al., 2012).
The low relative stomatal O 3 flux in the Scots pine stand in Brasschaat could be the result of the sparse canopy with low LAI.Although no relation between stomatal O 3 flux and LAI was found in a previous study on this site (Neirynck et al., 2012), inter-annual and seasonal variation in LAI is very small, rendering such a correlation analysis very difficult.

Ozone effects on GPP
A comparison of the frequency distributions of radiation, temperature, and VPD between the training dataset and the dataset with the days on which we expected an O 3 effect showed that the meteorological conditions in the latter dataset were fully represented in the training dataset.From the full overlap, we can rather safely assume that the GPP model did not include a biased response to these variables that could result in a GPP overestimation that we might wrongly interpret as an effect of O 3 .Also, O 3 as an input variable in the ANN did not have any explanatory power on GPP as it had the lowest MSE value close to the overall model MSE.Furthermore, a GPP model parameterised to include a carry-over effect of O 3 on GPP did not overestimate GPP at a statistically detectable level for days on which such an effect was not assumed to occur.From these results, we infer that carry-over effects of O 3 were unlikely to have occurred and that the assumption of the absence of (detectable) carry-over effects was valid.
The statistical tests on the datasets of the measured and the modelled GPP did not reveal a statistically significant model overestimation of daily GPP for the days on which we assumed an O 3 effect nor an overestimation of growing season GPP.Also, no significant correlations were found between growing season GPP residuals and stomatal O 3 flux, AOT40, and POD 1 , even though critical levels for AOT40 and POD 1 were exceeded in every single year of our study period.From these results and within the limits of the modelling approach applied in this study, we can infer that no significant effect of O 3 on GPP occurred.Some earlier studies have investigated the effect of O 3 on forest carbon uptake.A cumulative stomatal uptake of 27 mmol m −2 over the growing season did not result in any visible damage or a reduction in NEE on a poplar plantation in Belgium (Zona et al., 2014).Zapletal et al. (2011), on the other hand, reported that the CO 2 uptake of a Norway spruce forest in the Czech Republic increased with increasing stomatal O 3 flux followed by a sudden decrease in CO 2 uptake, suggesting that an O 3 flux threshold exists.Fares et al. (2013) showed a negative correlation between GPP and O 3 uptake in two Mediterranean ecosystems (a forest dominated by Pinus ponderosa and an orchard of cultivated Citrus sinensis, both in California, USA).A GPP reduction of 1-16 % in response to O 3 uptake under an ambient O 3 mixing ratio of 30-50 ppb was determined across vegetation types and environmental conditions in the United States by Yue and Unger (2014).The magnitude of the reduction depended on the sensitivity to O 3 of the species and on the biome types.
AOT40 is, at present, the European standard for forest protection (EEA, 2014) with a critical level of 5000 ppb h, equivalent to a growth reduction of 5 % (CLRTAP, 2015).In this study on Scots pine in Brasschaat, this value was far exceeded in all years (Fig. 9), yet no negative effect on GPP was observed in years with higher AOT40 values.
POD 1 is considered a more appropriate index for potential O 3 damage because it considers O 3 flux.The critical level of POD 1 is species-specific; a critical level of 8 mmol m −2 with a 2 % growth reduction is used for Norway spruce and a critical level of 4 mmol m −2 with a 4 % growth reduction is used for birch and beech (CLRTAP, 2015).A critical level for Scots pine has not yet been determined, and therefore the value of 8 mmol m −2 for Norway spruce is often adopted as a critical level for Scots pine.During our 14-year study period, this critical level was exceeded every single year, and again no significantly negative correlation between total GPP residuals and POD 1 was observed.In comparison to the AOT40 level, 2006 was not the year with the highest POD 1 .This difference between AOT40 and POD 1 in 2006 was due to stomatal closure; during high O 3 mixing ratio events, g st was rather low (Fig. S3).POD 1 was highest in the year 2002, when O 3 mixing ratios were relatively low, but g st was high.The low O 3 mixing ratios explain the lower AOT40 for 2002.
Notwithstanding the absence of a statistically significant positive correlation between GPP residuals and both AOT40 and POD 1 , critical levels for both AOT40 and POD 1 were exceeded every single year.AOT40 is based on O 3 mixing ratios, and these concentration-based indices have been shown to be weaker indicators of O 3 damage than flux-based indices (Karlsson et al., 2007;Simpson et al., 2007).The critical level of POD 1 for Scots pine was adopted from the critical level for Norway spruce (CLRTAP, 2015).Possibly, this critical level is too low for Scots pine.
Figure S4 shows a negative relationship between the measured growing season GPP and the O 3 dose, most notably and only significant between GPP and F st (Fig. S4a).These trends suggest a strong effect of O 3 on GPP, which contradicts the outcome of our modelling analysis.The relationships are negative because the steady GPP increase that can be observed from the year 2005 until the end of the study period coincides with a steady decrease in O 3 loads (Fig. S4d, e, f).To our judgement, this GPP increase is more likely to be the result of forest regrowth in response to decreased acidification at the site (Neirynck et al., 2008) than a response to decreased O 3 loads.This forest regrowth is accounted for in our modelling analyses by means of the LAI input in the ANN, allowing us to disentangle the effect of both factors on GPP.We furthermore believe that the observed trends do not reflect a causal relationship between GPP reduction and O 3 loads because the GPP decrease is the strongest at low O 3 loads but virtually absent at high O 3 loads (Fig. S4a, b).This is not what would be expected in the case of an O 3 effect.
Overall, no significant O 3 effects on daily and growing season GPP were found with our modelling approach.It can thus be concluded that O 3 did not affect the GPP of the pine forest, at least not if the assumptions we made in our approach to detect O 3 effects are valid.The most crucial assumption involves the distinction between days at which a GPP effect did and did not occur.It was not possible to identify these days with great precision due to lack of knowledge on the defensive capacity of the trees and their ability to repair O 3 damage.To overcome this, we repeated our analysis with three different peak thresholds for daily stomatal O 3 uptake rates above which an effect would occur and with three different delay periods over which an induced O 3 effect would last.The fact that all nine analyses produced the same outcome provides validity for our conclusions, despite the uncertainty involved in the identification of days with O 3 effects.
The lack of detected O 3 effects on GPP does not mean that O 3 did not negatively affect this Scots pine stand in Brasschaat.Stomatal O 3 uptake has been linked here to reductions in GPP only.As already stated in the Introduction, protective responses such as compensation and enhanced tolerance occur in trees (Skärby et al., 1998).It is likely that the trees at our study site were able to fully detoxify the O 3 taken up.The respiratory cost involved might have come at the expense of biomass production and growth, while gross C uptake remained unaffected.Future analyses, such as a tree ring analysis, may provide an answer to whether this is the case.

Summary
We parameterised a multiplicative stomatal model for a Scots pine stand in Brasschaat.This species-and site-specific parameterised model performed very well.With this model embedded in a resistance scheme, stomatal O 3 fluxes were calculated and used to test for O 3 effects on GPP.Only very small reductions in growing season GPP were calculated.Although the critical levels for AOT40 and POD 1 were exceeded in every single year, no significant correlations between total GPP residuals and stomatal O 3 flux, AOT40, and POD 1 were found.Within the limitations of the approach used in this study, we can thus conclude that O 3 did not affect the gross carbon uptake by the Scots pine stand in Brasschaat.

Figure 1 .
Figure 1.The fingerprint of air temperature (T air ), incoming global radiation (R g ), vapour pressure deficit (VPD), and measured gross primary productivity (GPP) averaged over the period 1998-2013.The day of the year is plotted on the y axis and the hour of the day on the x axis.

Figure 2 .
Figure 2. The time series of the weekly total precipitation and mean soil water potential (SWP).The precipitation and SWP data are averaged over the period 1998-2013.The error bars represent the 95 % confidence intervals.

Figure 3 .
Figure 3.The seasonal course of the LAI for each of the 14 growing seasons used in this study.

Figure 4 .
Figure 4.The measured versus the modelled stomatal conductance (g st ) for the parameterisation dataset (a) (n = 205) and the validation dataset (b) (n = 205).The black line is the 1 : 1 line.The red line is the linear fit for which the equation is given in the figure.Also shown are the p values of the tests for the slope being different from 1 (p a ) and the intercept different from 0 (p b ).

Figure 5 .
Figure 5.The measured stomatal conductance (g st ) as a function of the different variables used in the multiplicative model: photosynthetically active radiation (PAR), air temperature (T air ), vapour pressure deficit (VPD), and soil water potential (SWP).The red line represents the boundary line for which the functions are given in Appendix B (Eqs.B3-B6) (n = 205).

Figure 6 .
Figure 6.Histograms of the meteorological variables for the training dataset (red) and the high O 3 uptake dataset (blue).The subplots represent global radiation R g (a), minimum temperature T min (b), maximum temperature T max (c), mean temperature T mean (d), and vapour pressure deficit VPD (e).

Figure 7 .
Figure 7.The measured GPP plotted as a function of the modelled GPP for two different datasets: (a) only the days before the first major O 3 peak in every year and (b) the training dataset with the days after the first major O 3 peak in every year, excluding those with high O 3 fluxes plus 6 following days to train the network.The black line is the 1 : 1 line.The blue line is the regression fit including 95 % confidence intervals (in grey).

Table 3 .
The ranking of the parameters defining GPP in the ANN by replacing each input variable with a random permutation of its values.(a) The parameters with their mean square error (MSE; mol m −2 day −1 ) for the model without O 3 .(b) The parameters with their MSE for the model with O 3 .The overall model MSE without any random permutation is also shown.

Figure 8 .
Figure 8.The measured versus the modelled gross primary productivity (GPP) for days used for model training and testing (a, b), for days on which an O 3 effect was assumed (c, d), and for the entire growing season (e, f).GPP model 1 was trained without days with the highest stomatal O 3 uptake, whereas GPP model 2 was trained to test for possible lag effects of O 3 on GPP.The black lines are the fitted linear regression lines and the grey lines mark the 95 % confidence bands.Also shown are the p values for the tests of the slope and intercept from the regression y = ax + b being different from 1 and 0, respectively.

Figure 9 .
Figure 9.The residuals of growing season gross primary productivity (GPP) as a function of (a, b) the total stomatal O 3 flux over the growing season (F st ), (c, d) AOT40, and (e, f) POD 1 .PLA is the projected leaf area.The negative residuals indicate model overestimation of GPP.GPP model 1 was trained without days with the highest stomatal O 3 uptake, whereas GPP model 2 was trained to test for possible lag effects of O 3 on GPP.The black lines are the fitted linear regression lines and the grey lines mark the 95 % confidence bands.Also shown are the p values for the tests of the slope and intercept from the regression y = ax + b being different from 0 (n = 14).

Table 1 .
The optimised parameter values of the multiplicative stomatal model.

Table 2 .
Performance statistics for the multiplicative stomatal model: mean bias (MB), relative mean error (RME), systematic and unsystematic root mean square error (RMSE s/u ), Willmott's index of agreement (d), model efficiency (ME), and the coefficient of determination (R 2 ).