Large but decreasing effect of ozone on the European carbon sink

Abstract

The capacity of the terrestrial biosphere to sequester carbon and mitigate climate change is governed by the ability of vegetation to remove emissions of ${\mathrm{CO}}_{\mathrm{2}}$ through photosynthesis. Tropospheric ${\mathrm{O}}_{\mathrm{3}}$, a globally abundant and potent greenhouse gas, is, however, known to damage plants, causing reductions in primary productivity. Despite emission control policies across Europe, background concentrations of tropospheric ${\mathrm{O}}_{\mathrm{3}}$ have risen significantly over the last decades due to hemispheric-scale increases in${\mathrm{O}}_{\mathrm{3}}$ and its precursors. Therefore, plants are exposed to increasing background concentrations, at levels currently causing chronic damage. Studying the impact of ${\mathrm{O}}_{\mathrm{3}}$ on European vegetation at the regional scale is important for gaining greater understanding of the impact of${\mathrm{O}}_{\mathrm{3}}$ on the land carbon sink at large spatial scales. In this work we take a regional approach and update the JULES land surface model using new measurements specifically for European vegetation. Given the importance of stomatal conductance in determining the flux of ${\mathrm{O}}_{\mathrm{3}}$ into plants, we implement an alternative stomatal closure parameterisation and account for diurnal variations in ${\mathrm{O}}_{\mathrm{3}}$ concentration in our simulations. We conduct our analysis specifically for the European region to quantify the impact of the interactive effects of tropospheric ${\mathrm{O}}_{\mathrm{3}}$ and${\mathrm{CO}}_{\mathrm{2}}$ on gross primary productivity (GPP) and land carbon storage across Europe. A factorial set of model experiments showed that tropospheric${\mathrm{O}}_{\mathrm{3}}$ can suppress terrestrial carbon uptake across Europe over the period 1901 to 2050. By 2050, simulated GPP was reduced by 4 to 9 % due to plant ${\mathrm{O}}_{\mathrm{3}}$ damage and land carbon storage was reduced by 3 to 7 %. The combined physiological effects of elevated future ${\mathrm{CO}}_{\mathrm{2}}$ (acting to reduce stomatal opening) and reductions in ${\mathrm{O}}_{\mathrm{3}}$ concentrations resulted in reduced ${\mathrm{O}}_{\mathrm{3}}$ damage in the future. This alleviation of ${\mathrm{O}}_{\mathrm{3}}$ damage by ${\mathrm{CO}}_{\mathrm{2}}$-induced stomatal closure was around 1 to 2 % for both land carbon and GPP, depending on plant sensitivity to ${\mathrm{O}}_{\mathrm{3}}$. Reduced land carbon storage resulted from diminished soil carbon stocks consistent with the reduction in GPP. Regional variations are identified with larger impacts shown for temperate Europe (GPP reduced by 10 to 20 %) compared to boreal regions (GPP reduced by 2 to 8 %). These results highlight that ${\mathrm{O}}_{\mathrm{3}}$ damage needs to be considered when predicting GPP and land carbon, and that the effects of ${\mathrm{O}}_{\mathrm{3}}$ on plant physiology need to be considered in regional land carbon cycle assessments.

Introduction

The terrestrial biosphere absorbs around 30 % of anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions and acts to mitigate climate change (Le Quéré et al., 2015). Early estimates of the European carbon balance suggest a terrestrial carbon sink of between 135 and 205 TgC yr${}^{-\mathrm{1}}$ (Janssens et al., 2003). Schulze et al. (2009) determined a larger carbon sink of 274 TgC yr${}^{-\mathrm{1}}$, and more recent estimates suggest a European terrestrial sink of between 146 and 184 TgC yr${}^{-\mathrm{1}}$ (Luyssaert et al., 2012). The carbon sink capacity of land ecosystems is dominated by the ability of vegetation to sequester carbon through photosynthesis and release it back to the atmosphere through respiration. Therefore, any change in the balance of these fluxes will alter ecosystem source–sink behaviour.

In recent decades much attention has focussed on the effects of rising atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ on vegetation productivity (Ceulemans and Mousseau, 1994; Norby et al., 1999, 2005; Saxe et al., 1998). The Norby et al. (2005) synthesis of Free-Air ${\mathrm{CO}}_{\mathrm{2}}$ Enrichment (FACE) experiments suggests a median stimulation ($\mathrm{23}\pm \mathrm{2}$ %) of forest net primary production (NPP) in response to a doubling of ${\mathrm{CO}}_{\mathrm{2}}$. Similar average increases (20 %) were observed for ${\mathrm{C}}_{\mathrm{3}}$ crops, although this translated into smaller gains in biomass (17 %) and crop yields (13 %) (Long et al., 2006). Little attention, however, has been given to tropospheric ozone (${\mathrm{O}}_{\mathrm{3}}$), a globally abundant air pollutant recognised as one of the most damaging pollutants for forests (Karlsson et al., 2007; Royal-Society, 2008; Simpson et al., 2014b). Tropospheric ${\mathrm{O}}_{\mathrm{3}}$ is a secondary air pollutant formed by photochemical reactions involving carbon monoxide (CO), volatile organic compounds (VOCs), methane (${\mathrm{CH}}_{\mathrm{4}}$) and nitrogen oxides (${\mathrm{NO}}_x$) from both human-made and natural sources, as well as downward transport from the stratosphere and lightning, which is a source of ${\mathrm{NO}}_x$. The phytotoxic effects of ${\mathrm{O}}_{\mathrm{3}}$ exposure are shown to decrease vegetation productivity and biomass, with consequences for terrestrial carbon sequestration (Felzer et al., 2004; Loya et al., 2003; Mills et al., 2011b; Sitch et al., 2007). Few studies, however, consider the simultaneous effects of exposure to both gases, and few Earth system models (ESMs) currently explicitly consider the role of tropospheric ${\mathrm{O}}_{\mathrm{3}}$ in terrestrial carbon dynamics (IPCC, 2013), both of which are important in understanding the carbon sequestration potential of the land surface and future carbon dynamics regionally and globally (Le Quéré et al., 2016; Sitch et al., 2015).

Due to increased anthropogenic precursor emissions over the industrial period, background concentrations of ground-level ${\mathrm{O}}_{\mathrm{3}}$ have risen (Vingarzan, 2004). Background ${\mathrm{O}}_{\mathrm{3}}$ is generally defined as the ${\mathrm{O}}_{\mathrm{3}}$ pollution present in a region that is not attributed to local anthropogenic sources (Vingarzan, 2004). ${\mathrm{O}}_{\mathrm{3}}$ levels at the start of the 20th century are estimated to be around 10 ppb for the Montsouris Observatory site near Paris, data for Arkona on the Baltic coast increased from ca. 15 ppb in the 1950s to 20–27 ppb by the early 1980s, and the Irish coast site Mace Head shows around 40 ppb by the year 2000 (Logan et al., 2012; Parrish et al., 2012). Present-day annual average background ${\mathrm{O}}_{\mathrm{3}}$ concentrations reported in the review of Vingarzan (2004) show ${\mathrm{O}}_{\mathrm{3}}$ concentrations range between approximately 20 and 45 ppb, with the greatest increase occurring since the 1950s. Trends vary from site to site though, even on a decadal basis (Logan et al., 2012; Simpson et al., 2014b), depending, for example, on local/regional trends in precursor (especially ${\mathrm{NO}}_x$) emissions, elevation and exposure to long-range transport of ${\mathrm{O}}_{\mathrm{3}}$. Nevertheless, there is some indication that background ${\mathrm{O}}_{\mathrm{3}}$ levels over the mid-latitudes of the Northern Hemisphere have continued to rise at a rate of approximately 0.5–2 % per year, although not uniformly (Vingarzan, 2004). As a result of controls on precursor emissions in Europe and North America, peak ${\mathrm{O}}_{\mathrm{3}}$ concentrations in these regions have decreased or stabilised over recent decades (Cooper et al., 2014; Logan et al., 2012; Parrish et al., 2012; Simpson et al., 2014b). Nevertheless, climate change may increase the frequency of weather events conducive to peak ${\mathrm{O}}_{\mathrm{3}}$ incidents in the future (e.g. summer droughts and heatwaves; Sicard et al., 2013), and may increase biogenic emissions of the ${\mathrm{O}}_{\mathrm{3}}$ precursors isoprene and ${\mathrm{NO}}_x$, although such impacts are subject to great uncertainty (Simpson et al., 2014b; Young et al., 2009, 2013). Intercontinental transport of air pollution from regions such as Asia are thought to contribute substantially to rising background ${\mathrm{O}}_{\mathrm{3}}$ concentrations over the last decades (Cooper et al., 2010; Verstraeten et al., 2015). Northern Hemisphere background concentrations of ${\mathrm{O}}_{\mathrm{3}}$ are now close to established levels for impacts on human health and the terrestrial environment (Royal-Society, 2008). Therefore, although peak ${\mathrm{O}}_{\mathrm{3}}$ concentrations are in decline across Europe, plants are exposed to increasing background levels, at levels currently causing chronic damage (Mills et al., 2011b). Intercontinental transport means future ${\mathrm{O}}_{\mathrm{3}}$ concentrations in Europe will be partly dependent on how ${\mathrm{O}}_{\mathrm{3}}$ precursor emissions evolve globally (Auvray and Bey, 2005; Derwent et al., 2015).

Rising background ${\mathrm{O}}_{\mathrm{3}}$ concentrations impact agricultural yields and nutritional quality of major crops (Ainsworth et al., 2012; Avnery et al., 2011), with consequences for global food security (Tai et al., 2014). Increasing background levels of ${\mathrm{O}}_{\mathrm{3}}$ are damaging to ecosystem health and reduce the global land carbon sink (Arneth et al., 2010; Sitch et al., 2007). Reduced uptake of carbon by plant photosynthesis due to ${\mathrm{O}}_{\mathrm{3}}$ damage allows more ${\mathrm{CO}}_{\mathrm{2}}$ to remain in the atmosphere. This effect of ${\mathrm{O}}_{\mathrm{3}}$ on plant physiology represents an additional climate warming to the direct radiative forcing of ${\mathrm{O}}_{\mathrm{3}}$, a potent greenhouse gas (Collins et al., 2010; Sitch et al., 2007), the magnitude of which, however, remains highly uncertain (IPCC, 2013).

Dry deposition of ${\mathrm{O}}_{\mathrm{3}}$ to terrestrial surfaces, primarily uptake by stomata on plant foliage and deposition on external surfaces of vegetation (Fowler et al., 2001, 2009), is a large sink for ground-level ${\mathrm{O}}_{\mathrm{3}}$ (Wild, 2007; Young et al., 2013). On entry to sub-stomatal spaces, ${\mathrm{O}}_{\mathrm{3}}$ reacts with other molecules to form reactive oxygen species (ROS). Plants can tolerate a certain level of ${\mathrm{O}}_{\mathrm{3}}$ depending on their capacity to scavenge and detoxify the ROS (Ainsworth et al., 2012). Above this critical level, long-term chronic ${\mathrm{O}}_{\mathrm{3}}$ exposure reduces plant photosynthesis and biomass accumulation (Ainsworth, 2008, 2012; Matyssek et al., 2010a; Wittig et al., 2007, 2009), either directly through effects on photosynthetic machinery such as reduced Rubisco content (Ainsworth et al., 2012; Wittig et al., 2009) and/or indirectly by reduced stomatal conductance ($g_{\mathrm{s}})$ (Kitao et al., 2009; Wittig et al., 2007), alters carbon allocation to different pools (Grantz et al., 2006; Wittig et al., 2009), accelerates leaf senescence (Ainsworth, 2008; Nunn et al., 2005; Wittig et al., 2009) and changes plant susceptibility to biotic stress factors (Karnosky et al., 2002; Percy et al., 2002).

