1 No tropospheric ozone impact on the carbon uptake by 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 10 investigated the effect of high O3 events on gross primary production (GPP) for a Scots pine stand near Antwerp, Belgium over the period 1998-2013. Stomatal O3 fluxes were modelled using in situ O3 concentration 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 parameterised for days with low stomatal O3 uptake rates was used to simulate 15 GPP during periods of high stomatal O3 uptake. Eventual 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 PODy. Although the critical levels for both indices were exceeded in every single year, no significant negative effects were found of ozone on GPP and no correlations between GPP residuals and AOT40 and PODy were found. Overall, we conclude that no O3 effects were detected on the carbon 20 uptake by this Scots pine stand.


Introduction
Tropospheric ozone (O3) is a secondary air pollutant that has the potential to negatively affect vegetation, leading to reduced growth and carbon sequestration potential (ICP Vegetation, 2012;Middleton, 1956).Background concentrations of tropospheric O3 have increased with 36 % since pre-industrial times (IPCC, 2001) and are 25 projected to further increase considerably until about 2050 (IPCC, 2007).Depending on the scenarios, background O3 levels might either increase or decrease after 2050 (IPCC, 2007).
In recent years, many studies have been conducted to investigate the mechanisms underlying the O3 impacts on vegetation.Ozone reduces plant growth by altering photosynthetic rates, carbohydrate production, carbon sequestration, carbon allocation, and carbon translocation (Reich and Amundson, 1985;Andersen et al., 30 1997;Beedlow et al., 2004).Once O3 enters the leaves through the stomata, it can affect plant growth by direct cellular damage (Mauzerall and Wang, 2001), leading to visible leaf injury and reduced leaf longevity (Noble and Jensen, 1980).In response to O3, respiratory processes increase, which will also effect the tree's carbon balance (Darall, 1989).Skärby et al. (1987) proved that dark respiration of Scots pine shoots increased after long-term exposure to a low level of O3.Protective responses, such as compensation (e. g. repair of injured tissue), avoidance Biogeosciences Discuss., doi:10.5194/bg-2016-12, 2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.
(e. g. stomatal closure), and tolerance (e. g. alteration of metabolic pathways), all consume carbon and, hence, resistance to O3 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 O3, several indices have been created, e. g.AOT40 (ppb h), the cumulated O3 concentration in excess of a threshold of 40 ppb, and PODy, the accumulated O3 flux above a flux threshold y (nmol m -2 s -1 ).

40
Critical levels, quantitative estimates of exposure to O3 above which direct adverse effects may occur (CLRTAP, 2015), have been determined for these indices based on O3 dose-response relationships from fumigation experiments with enhanced O3 concentrations (Karlsson et al., 2004).The magnitude of the O3 impact on plants depends on the intensity of O3 exposure, environmental factors influencing both plant photosynthesis and the O3 flux to plant surfaces, and plant species-specific defensive mechanisms (Musselman and Massman, 1999).Because 45 of the variable plant responses to similar O3 concentrations, the question arises whether widely applicable tolerable limits of O3 concentration exist (Skärby et al., 1998).
While high stomatal O3 fluxes have been shown to affect crop yields and tree seedlings, it is not sure whether O3 uptake or O3 flux also negatively affects the carbon uptake by mature forest trees.Many studies determined the effect of O3 on seedlings and young trees (Buker et al., 2015), but little is known about the effect on mature trees.

50
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;Chappelka and Samuelson, 1998).In addition to the uncertainties related with the upscaling 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 O3 uptake on carbon uptake under 55 ambient O3 concentrations by trees has hardly been studied in situ.Some studies showed reductions in plant growth due to stomatal O3 uptake (Zapletal et al., 2011;Fares et al., 2013;Yue and Unger, 2013), while other studies did not show any effect (Zona et al., 2014;Samuelson, 1994).Whether or not an effect of stomatal O3 uptake was found was species-and site-specific, and there is a clear need for more studies investigating the effect of O3 on carbon uptake by mature trees in the field (Chappelka and Samuelson, 1998).

