Magnitude, trends, and impacts of ambient long-term ozone exposure in the United States from 2000 to 2015

Long-term exposure to ambient ozone (O₃) is associated with a variety of impacts, including adverse human-health effects and reduced yields in commercial crops. Ground-level O₃ concentrations for assessments are typically predicted using chemical transport models; however such methods often feature biases that can influence impact estimates. Here, we develop and apply artificial neural networks to empirically model long-term O₃ exposure over the continental United States from 2000 to 2015, and we generate a measurement-based assessment of impacts on human-health and crop yields. Notably, we found that two commonly used human-health averaging metrics, based on separate epidemiological studies, differ in their trends over the study period. The population-weighted, April–September average of the daily 1 h maximum concentration peaked in 2002 at 55.9 ppb and decreased by 0.43 [95 % CI: 0.28, 0.57] ppb yr⁻¹ between 2000 and 2015, yielding an ∼18 % decrease in normalized human-health impacts. In contrast, there was little change in the population-weighted, annual average of the maximum daily 8 h average concentration between 2000 and 2015, which resulted in a ∼5 % increase in normalized human-health impacts. In both cases, an aging population structure played a substantial role in modulating these trends. Trends of all agriculture-weighted crop-loss metrics indicated yield improvements, with reductions in the estimated national relative yield loss ranging from 1.7 % to 1.9 % for maize, 5.1 % to 7.1 % for soybeans, and 2.7 % for wheat. Overall, these results provide a measurement-based estimate of long-term O₃ exposure over the United States, quantify the historical trends of such exposure, and illustrate how different conclusions regarding historical impacts can be made through the use of varying metrics.

1 Introduction

Tropospheric ozone (O₃) is a secondary pollutant that is photochemically formed from precursor gases. Exposure to ambient O₃ is associated with adverse health effects in humans (U.S. EPA, 2013) and reduced yields in commercial crops (Chameides et al., 1994; Mauzerall and Wang, 2001). These impacts have driven efforts to reduce ground-level O₃ in the United States, specifically targeting peak levels of O₃ concentrations through regulations that control anthropogenic precursor emissions, such as nitrogen oxides (NO_x) and volatile organic compounds (VOCs). O₃ reduction efforts have been widely successful in reducing peak concentrations (Simon et al., 2015; Lefohn et al., 2017; Fleming et al., 2018), but impacts related to both human health and crop yields nonetheless persist (Cohen et al., 2017; Seltzer et al., 2018; Zhang et al., 2018; Shindell et al., 2019).

Quantifying impacts requires an estimate of exposure to O₃, which is most commonly accomplished through the use of chemical transport models (CTMs; e.g., Anenberg et al. 2010; Silva et al., 2013; Lelieveld et al., 2015; Malley et al., 2017; Shindell et al., 2018; Stanaway et al., 2018). CTMs apply state-of-the-science knowledge to simulate O₃ formation, termination, and transport, while also providing complete spatial and temporal coverage over a particular domain – a desired trait for impact assessments. However, estimates of exposure and impacts can vary substantially across CTM studies. For example, two CTM-based studies estimated 2005 respiratory-related premature mortalities in the USA using the same relative risk function (Jerrett et al., 2009), yet yielded results that differed by $\sim \mathrm{3}\times$ (i.e., 13 000 vs. 38 000; Zhang et al., 2018; Lelieveld et al., 2013). While CTMs accurately reproduce many features of atmospheric chemistry (Shindell et al., 2013; Hu et al., 2018), one important issue associated with CTM-based impact assessments is that CTMs are consistently biased high when predicting O₃ concentrations (e.g., Schnell et al., 2015; Travis et al., 2016; Yan et al., 2016; Seltzer et al., 2017, Porter et al., 2017; Guo et al., 2018). Such biases can influence estimates of impacts and are often amplified by nonlinear concentration–response functions (Seltzer et al., 2018). Measurement-based methods, including area-weighted average of nearby monitors, nearest monitor, inverse distance weighting, Kriging interpolation, and multiple-linear regression under a Bayesian framework, can also be used to estimate exposure (Bell, 2006; Brauer et al., 2008; Marshall et al., 2008; Chang et al., 2010; Seltzer et al., 2018). However, a notable limitation of such methods stems from the sparse spatial coverage of monitoring sites. While these limitations might be minor in areas with dense monitoring, such methods can become insufficient as the distance from monitors increases (Bell, 2006).

O₃ exposure trends are also of great interest to researchers and air quality managers. To accurately model trends of O₃ exposure, many dimensions of variability must be captured. For the annual average of the maximum daily 8 h average O₃ concentration (hereafter MDA8), a metric that has been used to quantify cause-specific long-term O₃ exposure associations in epidemiological studies (e.g., Turner et al., 2016; Lim et al., 2019), the O₃ diurnal and seasonal cycles must be accurately simulated over time. CTM evaluation studies also report the existence of seasonal, spatial, and diurnal variability in model performance (Cooper et al., 2014; Schnell et al., 2015; Seltzer et al., 2017; Lin et al., 2017; Guo et al., 2018; Strode et al., 2019; Young et al., 2018), which can lead to conflicting conclusions regarding trends in exposure. For example, Zhang et al. (2018) report a ∼9 % decrease in the population-weighted, daily maximum 1 h exposure concentration of O₃ in the US warm season between 1990 and 2010. Meanwhile, a separate study reported no change in the population-weighted, daily maximum 8 h exposure concentration of warm-season O₃ over those same 2 decades (Stanaway et al., 2018). Since monitoring data are sparse, quantification of trends using observations requires either continuous, long-term measurement data at a particular site or the aggregation of observations into regions (e.g., Southeast, Northeast, Great Plains) and/or urban–rural–suburban classifications. Many studies have indeed made use of such data to assess O₃ trends (Jaffe and Ray, 2007; Cooper et al., 2012, 2014; Parrish et al., 2012; Simon et al., 2015). The recent publication of the Tropospheric Ozone Assessment Report (TOAR) database (Schultz et al., 2017) has created a rich observational dataset and further expanded the number of such assessments (e.g., Chang et al., 2017; Gaudel et al., 2018; Lefohn et al., 2018; Fleming et al., 2018; Mills et al., 2018b).