The response of plants to ${\mathrm{O}}_{\mathrm{3}}$ is very wide-ranging as reported in the literature from different field studies. The Wittig et al. (2007) meta-analysis of temperate and boreal tree species showed that raised ${\mathrm{O}}_{\mathrm{3}}$ concentrations significantly reduced leaf-level light-saturated net photosynthetic uptake ($-\mathrm{19}$ %, range: $-\mathrm{3}$ to $-\mathrm{28}$ % at a mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 85 ppb) and $g_{\mathrm{s}}$ ($-\mathrm{10}$ %, range: $+\mathrm{5}$ to $-\mathrm{23}$ % at a mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 91 ppb) in both broadleaf and needleleaf tree species. In the Feng et al. (2008) meta-analysis of wheat, ${\mathrm{O}}_{\mathrm{3}}$ reduced aboveground biomass ($-\mathrm{18}$ % at a mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 70 ppb) photosynthetic rate ($-\mathrm{20}$ % at a mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 73 ppb) and $g_{\mathrm{s}}$ ($-\mathrm{22}$ % at a mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 79 ppb). One of few long-term field-based ${\mathrm{O}}_{\mathrm{3}}$ exposure studies (AspenFACE) showed that after 11 years of exposing mature trees to ${\mathrm{O}}_{\mathrm{3}}$ (mean ${\mathrm{O}}_{\mathrm{3}}$ concentration of 46 ppb), ${\mathrm{O}}_{\mathrm{3}}$ decreased ecosystem carbon content ($-\mathrm{9}$ %) and decreased NPP ($-\mathrm{10}$ %), although the ${\mathrm{O}}_{\mathrm{3}}$ effect decreased through time (Talhelm et al., 2014). Zak et al. (2011) showed this was partly due to a shift in community structure as ${\mathrm{O}}_{\mathrm{3}}$-tolerant species, competitively inferior in low-${\mathrm{O}}_{\mathrm{3}}$ environments, outcompeted ${\mathrm{O}}_{\mathrm{3}}$-sensitive species. Gross primary productivity (GPP) was reduced ($-\mathrm{12}$ to $-\mathrm{19}$ %) at two Mediterranean ecosystems exposed to ${\mathrm{O}}_{\mathrm{3}}$ (ranging between 20 and 72 ppb across sites and through the year) studied by Fares et al. (2013). Biomass of mature beech trees was reduced ($-\mathrm{44}$ %) after 8 years of exposure to ${\mathrm{O}}_{\mathrm{3}}$ ($\sim$ 150 ppb) (Matyssek et al., 2010a). After 5 years of ${\mathrm{O}}_{\mathrm{3}}$ exposure (ambient $+\mathrm{20}$ to $+\mathrm{40}$ ppb) in a semi-natural grassland, annual biomass production was reduced ($-\mathrm{23}$ %), and in a Mediterranean annual pasture ${\mathrm{O}}_{\mathrm{3}}$ exposure significantly reduced total aboveground biomass (up to $-\mathrm{25}$ %) (Calvete-Sogo et al., 2014). However, these were empirical studies at individual sites, and these focus on ${\mathrm{O}}_{\mathrm{3}}$ effects on plant physiology and productivity but do not quantify the impact on the land carbon sink. Modelling studies are needed to scale site observations to the regional and global scales. Models generally suggest that plant productivity and carbon sequestration will decrease with ${\mathrm{O}}_{\mathrm{3}}$ pollution, though the magnitudes vary. For example, based on a limited dataset to parameterise plant ${\mathrm{O}}_{\mathrm{3}}$ damage for a global set of plant functional types (PFTs), Sitch et al. (2007) predicted a decline in global GPP of 14 to 23 % by 2100. A second study by Lombardozzi et al. (2015) predicted a 10.8 % decrease in present-day (2002–2009) GPP globally. Here we take a regional approach and take advantage of the latest measurements showing changes in plant productivity with accumulated exposure to ${\mathrm{O}}_{\mathrm{3}}$ specifically for a range of European vegetation from different regions (CLRTAP, 2017) with which to calibrate the JULES model for plant sensitivity to ${\mathrm{O}}_{\mathrm{3}}$, and we conduct our analysis specifically for the European region.

Understanding the response of plants to elevated tropospheric ${\mathrm{O}}_{\mathrm{3}}$ is challenged by the large variation in ${\mathrm{O}}_{\mathrm{3}}$ sensitivity both within and among species (Karnosky et al., 2007; Kubiske et al., 2007; Wittig et al., 2009). Additionally, other environmental stresses that affect stomatal behaviour will affect the rate of ${\mathrm{O}}_{\mathrm{3}}$ uptake and therefore the response to ${\mathrm{O}}_{\mathrm{3}}$ exposure, such as high temperature, drought and changing concentrations of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ (Mills et al., 2016; Fagnano et al., 2009; Kitao et al., 2009; Löw et al., 2006). Increasing concentrations of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, for example, are suggested to provide some protection against ${\mathrm{O}}_{\mathrm{3}}$ damage by causing stomata to close (Harmens et al., 2007; Wittig et al., 2007); however the long-term effects of ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation on plant growth and carbon storage remain uncertain (Baig et al., 2015; Ciais et al., 2013). Further, in some studies, stomata have been shown to respond sluggishly, losing their responsiveness to environmental stimuli with exposure to ${\mathrm{O}}_{\mathrm{3}}$, which can lead to higher ${\mathrm{O}}_{\mathrm{3}}$ uptake, increased water loss and therefore greater vulnerability to environmental stresses such as drought (Mills et al., 2009, 2016; Paoletti and Grulke, 2010; Wilkinson and Davies, 2009).

Given the critical role $g_{\mathrm{s}}$ plays in the uptake of both ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$, we use an alternative representation and parameterisation of $g_{\mathrm{s}}$ in JULES by implementing the Medlyn et al. (2011) $g_{\mathrm{s}}$ formulation. This model is based on the optimal theory of stomatal behaviour and has advantages over the current JULES $g_{\mathrm{s}}$ formulation of Jacobs (1994) including (i) a single parameter ($g_{\mathrm{1}})$ compared to two parameters in Jacobs (1994), (ii) the $g_{\mathrm{1}}$ parameter is related to the water-use strategy of vegetation and is easier to parameterise with commonly measured leaf- or canopy-level observations of photosynthesis, $g_{\mathrm{s}}$, and humidity, and (iii) values of $g_{\mathrm{1}}$ are available for many different PFTs derived from a global dataset of leaf-level measurements (Lin et al., 2015).

The main objective of this work is to assess the impact of historical and projected (1901 to 2050) changes in tropospheric ${\mathrm{O}}_{\mathrm{3}}$ and atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration on predicted GPP and the land carbon sink for Europe. These are the two greenhouse gases that directly affect plant photosynthesis and $g_{\mathrm{s}}$. We use a factorial suite of model experiments, using the Joint UK Land Environment Simulator (JULES) (Best et al., 2011; Clark et al., 2011) and the land surface model of the UK Earth System Model (UKESM) (Collins et al., 2011) to simulate plant ${\mathrm{O}}_{\mathrm{3}}$ uptake and damage and to investigate the impact of both ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ on plant water use and carbon uptake. In this work, the JULES model is recalibrated using the latest observations of vegetation sensitivity to ${\mathrm{O}}_{\mathrm{3}}$, with the addition of a separate parameterisation for temperate and boreal regions versus the Mediterranean. The ${\mathrm{O}}_{\mathrm{3}}$ sensitivity of each PFT in JULES was recalibrated for both a high and low sensitivity to account for uncertainty in the ${\mathrm{O}}_{\mathrm{3}}$ response, in part due to the observed variation in ${\mathrm{O}}_{\mathrm{3}}$ sensitivity among species. This includes ${\mathrm{O}}_{\mathrm{3}}$ sensitivities for agricultural crops (wheat – high sensitivity) versus natural grassland (low sensitivity), with separate sensitivities for Mediterranean grasslands. For forests JULES is parameterised with ${\mathrm{O}}_{\mathrm{3}}$ sensitivities for broadleaf and needleleaf trees (with a high and low ${\mathrm{O}}_{\mathrm{3}}$ sensitivity for both), with separate sensitivities (high and low) for Mediterranean broadleaf species. We make a separate distinction for the Mediterranean region where possible because the work of Büker et al. (2015) showed that the sensitivity of dominant Mediterranean trees to ${\mathrm{O}}_{\mathrm{3}}$ is different from temperate species. In addition, we introduce an alternative $g_{\mathrm{s}}$ scheme into JULES as described above. JULES is forced with spatially varying daily ${\mathrm{O}}_{\mathrm{3}}$ concentrations from a high-resolution atmospheric chemistry model for Europe that are disaggregated to hourly concentrations; therefore our simulations account for diurnal variations in ${\mathrm{O}}_{\mathrm{3}}$ concentration and ${\mathrm{O}}_{\mathrm{3}}$ responses, allowing for improved estimates of ${\mathrm{O}}_{\mathrm{3}}$ uptake by vegetation. We do not attempt to make a full assessment of the carbon cycle of Europe, instead we target ${\mathrm{O}}_{\mathrm{3}}$ damage, which is currently a missing component in earlier carbon cycle assessments (Le Quéré et al., 2018; Sitch et al., 2015). To this end, we prescribe changing ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ concentrations from 1901 to 2050 but using a fixed pre-industrial climate. We acknowledge the use of a “fixed” pre-industrial climate omits the additional uncertainty of the interaction between climate change and $g_{\mathrm{s}}$, which will affect the rate of ${\mathrm{O}}_{\mathrm{3}}$ uptake and therefore ${\mathrm{O}}_{\mathrm{3}}$ concentrations. In addition, using uncoupled chemistry and climate is a further source of uncertainty. To understand the impact of these complex feedback mechanisms is an important area for future work, but in the current study our aim is to isolate the physiological response of plants to both ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ and determine the sensitivity of predicted GPP and the land carbon sink to this process, as the impact of ${\mathrm{O}}_{\mathrm{3}}$ on the land carbon sink currently remains largely unknown at large spatial scales for Europe.

Methods

Representation of ${\mathrm{O}}_{\mathrm{3}}$ effects in JULES