60
Here we investigate the effect of high O3 events on gross primary production (GPP) for a Scots pine stand in Flanders, Belgium.At current ambient O3 levels, critical levels for both AOT40 and POD1 are already being exceeded in this Scots pine stand (Neirynck et al., 2012).This indicates that even at current ambient O3 levels tree productivity might be affected.Ozone-induced GPP reduction is most likely to occur during or shortly after days with high stomatal O3 uptake.An effect of stomatal O3 uptake on GPP can be detected when a GPP model 65 parameterised for days with low stomatal O3 uptake rates is extrapolated to high stomatal O3 uptake eventsi.e., days where an effect on GPP is assumed -and the model overestimates GPP during these events.This study therefore tests the hypothesis that GPP of the studied pine forest is reduced during or shortly after high stomatal O3 uptake events.Discuss., doi:10.5194/bg-2016Discuss., doi:10.5194/bg- -12, 2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Study area
The study area consisted of a 2-ha Scots pine stand in a 150-ha coniferous/deciduous forest named 'De Inslag', situated in Brasschaat (+51° 18' 33'' N, +04° 31' 14'' E), northeast of the Antwerp agglomeration and eastnortheast 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).

75
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 trees were planted in 1929 (Neirynck et al., 2008).In 1995, tree density amounted to 542 trees ha -1 .In the 80 autumn of 1999, the forest was thinned, which resulted in 376 trees ha -1 in 2001.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 one-year-old needles (Op de Beeck et al., 2010).

