Interactive comment on “ Ozone vegetation damage effects on gross primary productivity in the United States ”

Yue and Unger provide an interesting article that estimates the impact of ozone on GPP in the US. The authors find that the inclusion of ozone reduces GPP in the US, primarily in the eastern half of the country. The ozone concentrations and GPP are evaluated against measurements, and the sensitivity of using different land cover datasets and meteorological datasets was also tested. Overall, this is a relatively clear and very thorough analysis of GPP responses and model sensitivities, though there are a few problems with methodology and analyses that need to be clarified.


Introduction
The effects of tropospheric ozone (O 3 ) damage on US forests have been studied for half a century (Karnosky et al., 2007), but the impacts of O 3 on the North American carbon balance are still relatively poorly understood (Felzer et al., 2004;Huntingford et al., 2011).O 3 is a secondary pollutant produced in the troposphere during the photochemical oxidation of carbon monoxide, methane, and volatile organic com-pounds (VOCs) by the major tropospheric oxidant, the hydroxyl radical, in the presence of sunlight and nitrogen oxides.Fossil-fuel, biofuel and biomass burning since the industrial and agricultural revolutions have greatly increased the emissions of O 3 precursors and led to an approximate doubling of O 3 levels over the US since the preindustrial.Deposition through stomatal uptake is an important sink for O 3 but damages photosynthesis, reduces plant growth and biomass accumulation, limits crop yields, and affects stomatal control over plant transpiration of water vapor between the leaf surface and atmosphere (Ainsworth et al., 2012;Hollaway et al., 2012).
Understanding the O 3 pollution influence on the North American forest sink is crucial to any effort to mitigate climate change by stabilizing atmospheric carbon dioxide (CO 2 ) concentrations.Currently, North America is acting as a net source of CO 2 to the atmosphere (King et al., 2012), mainly due to fossil-fuel combustion in the US and the deforestation in Mexico.Sequestration of atmospheric CO 2 by forest ecosystems is a major control on atmospheric CO 2 abundance and its growth rate (Pan et al., 2011).Terrestrial ecosystems of North America absorb the equivalent of about 35 % of North America's fossil fuel based CO 2 emissions, representing a source-to-sink ratio of nearly 3 : 1. Forest regrowth in the US is responsible for 30-70 % of this North American CO 2 sink, which varies significantly from year to year (Pacala et al., 2001;Goodale et al., 2002;Pan et al., 2011;King et al., 2012).However, O 3 damage may in part dampen the level of carbon sequestration by North American ecosystems (Felzer et al., 2004(Felzer et al., , 2005)).
Experimental studies that examine O 3 impacts on plant productivity are typically performed for individual vegetation types, on the scale of sites, and within a limited time Published by Copernicus Publications on behalf of the European Geosciences Union.
period (e.g., Wittig et al., 2007;Feng et al., 2008;Lombardozzi et al., 2013).For example, based on measurements reported from over 100 studies, Wittig et al. (2007) estimated that chronic O 3 exposure depressed photosynthesis by 11 % and stomatal conductance by 13 % for several tree species at the ambient O 3 level of ∼ 45 ppbv relative to that in O 3 -free air.The O 3 damage effect is strongest for crops.With data sets from ∼ 50 peer-reviewed studies, Feng et al. (2008) estimated that elevated O 3 levels significantly decrease wheat photosynthetic rates by 20 % and stomatal conductance by 22 %.Emerging research has found that the O 3 vegetation damage effects may result in a loss of plant stomatal control, and a consequent decoupling of the stomatal response from photosynthesis inhibition (Lombardozzi et al., 2012a(Lombardozzi et al., , b, 2013)).
Previous work has found that in the US region during 1989-1993, O 3 pollution reduced net primary productivity (NPP) by 3-7 % overall, and up to 13 % in hot spots including the southeast and in the Midwest agricultural lands (Felzer et al., 2004(Felzer et al., , 2005)).The indirect CO 2 radiative forcing due to the vegetation damage effects of anthropogenic O 3 increases since the industrial revolution may be as large as +0.4 Wm −2 (Sitch et al., 2007), which is 25 % of the magnitude of the direct CO 2 radiative forcing over the same period, and of similar magnitude to the direct O 3 radiative forcing.Through this perturbation of the carbon cycle, O 3 pollution affects the climate system on considerably longer timescales than its own atmospheric lifetime (Unger and Pan, 2012).Over the past decade since this previous assessment surface O 3 levels in most of the US have decreased (Lefohn et al., 2010) due to domestic emission reductions following the implementation of air quality control legislation (Bloomer et al., 2010).However, increasing O 3 concentration is observed over western US (Jaffe and Ray, 2007).Such a trend may in part be related to the inter-continental flow from Asia (Cooper et al., 2010) and the global increase in methane (Rigby et al., 2008).
The major goal of this study is to assess O 3 damage effects on gross primary productivity (GPP) in the US for the recent decade 1998-2007 using a data-constrained vegetation model.In this work, we describe the implementation of a semi-mechanistic O 3 damage function (Sitch et al., 2007) into the Yale Interactive Terrestrial Biosphere model (YIBs) that includes enzyme-kinetic biophysics (Unger et al., 2013).In the first stage of the study, we utilize eddy-derived GPP flux measurements at 40 sites across the US and Canada that have been collated for the North American Carbon Program (NACP) site-level interim synthesis (Huntzinger et al., 2012;Schaefer et al., 2012;Barr et al., 2013;Ricciuto et al., 2013) to evaluate an off-line version of the vegetation model's site level GPP simulation and to assess the impact of surface O 3 damage at those sites.In the second stage of the study, the impacts of O 3 damage on GPP throughout the entire US are quantified using a regional configuration of the vegetation model.