In this study, we applied artificial neural networks (ANNs) and the TOAR database to estimate a suite of O₃ impact metrics related to human health and crop yield over the contiguous United States from 2000 to 2015 at $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ resolution. Specifically, we took advantage of the improved long-term coverage afforded by the TOAR database to develop a framework that empirically estimates O₃ exposure with complete spatial and temporal coverage over the United States. ANNs have been previously used to make O₃ predictions (Ruiz-Suárez et al., 1995; Yi and Prybutok, 1996; Comrie, 1997; Gardner and Dorling, 2000; Dutot et al., 2007; Di et al., 2017), but generally at the monitor or city level. Our main goal was to better quantify the magnitude and trends of population-weighted and agriculture-weighted long-term (i.e., months, annual) O₃ exposure in the USA over many consecutive years and use those estimates to generate a measurement-based assessment of impacts and trends on human health and crop yields. In addition, we tested and applied the ANN to meteorologically adjust exposure predictions, thus eliminating a substantial proportion of the short-term variability and enabling a separate quantification of long-term O₃ exposure trends.

2 Methods

2.1 Observational dataset and impact metrics

Daily O₃ observations spanning 2000–2015 from the University of New Hampshire Air Quality and Climate Program (Airmap), the U.S. Air Quality System (AQS), the Canadian Air and Precipitation Monitoring Network (CAPMoN), the U.S. Clean Air Status and Trends Network (CASTNET), the Global Atmosphere Watch (GAW), and the Canada National Air Pollution Surveillance (NAPS) monitoring networks in North America were retrieved from the Tropospheric Ozone Assessment Report (TOAR) database (Schultz et al., 2017). The reader is referred to Schultz et al. (2017) for a detailed description of these networks, including variations in network area type (i.e., urban vs. suburban vs. rural) and number of monitors. These daily observations were used to calculate two human-health- and two crop-yield-relevant averaging metrics. The first human-health metric is from the Jerrett et al. (2009), hereafter J2009, long-term O₃ exposure epidemiology study. Using data from the American Cancer Society Cancer Prevention Study II (ACS CPS-II) cohort, J2009 estimated changes in cause-specific mortalities attributable to incremental changes in the April–September average of the daily 1 h maximum O₃ concentration (hereafter MDA1). The second human-health metric is from the Turner et al. (2016), hereafter T2016, long-term O₃ exposure epidemiology study. T2016, using an expanded version of the ACS CPS-II cohort that included more follow-up years, a larger population, and more events (i.e., deaths), reported changes in cause-specific mortalities attributable to incremental changes in the annual average of the maximum daily 8 h average O₃ concentration (hereafter MDA8). To elucidate the influence of the underlying seasonal trends on the MDA8 metric, we also subdivided this annual metric into 3-month seasonal windows (i.e., summer: June–August; spring: March–May). These seasonal divisions feature the following labels: MDA8-MAM (spring), MDA8-JJA (summer), MDA8-SON (fall), and MDA8-DJF (winter).

The two crop-loss metrics included here were the M12 (12 h mean) and AOT40 (accumulated amount of O₃ over the 40 ppb threshold) averaging metrics. Both have been used in a variety of crop loss assessments (e.g., Van Dingenen et al., 2009; Avnery et al., 2011; Shindell et al., 2019). The M12 metric, which can be used to calculate impacts on maize and soybean relative yields, is defined as the mean O₃ value for the local hours of 08:00–20:00, averaged over the 3 months prior to the start of the harvest period. The AOT40 metric, which can be used to calculate impacts on maize, soybeans, and wheat, is an accumulative index and defined as a summation of the hourly mean O₃ values over 40 ppb for the local hours of 08:00–20:00, also averaged over the 3 months prior to the start of the harvest period. We initialized the start of the harvest period to be consistent with Avnery et al. (2011). For maize and soybeans, the 3-month averaging period was initialized in July. Wheat features two varieties with separate initialization months for harvesting. One is initialized in March and the other is initialized in May. Exposure results of both varieties are included for illustrative and seasonal comparisons. It should be noted that long-term O₃ exposure also stunts the yields of a variety of other crops, such as rice (Mills et al., 2007; Van Dingenen et al., 2009; Shindell et al., 2019), but inclusion of these impacts was not considered here since they are not major commercial crops in the United States.

2.2 Artificial neural network