85
In this study, continuous measurements over the period 1998 -2013 were used, excluding 1999 and 2003 due to poor data quality or coverage.Different meteorological variables were measured on a tower with a height of 41 m, set up in 1996 (Gielen et al., 2013).A sonic anemometer (Model Solent 1012R2, Gill Instruments, Lymington, UK) measures the wind velocity (WV; m s -1 ).Meteorological data include vertical profiles of air temperature (Tair; °C) and humidity (RH; (%) (HMP 230 dew point transmitter and PT100, Vaisala, Finland) in aspirated radiation 90 shields at 2, 24 and 40 m height.Wind speed measurements (LISA cup anemometer, Siggelkow GMBH, Germany) are conducted at 24, 32 and 40 m height.At the top of the tower, ingoing and outgoing short-wave and long-wave radiation are measured by a CNR1-radiometer (pyranometer/pyrgeometer, Kipp and Zonen, the Netherlands) and a CMP6-pyranometer (Kipp and Zonen, the Netherlands).A wind vane (potentiometer W200P, Campbell, UK) is mounted on a tower rail.Rainfall is registered by a tipping bucket rain gauge (NINA precipitation pulse transmitter,
Soil water content (SWC; m 3 m -3 ) was measured at 25 cm below the soil surface using a TDR (Time Domain Reflectometer) sensor at three days to biweekly resolution and subsequently interpolated to obtain daily estimates, taking into account water inputs via precipitation (Gielen et al., 2010).Soil water potential (SWP; MPa) was 100 derived from SWC measurements (m 3 m -3 ) with the model of van Genuchten (van Genuchten, 1980).All climatic variables were measured every 10 seconds and half hourly means were stored on a data logger (Campbell CR1000, UK) in an air-conditioned cabinet adjacent to the tower.Data gaps were filled with data from nearby weather stations.
Vertical profiles of O3 concentrations are being measured at two inlets above the canopy (at 24 and 40 m) using measured using an infrared gas analyser (IRGA) (Model LI-6262, LI-COR Inc., Lincoln, NE, USA).The vertical CO2 flux between the forest and the atmosphere, net ecosystem exchange (NEE), was measured with the eddy covariance technique following standard data quality procedures (Carrara et al., 2003;Gielen et al., 2013;Carrara et al., 2004).Gross primary production (µmol C m -2 s -1 ) is derived from NEE data, by subtracting the modelled 110 total (autotrophic and heterotrophic) ecosystem respiration from the measured NEE.The ecosystem respiration or total carbon loss is modelled with standardised algorithms as presented in Falge et al. (2001).
The LAI time series for the Brasschaat forest was reconstructed based on the historical data.The general approach was to use the fragmentary LAImax data that were measured by Gond et al (1999) in 1997, by Konôpka et al. (2005) et al., 2005).All LAImax measurements were interpolated lineairly to derive LAImax values for the missing years.
The thinning event in 1999 was accounted for by subtracting the removed leaf biomass by using the allometric relations from Yuste et al. (2005) and specific leaf are measurements from Op de Beeck et al. (2010).The seasonal 120 pattern of the Op de Beeck et al. ( 2010) measurements was used and kept constant over the entire time series.

Stomatal conductance measurements
During the summers of 2007 and 2013, stomatal conductance to H2O (gst, H2O) of Scots pine needles was measured at the site.These gst,H2O measurements were based on gas exchange measurements (photosynthesis and transpiration), which were carried out with a leaf cuvette connected to an IRGA (LI-6400, LI-COR, Lincoln, 125 Nebraska, USA).
Stomatal measurements were carried out in highly different environmental conditions during 2007 (cold and wet summer) and 2013 (warm and dry summer).Diurnal stomatal measurements and stomatal responses to PAR, Tair, and VPD, both on one-year-old needles and on current-year needles, were carried out during the summers of 2007 (Op de Beeck et al., 2010) and 2013.All measurements in 2007 were carried out on the needles of the two trees 130 closest to the tower, both in the lower and the upper canopy.In 2013 only the tree closest to the tower was accessible for measurements.The total number of sets of needles on which measurements were carried out, amounts to 10 in 2007 and 16 in 2013.The LI-6400 instrument calculated the gst,H2O assuming the whole area of the cuvette (2x3 cm) was covered with the needles.To obtain the gst to O3 (from here on referred to as 'gst') of the six or eight needles in the cuvette, two corrections needed to be made: one for the needle area, because needles did 135 not completely cover the area of the cuvette, and a second one for converting the data from gst, H2O to gst.
The width of the needles was measured at 40 x magnification under a binocular microscope (M3 Wild, Wild Heerbrugg, Gais, Switzerland) in combination with an ocular equipped with a reticule (Leitz, Wetzlar, Germany, periplan, GW 10xm).The width per needle was measured at three places: at the top, in the middle, and at the base.
An average of those three measurements was multiplied by the length of the needle inside the cuvette, 3 cm.After 140 measuring the needle area with a microscope, the gst,H2O data were corrected for the lower needle area.These data were multiplied by 0.61, which is the ratio of the molecular diffusivities of water vapour and O3 in the air, to convert from gst,H2O to gst.

Multiplicative stomatal model: description
Stomatal conductance was modelled using the multiplicative gst model, first described by Jarvis (1976).The model 145 has been developed to calculate species-specific gst according to phenology and environmental conditions (Emberson et al., 2000) and is described in detail in Appendix A.
We modified the model to make it more applicable for Scots pine.In this modified model (Eq. 1) PAR, Tair, VPD, and SWP influence the range between gmax and gmin instead of gmax and zero.This modification was needed, because in the Brasschaat pine forest stomata never completely close, hence gst is never zero (Op de Beeck et al., 2010).
Here gst is the stomatal conductance to O3 and gmax is the maximal stomatal conductance to O3.The functions fPHEN, fPAR, fT, fVPD, and fSWP represent the modification of gmax by, respectively, phenology, PAR, Tair, VPD, and SWP.
The function fmin is the ratio of gmin and gmax where gmin is the minimal stomatal conductance to O3 (see Appendix A for more detailed information).Impaired stomatal aperture mechanisms (stomatal sluggishness) due to ozone 155 exposure (Paoletti and Grulke, 2010) were not included in the model development.

Multiplicative stomatal model: parameterisation and validation
For the optimisation of the parameters of the different functions in the model, we assumed that the phenology function was 1.This was deemed a fair assumption, because gst,H2O was measured on mature needles in the summer (July/August 2007 and 2013), in the middle of the growing season.

160
The data set, including measured gst, PAR, Tair, VPD, and SWP, was split into two subsets by grouping odd and even rows for data sorted by PAR.One set was then used for parameterisation, the other for validation.The stomatal model was parameterised using the computer program Matlab (version 2013a).The optimisation of all parameters was done with the function 'lsqcurvefit' in Matlab.It finds the best parameter values, starting with an initial value, to best fit the function of the stomatal model to measured gst and can thus be used to fit a nonlinear 165 function with more than two independent variables.All parameters of fPAR, fTair, fVPD, and fSWP were optimised separately, with initial values that were estimated visually from plots of the functions to the dataset.

Multiplicative stomatal model: model evaluation
The

Canopy model 180
We applied a canopy model to scale up gst, measured at leaf level, to the canopy level.The canopy model consists of different horizontal leaf layers and includes a radiation transfer model (Goudriaan, 1977), a solar elevation model (Campbell and Norman, 1998) and the modified multiplicative stomatal model (Emberson et al., 2000).The model is described in detail in Appendix C.
where gtot is the total conductance to O3 (mol (m² ground area -1 ) s -1 ); gaero is the aerodynamic conductance and is set to 1; gbl is the boundary layer conductance to O3; gcan is the canopy conductance.
The boundary layer conductance to O3 was calculated with the following formula (Baldocchi et al., 1987): where K is the von Karman constant (0.43); u* (m s -1 ) is the friction velocity, which is derived from the measured 195 momentum fluxes; Sc is the Schmidt number (1.07 for O3); Pr is the Prandtl number (0.72 for O3); 44.64 mol m -3 is the molar density of air, and is applied for converting the unit of gbl from m s -1 to mol m -2 s -1 .
The canopy conductance consisted of a stomatal and a non-stomatal component.Since the stomatal component varies throughout the canopy, the canopy was divided into eight sublayers so that the leaves were evenly distributed between the horizontal layers.Dividing the canopy into sufficient sublayers was necessary in order to model fluxes 200 well.Eight sublayers were considered to be sufficient, as indicated in a sensitivity test with more and less sublayers.
For each leaf layer, the model calculates gst for sunlit and shaded needles, taking the solar elevation angle into account.Non-stomatal conductance was assumed to be constant over the canopy and was set at 0.16.This value was derived from long-term O3 flux measurements in Brasschaat (Neirynck et al., 2012).To obtain an O3-damage free GPP model, data from days for which an O3 effect was expected were removed from the dataset.These were the days with stomatal O3 uptake values in the upper 2%, 5 %, and 10% of stomatal O3 220 flux.As the results of a 2% and 10% cut-off were equal to a 5% cut-off, we report only results of a 5% cut-off.
With 2/3 of the data from days with low stomatal O3 uptake, the artificial neural network was trained.The other 1/3 was used to test the model.This O3-damage free GPP model was then run with all data.The absolute and relative differences in GPP accumulated over the growing season between EC-derived and modelled values were calculated, to investigate whether or not there was a reduction of GPP.

225
The relation between the residuals of total GPP and both AOT40, POD1 and POD2 was examined.Therefore, a linear fit between the residuals and the indices was made.A significant negative correlation would exist if the slope is significant different from 0 (p < 0.05) and intercept is not significant different from 0 (p > 0.05).These yearly GPP residuals were also plotted to the stomatal O3 flux to investigate their relation and a linear fit was made of which the significance was tested.If GPP was increasingly overestimated in the presence of higher stomatal O3 230 fluxes, this would indicate a deleterious O3 effect.
Ozone effects possibly appear and last during a period of several days after the O3 peaks, and as a result they will not be detected in the above analyses.Due to these possibly lag effects of O3, the above analyses were repeated, but now excluding the days with high stomatal O3 uptake along with the two subsequent days removed from the training and testing datasets.

Multiplicative stomatal model and simulated O3 fluxes
The best fitting parameter values for the multiplicative stomatal model are presented in Table 1 and different statistics to evaluate the model performance are presented in Table 2.For the parameterisation dataset, the measured data were fitted against modelled gst and plotted in Fig. 4. The slope of the linear fit was not significantly 250 different from 1 (p > 0.05) and the intercept was not significantly different from 0 (p > 0.05).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.05) and the intercept was not significantly different from 0 (p > 0.05).
In Fig. 5 the measured gst was plotted against the different model input variables: PAR, Tair, VPD, and SWP, and for each plot the boundary function was fitted.
The average daily O3 fluxes for the different years are presented in Fig. S2.The daily total Fst ranges from 1.42 to 2.00 nmol O3 m -2 day -1 .In 2011 the daily total Fst was the lowest, while in 2002 the highest stomatal flux was reached.The annual average ratio Fst/Ftot varied between 24-28 % (Fig. S2).We observed the lowest ratios in the beginning and at the end of the growing season.Above-average ratios were observed at the peak of the growing 260 season.