Vegetation biophysics
Here, we apply an off-line version of the YIBs model that previously was implemented into the NASA Goddard Institute for Space Studies global chemistry-climate model (Unger et al., 2013).The off-line model can be run at the site-level or in regional mode for a designated region.The vegetation biophysics module computes the photosynthetic uptake of CO 2 coupled with the transpiration of water vapor at the 1 h physical integration time step of the off-line model.The vegetation biophysics calculates C3 and C4 photosynthesis using the well-established Michealis-Menten enzymekinetics leaf model of photosynthesis (Farquhar et al., 1980;von Caemmerer and Farquhar, 1981) and the stomatal conductance model of Ball and Berry (Collatz et al., 1991).The coupled photosynthesis/stomatal conductance leaf model has been widely used to project terrestrial biosphere responses to global change.The model is briefly summarized here for transparency and completeness.The leaf model assumes that the rate of net CO 2 assimilation (A net ) in the leaves of C3 and C4 plants is limited by one of three processes: (i) the capacity of the ribulose 1,5-bisphosphate (RuBP) carboxylaseoxygenase enzyme (Rubisco) to consume RuBP (J c ); (ii) the capacity of the Calvin cycle and the thylakoid reactions to regenerate RuBP supported by electron transport (J e ); (iii) the capacity of starch and sucrose synthesis to consume triose phosphates and regenerate inorganic phosphate for photophosphorylation in C3 and phosphoenolpyruvate (PEP) limitation in C4 (J s ).J c , J e , and J s are described as functions of the maximum carboxylation capacity (V cmax ) at the optimal temperature, 25 • C, and the internal leaf CO 2 concentration (C i ).The gross rate of carbon assimilation from photosynthesis (A) is given by the following: Net carbon assimilation is given by the following: where R d is the rate of dark respiration: Leaf stomata control the uptake of CO 2 vs. the loss of H 2 O.At equilibrium, the stomatal conductance of water vapor through the leaf cuticle (g s in mol [H 2 O] m −2 s −1 ) depends on the net rate of carbon assimilation: where m and b are the slope and intercept derived from empirical fitting to the Ball and Berry stomatal conductance equations, RH is relative humidity, c s is the CO 2 concentration at the leaf surface, and r s is the stomatal resistance  (Friend and Kiang, 2005) and the Community Land Model (Oleson et al., 2010) with updates from Bonan et al. (2011) (Table 1).In both the site-level and regional models, we apply these model PFTspecific photosynthesis parameters and do not tune or calibrate to the local vegetation properties.The model calculates evapotranspiration as a function of the stomatal conductance.However, we do not consider the feedback of the changes in evapotranspiration to the boundary-layer meteorology because we use prescribed meteorological variables from reanalysis in the simulations.The canopy radiative transfer scheme assumes a closed canopy and layers the canopy for light stratification using an adaptive number of layers (typically 2-16) (Friend and Kiang, 2005).Each canopy layer distinguishes sunlit and shaded regions for which the direct and diffuse photosynthetically active radiation (PAR) is computed (Spitters et al., 1986).The coupled photosynthesis and stomatal conductance equations are solved analytically using a cubic function of A net .C i is calculated explicitly at the leaf level.Scaling of the leaf to canopy level is through stratification of canopy light levels and leaf area profiles.The photosynthetic uptake of CO 2 is accumulated into a carbon reserve pool, from which other processes may allocate uses.