We utilized feed-forward artificial neural networks (ANNs), which are also referred to as multilayer perceptrons, to model the four metrics considered here, with a unique network for each metric. ANNs were constructed using the Keras API (https://keras.io, last access: 5 February 2020; Chollet, 2015) and TensorFlow machine-learning library (https://www.tensorflow.org, last access: 5 February 2020; Abadi et al., 2015). Broadly, ANNs consist of several interconnected layers, beginning with an input data layer, ending with an output data layer, and having at least one “hidden” layer between the input and output that models the nonlinear relationships of the system. Each layer is connected via a set of coefficients at individual “nodes” that are optimized through model training, similar to a multiple linear regression (MLR) model. In contrast to a MLR, a layer in an ANN may have multiple nodes, and the output from each node proceeds through an “activation function”. An ANN activation function can take many shapes, but the two most common are a sigmoidal function (which converts the node output to a probability) and a rectified linear (ReLu) function (which applies a threshold to a linear function). The ANNs used here consisted of one input, three hidden, and one output layer. All nodes in each hidden layer featured a ReLu activation function, including the output layer to ensure all predictions were non-negative. The three hidden layers, each of which included a bias term, consisted of 32 nodes each. This particular architecture was selected following the testing of various configurations (i.e., differences in the number of nodes and layers), with the ultimate goal to prevent over-fitting of model parameters and maximizing model generalization (see Sect. 3.1 for added discussion).

Table 1 Variable input parameters (i.e., non-fixed effect parameters) for each artificial neural network (ANN). All 24 h and 12 h (08:00–20:00) periods were adjusted to local times.

Daily observations from the TOAR dataset spanning 2000–2015 were paired with MERRA-2 meteorological reanalysis data (Gelaro et al., 2017), anthropogenic emissions data from the Community Emissions Data System (CEDS) inventory (Hoesly et al., 2018), monthly East Asian anthropogenic emissions (Hoesly et al., 2018, Zheng et al., 2018), and monthly methane concentrations (GLOBALVIEW-CH4, 2009). Details regarding these parameters are provided in Table 1. Meteorological variables that were considered O₃ covariates largely follow Li et al. (2019). Local anthropogenic emissions (Hoesly et al., 2018) included nitrogen oxides (NO_x), non-methane volatile organic carbon (NMVOC; includes total weight of all species), and carbon monoxide (CO). Since emissions from East Asia have a large impact on North American ground-level O₃ concentrations (Liang et al., 2018) and have dramatically changed in recent decades (Zheng et al., 2018), monthly total emissions from all East Asian countries (i.e., China, Japan, South Korea, North Korea, and Mongolia) were included as an input. Emissions from all East Asian countries were retrieved from the CEDS inventory, with the exception of Chinese emissions, which were retrieved from the Multi-resolution Emission Inventory for China (MEIC) inventory (Zheng et al., 2018). As the last year included in the CEDS inventory is 2014, anthropogenic emissions for 2014 were repeated for 2015. To incorporate geographical differences and long-term drivers not included as input, several fixed-effect parameters were also used as input, including latitude, longitude, and year. A single input into each ANN consists of all the variables described above, which are paired in space and time to an observation retrieved from the TOAR database. Finally, all input data were normalized by subtracting the mean and dividing by the standard deviation of the training dataset.

Prior to model training, the complete dataset was divided into three components – training, validation, and testing. The training dataset is used to iteratively tune the coefficients in the ANN, the validation dataset is used to ensure the training process does not over-fit the ANN parameters to match the training dataset, and the testing dataset is used to evaluate how well the trained model performs. To compile these components, all available data in a given month were collected and 3 random, consecutive days were removed for validation, and 4 random, consecutive days were removed for testing. The remaining days became the training dataset. Overall, the size of the training dataset (i.e., the number of compiled inputs) eclipsed 5 million values; therefore, the number of trainable parameters was nearly 4 orders of magnitude smaller. The optimization of all coefficients at each node in the ANN is accomplished through stochastic gradient descent (SGD) optimization. SGD consists of (a) taking mini-batches of the training dataset, (b) estimating the gradient of all coefficients relative to the known output, (c) taking a small iterative step towards an optimal solution, (d) repeating with a new mini-batch of the training dataset, and (e) repeating steps (a)–(d) until the entire training dataset has been fed through the network. Proceeding through steps (a)-(e) is referred to as an epoch, and network training proceeds through multiple epochs. In total, we used the Adam optimizer (Kingma and Ba, 2015), with a learning rate of 0.001 (i.e., the size of step (c)), a decay factor of 0.9 (i.e., a shrinking of the step (c) size), and a mean-squared error target cost-function. Each ANN was trained for 3000 epochs, with a shuffling of the training data between each epoch. Through monitoring of the model training using the validation dataset, it was determined that 3000 epochs were sufficient to optimize the system without over-fitting.

To quantify the added benefit of the ANN over a simplified model, a comparison with results from a MLR is included. In addition, since exposure mostly occurs at unobserved locations, and all of the model training explained thus far is only evaluated at observed locations (i.e., from the TOAR database), we added an additional step to test our methods. In short, we performed several CTM simulations and sampled the daily-level CTM predictions of each metric at all available monitoring locations, generating what we refer to as a “pseudo-observational dataset”. We then followed the same machine learning process described above, except using the pseudo-observational dataset and four newly trained ANNs, to predict the population-weighted (MDA1/MDA8) and agriculture-weighted (M12/AOT40) exposure values estimated by the CTM. Through this process, we can assess the network's ability to predict total exposure through the exclusive use of sparse measurements.

2.3 Chemical-transport modeling

GEOS-Chem was used to generate ground-level CTM predictions of O₃ (v11-01; http://www.geos-chem.org, last access: 5 February 2020; Bey et al., 2001). A nested version of the model at $\mathrm{0.5}{}^{\circ }\times \mathrm{0.625}{}^{\circ }$ horizontal resolution, driven by native-resolution MERRA-2 meteorology and fed varying $\mathrm{2.0}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ boundary conditions, was utilized to simulate O₃ throughout the continental United States for the years 2000, 2003, 2005, 2007, 2010, 2012, and 2014. The model includes comprehensive ${\mathrm{HO}}_x-{\mathrm{NO}}_x-\mathrm{VOC}-{\mathrm{O}}_x$ gas chemistry, coupled to an aerosol module that includes sulfate–nitrate–ammonium chemistry (Park et al., 2004; Pye et al., 2009), primary carbonaceous aerosols (Park et al., 2003), mineral dust (Fairlie et al., 2007), and sea salt (Jaegle et al., 2011), with aerosol thermodynamics simulated using ISORROPIA II (Fountoukis and Nenes, 2007). Global anthropogenic emissions come from the CEDS inventory (Hoesly et al., 2018) and were processed through the Harvard–NASA Emission Component (HEMCO; Keller et al., 2014). All nested simulations featured a 2-month spin-up, each O₃ metric was calculated at local time, and ground-level concentrations (10 m) from the first level of the GEOS-Chem output were calculated using the methods described in Sect. 3 of Zhang et al. (2012).

2.4 Calculation of metric trends

Trends of all metrics are presented both spatially and weighted towards the subject of interest (i.e., population-weighted or agriculture-weighted). Since the metrics considered here are based on long-term exposure, all trends were assessed at the annual timescale (i.e., one data point, either grid cell or a population-/agriculture-weighted value) using a linear least-squares regression. To calculate population-weighted exposure concentrations, we used population data from the 2017 revisions to the UN Population Division (United Nations, 2017), distributed to grid cells using population density data from the Gridded Population of the World (GPW) version 4 (CIESIN, 2016). Agriculture-weighted exposure concentrations were calculated using crop production and density data from the Food and Agricultural Organization datasets (FAO, 2010).

We also accounted for short-term variability in metric trends by modeling meteorologically adjusted predictions of each metric. To evaluate the ANNs' ability to complete this task, we performed CTM simulations of 2003, 2005, 2007, 2010, 2012, and 2014 using meteorological conditions from each respective year, but frozen anthropogenic emissions and methane concentrations from 2000. We then used the previously trained ANNs (i.e., the ANNs generated using the pseudo-observational data) to predict the population-weighted and agriculture-weighted exposure metrics from these CTM sensitivity simulations. Finally, we compared the CTM- and ANN-predicted trends attributed to meteorology between 2000 and 2015. This enabled us to evaluate how well the ANN can meteorologically adjust exposure trends. From there, the same methods were applied using the ANNs trained with the TOAR data to estimate meteorologically adjusted trends of the population-weighted and agriculture-weighted exposure metrics. Specifically, all variables were held frozen at 2000 values, except for the MERRA-2 meteorological conditions.

2.5 Calculation of human-health and crop-yield impacts

Human-health impacts were quantified using the exposure–response relationships and averaging metrics reported by J2009 and T2016. Both epidemiological studies found a significant relationship between exposure to long-term O₃ and premature respiratory mortality. Respiratory impacts are the lone end point considered here since they are the most common impact reported by the community. However, T2016 and several other studies (Jerrett et al., 2013; Crouse et al., 2015; Cakmak et al., 2016; Lim et al., 2019) reported a significant relationship between long-term O₃ exposure and other mortality end points, such as cardiovascular disease.

Impact assessments for human health generally report results as the estimated number of premature mortalities attributable to long-term exposure. However, these results can often be driven by non-exposure variables, such as changes in population count (Cohen et al., 2017), baseline mortality rates (Cohen et al., 2017), and population aging (Apte et al., 2018). To eliminate the influence of changes in the total population count on net impacts, we normalized our results and report estimated health impacts as premature mortalities per 100 000 people attributable to long-term O₃ exposure. We also illustrate the percent contributions of each variable (i.e., population aging, changes in baseline mortality rates, and exposure) on the net health impact calculations.

Normalized premature mortalities attributable to long-term O₃ exposure were calculated as follows.

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }X=\left\{\begin{array}{ll}\mathrm{0}& \text{if}\left[{\mathrm{O}}_{\mathrm{3}}\right]\le \mathrm{TMREL}\\ \left[{\mathrm{O}}_{\mathrm{3}}\right]-\mathrm{TMREL}& \text{if}\left[{\mathrm{O}}_{\mathrm{3}}\right]>\mathrm{TMREL}\end{array}\right.\text{(2)}& {\displaystyle }\mathrm{HR}={\exp }^{\mathit{\beta }\mathrm{\Delta }Y}\text{(3)}& {\displaystyle }\mathrm{AF}=\mathrm{1}-{\exp }^{-\mathit{\beta }\mathrm{\Delta }X}\text{(4)}& {\displaystyle }\mathrm{\Delta }{\mathrm{Mort}}_i=y_{\mathrm{0}i}\times \mathrm{AF}\times {\mathrm{population}}_i\text{(5)}& {\displaystyle }\begin{array}{rl}\mathrm{Normalized}\mathrm{Mort}& =\left({\sum }_{i=\mathrm{1}}^n\mathrm{\Delta }{\mathrm{Mort}}_i/{\sum }_{i=\mathrm{1}}^n{\mathrm{population}}_i\right)\\ & \times \mathrm{100}\mathrm{000}\end{array}\end{array}$$

Here TMREL is the theoretical minimum risk exposure level, ΔX is the predicted long-term O₃ exposure concentration above the TMREL, β is the exposure–response factor, HR is the hazard ratio reported by the epidemiological study, ΔY is 10 ppb in both epidemiological studies, AF is the attributable fraction of the disease burden attributable to long-term O₃ exposure, y₀ is the cause-specific, age-binned baseline mortality rate, population is the age-binned population count, i is the age bin index, ΔMort is the estimated number of cause-specific, age-binned premature mortalities, n is the number of age bins, and Normalized Mort is the estimated number of cause-specific premature mortalities per 100 000 people attributable to long-term O₃ exposure. Baseline mortality rates were retrieved from the 2017 GBD (Global Burden of Disease) project (Stanaway et al., 2018) and mapped to best match the ICD-10 (International Statistical Classification of Diseases and Related Health Problems) codes reported in T2016. The hazard ratio for respiratory diseases was 1.040 (95 % CI: 1.013, 1.067) and 1.12 (95 % CI: 1.08, 1.16) in J2009 and T2016, respectively. The TMRELs used were 33.3 ppb when using the J2009 averaging metric and 26.7 ppb when using the T2016 averaging metric, as reported by each epidemiological study.

We report agriculture (maize, soybean, and wheat) impacts in terms of a national relative yield loss (RYL) due to long-term O₃ exposure. We utilized the concentration–response function and RYL methods outlined in Van Dingenen et al. (2009), as summarized below.

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}\mathrm{Maize}\mathrm{RYL}\left[M\mathrm{12}\right]& =\mathrm{1}-\left(\exp \left[-{\left({\displaystyle \frac{M\mathrm{12}}{\mathrm{124}}}\right)}^{\mathrm{2.83}}\right]/\right.\\ & \left.\exp \left[-{\left({\displaystyle \frac{\mathrm{20}}{\mathrm{124}}}\right)}^{\mathrm{2.83}}\right]\right)\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}\mathrm{Soybean}\mathrm{RYL}\left[M\mathrm{12}\right]& =\mathrm{1}-\left(\exp \left[-{\left({\displaystyle \frac{M\mathrm{12}}{\mathrm{107}}}\right)}^{\mathrm{1.58}}\right]/\right.\\ & \left.\exp \left[-{\left({\displaystyle \frac{\mathrm{20}}{\mathrm{107}}}\right)}^{\mathrm{1.58}}\right]\right)\end{array}\text{(8)}& {\displaystyle }\mathrm{Maize}\mathrm{RYL}\left[\mathrm{AOT}\mathrm{40}\right]=\mathrm{AOT}\mathrm{40}\times \mathrm{0.00356}\text{(9)}& {\displaystyle }\mathrm{Soybean}\mathrm{RYL}\left[\mathrm{AOT}\mathrm{40}\right]=\mathrm{AOT}\mathrm{40}\times \mathrm{0.0113}\text{(10)}& {\displaystyle }\mathrm{Wheat}\mathrm{RYL}\left[\mathrm{AOT}\mathrm{40}\right]=\mathrm{AOT}\mathrm{40}\times \mathrm{0.0163}\end{array}$$