Ozone effects on GPP
Total GPP (mol C m -2 day -1 ) was calculated for days with low stomatal O3 uptake, high stomatal O3 uptake and for the entire growing season, using both the EC-derived GPP data and the modelled GPP data (Table 3).For days with low stomatal O3 uptake, the average daily total GPP was 0.48 mol C m -2 day -1 , and the models reproduced

265
GPP very well (Table 3).When we calculated total GPP for days with high stomatal O3 uptake, the EC-derived fluxes were much higher than for the days with low stomatal O3 uptake.This was probably due to the higher irradiation that typically occurs during high O3 events and stimulates GPP.The higher GPP, however, also suggests that negative O3 effects on GPP were highly unlikely.This is exacerbated by the fact that our models almost consistently underestimate GPP during high O3 events (Table 3), whereas we hypothesised the exact opposite, 270 namely that the models would overestimate GPP during these events because they were parameterised for low O3 days.We also observed no differences between both models, suggesting no lagged O3 effects on GPP (Table 3).
A weak, negative correlation between total GPP residuals and Fst exists for the GPP model trained without days with high stomatal O3 uptake (Fig. 7, A), while a small positive correlation is shown for the GPP model which tested for lag effects of O3 (Fig. 7, B).However, these differences were not statistically significant at p<0.05.For 275 both models, correlations between total GPP residuals and AOT40, and between total GPP residuals and both POD1 and POD2 existed.These correlations were also not statistically significant at p<0.05 (Fig. 7, C, D, E, F, G, and H).