O 3 damage effect on photosynthesis
O 3 oxidizes cellular membranes and photosynthetic tissues when it enters leaves through stomata, leading to reductions in photosynthesis and GPP.O 3 damage inhibits stomatal conductance, which is closely related to the photosynthetic rate, resulting in a reduction in transpiration.A semi-mechanistic parameterization is employed to estimate the O 3 damage effects to both photosynthesis and stomatal conductance (Sitch et al., 2007).The exposure to O 3 leads to reductions in photosynthesis: where F is the reduction fraction calculated as where a is the O 3 sensitivity coefficient derived from observations.Two cases are examined: high and low O 3 sensitivity following Sitch et al. (2007).U >O3T is the instantaneous leaf uptake of O 3 flux above a plant function type (PFT)-specific threshold of O3T (Table 1), Here F O 3 is the O 3 flux entering the leaf through the stomata, where [O 3 ] is the O 3 concentration at the top of the canopy, and r b is the boundary layer resistance.The stomatal resistance to O 3 is calculated based on stomatal resistance to water r s with a ratio constant κ = 1.67.From Eq. ( 4), the decrease in A net reduces the stomatal conductance g s proportionally, The r s and g s are the O 3 -damaged stomatal resistance and conductance, respectively.When the plant is exposed to [O 3 ] (Eq.8), the excess O 3 flux entering leaves (Eq.7) causes F < 1 (Eq.6), decreasing A net (Eq. 5) while increasing the stomatal resistance (Eq.9).The latter will act to reduce the O 3 uptake flux (Eq.8) to protect the plant.Thus, the scheme considers associated changes in both photosynthetic rate and stomatal conductance.When photosynthesis is inhibited by O 3 , the stomatal conductance decreases accordingly to resist more air passing through the stomata, resulting in a decline of the oxidant fluxes inside leaves, as described through Eqs. ( 5

Vegetation structure
The YIBs vegetation model simulates eight PFTs, using either C3 or C4 photosynthesis (Table 1).We apply two different sets of land cover and leaf area index (LAI) in the simulations.The first set is the PFT-specific vegetation cover fraction and LAI retrieved by the Moderate Resolution Imaging Spectroradiometer (MODIS, Knyazikhin et al., 1998).
The value on a specific day is linearly interpolated from the monthly means of the nearest two months based on the distance of this day to the middle dates of those two months.a Plant function types (PFTs) are tundra (TDA), C3 grassland (GRAC3), C4 savanna/grassland (GRAC4), shrubland (SHR), deciduous broadleaf forest (DBF), evergreen needleleaf forest (ENF), tropical rainforest (TRF), and cropland (CRO).b For site-level simulations, we consider CRO to be a C4 plant.For regional simulation, we consider half CRO as C3 plants (soybean) and the rest C4 plant (corn).

Meteorological forcing
For the site-level simulations, we use hourly in situ measurements of surface meteorological variables, including surface air temperature, specific humidity, wind speed, surface pressure, and CO 2 concentrations.There are some missing values in the measurements due to occasional instrument failure.We gap-fill the site-based observations with that from the MERRA-land data (Reichle et al., 2011), which is interpolated to each site based on the site location.
For the regional simulations, the off-line YIBs model uses hourly MERRA-land data climatic variables including the following: surface air temperature, specific humidity, wind speed, surface pressure, precipitation, direct PAR, and diffuse PAR, and soil temperature and soil moisture at six soil depths.The original data resolution of 0.5 • × 0.667 • by latitude and longitude is degraded to 1 • × 1.333 • due to current disk space limitation.

Surface [O 3 ]
Hourly and daily maximum 8 h average surface [O 3 ] representative of the present day climate (∼ 2005) are taken from previous simulations using NASA Model-E2 (Shindell et al., 2013).The global model has 2 • × 2.5 • latitude by longitude horizontal resolution with 40-vertical layers extending to 0.1 hPa.The gas-phase chemistry and aerosol modules are fully integrated, so that these components interact with each other and with the physics of the climate model (Bell et al., 2005;Shindell et al., 2006Shindell et al., , 2013;;Unger, 2011).The model surface O 3 is validated using measurements from 73 Clean Air Status and Trends Network (CASTNET) sites operated by the United States Environmental Protection Agency (EPA) (http://epa.gov/castnet/javaweb/ozone.html) and ∼ 1200 monitor sites managed by the EPA AIRDATA (http://www.epa.gov/airdata/).These sites are operated on the county level scale.The CASTNET provides hourly [O 3 ] at rural sites from 1996-2005.The AIRDATA network provides daily maximum 8 h average (MDA8) [O 3 ], covering both urban and rural regions.We use AIRDATA data for the year 2005.

Site-level runs
We configure a site-level version of the YIBs model for the 40 eddy covariance flux tower sites described in detail in the NACP synthesis (Fig. S1 in the Supplement and Appendix Table A1, Schaefer et al., 2012).Meteorological measurements are available for a wide range of time periods across the different sites ranging from the minimum of 1 year at Fermi Lab (US-IB1) and the maximum of 15 years at Harvard Forest (US-HA1).These sites cover a range of different vegetation types including the following: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands, croplands, closed shrublands, mixed forests, permanent wetlands, and woody savannas.Table S1 in the Supplement summarizes how the NACP vegetation types are mapped onto the eight model PFTs.For the site-level simulations, we assume C4 photosynthetic pathway for all cropland sites, which are mainly corn (Schaefer et al., 2012).The local site LAI values are not available.As a result, we use the MERRA or MODIS LAI for the simulations.
For each site, a group of six sensitivity simulations are performed (Table 2).We conduct the first four runs using different combinations of meteorological and vegetation forcings,  1.
To quantify the performance of the vegetation model, we estimate the χ 2 for each site following the method described in Schaefer et al. (2012), where is the difference between the pair of simulated and observed GPPs.ε i are the observational uncertainties resulting from turbulence, gap-filling, flux partitioning, and u * threshold determination (Barr et al., 2013).n is the length of observations (e.g., the number of days for the daily variables).The lower the χ 2 , the smaller the model biases.If χ 2 < 1, the simulation bias is on average smaller than the measurement uncertainty, indicating a good performance of the model.Here, we define a reasonable performance of χ 2 < 4, when the residual is less than twice the measurement uncertainty.We also calculate the root mean square error (RMSE) as follows: We validate the simulated O 3 damage effect with measurements from literature.Ishii et al., 2004;Zhang et al., 2012).Second, the diurnal cycles and seasonality of [O 3 ] are very different for different sites (Bloomer et al., 2010), making it difficult to apply a uniform temporal cycle for all the NACP sites.The reductions in GPP at these simulations are compared with results from field measurements at the corresponding [O 3 ] level.

Regional run over US
A gridded version of the YIBs model at 1 • × 1.333 • latitude by longitude horizontal resolution for the US region is driven with MERRA meteorological forcings for the period 1998-2007.In the regional model, vegetation cover types are from the ISLSCP and LAI is form the MERRA-land reanalysis.We assign the MERRA LAI to the corresponding PFT types defined by ISLSCP (Fig. S2 in the Supplement).The 18 ISLSCP land types are converted to 8 PFTs used in the model (Table S1 in the Supplement).Some of the ISLSCP land types, such as the deciduous needleleaf forest, are not represented in the YIBs model.However, the coverage of these types is very small in the US (Fig. S2 in the Supplement) and will not influence the regional simulation after the conversion to the model types.For the regional simulation, we assume that the total crop area in each crop grid cell is split 50 % C3 and 50 % C4 to account for the dominance of both soybean (C3) and corn (C4) crops in the central and northern US agricultural regime.We perform two simulation cases with high and low O 3 damage sensitivity.Finally, to understand how the O 3 vegetation damage effect may respond to possible future changes in [O 3 ], we perform four additional sensitivity experiments with ±25 % changes in [O 3 ] for each O 3 sensitivity case.

Evaluation of O 3 -free GPP at NACP sites
We first compare the monthly mean LAI from MERRA and MODIS at each NACP site (Fig. S3 in the Supplement).For each site, the MERRA LAI is averaged for the period when GPP measurements are available.The two data sets show similar annual cycles at several sites but are inconsistent for 7 out of 20 ENF sites (CA-Ca1, CA-Ca2, CA-Ca3, CA-NS1, US-Me2, US-Me3, and US-Me5) and 2 out of 5 shrubland sites (US-SO2 and US-Ton).In addition, for grasslands and croplands, the data sets exhibit different seasonality.It must be emphasized that the MERRA and MODIS LAI represent the average state in the retrieval product grid cells and as such  may not represent the local LAI for the actual PFT species at each site.
The long-term monthly mean O 3 -free GPP from the simulation METsite_LAImerra is compared with observations at NACP sites (Figs. 1 and S4 in the Supplement).The simulations capture the magnitude and seasonality of GPP for most sites especially for deciduous broadleaf and cropland.The largest model overestimate (factor of 3-8) occurs at two ENF sites, CA-SJ1 and CA-SJ2 in North Central Canada (Fig. S4 in the Supplement).For the grassland sites, the model-observation correlation is low because the seasonality is not well captured, especially for US-ARM (in Great Plains) and US-Var (in California), where the modeled maximum GPP occurs in summer (July), 2-3 months later than in the measurements (April) (Fig. S4 in the Supplement).This incorrect model seasonality is a result of the MERRA LAI (compare Fig. S3 in the Supplement) that does not begin to increase rapidly until after May and is not consistent with the local LAI at the site.In reality, California grasslands exhibit rapid growth in spring then mature and die after April or May (Chiariello, 1989).The grasslands in the Great Plains may have up to six different phenological groups, including some species active in spring (e.g., in US-ARM) while some others peak in summer (e.g., in US-Shd) (Henebry, 2003).On the annual mean basis, the correlation coefficient between simulated GPP and observations at all 40 sites is 0.65.The correlation is higher (0.81) for summer (June-August, Fig. S5 in the Supplement).The annual GPP averaged over all 40 sites is 3.8 g C m −2 day −1 , 27 % higher than the observational average (3.0 g C m −2 day −1 ).
Among the 40 NACP sites, 22 have reasonable performance with χ 2 < 4 for the simulation METsite_LAImerra (Fig. 2a).For these sites, 12 are ENF with the best (χ 2 = 1.2) at site US-Dk3 (Fig. S4 in the Supplement).The ENF sites usually have multiple years of measurements and provide good samples for testing the consistency between simulations and observations.Simulations at other four DBF, three cropland, and three shrubland sites have χ 2 < 4. Compared with the 24 land surface models in Schaefer et al. (2012), the YIBs model shows significant improvement at the crop PFT sites (χ 2 < 4.1 vs. χ 2 > 6).YIBs simulates GPP with χ 2 < 4 at 22 sites in total compared to 16 sites for the ensemble simulations in Schaefer et al. (2012).YIBs GPP simulation is weaker (χ 2 > 4) at 18 sites including eight ENF sites, two DBF, and two shrubland PFT sites.The common feature of the biases at these sites is the overestimation of peak GPP during summer (e.g., CA-SJ1, CA-SJ2, CA-Mer, Fig. S4 in the Supplement).It is possible that the model does not represent the full realism of the biophysical processes accurately.However, we assert that the most likely cause of the model overestimate is the uniform application of model For example, sites CA-SJ1, CA-SJ2, and CA-SJ3 are located close to each other.Simulations at these sites have similar magnitude in simulated GPP while observations show distinct variability between the sites.We compare R 2 , RMSE, and χ 2 for the different sensitivity experiments in order to ascertain which combination of meteorological and LAI forcings best reproduces the measured GPP over North America (Table 3 and Fig. 3).The sites CA-Let, CA-NS1, US-Var, CA-SJ1, and CA-SJ2 are excluded from the analysis because of excessive bias(Fig.S6a in the Supplement).The average R 2 increases while RMSE decreases when MERRA reanalyses are substituted with sitebased meteorology, or the MERRA LAI is used instead of MODIS LAI (Table 3).The choice of LAI forcing has the most significant impact on YIBs simulation performance, consistent with previous work that showed the dominance of phenology over meteorology in controlling local terrestrial carbon exchange (Desai et al., 2008;Puma et al., 2013).Using MODIS LAI, YIBs has a total χ 2 of 9.2 that shows an average reduction of 4.7 (52 %) with MERRA LAI (Table 3 and Fig. 3).Applying the site meteorology relative to MERRA meteorological forcings offers smaller improvements.For example, the total χ 2 value decreases by 5 % in MET-site_LAImodis compared with that in METmerra_LAImodis and 15 % in METsite_LAImerra relative to that in MET-merra_LAImerra (Table 3).

Evaluation of modeled surface [O 3 ]
We validate summertime surface O 3 simulated by the NASA Model-E2 chemistry-climate model with observations from the CASTNET and AIRDATA (Fig. 4).High O 3 level appears in the eastern US due to anthropogenic emissions and in the mountainous western US due to high elevation.The model generally captures this spatial pattern with a correlation coefficient of 0.39 against observations over the selected 73 CASTNET sites (Fig. 4a-b).The simulation overestimates the O 3 level by ∼ 4 ppbv (12 %) in the east and ∼ 1 ppbv (3 %) in the west.The CASTNET sites are located in rural sites, which usually have lower [O 3 ] than that in urban areas, except for some megacities where the excessive NO x emissions result in lower O 3 level (Gregg et al., 2003).Therefore, we also compare the simulated MDA8 [O 3 ] with monitored at ∼ 1200 AIRDATA sites, which covers both urban and rural regions (Fig. 4c).In the eastern US, the model captures high [O 3 ] centers around Michigan, Indiana, and Ohio states and that along the northeast coast.In the west, the simulation reproduces high [O 3 ] over mountain regions and in California.On average, the simulated MDA8 [O 3 ] is lower by ∼ 0.5 ppbv (1 %) in the east and ∼ 3.5 ppbv (7 %)   S6a in the Supplement).We also calculate the mean values of R 2 and RMSE for these sites.We calculate the total χ 2 (shown as red bars in Fig. 3) using all the available observations over all sites.in the west.The correlation coefficient between simulations and observations is as high as 0.51 (Fig. 4d).