Table 2 Daily-level training, validation, and testing performance metrics (RMSE) of the ANN using GEOS-Chem sampled data at TOAR locations compared to a multiple linear regression model. Note ppbh denotes parts per billion-hour.

3 Results

3.1 Artificial neural network training and evaluation

We began our evaluation using the pseudo-observational dataset derived from daily GEOS-Chem output sampled at all available monitoring locations. From there, we used this dataset to train four ANNs (i.e., one for each metric) and attempted to recreate the original GEOS-Chem output. Through this process, we attempted to determine the strength of an ANN in reconstructing complete exposure maps using sparse observation data. The RMSE results from the initial training and validation datasets were similar (Table 2), indicating that the network was not over-fitting and generalized the system well. When compared to a MLR, the RMSE testing results were ∼33 % lower, demonstrating the added benefit of the ANN (Table 2). Population-weighted and agriculture-weighted exposure estimates from the ANN closely matched the predictions from GEOS-Chem (red vs. blue in Fig. 1) for all metrics and featured a high coefficient of correlation (insets in Fig. 1). An exception was the marginal high bias of the AOT40 metrics for wheat early in the time series. These small deviations were due to a few factors. First, regions of dense agriculture production are limited and generally located in areas with fewer monitors, thus limiting the extent of model training. Second, the AOT40 metric is an accumulation index, which can lead to the amplification of small biases. Separately, we also found the ANN to perform well when meteorologically adjusting the predicted exposure trends (i.e., the short-term trends attributable to meteorology; see insets and green vs. yellow lines in Fig. 1). In total, the ANN was able to reproduce the complete exposure predictions with high fidelity, as estimated by GEOS-Chem, using information strictly from monitoring locations.