JULES calculates the land–atmosphere exchanges of heat, energy, mass, momentum and carbon on a sub-daily time step, and includes a dynamic vegetation model (Best et al., 2011; Clark et al., 2011; Cox, 2001). This work uses JULES version 3.3 (http://www.jchmr.org, last access: 10 September 2017) at a 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ spatial resolution and hourly model time step; the spatial domain is shown in Fig. S1 in the Supplement. JULES has a multilayer canopy radiation interception and photosynthesis scheme (10 layers in this instance) that accounts for direct and diffuse radiation, sun fleck penetration through the canopy, inhibition of leaf respiration in the light and change in photosynthetic capacity with depth into the canopy (Clark et al., 2011; Mercado et al., 2009). Soil water content also affects the rate of photosynthesis and $g_{\mathrm{s}}$. It is modelled using a dimensionless soil water stress factor, $\mathit{\beta }$, which is related to the mean soil water concentration in the root zone, and the soil water contents at the critical and wilting points (Best et al., 2011).

To simulate the effects of stomatal ${\mathrm{O}}_{\mathrm{3}}$ deposition on vegetation productivity and water use, JULES uses the flux-gradient approach of Sitch et al. (2007), modified to include non-stomatal deposition following Tuovinen et al. (2009). A similar approach is taken by Franz et al. (2017) in the OCN model; however plant ${\mathrm{O}}_{\mathrm{3}}$ damage is a function of accumulated ${\mathrm{O}}_{\mathrm{3}}$ exposure over time. In JULES, plant ${\mathrm{O}}_{\mathrm{3}}$ damage is instantaneous because the impact of cumulative ${\mathrm{O}}_{\mathrm{3}}$ exposure on plant productivity has already been calibrated with observations (described below). JULES uses a coupled model of $g_{\mathrm{s}}$ and photosynthesis; the potential net photosynthetic rate ($A_{\mathrm{p}}$, mol ${\mathrm{CO}}_{\mathrm{2}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is modified by an “${\mathrm{O}}_{\mathrm{3}}$ uptake” factor ($F$ the fractional reduction in photosynthesis), so that the actual net photosynthesis ($A_{\mathrm{net}}$, mol ${\mathrm{CO}}_{\mathrm{2}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is given by Eq. (1) (Clark et al., 2011; Sitch et al., 2007). Because of the relationship between these two fluxes, the direct effect of ${\mathrm{O}}_{\mathrm{3}}$ damage on photosynthetic rate also leads to a reduction in $g_{\mathrm{s}}$. An alternative approach was taken by Lombardozzi et al. (2012) in the CLM model in which photosynthesis and $g_{\mathrm{s}}$ are decoupled, so that ${\mathrm{O}}_{\mathrm{3}}$ exposure affects carbon assimilation and transpiration independently. In JULES, changes in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration also affect photosynthetic rate and $g_{\mathrm{s}}$; consequently the interactive effects of changing concentrations of both ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ are allowed for. $$A_{\mathrm{net}}=A_{\mathrm{p}}F$$

The ${\mathrm{O}}_{\mathrm{3}}$ uptake factor ($F$) is defined as $$F=\mathrm{1}-a\cdot \max \left[F_{{\mathrm{O}}_{\mathrm{3}}}-F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}},\mathrm{0.0}\right].$$

$F_{{\mathrm{O}}_{\mathrm{3}}}$ is the instantaneous leaf uptake of ${\mathrm{O}}_{\mathrm{3}}$ (nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$, $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$ is a PFT-specific threshold for ${\mathrm{O}}_{\mathrm{3}}$ damage (nmol m${}^{-\mathrm{2}}$ PLA s${}^{-\mathrm{1}}$, projected leaf area) and “$a$” is a PFT-specific parameter representing the fractional reduction of photosynthesis with ${\mathrm{O}}_{\mathrm{3}}$ uptake by leaves. Following Tuovinen et al. (2009), the flux of ${\mathrm{O}}_{\mathrm{3}}$ through stomata, $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$, is represented as follows: $$F_{{\mathrm{O}}_{\mathrm{3}}}={\mathrm{O}}_{\mathrm{3}}\left({\displaystyle \frac{g_{\mathrm{b}}\left(\frac{g_{\mathrm{l}}}{K_{{\mathrm{O}}_{\mathrm{3}}}}\right)}{g_{\mathrm{b}}+\left(\frac{g_{\mathrm{l}}}{K_{{\mathrm{O}}_{\mathrm{3}}}}\right)+g_{\mathrm{ext}}}}\right).$$

${\mathrm{O}}_{\mathrm{3}}$ is the molar concentration of ${\mathrm{O}}_{\mathrm{3}}$ at reference (canopy) level (nmol m${}^{-\mathrm{3}})$, $g_{\mathrm{b}}$ is the leaf-scale boundary layer conductance (m s${}^{-\mathrm{1}}$, Eq. 3b), $g_{\mathrm{l}}$ is the leaf conductance for water (m s${}^{-\mathrm{1}})$, $K_{{\mathrm{O}}_{\mathrm{3}}}$ accounts for the different diffusivity of ozone to water vapour and takes a value of 1.51 after Massman (1998), and $g_{\mathrm{ext}}$ is the leaf-scale non-stomatal deposition to external plant surfaces (m s${}^{-\mathrm{1}})$, which takes a constant value of 0.0004 m s${}^{-\mathrm{1}}$ after Tuovinen et al. (2009). The leaf-level boundary layer conductance ($g_{\mathrm{b}})$ is calculated as in Tuovinen et al. (2009). $$g_{\mathrm{b}}=\mathit{\alpha }{\mathrm{Ld}}^{-\mathrm{1}/\mathrm{2}}{U}^{-\mathrm{1}/\mathrm{2}}$$

$\mathit{\alpha }$ is a constant (0.0051 m s${}^{-\mathrm{1}/\mathrm{2}})$, Ld is the cross-wind leaf dimension (m) defined per PFT as 0.05 for trees, 0.02 for grasses (${\mathrm{C}}_{\mathrm{3}}$ and ${\mathrm{C}}_{\mathrm{4}}$), and 0.04 for shrubs and $U$ is wind speed at canopy height (m s${}^{-\mathrm{1}})$. The rate of ${\mathrm{O}}_{\mathrm{3}}$ uptake is dependent on $g_{\mathrm{s}}$, which is dependent on photosynthetic rate. Given $g_{\mathrm{s}}$ is a linear function of photosynthetic rate in JULES (Clark et al., 2011), from Eq. (1) it follows that $$g_{\mathrm{s}}=g_{\mathrm{l}}F.$$

The ${\mathrm{O}}_{\mathrm{3}}$ flux to stomata, $F_{{\mathrm{O}}_{\mathrm{3}}}$, is calculated at leaf level and then scaled to each canopy layer differentiating sunlit and shaded leaf photosynthesis, and is finally summed up to the canopy level. Because the photosynthetic capacity, photosynthesis and therefore $g_{\mathrm{s}}$ decline with depth into the canopy, this in turn affects ${\mathrm{O}}_{\mathrm{3}}$ uptake, with the top leaf level contributing most to the total ${\mathrm{O}}_{\mathrm{3}}$ flux and the lowest level contributing least.

Calibration of ${\mathrm{O}}_{\mathrm{3}}$ uptake model

Here we use the latest literature on flux-based ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships derived from observed field data across Europe (CLRTAP, 2017) to determine the key PFT-specific ${\mathrm{O}}_{\mathrm{3}}$ sensitivity parameters in JULES ($a$ and $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$). Synthesis of information expressed as ${\mathrm{O}}_{\mathrm{3}}$ flux-based dose–response relationships derived from field experiments is carried out by the United Nations Convention on Long-Range Transboundary Air Pollution (CLRTAP), this information is then used as a policy tool to inform emission reduction strategies in Europe to improve air quality (CLRTAP, 2017; Mills et al., 2011a). Derivation of ${\mathrm{O}}_{\mathrm{3}}$ flux-based dose–response relationships for different vegetation types uses the accumulated stomatal ${\mathrm{O}}_{\mathrm{3}}$ flux above a threshold (often referred to as the phytotoxic ${\mathrm{O}}_{\mathrm{3}}$ dose above a threshold of “$y$” i.e. POD${}_y$) as the dose metric, and the percentage change in biomass as the response metric (Emberson et al., 2007; Karlsson et al., 2007). We use these observation-based ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships to calibrate each JULES PFT for sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ using available relationships for the closest matching vegetation type. For JULES, $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$ is the threshold for ${\mathrm{O}}_{\mathrm{3}}$ damage, and values for this parameter are taken from the ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships as the POD${}_y$ value (see CLRTAP, 2017, and Büker et al., 2015, for derivation of POD${}_y$ values). The actual sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ is determined by the slope of the ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationship, i.e. how much biomass changes with accumulated stomatal uptake of ${\mathrm{O}}_{\mathrm{3}}$ above the damage threshold; this relates to the parameter $a$ in JULES. The parameter $a$ is a PFT-specific parameter representing the fractional reduction of photosynthesis with ${\mathrm{O}}_{\mathrm{3}}$ uptake by leaves. Values for this parameter are found for each PFT by running JULES with different values of $a$, which alter the instantaneous photosynthetic rate, but then calculating the accumulated stomatal flux of ${\mathrm{O}}_{\mathrm{3}}$ and the change in productivity until the slope of this relationship produced by the JULES simulations matches that of the ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships derived from observations. Essentially we calibrate each JULES PFT for sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ by reproducing the observation-based ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships.

Each PFT was calibrated for high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity to account for uncertainty in the sensitivity of different plant species to ${\mathrm{O}}_{\mathrm{3}}$, using the approach of Sitch et al. (2007). Therefore, when using our results to assess the impact of ${\mathrm{O}}_{\mathrm{3}}$ at the land surface, we are able to provide a range in our estimates to help address some of the uncertainty in the ${\mathrm{O}}_{\mathrm{3}}$ response of different vegetation types. In addition, where possible owing to available data, a distinction was made for Mediterranean regions. This was because the work of Büker et al. (2015) showed that different ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships are needed to describe the ${\mathrm{O}}_{\mathrm{3}}$ sensitivity of dominant Mediterranean trees. For the ${\mathrm{C}}_{\mathrm{3}}$ herbaceous PFT, the dominant land cover type across the European domain in this study (Fig. S2), the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity was calibrated against observations for wheat to give a representation of agricultural regions and wheat is one of the most sensitive grasses to ${\mathrm{O}}_{\mathrm{3}}$ (Fig. S3, Table S1). For the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity JULES was calibrated against the dose–response function for natural grassland to give a representation of natural grassland and this vegetation has a much lower sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ damage; for the Mediterranean region we used a function for Mediterranean natural grasslands, all taken from CLRTAP (2017) (Fig. S3, Table S1 in the Supplement). Tree–shrub PFTs were calibrated against observed ${\mathrm{O}}_{\mathrm{3}}$ dose–response functions for the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity: broadleaf trees (temperate–boreal): birch–beech dose–response relationship; broadleaf trees (Mediterranean): deciduous oak dose–response relationship; needleleaf trees: Norway spruce dose–response relationship; shrubs: birch–beech dose–response relationship; all from CLRTAP (2017) (Fig. S3, Table S1). Data on ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships for different vegetation types is very limited; therefore for the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity calibration for trees–shrubs we assumed a 20 % decrease in sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ based on the difference in sensitivity between high- and low-sensitivity tree species in the Karlsson et al. (2007) study. Due to limitations in data availability, the shrub parameterisation uses the observed dose–response functions for broadleaf trees. Similarly, the parameterisation for ${\mathrm{C}}_{\mathrm{4}}$ herbaceous uses the observed dose–responses for ${\mathrm{C}}_{\mathrm{3}}$ herbaceous; however the fractional cover of ${\mathrm{C}}_{\mathrm{4}}$ herbs across Europe is low (Fig. S2), so this assumption affects a very small percentage of land cover.

To calibrate the JULES ${\mathrm{O}}_{\mathrm{3}}$ uptake model, JULES was run across Europe forced using the WFDEI observational climate dataset (Weedon, 2013) at 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ spatial and 3 h temporal resolution. JULES uses interpolation to disaggregate the forcing data down from 3 h to an hourly model time step. The model was spun up over the period from 1979 to 1999 with a fixed atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration of 368.33 ppm (1999 value from Mauna Loa observations; Tans and Keeling, 2014). Zero tropospheric ozone concentration was assumed for the control simulation. For the simulations with ${\mathrm{O}}_{\mathrm{3}}$, spin-up used spatially explicit fields of present-day ${\mathrm{O}}_{\mathrm{3}}$ concentration produced using the UK Chemistry and Aerosol (UKCA) model with standard chemistry from the run evaluated by O'Connor et al. (2014). A fixed land cover map was used based on IGBP (International Geosphere-Biosphere Programme) land cover classes (IGBP-DIS). Therefore as the vegetation distribution was fixed and the calibration was not looking at carbon stores, a short spin-up time was adequate to equilibrate soil temperature and soil moisture. JULES was then run for the year 2000 with a corresponding ${\mathrm{CO}}_{\mathrm{2}}$ concentration of 369.52 ppm (from Mauna Loa observations; Tans and Keeling, 2014) and monthly fields of spatially explicit tropospheric ${\mathrm{O}}_{\mathrm{3}}$ (O'Connor et al., 2014) as necessary.

Calibration was performed using four simulations: with (i) zero tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentration, as the control simulation (control), (ii) tropospheric ${\mathrm{O}}_{\mathrm{3}}$ at the current ambient concentration (${\mathrm{O}}_{\mathrm{3}}$), (iii) ambient $+\mathrm{20}$ ppb (${\mathrm{O}}_{\mathrm{3}}$ $+$ 20) and (iv) ambient $+\mathrm{40}$ ppb (${\mathrm{O}}_{\mathrm{3}}$ $+$ 40). The different ${\mathrm{O}}_{\mathrm{3}}$ simulations (i.e. ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{O}}_{\mathrm{3}}$ $+$ 20 and ${\mathrm{O}}_{\mathrm{3}}$ $+$ 40) were used to capture the range of ${\mathrm{O}}_{\mathrm{3}}$ conditions in the data used in the observation-based ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships used in this study for calibration. Often data were from experiments using artificially manipulated conditions of ambient $+\mathrm{40}$ ppb ${\mathrm{O}}_{\mathrm{3}}$, for example. For each JULES ${\mathrm{O}}_{\mathrm{3}}$ simulation, the value of $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$ was taken from the vegetation-specific ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationship as the threshold ${\mathrm{O}}_{\mathrm{3}}$ concentration above which damage to vegetation occurs. An initial estimate of the parameter $a$ was used. Then for each PFT and each simulation, hourly estimates of NPP (our proxy for biomass – although not identical, they are related) and ${\mathrm{O}}_{\mathrm{3}}$ uptake in excess of $F_{{\mathrm{O}}_{\mathrm{3}}\mathrm{crit}}$ were accumulated over a PFT-dependent accumulation period. The accumulation periods were $\sim$ 6 months for broadleaf trees and shrubs, all year for needleleaf trees and $\sim$ 3 months for herbaceous species, through the growing season, following guidelines in CLRTAP (2017). Additionally, in accordance with the methods used in CLRTAP (2017) that describe how the ${\mathrm{O}}_{\mathrm{3}}$ dose–response relationships are derived from observations, we use the stomatal ${\mathrm{O}}_{\mathrm{3}}$ flux per projected leaf area to top canopy sunlit leaves. The percentage change in total NPP was calculated for each ${\mathrm{O}}_{\mathrm{3}}$ simulation and plotted against the cumulative uptake of ${\mathrm{O}}_{\mathrm{3}}$ over the PFT-specific accumulation period. The linear regression of this relationship was calculated, and slope and intercept were compared against the slope and intercept of the observed dose–response relationships. Values of the parameter $a$ were adjusted, and the procedure was repeated until the linear regression through the simulation points matched that of the observations (Fig. S3, Table S1).

Representation of stomatal conductance and site-level evaluation

In JULES, $g_{\mathrm{s}}$ (m s${}^{-\mathrm{1}})$ is represented following the closure proposed by Jacobs (1994): $$g_{\mathrm{s}}=\mathrm{1.6}R{T}_{\mathrm{l}}{\displaystyle \frac{A_{\mathrm{net}}\mathit{\beta }}{c_a-c_i}}.$$ In this parameterisation, $c_i$ is unknown and in the default JULES model is calculated as in Eq. (6), hereafter called JAC: $$c_i=\left(c_a-c_{\ast }\right)f_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}{q}_{\mathrm{crit}}}}\right)+c_{\ast }.$$

$\mathit{\beta }$ is a soil moisture stress factor, the factor 1.6 accounts for $g_{\mathrm{s}}$ being the conductance for water vapour rather than ${\mathrm{CO}}_{\mathrm{2}}$, $R$ is the universal gas constant (J K${}^{-\mathrm{1}}$ mol${}^{-\mathrm{1}})$, $T_{\mathrm{l}}$ is the leaf surface temperature (K), $c_a$ and $c_i$ (both Pa) are the leaf surface and internal ${\mathrm{CO}}_{\mathrm{2}}$ partial pressures, respectively, $c_{\ast }$ (Pa) is the ${\mathrm{CO}}_{\mathrm{2}}$ photorespiration compensation point, $\mathrm{d}q$ is the humidity deficit at the leaf surface (kg kg${}^{-\mathrm{1}})$, $\mathrm{d}{q}_{\mathrm{crit}}$ (kg kg${}^{-\mathrm{1}})$ and $f_{\mathrm{0}}$ are PFT-specific parameters representing the critical humidity deficit at the leaf surface, and the leaf internal-to-atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ ratio ($c_i/c_a)$ at the leaf specific humidity deficit (Best et al., 2011); values are shown in Table S1.

In this work, we replace Eq. (6) with the closure described in Medlyn et al. (2011), using the key PFT-specific model parameter $g_{\mathrm{1}}$ (kPa${}^{\mathrm{0.5}}$), and d$q$ is expressed in kilopascals, shown in Eq. (7), hereafter called MED: $$c_i=c_a\left({\displaystyle \frac{g_{\mathrm{1}}}{g_{\mathrm{1}}+\sqrt{\mathrm{d}q}}}\right).$$