O 3 damage effects at NACP sites
We must apply the simulated O 3 to quantify the O 3 vegetation damage at the NACP sites because the sites do not monitor local [O 3 ].The summer average [O 3 ] is 30-50 ppbv at 24 US sites (Fig. 5a).The O 3 damage effect is relatively stronger at sites with both high O 3 -free GPP and ambient [O 3 ] (Fig. 5d).The most significant damages are predicted at US-MMS (DBF) and US-Dk3 (ENF) sites, where the GPP reductions are 5-14 % depending on the low or high O 3 sensitivity (Fig. 5d).At these two sites, the high stomatal conductance (4.0 and 3.4 mm s −1 , Fig. 5b) and ambient [O 3 ] (both 43 ppbv, Fig. 5a) result in the largest O 3 stomatal flux (both ∼ 0.3 mmol m −2 d −1 , Fig. 5c) among the 24 sites.The lowest O 3 damage (1-2 % GPP reduction) appears in the three shrub sites, US-Ton, US-SO2, and US-Los, although mean [O 3 ] there is as high as 43 ppbv.The main reason for the limited O 3 damage is the low stomatal conductance (average 1.4 mm s −1 , Fig. 5b) related to the small O 3 -free GPP (average 4.6 g C m −2 d −1 , Fig. 5d).Similarly, the O 3 damage for C3 grass is as low as 1-2 %, although the GPP of this plant is highly sensitive to O 3 (Table 1).For ENF and DBF sites, the average site-level O 3 damage effects are estimated to be 2-5 and 3-9 % respectively with differences between these ecosystem types predominantly driven by differences in sensitivity to O 3 .The four C4 crop sites, US-Ne1, US-IB1, US-Ne2, and US-Ne3, exhibit the highest O 3 -free GPP but show only moderate O 3 damage effects (GPP reductions of 4-6 %, Fig. 5d).This result is driven by low ambient [O 3 ] at the C4 crop sites (average 32 ppbv, Fig. 5a) in combination with the reduced C4 stomatal conductance (higher water use efficiency) relative to C3 plants (average 3.2 mm s −1 , Fig. 5b).Indeed, the C4 photosynthetic pathway has been observed to offer protection against O 3 damage (Heagle et al., 1989;Rudorff et al., 1996).Inclusion of O 3 damage effect improves the site-level simulations (Fig. 2b-c).For 36 out of the 40 sites, the χ 2 of simulated GPP decreases when considering vegetation responses to O 3 , and the improvement is better when higher O 3 sensitivity is applied.At these sites, for example, CA-TP4, US-Dk3, US-MMS, and US-PFa, the reduced GPP at peak seasons is closer to measurements (not shown), leading to smaller biases for simulations.On average, the χ 2 decreases by 3-8 % at these sites, depending on the O 3 sensitivity in the simulation (Table 3 and Fig. 3).Finally, the simulated annual GPP averaged over all NACP sites changes from 3.8 g C m −2 day −1 to 3.6 g C m −2 day −1 with the high O 3 sensitivity simulation case, closer to the observations of 3.0 g C m −2 day −1 .The bias-correction from O 3 damage is much smaller relative to the effect of phenology (Fig. 3).Moreover, the O 3 -induced damage does not improve the GPP correlation between observations and simulations, which remains similar at ∼ 0.8 (for 40 sites) with and without O 3 effects (Fig. S5 in the Supplement).(1996) for colluna vulgaris, Mulchi et al. (1992) for soybean, and Ishii et al. (2004) and Ainsworth (2008) for rice.Values for rice are denoted in green and others in red.The author initials are indicated for the corresponding studies.