Figure 1 Red – GEOS-Chem (GC) simulated values of all metrics. Blue – ANN predictions of GC using daily samples from the GEOS-Chem simulations at all available monitoring locations. Green – meteorological trend simulated by GEOS-Chem with all input frozen at 2000 levels, with the exception of meteorological variables. Yellow – ANN prediction of GC-MET (green) using the previously trained (i.e., blue) neural networks. Inset within each panel is the coefficient of correlation for the GC–ANN and GC-MET–ANN-MET time series.

Table 3 Daily-level training, validation, and testing performance metrics (RMSE) of the ANN using TOAR observations compared to a multiple linear regression model. Note that ppbh denotes parts per billion-hour.

With confidence in the overall framework, we then trained new ANNs using daily 2000–2015 observations from the TOAR database. Little difference between the training, validation, and testing performance metrics indicated that each ANN was not over-fitting to the training dataset (Table 3). In addition, we again found the ANN to perform ∼30 % better than a MLR model (Table 3). When compared to the original TOAR database, we found high accuracy between each ANN-predicted long-term metric and the original observations (Figs. S1–S4; Table S1 in the Supplement). The RMSE of the MDA1 and MDA8 predictions ranged from 3.1 to 4.4 and 2.3 to 3.9 ppb, respectively (Table S1). The r² of the two metrics ranged from 0.77 to 0.84 and 0.74 to 0.82. Similar levels of bias (RMSE) and correlation (r²) were found when comparing the long-term agriculture metrics (Table S1).

Figure 2 Trends of the MDA1 (a; ppb yr⁻¹) and MDA8 (b; ppb yr⁻¹) health metrics from 2000 to 2015. Trends of MDA8-JJA (c; ppb yr⁻¹) and MDA8-DJF (d; ppb yr⁻¹) from 2000 to 2015. The p values from these trends can be found in Fig. S5.

3.2 Magnitude and trends of long-term O₃ exposure metrics

3.2.1 Human-health-relevant metrics

The MDA1 metric featured large reductions throughout the study period, with downward trends exceeding 1 ppb yr⁻¹ in the Southeast and in portions of California (Fig. 2). As a result, exposure throughout this period simultaneously decreased. The national population-weighted exposure concentration peaked in 2002 at 55.9 ppb, reached a minimum of 48.2 ppb in 2014, and featured sizable year-to-year fluctuations due to inter-annual variation (Fig. 3). From 2000 to 2015, the national population-weighted exposure concentration of the MDA1 metric featured an annual decrease of 0.43 [95 % CI: 0.28, 0.57] ppb yr⁻¹ (Table S4). After adjusting for meteorology, the trend changed to −0.41 [95 % CI: −0.35, −0.47] ppb yr⁻¹. The similar mean values of these two trends suggest that nearly all of the MDA1 reductions are due to non-meteorological drivers (i.e., emission changes, intercontinental transport, methane). Changes in exposure featured an east–west divide, with population-weighted exposure concentrations decreasing by 0.49 [95 % CI: 0.28, 0.69] ppb yr⁻¹ in the east and 0.31 [95 % CI: 0.21, 0.41] ppb yr⁻¹ in the west (Fig. 3, Table S4).

Figure 3 (a) Population-weighted exposure concentrations of the MDA1 human-health metric from 2000 to 2015. The meteorologically adjusted trend is in black with the slope in the inset. (b) The 2000–2015 population-weighted trends (ppb yr⁻¹) of the MDA1 (green) and MDA8 (red) metrics. The west–east divide is made along the 95^∘ W meridian and the whiskers span the 95 % confidence interval. (c) Population-weighted exposure concentrations of the MDA8 human-health metric from 2000 to 2015. The meteorologically adjusted trend is in black with the slope in the inset. Tabulated values of these plots can be found in Table S2 and S4.

In contrast, the MDA8 metric featured more modest decreases in the Southeast USA and scattered areas with increasing trends (Fig. 2). This divergence between the two human-health metrics is due to the different averaging periods (i.e., the traditional “ozone season” vs. an annual average). If only summer months were considered when calculating the MDA8 metric (i.e., MDA8-JJA), the two trends would be spatially and quantitatively consistent (Fig. 2). However, O₃ increases during the winter months (i.e., MDA8-DJF) partially compensated for the summer decreases, resulting in no discernable trend for the national population-weighted MDA8 metric (Fig. 3). After adjusting for meteorology, the national population-weighted MDA8 trend from 2000 to 2015 is −0.02 [95 % CI: 0.01, −0.04] ppb yr⁻¹ (Fig. 3). Similar to the MDA1 metric, trends featured an east–west divide. It is interesting to note that the western MDA8 trends were slightly positive and the eastern MDA8 trends were slightly negative (Table S4). Separately, prior studies (e.g., Bloomer et al., 2010; Cooper et al., 2012; Parrish et al., 2012; Clifton et al., 2014; Simon et al., 2015; Strode et al., 2015; Fleming et al., 2018) have highlighted the existence of a seasonal shift in the distribution of O₃ concentrations throughout the United States. We find that these shifts have not only manifested in contrasting seasonal trends (i.e., summer decreases vs. winter increases), but have also led to changes in the dominant months of O₃ exposure. For example, the population-weighted exposure concentrations during the spring and summer months (MDA8-MAM vs. MDA8-JJA) were nearly equivalent from 2013 to 2015 (Table S2).