Multiplicative stomatal model 280
All statistics shown in Table 2 clearly indicated that the fitted multiplicative stomatal model performed well.For both parameterisation and validation datasets, the model explained 72 % of the variance in gst.For both datasets, slope and intercept of the linear regression lines of measured versus modelled gst were not significantly different from 1 and 0, respectively (Fig. 4).Moreover, the model efficiency (ME in  (Willmott et al., 1985), which was the case for this model.Low mean bias (MB) and low mean relative error (MRE) further indicated very good performance.
The good performance of the model can also be observed in Fig. 5, in which the boundary lines represented the response of gst to the independent variables when other variables were not limiting (Martin et al., 1997).The 290 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).
Multiplicative stomatal models based on Jarvis (1976) have been parameterised earlier for generic Scots pine forests in Europe (Mills et al., 2011;Buker et al., 2015) and used to estimate critical levels for this species.

295
However, the empirical the dose-response relationship for Scots pine is based on only one two-year fumigation study on small seedlings and, therefore, high uncertainty exists in the modelled O3 impact on Scots pine growth.
The parameterisation of Mills et al. (2011) andBüker et al. (2015) differ from that of this study in a number of parameters.First, the needles of the Scots pine stand in Brasschaat had a higher night-time gst (gmin) and will therefore take up more O3 at night.Maximal gst, in contrast, is lower in Brasschaat than estimated for other Scots 300 pine forests, implying that during episodes of high O3 concentrations, the Brasschaat site is unlikely to take up very high amounts of O3.This may have contributed to the absence of a clear O3 response at our study site.Also the Scots pine stand in Brasschaat is less sensitive to drought stress than the generic model, due to a higher VPDmax and a wider SWP range.The wider SWP range is mainly due to a clearly lower SWPmax.These differences between the parameter values and, hence, in gst for generic Scots pine forests and for the Scots pine stand in Brasschaat will lead to different critical levels and under-or overestimation of possible O3 damage.Species-specific parameterisation is important, but site-specific parameterisation is clearly important as well.