PFT-specific values of the $g_{\mathrm{1}}$ parameter were derived for European vegetation from the database of Lin et al. (2015) and are shown in Table S1. The $g_{\mathrm{1}}$ parameter represents the sensitivity of $g_{\mathrm{s}}$ to the assimilation rate, i.e. plant water use efficiency, and was derived as in Lin et al. (2015) by fitting the Medlyn et al. (2011) model to observations of $g_{\mathrm{s}}$, photosynthesis, and vapour pressure deficit (VPD), assuming an intercept of zero.

The impact of $g_{\mathrm{s}}$ model formulation (JAC versus MED) on simulated water, ${\mathrm{O}}_{\mathrm{3}}$, carbon and energy fluxes is compared for two contrasting grid points – wet (low soil moisture stress) and dry (high soil moisture stress) in the European domain. JULES was spun up for 20 years (1979–1999) at two grid points in central Europe representing a wet (low soil moisture stress; lat: 48.25, long: 5.25) and a dry site (high soil moisture stress; lat: 38.25, long: $-\mathrm{7.75}$). The modelled soil moisture stress factor (fsmc) at the wet site ranged from 0.8 to 1.0 over the year 2000 (1.0 indicates no soil moisture stress), and at the dry site fsmc steadily declined from 0.8 at the start of the year to 0.25 by the end of the summer. The WFDEI meteorological forcing dataset (Weedon, 2013), along with atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration for the year 1999 (368.33 ppm), was used and either no ${\mathrm{O}}_{\mathrm{3}}$ (i.e. the ${\mathrm{O}}_{\mathrm{3}}$ damage model was switched off) for the control simulations or spatially explicit fields of present-day ${\mathrm{O}}_{\mathrm{3}}$ concentration produced using the UK Chemistry and Aerosol (UKCA) model from the run evaluated by O'Connor et al. (2014) for the simulations with ${\mathrm{O}}_{\mathrm{3}}$ were used. Following the spin-up period, JULES was run for 1 year (2000) with corresponding atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration, and tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentrations as described above. The control and ${\mathrm{O}}_{\mathrm{3}}$ simulations were performed for both JAC and MED model formulations. Land cover for the spin-up and main run was fixed at 20 % for each PFT. For the simulations including ${\mathrm{O}}_{\mathrm{3}}$ damage, the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity parameterisation was used. The difference between these simulations was used to assess the impact of $g_{\mathrm{s}}$ model formulation on the leaf-level fluxes of carbon and water. We calculate and report (results Sect. 3.1) the difference in mean annual water use that results from the above simulations using the different $g_{\mathrm{s}}$ models. For each day of the simulation we calculate the percentage difference in water use between the two simulations; we then calculate the mean and standard deviation over the year to give the annual mean leaf-level water use.

Site-level evaluation of the two $g_{\mathrm{s}}$ models compared to FLUXNET observations was carried out to evaluate the seasonal cycles of latent and sensible heat using the two $g_{\mathrm{s}}$ models JAC and MED compared to observations. Seven FLUXNET towers were selected to represent a range of land cover types as shown in Table S2. JULES was set up for each site using observed site-level hourly meteorology, and the vegetation cover was prescribed according to the fractional covers of the different JULES surface types shown in Table S2. Following a spin-up period, simulations were run at each site for the years shown in Table S2.

Model simulations for Europe

Forcing datasets

We used the WATCH meteorological forcing dataset (Weedon et al., 2010, 2011) at 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ spatial and 3 h temporal resolution for our JULES simulations. JULES interpolates this down to an hourly model time step. For this study, the climate was kept constant by recycling over the period from 1901 to 1920 to allow us to focus on the impact of ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{CO}}_{\mathrm{2}}$ and their interactive effects.

JULES was run with prescribed annual mean atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentrations. Pre-industrial global ${\mathrm{CO}}_{\mathrm{2}}$ concentrations for 1900 to 1960 were taken from Etheridge et al. (1996), for 1960 to 2002 were from Mauna Loa (Keeling and Whorf, 2004), as calculated by the Global Carbon Project (Le Quéré et al., 2016), and for 2003–2050 were based on the IPCC SRES A1B scenario and were linearly interpolated to gap-fill missing years (Fig. 1).

Regional time series of canopy height ${\mathrm{O}}_{\mathrm{3}}$ (ppb) forcing from EMEP (a–c) and (d) global atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ (ppm) concentration (this does not vary regionally; black dots show data points; the black line shows interpolated points). Each panel for the ${\mathrm{O}}_{\mathrm{3}}$ forcing shows the regional annual average (woody PFTs, black solid line; herbaceous PFTs, black dashed line) and the annual maximum ${\mathrm{O}}_{\mathrm{3}}$ concentration above: woody PFTs (red) and herbaceous PFTs (blue).

[Figure omitted. See PDF]

JULES was run including dynamic vegetation with a land cover mask giving the fraction of agriculture in each 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ grid cell based on the Hurtt et al. (2011) land cover database for the year 2000. The agricultural mask is fixed and does not change over the simulation period. This means that whilst the model is allowed to evolve its own vegetation cover outside of the agricultural mask, within the mask only ${\mathrm{C}}_{\mathrm{3}}$ and ${\mathrm{C}}_{\mathrm{4}}$ herbaceous PFTs are allowed to grow, with no competition from other PFTs. Therefore, through the simulation period, regions of agriculture are maintained as such and not outcompeted by forests, for example, allowing for a more accurate representation of the land cover of Europe in the model. No form of land management is simulated (i.e. no crop harvesting, ploughing, rotation or grazing); growth and leaf area index (LAI) are determined by resource availability and phenology. Outside of the agricultural mask, dynamic vegetation means that grid cell PFT coverage and LAI are the result of resource availability, phenology and simulated competition. Across the model domain, simulated mean annual LAI was dominantly within the range of 2 to 5 m${}^{\mathrm{2}}$ m${}^{-\mathrm{2}}$ (Figs. S4 and S5). Following a full spin-up period (to ensure equilibrium vegetation, carbon and water states), there was no significant change in the fractional cover of each PFT over the simulation period (1901–2050). By 2050, increases in boreal forest cover occurred, but this was less than 2 % and limited to very small areas; given this small change we show just the land cover for 2050 in Fig. S2.

Tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentration was produced by the EMEP MSC-W model at 0.5${}^{\circ }$ $\times$ 0.5${}^{\circ }$ resolution (Simpson et al., 2012), driven with meteorology from the regional climate model RCA3 (Kjellström et al., 2011; Samuelsson et al., 2011), which provides a downscaling of ECHAM A1B-r3 (simulation 11 of Kjellström et al., 2011). This set-up (EMEP $+$ RCA3) is also used by Langner et al. (2012a), Simpson et al. (2014a), Tuovinen et al. (2013), Franz et al. (2017) and Engardt et al. (2017), in which further details and model evaluation can be found. Unfortunately, the three-dimensional RCA3 data needed by the EMEP model were not available prior to 1960, but as in Engardt et al. (2017) the meteorology of 1900–1959 had to be approximated by assigning random years from 1960 to 1969. This procedure introduces some uncertainty of course, although Langner et al. (2012b) show that for the period from 1990 to 2100 it is emissions change, rather than meteorological change, that drives modelled ${\mathrm{O}}_{\mathrm{3}}$ concentrations. The emissions scenarios for 1900–2050 merge data from the International Institute of Applied System Analysis (IIASA) for 2005–2050 (the so-called ECLIPSE 4a scenario), recently revised EMEP data for 1990 and a scaling back from 1990 to 1900 using data from Lamarque et al. (2013). The trend in emissions of the major ${\mathrm{O}}_{\mathrm{3}}$ precursors ${\mathrm{NO}}_x$, non-methane volatile organic compounds (NMVOCs) and isoprene are shown from 1900 to 2050 over Europe in Fig. S6. Isoprene emissions are not inputs to the EMEP model, but rather calculated at each time step using temperature, radiation and land-cover-specific emission factors (Simpson et al., 2012). Changes in the assumed background concentration of ${\mathrm{CH}}_{\mathrm{4}}$ (from RCP6.0) (van Vuuren et al., 2011) are also shown in Fig. S6. Engardt et al. (2017) show the trend in emissions of ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{NH}}_{\mathrm{3}}$ from 1900 to 2050 over Europe. The EMEP model accounts for changes in biogenic volatile organic compound (BVOC) emissions as a result of predicted ambient temperature changes.

${\mathrm{O}}_{\mathrm{3}}$ concentrations from EMEP MSC-W were calculated at canopy height for two land cover categories: forest and grassland (Figs. S7 and S8), which are taken as surrogates for high and low vegetation, respectively. These canopy-height-specific concentrations allow for the large gradients in ${\mathrm{O}}_{\mathrm{3}}$ concentration that can occur in the lowest tens of metres, giving lower ${\mathrm{O}}_{\mathrm{3}}$ for grasslands than seen at 20 m in a forest canopy, for example (Gerosa et al., 2017; Simpson et al., 2012; Tuovinen et al., 2009). These canopy-level ${\mathrm{O}}_{\mathrm{3}}$ concentrations are used as input to JULES, using the EMEP ${\mathrm{O}}_{\mathrm{3}}$ concentrations for forest for the forest JULES PFTs (broadleaf or needleleaf, tree and shrub), and the EMEP ${\mathrm{O}}_{\mathrm{3}}$ concentrations for grassland for the grass–herbaceous JULES PFTs (${\mathrm{C}}_{\mathrm{3}}$ and ${\mathrm{C}}_{\mathrm{4}}$). This study used daily mean values of tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentration from EMEP disaggregated down to the hourly JULES model time step. The daily mean ${\mathrm{O}}_{\mathrm{3}}$ forcing was disaggregated to follow a mean diurnal profile of ${\mathrm{O}}_{\mathrm{3}}$, which was generated from hourly ${\mathrm{O}}_{\mathrm{3}}$ output from EMEP MSC-W for the two land cover categories (forest and grassland as described above) across the same model domain. ${\mathrm{O}}_{\mathrm{3}}$ concentrations follow a diurnal cycle and peak during the day; therefore accounting for the diurnal variation in ${\mathrm{O}}_{\mathrm{3}}$ concentrations allows for a more realistic estimation of ${\mathrm{O}}_{\mathrm{3}}$ uptake.

Figure 1 shows large increases in tropospheric ${\mathrm{O}}_{\mathrm{3}}$ from pre-industrial to present day (2001), which is in line with modelling studies (Young et al., 2013) and site observations (Derwent et al., 2008; Logan et al., 2012; Parrish et al., 2012), and is predominantly a result of increasing anthropogenic emissions (Young et al., 2013). Figures S7 and S8 show this large increase in ground-level ${\mathrm{O}}_{\mathrm{3}}$ concentrations from 1901 to 2001 occurs in all seasons. Present-day ${\mathrm{O}}_{\mathrm{3}}$ concentration shows a strong seasonal cycle, with a spring–summer peak in concentrations in the mid-latitudes of the Northern Hemisphere (Derwent et al., 2008; Parrish et al., 2012; Vingarzan, 2004). Seasonal cycles have been changing over the past decades, however, attributed to changes in ${\mathrm{NO}}_x$ and other emissions, as well as changes in transport patterns (Parrish et al., 2013). These changes will likely continue in the future as emissions and meteorological factors impact photochemical ${\mathrm{O}}_{\mathrm{3}}$ production and transport patterns. Indeed, the ${\mathrm{O}}_{\mathrm{3}}$ concentrations used in the simulations in this study show increased ${\mathrm{O}}_{\mathrm{3}}$ levels in winter and in some regions in autumn and spring in 2050 compared to present day, which may be due to reduced titration of ${\mathrm{O}}_{\mathrm{3}}$ by NO as a result of reduced ${\mathrm{NO}}_x$ emissions in the future (Royal Society, 2008). Summer ${\mathrm{O}}_{\mathrm{3}}$ concentrations are lower in 2050, however, compared to 2001.

Spin-up and factorial experiments

JULES was spun up by recycling the climate from the early part of the twentieth century (1901 to 1920) using atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ (296.1 ppm) and ${\mathrm{O}}_{\mathrm{3}}$ concentrations from 1901 (Figs. S7 and S8). Model spin-up was 2000 years by which point the carbon pools and fluxes were in steady state with zero mean net land–atmosphere ${\mathrm{CO}}_{\mathrm{2}}$ flux. We performed the following transient simulations for the period from 1901 to 2050 with continued recycling of the climate as used in the spin-up, for both high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities:

• run_${\mathrm{O}}_{\mathrm{3}}$: fixed 1901 ${\mathrm{CO}}_{\mathrm{2}}$, varying ${\mathrm{O}}_{\mathrm{3}}$;

• run_${\mathrm{CO}}_{\mathrm{2}}$: varying ${\mathrm{CO}}_{\mathrm{2}}$, fixed 1901 ${\mathrm{O}}_{\mathrm{3}}$;

• run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$: varying ${\mathrm{CO}}_{\mathrm{2}}$, varying ${\mathrm{O}}_{\mathrm{3}}$.

