Synthetic ozone deposition and stomatal uptake at flux tower sites

. We develop and evaluate a method to estimate O 3 deposition and stomatal O 3 uptake across networks of eddy covariance ﬂux tower sites where O 3 concentrations and O 3 ﬂuxes have not been measured. The method combines standard micrometeorological ﬂux measurements, which constrain O 3 deposition velocity and stomatal conductance, with a gridded dataset of observed surface O 3 concentrations. Measurement errors are propagated through all calculations to quantify O 3 ﬂux uncertainties. We evaluate the method at three sites with O 3 ﬂux measurements: Harvard Forest, Blodgett Forest, and Hyytiälä Forest. The method reproduces 83 % or more of the variability in daily stomatal uptake at


Introduction
Surface ozone (O 3 ) is toxic to both people and plants. Present-day and recent historical O 3 levels reduce carbon sequestration in the biosphere (Reich and Lassoie, 1984;Guidi et al., 2001;Sitch et al., 2007;Ainsworth et al., 2012), perturb the terrestrial water cycle (Lombardozzi et al., 2012(Lombardozzi et al., , 2015, and cause around $ 25 billion in annual crop losses (Reich and Amundson, 1985;Van Dingenen et al., 2009;Avnery et al., 2011;Tai et al., 2014). The basic plant responses to O 3 injury are well established from controlled exposure experiments (e.g., Wittig et al., 2009;Ainsworth et al., 2005Ainsworth et al., , 2012Hoshika et al., 2015), but few datasets are available to quantify O 3 fluxes and responses for whole ecosys-Published by Copernicus Publications on behalf of the European Geosciences Union.
tems or plant functional types that are represented within regional and global biosphere and climate models. The eddy covariance method has been widely used to measure landatmosphere fluxes of carbon, water, and energy and evaluate their representation in models (Baldocchi et al., 2001;Bonan et al., 2011), but few towers measure O 3 fluxes Fowler et al., 2001;Keronen et al., 2003;Gerosa et al., 2004;Lamaud et al., 2009;Fares et al., 2010;Stella et al., 2013;Zona et al., 2014). A recent review identified just 78 field measurements of O 3 fluxes over vegetation during the last 4 decades, many lasting just a few weeks (Silva and Heald, 2018). This paper demonstrates a reliable method to estimate O 3 fluxes at 103 eddy covariance flux towers spanning over 2 decades to enable O 3 impact studies on ecosystem scales.
The land surface is a terminal sink for atmospheric O 3 due to the reactivity of O 3 with unsaturated organic molecules and the modest solubility of O 3 in water. Surface deposition is 20 % of the total loss in tropospheric O 3 , making it an important control on air pollution (Wu et al., 2007;Young et al., 2013;Kavassalis and Murphy, 2017). This O 3 deposition flux includes stomatal uptake into leaves, where O 3 can cause internal oxidative damage, and less harmful non-stomatal deposition to plant cuticles, stems, bark, soil, and standing water (Fuhrer, 2000;Zhang et al., 2002;Ainsworth et al., 2012). O 3 can also react with biogenic volatile organic compounds, particularly terpenoid compounds, in the plant canopy air, and this process is commonly included in non-stomatal deposition (Kurpius and Goldstein, 2003). The deposition flux (mol O 3 m −2 s −1 ) can be described as where χ and χ 0 are the O 3 mole fractions (mol mol −1 ) in the atmosphere and at the surface, respectively, n is the molar density of air (mol m −3 ), and v d is a deposition velocity (m s −1 ) that expresses the net vertical O 3 transport between the height where χ is measured and the surface. F O 3 is defined positive for flux towards the ground. Equation (1) reasonably assumes that χ 0 = 0 because terrestrial surfaces have abundant organic compounds that react with and destroy O 3 . The deposition velocity can be decomposed into resistances (s m −1 ) for aerodynamic transport (r a ), diffusion in the quasilaminar layer (r b ), stomatal uptake (r s ), and non-stomatal deposition (r ns ) (Wesely, 1989): For stomatal and non-stomatal processes, the rates are often expressed as conductances (m s −1 ), which are the inverse of the resistances: g s = r −1 s and g ns = r −1 ns . The sum of stomatal and non-stomatal conductances is the vegetation canopy conductance, g c = g s + g ns . The stomatal O 3 flux is the portion of F O 3 that enters the stomata, and can be described as F s,O 3 = F O 3 g s (g s + g ns ) −1 = v d nχ g s (g s + g ns ) −1 .
To construct the synthetic O 3 flux, or SynFlux, we use measurements of O 3 concentration and standard eddy covariance flux measurements to derive nearly all of the terms in Eqs.
(1)-(3) from surface observations, using some additional information from remote sensing and models. This enables the estimation of F O 3 and F s,O 3 , as described in Sect. 2. Section 3 evaluates the method against observations at three sites that measure F O 3 and examines the importance of stomatal and non-stomatal deposition. Section 4 uses SynFlux to assess the spatial patterns of O 3 uptake to vegetation and to compare flux-based metrics of O 3 damage with concentration-based metrics. Finally, we discuss the strengths, limitations, and implications of our approach in Sect. 5.
2 Data sources and methods

SynFlux: synthetic O 3 flux
The FLUXNET2015 dataset (Pastorello et al., 2017) aggregates measurements of land-atmosphere fluxes of CO 2 , H 2 O, momentum, and heat at sites around the world (http://fluxnet. fluxdata.org/data/fluxnet2015-dataset; last access: 24 February 2017). Measurements are made with the eddy covariance method on towers above vegetation canopies (Baldocchi et al., 2001;Anderson et al., 1984;Goldstein et al., 2000) with consistent gap filling (Reichstein et al., 2005;Vuichard and Papale, 2015) and quality control across sites (Pastorello et al., 2014). Flux and meteorological quantities are reported in half-hour intervals. We analyze data from all sites in the United States and Europe in the FLUXNET2015 Tier 1 dataset. This analysis is restricted to the US and Europe because these regions have dense O 3 monitoring networks, described below. There are 103 sites meeting these criteria, all listed in Table S1 in the Supplement with references to full site descriptions. Three of these sites -Blodgett Forest, Harvard Forest, and Hyytiälä Forest -measure O 3 flux with the eddy covariance method, which we will use in Sect. 3 to evaluate our methods. SynFlux aims to constrain O 3 deposition and stomatal uptake as much as possible from measured water, heat, and momentum fluxes, in contrast to other methods (Finkelstein et al., 2000;Schwede et al., 2011;Yue et al., 2014) that rely more heavily on atmospheric models or parameterizations of stomatal conductance. From the eddy covariance measurements, we derive the resistance components of Eq. (2) using methods similar to past studies (Kurpius and Goldstein, 2003;Gerosa et al., 2005;Fares et al., 2010). The aerodynamic and quasi-laminar layer resistances (r a and r b , respectively) are derived from measured wind speed, friction velocity, and fluxes of sensible and latent heat every half hour using Monin-Obukhov similarity theory (Foken, 2017). The stomatal conductance for O 3 (g s ) is derived from the measured water vapor flux and meteorological data every half  teith, 1981;Gerosa et al., 2007). Supplement S1 provides further details of the resistance and conductance calculations. Some studies instead calculate g s from gross primary productivity (Lamaud et al., 2009;El-Madany et al., 2017), but that method is less widely used than the Penman-Monteith approach adopted here. The Penman-Monteith method of calculating stomatal conductance has been successfully applied across FLUXNET sites previously (Medlyn et al., 2017;Novick et al., 2016;Knauer et al., 2017;Lin et al., 2018). Those studies and others caution that, since evapotranspiration measurements include evaporation from ground, the stomatal conductance could be overestimated. While there are methods for quantifying and removing the evaporative fraction of evapotranspiration from eddy covariance data (Wang et al., 2014;Zhou et al., 2016;Scott and Biederman, 2017), a more common approach is to restrict analysis to conditions when transpiration dominates. We follow this second approach, analyzing only daytime data during the growing season, and use filtering criteria similar to Knauer et al. (2017). We define daytime as Sun elevation angle above 4 • and the growing season as days when gross primary productivity (GPP) exceeds 20 % of the annual maxima in GPP.
To avoid complications to the Penman-Monteith equation from wet canopies, we exclude times when dew may be present (RH > 80 %), and days with precipitation (> 5 mm). We also exclude the top and bottom 1 % of g s values, which include many unrealistic outliers (e.g., |g s | > 0.5 m s −1 ).  (Schnell et al., 2014). This dataset has 1 • spatial resolution, so some differences from measured O 3 abundances at individual sites are inevitable. Schnell et al. (2014) estimated these errors to be 6-9 ppb (rms) or about 15 % of summer mean O 3 in the US and similar in Europe. Figure 2 shows that the daytime gridded O 3 concentrations correlate well with observations at three flux tower sites where O 3 was measured (R 2 = 0.63-0.87) and have modest negative bias (5-10 ppb, −12 % to −28 %), consistent with the accuracy reported by Schnell et al. (2014). We use the Zhang et al. (2003) parameterization of non-stomatal conductance, which accounts for O 3 deposition to leaf cuticles and ground and was developed from measurements in the eastern United States. The parameterization requires leaf-area index, which we take from satellite remote sensing (Claverie et al., 2014(Claverie et al., , 2016, snow depth, which we take from MERRA2 reanalysis (GMAO, 2015;Gelaro et al., 2017), and standard meteorological data provided by FLUXNET2015. Uncertainties in these variables are described in Sect. 2.4. Performance of the non-stomatal parameterization is examined in Sect. 3.2. Figure 3 shows the stomatal O 3 flux at each site calculated with Eq. (3), and then averaged over the growing season. Figure S1 shows the corresponding total O 3 flux (Eq. 1). We refer to these products as the "synthetic" total O 3 flux (F

Observed O 3 flux
We evaluate SynFlux and its inputs at three sites where O 3 flux measurements are available: Harvard Forest, Massachusetts, United States ; Blodgett Forest, California, United States (Fares et al., 2010); and Hyytiälä Forest, Finland (Keronen et al., 2003;   Gridded and observed daily daytime O 3 concentrations at Blodgett, Harvard, and Hyytiälä forests. Inset numbers provide the coefficient of determination (R 2 ), mean and median bias, the standard major axis (SMA) slope, the Thiel-Sen (Sen) slope, and the 68 % confidence interval of the slopes. The black arrow points towards outliers that are not shown. Rannik et al., 2009). These forest sites sample a range of environmental and ecosystem conditions summarized in Table 1. All three sites have at least 6 years of halfhourly or hourly flux measurements. Two sites are evergreen needleleaf forests (Blodgett and Hyytiälä), while one is a deciduous broadleaf forest containing some scattered stands of evergreen needleleaf trees (Harvard). Climate also differs across these sites. Blodgett Forest has a Mediterranean climate with cool, wet winters and hot, dry summers. Hyytiälä and Harvard forests have cold winters and wetter summers, with Harvard Forest being the warmer of the two.
Harvard Forest water vapor flux measurements were recalibrated for this work based on matching water vapor mixing ratio measured by the flux sensor to levels calculated from ambient relative humidity and air temperature, resulting in a 30 % increase in evapotranspiration during the 1990s and no change since 2006. In addition, we remove sub-canopy evaporation from the measured water vapor flux before the Penman-Monteith calculation. Based on past measurements at these sites, the sub-canopy fraction of evapotranspiration is 20 % at Hyytiälä Forest and 10 % at Harvard Forest in summer Launiainen et al., 2005). We are unable to make this correction at all FLUXNET sites since water vapor flux is typically measured only above canopy.
At these three sites, observation-derived v d , g ns , and F s,O 3 can be derived from the F O 3 measurements with methods that differ slightly from Sect. 2.1. O 3 deposition velocity is inferred from measurements of O 3 concentration and flux via v d = F O 3 (nχ ) −1 . Resistance or conductance terms r a , r b , and g s are calculated as described in Sect. 2.1, and then both canopy and non-stomatal conductance are derived from ob- and g ns = g c − g s , respectively. With those values, Eq. (3) and use the same observation-derived g s , r a , and r b but different values of g ns , v d , and O 3 mole fraction.

Gap filling for friction velocity
The FLUXNET2015 dataset uses gap filling for most flux and meteorological measurements (Vuichard and Papale, 2015), but not for friction velocity (u * ), which is required to calculate v d and F syn s,O 3 . Filling this one variable would significantly reduce the fraction of missing data in our analysis. Monin-Obukhov similarity theory predicts that friction velocity will be proportional to wind speed in the surface layer, for a given roughness length and stability regime (Foken, 2017). On this basis, we regress the available friction velocity measurements against wind speed and net radiation (a proxy for stability) separately for each site and month (a proxy for vegetation roughness). This gap filling was possible at 91 sites that report net radiation measurements.
The predicted friction velocities from the regression model are correlated with available observations (R 2 > 0.5) and have minimal mean bias (±10 %) at 85 out of 91 eligible sites ( Fig. S3 in the Supplement), with most sites (63 out of 91) showing strong correlations (R 2 > 0.7). At the remaining six sites with lower regression model performance (R 2 < 0.5) we do not use u * gap filling. The u * gap filling increases the number of F syn s,O 3 estimates by 1 %-20 %. Time periods with u * gaps have no significant bias in meteorological conditions (e.g., mean wind speed, radiation, energy fluxes) compared to periods with u * measurements. As a result, the differences in monthly mean F syn s,O 3 with and without gap filling are small (10 % rms). So, although the u * gap filling is a potential source of uncertainty, the F syn s,O 3 estimates are robust. The following analysis will use the gap-filled data, but our results do not change in any meaningful way if we use the unfilled data.

Error analysis, averaging, and numerical methods
We quantify the errors in F syn O 3 , F syn s,O 3 , and all other calculated variables from the measurement uncertainties using standard techniques for propagation of errors through all equations (see Supplement S2). This method provides the uncertainty, quantified as the standard deviation, of each variable in each half-hour interval. The error analysis reveals that F syn s,O 3 and other derived quantities have uncertainties that change from hour to hour by 2 orders of magnitude (Fig. S2). In addition, many extreme values of F syn s,O 3 , g s , and other variables have very large uncertainties. We retain these outliers in our analysis and use the error analysis to appropriately reduce their influence on averages and other statistics, as described below, without discarding data.
The FLUXNET2015 dataset contains error estimates for sensible and latent heat measurements. We use these reported values in the error analysis. Where uncertainties in these fluxes are missing, we fill the gaps using a linear regression of available flux errors against flux values for that site. For friction velocity, the uncertainty is the prediction error in the linear model used for gap filling (Sect. 2.3). Based on expert judgment, the standard deviation of O 3 mole fraction is set to 20 %, pressure to 0.5 hPa, temperature to 0.5 K, relative humidity to 5 %, and canopy height to the lesser of 15 % or 2 m. For remotely sensed leaf-area index, the uncertainty is 1.1 m 2 m −2 for all vegetation types (Claverie et al., 2013(Claverie et al., , 2016. Snow depth uncertainty in MERRA2 is 0.08 m (Reichle et al., 2017). The Zhang et al. (2003) g ns parameterization has five vegetation-specific parameters and all are assigned 50 % standard deviation. Zero error is assumed for the flux tower height. Based on these inputs, the median relative uncertainty in F syn s,O 3 is 44 %, but it rises to several hundred percent for some half-hour intervals. The error analysis shows that most of the uncertainty in F syn s,O 3 derives from uncertainty in the latent heat flux measurement.
Daily and monthly averages of F syn s,O 3 and other quantities are constructed in stages. We first calculate a mean diurnal cycle for the day or month by pooling measurements during each hour in a maximum likelihood estimate, a weighted average that accounts for the uncertainty in each measurement. The maximum likelihood estimate is appropriate when combining values from the same distribution, which is expected to apply for measurements within a particular hour, but not across hours of the day. We then average across hours with an unweighted mean to calculate the daily or monthly value. For the daily averages, there are one to two observations within each hour. For the monthly averages, there are typically 30 to 60 in each hour of the day. We calculate seasonal averages with an unweighted mean of monthly values. Uncertainties are propagated through each stage of these averages, as detailed in Supplement S2. We compared averages with and without uncertainty weighting. The uncertainty-weighted averages tend to be smaller and less variable than unweighted averages because the error propagation identifies when outliers and large values have greater uncertainty. For example, the monthly values of g c derived from observations at Harvard Forest are 0.57 ± 0.11 cm s −1 with uncertainty weighting and 0.68 ± 0.17 cm s −1 without. Our discussion focuses on uncertainty-weighted daily averages of daytime data.

Data availability
The SynFlux dataset produced in this work is available at https://doi.org/10.5281/zenodo.1402054 (last access: 30 August 2018). The dataset includes synthetic stomatal and total O 3 fluxes, O 3 concentrations, O 3 deposition velocity, canopy conductance, stomatal conductance, and all of their propagated uncertainties. Monthly mean values are provided with and without u * gap filling, for 103 sites totaling 926 site years.   (Fig. S2), but the error analysis reveals that many of the outlying points have large uncertainties. For 98 % of points, the differences between F syn s,O 3 and F obs s,O 3 are less than the 95 % confidence interval derived from the error analysis (two-sided t-test). Thus, the errors in F syn s,O 3 are consistent with the propagated uncertainty in the observations. The half-hourly F syn s,O 3 values perform similarly well against observations (Fig. S4), but our analysis focuses on averages. The performance of daily F syn s,O 3 is partially due to resolving the seasonal cycle. If we subtract the mean seasonal cycle from both synthetic and observation-derived daily F s,O 3 , the residual correlation is R 2 = 0.5-0.7 (vs. 0.9 with the seasonal cycle included). This represents the skill of SynFlux at reproducing within-month and interannual variability. Overall, these results suggest that synthetic F syn s,O 3 is a reliable estimate of stomatal O 3 uptake into plants that can be used at flux tower sites without O 3 measurements.
The measurements also enable us to evaluate synthetic total deposition, F (3). The canopy resistance for O 3 is normally much greater than the quasi-laminar layer and aerodynamic resistances, meaning r c r a and r c r b , often by a factor of 3-10. Therefore, the O 3 deposition velocity is approximately v d ≈ r −1 c = g c . Under these conditions, Eq. (1) simplifies to F O 3 ≈ nχ (g s +g ns ) and Eq. (3) simplifies to F s,O 3 ≈ nχ g s . While g s is calculated from measured H 2 O fluxes, g ns comes from a parameterization, which inevitably introduces error into g ns and F syn  Figure 5 shows the seasonal cycles of observation-derived O 3 deposition velocity and its important components at the three study sites with O 3 flux measurements. For low or moderately reactive gases like O 3 , canopy resistance is typically greater than aerodynamic or quasi-laminar layer resistance, so it controls the overall deposition velocity. At these three sites, deposition velocity is lowest in winter (0.1-0.2 cm s −1 ) and highest in summer (0.5-0.6 cm s −1 ). Stomatal conductance peaks during warm and wet months, which explains most of this seasonal variation, except at Blodgett Forest as discussed below. Traditionally, stomatal conductance was thought to exceed non-stomatal conductance during the growing season at most vegetated sites (Wesely, 1989;Zhang et al., 2003), although this has been challenged more recently (Altimir et al., 2006;Stella et al., 2011a;Wolfe et al., 2011;Plake et al., 2015). At both Harvard and Hyytiälä forests, the mean stomatal conductance (0.2-0.6 cm s −1 ) is 1.5-6 times larger than non-stomatal conductance (0.08-0.2 cm s −1 ) during the growing season, so about 60 %-90 % of O 3 deposition occurs through stomatal uptake. At Blodgett, nonstomatal conductance slightly exceeds stomatal conductance in summer (0.4 vs. 0.3 cm s −1 ). The fast non-stomatal deposition is explained by O 3 reacting with biogenic terpenoid emissions below the flux measurement height (Kurpius and Goldstein, 2003;Fares et al., 2010). As documented in past work, these biogenic emissions depend strongly on temperature and light and have a large seasonal cycle with maxima in summer and minima in winter, so stomatal uptake is generally < 50 % of O 3 deposition at Blodgett in the summer but > 70 % in winter (Kurpuis and Goldstein, 2003;Fares et al., 2010;Wolfe et al. 2011).

Stomatal and non-stomatal deposition
A recent analysis of O 3 flux measurements at Harvard Forest suggests that non-stomatal deposition averages 40 % of daytime O 3 deposition during summer months, with a range of 20 %-60 % across years (Clifton et al., 2017). Our analysis of the same site does not support such a large role for non-stomatal deposition at this site in summer. For each year, we calculate summer daytime means of g s and g c by averaging the June-September values, and then calculate the nonstomatal fraction of deposition (1 − g s /g c ). Averaged across years 1993-2000, we find that 8 % of daytime O 3 deposition is non-stomatal during the summer, with a range of −33 % to 34 % across years. Negative fractions mean that stomatal conductance is large enough to explain all O 3 deposition. A large negative non-stomatal fraction (−33 %) occurs in only one year (1996) and no other year is less than −11 %, which is within uncertainty of 0 % (2σ ) according to the error propagation. Despite the small or zero non-stomatal fraction found here, our results continue to support the large year-to-year variability of this fraction reported by Clifton et al. (2017). The re-calibrated latent heat flux measurements are the main reason that our results differ from prior work and Supplement S3 provides further details. At Hyytiälä Forest, our results are consistent with prior work that found that the non-stomatal deposition is 26 % to 44 % of daytime O 3 deposition during the growing season (Rannik et al., 2012). Nevertheless, non-stomatal deposition equals or exceeds stomatal uptake where there are large terpene emissions (e.g., Blodgett) and at some other temperate sites that probably lack large biogenic emissions (Fowler et al., 2001;Cieslik, 2004;Lamaud et al., 2009;Stella et al., 2011b;El-Madany et al., 2017). We also examined interannual variation in O 3 deposition velocity. We find that the mean summer daytime At Blodgett Forest, the parameterized g ns is about half of observation-derived g ns in summer, but this is not surprising since the parameterization does not account for O 3 reactions with biogenic volatile organic compounds (BVOCs), which are known to be important at this site (Fares et al., 2010). In winter, however, the parameterized g ns values at Blodgett Forest are similar to observations (0.10 vs. 0.08 cm s −1 ). The parameterization is therefore able to roughly predict mean non-stomatal conductance in the absence of major BVOC emissions. Nevertheless, the parameterization reproduces almost none of the daily variability of g ns at any site (R 2 < 0.1, Fig. S7). This corroborates the recent field assessment that non-stomatal conductance is a weak point of most current dry deposition algorithms (Wu et al., 2018). We attempted, unsuccessfully, to use BVOC emissions from the MEGAN biogenic emission model (Guenther et al., 2012) to improve the g ns parameterization, but the correlations between compounds that react fastest with O 3 (monoterpenes and sesquiterpenes) and the observation-derived daily mean g ns were poor (R 2 ≤ 0.15). On that basis, F syn O 3 may also underestimate total O 3 deposition at other sites with high monoterpene and sesquiterpene emissions, such as warm-weather pine forests, but F Michigan, Nebraska, and Ohio) due to its moderate O 3 concentrations (Fig. S6) and moisture levels, which promotes stomatal conductance (Fig. 1) compared to surrounding sites due to irrigation and naturally wet soil in the California Delta. A combination of topography and climate is also an important factor in California: forest sites in the Sierra Nevada have lower g s and F syn s,O 3 than the lowland crops and wetland grasses. In Oregon, an evergreen needleleaf site regrowing after a fire has higher g s and F syn s,O 3 than two older forest stands nearby. The differences between nine Wisconsin forest sites, however, are mostly due to different years of data at each site combined with interannual variability in F syn s,O 3 ; fluxes at these sites are similar in overlapping years.
Variability across the 60 sites in Europe is controlled by similar factors. Stomatal uptake ranges from 1.4 to 9.6 nmol O 3 m −2 s −1 , with an average of 4.7 nmol O 3 m −2 s −1 (Fig. 3). The Mediterranean region has high O 3 concentrations ( Fig. S8) but generally low stomatal conductance due to the dry climate (Fig. 1). Within this region, vegetation type explains broad patterns. Shrub sites in Spain, France, and Sardinia have very low g s (∼ 0.15 cm s −1 ), so (3-6 nmol O 3 m −2 s −1 ), despite similar climate and O 3 . In central and northern Europe, temperate climate promotes higher stomatal conductance, while O 3 concentrations remain modest throughout the growing season. The largest F syn s,O 3 is 9.8 nmol O 3 m −2 s −1 at a deciduous broadleaf forest in Switzerland, while nearby evergreen forests, cereal crops, and grasslands all have lower fluxes (6-8 nmol O 3 m −2 s −1 ). While Finland has a generally low F syn s,O 3 of 2-5 nmol O 3 m −2 s −1 , the high end of this range is similar to rural sites in Germany, illustrating that O 3 can impact remote ecosystems with high stomatal conductance, even where O 3 concentrations are low. Table 2 quantifies SynFlux, O 3 deposition velocity, and conductance for each plant functional type. Wetlands, crops, and forests have the highest average F syn s,O 3 , which is about 2 times higher than woody savanna or shrublands, the vegetation types with the lowest F syn s,O 3 . At wetland sites, g s and F syn s,O 3 could be overestimated due to evaporation of surface water (Sect. 2.1), but any error is likely modest because our estimates of stomatal conductance at these sites (0.48±0.16 cm s −1 ; Table 2) are reasonable for wetland vegetation (up to 1 cm s −1 ; Drake et al., 2013). The vegetation types rank in the same order for stomatal conductance, again showing stomata as the main control on O 3 uptake into vegetation. Stomatal uptake exceeds non-stomatal uptake for all plant functional types except woody savanna and shrubland. O 3 deposition velocities reported in Table 2 fall within the ranges of past literature, as reviewed by Silva and Heald (2017). However, while Silva and Heald found that the mean deposition velocity was greater over deciduous forests than coniferous forests, crops, or grass, we do not. Rather, we find that variability between sites within each of these categories is large, having a standard deviation of about 30 % of the multi-site mean.

Metrics for O 3 damage to plants
Since O 3 injures plants mainly by internal oxidative damage after entering the leaves through stomata, the most physiological predictor of plant injuries is the cumulative uptake of O 3 (CUO, Reich, 1987;Fuhrer, 2000;Karlsson et al., 2004;Cieslik, 2004;Matyssek et al., 2007). CUO is defined as the cumulative stomatal O 3 flux exceeding a threshold flux Y that can be detoxified by the plant, integrated over a period of time: Here, H (x) is the Heaviside step function and t i is the time elapsed during measurement of F s,O 3 ,i . The sum is carried out over time i in the growing season, which we define based on GPP (Sect. 2.1). The detoxification threshold varies across vegetation types, even among related species (Karlsson et al., 2004;Büker et al., 2015), and thresholds for specific FLUXNET sites are generally unknown. As a compromise, we calculate CUO, with Y = 0, and also CUO3, with Y = 3 nmol O 3 m −2 s −1 , which has been suggested as a reasonable generic threshold . CUO is always greater than CUO3, but the sites with high CUO tend to also have high CUO3, so their spatial patterns are similar (Fig. S8).
While CUO is a physiological dose, concentration-based metrics remain common for assessing ozone impacts because they are easier to measure. Concentration-based metrics quantify O 3 in ambient air irrespective of whether that O 3 enters leaves. These metrics follow the general form where w(χ ) is a weighting function applied to the O 3 mole fraction χ, and χ c is a constant. Like CUO, the sum is usually over time i during the growing season. Three of the most common concentration-based O 3 metrics are the mean O 3 concentration, the accumulated concentration over a threshold of 40 ppb (AOT40; UNECE, 2004), and the sigmoidal-weighted index (W126; Lefohn and Runeckles, 1987). For mean, w (χ) = t i −1 and χ c = 0. For AOT40, w(χ ) = H (χ − χ c ) and χ c = 40 ppb. For W126, w (χ) = 1 + 4403 exp − 126 ppb −1 χ −1 and χ c = 0. Both AOT40 and W126 use only daytime (08:00-20:00) measurements and W126 also takes the maximum value over all 3-month periods during the growing season. The weighting functions for AOT40 and W126 give little or no weight to O 3 concentrations below 40 ppb. In addition, W126 gives increasing weight to concentrations up to about 110 ppb and full weight for higher concentrations based on the understanding that exposure to high O 3 concentrations is more injurious than moderate or low concentrations. Other concentration-based metrics (e.g., SUM60) use other thresholds or weighting functions, but many are strongly correlated with AOT40 or W126 or otherwise qualitatively similar (Paoletti et al., 2007). The spatial patterns of AOT40 and W126 closely resemble that of mean O 3 concentration in the US and Europe despite their different weighting functions (Fig. S9). AOT40 and W126 are well correlated with each other across sites (R 2 = 0.87) and with mean O 3 mole fraction (R 2 = 0.76 and R 2 = 0.52 for mean O 3 vs. AOT40 and W126, respectively) despite their different weighting functions. As a result, all of these concentration-based metrics have similar spatial patterns in the US and Europe. The CUO and CUO3 spatial patterns, however, are similar to F syn s,O 3 and distinct from the concentration-based metrics. This illustrates that locations with high AOT40 or W126, like the southwestern US or Mediterranean Europe, can have low CUO.
Even though concentration-based metrics do not measure the physiological O 3 dose to plants, they can be useful if the metric is proportional to the flux-based dose and injuries. Indeed, many controlled experiments and observational studies have documented correlations between both AOT40 and W126 and either uptake or plant injuries (e.g., Fuhrer et al., 1997;Cieslik, 2004;Musselman et al., 2006;Matyssek et al., 2010). However, many of these studies were carried out at a single site or under conditions where stomatal conductance was relatively steady while O 3 concentrations varied, for example by maintaining well-watered soil. When stomatal conductance varies widely, such as between arid and humid climates or seasons, concentration-based metrics may not correlate with stomatal O 3 flux . Figure 6 shows that all of the concentration-based metrics are poorly correlated with CUO across the sites (AOT40: R 2 = 0.05, W126: R 2 = 0.03, mean O 3 : R 2 = 0.04). Humidity helps explain some of the scatter in Fig. 6. The sites with high concentration-based metrics and low CUO have high vapor pressure deficit (VPD) and low stomatal conductance, and are mostly in the western US and Mediterranean Europe. Restricting the analysis to humid sites (VPD < 1.5 kPa) does not improve the correlation (R 2 ≈ 0.05) and at the arid sites (VPD > 1.6 kPa) the concentration-based metrics are modestly anti-correlated with CUO (AOT40: R 2 = 0.19, W126: R 2 = 0.05, mean O 3 : R 2 = 0.37). This result reinforces that concentration-based metrics can misrepresent CUO and plant injuries .
From the CUO values in Table 2, we can estimate the range of O 3 impacts on biomass production at the FLUXNET sites. Although species vary in their sensitivity to O 3 (Lombardozzi et al., 2013), several studies suggest that the biomass production of broadleaf and needleleaf trees decreases by 0.2 % to 1 % per mmol O 3 m −2 of CUO (Karlsson et al., 2004;Wittig et al., 2007;Hoshika et al., 2015). Combining the mean CUO for each plant functional type (Table 2) with these sensitivities, our work implies that O 3 reduces the biomass production at these FLUXNET sites by 6 %-29 % for deciduous broadleaf forests and 4 %-20 % for needleleaf forests. The range represents the spread of reported doseresponse sensitivities within each plant type, meaning the least and most O 3 -sensitive species. Several broadleaf crops are more sensitive to O 3 , with biomass reductions of 1.3 %-1.6 % per mmol O 3 m −2 of CUO3 . That sensitivity implies a 20 %-24 % drop in biomass production at FLUXNET crop sites. Some studies have quantified O 3 dose-response relationships with other thresholds Y = 1.6 to 6 nmol O 3 m −2 s −1 (e.g., Karlsson et al., 2007;Pleijel et al., 2004Pleijel et al., , 2014, but the sensitivities have a similar magnitude.  Fares et al. (2013) also demonstrated 12 %-19 % reduction in gross primary production due to O 3 at some of the same crop and forest FLUXNET sites. Using prognostic models of O 3 concentrations and stomatal uptake, several past studies have also suggested that O 3 reduces biomass production and CO 2 sequestration by 4 %-20 % in the US and Europe (Sitch et al., 2007;Wittig et al., 2007;Yue et al., 2014Yue et al., , 2016Lombardozzi et al., 2015). Our results support this range of impacts, although some FLUXNET sites and species likely experience greater O 3 injury, but here the CUO is highly constrained from observations and therefore avoids the additional uncertainties of atmosphere-biosphere models.

Conclusions
We have demonstrated a method to estimate O 3 fluxes and stomatal O 3 uptake at eddy covariance flux towers wherever regional O 3 monitors exist. , is sensitive to errors in the parameterized non-stomatal conductance, but mean values are still with a factor of 2 of observations. The errors in this dataset are modest compared with differences between observations and regional and global atmospheric chemistry models that are frequently a factor of 2 or more (Zhang et al., 2003;Hardacre et al., 2015;Clifton et al., 2017;Silva and Heald, 2017), illustrating the utility of this dataset for evaluating models and O 3 impacts.
Across flux tower sites in the US and Europe, F syn s,O 3 ranges from 0.5 to 11.0 nmol O 3 m −2 s −1 during the summer growing season. The spatial pattern of F syn s,O 3 is mainly controlled by stomatal conductance rather than O 3 concentration. Patterns of stomatal conductance and F syn s,O 3 in turn are explained by climate, especially atmospheric and soil moisture, vegetation types, and land management, such as irrigation. O 3 concentration-based metrics (AOT40, W126, mean O 3 ) have been widely used to evaluate O 3 damages to plants because they are easier and cheaper to measure than the cumulative uptake of O 3 (CUO) into leaves. However, these metrics have very little correlation with CUO (R 2 ≤ 0.05) across FLUXNET sites. Using dose-response relationships between CUO and biomass reduction, we estimate that O 3 reduces biomass production and carbon uptake by 4 %-29 %, depending on the site and plant type. Unlike most past estimates, which have used prognostic models of O 3 uptake, our assessment of biomass reduction is based on O 3 fluxes that are tightly constrained by observations. To promote further applications in ecosystem monitoring and modeling, the SynFlux dataset is publicly available as monthly averages of F syn s,O 3 , F syn O 3 , O 3 deposition velocity, stomatal conductance, and related variables.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This work was supported by the Winchester Fund and by the Council on Research Creativity at Florida State University. Eddy covariance data used here were acquired and shared by the FLUXNET community, including the AmeriFlux and CarboEuropeIP networks. The FLUXNET eddy covariance data processing and harmonization were 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 Center, and the OzFlux, ChinaFlux, and AsiaFlux offices. Trevor