Stomatal O3 fluxes at canopy scale
The stomatal O3 flux contributed on average for 26 % to the total O3 flux over the study period (Fig. S2).This fraction is similar to the 21 % stomatal O3 flux in a Danish Norway spruce stand (Mikkelsen et al., 2004) and the 310 30 % stomatal O3 flux in Quercus ilex in Italy (Vitale et al., 2005;Gerosa et al., 2005).Cieslik (2004) showed that in Southern Europe stomatal O3 flux of different vegetation types, such as pine forest and Mediterranean shrubs, is typically less than 50 % of the total O3 flux.A five-year study on a Mediterranean Pinus ponderosa stand showed a stomatal O3 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 O3 uptake (Neirynck et al., 2012).

315
The low relative stomatal O3 flux in the Scots pine stand in Brasschaat could be the result of the sparse canopy with low LAI.Although no relation between stomatal O3 flux and LAI was found in a previous site study on this site (Neirynck et al., 2012), interannual and seasonal variation in LAI is very small, rendering such a correlation analysis very difficult.

Ozone effects on GPP 320
Although our models reproduced GPP very well, we did not observe immediate or lagged effects of high stomatal O3 uptake on GPP (Table 3; Fig. 7, A, B).Some earlier studies have investigated the effect of O3 on forest carbon uptake.Cumulative stomatal uptake of 27 mmol m -2 over the growing season did not result in any visible damage or a reduction in NEE of a poplar plantation in Belgium (Zona et al., 2014).Zapletal (2011), on the other hand, reported that CO2 uptake of a Norway spruce forest in the Czech Republic increased with increasing stomatal O3 325 flux, followed by a sudden decrease in CO2 uptake, suggesting that an O3 flux threshold exists.Fares (2013) showed a negative correlation between GPP and O3 uptake at two Mediterranean ecosystems (a forest dominated by Pinus ponderosa in California, USA and an orchard site of Citrus sinensis cultivated in California, USA).A GPP reduction of 1-16% in response to O3 uptake under ambient O3 concentrations of 30-50 ppb was determined across vegetation types and environmental conditions in the United States by Yue and Unger (2013).The 330 magnitude of reduction depended on the sensitivity to O3 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 % (Mills et al., 2011).In this study on Scots pine in Brasschaat, this value was far exceeded in all years (Fig. 7), yet no negative effect on GPP was observed in years with higher AOT40 values.Particularly noteworthy is the extreme high AOT40-value of 2006, which was due to the high O3 335 concentrations during that year, which, nevertheless, did not result in GPP reductions (Table 3).
PODy is considered a more appropriate index for potential O3 damage because it considers O3 flux.The critical level of POD1 is species-specific; a critical level of 8 mmol m -2 with 2 % growth reduction is used for Norway spruce and a critical level of 4 mmol m -2 with 4 % growth reduction is used for birch and beech (Mills et al., 2011).spruce is often adopted as critical level for Scots pine.During this study, this critical level was exceeded every single year, and again no negative correlation between total GPP residuals and POD1 was observed.In comparison to the AOT40 level, 2006 was not the year with the highest POD1.This difference between AOT40 and POD1 during 2006 was due to stomatal closure; during high O3 concentration events, gst was rather low (Fig. 6).POD1 was highest in the year 2002, when O3 concentrations were relatively low, but gst was high.The low O3 345 concentrations explain the lower AOT40 for 2002.
Notwithstanding the absence of a statistically significant negative correlation between GPP residuals and both AOT40 and PODy, critical levels for both AOT40 and PODy were exceeded every single year.AOT40 is based on O3 concentration and these concentration-based indices have been shown to be weaker indicators for O3 damage than flux-based indices (Karlsson et al., 2007;Simpson et al., 2007).The critical level of PODy for Scots pine was 350 adopted from the critical level for Norway spruce (Mills et al., 2011).Possibly this critical level is too low for Scots pine.As shown by Reich (1987), pines are less sensitive to O3 compared to hardwoods and crops.This supports the idea of a too low critical level.
Overall, no significant O3 effects on GPP accumulated over the growing season were found.Although no significant O3 effects on GPP were found in this study, it is still possible that O3 negatively affected this Scots pine 355 stand in Brasschaat.Stomatal O3 uptake was linked 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 trees at our study site were able to fully detoxify the incorporated O3.As a result, no O3 effects on carbon uptake were detectable.However, this protection may have come at a respiratory cost, which may have reduced the NPP/GPP ratio of this forest.The NPP/GPP ratio of our study site was very low (Gielen et al., 2013).In 360 addition to the poor nutrient status (limitation by P and Mg, extremely high N deposition; (Neirynck et al., 2008)), O3 uptake may partly be responsible.This can, however, not be tested because pine forest NPP data were not available at annual timescale.