We use these simulations to investigate the direct effects of changing atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ concentrations, individually and combined, on plant water use, GPP and the land C sink through the twentieth century and into the future, specifically over three time periods: historical (1901–2001), future (2001–2050) and over the full time series (1901–2050). For each time period we calculate the difference between the decadal means calculated at the start and end of the analysis period for each variable of interest. Therefore our results report the change in GPP, for example, over the analysis period. For each variable analysed (GPP, NPP, vegetation carbon, soil carbon, total land carbon and $g_{\mathrm{s}}$), we use the mean over 10 years to represent each time period, e.g. the mean over 2040 to 2050 is what we call 2050, and 1901 to 1910 is what we refer to as 1901. The difference among the simulations gives the effect of ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ either separately or in combination over the different time periods. We look at the percentage change due to either ${\mathrm{O}}_{\mathrm{3}}$ at pre-industrial ${\mathrm{CO}}_{\mathrm{2}}$ concentration (i.e. without the additional effect of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ on stomatal behaviour – run_${\mathrm{O}}_{\mathrm{3}}$), ${\mathrm{CO}}_{\mathrm{2}}$ (at fixed pre-industrial ${\mathrm{O}}_{\mathrm{3}}$ concentration, run_${\mathrm{CO}}_{\mathrm{2}}$) or the combined effect of both gases (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$), e.g. 100 $\cdot$ (var${\mathrm{O}}_{\mathrm{3}}$[2050]–var${\mathrm{O}}_{\mathrm{3}}$[1901])$/$var${\mathrm{O}}_{\mathrm{3}}$[1901] gives the ${\mathrm{O}}_{\mathrm{3}}$ effect (at fixed ${\mathrm{CO}}_{\mathrm{2}}$) over the full experimental period. The meteorological forcing is prescribed in these simulations and is therefore the same among the model runs. Other climate factors, such as VPD, temperature and soil moisture availability are accounted for in our simulations, but our analysis isolates the effects of ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ $+$ ${\mathrm{CO}}_{\mathrm{2}}$. We also use a paired $t$ test to determine statistically significant differences between the different (high and low) plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities.

Evaluation

To evaluate our JULES simulations we compare mean GPP from 1991 to 2001 for each of the JULES scenarios and both high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities against the observation-based globally extrapolated Flux Network model tree ensemble (MTE) (Jung et al., 2011). We use a paired $t$ test to determine statistically significant differences in the mean responses.

Results

Impact of $g_{\mathrm{s}}$ model formulation and site-level evaluation

The impact of the $g_{\mathrm{s}}$ model on simulated $g_{\mathrm{s}}$ is shown for the site with low soil moisture stress (wet site, Fig. 2). For the broadleaf tree and ${\mathrm{C}}_{\mathrm{3}}$ herbaceous PFT, the MED model simulates a larger conductance compared to the JAC model. In other words, with the MED model these two PFTs are parameterised with a less conservative water use strategy, which, for the grid point shown in Fig. 2, increased the annual mean water use by 35 % ($\pm$29 %) and 45 % ($\pm$32 %), respectively. In contrast, the needleleaf tree, ${\mathrm{C}}_{\mathrm{4}}$ herbaceous and shrub PFTs are parameterised with a more conservative water use strategy with the MED model, and the mean annual $g_{\mathrm{s}}$ was decreased by 13 % ($\pm$12 %), 27 % ($\pm$10 %) and 36 % ($\pm$13 %), respectively, compared to the JAC model. This comparison was also performed for a dry site (high soil moisture stress), and similar results were found (Fig. S9). The effect of $g_{\mathrm{s}}$ formulation on simulated photosynthesis was much smaller because of the lower sensitivity of the limiting rates of photosynthesis to changes in $c_i$ in the model compared to the effect of the same change in $c_i$ on modelled $g_{\mathrm{s}}$ (Figs. S10 and S11). Changes in $g_{\mathrm{s}}$ impact the partitioning of simulated energy fluxes. In general, increased $g_{\mathrm{s}}$ results in increased latent heat and thus decreased sensible heat flux, and vice versa where $g_{\mathrm{s}}$ is decreased (Figs. S10 and S11). Also shown is the effect of the MED model on ${\mathrm{O}}_{\mathrm{3}}$ flux into the leaf (Figs. S12 and S9 bottom panel). For the broadleaf tree and ${\mathrm{C}}_{\mathrm{3}}$ herbaceous PFTs, the MED model simulates a larger conductance and therefore a greater flux of ${\mathrm{O}}_{\mathrm{3}}$ through stomata compared to JAC, and this is indicative of the potential for greater reductions in photosynthesis (Figs. S10 and S11 top row). The reverse is seen for the needleleaf tree, ${\mathrm{C}}_{\mathrm{4}}$ herbaceous and shrub PFTs.

Comparison of simulated $g_{\mathrm{s}}$ with MED ($y$ axis) versus JAC ($x$ axis) for all five JULES PFTs at one grid point (lat: 48.25, long: 5.25); shown are hourly values for the year 2000.

[Figure omitted. See PDF]

Site-level evaluation of the seasonal cycles of latent and sensible heat with both JAC and MED models compared to FLUXNET observations showed in general the MED model improved the seasonal cycle of both fluxes (lower RMSE), but the magnitude of this varied from site to site (Fig. S13). At the deciduous broadleaf site, US-UMB, MED resulted in improvements of the simulated seasonal cycle particularly in the summer months for both fluxes (RMSE decreased from 42.7 and 31.5 to 38.5 and 28.0 W m${}^{-\mathrm{2}}$ for latent and sensible heat, respectively). At the second deciduous broadleaf site IT-CA1, however, there was almost no difference between the two $g_{\mathrm{s}}$ models. Both evergreen needleleaf forest sites (FI-Hyy and DE-Tha) saw improvements in the simulated seasonal cycles of latent and sensible heat with the MED model, primarily as a result of lower latent heat flux in the spring and summer months and higher sensible heat flux over the same period. At FI-Hyy, RMSE decreased from 10.1 and 7.4 to 6.7 and 6.7 W m${}^{-\mathrm{2}}$ for latent and sensible heat, respectively, and at DE-Tha, RMSE decreased from 16.0/11.9 to 10.5/10.6 W m${}^{-\mathrm{2}}$ for latent/sensible heat, respectively. With the MED model the monthly mean latent heat flux was improved at the ${\mathrm{C}}_{\mathrm{3}}$ grass site (CH-Cha) as a result of increased flux in the summer months (RMSE decreased from 15.7 to 13.8 W m${}^{-\mathrm{2}})$; however there was no improvement in the sensible heat flux and RMSE with MED was increased (from 3.9 to 4.9 W m${}^{-\mathrm{2}})$. At the ${\mathrm{C}}_{\mathrm{4}}$ grass site (US-SRG), small improvements were made in the seasonal cycle of both latent and sensible heat with the MED model. At the deciduous savannah site (CG-Tch), which included a high proportion of shrub PFT in the land cover type used in the site simulation, large improvements in the seasonal cycle of both fluxes were simulated with the MED model, as a result of a decrease in the latent heat flux and an increase in the sensible heat flux (RMSE decreased from 39.5 and 31.6 to 30.4 and 24.4 W m${}^{-\mathrm{2}}$ for latent and sensible heat, respectively).

Evaluation of the JULES ${\mathrm{O}}_{\mathrm{3}}$ model

For all JULES scenarios, similar spatial patterns of GPP are simulated compared to MTE (Figs. 3 and S14). MTE estimates a mean GPP for present day in Europe of 938 gC m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (Fig. 3). JULES tends to under-predict GPP relative to the MTE product. Estimates of GPP from JULES with both transient ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$) give a mean across Europe of 813 gC m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity) to 881 gC m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity), both of which are significantly different from the MTE product ($t=\mathrm{27}$, $d.f.=\mathrm{5750}$, $p$ < 2.2e${}^{-\mathrm{16}}$ (high); $t=\mathrm{4.3}$, $d.f.=\mathrm{5750}$, $p$ < 1.5e${}^{-\mathrm{05}}$ (low); Fig. 3). Forcing with ${\mathrm{CO}}_{\mathrm{2}}$ alone (run_${\mathrm{CO}}_{\mathrm{2}}$) gives a mean GPP across Europe of 900 to 923 gC m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, respectively), and ${\mathrm{O}}_{\mathrm{3}}$ alone (run_${\mathrm{O}}_{\mathrm{3}}$ – without the protective effect of ${\mathrm{CO}}_{\mathrm{2}}$) reduces estimated GPP to 732 to 799 gC m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (Fig. S14). At latitudes > 45${}^{\circ }$ N JULES has a tendency to under-predict MTE GPP, and at latitudes < 45${}^{\circ }$ N JULES tends to over-predict MTE GPP (Fig. S15). These regional differences are highlighted in Fig. S16, in which in the Mediterranean region, JULES tends to over-predict compared to MTE GPP, so simulations with ${\mathrm{O}}_{\mathrm{3}}$ reduce the simulated GPP, bringing it closer to MTE. In the temperate region however, JULES tends to underestimate MTE GPP, so the addition of ${\mathrm{O}}_{\mathrm{3}}$ reduces simulated GPP further (Fig. S16). In the boreal region, JULES under-predicts GPP, but to a lesser extent than in the temperate region, and the addition of ${\mathrm{O}}_{\mathrm{3}}$ has less impact on reducing the GPP further (Fig. S16).

Mean GPP (g C m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}})$ from 1991 to 2001 for (a) the observationally based globally extrapolated Flux Network model tree ensemble (MTE) (Jung et al., 2011); (b, c) model simulations with transient ${\mathrm{CO}}_{\mathrm{2}}$ and transient ${\mathrm{O}}_{\mathrm{3}}$ (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$), with high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively.

[Figure omitted. See PDF]

Simulated percentage change in total carbon stocks (Land C) and gross primary productivity (GPP) due to ${\mathrm{O}}_{\mathrm{3}}$ effects at fixed pre-industrial atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (run_${\mathrm{O}}_{\mathrm{3}}$). Changes are shown for the periods 1901 to 2001 and 2001 to 2050 for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities.

[Figure omitted. See PDF]

European simulations – historical period: 1901–2001

Over the historical period (1901–2001), run_${\mathrm{O}}_{\mathrm{3}}$ reduced GPP under both the low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity parameterisations by $-\mathrm{3}$ to $-\mathrm{9}$ %, respectively (Table 1), and this difference in simulated GPP was significant ($t=\mathrm{102.2}$, $d.f.=\mathrm{6270}$, $p$ < 2.2e${}^{-\mathrm{16}})$. Figure 4 highlights regional variations; simulated reductions in GPP are up to 20 % across large areas of Europe, and up to 30 % in some Mediterranean regions under the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity. Some boreal and Mediterranean regions show small increases in GPP over this period, associated with ${\mathrm{O}}_{\mathrm{3}}$-induced stomatal closure enhancing water availability in these drier regions (Fig. 5). This allows for greater stomatal conductance later in the year when soil moisture may otherwise have been limiting to growth (up to 10 %, Fig. 5), and therefore caused higher GPP, but these regions comprise only a small area of the entire domain. Indeed, over much of Europe, ${\mathrm{O}}_{\mathrm{3}}$-induced stomatal closure led to reduced $g_{\mathrm{s}}$ (up to 20 %) across large areas of temperate Europe and the Mediterranean, and even greater reductions in some smaller regions of the southern Mediterranean (Fig. 6), and these are not associated with notable increases in soil moisture availability (Fig. 5), resulting in depressed GPP over much of Europe as described above. Under the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, similar spatial patterns occur, but the magnitude of GPP change (up to $-\mathrm{10}$ % across much of Europe) and $g_{\mathrm{s}}$ change ($-\mathrm{5}$ to $-\mathrm{10}$ %) is lower compared to the high sensitivity. Over the twentieth century the land carbon sink is suppressed ($-\mathrm{2}$ to $-\mathrm{6}$ %, Table 1). Large regional variation is shown in Fig. 4, with temperate and Mediterranean Europe seeing a large reduction in land carbon storage, particularly under the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity (up to $-\mathrm{15}$ %).

Simulated changes in the European land carbon cycle due to changing ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ concentrations (independently and together). Shown are changes in total carbon stocks (Land C) and gross primary productivity (GPP) over three different periods (historical: 1901 to 2001; future: 2001 to 2050; and full time series: 1901 to 2050). Absolute (top) and relative (bottom) differences are shown. For 2001 to 2050, please refer to Table S4 for the initial value for each run. See the Supplement for details of the estimation of the ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ effects and their interaction.

	High plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity
	1901–2001		2001–2050		1901–2050
	GPP	Land C	GPP	Land C	GPP	Land C
	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)
Value in 1901	9.05	167	–	–	9.05	167
Absolute change
${\mathrm{O}}_{\mathrm{3}}$	$-\mathrm{0.8}$1	$-\mathrm{9.21}$	0.01	$-\mathrm{2.44}$	$-\mathrm{0.80}$	$-\mathrm{11.65}$
${\mathrm{CO}}_{\mathrm{2}}$	1.16	4.24	1.42	12.98	2.58	17.22
${\mathrm{CO}}_{\mathrm{2}}$ $+$ ${\mathrm{O}}_{\mathrm{3}}$	0.13	$-\mathrm{3.28}$	1.66	11.11	1.79	7.83
% change
${\mathrm{O}}_{\mathrm{3}}$	$-\mathrm{8.95}$	$-\mathrm{5.51}$	0.12	$-\mathrm{1.55}$	$-\mathrm{8.84}$	$-\mathrm{6.98}$
${\mathrm{CO}}_{\mathrm{2}}$	12.82	2.54	13.91	7.58	28.51	10.31
${\mathrm{CO}}_{\mathrm{2}}$ $+$ ${\mathrm{O}}_{\mathrm{3}}$	1.44	$-\mathrm{1.96}$	18.08	6.79	19.78	4.69
	Low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity
	1901–2001		2001–2050		1901–2050
	GPP	Land C	GPP	Land C	GPP	Land C
	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)	(Pg C yr${}^{-\mathrm{1}})$	(Pg C)