Evaluation of simulated O 3 vegetation damage against field and laboratory data
We compare the simulated O 3 damage effect with field and laboratory measurements from the published literature (Fig. 6).In total, 14 additional sensitivity experiments are performed with different levels of [O 3 ] at each NACP site (see Sect. 2.2.1).GPP reductions increase accordingly with the increasing [O 3 ] (Fig. 6).For a given [O 3 ], the O 3 damage effect is strongest for C4 crops (despite the lower g s : A net ratio) but weakest for shrubland.YIBs simulates reasonable O 3 damage to GPP for all model PFTs compared to the meta-analyses of Wittig et al. (2007) and Lombardozzi et al. (2013).Field studies in shrubland are limited.Zhang et al. (2012) investigated the responses of four shrub species to [O 3 ] = 70 ppbv and found large reductions in net photosynthesis of 50-60 %.The average O 3 -free A net of those shrub species was 8-16 g [C] m −2 s −1 , much higher than even the gross photosynthesis (A) of 6 g [C] m −2 s −1 at the shrub NACP sites, likely because the latter are located in dry and/or cold areas (Fig. S1 in the Supplement).The YIBs simulated O 3 vegetation damage effects for C4 plants are in good agreement with field measurements from Taylor et al. (2002) and Grantz et al. (2012).In the case of C3 grass and C3 crop, the model simulates consistent GPP reduction percentages with observations from Feng et al. ( 2008) for wheat, Foot et al. (1996) for colluna vulgaris, and Mulchi et al. (1992) for soybean.However, these O 3 damage results are all > 50 % less than for available measurements in rice crops (Ishii et al., 2004;Ainsworth et al., 2008), suggesting that rice may have much higher O 3 sensitivity than other C3 plants.In the US rice plantation area is much smaller than that of soybean and corn.Therefore, we adopt the O 3 sensitivity parameters for C3/C4 plants shown in Table 1 for the regional simulations.