Summary
We parameterised a multiplicative stomatal model for a Scots pine stand in Brasschaat.This species-and site-365 specific parameterised model performed very well.With this model, stomatal O3 fluxes were calculated and used to test for O3 effects on GPP.Only very small reductions in growing season GPP were calculated.Although critical levels for AOT40 and PODy were exceeded in every single year, no significant negative correlations between total GPP residuals and stomatal O3 flux, AOT40, and PODy were found.In general, we can thus conclude that no O3 effects were detected on the carbon uptake by the Scots pine stand in Brasschaat.First, the canopy is divided into different horizontal layers, with each a sunlit and shaded fraction.For each layer fraction the incoming PAR is calculated with the radiation transfer submodel.With the multiplicative stomatal model gst is calculated for each layer fraction.460 For each layer fraction the total leaf conductance (mol m -2 leaf area s -1 ) is calculated by summing gst and gns, the non-stomatal conductance.This value is multiplied by LAI of the layer fraction and the values for both the sunlit and the shaded layer fraction are summed to obtain the total layer conductance (mol m -2 ground area s -1 ).All layer conductances can be summed to obtain the conductance of the canopy as a whole (gcan).465 The total conductance is a function of gbl and gcan based on an electrical conductance analogy model.
where gaero is the aerodynamic conductance; gbl is the boundary layer conductance; gcan is the canopy conductance.
Total O3 flux (nmol m -2 ground area s -1 ) is the O3 flux for the whole canopy and is then calculated by: 470 where [O3] is the O3 concentration (ppb).where gst is the stomatal conductance (mol m -2 ground area s -1 ); gns is the non-stomatal conductance (mol 475 m -2 ground area s -1 ).
Non-stomatal O3 flux (Fns) is the difference between total O3 flux and stomatal O3 flux:

Part 2 The solar elevation submodel
This submodel calculates the solar elevation angle, β (radians), at each time step (Campbell and Norman, 480 1998).First the sunlit LAI fraction of each horizontal leaf layer i is calculated with Beer's law: where kb is the direct radiation extinction coefficient; Ω is a factor accounting for inter-and intra-crown foliage clumping; LAIc(i) is the cumulative LAI above a horizontal leaf layer i. 500 A spherical needle angle distribution is assumed, hence   = 0.5/ sin  (de Pury, 1997).
The shaded LAI fraction of each horizontal leaf layer i is calculated as follows: where fsun(i) is the sunlit LAI fraction.
The intensity of direct radiation does not decline through the canopy, but the diffuse radiation does and 505 is calculated with Beer's law: where Id0 is the incoming diffuse radiation.
in 2003 and by Op de Beeck et al. (2010) in 2007.The latter was done with hemispherical pictures while 115 the first in 1997 and 2003 were done with the LAI-2050 instrument (LI-COR, Lincoln, Nebraska, USA).To assure consistency across the time series, measurements were corrected for clumping by using the factor 0.83 (Jonckheere Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.
parameterised multiplicative stomatal model was then tested against the validation dataset.Measured gst values were plotted against the modelled gst values.A linear function  =  +  was fitted, where 'a' should be not 170 significantly different from one (p > 0.05) and 'b' should be not significantly different from zero (p > 0.05) for both parameterisation and validation dataset.We evaluated the model performance with the following statistics: the coefficient of determination or R squared (R²) as a goodness-of-fit measure and error-based measures including Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.mean bias (MB), relative mean error (RME), Willmott's index of agreement (d), model efficiency (ME), root mean squared error (RMSE), and its systematic (RMSEs) and unsystematic component (RMSEu).In Appendix B these 175 error-based statistics are explained.The measured gst was plotted in function of the different input variables (PAR, Tair, VPD, and SWP) and the boundary function of each plot was fitted.This was done in order to test how well the obtained parameter values were estimated in function of the measured gst.
Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.The stomatal and non-stomatal O3 fluxes (nmol m -2 s -1 ) were calculated by multiplying the proportion of gst and 205 gns of the canopy per ground area with the O3 concentration.These obtained half-hourly fluxes were aggregated to daily fluxes.These daily fluxes were averaged in order to know the average daily O3 uptake by the canopy for the different years.The ratio Fst/Ftot was calculated and this gives an indication of the contribution of the stomatal O3 flux to the total O3 flux.2.6 Ozone effects210A feed-forward back propagation Artificial Neural Network (ANN) in Matlab (Matlab Matworks R2013a, The MathWorks Inc., Natick, Massachusetts, USA) was used to simulate GPP of the Scots pine forest.The ANN, which contained 10 nodes organised in 1 layer, was trained with 70% data random selected data measured data and validated based on the remaining 30% of data set (R²=0.72).The daily GPP data of the growing seasons between1998-2013, except 1999 and 2003, were used as dependent target variable, whereas year, day of year, 215 Tmin, Tmax, Tmean, average VPD, SWC, Rg, average Tsoil, and average WV of these periods were used as independent input variables.Daily totals of the variables were used, with the exception of VPD, Tsoil, and WV for which daily averaged values were used.
figure gives a good overview of how meteorology and GPP typically changed over time in this forest; interannual anomalies from the average patterns can be found in Fig.S1.Distinct daily and seasonal patterns can be observed 240 Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.

A
critical level for Scots pine has not yet been determined and therefore the value of 8 mmol m -2 for Norway Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.
Stomatal O3 flux is the fraction of the total O3 flux taken up by the stomata and is calculated by: Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.
The total received radiation by a sunlit fraction (Isun(i)) is the sum of direct and diffuse radiation.Shaded leaves only receive diffuse radiation: 510  ℎ() =  () the averaged leaf angle for a uniform needle angle distribution; Ib0 is the incoming direct radiation at top of the canopy; Id(i) is the diffuse radiation for a horizontal leaf layer i.Total received irradiance is now converted to total received PAR and split into PAR per horizontal leaf 515 layer.Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Fig. 1 .
Fig. 1.Fingerprint of the meteorological and eddy-flux derived gross primary productivity (GPP) measurements averaged over the period 1998-2013.Day of year is plotted on the y-axis and hour of day on the x-axis.Air temperature (Tair), incoming global radiation (Rg), vapour pressure deficit (VPD), and GPP are plotted.

Fig. 2 .
Fig. 2.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.Error bars represent the 95% confidence intervals.

Fig. 3 .
Fig. 3. Seasonal changes of LAI over the 15 year study period at the Brasschaat Scots pine site.

Fig. 5 .
Fig. 5. Measured stomatal conductance (gst) in function of the different variables used in the multiplicative model: photosynthetically active radiation (PAR), air temperature (Tair), vapour pressure deficit (VPD), and soil water potential (SWP).The red line represents the boundary line for which the functions are given in Appendix A (A3-A6). (n = 205) 715

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 285 provides low root-mean-square error (RMSE), while the systematic component (RMSEs) should approach zero and the unsystematic component (RMSEu) should approach RMSE

Table 2 .
Statistics of the model evaluation.The statistics used to evaluate the model performance are mean bias Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Table 3 .
GPP values (mol C m -2 day -1 ) for the days with low stomatal O3 uptake, high stomatal O3 uptake, and the growing season total (GS) of different years as well as the multi-annual mean.One GPP model (GPPmod1) was trained without days with high stomatal O3 uptake, whereas a second GPP model (GPPmod2) also excluded days 685 following high O3 events, and was thus trained to test for lag effects of O3.The relative difference between modelled and EC-derived GPP estimates is presented between brackets.Positive values indicate an overestimation by the model and therefore a potential O3 effects on GPP.Biogeosciences Discuss., doi:10.5194/bg-2016-12,2016 Manuscript under review for journal Biogeosciences Published: 23 February 2016 c Author(s) 2016.CC-BY 3.0 License.