Value in 1901	9.34	167.5	–	–	9.34	167.5
Absolute change
${\mathrm{O}}_{\mathrm{3}}$	$-\mathrm{0.30}$	$-\mathrm{3.59}$	0.02	$-\mathrm{1.07}$	$-\mathrm{0.40}$	$-\mathrm{4.66}$
${\mathrm{CO}}_{\mathrm{2}}$	1.15	6.43	1.35	13.14	2.50	19.57
${\mathrm{CO}}_{\mathrm{2}}$ $+$ ${\mathrm{O}}_{\mathrm{3}}$	0.65	2.50	1.50	12.35	2.15	14.85
% change
${\mathrm{O}}_{\mathrm{3}}$	$-\mathrm{3.21}$	$-\mathrm{2.14}$	0.22	$-\mathrm{0.65}$	$-\mathrm{4.28}$	$-\mathrm{2.78}$
${\mathrm{CO}}_{\mathrm{2}}$	12.31	3.84	12.87	7.55	26.77	11.68
${\mathrm{CO}}_{\mathrm{2}}$ $+$ ${\mathrm{O}}_{\mathrm{3}}$	6.96	1.49	15.02	7.26	23.02	8.87

Combined, the physiological response to changing ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ concentrations (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$) results in a net loss of land carbon over the twentieth century under the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity ($-\mathrm{2}$ %, Table 1), dominated by loss of soil carbon (Table S3). This reflects the large increases in tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentration observed over this period (Fig. 1). Under low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, the land carbon sink has started to recover by 2001 ($+\mathrm{1.5}$ %) owing to the recovery of the soil carbon pool beyond 1901 values over this period (Table S3).

Simulated percentage change in plant available soil moisture (fsmc) due to ${\mathrm{O}}_{\mathrm{3}}$ effects at fixed pre-industrial atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (run_${\mathrm{O}}_{\mathrm{3}}$). Changes are shown for the periods 1901 to 2001 and 2001 to 2050 for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities.

[Figure omitted. See PDF]

To gain perspective on the magnitude of the ${\mathrm{O}}_{\mathrm{3}}$-induced flux of carbon from the land to the atmosphere we relate changes in total land carbon to carbon emissions from fossil fuel combustion and cement production for the EU28-plus countries (EU-28 plus Albania, Bosnia and Herzegovina, Iceland, Belarus, Serbia, Moldova, Norway, Turkey, Ukraine, Switzerland and Macedonia) from the data of Boden et al. (2013). We recognise that our simulation domain is slightly larger than the EU28-plus as it includes a small area of western Russia so direct comparisons cannot be made, but this still provides a useful measure of the size of the carbon flux. For the period 1970 to 1979 the simulated loss of carbon from the European terrestrial biosphere due to ${\mathrm{O}}_{\mathrm{3}}$ effects on vegetation physiology was on average 1.32 Pg C (high vegetation sensitivity) and 0.71 Pg C (low vegetation sensitivity) (Table 2). This ${\mathrm{O}}_{\mathrm{3}}$-induced reduced C uptake of the land surface is equivalent to around 8 to 16 % of the emissions of carbon from fossil fuel combustion and cement production over the same period for the EU28-plus countries (Table 2). Currently the emissions data availability goes up to 2011. Over the last observable decade (2002 to 2011) the simulated reduction in land carbon due to ${\mathrm{O}}_{\mathrm{3}}$ has declined, but is still equivalent to 2 to 4 % of the emissions of carbon from fossil fuels and cement production for the EU28-plus countries (Table 2). By comparison with one of the largest anthropogenic emissions of carbon for Europe, we show here that the potential effect of ${\mathrm{O}}_{\mathrm{3}}$ on reducing the size of the European land carbon sink is notable.

Simulated change in total land carbon due to ${\mathrm{O}}_{\mathrm{3}}$ damage with changing atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration for the two vegetation sensitivities. The sum of carbon emissions for each decade from fossil fuel combustion and cement production for the EU28 plus Albania, Bosnia and Herzegovina, Iceland, Belarus, Serbia, Moldova, Norway, Turkey, Ukraine, Switzerland and Macedonia (EU28-plus) is shown; the data are from Boden et al. (2013). The simulated change in land carbon as a result of ${\mathrm{O}}_{\mathrm{3}}$ damage is depicted as a percentage of the EU28-plus emissions to demonstrate the magnitude of the additional source of carbon to the atmosphere from plant ${\mathrm{O}}_{\mathrm{3}}$ damage.

	Mean (Pg C)
	1970–1979	1980–1989	1990–1999	2000–2009	2002–2011
Modelled ${\mathrm{O}}_{\mathrm{3}}$ effect on land C sink
Higher sensitivity	$-\mathrm{1.32}$	$-\mathrm{1.01}$	$-\mathrm{0.97}$	$-\mathrm{0.53}$	$-\mathrm{0.50}$
Low sensitivity	$-\mathrm{0.71}$	$-\mathrm{0.58}$	$-\mathrm{0.50}$	$-\mathrm{0.29}$	$-\mathrm{0.26}$
Sum of C emissions from fossil fuel	8.39	8.63	12.26	12.83	12.75
combustion and cement production (Pg C)
C lost from ${\mathrm{O}}_{\mathrm{3}}$ effect as a percentage of fossil fuel
and cement emissions (%)
Higher sensitivity	$-\mathrm{15.73}$	$-\mathrm{11.70}$	$-\mathrm{7.91}$	$-\mathrm{4.13}$	$-\mathrm{3.92}$
Low sensitivity	$-\mathrm{8.46}$	$-\mathrm{6.72}$	$-\mathrm{4.08}$	$-\mathrm{2.26}$	$-\mathrm{2.04}$

European simulations – future period: 2001–2050

Over the 2001 to 2050 period, region-wide GPP with ${\mathrm{O}}_{\mathrm{3}}$ only changing (run_${\mathrm{O}}_{\mathrm{3}}$) increased marginally ($+\mathrm{0.1}$ to $+\mathrm{0.2}$ %, high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities; Table 1, with a significant difference between the two plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities; $t=\mathrm{57}$, $d.f.=\mathrm{6270}$ $p$ < 2.2e${}^{-\mathrm{16}})$, although with large spatial variability as discussed below (Fig. 4g and h). Figures S7 and S8 show that despite decreased tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentrations by 2050 in summer compared to 2001 levels, all regions are exposed to an increase in ${\mathrm{O}}_{\mathrm{3}}$ over the wintertime, and some regions of Europe, particularly temperate and Mediterranean regions experience increases in ${\mathrm{O}}_{\mathrm{3}}$ concentration in spring and autumn. Therefore, although in the ${\mathrm{O}}_{\mathrm{3}}$ simulation overall simulated GPP for Europe shows a small increase, large spatial variability is shown in Fig. 4g and h because of the variability in ${\mathrm{O}}_{\mathrm{3}}$ concentration with region and season. Increased GPP (dominantly 10, but up to 20 % in some areas) on 2001 levels is simulated across areas of Europe; however, decreases of up to 21 % are simulated in some areas of the Mediterranean, up to 15 % in some areas of the boreal region and up to 27 % in the temperate zone (Fig. 4g and h).

When ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ effects are combined (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$), simulated GPP increases ($+\mathrm{15}$ to $+\mathrm{18}$ %, high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively, Table 1). This increase is greater than the enhancement simulated when ${\mathrm{CO}}_{\mathrm{2}}$ affects plant growth independently (run_${\mathrm{CO}}_{\mathrm{2}}$) because additional ${\mathrm{O}}_{\mathrm{3}}$-induced stomatal closure increases soil water availability in some regions, which enhances growth more in run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$, compared to run_${\mathrm{CO}}_{\mathrm{2}}$. Nevertheless, although the percentage gain is larger, the absolute value of GPP by 2050 remains lower in run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$ compared to GPP in run_${\mathrm{CO}}_{\mathrm{2}}$, highlighting the negative impact of ${\mathrm{O}}_{\mathrm{3}}$ at the land surface (Table S4).

Despite small increases in GPP in run_${\mathrm{O}}_{\mathrm{3}}$, the land carbon sink continues to decline from 2001 levels ($-\mathrm{0.7}$ to $-\mathrm{1.6}$ %, low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively; Table 1). This is because the soil and vegetation carbon pools continue to lose carbon as they adjust slowly to small changes in input (GPP); i.e. the soil carbon pool is not in equilibrium in 2001, and is declining in response to reduced litter input as a result of twentieth century ${\mathrm{O}}_{\mathrm{3}}$ impacts on GPP. Nevertheless, the negative effect of ${\mathrm{O}}_{\mathrm{3}}$ on the future land sink is markedly reduced relative to the historical period. Figure 4e and f, however, highlight regional differences. Boreal regions and parts of central Europe see minimal ${\mathrm{O}}_{\mathrm{3}}$ damage, whereas some areas of southern and northern Europe see further losses of up to 8 % on 2001 levels. The run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$ simulation is dominated by the physiological effects of changing ${\mathrm{CO}}_{\mathrm{2}}$, with land carbon sink increases of up to 7 % (Table 1).

Simulated percentage change in stomatal conductance ($g_{\mathrm{s}})$ due to ${\mathrm{O}}_{\mathrm{3}}$ effects at the fixed pre-industrial atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (run_${\mathrm{O}}_{\mathrm{3}}$). Changes are shown for the periods 1901 to 2001 and 2001 to 2050 for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities.

[Figure omitted. See PDF]

European simulations – full experimental period: 1901–2050

From 1901 to 2050, run_${\mathrm{O}}_{\mathrm{3}}$ reduces GPP ($-\mathrm{4}$ to $-\mathrm{9}$ %, with a significant difference between the low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities ($t=\mathrm{95}$, $d.f.=\mathrm{6270}$ $p$ < 2.2e${}^{-\mathrm{16}}\left)\right)$ and land carbon storage ($-\mathrm{3}$ to $-\mathrm{7}$ %, Table 1). Regionally, ${\mathrm{O}}_{\mathrm{3}}$ damage is lowest in the boreal zone. GPP decreases are largely between 5–8% and 2–4 % for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively, with large areas minimally affected by ${\mathrm{O}}_{\mathrm{3}}$ damage (Fig. 7), consistent with lower $g_{\mathrm{s}}$ values of needleleaf trees that dominate this region, and so lower ${\mathrm{O}}_{\mathrm{3}}$ uptake (Figs. S17 and S18). In the temperate region, ${\mathrm{O}}_{\mathrm{3}}$ damage is extensive, with reductions in GPP dominantly from 10 to 15 % for the low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively. Across significant areas of this region reductions in GPP are up to 20 % under high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity (Fig. 7). In the Mediterranean region, ${\mathrm{O}}_{\mathrm{3}}$ damage reduces GPP by 5 to 15 % and 3 to 6 % for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively, with some areas seeing greater losses of up to 20 % under the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, but this is less extensive than that seen in the temperate zone (Fig. 7). In these drier regions, ${\mathrm{O}}_{\mathrm{3}}$-induced stomatal closure can increase available soil moisture (Fig. S17 and S18).

The run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$ simulation shows that ${\mathrm{CO}}_{\mathrm{2}}$-induced stomatal closure can help alleviate ${\mathrm{O}}_{\mathrm{3}}$ damage by reducing the uptake of ${\mathrm{O}}_{\mathrm{3}}$ (Table S6). In these simulations, ${\mathrm{CO}}_{\mathrm{2}}$-induced stomatal closure was found to offset ${\mathrm{O}}_{\mathrm{3}}$ suppression of GPP, such that GPP by 2050 is 3 to 7 % lower due to ${\mathrm{O}}_{\mathrm{3}}$ exposure (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$), rather than 4 to 9 % lower in the absence of increasing ${\mathrm{CO}}_{\mathrm{2}}$ (run_${\mathrm{O}}_{\mathrm{3}}$, Table S6). Figure 6 shows this spatially, ${\mathrm{O}}_{\mathrm{3}}$ damage is reduced when the effect of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ on stomatal closure is accounted for; however, despite this, the land carbon sink and GPP remain significantly reduced due to ${\mathrm{O}}_{\mathrm{3}}$ exposure.

From 1901 to 2050, run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$ results in an increase in European land carbon uptake ($+\mathrm{5}$ to $+\mathrm{9}$ %) and an increase in GPP ($+\mathrm{20}$ to $+\mathrm{23}$ %) by 2050 for the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively (Table 1). Nevertheless, despite this increase there remains a large negative impact of ${\mathrm{O}}_{\mathrm{3}}$ on the European land carbon sink (Fig. S19). By 2050 the simulated enhancement of land carbon and GPP in response to elevated ${\mathrm{CO}}_{\mathrm{2}}$ alone (run_${\mathrm{CO}}_{\mathrm{2}}$) is reduced by 3 to 6 % (land carbon) and 4 to 9 % (GPP) for the low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively, when ${\mathrm{O}}_{\mathrm{3}}$ is also accounted for (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$, Table 1). This is a large reduction in the ability of the European terrestrial biosphere to sequester carbon.

Discussion

Evaluation of $g_{\mathrm{s}}$ models and JULES ${\mathrm{O}}_{\mathrm{3}}$ model