It should be noted that comparing previously reported seasonal trends of O₃ is difficult due to varying study periods, averaging metrics, and selection of monitoring networks. Oftentimes, rural locations are highlighted, enabling the isolation of trends in background O₃ concentrations or the minimization of the influence of nearby changes in anthropogenic emissions (e.g., Jaffe and Ray, 2007; Cooper et al., 2012; Jaffe et al., 2018). These study-specific choices can effect conclusions. For example, Simon et al. (2015) report that rural O₃ monitors more often feature statistically significant decreases in national mean MDA8 O₃ during summer months and urban O₃ monitors more often feature statistically significant increases in national mean MDA8 O₃ during winter months. For this study, since our focus is on changes in exposure, we incorporate all available observational data, including monitors in urban cores. As such, when compared to prior studies, our conclusions regarding O₃ trends may be different.

Cooper et al. (2012), using rural monitoring data spanning 1990–2010, reported a −0.45 and a +0.10 ppb yr⁻¹ trend in daytime O₃ during summer months for the eastern and western USA, respectively. However, both trends featured wide ranges. Jaffe et al. (2018), using a limited number of high-elevation, rural monitoring sites, reported decreasing trends of median summertime O₃ between 2000 and 2016 at most analyzed locations, with stronger decreases in the east than west (∼1 ppb yr⁻¹ vs. ∼0.5 ppbẏr). Lin et al. (2017) also used rural monitoring data, but increased the coverage to include 1988–2014 and found a 0.4–0.8 ppb yr⁻¹ decreasing trend of median MDA8-JJA concentrations in the eastern USA and mixed trends in the west. Fleming et al. (2018) incorporated both urban and nonurban monitors and showed that the observed magnitude of several warm season human-health ozone metrics is similar for North American urban and nonurban sites and that the trends are only slightly smaller for the urban areas. Broadly, we also found a dramatic divide between east and west summertime O₃ exposure trends, but our results did feature some differences from prior studies. For example, our exposure-focused (i.e., population-weighted) estimates of eastern USA trends are similar to the mean reported by Cooper et al. (2012) and on the low end of that reported by Lin et al. (2017). We also found a consistent decreasing trend in western MDA8-JJA exposure, as well as smaller levels of trend uncertainty (Table S4).

Cooper et al. (2012) also reported a uniform east and west increase in rural wintertime O₃ concentrations of 0.12 ppb yr⁻¹. However, the exclusive selection of rural monitors precludes the extrapolation of those results to estimate exposure trends. This is well illustrated by Simon et al. (2015), who used an extensive network of 1998–2013 observations to show that there was a strong rural–urban divide in mean winter O₃ trends, with increasing trends more prevalent in urban areas. Indeed, when compared to Cooper et al. (2012), we found a much stronger trend increase in MDA8-DJF exposure (Table S4). Our results indicate that the national trend in MDA8-DJF exposure was +0.33 [95 % CI: 0.37, 0.28] ppb yr⁻¹ (Table S4), with a near-uniform increase in both the east and west.

Figure 4 Trends of the M12 (a; ppb yr⁻¹) and AOT40 (b, c, d; ppmh yr⁻¹) agriculture metrics from 2000 to 2015. Note: JAS: July–September; MAM: March–May; MJJ: May–July. The p values from these trends can be found in Fig. S6.

3.2.2 Crop-loss-relevant metrics

Since the averaging months for crop-loss metrics are dependent on crop variety, the magnitude and trends can feature distinct patterns (Table S3). All maize and soybean metrics envelop the months of July–September. As such, and consistent with the MDA1 results, this averaging period yielded widespread decreases for both the M12 and AOT40 metrics, with the strongest reductions in the southeast and California (Fig. 4a, b). However, both of these crops are predominantly grown in the Midwest and Great Plains (Fig. S7). These regions generally experienced smaller trend reductions. Nationally, the agriculture-weighted trends of the M12 metric for maize and soybeans were −0.35 [95 % CI: −0.17, −0.54] ppb yr⁻¹ and −0.39 [95 % CI: −0.19, −0.59] ppb yr⁻¹ (Fig. 5; Table S4). The agriculture-weighted trends of the AOT40 metric for maize and soybeans were −0.35 [95 % CI: −0.18, −0.51] ppmh yr⁻¹ and −0.39 [95 % CI: −0.21, −0.56] ppmh yr⁻¹ (Fig. 5; Table S4). After adjusting for meteorology, the mean trend for both metric and crop pairings was reduced marginally, suggesting that meteorological factors played a small role in the net trends from 2000 to 2015 (Fig. 5b).

Figure 5 (a) Agriculture-weighted exposure concentrations of the M12 agriculture metrics from 2000 to 2015. The meteorologically adjusted trend of each metric is in black with the slope in the inset. (b) The 2000–2015 agriculture-weighted trends (ppb yr⁻¹ for M12, ppmh yr⁻¹ for AOT40) of the M12 and AOT40 metrics. The whiskers span the 95 % confidence interval. (c) Agriculture-weighted exposure concentrations of the AOT40 agriculture metrics from 2000 to 2015 (trends of AOT40 soybean and AOT40 maize are nearly equivalent). The meteorologically adjusted trend of each metric is in black with the slope in the inset. Note: variable averaging periods are considered, reflecting differences in crop harvest seasons. Tabulated values of the left and right plots can be found in Tables S3 and S4.

Both agriculture-weighted AOT40 averaging periods for wheat (MAM and MJJ) featured decreasing, but considerably different, trends (Fig. 5, bottom panels). These trend differences again highlight the seasonal shift in O₃ concentrations. From 2000 to 2006, the AOT40-MJJ wheat metric was ∼40 %–60 % higher than the AOT40-MAM wheat metric (Table S3). However, by 2014, both metrics were nearly equal. The 40 ppb accumulation threshold applied in the AOT40 calculation also amplifies this convergence. Towards the end of the study period, mean daytime O₃ concentrations in the Midwest and Great Plains had decreased sufficiently for the two metrics to nearly intersect.

Figure 6 Annual trends in the daytime 2 m temperature and daytime 2 m specific humidity from 2000 to 2015.

