Development and evaluation of an ozone deposition scheme for coupling to a terrestrial biosphere model

. Ozone ( O 3 ) is a toxic air pollutant that can damage plant leaves and substantially affect the plant’s gross primary production (GPP) and health. Realistic estimates of the effects of tropospheric anthropogenic O 3 on GPP are thus potentially important to assess the strength of the terrestrial biosphere as a carbon sink. To better understand the impact of ozone damage on the terrestrial carbon cycle, we developed a module to estimate O 3 uptake and damage of plants for a state-of-the-art global 5 terrestrial biosphere model called OCN. Our approach accounts for ozone damage by calculating (a) O 3 transport from 45 m height to leaf level, (b) O 3 ﬂux into the leaf, and (c) ozone damage of photosynthesis as a function of the accumulated O 3 uptake over the life-time of a leaf. A comparison of modelled canopy conductance, GPP, and latent heat to FLUXNET data across European forest and grassland sites shows a general good performance of OCN including ozone damage. This comparison provides a good baseline on 10 top of which ozone damage can be evaluated. In comparison to literature values, we demonstrate that the new model version produces realistic O 3 surface resistances, O 3 deposition velocities, and stomatal to total O 3 ﬂux ratios. A sensitivity study reveals that key metrics of the air-to-leaf O 3 transport and O 3 deposition, in particular the stomatal O 3 uptake, are reason-ably robust against uncertainty in the underlying parameterisation of the deposition scheme. Nevertheless, correctly estimating canopy conductance plays a pivotal role in the estimate of cumulative O 3 uptake. We further ﬁnd that accounting for stomatal 15 and non-stomatal uptake processes substantially affects simulated plant O 3 uptake and accumulation, because aerodynamic resistance and non-stomatal O 3 destruction reduce the predicted leaf-level O 3 concentrations. Ozone impacts on GPP and transpiration in a Europe-wide simulation indicate that tropospheric O 3 impacts the regional carbon and water cycling less than expected from previous studies. This study presents a ﬁrst step towards the integration of atmospheric chemistry and ecosystem dynamics modelling, which would allow to assess the wider feedbacks between vegetation ozone uptake and tropospheric ozone burden.