Comparison of the new $g_{\mathrm{s}}$ model implemented in this study (MED) with the $g_{\mathrm{s}}$ model currently used as a standard in JULES (JAC) revealed large differences in $g_{\mathrm{s}}$ for each PFT, principally as a result of the data-based parameterisation of the new model. Water use increased for the broadleaf tree and ${\mathrm{C}}_{\mathrm{3}}$ herbaceous PFTs using the MED model compared to JAC, but decreased for the needleleaf tree, ${\mathrm{C}}_{\mathrm{4}}$ herbaceous and shrub PFTs, which displayed a more conservative water use strategy compared to JAC. These changes are in line with the work of De Kauwe et al. (2015), who found a reduction in annual transpiration for evergreen needleleaf, tundra and ${\mathrm{C}}_{\mathrm{4}}$ grass regions when implementing the Medlyn $g_{\mathrm{s}}$ model into the Australian land surface scheme CABLE. Site-level evaluation of the models against FLUXNET observations showed that in general the MED model improved simulated seasonal cycles of latent and sensible heat. The magnitude of the improvement varied with site; improvements were seen at the deciduous savanna site, and at the needleleaf tree sites and broadleaf tree site (US_UMB) in the spring and summer. However, much smaller improvements were seen at the grass sites. Changes in $g_{\mathrm{s}}$ in this study resulted in differences in latent and sensible heat fluxes. Changes in the partitioning of energy fluxes at the land surface could have consequences for the intensity of heatwaves (Cruz et al., 2010; Kala et al., 2016), run-off (Betts et al., 2007; Gedney et al., 2006) and rainfall patterns (de Arellano et al., 2012), although fully coupled simulations would be necessary to detect these effects. The differences in simulated $g_{\mathrm{s}}$ led to differences in uptake of ${\mathrm{O}}_{\mathrm{3}}$ between the two models because the rate of $g_{\mathrm{s}}$ is the predominant determinant of the flux of ${\mathrm{O}}_{\mathrm{3}}$ through stomata. Higher ${\mathrm{O}}_{\mathrm{3}}$ uptake is indicative of greater damage. Therefore, given that ${\mathrm{C}}_{\mathrm{3}}$ herbaceous vegetation is the dominant land cover class across the European domain used in this study, this suggests a greater ${\mathrm{O}}_{\mathrm{3}}$ impact for Europe would be simulated with the MED model compared to JAC in our simulations in which chemistry is uncoupled from the land surface.