We posit that the influence of meteorology on the agriculture-weighted trends, as indicated by the marginal difference in the mean of the meteorologically adjusted and non-adjusted trends, is primarily due to two factors. Prior analysis has shown that two important meteorological variables influencing O₃ include temperature and humidity (Camalier et al., 2007; Jacob and Winner, 2009). The temperature–O₃ mechanism is a function of increasing temperatures promoting peroxyacetyl nitrate decomposition (leading to ozone increases near NO_x source regions, but decreases in remote areas; Doherty et al., 2013) and increases in isoprene emissions. The humidity–O₃ mechanism is a function of increasing water vapor concentrations, promoting O₃ chemical destruction. According to the MERRA-2 reanalysis product, the Midwest and Great Plains regions featured both decreasing trends in daytime 2 m temperature and increasing trends in daytime 2 m specific humidity (Fig. 6). In addition, acute O₃ episodes are notably sensitive to particular meteorological variables (Russell et al., 2016; Fix et al., 2018), such as temperature, providing an environment where meteorological variability can disproportionately influence the magnitude of AOT40 values.

Figure 7 (a) Annual estimates of normalized premature mortalities (number per 100 000 people per year) attributable to long-term O₃ exposure using the MDA1 and MDA8 averaging metrics and exposure–response function from J2009 and T2016, respectively. The shaded region reflects the confidence interval reported in each underlying epidemiological study. (b) The 2000 vs. 2015 percent contributions of population aging, changing baseline mortality rates, and long-term O₃ exposure to net normalized premature mortalities using the MDA1 and MDA8 averaging metrics and exposure–response function from J2009 and T2016, respectively. Tabulated values of (a) can be found in Table S5.

3.3 Estimates of long-term O₃ exposure impacts

3.3.1 Human health

Human-health impacts, reported as the estimated number of premature respiratory mortalities attributable to long-term O₃ exposure per 100 000 people, were strongly dependent on the choice of exposure–response relationship (Fig. 7). First, the T2016 results reported nearly double the estimated human-health impacts attributable to long-term O₃ exposure. For example, in 2010, the J2009 and T2016 estimated impacts were ∼5.4 [95 % CI: 1.8, 8.7] and ∼11.3 [95 % CI: 7.9, 14.5] premature mortalities per 100 000 people, respectively. Second, the diverging trends of the two exposure metrics (Fig. 3) are reflected in the estimated impacts (Fig. 7). Between 2000 and 2015, the MDA1 population-weighted exposure concentration decreased from ∼53.7 to ∼48.3 ppb (Table S2). As a result, the estimated human-health impacts using the J2009 parameters decreased from ∼6.0 [95 % CI: 2.0, 9.7] to ∼5.0 [95 % CI: 1.7, 8.0] premature mortalities per 100 000 people (Table S5). In contrast, the MDA8 population-weighted exposure concentration decreased from ∼39.9 to ∼39.1 ppb, yet the impacts using the T2016 parameters increased from ∼10.8 [95 % CI: 7.6, 13.8] to ∼11.3 [95 % CI: 7.9, 14.5] premature mortalities per 100 000 people (Fig. 7 and Table S5). These differences in estimated impacts are not only due to changes in exposure. Over this period, an aging population structure promoted increased susceptibility to O₃ impacts. In addition, depending on the age bin, baseline mortality rates for respiratory diseases either marginally decreased or remained approximately stable.

While impacts due to changes in exposure for both metrics decreased between these end points, albeit by different magnitudes (blue bars in Fig. 7b; −25.5 % vs. −5.7 %), these other determinants played a strong role in modulating the estimated impact trends (Fig. 7b). The net changes in 2015 vs. 2000 normalized human-health impacts using the J2009 and T2016 exposure–response relationships and averaging metrics were −17.8 % and +4.7 %, respectively (black bars in Fig. 7b). In both calculations, an aging population structure substantially eliminated much of the gains from exposure decreases (+15.5 %). Changing baseline mortality rates were more modest, decreasing both calculations by 4.7 %.

The differences in estimated human-health impacts when using the J2009 and T2016 exposure–response relationship and averaging metrics reported here are consistent with prior studies (Malley et al., 2017; Seltzer et al., 2018; Shindell et al., 2018). That is, the estimated human-health impacts when using the T2016 exposure–response relationship and averaging metric are considerably higher than the results computed when using the J2009 parameters. However, to our knowledge, the historical evolving differences between the two have yet to be shown. For example, the T2016 results were ∼80 % higher than the J2009 results in 2000 (Table S5). By 2008, the T2016 results were nearly double the J2009 results, and this difference continued to grow over time (∼130 % in 2015).

Between 2000 and 2015, our net estimated premature mortalities attributable to long-term O₃ exposure in the USA ranged from ∼14 500 to 19 200 when using the J2009 parameters and from ∼29 800 to 37 600 when using the T2016 parameters. These results are lower than analogous prior studies that are based solely on CTM estimates of O₃ exposure. An exception is Zhang et al. (2018), who found comparable results when using the J2009 epidemiological study. However, Zhang et al. (2018) report a 13 % increase in premature mortalities attributable to long-term O₃ exposure in the United States between 1990 and 2010, despite O₃ decreases. We find a ∼6.7 % decrease in premature mortalities attributable to long-term O₃ exposure, albeit over 2000–2015. This is likely due to the dramatic decreases in O₃ precursor emissions that occurred post-2000 (Xing et al., 2013; Simon et al., 2015).

3.3.2 Crop loss

Agriculture impacts for each of the crop varieties considered here decreased from 2000 to 2015 (Fig. 8). When using the M12 metric, the estimated national RYL values for maize and soybeans in 2000 were 4.6 % and 16.3 %, respectively (Fig. 8 and Table S5). These values decreased to 2.9 % and 11.2 % in 2015. When using the AOT40 metric, the estimated national RYL values for maize, soybeans, and wheat for the year 2000 were 3.4 %, 11.9 %, and 12.1 %, respectively. By 2015, these RYL values dropped to 1.6 %, 4.8 %, and 9.4 %, respectively. Broadly, these estimated agriculture yield impacts are comparable to the global “ozone yield gaps” (i.e., RYL) modeled by Mills et al. (2018a), who considered the flux-based, stomatal uptake of O₃ for each crop.

Figure 8 Estimates of the national relative yield loss for a variety of commercial crops using ANN-calculated exposure metrics. Tabulated values of this plot can be found in Table S5.