O 3 vegetation damage effect on GPP in US region
High values of simulated summertime GPP (including O 3 damage effect) appear east of 95 • W in the US (Fig. 7a),  because the land surface there is covered by crops and forests.A high center of GPP (> 10 g C m −2 day −1 ) appears over cropland in the north central US In the western US, the coverage of grass and shrub and the low water availability (low precipitation and soil moisture) over semi-arid regions lead to a low carbon-assimilation rate.The regional gridded simulated GPP reproduces the JJA growing season average NACP site-level fluxes with a correlation coefficient of 0.62 for 32 sites below 50 • N (points in Fig. 7a).The correlation is lower than the 0.84 estimated for the site-level simulation METsite_LAImerra at the same sites and the same season.

Atmos
Since the meteorological forcings and LAI are similar, the difference in land cover, ISLSCP vs. site definitions (Figs.S1 and S2 in the Supplement), accounts for the discrepancy between regional and site-level simulations.
On average, the simulated summer GPP (including the high O 3 damage effect) is 9.5 g C m −2 day −1 in the eastern US (east of 95 • W) and 3.9 g C m −2 day −1 in the western US (west of 95 • W), giving a mean value of 6.1 g C m −2 day −1 for the US region.The total carbon uptake is estimated to be 4.43 ± 0.18 Pg C during the summer growing season, accounting for 57-60 % of the annual average value of 7.59 ± 0.25 Pg C over the 1998-2007 period.Our estimate of annual carbon uptake is consistent with previous published US Similar reduction fractions are predicted for the annual GPP.
US surface [O 3 ] exhibits a decreasing trend over the past 2 decades, especially in the eastern US, due to precursor emission controls (Lefohn et al., 2010).However, the community continues to debate how surface [O 3 ] will respond to future emissions and climate change.On the one hand, surface [O 3 ] may decline by the mid 21st century due to large reductions in regional anthropogenic precursor emissions (Wu et al., 2008).On the other hand, climate change effects alone may increase local surface [O 3 ] due to the warmer, drier, and more stable environment (Leibensperger et al., 2008;Wu et al., 2008).Due to the uncertainty in future surface [O 3 ] projections, our strategy here is to perform four additional sensitivity experiments with ±25 % changes in [O 3 ] for each O 3 sensitivity case.Increases of 25 % in [O 3 ] may reduce GPP in the eastern US by 6-11 %, with a maximum local reduction of 25 % for the high O 3 sensitivity case (Fig. 9d).The damage magnitude with low O 3 sensitivity (Fig. 9b) mimics the present-day estimate with high O 3 sensitivity (Fig. 8b).In contrast, the O 3 damage to the eastern US GPP is as low as 2-4 % in response to 25 % decreases in [O 3 ] (Fig. 9a, c), suggesting a substantial co-benefit to ecosystem-health of O 3 precursor emissions control.