Simulated percentage change in total carbon stocks (Land C) and gross primary productivity (GPP) due to (i) (a, e, i, m) ${\mathrm{O}}_{\mathrm{3}}$ effects at fixed pre-industrial atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (run_${\mathrm{O}}_{\mathrm{3}}$), (ii) (b, f, j, n) ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation at fixed pre-industrial ${\mathrm{O}}_{\mathrm{3}}$ concentration (run_${\mathrm{CO}}_{\mathrm{2}}$), (iii) (c, g, k, o) the interaction between ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ effects (run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$), and (iv) (d, h, l, p) ${\mathrm{O}}_{\mathrm{3}}$ effects with changing atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (i.e. ${\mathrm{O}}_{\mathrm{3}}$ damage accounting for the effect of ${\mathrm{CO}}_{\mathrm{2}}$-induced stomatal closure; run_both_${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$ – run_${\mathrm{CO}}_{\mathrm{2}}$). Changes are depicted for the period from 1901 to 2050 for high and low ozone plant sensitivities.

[Figure omitted. See PDF]

We evaluated the JULES ${\mathrm{O}}_{\mathrm{3}}$ model by comparing modelled GPP against the Jung et al. (2011) MTE product. Similar spatial patterns of GPP were simulated by JULES compared to MTE. Zonal means also showed similar patterns of GPP, although JULES under-predicted GPP compared to MTE at latitudes > 45${}^{\circ }$ N (temperate and boreal regions; all simulations) and over-predicted GPP at latitudes < 45${}^{\circ }$ N (Mediterranean region; all simulations). The simulations with transient ${\mathrm{O}}_{\mathrm{3}}$ (i.e. ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$$+$${\mathrm{O}}_{\mathrm{3}}$) showed large differences in GPP between the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity simulations, which is to be expected given that the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity simulations were parameterised to be “damaged” more by ${\mathrm{O}}_{\mathrm{3}}$, i.e. greater reduction of photosynthesis and $g_{\mathrm{s}}$ with ${\mathrm{O}}_{\mathrm{3}}$ exposure compared to the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity simulations. This difference was largest in the temperate zone, largely because of ${\mathrm{C}}_{\mathrm{3}}$ grass cover being the dominant land cover here and the difference in the sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ between the high and low calibrations is significantly larger for ${\mathrm{C}}_{\mathrm{3}}$ grasses compared to the needleleaf trees that dominate in the boreal region. Additionally, a longer growing season in the temperate region may allow for greater uptake of ${\mathrm{O}}_{\mathrm{3}}$ into vegetation. ${\mathrm{C}}_{\mathrm{3}}$ grass is also the dominant land cover in the Mediterranean region with a different calibration used for Mediterranean grasses for the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, which is less sensitive to ${\mathrm{O}}_{\mathrm{3}}$ than the temperate ${\mathrm{C}}_{\mathrm{3}}$ grasses, but high soil moisture stress is common throughout the growing season in the Mediterranean, limiting the uptake of ${\mathrm{O}}_{\mathrm{3}}$ through stomata, which likely diminishes the difference between the high and low calibrations. In general, incorporating plant ${\mathrm{O}}_{\mathrm{3}}$ damage into JULES leads to worse agreement with the MTE GPP product; however, this is expected to some degree as we are adding an explicit representation of ${\mathrm{O}}_{\mathrm{3}}$ damage to a model calibrated to reproduce present-day GPP and drawdown of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. Inevitably this implicitly includes ${\mathrm{O}}_{\mathrm{3}}$ damage to vegetation. Explicit representation of plant ${\mathrm{O}}_{\mathrm{3}}$ damage is important to investigate how ${\mathrm{O}}_{\mathrm{3}}$ damage changes through time under different emissions scenarios and the interactive effects with other gases (such as ${\mathrm{CO}}_{\mathrm{2}}$) and with climate change. The percentage changes we simulate are therefore important to demonstrate the sensitivity of modelled GPP and land carbon to this process.

Comparison of modelled estimates of ${\mathrm{O}}_{\mathrm{3}}$ damage

Our estimates suggest ${\mathrm{O}}_{\mathrm{3}}$ (run_${\mathrm{O}}_{\mathrm{3}}$) reduced GPP by 2001 by 3 to 9 % on average across Europe and NPP by 5 to 11 % for the low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively (Table S3). Anav et al. (2011) simulated a 22 % reduction of GPP across Europe for 2002 using the ORCHIDEE model. Present-day ${\mathrm{O}}_{\mathrm{3}}$ exposure reduced GPP by 10 to 25 % in Europe and 10.8 % globally in the study by Lombardozzi et al. (2015) using the Community Land Model (CLM). ${\mathrm{O}}_{\mathrm{3}}$ reduced NPP by 11.2 % in Europe from 1989 to 1995 using the Terrestrial Ecosystem Model (TEM) (Felzer et al., 2005). Globally, concentrations of ${\mathrm{O}}_{\mathrm{3}}$ predicted for 2100 reduced GPP by 14 to 23 % using a former parameterisation of ${\mathrm{O}}_{\mathrm{3}}$ sensitivity in JULES (Sitch et al., 2007). The recent study by Franz et al. (2017) showed mean GPP declined by 4.7 % over the period 2001 to 2010 using the OCN model over the same European domain and using the same ${\mathrm{O}}_{\mathrm{3}}$ forcing produced by EMEP MSC-W as used in this study. Our estimates of changes in present-day GPP and NPP are at the lower end of previously modelled estimates. Simulated ${\mathrm{O}}_{\mathrm{3}}$ impacts will be dependent on model ${\mathrm{O}}_{\mathrm{3}}$ concentrations, meteorology, plant sensitivity to ${\mathrm{O}}_{\mathrm{3}}$ and process representation of ${\mathrm{O}}_{\mathrm{3}}$ damage. It is therefore difficult to hypothesise exactly why modelled estimates differ, but it is suggested that an ensemble approach to modelling ${\mathrm{O}}_{\mathrm{3}}$ impacts on the terrestrial biosphere would be beneficial to understand some of these differences and provide estimates of ${\mathrm{O}}_{\mathrm{3}}$ damage with uncertainties.

Impacts of ${\mathrm{O}}_{\mathrm{3}}$ at the land surface

In this study, ${\mathrm{O}}_{\mathrm{3}}$ has a detrimental effect on the size of the land carbon sink for Europe. This is primarily through a decrease in the size of the soil carbon pool as a result of reduced litter input to the soil, consistent with reduced GPP and NPP. Field studies show that in some regions of Europe, soil carbon stocks are decreasing (Bellamy et al., 2005; Capriel, 2013; Heikkinen et al., 2013; Sleutel et al., 2003). The study of Bellamy et al. (2005), for example, showed that carbon was lost from soils across England and Wales between 1978 and 2003 at a mean rate of 0.6 % per year with little effect of land use on the rate of carbon loss, suggesting a possible link to climate change. It is understood that climate change is likely to affect soil carbon turnover. Increased temperatures increase microbial decomposition activity in the soil and therefore increase carbon losses through higher rates of respiration (Cox et al., 2000; Friedlingstein et al., 2006; Jones et al., 2003). However, some studies have found that ${\mathrm{O}}_{\mathrm{3}}$ can decrease soil carbon content. Talhelm et al. (2014), for example, found ${\mathrm{O}}_{\mathrm{3}}$ reduced carbon content in near surface mineral soil of forest soils exposed to 11 years of ${\mathrm{O}}_{\mathrm{3}}$ fumigation. Hofmockel et al. (2011) found elevated ${\mathrm{O}}_{\mathrm{3}}$ reduced the carbon content in more stable soil organic matter pools, and Loya et al. (2003) showed that the fraction of soil carbon formed in forest soils over a 4-year experimental period when fumigated with both ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ was reduced by 51 % compared to the soil fumigated with ${\mathrm{CO}}_{\mathrm{2}}$ alone. It is agreed that amongst other factors that change with ${\mathrm{O}}_{\mathrm{3}}$ exposure such as litter quality and composition, reduced litter quantity also has significant detrimental consequences for soil carbon stocks (Andersen, 2003; Lindroth, 2010; Loya et al., 2003). Results from this study therefore suggest that increasing tropospheric ${\mathrm{O}}_{\mathrm{3}}$ may be a contributing factor to the declining soil carbon stocks observed across Europe as a result of reduced litter input to the soil carbon pool consistent with reduced NPP.

We acknowledge, however, that our model simulations do not include coupling of nitrogen and carbon cycles or land management practices. We include a representation of agricultural regions through the model calibration against the wheat ${\mathrm{O}}_{\mathrm{3}}$ sensitivity function (CLRTAP, 2017), and in our simulations the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity scenario uses this calibration against wheat for all ${\mathrm{C}}_{\mathrm{3}}$$/$${\mathrm{C}}_{\mathrm{4}}$ land cover, which dominates our model domain. Wheat is known to be one of the most ${\mathrm{O}}_{\mathrm{3}}$-sensitive crop species, however, so it is possible that our simulations overestimate the ${\mathrm{O}}_{\mathrm{3}}$ impact at the land surface. However, the low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity calibration against natural grasslands provides a counter estimate of the impact of ${\mathrm{O}}_{\mathrm{3}}$ at the land surface; therefore it is important to consider the range our results provide (i.e. both the high and low plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities) as an indicator of the impact of ${\mathrm{O}}_{\mathrm{3}}$ on the land surface. As with all uncoupled modelling studies, a change in $g_{\mathrm{s}}$ and flux will impact the ${\mathrm{O}}_{\mathrm{3}}$ concentration itself. Therefore adopting the Medlyn formulation with a higher $g_{\mathrm{s}}$ values and subsequently higher ${\mathrm{O}}_{\mathrm{3}}$ flux for broadleaf and ${\mathrm{C}}_{\mathrm{3}}$ PFTs (Fig. 2) would lead to reduced ${\mathrm{O}}_{\mathrm{3}}$ concentration, which in turn may dampen the effect of higher $g_{\mathrm{s}}$ on ${\mathrm{O}}_{\mathrm{3}}$ flux. Although the higher uptake of ${\mathrm{O}}_{\mathrm{3}}$ by vegetation may lead to more damage and increase ${\mathrm{O}}_{\mathrm{3}}$ concentrations, in an uncoupled chemistry–land modelling system such as this it is not possible to predict which process would dominate. Additionally, this version of JULES does not have a crop module; it has no land management practices such as harvesting, ploughing or crop rotation – processes which may have counteracting effects on the land carbon sink. Further, without a coupled carbon and nitrogen cycle, it is likely that the ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation response of GPP and the land carbon sink is overestimated in some regions of our simulations since nitrogen availability limits terrestrial carbon sequestration of natural ecosystems in the temperate and boreal zone (Zaehle, 2013). This would have consequences for our modelled ${\mathrm{O}}_{\mathrm{3}}$ impact, particularly into the future when the large ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation effect was responsible for partly offsetting the negative impact of ${\mathrm{O}}_{\mathrm{3}}$. However, in our simulations a high fraction of land cover is agricultural, which we assume would be optimally fertilised. Our simulations also use a fixed climate, so we do not include the effect of climate change on shifting plant phenology. Therefore, our results may underestimate plant ${\mathrm{O}}_{\mathrm{3}}$ damage since if the growing season started earlier or finished later plants in some regions would be exposed to higher ${\mathrm{O}}_{\mathrm{3}}$ concentrations. Nevertheless, we emphasise that this study provides a sensitivity assessment of the impact of plant ${\mathrm{O}}_{\mathrm{3}}$ damage on GPP and the land carbon sink.

Another caveat we fully acknowledge is that at the leaf level JULES is parameterised to reduce $g_{\mathrm{s}}$ with ${\mathrm{O}}_{\mathrm{3}}$ exposure. Whilst this response is commonly observed (Wittig et al., 2007; Ainsworth et al., 2012), there is evidence to suggest that ${\mathrm{O}}_{\mathrm{3}}$ impairs stomata in some species, making them non-responsive to environmental stimuli (Hayes et al., 2012; Hoshika et al., 2012a; Mills et al., 2009; Paoletti and Grulke, 2010). In drought conditions the mechanism is thought to involve ${\mathrm{O}}_{\mathrm{3}}$-stimulated ethylene production, which interferes with the stomatal response to abscisic acid (ABA) signalling (Wilkinson and Davies, 2009, 2010). Such stomatal sluggishness can result in higher ${\mathrm{O}}_{\mathrm{3}}$ uptake and injury, increased water loss, and therefore greater vulnerability to environmental stresses (Mills et al., 2016). McLaughlin (2007a, b) and Sun et al. (2012) provide evidence of increased transpiration and reduced streamflow in forests at the regional scale in response to ambient levels of ${\mathrm{O}}_{\mathrm{3}}$ and suggest this could increase the frequency and severity of droughts. Hoshika et al. (2012b), however found that despite sluggish stomatal control in ${\mathrm{O}}_{\mathrm{3}}$-exposed trees, whole tree water use was lower in these trees because of lower gas exchange and premature leaf shedding of injured leaves. To our knowledge, the study of Hoshika et al. (2015) is the first to include an explicit representation of sluggish stomatal control in a land–atmosphere model. They show that sluggish stomatal behaviour has implications for carbon and water cycling in ecosystems. However, it is by no means a ubiquitous response, and it is not fully understood which species respond this way and under what conditions (Mills et al., 2016; Wittig et al., 2007). Nevertheless, this remains an important area of future work.

In this work we implement the stomatal closure proposed in Medlyn et al. (2011), which uses the parameter $g_{\mathrm{1}}$. Hoshika et al. (2013) show a significant difference in the $g_{\mathrm{1}}$ parameter (higher in elevated ${\mathrm{O}}_{\mathrm{3}}$ compared to ambient) in Siebold's beech in June of their experiment. However, this is only at the start of the growing season; further measurements show no difference in this parameter among ${\mathrm{O}}_{\mathrm{3}}$ treatments. Quantifying an ${\mathrm{O}}_{\mathrm{3}}$ effect directly on $g_{\mathrm{1}}$ would require a detailed meta-analysis of empirical data on photosynthesis and $g_{\mathrm{s}}$ for different PFTs, which is currently lacking in the literature.

A further caveat of this study is that the ${\mathrm{O}}_{\mathrm{3}}$ concentrations used to force the model are offline, in this case generated by the EMEP MSC-W model. This means the depositional sink is different in JULES (Medlyn formulation) compared to the EMEP model, which uses the $g_{\mathrm{s}}$ formulation presented in Emberson et al. (2000, 2001). Because we link two different model systems, the $g_{\mathrm{s}}$ values in the EMEP model differ from those obtained using the Medlyn formulation, which would ultimately lead to different ${\mathrm{O}}_{\mathrm{3}}$ concentrations. The role of EMEP in this study is to provide ${\mathrm{O}}_{\mathrm{3}}$ concentrations at the top of the vegetation canopy to force JULES and not $g_{\mathrm{s}}$; how the different depositional sinks would affect simulated ${\mathrm{O}}_{\mathrm{3}}$ concentrations at canopy height has not been investigated.

These offline simulations show the sensitivity of GPP and the land carbon sink to tropospheric ${\mathrm{O}}_{\mathrm{3}}$, suggesting that ${\mathrm{O}}_{\mathrm{3}}$ is an important predictor of future GPP and the land carbon store across Europe. There are uncertainties in our estimates, however, from the use of uncoupled tropospheric chemistry, meteorology and stomatal function. For example, increased frequency of drought in the future would reduce stomatal conductance (assuming no sluggish stomatal response) and thus ${\mathrm{O}}_{\mathrm{3}}$ uptake. Since our offline simulations do not include this feedback, it is possible the ${\mathrm{O}}_{\mathrm{3}}$ effect is overestimated here. Given the complexity of potential interactions and feedbacks, it remains difficult to diagnose the importance of individual factors (e.g. the direct physiological response) in a fully coupled simulation. Once the importance of a process is demonstrated offline, it provides evidence of the need to incorporate such a process in coupled regional and global simulations.

${\mathrm{O}}_{\mathrm{3}}$ as a missing component of carbon cycle assessments?

Comprehensive analyses of the European carbon balance suggest a large biogenic carbon sink (Janssens et al., 2003; Luyssaert et al., 2012; Schulze et al., 2009). However, estimates are hampered by large uncertainties in key components of the land carbon balance, such as estimates of soil carbon gains and losses (Ciais et al., 2010; Janssens et al., 2003; Schulze et al., 2009, 2010). We suggest that the effect of ${\mathrm{O}}_{\mathrm{3}}$ on plant physiology is a contributing factor to the decline in soil carbon stores observed across Europe, and as such this ${\mathrm{O}}_{\mathrm{3}}$ effect is a missing component of European carbon cycle assessments. Over the full experimental period (1901 to 2050), our results show elevated ${\mathrm{O}}_{\mathrm{3}}$ concentrations reduce the amount of carbon that can be stored in the soil by 3 to 9 % (low and high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivities, respectively), which almost completely offsets the beneficial effects of ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation on soil carbon storage under the high plant ${\mathrm{O}}_{\mathrm{3}}$ sensitivity . This would contribute to a change in the size of a key carbon sink for Europe and is particularly important when we consider the evolution of the land carbon sink into the future given the impact of ${\mathrm{O}}_{\mathrm{3}}$ on soil carbon sequestration and the high uncertainty of future tropospheric ${\mathrm{O}}_{\mathrm{3}}$ concentrations. Schulze et al. (2009) and Luyssaert et al. (2012) extended their analysis of the European carbon balance to include additional non-${\mathrm{CO}}_{\mathrm{2}}$ greenhouse gases (${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$). Both studies found that emissions of these offset the biogenic carbon sink, reducing the climate mitigation potential of European ecosystems. This highlights the importance of accounting for all fluxes and stores in carbon and greenhouse gas balance assessments, of which ${\mathrm{O}}_{\mathrm{3}}$ and its indirect effect on the ${\mathrm{CO}}_{\mathrm{2}}$ flux via direct effects on plant physiology are currently missing.

Interactive effects of ${\mathrm{O}}_{\mathrm{3}}$ and ${\mathrm{CO}}_{\mathrm{2}}$

We looked at the interactive effects of ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$. Our results support the hypothesis that elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ provides some protection against ${\mathrm{O}}_{\mathrm{3}}$ damage because of lower $g_{\mathrm{s}}$ that reduces uptake of ${\mathrm{O}}_{\mathrm{3}}$ through stomata (Harmens et al., 2007; Wittig et al., 2007). In the present study, reductions in GPP and the land carbon store due to ${\mathrm{O}}_{\mathrm{3}}$ exposure were lower when simulated with concurrent changes in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. Despite acclimation of photosynthesis after long-term exposure to elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ of field-grown plants (Ainsworth and Long, 2005; Medlyn et al., 1999), there is no evidence to suggest that $g_{\mathrm{s}}$ acclimates (Ainsworth et al., 2003; Medlyn et al., 2001). This suggests the protective effect of elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ against ${\mathrm{O}}_{\mathrm{3}}$ damage will be sustained in the long term. However, although meta-analysis suggests a general trend of reduced $g_{\mathrm{s}}$ with elevated ${\mathrm{CO}}_{\mathrm{2}}$ (Ainsworth and Long, 2005; Medlyn et al., 1999), this is not a universal response. Stomatal responses to exposure to elevated ${\mathrm{CO}}_{\mathrm{2}}$ with FACE treatment varied with genotype and growth stage in a fast-growing poplar community (Bernacchi et al., 2003; Tricker et al., 2009). In other mature forest stands, limited stomatal response to elevated ${\mathrm{CO}}_{\mathrm{2}}$ was observed after canopy closure (Ellsworth, 1999; Uddling et al., 2009). Also, some studies found that stomatal responses to ${\mathrm{CO}}_{\mathrm{2}}$ were significant only under high atmospheric humidity (Cech et al., 2003; Leuzinger and Körner, 2007; Wullschleger et al., 2002). These examples illustrate that stomatal responses to elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ are not universal, and as such the protective effect of ${\mathrm{CO}}_{\mathrm{2}}$ against ${\mathrm{O}}_{\mathrm{3}}$ injury cannot be assumed for all species at all growth stages under wide-ranging environmental conditions.

Conclusions

What is abundantly clear is that plant responses to both ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ are complicated by a host of factors that are only partly understood, and it remains difficult to identify general global patterns given that effects of both gases on plant communities and ecological interactions are highly context and species specific (Ainsworth and Long, 2005; Fuhrer et al., 2016; Matyssek et al., 2010b). This study quantifies the sensitivity of the land carbon sink for Europe and GPP to changing concentrations of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ from 1901 to 2050. We have used a state-of-the-art land surface model calibrated for European vegetation to give our best estimates of this sensitivity within the limits of data availability to calibrate the model for ${\mathrm{O}}_{\mathrm{3}}$ sensitivity, current knowledge and model structure. In summary, this study has shown that potential gains in terrestrial carbon sequestration over Europe resulting from elevated ${\mathrm{CO}}_{\mathrm{2}}$ can be partially offset by concurrent rises in tropospheric ${\mathrm{O}}_{\mathrm{3}}$ over 1901–2050. Specifically, we have shown that the negative effect of ${\mathrm{O}}_{\mathrm{3}}$ on the land carbon sink was greatest over the twentieth century, when ${\mathrm{O}}_{\mathrm{3}}$ concentrations increased rapidly from pre-industrial levels. Over this period soil carbon stocks were diminished over agricultural areas, consistent with reduced NPP and litter input. This loss of soil carbon was largely responsible for the decrease in the size of the land carbon sink over Europe. The ${\mathrm{O}}_{\mathrm{3}}$ effect on the land carbon store and flux was reduced into the future as ${\mathrm{CO}}_{\mathrm{2}}$ concentration rose considerably and changes in ${\mathrm{O}}_{\mathrm{3}}$ concentration were less pronounced. However, there remained a large cumulative negative impact on the land carbon sink for Europe by 2050. The interaction between the two gases was found to reduce ${\mathrm{O}}_{\mathrm{3}}$ injury owing to reduced stomatal opening in elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. However, primary productivity and land carbon storage remained suppressed by 2050 due to plant ${\mathrm{O}}_{\mathrm{3}}$ damage. Expressed as a percentage of the emissions from fossil fuel and cement production for the EU28-plus countries, the carbon emissions from ${\mathrm{O}}_{\mathrm{3}}$-induced plant injury are a source of anthropogenic carbon previously not accounted for in carbon cycle assessments. Our results demonstrate the sensitivity of modelled terrestrial carbon dynamics to the direct effect of tropospheric ${\mathrm{O}}_{\mathrm{3}}$ and its interaction with atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ on plant physiology, demonstrating this process is an important predictor of future GPP and trends in the land carbon sink. Nevertheless, this process remains largely unconsidered in regional and global climate model simulations that are used to model carbon sources and sinks and carbon–climate feedbacks.

The JULES model can be downloaded from the Met Office Science Repository Service (https://code.metoffice.gov.uk/trac/jules, last access: 11 September 2017 – see here for a helpful how-to http://jules.jchmr.org/content/getting-started, last access: 11 September 2017). Model output data presented in this paper and the exact version of JULES with name lists are available upon request from the corresponding author.

The Supplement related to this article is available online at https://doi.org/10.5194/bg-15-4245-2018-supplement.

RJO performed the JULES model simulations and analysis. DS provided regional tropospheric ozone forcing from the EMEP MSC-W model. All authors wrote the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

Rebecca J. Oliver and Lina M. Mercado were supported by the EU FP7 (ECLAIRE, 282910) and JWCRP (UKESM, NEC05816). This work was also supported by EMEP under UNECE. Stephen Sitch and Lina M. Mercado acknowledge the support of the NERC SAMBBA project (NE/J010057/1). The UK Met Office contribution was funded by BEIS under the Hadley Centre Climate Programme (GA01101). Gerd A. Folberth also acknowledges funding from the EU's Horizon 2020 research and innovation programme (CRESCENDO, 641816). We also thank Magnuz Engardt of SMHI for providing the RCA3 climate dataset. This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CARBOMONT, ChinaFLUX, FLUXNET Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, TERN OzFlux, TCOS Siberia and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonisation was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project and Fluxdata project of FLUXNET, with the support of CDIAC and the ICOS Ecosystem Thematic Centre and the OzFlux, ChinaFlux, and AsiaFlux offices. We also thank the two anonymous reviewers, who helped to improve this paper. Edited by: Xinming Wang Reviewed by: two anonymous referees