Several other characteristics are consistent among all of the crop varieties and metric combinations considered here. For one, estimated RYL featured sizable inter-annual variability, indicating that the impacts calculated from a single year might not be representative of a particular period. For example, the RYL for soybeans, when using the AOT40 metric, increased from 7.8 % in 2004 to 11.6 % in 2005 – a nearly ∼50 % increase. Second, similar to Van Dingenen et al. (2009) and Lapina et al. (2016), impacts were consistently higher when utilizing the M12 metric and the associated concentration–response functions. These differences also became amplified over time. The RYL for soybeans in 2000 using the M12 metric (16.3 %) was ∼37 % higher than the RYL using the AOT40 metric (11.9 %). This difference increased to ∼135 % (11.2 % vs. 4.8 %) by 2015 (Table S5). These diverging trends occur for two reasons. First, the daytime O₃ concentrations approached the AOT40 threshold of 40 ppb post-2007 (Table S3). This decrease precipitated disproportional improvements in AOT40-calculated RYL. Second, and to a lesser degree, the slopes of the two soybean concentration–response functions are different (Fig. S8).

4 Uncertainties, limitations, and additional remarks

Studies quantifying the health impacts attributable to long-term PM_2.5 exposure oftentimes use higher-resolution products (e.g., $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$) that harness satellite data (e.g., Apte et al., 2015; Cohen et al., 2017; van Donkelaar et al., 2019). However, a number of complications prevent such products for surface O₃ (Duncan et al., 2014). Regardless, we believe this $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ product is of sufficient resolution to estimate long-term O₃ exposure for a number of reasons. First, O₃ features a residence time on the order of hours to days in the lower troposphere and in urban environments (Parrish et al., 2012; Monks et al., 2015), providing sufficient time for localized and regional mixing. Second, unlike short-term O₃ exposure, long-term O₃ exposure is less sensitive to singular events that are more heterogeneous in space and time. Third, regional CTM studies report only marginal differences in O₃ concentrations and estimated impacts when scaling from 12 km to resolutions comparable to those in this analysis (Punger and West, 2013; Gan et al., 2016).

In terms of impacts, there is evidence that O₃ affects more than what was presented here. For example, several epidemiological studies suggest that human-health impacts may extend to cardiovascular mortality (Jerrett et al., 2013; Crouse et al., 2015; Cakmak et al., 2016; Turner et al., 2016; Lim et al., 2019). Separately, our analysis applied a log-linear exposure–response function when performing the human-health calculations since it is most common method applied in the community. There is evidence that this relationship may instead be linear (Di et al., 2017). For agriculture, the exposure–response functions utilized here are “pooled” from studies featuring a limited number of cultivars grown in the USA and Europe (Van Dingenen et al., 2009). While considered reliably representative of the commonly grown cultivar population in these regions, extrapolation of these relationships to the national level may introduce additional uncertainty. In addition, the methodology selected here does not take into account changes in plant conditions that may limit or exacerbate conditions which influence the opening of stomata and the ability of a plant to uptake O₃, such as temperature and soil moisture. The results presented here also demonstrate the need for additional epidemiological studies to test the utility of common averaging metrics that are used when estimating health impacts. Specifically, clarity is needed regarding whether long-term O₃ health impacts are more sensitive to peak averaged (i.e., the MDA1 metric) or annually averaged (i.e., the MDA8 metric) O₃ exposure.

Finally, long-term trends of O₃ are driven by a number of mechanisms, including intercontinental transport (Fiore et al., 2009; Lin et al., 2012, 2017) and methane concentrations (Fiore et al., 2002; Shindell et al., 2017; Lin et al., 2017). For example, Lin et al. (2015) conclude that rising Asian emissions and global methane have played a key role in the increase in western USA springtime O₃ from 1995 to 2014. These drivers merit additional study, with an emphasis on exploring seasonal differences that influence impact metrics. Furthermore, inclusion of observations from the most recent years (i.e., 2015–2017) should be targeted. Since Chinese emissions of NO_x peaked in 2012 (Zheng et al., 2018), our current and future estimates of intercontinental transport influences on background O₃ might warrant revisiting.

5 Conclusions

Through the application of artificial neural networks, we empirically model the magnitude and trends of long-term (i.e., seasonal, annual) ambient O₃ over the continental United States from 2000 to 2015. We then used these estimates of long-term O₃ exposure to generate a measurement-based assessment of impacts on human health and crop yields. All metrics with averaging periods spanning the traditional O₃ season (i.e., warm months) featured peak exposure in 2002, with net decreases over the course of the study period. For example, the population-weighted, April–September average of the daily 1 h maximum O₃ concentration (i.e., MDA1 from Jerrett et al., 2009) decreased by 0.43 [95 % CI: 0.28, 0.57] ppb yr⁻¹ between 2000 and 2015. In contrast, there was little change in the population-weighted, annual average of the maximum daily 8 h average O₃ concentration (i.e., MDA8 from Turner et al., 2016) between 2000 and 2015. There were compensating seasonal effects, with wintertime O₃ increases and summertime O₃ decreases, yielding a net population-weighted trend of −0.03 [95 % CI: 0.04, −0.10] ppb yr⁻¹. Human-health metric trends also featured an east–west divide, with stronger decreases in the eastern USA. All agriculture-weighted crop-loss metrics featured decreasing trends over the study period.

Human-health impacts were quantified in terms of the estimated number of premature respiratory mortalities attributable to long-term O₃ exposure per 100 000 people. Crop-loss impacts were quantified in terms of the estimated national relative yield loss for a variety of commercial crops. Normalized human-health impact estimates decreased by ∼18 % and increased by ∼5 % when using the Jerrett et al. (2009) and Turner et al. (2016) averaging metrics and parameters, respectively. In both cases, exposure changes and an aging population structure played a substantial role in modulating these trends. When using the M12 metric, the relative yield loss (RYL) due to O₃ exposure for maize and soybeans improved by 1.7 % and 5.1 %. When using the AOT40 metric, the net benefits were greater, with the RYL for maize, soybeans, and wheat improving by 1.9 %, 7.1 %, and 2.7 %, respectively. These different responses are mainly due to the daylight O₃ concentrations approaching the 40 ppb AOT40 threshold by the end of the study period. Overall, these results provide a measurement-based estimate of long-term O₃ exposure over the United States, quantify the historical trends of such exposure, and illustrate how different conclusions regarding historical impacts can be made through the use of varying metrics.