Conclusions
We have performed an updated assessment of O 3 vegetation damage effects on GPP in the US for the period using the YIBs vegetation model.The semi-mechanistic parameterization of O 3 inhibition on photosynthesis proposed by Sitch et al. (2007) has been implemented into this processbased vegetation model.The simulated O 3 damage effects are consistent with laboratory and field measurements reported in previously published studies.We evaluated the simulated O 3 -free and O 3 -damaged GPP with in situ measurements from 40 NACP sites.The O 3 -free and O 3 -damaged GPP simulations capture the seasonality and interannual variability of GPP at most sites.The model GPP biases are lowest at forest and cropland sites but highest at grassland sites.Model GPP is highly sensitive to choice of LAI forcing.Simulations that apply MERRA LAI generally perform better (show lower biases) than those with MODIS LAI.In response to the simulated ambient [O 3 ] of 30-50 ppbv, simulated GPP decreases by 1-14 % at the NACP sites, depending on the O 3 sensitivity and PFT type.Maximum reductions of 5-14 % occur in two forest sites, where both O 3 -free GPP and ambient [O 3 ] are relatively high.Inclusion of the O 3 damage offers only a small improvement to the simulated annual average GPP at NACP sites (from 3.8 g C m −2 day −1 to 3.6 g C m −2 day −1 ) such that the model still overestimates the observational average of 3.0 g C m −2 day −1 .The model GPP overestimate is most likely related to the use of generic PFT-specific photosynthesis parameters and the satellite prescribed LAI that may not represent the local site LAI.In this work, we assumed a coupled response between photosynthesis and stomatal conductance.Emerging research has found that the O 3 vegetation damage effects can actually result in a loss of plant stomatal control, and a consequent Atmos.Chem.Phys., 14, 9137-9153, 2014 www.atmos-chem-phys.net/14/9137/2014/decoupling of the stomatal response from photosynthesis inhibition (Lombardozzi et al., 2012a, b).Treatment of this decoupled response in the YIBs model would lead to a higher level of O 3 flux entering leaves, thus causing stronger damage.Interestingly, this mechanism would provide a way to improve the simulated GPP overestimates.However, other studies have suggested that the O 3 damage to GPP may be offset by the benefits of co-located nitrogen deposition (Ollinger et al., 2002;Felzer et al., 2007), or even limited by carbon-nitrogen interactions (Kvalevag and Myhre, 2013).
Regional simulations for the US yield a summertime GPP (with high sensitivity O 3 damage) of 6.1 g C m −2 day −1 (9.5 g C m −2 day −1 in the eastern US and 3.9 g C m −2 day −1 in the western US).The total carbon uptake was estimated to be 4.43 ± 0.18 Pg C for the summer, accounting for 57-60 % of the annual value of 7.59 ± 0.25 Pg C over the 1998-2007 period.Carbon assimilation rate is suppressed by 4-8 % on average in the summertime eastern US with maximum local damage of 11-17 % in states close to the Great Lakes and along the eastern coast.When [O 3 ] is decreased by 25 %, O 3 damage to GPP is only 2-4 % in the eastern US, indicating substantial improvements to vegetation health and carbon assimilation rate.Previously, Felzer et al. (2004) found annual average O 3 -induced NPP reductions of 3-7 % over the US for 1989-1993 and simulated the largest reductions in states close to the Great Lakes and along the East Coast, where the high O 3 sensitivity of crops makes the dominant contribution.Our study examined O 3 damage effects a decade later than Felzer et al. (2004) but gives consistent results.Qualitatively, this consistency between decades may be explained by the offsetting influences of (i) surface O 3 reductions due to air quality control legislation and (ii) GPP increases due to CO 2 -fertilization and rising temperatures.Felzer et al. (2004) estimated a maximum local NPP reduction of 34 %, which is double the maximum of 17 % in our analyses.Furthermore, Felzer et al. (2004) found widespread reductions of > 6 % in the Midwest where there is almost no O 3 impact in this study (Fig. 8).Differences between the studies are mostly likely driven by the use of different vegetation cover and LAI data sets, and the use of a semi-mechanistic flux-based uptake in this study vs. the concentration-based uptake method elsewhere.
The current work has used an off-line approach.Yet, the O 3 -vegetation-meteorology system is strongly coupled.For instance, plant productivity itself controls the emission of isoprene, a major O 3 precursor.The O 3 -induced modification to stomatal conductance may inhibit evapotranspiration, leading to changes in canopy temperature, precipitation, soil moisture, and other surface hydrology and meteorology (Bernacchi et al., 2007;VanLoocke et al., 2012).In future work, we will study O 3 vegetation damage effects using YIBs embedded within a fully coupled global chemistryclimate model framework in order to account for these feedbacks including altered canopy energy fluxes and partitioning between latent and sensible heat that drive regional climate and hydrology.In addition, the O 3 damage algorithm parameters were calibrated using limited measurements for a few plant species, and were based on biomass yield not photosynthetic rate (Sitch et al., 2007).Future work will exploit recent extensive meta-data analyses (Lombardozzi et al., 2013;Wittig et al., 2007) to refine the ozone damage parameterization in YIBs including the decoupled modification of photosynthesis and stomatal conductance.a Site information is adopted from Schaefer et al. (2012), except that the operational time span listed here is only for the period when measurements of GPP are available.b PFT names are as follows: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), closed shrublands (CSH), mixed forests (MF), permanent wetlands (WET), and woody savannas (WSA).c The land type at US-ARM is cropland in Schaefer et al. (2012).However, the site is covered by cattle pasture and wheat fields (https://www.arm.gov/sites/sgp), which are more like C3 grassland.
The Supplement related to this article is available online at doi:10.5194/acp-14-9137-2014-supplement.

Figure 1 .
Figure 1.Comparison between monthly average GPP in the YIBs simulations and NACP observations grouped by PFT type (where in situ measurements are available).The red lines indicate linear regression between observations and simulation results.The regression fits and correlation coefficients are shown on each panel.The land types include evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), shrublands (SHR), grasslands (GRA), and croplands (CRO).The model-observation comparison for each site is shown in Fig. S4 in the Supplement.

Figure 2 .
Figure2.Histogram of (a) χ 2 for O 3 -free GPP and changes in χ 2 after the inclusion of O 3 damage impact with (b) low and (c) high O 3 sensitivity.Each bar represents the number of sites where the χ 2 or χ 2 of simulations fall between the specific thresholds as defined by the x axis intervals.The minimum and maximum of χ 2 or χ 2 are indicated as the two ends of x axis in the plots.The land cover definitions are as follows: ENF, Evergreen Needleleaf Forest; DBF, Deciduous Broadleaf Forest; SHR, Shrubland; GRA, Grasslands; CRO, Croplands.Results for each site are detailed in Fig.S6in the Supplement.