Introduction
Tropospheric ozone (O 3 ) is a highly reactive and toxic gas. It enters the plants mainly through the stomata of the leaf, where it forms reactive oxygen species (ROS) which have the potential to damage the leaf. While leaves possess physiological pathways to produce compounds like ascorbate and polyamines, which help to neutralise the oxidising power of ROS (Kronfuß et al., 1998;Kangasjärvi et al., 1994;Tausz et al., 2007), ozone injury may occur when the leaf's anti-oxidant system becomes overwhelmed (Wieser and Matyssek, 2007).
In Western Europe, tropospheric O 3 levels have increased approximately by a factor 2 to 5 from pre-industrial values to the 1990s (Cooper et al., 2014;Marenco et al., 1994;Staehelin et al., 1994) (although the low values at the start of this period are very uncertain) and approximately doubled between 1950 and 1990s in the northern hemisphere (Parrish et al., 2012;Cooper et al., 2014). The major causes for this increased O 3 formation are the increased emission of O 3 precursor trace gases such 10 as NO x and CO, primarily from combustion sources, non-methane volatile organic compounds from anthropogenic sources (combustion, solvents), and methane emissions from agriculture and industry (Fusco and Logan, 2003;Vingarzan, 2004). For instance, in Western Europe, NO x emissions have risen by a factor of 4.5 between 1955 and 1985 (Staehelin et al., 1994).
In addition, downward transport of O 3 from the stratosphere to the troposphere (Vingarzan, 2004;Young et al., 2013) and intercontinental transport (Vingarzan, 2004;Jenkin, 2008;Fiore et al., 2009) can increase local and regional O 3 concentrations. 15 A commonly observed consequence of elevated levels of O 3 exposure is a decline in net photosynthesis (Morgan et al., 2003;Wittig et al., 2007), which may result from the damage of the photosynthetic apparatus or increased respiration due to the production of defence compounds and investments in injury repair (Wieser and Matyssek, 2007;Ainsworth et al., 2012).
The reduction in net photosynthesis results in reduced growth and hence a reduced leaf area, and plant biomass (Morgan et al., 2003;Lombardozzi et al., 2013;Wittig et al., 2009). The tight coupling between photosynthesis and stomatal conductance 20 further affects canopy conductance, and thereby transpiration rates (Morgan et al., 2003;Wittig et al., 2009;Lombardozzi et al., 2013), likely affecting the ecosystem water balance.
Due to its phytotoxic effect, elevated O 3 levels as a consequence of anthropogenic air pollution may affect the land carbon cycle, and potentially reduce the net land carbon uptake capacity (Sitch et al., 2007;Arneth et al., 2010;Simpson et al., 2014a), which currently corresponds to about a quarter of the anthropogenic fossil fuel emissions as a result of a sustained imbalance 25 between photosynthetic carbon uptake and carbon loss through respiration and disturbance processes (Le . However, the extent to which O 3 affects plant health regionally and thereby alters terrestrial biogeochemistry and the terrestrial water balance is still subject of large uncertainty (Simpson et al., 2014a).
A number of O 3 exposure indices have been proposed to assess the potential detrimental effect of tropospheric O 3 on the plants (LRTAP-Convention, 2010;Mills et al., 2011b). In Europe, the standard method of these indices is the concentration- 30 based AOTX [ppb h] (accumulated O 3 concentration over a threshold of X ppb), which relates the free-air O 3 concentration to observed plant damage. Models assessing ozone damage to gross or net primary production based on AOTX have been used for many years and indicate that substantial reduction in plant growth and carbon sequestration occurs globally and may reach reductions of more than 40 % at O 3 hot spots (Felzer et al., 2004(Felzer et al., , 2005Ren et al., 2011;Anav et al., 2011). A significant caveat of concentration-based assessments of ozone toxicity effects is that species differ vastly in their canopy conductance as well as regional provenances of species. Stomatal control of the leaf gas exchange regulates photosynthesis, and varies inter alia with plant specific photosynthetic capacity and intrinsic water-use efficiency of photosynthesis, phenology, as well as environmental factors such as incident light, atmospheric vapour pressure deficit (VPD), air temperature. The consequent differences in stomatal conductance implies that the actual O 3 dose, and thus the level of ozone-related damage, 5 differs between species exposed to similar atmospheric O 3 concentrations (Wieser and Havranek, 1995). The O 3 dose, that is the integral of the instantaneous O 3 stomatal flux over a given period of time, has been observed to strongly correlate with the amount of injury of a plant, suggesting that plants with higher stomatal conductance are subject to higher doses and hence more susceptible to injury (Reich, 1987;Wittig et al., 2009).
Accounting for the O 3 dose rather than the O 3 exposure in assessments of ozone damage results in diverging regional 10 patterns of ozone damage, as regions with the highest exposure (O 3 concentrations) do not always coincide with regions of high uptake (Emberson et al., 2000;Mills et al., 2011a;Simpson et al., 2007). Regions with low AOT40 (AOTX above a threshold of 40 ppb) values might show moderate to high values of O 3 uptake because the flux approach accounts for climatic conditions that enable high stomatal conductances and hence high values of O 3 uptake (Mills et al., 2011a). Observed ozone damage in the field seems to be better correlated to flux-based risk assessment compared to concentration based methods (Mills 15 et al., 2011a). Following this the LRTAP Convention recommends flux based methods as the preferred tool for risk assessment (LRTAP-Convention, 2010).
When calculating the O 3 uptake into the plants, it is important to consider that stomatal uptake is not the only surface sink of O 3 . O 3 destruction also occurs at non-stomatal surfaces such as the leafs cuticle and soil surface. The stomatal flux represents approximately half of the total O 3 flux to the surface (Gerosa et al., 2004;Fowler et al., 2009;Cieslik, 2004;Simpson et al., 20 2003). Accounting for this non-stomatal O 3 deposition reduces the amount of O 3 uptake into the plants by reducing the surface O 3 concentration (Tuovinen et al., 2009) and thus has the potential to affect flux-based ozone damage estimates.
A further challenge in estimating plant damage related to O 3 uptake is that plants differ in their ability to remove any ROS from the leaf before damage of leaf cellular organs is incurred (Luwe and Heber, 1995). Conceptionally, one can describe the capacity as a plant-specific O 3 dose, with which the anti-oxidant system of the leaves can cope such that no damage is observed 25 (Musselman et al., 2006). The production of defence compounds increases respiration costs and following this reduces net primary production what may result in reduced growth and biomass (Ainsworth et al., 2012). Ozone damage is only incurred, once the O 3 flux into the leaf exceeds this dose. A commonly used index to assess flux-based damage to plants is the PODy [Phytotoxic Ozone Dose, nmol m −2 s −1 ], which gives the accumulated O 3 flux above a threshold of Y nmol m −2 s −1 for all daylight hours and a given time period. Common threshold values for PODy range from 1-6 nmol m −2 s −1 (Pleijel et al., 2007;30 LRTAP-Convention, 2010;Mills et al., 2011b), depending on the specific species sensitivity to O 3 .
Only a few terrestrial biosphere models have adopted the flux approach to relate O 3 exposure to plant damage and thus estimate O 3 induced reductions in terrestrial carbon sequestration in a process-based manner. Sitch et al. (2007) developed a version of the JULES model in which stomatal O 3 uptake directly affects net primary production (NPP), thereby ignoring the effect of reduced photosynthesis under elevated levels of O 3 on water fluxes. Lombardozzi et al. (2015) proposed a revised 35 version of the CLM model, in which O 3 imposes fixed reductions to net photosynthesis for two out of three modelled plant types. Atmospheric O 3 concentrations and the amount of cumulated O 3 uptake directly affect net photosynthesis only for one plant type.
In this paper, we present a new, globally applicable model to calculate O 3 uptake and damage in a process-oriented manner, coupled to the terrestrial energy, water, carbon and nitrogen budget of the OCN terrestrial biosphere model (Zaehle and Friend,5 2010).
In this model, the canopy O 3 abundance is calculated using aerodynamic resistance and surface resistances to soil surface, vegetation surfaces and stomatal cavities to take account of non-stomatal O 3 destruction. Canopy O 3 abundance is used to simulate stomatal O 3 uptake given instantaneous values of net photosynthesis and stomatal conductance. O 3 uptake and its effect on net photosynthesis is then calculated based on an extensive meta-analysis across 28 tree species by Wittig et al. (2007) 10 considering the ability of plants to detoxify a proportion of the O 3 dose (Sitch et al., 2007).
We first give a detailed overview of the ozone scheme (Section 2.1); evaluate modelled gross primary production (GPP), canopy conductance, latent heat fluxes and LAI against data from the FLUXNET database (Baldocchi et al., 2001) to test the ability of the model to simulate observed values of key components affecting calculate O 3 uptake (Section 3.1); evaluate the simulated O 3 metrics against reported values in the literature (Section 3.2); provide a sensitivity analysis of critical variables 15 and parameters of the deposition model to evaluate the reliability of simulated values of O 3 uptake (Section 3.3); give an estimate of the effect of the present-day O 3 burden on European GPP and transpiration(Section 3.4); and estimate the impact of using the O 3 deposition scheme on O 3 uptake and cumulated uptake (Section 3.5).

Methods
We developed an ozone deposition and leaf-uptake module for the terrestrial biosphere model OCN (Zaehle and Friend, 2010). 20 OCN is a further development of the land-surface-scheme ORCHIDEE (O) (Krinner et al., 2005), and simulates the terrestrial coupled carbon (C), nitrogen (N) and water cycles for twelve plant functional types driven by climate data, atmospheric composition (N deposition, as well as atmospheric CO 2 and O 3 burden), and land use information (land cover and fertiliser application).
In OCN net photosynthesis is calculated for shaded and sun-lit leaves in a multi-layer canopy with up to 20 layers (each 25 with a thickness of up to 0.5 leaf area index) following a modified Farquhar-scheme and considering the light profiles of diffuse and direct radiation (Zaehle and Friend, 2010). Photosynthetic capacity depends on leaf nitrogen concentration and leaf area, which are both affected by ecosystem available N. Increases in leaf nitrogen content enable higher net photosynthesis and higher stomatal conductance per unit leaf area. This in turn affects transpiration as well as O 3 uptake and ozone damage estimates. Leaf N is highest in the top canopy and monotonically decreases with increasing canopy depth. Following this, 30 stomatal conductance and O 3 uptake is generally highest in the upper canopy and lowest in the bottom of the canopy.
The O 3 and N-deposition data used for this study are provided by the EMEP MSC-W (European Monitoring and Evaluation Programme Meteorological Synthesising Centre -West) chemical transport model (CTM) (Simpson et al., 2012). The O 3 flux and deposition modules used in the EMEP model are advanced compared to most CTMs, and have been documented in a number of papers (Emberson et al., 2001;Tuovinen et al., 2004Tuovinen et al., , 2009Simpson et al., 2007Simpson et al., , 2012Klingberg et al., 2008).
The ozone deposition scheme for OCN is adapted from the model used by EMEP MSC-W (Simpson et al., 2012) to fit the land surface characteristics and process descriptions of the ORCHIDEE model. The leaf-level ozone concentrations computed by EMEP can not directly be used by OCN, since EMEP and OCN differ in a number of properties, as for instance in the 5 number of simulated plant functional types, and importantly their ecophysiological process representation. Both models differ in the simulation of various ecosystem processes (e.g. phenology, canopy processes, biogeochemical cycles, and vegetation dynamics, which are more explicitly represented in OCN), which in sum impact stomatal and non-stomatal ozone deposition and through this the leaf-level ozone concentration. A possible further development of the new OCN is the coupling to a CTM, to allow for a consistent simulation of tropospheric O 3 burden and vegetation O 3 uptake.

Ozone module
The ozone deposition scheme calculates O 3 deposition to the leaf surface from the free atmosphere, represented by the O 3 concentration at the lowest level of the atmospheric chemistry transport model (CTM), taken to be at 45 m above the surface.
The total O 3 dry deposition flux (F g ) to the ground surface is calculated as (1) 15 where χ O3 atm is the O 3 concentration at 45 m and V g is the deposition velocity at that height. In OCN V g is taken to be dependent on the aerodynamic resistance (R a ), canopy-scale quasi-laminar layer resistance (R b ) and the compound surface resistance (R c ) to O 3 deposition.
R b is calculated from the friction velocity (u * ) as 5 The R a between 45 m height and the canopy is not computed by OCN and is inferred from the logarithmic wind profile (for more details see Appendix A). R c is calculated as the sum of the parallel resistances to stomatal/canopy (1/G O3 c ) and non-stomatal O 3 uptake (1/G ns ) (Simpson et al., 2012, eq. 55) The stomatal conductance to O 3 G O3 st (m s −1 ) is computed by OCN (Zaehle and Friend, 2010) as: where G O3 st is calculated as a function of net photosynthesis at saturating C i (A n,sat ) where g 1 is the intrinsic slope between A n and G st . It further depends on a number of scalars to account for the effect of soil moisture (f (Θ)), water transport limitation with canopy height (f (height)), and atmospheric drought (f ( q air )), as well as an empirical non-linear sensitivity to the internal leaf CO 2 concentration (f (C i )), all as described in Friend and Kiang (2005). The factor 1.51 accounts for the different diffusivity of O 3 from water vapour (Massman, 1998). The canopy conductance to O 3 G O3 c is calculated by summing the G O3 st of all canopy layers. To yield reasonable conductance values in OCN compared to FLUXNET data (see Sect. 3.1), the original intrinsic slope between A n and G c called α in Friend and Kiang (2005) is adapted such that g 1 = 0.7α.
The non-stomatal conductance G ns follows the EMEP approach (Simpson et al., 2012, eq. 60) and represents the O 3 fluxes between canopy air space and surfaces other than the stomatal cavities. The model accounts for O 3 destruction on the leaf 15 surface (r ext ), within-canopy resistance to O 3 transport (R inc ), and ground surface resistance (R gs ) where the surface area index SAI is equal to the leaf area index LAI for herbaceous PFTs (grasses and crops), and SAI = LAI + 1 for tree PFTs according to Simpson et al. (2012), to account for O 3 destruction on branches and stem. Unlike EMEP, we do not apply a day of the growing season constraint for crop exposure to O 3 , which in OCN is accounted for by the 20 simulated phenology and seasonality of photosynthesis. The external leaf-resistance (r ext ) per unit surface area is calculated as where the base external leaf-resistance (r ext,b ) of 2500 m s −1 is scaled by a low temperature correction factor F T and with 1 ≤ F T ≤ 2 and T s the 2 m air-temperature ( • C Simpson et al., 2012, eq. 60). For temperatures below -1 • C non-stomatal resistances are increased up to two times (Simpson et al., 2012;Zhang et al., 2003). The within-canopy resistance (R inc ) is calculated as where b is an empirical constant (set to 14 s −1 ) and h is the canopy height in m. The ground-surface resistance R gs is calculated with the constraint of 0 ≤ f snow ≤ 0.5, to prevent negative values in the first fraction of eq. 10. s d,max is taken to be 10 kgm −2 15 (Ducoudré et al., 1993).
Given these resistances, the canopy O 3 concentration (χ O3 c , nmol m −3 ) is then calculated based on a constant flux assumption χ O3 c and the stomatal conductance to O 3 (G O3 st in m s −1 ) are used to calculate the O 3 flux into the leaf cavities (F st , nmol m −2 s −1 ): According to Laisk et al. (1989) the leaf internal O 3 concentration (χ O3 i ) is assumed to be zero. The OCN implementation of deposition and flux described above is a simplification of the deposition system used by EMEP 5 in order to fit the process representation of ORCHIDEE, from which OCN has inherited its biophysical modules. The external leaf resistance is not included in the calculation of F st (Tuovinen et al., , 2009) what results in an overestimation of stomatal O 3 uptake. Further, OCN's calculation of R a is based upon neutral stability conditions (see Appendix), whereas the EMEP model makes use of rather detailed stability correction factors. However, a series of calculations with the full EMEP model have shown that the uncertainties associated with these simplifications are small, typically 0.5-5 mmol m −2 . As base-10 case values of POD0 are typically ca. 30-50 mmol m −2 in EU regions, these approximations do not seem to be a major cause of error, at least in regions with substantial ozone (and carbon) uptake. The full coupling of OCN to a CTM would be desirable to eliminate this bias and allow for a consistent calculation of tropospheric and surface near O 3 burdens.

Relating stomatal uptake to leaf damage
An accumulation of F st over time gives the accumulated uptake of O 3 for a particular canopy layer (CU O l , mmol m −2 ), or 15 for l = 1 (top canopy layer) the phytotoxic O 3 dose, (P OD , mmol m −2 ) where c = 10 −6 converts from nmol to mmol and the integration time step is 1800 seconds.
The CU O l is used to approximate the damage to net photosynthesis (A n ) by using the damage relationship of Wittig et al. (2007): where the factor 100 scales the percentage values of damage to fractions. Net photosynthesis accounting for ozone damage (A O3 n ) is then calculated by subtracting the damage fraction from the undamaged value of A n : Since G st and A n are tightly coupled (see eq. 5), a damage of A n results in a simultaneous reduction in G st . The canopy- Note that CUO and POD can be directly compared to estimates according to the LRTAP-Convention (2010) notation, when analysing only the top canopy layer (Mills et al., 2011b).

Sensitivity analysis
A sensitivity analysis is conducted to estimate the sensitivity of the modelled plant O 3 uptake to the parameterisation of 15 the model, to establish the robustness of the model, and to identify the most influential parameters. Three parameters (groundsurface resistance (R gs ), external leaf-resistance (r ext ), empirical constant (b), see eq. 10, 6, 9, respectively) and three modelled quantities (canopy conductance (G c ), aerodynamic resistance (R a ), and canopy-scale quasi-laminar layer resistance (R b ), see eq. 5, 2), with considerable uncertainty due to the underlying parameters used to calculate these quantities, are perturbed within ±20% of their central estimate. R c (eq. 4), V g (eq. 2) and the O 3 flux ratio (F R ) calculated as the ratio of F stC and the total O 3 flux to the surface (F g , eq. 1).
The summer months June, July, and August (JJA) are selected from the simulation output and used for further analysis. For , the sensitivity to changes in all six perturbed parameters/variables is estimated by calculating partial correlation coefficients (PCC) and partial ranked correlation coefficients (PRCC) (Helton and Davis, 2002).
PCC's record the linear relationship between two variables where the linear effects of all other variables in the analysis are removed (Helton and Davis, 2002). In case of nonlinear relationships, RPCC can be used, which implies a rank transformation to linearise any monotonic relationship, such that the regression and correlation procedures as in the PCC can follow (Helton and Davis, 2002). We estimate the magnitude of the parameter effect by creating mean summer values of the four prognostic variables for each sensitivity run, and regressing these values against the corresponding parameter/variable scaling values of the respective model run. NEE measurements are used to estimate gross primary production (GPP) by the flux-partitioning method according to (Reichstein et al., 2005). Canopy conductance (G c ) is derived by inverting the Penmen-Monteith equation given the observed LE and atmospheric conditions as described in Knauer et al. (2015).
The half-hourly FLUXNET and model fluxes are filtered prior to deriving average growing-season fluxes (bud break to litter 25 fall) to reduce the effect of model biases on the model-data comparison. Night-time and morning/evening hours are excluded by removing data with lower than 20 % of the daily maximum short-wave downward radiation. To avoid any biases associated with the soil moisture or atmospheric drought response of OCN, we further exclude data points with a modelled soil moisture constraint factor (range between 0-1) below 0.8 and an atmospheric vapour pressure deficit larger than 0.5 kPa.
Daily mean values are calculated from the remaining time steps only where both modelled and observed values are present. 30 The derived daily values are furthermore constrained to the main growing season by excluding days where the daily GPP is less than 20 % of the yearly maximum daily GPP.
To derive representative diurnal cycles, data for the month July are filtered for daylight hours (taken as incoming short-wave radiation ≥ 100 W m −2 ), and excluding periods of soil or atmospheric drought stress as above. This is done for modelled F stC , R c , V g , F R and for both modelled and FLUXNET observed GP P and G c .

Modelling protocol and data for regional simulations
For the regional simulations, OCN is run at a spatial resolution of 0.5 • x 0.5 • on a spatial domain focused on Europe. Daily 5 meteorological forcing (temperature, precipitation, short-wave and long-wave downward radiation, atmospheric specific humidity and wind speed) for the years 1961 to 2010 is obtained from RCA3 regional climate model (Samuelsson et al., 2011;Kjellstrom et al., 2011), nested to the ECHAM5 model (Roeckner et al., 2006), and has been bias corrected for temperatures and precipitation using the CRU climatology (New et al., 1999). Reduced and oxidised nitrogen deposition in wet and dry forms and O 3 concentrations at 45 m height for the same years are obtained from the EMEP model, which is also run with 3 Results

Evaluation against daily eddy-covariance data
Figure 1 a shows that, for most sites, modelled and observation-based GPP agree well (see Tab. 2 for R 2 and RMSE values).
The standard deviation is larger for the observation-based estimates because of the high level of noise in the eddy-covariance data. For sites dominated by needle-leaved trees, the modelled and observation-based GPP values are very close, with only When comparing modelled and observed latent heat fluxes (LE), the model fits the observations best at the needle-leaved forest sites (Fig. 1 c). However, LE is overestimated at nine out of ten broadleaved forest sites, but remains within the range of the large observational standard deviation. At sites dominated by C3 grasses the modelled LE differs considerably from observed value, at two sites overestimating and two underestimating the fluxes, again within the observational-standard deviation.
In agreement with the comparison of GPP and LE, the comparison of modelled to observation-based canopy conductance 5 (G c ) shows the best agreement for sites dominated by needle-leaved trees (Fig. 1 b). At sites dominated by broadleaved trees, the modelled G c varies more widely from the FLUXNET G c . The modelled G c at sites dominated by C3 grasses is in very good agreement to FLUXNET G c with slightly overestimating G c at 2 out of 3 sites except for the DE-Meh site, where means differ outside the standard deviation (see Fig. 10b).
The comparison of the average modelled summertime LAI and point measurements at the FLUXNET illustrates that the 10 variability in the measured LAI is much greater than that of OCN ( Fig. 1 d) For further evaluation of the modelled O 3 uptake, we analysed the diurnal cycles of four key O 3 variables (O 3 uptake (F stC ), 20 O 3 surface resistance (R c ), O 3 deposition velocity (V g ), and flux ratio (F R )) as well as GP P and G c at three sites, one of the three categories broadleaved, needle-leaved and C3 grass sites respectively. The selection criteria are that modelled and FLUXNET GPP and LAI agree well and a minimum of five observation years is available to reduce possible biases from the inability of the model to simulate short-term variations from the mean. The selected sites are a temperate broadleaved summer green forest (IT-Ro1), a boreal needle-leaved evergreen forest (FI-Hyy), as well as a temperate C3 grass land (CH-25 Oe1). We evaluate modelled GP P and G c against observations from the FLUXNET sites. The modelled mean diurnal cycles are compared to reported values in the literature since we did not have access to site-specific observations.
Modelled and observed mean diurnal cycles of GP P and G c are in general agreement at the three selected FLUXNET sites (see Fig. 2 a,g,m and b,h,n) with particularly good agreement for the mean diurnal cycle of GP P at the needle-leaved site FI-

13
The mean diurnal cycles of G c derived from the FLUXNET data are again best matched at the site FI-Hyy, whereas the model generally overestimates the diurnal cycle of G c slightly at the site IT-Ro1, and overestimates peak G c at the CH-Oe1 site.
The fact that OCN does not always simulate the observed midday depression of G c , suggests that the response of stomata to atmospheric and soil drought in OCN requires further evaluation and improvement. Similar to the daily mean values (see Fig. 1 a,b) the mean hourly values show the best match of GP P and G c for the needle-leaved tree site and stronger deviations for the 5 sites covered by broadleaved trees and C3 grasses.
The stomatal O 3 uptake F stC (Fig. 2 c,i,o) is close to zero during night time when the stomata are assumed to be closed, because gross photosynthesis is zero. At FI-Hyy and CH-Oe1, peak uptake occurred at noon at values between 8-9 nmol m −2 s −1 , when photosynthesis (Fig. 2 g,m) and stomatal conductance (Fig. 2 h,n) are highest. At the Italian site IT-Ro1, maximum uptake occurs in the afternoon hours around 15 h, with much larger standard deviation compared to the other two sites (Fig. 2 c)). 10 The magnitude of stomatal O 3 uptake corresponds well to some values reported e.g. for crops (Gerosa et al., 2003(Gerosa et al., , 2004 conditions  and 1-6 nmol m −2 s −1 for diverse southern European vegetation types (Cieslik, 2004). Much higher values are reported for Picea abies (50-90 nmol m −2 s −1 ), Pinus cembra (10-50 nmol m −2 s −1 ) and Larix decidua (10)(11)(12)(13)(14)(15) 40 nmol m −2 s −1 ) at a site near Innsbruck Austria (Wieser et al., 2003), where canopy O 3 uptake was estimated by sapflow measurements in contrast to the studies mentioned before where the eddy covariance technique was applied. The much higher F stC values in that study result from much higher canopy conductances to O 3 (G O3 c ), which are up to 12 times higher than the modelled G O3 c values in our study (see Fig. 2, G O3 c = Gc 1.51 ). The ratio between the stomatal O 3 uptake and the total surface uptake (F R ) is close to zero during night time hours and 20 increases steeply in the morning hours (Fig. 2 d,j,p). The 24 h average is approximately 0.3 for IT-Ro1 and 0.4 for FI-Hyy and CH-Oe1 (Fig. 2 d,j,p). Peak hourly mean values are close to 0.6 at IT-Ro1, around 0.7 at FI-Hyy and close to 0.8 at CH-Oe1.
These values are comparable to the ratios reported for crops (0.5-0.6 Gerosa et al., 2004;Fowler et al., 2009), Norway spruce (Mikkelsen et al., 2004, 0.3-0.33) and diverse southern European vegetation types (Cieslik, 2004, 0.12 -0.69 (Mikkelsen et al., 2004;Tuovinen et al., 2001). Mean daytime deposition velocities of 0.006 m s −1 (range 0.003-0.008 m s −1 ) are reported at a finish mountain birch site (Tuovinen et al., 2001). Simulated monthly mean values of V g differ substantially between the sites (see Fig. 11). When comparing the monthly means over all sites (Fig. 11 dashed line) of a functional group (broadleaved, needle-leaved, C3 grasses) to the ensemble mean of 15 CTM's (Hardacre et al., 2015) the values simulated here are higher for needle-leaved tree sites. For broadleaved tree sites and and grassland sites higher value but still within the observed ensemble range are found for the summer months.
The modelled hourly mean O 3 surface resistance R c is highest with approximately 400 sm −1 during night time and decreases during daytime to values of 100-180 sm −1 , where the lowest surface resistance of approximately 100 sm −1 is modelled at the 5 grassland site CH-Oe1 (Fig. 2 f,l,r). These values are slightly higher than independent estimates (for grasses and crops obtained for other sites) of noon surface resistances ranging 50-100 sm −1 (Padro, 1996;Coyle et al., 2009;Gerosa et al., 2004;Tuovinen et al., 2004). Tuovinen et al. (2004) reported noon values of approximately 140 sm −1 for a Scots pine forest and 70-140 sm −1 for a Norway spruce forest site (Tuovinen et al., 2001), which compares well with the modelled R c values at the needle-leaved forest site (FI-Hyy; Fig. 2 l). Higher noon values of approximately 250 sm −1 are reported at a Danish Norway spruce site 10 (Mikkelsen et al., 2004). For a Mountain Birch forest noon values of 110-140 sm −1 (Tuovinen et al., 2001) are observed which is slightly lower than the modelled value at the IT-Ro1 site (dominated by broadleaved tree PFT).

Sensitivity analysis
We assess the sensitivity of the modelled O 3 uptake and deposition, represented by F g , F stC , V g , and R c to uncertainty in six weakly constrained variables and parameters of the O 3 deposition scheme (R a , b, r ext ,R gs , G c , and R b ). Fig. 3 a shows   15 for example the results for the boreal needle-leaved forest FI-Hyy. As expected, all uptake/deposition variables, except for the flux ratio (F R ) are negatively correlated to the aerodynamic resistance R a , which describes the level of decoupling of the atmosphere and land surface. Increasing R a decreases the canopy internal O 3 concentration and hence stomatal (F stC ) and total (F g ) deposition as well as the deposition velocity (V g ). The flux ratio F R is slightly positively correlated to changes in R a due to the stronger negative correlation of F stC relative to F g . 20 In decreasing order, but as expected, the level of external leaf-resistance (r ext ), the scaling factor b (eq. 9), the soil resistance (R gs ), and the canopy-scale quasi-laminar layer resistance (R b ) increase R c and consequently reduce F g and V g . Reducing the non-stomatal deposition by increasing r ext , b,R gs , and R b increases the canopy internal O 3 concentration and thus stomatal O 3 uptake (F stC ). The combined effects of a reduction of total deposition F g and an increase of F stC cause a positive correlation of F R to r ext , b,R gs , and R b . 25 Increasing canopy conductance (G c ) increases stomatal O 3 uptake (F stC ) and thereby also increases V g and F g . The increased total O 3 uptake (F g ) decreases the surface resistance to O 3 uptake R c resulting in a negative correlation of R c with G c . The stronger increase in F stC relative to F g results in a positive correlation of F R .
Despite these partial correlations, only changed values for r ext and G c have a notable effect on the predicted fluxes (Fig. 3 b), whereas for the other factors (R a , b, andR gs ), the impact on the simulated fluxes is less than 0.1 % due to a 1 % change in 30 the variables/parameters of the deposition scheme.
The flux ratio F R is very little affected by varying r ext and G c .
Notwithstanding the perturbations, all four O 3 related flux variables show a fairly narrow range of simulated values (Fig. 4).
For all four variables the unperturbed model and the ensemble mean lie on top of each other (see dashed red and yellow line in Fig. 4 a-d). The seasonal course of the surface resistances and fluxes are maintained. The simulations show a strong day to day variability of F stC , which is conserved with different parameter combinations, and which is largely driven by the day-to-day variations in G c and the atmospheric O 3 concentration (see Fig. 4 f and e respectively). Ozone uptake by the leaves reduces the O 3 surface resistance during the growing season such that R c becomes lowest. The cumulative uptake of O 3 (CUO) is lowest at the beginning of the growing season but not zero because the evergreen pine at the Hyytiälä site accumulates O 3 over several 5 years (Fig. 4 f). The CUO increases during the growing season and declines in autumn when a larger fraction of old needles are shed.
The little impact of the perturbations on the simulated O 3 uptake and deposition variables suggests that he calculated O 3 uptake is relatively robust against uncertainties in the parameterisation of some of the lesser known surface properties.
3.4 Regional simulations 10 We used the model to simulate the vegetation productivity, O 3 uptake, and associated ozone damage of plant production over North of 60 • N OCN has the tendency to produce lower estimates of GPP than inferred from the observation-based product, which is particularly pronounced in low productivity mountain regions of Norway and Sweden. It is unclear whether this bias is indicative of a too strong N limitation in the OCN model. Average decadal O 3 concentrations generally increase from Northern to Southern Europe (Fig. 6 a) and with increasing 25 altitude, with local deviations from this pattern in centres of substantial air pollution. The pattern of foliar O 3 uptake differs distinctly from that of the O3 concentrations, showing highest uptake rates in Central, Eastern and parts of Southern Europe ( Fig. 6 b), associated with centres of high rates of simulated gross primary production (Fig. 5 a) and thus canopy conductance.
The cumulative O 3 uptake reaches values of 40-60 mmol m −2 in large parts of Central Europe (Fig. 6 c). The highest accumulation rates of 80-110 mmol m −2 are found in Eastern Europe and parts of Scandinavia as well as in Italy, the Alps and the Bordeaux region. The concentration based exposure index AOT40 (Fig. 6 d) shows a strong north south gradient similar to the O 3 concentration (Fig. 6 a) and is distinctly different to the flux based CUO pattern (Fig. 6 c).
Simulated reduction of mean decadal GPP due to O 3 range from 80-160 g C m −2 yr −1 over large areas of Central, Eastern, and South-eastern Europe (Fig. 7 a) and is generally largest in regions of high productivity. The relative reduction of GPP is fairly consistent across large areas in Europe and averages 6-10 % (Fig. 7 b). Higher reductions in relative terms are found in regions with high cover of C4 PFTs, e.g. Black Sea area. Lower relative reductions are found in Northern and parts of Southern Europe where productivity is low and stomatal O 3 uptake is reduced by e.g. low O 3 concentrations or drought control on stomatal fluxes respectively. Slight increases or strong decreases in relative terms are found in regions with very small productivity like in Northern Africa and the mountainous regions of Scandinavia. A slight increase in GPP might be 5 caused by feedbacks of GPP damage on LAI, canopy conductance and soil moisture content such that e.g. water savings enable a prolonged growing season and thus a slightly higher GPP. Overall, simulated European productivity has been reduced from 10.6 Pg C yr −1 to 9.8 Pg C yr −1 corresponding to a 7.6 % reduction.
The O 3 induced reductions in GPP are associated with a reduction in mean decadal transpiration rates of 8-15 mm yr −1 over large parts of Central and Eastern Europe (Fig. 7 c). These reductions correspond to 3-6 % of transpiration in Central 10 Europe and 6-10 % in Northern Europe. As expected, the relative reductions in transpiration rates are therefore slightly less than for GPP due to the role of aerodynamic resistance in controlling water fluxes in addition to canopy conductance. Very high reductions in transpiration are found in the Eastern Black Sea area associated with strong reductions in GPP and in the mountainous regions of Scandinavia where absolute changes in transpiration are very small. Regionally (in particular in Eastern Spain, Northern Africa and around the Black Sea) lower reductions in transpiration or even slight increases are found (Fig. 7 d). 15 These are related to O 3 induced soil moisture savings during the wet growing season, leading to lower water stress rates during the drier season. The very strong reduction in transpiration West of the Crimean Peninsula are related to the strong reductions in GPP mentioned above. Overall, simulated European mean transpiration has been reduced from 170.4 mm to 163.3 mm corresponding to a 4.2 % reduction.
3.5 Impacts of using the ozone deposition scheme 20 At the FI-Hyy site the canopy O3 concentration, uptake and accumulated uptake (CUO) increases approximately 10-15 % for the D-STO model (non-stomatal depletion of O 3 is zero) and 20-25 % for the ATM model version (canopy O 3 concentration is equal to the atmospheric concentration) compared to the standard deposition scheme (D) used here (Fig. 8a-c and Fig. 12).
The exact values however are site and PFT specific (see Fig. 12 for the CH-Oe1 and IT-Ro1 site).
The regional impact of using the ozone deposition scheme on CUO is shown in Fig. 9. CUO substantially decreases for the 25 D-STO (Fig. 9b) compared to the ATM model (Fig. 9a). Using the standard deposition model D (Fig. 9c)  We extended the terrestrial biosphere model OCN by a scheme to account for the atmosphere-leaf transfer of O 3 in order to better account for air pollution effects on net photosynthesis and hence regional to global water, carbon, and nitrogen cycling. This ozone deposition scheme calculates canopy O 3 concentrations and uptake into the leaves depending on surface conditions and vegetation carbon uptake.
Estimates of the regional damage to annual average GPP (-7.6 %) and transpiration (-4.2 % ) simulated by OCN for 2001-2010 are lower than previously reported estimates. Meta-analyses suggest on average a 11 % (Wittig et al., 2007) and a 21 % (Lombardozzi et al., 2013) reduction of instantaneous photosynthetic rates. However because of carry-over effects this does 5 not necessarily translate directly into reductions in annual GPP. Damage estimates using the Community Land Model (CLM) suggest GPP reductions of 10-25 % in Europe and 10.8 % globally (Lombardozzi et al., 2015). Reductions in transpiration have been estimated as 5-20 % for Europe and 2.2 % globally (Lombardozzi et al., 2015). Lombardozzi et al. (2015) however used fixed reductions of photosynthesis (12-20 %) independent of cumulative O 3 uptake for 2 out of 3 simulated plant types.
Damage was only related to cumulative O 3 uptake for one plant type with a very small slope and hence little increase in damage

Atmosphere-leaf transport of ozone
The sensitivity analysis in Section 3.3 demonstrates that the estimate of canopy conductance (G c ) is crucial for calculating plant ozone uptake, therefore reliable observations to constrain modelled canopy conductance are highly important. The site-level evaluation shows that OCN produces reasonable estimates of simulated gross primary productivity (GPP), canopy conductance, and latent heat fluxes (LE) compared to FLUXNET observations. This agreement has to be seen in the light of the a diverse set 25 of random and systematic errors in the eddy covariance measurements, and derived flux and conductance estimates (Richardson et al., 2012;Knauer et al., 2016). Next to uncertainties about the strength of the aerodynamic coupling between atmosphere and canopy, problems exist at many sites with respect to the energy balance closure (Wilson et al., 2002). Failure to close the energy balance can cause underestimation of sensible and latent heat, as well as an overestimation of available energy, with mean bias of 20 % where the imbalance is greatest during nocturnal periods (Wilson et al., 2002). This imbalance propagates 30 to estimates of canopy conductance, which is inferred from latent and sensible heat fluxes. The energy imbalance furthermore appears to affect estimates of CO 2 uptake and respiration (Wilson et al., 2002). Flux partitioning algorithms which extrapolate night-time ecosystem respiration estimates to daytime introduce an additional potential for bias in the estimation of GPP (Reichstein et al., 2005). Nevertheless, the general good agreement of G c compared to FLUXNET estimates together with the finding that modelled values of key ozone variables are within observed ranges, supports the use of the extended OCN model for determining the effect of air pollution on terrestrial carbon, nitrogen, and water cycling.
A key difference from previous studies is our use of the use of the ozone deposition scheme, which reduces O 3 surface concentrations, and hence also the estimated O 3 uptake and accumulation (see Fig. 9). Accounting for stomatal and nonstomatal deposition in the calculation of the surface O 3 concentrations considerably impacts the estimated plant uptake of O 3 .  (Tuovinen et al., 2009). Flux-based ozone-damage assessment models may overestimate ozone-related damage unless they properly account for non-stomatal O 3 uptake at the surface.
We note that vegetation type and dynamics also impact the stomatal and non-stomatal deposition of O 3 , and hence the 15 calculation of the leaf-level O 3 concentrations. This impedes the use of CTM-derived leaf-level O 3 concentration, as CTM and vegetation specifications may differ strongly. Using the O 3 from the lowest level of the atmosphere reduces this problem, but running a terrestrial biosphere with a fixed atmospheric boundary condition (and not coupled to a atmospheric chemistrytransport model) is still a simplification that prevents biosphere-atmosphere feedbacks and therefore to potential discrepancies between vegetation and CTM model. Not accounting for this feedback and stomatal and non-stomatal O 3 deposition might 20 result in an overestimation of O 3 uptake and hence potential damage in the vegetation model. The deposition scheme in OCN offers the potential to couple vegetation and CTM modelling and is thus a step forward towards coupled atmosphere-vegetation simulations.

Estimating vegetation damage from ozone uptake
A key aspect of ozone damage estimates are the assumed dose-response-relationships, which relate O 3 uptake to plant damage. 25 The use of flux-based relationships is generally thought to improve damage estimates compared to concentration based metrics (e.g. AOT40), since stomatal constraints on O 3 uptake are taken into account, yielding very different spatial patterns of exposure hot spots . Similar to Simpson et al. (2007), we find strongly differing patterns between cumulative O 3 uptake (CUO) and AOT40 in our simulations here (see Fig. 6), where highest exposure is not only found in southern Europe where the O 3 concentration is highest but also in eastern Europe. 30 Several dose-response-relationships exist for biomass or yield damage (LRTAP-Convention, 2010, for an overview), there are few estimates of the likely cause of this damage, i.e. the reduction in net photosynthesis. In this study, the damage relationship to net photosynthesis proposed by Wittig et al. (2007) is used. The major advantage of this relationship is that it has been obtained by meta-analysis of many different tree species and thus might indicate an average response. This relationship is therefore used for all modelled plant functional types. However, a substantial disadvantage is that the meta-analysis implies a damage of 6.16 % at zero accumulated O 3 uptake with a rather minor increase in damage with increasing O 3 uptake. This might be an important factor explaining the lower ozone damage estimates of OCN compared to other terrestrial biosphere models.
In Lombardozzi et al. (2015) also a damage relationship derived from a meta-analysis is used however the disadvantage of predicted ozone damage at zero accumulated O 3 uptake there is even stronger compared to Wittig et al. (2007). Two out of 5 three modelled plant functional types assume -12.5 % and -16,1 % ozone damage at zero accumulated O 3 uptake (broadleaved and needle-leaved species respectively) and the third plant functional type (grass and crop) assumes 19.8 % at zero accumulated O 3 uptake together with a small increase in damage with increasing O 3 uptake (Lombardozzi et al., 2015). An evaluation of the different proposed damage functions implemented in terrestrial biosphere models (e.g. Wittig et al. (2007); Lombardozzi et al. (2015); Sitch et al. (2007)) is necessary to elucidate which are able to e.g. reproduce observed patterns of biomass damage and 10 hence might be suitable to predict regional or global damage estimates. Furthermore new damage relationships for different plant groups would be desirable for use in dynamic vegetation models to improve the ozone damage estimates, for example by ensuring an intercept close to one (zero damage at zero accumulated O 3 ).
The use of a (possibly PFT specific) flux threshold and its magnitude naturally also impacts the CUOY (canopy cumulative O 3 uptake above a threshold of Y nmol m −2 s −1 ) and possible damage estimates . The included damage 15 function (Wittig et al., 2007) is designed for the CUO without a flux threshold (Y = 0). The impacts of using different flux thresholds on regional estimates of O 3 uptake, accumulation and damage are still poorly understood and need further research.
It should be noted that using plant O 3 uptake based on leaf-level O 3 concentrations, as done here, together with empirical ozone-damage functions, where O 3 uptake is calculated from atmospheric O 3 concentrations, introduces a discrepancy. The O 3 uptake rates of the experiments forming the damage relationship however are calculated from mean ozone concentrations 20 e.g. over the exposure period and the respective average stomatal conductances (Wittig et al., 2007) such that the estimated O 3 uptake and cumulated uptake used to derive the damage relationship are coarse approximations and underlie considerable uncertainty. The error introduced in OCN by using leaf-level O 3 concentrations instead of atmospheric concentrations seems small, especially since the use of the leaf-level O 3 concentration is the physiological more appropriate approach.
In the current version of OCN only ozone damage to net photosynthesis is accounted for. Other processes like detoxification 25 of O 3 and injury repair (Wieser and Matyssek, 2007;Ainsworth et al., 2012), stomatal sluggishness (Paoletti and Grulke, 2010) and early senescence (Gielen et al., 2007;Ainsworth et al., 2012) are not accounted for. Decoupling of photosynthesis and stomatal conductance (e.g. through stomatal sluggishness) might impact GPP and transpiration damage estimates and requires further analysis. Accounting for direct impairment of the stomata might reduce the reported reductions in transpiration or even cause an increase compared to simulations with no ozone damage. Reduced carbon gain due to early senescence might impact 30 the growth and biomass accumulation of plants (Gielen et al., 2007;Ainsworth et al., 2012) and ought also be included in terrestrial biosphere models.

5 Conclusion
Estimates of O 3 impacts on plant gross primary productivity vary substantially. This uncertainty in the magnitude of damage and hence the potential impact on the global carbon budget is related to different approaches to model ozone damage. The use of a comparatively detailed ozone deposition scheme that accounts for non-stomatal as well as stomatal deposition, when calculating surface O 3 concentrations substantially affects O 3 uptake in our model. We therefore recommend that non-stomatal 5 O 3 uptake is routinely included in model assessments of ozone damage to obtain a better estimate of ozone uptake and accumulation. We show that O 3 uptake into the stomata is mainly determined by the canopy conductance in the ozone deposition scheme used here. This highlights the importance of reliable modelling of canopy conductances as well as realistic surface O 3 concentrations to obtain as accurate as possible estimates of O 3 uptake which are the basis for plant damage estimates.
Suitable ozone damage relationships to net photosynthesis for different plant groups are essential to relate the accumulated 10 O 3 uptake to plant damage in a model. Mean responses of plant groups similar to commonly modelled plant functional types are also desirable. Only few relationships exist which indicate mean responses of several species e.g. Wittig et al. (2007) and Lombardozzi et al. (2013) which however propose very different relationships. Furthermore, the impact of the plants ability to detoxify O 3 should be considered e.g. by using flux thresholds, as well as the combined effects of O 3 with air pollution (nitrogen deposition) and climate change (elevated CO 2 ) on the plants carbon uptake. where u * is calculated from the wind speed at 10 m height (u 10 ) using the atmospheric resistance calculations of the ORCHIDEE model (Krinner et al., 2005). The wind at 45 m (u 45 ) is approximated by assuming the logarithmic wind profile for neutral atmospheric conditions (Monteith and Unsworth, 2007) due to the lack of information on any other relevant atmospheric properties at 45 m height: where z 0 is the roughness length.

Appendix B: Emissions inventory
Emissions for the EMEP model were derived by merging data from three main sources. Firstly, emissions for 2005 and 2010 were taken from the so-called ECLIPSE database produced by IIASA for various EU Projects and the Task Force on Hemispheric Transport of Air Pollution (Amann et al., 2013;Stohl et al., 2015), although with improved spatial resolution over Europe by making use of the 7 km resolution MACC-2 emissions produced by TNO (Kuenen et al., 2011). For 1990 sions from land-based sources were taken directly from the EMEP database for that year, since 1990 had been the subject of recent review and quality-control (e.g. Mareckova et al., 2013). Emissions between 1990 and 2005 were estimated via linear interpolation between these 2005 and EMEP 1990 values. Emissions prior to 1990 were derived by scaling the EMEP 1990 emissions by the emissions ratios found in the historical data-series of Lamarque et al. (2010).
Emissions of the biogenic hydrocarbon isoprene from vegetation are calculated using the model's land cover and meteoro-10 logical data (Simpson et al., 2012(Simpson et al., , 1999. Emissions of NO from biogenic sources (soils, forest-fires, etc) were set to zero given both their uncertainty and sporadic occurrence. Tests have shown that this approximation has only a small impact on annual deposition totals to the EU area, even for simulations at the start of the 20th century. Volcanic emissions of sulfur dioxide (SO 2 ) were set to a constant value from the year 2010.