Figure 3 .
Figure 3.The calculated average χ 2 of GPP over NACP sites for six different simulations as listed in Table 2.The blue bars are results for all 40 NACP sites.The red bars are results excluding sites CA-Let, CA-NS1, US-Var, CA-SJ1, and CA-SJ2, where the simulated site-level χ 2 is larger than 16 as shown in Fig. S6a in the Supplement.

Figure 4 ..Figure 5 .Figure 6 .
Figure 4. Validation of simulated June-July-August (JJA) summertime average surface (a, b) diurnal mean and (c, d) daily maximum 8 h average O 3 with in situ measurements from (a, b) the EPA Clean Air Status and Trends Network (CASTNET) and (c, d) the AIRDATA.For (b) and (d), the blue points indicate sites east of 95 • W and the red points indicate sites west of 95 • W. The correlation coefficients are shown in (b) and (d).Each point in (b) and (d) represents the mean value for JJA at one specific site.Results for the model and observations are separated out in Fig. S7 in the Supplement.10941095 1096 1097 Simulated summertime (a) O 3 -exposed GPP and (b) O 3 stomatal flux over the U.S. simulated GPP is overlaid with in situ measurements from NACP.The simulations performed with land cover from ISLSCP and meteorological forcings from MERRA alysis.Fig. S8 separates results for model and measurements.

Figure 7 .
Figure 7. Simulated summertime (a) O 3 -exposed GPP and (b) O 3 stomatal flux over the US.The simulated GPP is overlaid with in situ measurements from NACP.The simulations are performed with land cover from ISLSCP and meteorological forcings from MERRA reanalysis.Figure S8 in the Supplement separates results for model and measurements.

Figure 8 .Figure 9 .
Figure 8. Simulated reduction fraction in summer GPP in the US due to (a) low and (b) high O 3 sensitivity for 1998-2007.
to assess the sensitivity of the results to local vs.reanalysis meteorological forcing and LAI (Table2).Two, MET-merra_LAImodis and METmerra_LAImerra, use hourly meteorology from MERRA-land reanalyses alone.The other two, METsite_LAImodis and METsite_LAImerra, use sitebased meteorology with gap-filled MERRA reanalysis.Simulations use two data sets of LAI: (1) METmerra_LAImerra and METsite_LAImerra use LAI from the MERRA-land reanalyses, which provide non PFT-specific LAI that we assign to the local PFT type at each site (TableA1), while (2) MET-merra_LAImodis and METsite_LAImodis use PFT-specific LAI retrieved by the MODIS.Later analyses show that MET-site_LAImerra has the lowest biases relative to other O 3free simulations.We perform two additional site-level simulations, which use the same forcings as that for MET-site_LAImerra but with the impact of O 3 uptake on photosynthesis.These two experiments, METsite_LAImerra_HO3 and METsite_LAImerra_LO3, use either high or low O 3 sensitivity as defined by the coefficient a in Table . Chem.Phys., 14, 9137-9153, 2014 www.atmos-chem-phys.net/14/9137/2014/ , 40, 60, 80, 100, 120, 140 ppbv, respectively.We do not include diurnal and seasonal variations of[O 3] in these sensitivity simulations as that in METsite_LAImerra for two reasons.First, field measurements for the O 3 vegetation damage are usually performed with fixed [O 3 ] during the growth season (e.g., Field and laboratory experiments may have different [O 3 ] compared to the ambient level we used complicating the validation.As a result, we perform 14 additional sensitivity simulations for each of NACP sites.All tests use meteorological and vegetation forcings the same as METsite_LAImerra (Table 2), except for the different [O 3 ] and O 3 sensitivity.These experiments are divided into two groups, seven in each, using either low or high O 3 sensitivity.In each group, simulations are performed with constant [O 3 ] at 20

Table 2 .
Description of the site-level simulations.The name of each simulation is composed of at least two words.The prefix indicates the source of meteorological forcings.The suffix or the second word indicates the sources of vegetation forcings.bForsimulationswithprefix "SITE" use site-based meteorological forcings, which are gap-filled with MERRA-land reanalyses.cAmbient[O3] is applied at each site.dLowand high indicate the sensitivity of GPP to [O 3 ] defined by the coefficient a in Table1. a

Atmos. Chem. Phys., 14, 9137-9153, 2014 www
Schaefer et al. (2012)hesis parameters that are not tuned to local site level vegetation parameters and, for instance, do not take into account plant species and age.Similar to the multimodel results inSchaefer et al. (2012), YIBs performance is weakest at the five grassland sites.In this case, the bias is mainly due to the delayed LAI seasonality in the MERRA satellite data set (Figs. S3 and S4 in the Supplement).In general, application of the remotely sensed LAI is a source of error because the gridded satellite data may not represent the local site changes in plant growth and phenology, especially for vegetation types with low biomass.The limitation of the satellite LAI spatial resolution implies that the model is unable to resolve GPP variability for sites in close proximity. .atmos-chem-phys.net/14/9137/2014/

Table 3 .
Statistics * for site-level simulations.Statistics include minimum and maximum values of R 2 , RMSE, and χ 2 for 35 NCAP sites with χ 2 < 16 (Fig.

Table A1 .
Descriptions of NACP sites in Canada (CA-) and US (US-) a .