Extremal dependence between temperature and ozone over the continental US

The co-occurrence of heat waves and pollution events and the resulting high mortality rates emphasize the importance of the co-occurrence of pollution and temperature extremes. Through the use of extreme value theory and other statistical methods, tropospheric surface ozone and temperature extremes and their joint occurrence are analyzed over the United States during the summer months (JJA) using measurements and simulations of the present and future climate and chemistry. Five simulations from the Chemistry-Climate Model Initiative (CCMI) reference experiment using specified dynamics (REFC1SD) were analyzed: the CESM1 CAM4-chem, CHASER, CMAM, MOCAGE and MRI-ESM1r1 simulations. In addition, a 25-year present-day simulation branched off the CCMI REFC2 simulation in the year 2000 and a 25-year future simulation branched off the CCMI REFC2 simulation in 2100 were analyzed using CESM1 CAM4-chem. The last two simulations differed in their concentration of carbon dioxide (representative of the years 2000 and 2100) but were otherwise identical. In general, regions with relatively high ozone extremes over the US do not occur in regions of relatively high temperature extremes. A new metric, the spectral density, is developed to measure the joint extremal dependence of ozone and temperature by evaluating the spectral dependence of their extremes. While in many areas of the country ozone and temperature are highly correlated overall, the correlation is significantly reduced when examined on the higher end of the distributions. Measures of spectral density are less than about 0.35 everywhere, suggesting that at most only about a third of the time do extreme temperatures coincide with extreme ozone. Two regions of the US have the strongest measured extreme dependence of ozone and temperature: the northeast and the southeast. The simulated future increase in temperature and ozone is primarily due to a shift in their distributions, not to an increase in their extremes. The locations where the right-hand side of the temperature distribution does increase (by up to 30 %) are consistent with locations where soil–moisture feedback may be expected. Future changes in the right-hand side of the ozone distribution range regionally between +20 % and −10 %. The location of future increases in the high-end tail of the ozone distribution are weakly related to those of temperature with a correlation of 0.3. However, the regions where the temperature extremes increase are not located where the extremes in ozone are large, suggesting a muted ozone response.

1 Introduction

The European heat wave of 2003, the Russian heat wave of 2010, and the extreme pollution and mortality increase that accompanied both events underline the danger of heat waves and the accompanying air pollution. Summertime increases in temperature are expected in the next century in all climate scenarios (Collins et al., 2013), with future heat waves expected to be more intense, more frequent and longer lasting (e.g., Meehl and Tebaldi, 2004). Here we examine the relationship between temperature extremes and ozone extremes in measurements and in current and future model simulations. An analysis of the joint extremes in ozone and temperature together may be particularly important as their joint impact on mortality is likely to be nonlinear (Wilson et al., 2014; Dear et al., 2005; Ren et al., 2008).

Over most of the US temperature is the first meteorological covariate with ozone (Porter et al., 2015). The relation between ozone and temperature is complex: it is determined not only by temperature-dependent ozone chemistry (Pusede et al., 2015) but also by other processes that correlate with temperature: for example, through meteorological factors such as stagnation events or cloud cover (e.g., see Jacob and Winner, 2009) or through temperature-dependent emissions (e.g., Weaver et al., 2009). The ozone–temperature relationship is often measured with a linear slope (e.g., Steiner et al., 2010). Increases in ozone with temperature have been reported in the range from 0 to 6 ppbv ^∘C⁻¹ depending on details of the analysis (Brown-Steiner et al., 2015). The mechanisms accounting for variations in the ozone–temperature slope are still uncertain but can be at least partially explained by the emission regime: the slope generally increases as ozone precursor emissions increase (e.g., Pusede et al., 2015).

When it comes to extreme values, the relationship between temperature and ozone becomes more complicated (e.g., Steiner et al., 2010; Shen et al., 2016), such that an overall linear slope fit does not necessarily capture the relationship. The extremal dependence between ozone and temperature has been explored using various methods. Sun et al. (2017) calculated the conditional probability of a high ozone day (ozone above the 90th percentile) given a high temperature day (temperature above the 90th percentile). They found probabilities that range from approximately 50 % in the northeastern US to somewhat less than 20 % in the western US. Schnell and Prather (2017) and Zhang et al. (2017) calculate the joint probability that ozone and temperature are extreme (above the 95th percentile) compared to the probability that either one of them is high. Schnell and Prather (2017) find that high temperature events and high ozone events co-occur up to 50 % of the time over the northeastern US between 1 April and 30 September. Zhang et al. (2017) obtain a qualitatively similar geographic pattern in the joint extremes of temperature and ozone. Here we propose a new method to measure the joint extremes of temperature and ozone based on the spectral dependence of the extremes. Changes in the future relation between ozone extremes and temperature extremes depend on (i) changes in the nature of future temperature extremes and (ii) the impact of these extremes on ozone. Globally the future increase in extreme temperatures (temperatures at the 20-year return period) in CMIP3 are similar to the increase in mean temperature (Seneviratne et al., 2012), although there are some important regional exceptions. This would suggest that for the most part the future temperature probability distribution simply shifts to higher temperatures but does not change in shape consistent with measured trends (McKinnon et al., 2016). Conversely, many studies suggest that the ozone distribution will increase predominantly on the high end due to changes in climate (e.g., see Weaver et al., 2009), although not all (e.g., Rieder et al., 2015). Wu et al. (2008) find increases in the high end of the probability distribution of both temperature and ozone in the midwestern US in 2050, attributing this to an increase in stagnation episodes with soil–moisture feedbacks impacting the temperature distribution. At high enough temperatures (>312 K), Steiner et al. (2010) find that the ozone increase with temperature is suppressed. Steiner et al. (2010) hypothesize that this is due to the diminished role of peroxyacetyl nitrate (PAN) chemistry and isoprene emissions at high temperatures. Shen et al. (2016) show that ozone suppression at high temperature occurs at 23 % of the CASTNET sites but hypothesize that the suppression is meteorologically induced. Based on the statistical relationship between ozone and temperature in the present climate, Shen et al. (2016) predict that future temperatures will lead to an average increase of 2.6 ozone violations per year in 2060 across the US. Meehl et al. (2018) examine the impact of an increase in future heat waves on ozone in two sets of future simulations: one with changing anthropogenic precursor emissions following Representative Concentration Pathway 6 (RCP6) and one in which the anthropogenic emissions remain fixed.

This article is organized as follows: in Sect. 2, we describe the data sets and model simulations used; in Sect. 3, we introduce the statistical procedures used to quantify the relationship between ozone and temperature. In Sect. 4, we present the results and then discuss these results in Sect. 5. Section 6 gives the conclusions.

2 Data and model descriptions

In this study we examine the simulated and measured relationship between ozone and temperature extremes over the US. In particular we analyze a number of specified dynamics REFC1SD simulations from the Chemistry-Climate Model Initiative (CCMI) (see Eyring et al., 2013b) for the period from 1992 to 2010. This allows a robust evaluation of simulated ozone and temperature extremes using analyzed meteorological fields and changing emissions against measurements. We also examine the impact of climate change on ozone and temperature extremes, comparing simulations of the current and future climate with fixed emissions. These latter simulations are free-running simulations in that the meteorology, sea surface temperatures (SSTs) and sea ice are calculated internally within the simulations. The free-running simulation of the current climate is compared with the REFC1SD simulations and the available measurements.

Most of the analysis in this paper emphasizes simulations with the Community Atmospheric Model with chemistry (CAM4-chem) within the Community Earth System Model (CESM1) (Lamarque et al., 2012). Brown-Steiner et al. (2015) evaluate both the specified dynamics and free-running model configurations of CAM4-chem against measurements over the US, including comparisons of ozone return periods. The horizontal grid resolution in the CESM1 simulations analyzed here is $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$; the free-running simulations have 26 vertical levels while the CESM1 REFC1SD simulation has 56 vertical levels. The CESM1 REFC1SD simulation (see Tilmes et al., 2016) uses analyzed meteorological data from Modern-Era Retrospective Analysis for Research and Applications (MERRA) from 1992 to 2010 and time-changing anthropogenic and biomass burning emissions as specified in Table 1. In the Supplement we include results from the REFC1SD simulations in an additional four models: the CHASER, CMAM, MOCAGE and MRI models. The number of REFC1SD models analyzed is limited to those with sufficient output to derive the maximum daily temperature and the maximum daily 8 h average ozone concentrations (MDA8). In the CMAM and MRI simulations both MDA8 ozone and daily maximum temperature are available daily; in the CHASER and MOCAGE simulations only daily MDA8 ozone data are available. Details on the additional model simulations are given in Morgenstern et al. (2017).

Table 1 Details and descriptions for each model.

¹ Greenhouse gas. ² Sea surface temperature. ³ Coupled Model Intercomparison Project. ⁴ Granier et al. (2011). ⁵ Guenther et al. (2012). ⁶ Hadley Centre Sea Ice and Sea Surface Temperature data set (Titchner and Rayner, 2014). ⁷ Modern-Era Retrospective Analysis for Research and Applications (Rienecker et al., 2011). ⁸ Assessment Report 5 (Eyring et al., 2013a).

The analysis of how extremes change with climate is limited to the CESM1 simulations. In the present-day free-running CESM1 simulation (the GCM2000 simulation) the CO₂ concentration is specified at 369 ppm, representative of the year 2000; in the future simulation (the GCM2100 simulation) the CO₂ concentration is specified at 669 ppm, representative of the 2100 concentration of CO₂ in the Representative Concentration Pathway 6 (RCP6) (see Table 1). The concentrations of all other greenhouse gases including methane are fixed at their year 2000 concentrations in both these simulations. Biogenic emissions are also fixed and are representative of the year 2000. Both the GCM2000 and GCM2100 simulations are 25-year simulations branched off of the CCMI CESM REFC2 simulations in the year 2000 and the year 2100, respectively (Tilmes et al., 2016). The first 5 years of each simulation are used as spin-up, with the latter 20 years analyzed. The global mean temperature change over the continental US between GCM2000 and GCM2100 is 2.1 ^∘C, while the temperature difference in the parent CCMI REFC2 simulations following RCP6 is 2.8 ^∘C. The smaller temperature increase between the GCM2000 and GCM2100 simulations is likely due to the fact that the emissions of greenhouse gases and short-lived forcing agents are held constant at the year 2000 levels in both simulations. In particular the aerosol emissions remain the same.

Hourly measured ozone and temperature data are taken from 23 CASTNET (Clean Air Status and Trends Network) stations with a nearly continuous data record during the period from 1992 to 2013 for the months of June, July and August (92 days each summer). In addition, to enhance the data record, we included two additional stations (Beaufort, NC, and Lassen Volcanic, CA) at which the first 2 or 3 years of data were missing, respectively. See Fig. 2 for station locations. CASTNET sites are situated to sample regional ozone concentrations so as to minimize the more local impact of urban areas. We supplement the CASTNET data with temperature and ozone measurements from the Environmental Protection Agency (EPA) Air Quality System (AQS) Data Mart for the years 1992–2010. This gives an additional 124 stations with nearly complete ozone and temperature data (see Supplement). Ozone data from the first model level provide a good estimate of 10 m ozone concentrations as measured by CASTNET (Brown-Steiner et al., 2015). The maximum daily 2 m temperature is used in both the CASTNET measurements and the simulations.

To render the data approximately stationary on both an interannual and seasonal basis, we adopt the procedures in Phalitnonkiat et al. (2016). Formally, let x_y,d represent the data on day d in year y, where x refers to either daily maximum temperature or MDA8 ozone. Since there are 91 days included in each summer period and 20 years (19 years for REFC1SD), d=1 refers to 1 June and d=91 refers to 30 August.

To minimize year-to-year variability so as to minimize any ozone trends while still keeping extreme data relevant, for each year y, we take the average of the data over that year but omit a number (a) of the highest values. That is, for a fixed year y, the resulting average is

$$\begin{array}{}\text{(1)}& m_{y,a}:={\displaystyle \frac{\mathrm{1}}{D-a}}\sum _{i=\mathrm{1}}^{D-a}{x}_{y,\left(i\right)},\end{array}$$

where D=91 is the total number of days for each year, and x_y,(i) is the order statistic of the fixed year y: $x_{y,\left(\mathrm{1}\right)}\le x_{y,\left(\mathrm{2}\right)}\le \mathrm{\dots }\le x_{y,\left(D\right)}$. Then, we calculate a daily ozone deviation:

$$\begin{array}{}\text{(2)}& {\widehat{x}}_{y,d}^G=x_{y,d}-m_{y,a}.\end{array}$$

In our analysis, we use a=10 as the default value, which preserves about 11 % of the extreme data. Sensitivity tests at a number of stations suggest the result is not sensitive to a. To eliminate seasonal effects, we average ${\widehat{x}}_{y,d}^G$ for each day d over all years Y=20 (or Y=19 for REFC1SD). That is, for each day d, we calculate

$$\begin{array}{}\text{(3)}& M_d={\displaystyle \frac{\mathrm{1}}{Y}}\sum _{y=\mathrm{1}}^Y{\widehat{x}}_{y,d}^G.\end{array}$$

M_d is then smoothed using a local polynomial regression since our sample size is rather small. In order not to overburden the notation, we will still use the notation M_d for the smoothed values of the estimates. Then we normalize the data by

$$\begin{array}{}\text{(4)}& {\widehat{x}}_{y,d}^{\mathrm{DS}}={\displaystyle \frac{{\widehat{x}}_{y,d}^G-M_d}{{\mathrm{SD}}_d}},\end{array}$$

where ${\mathrm{SD}}_d=\sqrt{\frac{\mathrm{1}}{Y}{\sum }_{y=\mathrm{1}}^Y({\widehat{x}}_{y,d}^G-M_d{)}^{\mathrm{2}}}$ is the standard deviation of day d. Later in the text, we refer to Eq. (4) as a normalized scale.

In addition to the transformations from Phalitnonkiat et al. (2016), we add another procedure to revert the normalized scale data back to its original scale while keeping the stationarity. That is, we rescale ${\widehat{x}}_{y,d}^{\mathrm{DS}}$ back to its original scale $\left({\widehat{x}}_{y,d}^{\mathrm{res}}\right)$ by using the formula

$$\begin{array}{ll}{\displaystyle }{\widehat{x}}_{y,d}^{\mathrm{res}}=& {\displaystyle }{\widehat{x}}_{y,d}^{\mathrm{DS}}\times \left({\displaystyle \frac{\mathrm{1}}{D}}\sum _{d^{\prime }=\mathrm{1}}^D{\mathrm{SD}}_{d^{\prime }}\right)+{\displaystyle \frac{\mathrm{1}}{D}}\sum _{d^{\prime }=\mathrm{1}}^D{M}_{d^{\prime }}\\ \text{(5)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{Y}}\sum _{y^{\prime }=\mathrm{1}}^Y{m}_{y^{\prime },a},\end{array}$$

where Y=20 (or Y=19 for REFC1SD).

3 Methodology

In this study in addition to using conventional methods, such as correlations, to quantify the relationship between temperature and ozone, we also propose a novel metric to capture the relationship between ozone and temperature extremes. Correlation coefficients are inadequate for capturing the relationship between the extremes of two variables since they are estimated from all observations and extremes represent a small percentage of these observations. An alternative metric is proposed using only the largest values of two variables. After some transformations which act to normalize the two variables (see Appendix B), their extreme dependency is characterized by a probability density function (pdf) that measures the angular density when the variables are plotted against each other. The area under the pdf is 1 by definition and the range of the pdf is $\left[\mathrm{0},\frac{\mathit{\pi }}{\mathrm{2}}\right]$. If the mass of the pdf is concentrated near 0 or near $\frac{\mathit{\pi }}{\mathrm{2}}$, extremes of the two variables are unlikely to be significant at the same time, which points to an independence of the extremes. Conversely, if the mass of the pdf is concentrated away from the end points $\mathrm{0},\frac{\mathit{\pi }}{\mathrm{2}}$, then simultaneous extremes of the two variables are likely. We refer to the procedure which normalizes the tails of the data so that the method described above works as the ranks method (see Appendix B).

Since the area under the curve from 0 to π∕2 is 1, we can consider only the area of the “middle” part, which we define to be between $\frac{\mathit{\pi }}{\mathrm{8}}$ and $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{8}}$, to represent the extreme dependence between two variables. Denote this amount by φ:

$$\begin{array}{}\text{(6)}& \mathit{\phi }:=\text{area}\left[{\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}},{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{8}}}\right].\end{array}$$

See the detailed explanation in Appendix C. Note that the range of φ is [0,1], where φ=1 refers to extreme dependence and φ=0 implies extreme independence.

Table 2 Rescaled data. MDA8 ozone averages (ppb) and daily maximum temperature averages (^∘C) in different regions over the US from CASTNET data and corresponding grid points from the CESM1 REFC1SD, the GCM2000 and the GCM2100 simulations. The averages are calculated from MDA8 (for ozone) and daily maximum (for temperature). Standard deviation (SD) is calculated among the stations in each region. The averages of ozone and temperature from each region are reported in italics (including all continental grid points in each box in Fig. 2a). The italics under “All” are averages of all points over the continental US.

Figure 1 Examples that show the correlation and extremal measure of dependence between two variables are not necessarily the same. The plots in the left column (a, c) use the data generated by Gaussian random vectors with correlation ρ=0.7 and each component has n=10 000 points sampled from N(10,1), where N(μ,σ²) is the normal distribution with mean μ and variance σ². The data are moderately correlated, while they have low extreme dependence (true φ=0; estimated φ=0.223). The plots in the right column (b, d) use the data generated by $({\mathrm{Var}}_{\mathrm{1}},{\mathrm{Var}}_{\mathrm{2}})=(Y_{\mathrm{1}},Y_{\mathrm{2}})$ with probability 0.8 and $({\mathrm{Var}}_{\mathrm{1}},{\mathrm{Var}}_{\mathrm{2}})=(Z,Z)$ with probability 0.2, where $(Y_{\mathrm{1}},Y_{\mathrm{2}})\sim N(\mathit{\mu }=[\mathrm{10},\mathrm{10}{]}^T,\mathbf{\Sigma }=[\mathrm{1},-\mathrm{0.9};-\mathrm{0.9},\mathrm{1}])$ follows a bivariate normal, μ is the mean vector and Σ is the covariance matrix, and $Z\sim N(\mathit{\mu }=\mathrm{10},{\mathit{\sigma }}^{\mathrm{2}}=\mathrm{9})$. The sample size is also n=10 000. The plots show the existence of tail dependence by having high angular density near $\frac{\mathit{\pi }}{\mathrm{4}}$ (φ=0.75); however, there is low correlation (true ρ=0; estimated $\mathit{\rho }=-\mathrm{0.05}$).

Figure 1 shows different scenarios of correlation and extreme dependence. Figure 1a gives a scenario of data with high correlation, yet the extremes of the data are only moderately dependent (Fig. 1c). In contrast, Fig. 1b gives an example of data with low correlation but highly dependent extremes (Fig. 1d).

4 Results

In this section we compare measured and simulated temperature and ozone records separately (Sect. 4.1) and then their joint dependence is analyzed in Sect. 4.2. The extremes of ozone and temperature and their extremal dependence is emphasized. Simulated ozone and temperature records from the REFC1SD, GCM2000 and GCM2100 CESM1 simulations are given in the main body of the paper. Simulated results from the REFC1SD simulations for the CHASER, CMAM, MOCAGE and MRI models are given in the Supplement. In the main body of the paper all the measurements shown are from CASTNET. Additional measurements at the AQS sites are given in the Supplement. For any given simulation, all percentiles are given with respect to that particular simulation. In particular, percentiles for the future simulations are given in terms of the future distributions. Note, in addition, that all quantities shown have been rescaled following Eq. (5).

4.1 Separate evaluation of temperature and ozone

The highest rescaled average daily maximum temperatures naturally occur in the south with local maximum in the southwestern US, the midwestern region and the east coast (Figs. 2a, c, e; S1 in the Supplement). The simulations do not represent the topography with the accuracy adequate to simulate temperatures in regions of large topographic relief characteristic of the western US. Overall, when evaluated at the CASTNET sites, temperature is slightly underestimated in the REFC1SD simulations and slightly overestimated in the CESM1 GCM2000 simulations (see Table 2). The CMAM and MRI REFC1SD simulations have large positive biases in mean temperature (Fig. S1). The spatial correlation between measured and simulated rescaled temperature in the CESM1 REFC1SD and GCM2000 simulations is between 0.57 and 0.53, respectively. In the GCM2100 simulation, rescaled maximum daily temperature increases by 2.43 ^∘C on average over the US compared with the GCM2000 simulation (Fig. 2c, e and Table 2), where the regions of high temperatures in the GCM2100 simulation expand prominently with the highest temperatures extending throughout the Midwest.

Table 3 Rescaled data. The 20-year return levels for MDA8 ozone (ppb) and daily maximum temperature (^∘C) (first and third columns) at the CASTNET sites and for the CESM1 REFC1SD simulation, the GCM2000 simulation and the GCM2100 simulation. The models are sampled only at the CASTNET stations. The 20-year return levels for ozone and temperature minus their averages (second and fourth columns).

In all simulations the width of the high end of the maximum daily temperature distribution, calculated as the difference between the 90th percentile and the average maximum daily temperatures (i.e., $\mathrm{mean}\left(T\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}(T)$), maximizes in the northern part of the domain, but with an extension of the high temperature differences extending southwards through the Midwest (Figs. 2b, d, f; S1b, d). Both the REFC1SD and GCM2000 simulations underestimate $\mathrm{mean}\left(T\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}(T)$, whereas the CMAM and MRI simulations are relatively unbiased (Fig. S1b, d). There is some evidence of a similar pattern to that simulated in the CASTNET measurements and AQS measurements. Overall, the conditional maximum temperature differences show little response to climate change (e.g., compare Fig. 2d and f), suggesting the high end of the future temperature distribution does not change markedly with respect to the mean. This is consistent with the historical changes in the temperature distributions (McKinnon et al., 2016).

In all simulations, rescaled MDA8 ozone is highest in the southwestern US and in the middle Atlantic regions extending towards the central Midwest (Figs. 3a, c, e; S3a, c). The westward extent of this ozone maximum is not reflected in the CASTNET data. Consistent with many general circulation models (GCMs) (e.g., Lamarque et al., 2012; Rieder et al., 2015), all the simulations have high ozone biases (Figs. 3a, c; S2a, c; S3a, c; Table 2). Averaged over all CASTNET stations, simulated surface ozone is biased high by approximately 12 ppb in the CESM1 REFC1SD simulations and 21 ppb in the GCM2000 simulation. The spatial correlation between measured and simulated ozone in the CESM1 REFC1SD simulation and the GCM2000 simulation is 0.24 and 0.23, respectively (with a p value at 0.25 and 0.26, respectively, for the alternative hypothesis of the correlation not being 0). In the GCM2100 simulation, ozone increases by approximately 2 ppb averaged over the US with respect to the GCM2000 simulation (Fig. 3e and Table 2).

Despite the simulated positive bias in average ozone, the simulated difference between the 90th percentile and average MDA8 ozone is biased low in the CESM1 simulations (Fig. 3b, d) with average biases of −0.79 and −4.28 ppb in the CESM1 REFC1SD and GCM2000 simulations, respectively. Thus the CESM1 simulations underestimate the width of the high end of the ozone distribution. Of the other REFC1SD simulations examined, only the MOCAGE simulation shows a high bias in the width of the high end of the MDA8 ozone distribution (Fig. S3d). In all simulations except CMAM in the southeastern US the overall simulated pattern is similar to the CASTNET measurements with the largest differences in the eastern part of the domain. The geographic pattern for the high-end width of the MDA8 ozone distribution ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|{\mathrm{O}}_{\mathrm{3}}>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$) is significantly different from the equivalent quantity for temperature. While the width of the maximum temperature distribution ($\mathrm{mean}\left(T\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}(T)$) maximizes in the central US (Figs. 2b, d, f; S1), the width of the MDA8 ozone distribution maximizes in the eastern US (Figs. 3b, d, f; S2, S3). On average the difference between 90th percentile MDA8 ozone and average MDA8 ozone in the CESM1 simulations increases only by 0.26 ppb in the future simulation (Fig. 3f).

Figure 2 Rescaled data. Average daily maximum temperature (^∘C) (a, c, e), and average daily maximum temperature (^∘C) conditioned on maximum temperature greater than the 90th percentile minus average daily maximum temperature (b, d, f) for the CESM1 REFC1SD simulation (1992–2010) (a, b), the GCM2000 simulation (2006–2025) (c, d) and the GCM2100 simulation (2106–2125) (e, f). CASTNET measurements (1992–2011) of each quantity are shown as filled diamonds in the first two rows. In the first two rows we also give the average bias as the model average minus the CASTNET average for each quantity, and the correlation as the spatial correlation between the model and the CASTNET measurements. In the last row we give the difference as the mean difference between GCM2100 and GCM2000 over the continental area between 21–51^∘ N and 230–300^∘ E and the correlation as the correlation between GCM2100 and GCM2000 over the continental area between 21–51^∘ N and 230–300^∘ E. The boxes in (a) show the division of the country into various regions: the northeast, the southeast, the Midwest and the west.

Figure 3 Rescaled data. Average MDA8 ozone (ppb) (a, c, e) and average MDA8 ozone (ppb) conditioned on MDA8 ozone greater than the 90th percentile minus average MDA8 ozone (b, d, f) for the CESM1 REFC1SD simulation (1992–2010) (a, b), the GCM2000 simulation (2006–2025) (c, d) and the GCM2100 simulation (2106–2125) (e, f). The biases, differences and correlations are defined similarly to in Fig. 2.

The relative difference between changes in extreme values and the change in median values in the future simulation compared to the present-day simulation can be expressed as the quantity Ψ.

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{\Psi }(X,Y):={\displaystyle }& {\displaystyle \frac{\text{Mean\_GCM2100}\left(X\right|Y>\mathrm{90}\mathit{\%})-\text{Mean\_GCM2100}(X|\mathrm{45}\mathit{\%}<Y<\mathrm{55}\mathit{\%})}{\text{Mean\_GCM2000}\left(X\right|Y>\mathrm{90}\mathit{\%})-\text{Mean\_GCM2000}(X|\mathrm{45}\mathit{\%}<Y<\mathrm{55}\mathit{\%})}},\end{array}$$

where X,Y are ozone or temperature. If the change in extreme increments in X given Y in the GCM2100 and GCM2000 simulations are the same, we expect the ratio to be 1.

The high-end width of the future maximum daily temperature distribution is projected to increase relative to the present-day temperature distribution by up to 30 % in the southeastern US, extending northwards through the eastern Midwest (see Fig. 4b). In contrast, Ψ(O₃,O₃) is less than 1 over much of the domain (Fig. 4a). Note, however, that the region where Ψ(O₃,O₃) is slightly greater than 1, extending from the southeastern US northwestward to the Midwest, corresponds quite closely to where Ψ(T,T) is greater than 1. The overall correlation between these quantities is 0.3, significant but weak. There have been varying predictions for whether future ozone increases in the extreme (e.g., Sun et al., 2017). Figure 4a suggests that in only a few locations in a future climate does the 90th percentile ozone concentration increase by at least 10 % over the increase in the median.

As an alternative way of viewing the data we also present the 20-year return levels to describe the marginal extremes (Fig. 5). For a stationary independent series the return level is simply related to the value at a given percentile of the distribution. We note, however, that the 20-year return level represents a value considerably further out on the high end of the distribution than the 90th percentile (compare Fig. 2b, d, f with Fig. 5b, d, f for temperature and Fig. 3b, d, f and Fig. 5a, c, e for ozone). Return levels are calculated using the procedure given in Phalitnonkiat et al. (2016) (see Appendix A for more detail). Differences between the 20-year return ozone MDA8 concentration and the mean concentration are generally higher in the eastern US than in the western US (Fig. 5), consistent with Fig. 3b, d, and f. This difference is underestimated in both the CESM1 REFC1SD and GCM2000 simulations, although more dramatically so in GCM2000 and particularly in the mid-Atlantic region (Table 3). Overall, the difference between the simulated 20-year return period ozone and the mean ozone are biased low by approximately 3.3 ppb in the CESM1 REFC1SD simulation and 9.4 ppb in the GCM2000 simulation (Fig. 5 and Table 3). This underestimation suggests that the simulations do not capture the width of the high-end MDA8 ozone distribution as measured by the 20-year return period minus the mean. Note that while the simulations underestimate the differences between the 20-year return MDA8 ozone concentration and the mean concentration, the 20-year return levels are biased high (Table 3). Of the four other REFC1SD simulations examined (CHASER, CMAM, MOCAGE, MRI simulations) all except the MOCAGE simulations underestimate the high-end tail of the ozone distribution as measured by the 20-year return period. All simulations except the MOCAGE show a relative minimum in width of the tail in the mid-Atlantic region, a minimum not captured in the measurements.

Future changes in the 20-year return period MDA8 ozone concentration and the difference between the 20-year return MDA8 concentration and the mean concentration between the GCM2100 and GCM2000 simulations are 2.6 and −0.28 ppb, respectively, as measured at the CASTNET sites. This result is consistent with Rieder et al. (2015) but is at odds with a number of studies that suggest future ozone levels will increase primarily at the high end due to the impact of climate (e.g., Wu et al., 2008). The only region in which we find an increase at the high end is the Midwest (Table 3).

Figure 4 Rescaled data. (a) Ψ(O₃,O₃); (b) Ψ(T,T). (See the definition of Ψ in Eq. 7).

Simulated differences between 20-year maximum daily return temperatures and mean temperature (Figs. 5, S4) are largest in the northern part of the domain and extend southwards through the Midwest, consistent with Figs. 2 and S1. The GCM2000 simulation generally captures the measured 20-year return maximum temperature level while the CESM1 REFC1SD is biased low by almost 3.5 ^∘C (Table 3). The CMAM simulation captures the width of the temperature distribution as measured by the 20-year return period of temperature while the MRI simulation is biased low (Fig. S4).

The maximum daily temperatures (Fig. 2) and the 20-year maximum daily temperature return levels increase in the GCM2100 simulation (Table 3); however, there are only relatively small increases in the temperature difference between the 20-year return value and the mean temperature.

4.2 Joint dependence of temperature and ozone

In this section, we examine the joint dependence of ozone and temperature in the simulations and in the data. In particular, we are interested in how high ozone events are related to high temperature events in the present and future climates. We use three measures to quantify this dependence and to compare it between the future and present climates and with the measurements: the ozone temperature correlation and conditional correlation and the metrics: Ψ and φ.

• We analyze the correlation between ozone and temperature to measure the overall linear correlation between these fields. We also analyze the correlation between MDA8 ozone and maximum daily temperature conditioned on maximum daily temperature greater than the 90th percentile to measure the relationship between ozone and temperature at higher temperatures.

• The quantity Ψ(O₃,T) measures the relative response (against the mean response) of MDA8 ozone at the 90th percentile level to daily maximum temperature at the 90th percentile level in the future versus present climate (Eq. 7).

• The quantity φ gives an explicit relationship between ozone and temperature extremes (Eq. 6).

The various simulations of the MDA8 temperature and maximum daily ozone correlation show some similarities but also distinct geographical differences (Figs. 6a, c, e; S6). All simulations have a region of low correlation within the middle of the country. The three REFCS1D simulations with high-frequency temperature and ozone output (the CESM1, CMAM and MRI simulations) show this region extending inland from the Gulf of Mexico, although the CESM1 REFC1SD simulation displaces this region of low correlations further to the east than the other two. The GCM2000 and the GCM2100 simulations displace the region of low correlations further to the north without the obvious connection to the Gulf of Mexico. Differences in the correlations among the simulations are also apparent in the western and eastern thirds of the country. The MRI and CMAM simulations show strong positive correlations in the northeastern US extending westward and southward, while the CESM1 simulations have weaker correlations throughout the east. All three simulations using CESM1 have a correlation maximum over the southeastern states, with the GCM2000 and GCM2100 simulations showing a relative minimum over the northeastern states. In contrast, the CMAM and MRI simulations (Fig. S6) show a maximum correlation over the northeastern states extending to the northwest. All simulations show a band of high correlations over the western states, with all simulations except for the CMAM simulation showing regionally high correlations over the Rockies. Based on the rather sparse CASTNET and AQS measurements, it is difficult to determine which simulation better captures the true correlation pattern.

Figure 5 Rescaled data. The 20-year return level minus average MDA8 ozone (ppb) (a, c, e) and 20-year return level minus average daily maximum temperature (^∘C) (b, d, f) for the CESM1 REFC1SD simulation (1992–2010) (a, b), the GCM2000 simulation (2006–2025) (c, d) and the GCM2100 simulation (2106–2125) (e, f). The 20-year return levels from CASTNET measurements (1992–2011) for each quantity are shown as filled diamonds in the first two rows. The bias and correlation in the first two rows and the difference and correlation in the last row are defined as in Fig. 2.

The conditional correlations between MDA8 ozone and maximum daily temperature when maximum daily temperature is greater than the 90th percentile are significantly reduced across the country in comparison with the unconditional correlations. Measured conditional correlations are, in all cases, marginally positive or negative. The simulated conditional correlations in the CESM1 are somewhat higher than measured, with a maximum in the Gulf of Mexico coast states. The conditional correlations in the CMAM and MRI simulations are distinctly lower than in CESM1 (Figs. 6b, d, f; S6). Shen et al. (2016) show a suppression of ozone at high temperatures at many sites across the US.

Figure 6 Normalized data. Unconditional correlations between maximum daily temperature and MDA8 ozone (a, c, e); correlations between maximum daily temperature and MDA8 ozone conditional on maximum daily temperature greater than the 90th percentile (b, d, f), for the CESM1 REFC1SD simulation (1992–2010) (a, b), the GCM2000 simulation (2006–2025) (c, d) and the GCM2100 simulation (2106–2125) (e, f). The unconditional and conditional correlations from CASTNET measurements (1992–2011) are shown as filled diamonds in the first two rows. The black dots in the right panels indicate the significant changes from the unconditional to conditional correlations.

A metric for the response of MDA8 ozone to high temperatures can be defined as MDA8 ozone (ppb) conditioned on daily maximum temperature greater than the 90th percentile minus average MDA8 ozone (Fig. 7a, b, c). While the geographic extent of the measurements is somewhat limited, the measured response to this metric appears to be high in an arc extending from the northeastern US along the eastern seaboard into the southeastern US. The southeastern US is also a region where the right-hand side (rhs) of the temperature distribution is rather narrow (see Figs. 2b, 5b). This suggests that ozone in the southeastern US is particularly sensitive to comparatively small changes in temperature. Note, however, the results do seem somewhat at odds with those of Shen et al. (2016), who find that temperature in the southeast does not improve their statistical model of ozone exceedances.

The MRI simulation (see Fig. S7b) captures this measured pattern the best of all the model simulations with a model–measurement correlation coefficient of 0.83. The CESM1 simulations miss the high response over the northeastern US: the largest simulated response in the CESM1 simulations extends off the eastern seaboard into the southeastern US (Fig. 7a, b, c) but does not extend to the northeastern US itself. In contrast, the CMAM simulation (see Fig. S7a) shows a high sensitivity of ozone to temperature extremes over the northeastern US but misses the extension of the response along the eastern seaboard of the US. Note that the CESM1 REFC1SD and the GCM2000 simulations show similar responses in the southeast even though the CESM1 REFC1SD simulation includes interactive isoprene emissions, while the GCM2000 simulations do not.

Figure 7 Rescaled data. Average MDA8 ozone (ppb) conditioned on daily maximum temperature greater than the 90th percentile minus average MDA8 ozone for (a) the CESM1 REFC1SD simulation (1992–2010), (b) the GCM2000 simulation (2006–2025) and (c) the GCM2100 simulation (2106–2125). The biases, differences and correlations are defined similarly to in Fig. 2. (d) Ψ(O₃,T) (See the definition of Ψ in Eq. 7).

Averaged over the continental US, MDA8 ozone conditioned on a maximum daily temperature greater than the 90th percentile minus average MDA8 ozone (Fig. 7c) decreases modestly by 0.26 ppb between GCM2000 and GCM2100. Given the 2.07 ppb future increase in mean ozone (see Fig. 3), this implies ozone conditioned on the 90th percentile of mean daily maximum temperature increases by 1.81 ppb. This is consistent with a suppression of ozone at high temperatures at many sites across the US. The comparative sensitivity of ozone to temperature increases in GCM2100 versus GCM2000 can be assessed with Ψ(O₃,T) (Eq. 7). While most of the country shows future decreases in temperature sensitivity, a number of regions, including the Gulf of Mexico coast states and the Pacific Northwest, show an increase in sensitivity by over 30 % (Fig. 7d). Both of these regions also show an increase in Ψ(T,T) (Fig. 4b).

Measured and simulated scatter plots of normalized MDA8 ozone versus maximum daily temperature are shown for three CASTNET sites in Fig. 8. In each plot, the extreme points after the normalization by the ranks method (see Appendix B for the detailed procedure) are shown in red. As described above (Sect. 3) the extremal dependence between the two variables is characterized by φ, where φ gives the proportion of the extreme points at which both variables are simultaneously extreme. These sites are selected to show a range of behavior in measured and simulated φ: at one site measured φ is larger than simulated φ (Ashland, ME), at one site it is less than that simulated (Sand Mountain, AL) and at one site the measured and simulated values are about the same (Beaufort, NC). At Ashland, Maine (Fig. 8a), the CESM1 REFC1SD, the GCM2000 and the GCM2100 simulations underestimate the extreme dependence of ozone on temperature, where about 25 % of the measured points have simultaneous ozone and temperature extremes; at Sand Mountain, Alabama (Fig. 8b), the model simulations overestimate the extreme dependence, where approximately 9 % of the measured data have simultaneous ozone and temperature extremes; at Beaufort, North Carolina (Fig. 8c), about 20 % of the simulated and measured extremes occur simultaneously for temperature and ozone.

The measured sites where ozone and temperature extremes tend to co-occur (Fig. 9) in the northeastern US and in the southeastern US are related to those sites where ozone shows the most response to high temperatures (Fig. 7a, b). Previous studies have used different methodologies to capture the extremal dependence in measurements between ozone and temperature (Sun et al., 2017; Schnell and Prather, 2017; Zhang et al., 2017). Sun et al. (2017) find that the conditional probability of a high ozone day given a high temperature day is approximately 50 % in the northeastern US and 30 % in the southeast and mid-Atlantic regions. Schnell and Prather (2017) also find that the co-occurrence of temperature and ozone extremes maximizes over the northeastern US (occurring 50 % or more of the time in their analysis) but decrease towards the Midwest (where joint occurrences occur 25 % or less of the time). They also find a secondary maximum of less amplitude in the joint occurrence of extremes over the southeastern US consistent with our analysis. However, Schnell and Prather (2017) do not find the spine of low co-occurrences clearly seen in the CASTNET data from northern Alabama to Pennsylvania (also see Fig. S7f for the AQS measured data) (Zhang et al., 2017). Zhang et al. (2017) find the co-occurrence of extreme ozone and temperature occurs 32 % of the time averaged over the US with a maximum over the northeastern US during JJA and indications of a possible secondary maximum over the southeastern US.

None of the simulations using CESM1 capture the measured high co-occurrences of ozone and temperature extremes in the northeastern US. The CMAM and MRI simulations (Fig. S7b, d) do better in this regard, although the CMAM simulation does not capture the maximum in the southeastern US. There are also discrepancies among the simulations in the midwestern US. Student's t test suggests that GCM2000 fails to capture the extreme dependence between temperature and ozone at the 95 % level; however, in the CESM1 REFC1SD simulation we cannot reject with a 95 % confidence interval the null hypothesis that the simulated and measured φ values are the same. Consistent with measurements, all simulations using CESM1 (including the GCM2100 simulation) show the maximum co-occurrence of temperature and ozone maximum in the southeastern US. It is in this region that the co-occurrence of ozone and temperature maxima increases in the future (Fig. 9d).

Figure 8 Normalized data. The scatter plots of temperature and ozone from selected CASTNET sites (or the corresponding grid point) for the CESM1 REFC1SD, the GCM200 or the GCM2100 simulation. (a) Ashland (ME), (b) Sand Mountain (AL), (c) Beaufort (NC). Extreme points are shown in red if the transformed points by using the ranks method are outside the unit circle (see the ranks method in Appendix B).

5 Discussion

All the CCMI REFC1SD simulations (the CESM1 REFCS1D, CHASER, CMAM, MOCAGE and MRI) and the GCM2000 simulation show a fundamental mismatch between the locations where the width of the rhs of the temperature distribution is large and those locations where rhs of the ozone distribution is large. As measured by the difference between the 20-year return temperature and the mean temperature (Figs. 5, S4) or the difference between the 90th percentile temperature and the mean temperature (Figs. 2, S1), the width of the rhs of the temperature distribution is highest in the northern portion of the domain with a southward extension through the midwestern states and into the northwestern states. This pattern is consistent with increased temperature variability at higher latitudes (e.g., Deser et al., 2012, and references therein) and a higher temperature variance in the interior of the country due to its greater continentality. However, ozone is most sensitive to temperature changes, as measured by the slope of ozone versus temperature, where ozone precursor emissions are large (Pusede et al., 2015). This is consistent with the fact that the width of the rhs of the ozone distribution is widest (Figs. 5, 3, and S2–S5) in the eastern third of the country, where emissions of ozone precursors are generally the largest. Ozone is also most sensitive to high temperatures in the eastern part of the US (Figs. 7 and S7). Geographical differences between the shape of the ozone and temperature distributions over the US impact the relationship between ozone and temperature extremes.

The response of ozone to changes in temperature is in part determined by the temperature–ozone correlation. Details of the temperature–ozone correlation are different in all simulations (see Figs. 6 and S6). Important differences include the location of the region of low correlations in the south central part of the US, the relative strength of the correlation in the northeastern and southeastern US and the pattern of correlations in the western US. Thus we might expect the ozone response to temperature extremes to differ in the different simulations. The CASTNET and AQS measurements (Figs. 6 and S6) generally support high temperature–ozone correlations in the northeastern US, low correlations along the Gulf of Mexico coast and generally high correlations in the Rockies and west coast states (also see Shen et al., 2016).

Some of the differences between the CESM1 REFC1SD and GCM2000 simulations are likely due to meteorological differences: while the CESM1 REFC1SD simulation is driven with analyzed meteorology, the GCM2000 simulation is driven with model-calculated meteorology. In general, the CESM1 REFC1SD simulation captures the measured relation between ozone and temperature better than the GCM2000 simulation over the northeastern US, although it does not fully capture their strong measured correlation (e.g., see Fig. 6). In the southeast the measured response appears to be generally well simulated in both simulations. Overall the GCM2000 fails to capture the extreme dependence as measured by φ between temperature and ozone at the 95 % level; in the CESM1 REFC1SD simulation we cannot reject the null hypothesis that the simulated and measured values are the same.

Figure 9 Normalized data. φ from (a) CESM1 REFC1SD simulation (1992–2010), (b) GCM2000 simulation (2006–2025) and (c) GCM2100 simulation (2106–2125). φ from CASTNET measurements (1992–2011) are shown as filled diamonds in (a) and (b). The bias and correlation in (a) and (b) and the difference and correlation in (c) are defined as in Fig. 2. (d) φ from the GCM2100 simulation (2106–2125) minus φ from the GCM2000 simulation (2006–2025).

A poor simulation of the Bermuda high in the CGM2000 simulation may be important in explaining some of the differences between the CESM1 REFC1SD and GCM2000 simulations. The position of the Bermuda high strongly impacts ozone distribution over the US (Zhu and Liang, 2013), with the second empirical orthogonal function of ozone variability strongly correlated with the location of the Bermuda high (Shen et al., 2016). A westward extension of the Bermuda high is correlated with high temperatures and low ozone over much of the southeastern US (Zhu and Liang, 2013) consistent with the low correlation between ozone and temperature in all the REFC1SD simulations extending northward from the Gulf of Mexico (Figs. 6, S6). We note this region of low correlations in the CMAM and MRI REFC1SD simulations is to the west of that simulated in CESM1 REFC1SD. In the GCM2000 simulation the Bermuda high is simulated too far to the west (not shown). Consistent with this, a maximum covariance analysis shows that the mode of variability associated with the Bermuda high is also displaced too far to the west (not shown). Thus, it is likely the pattern of variability associated with the Bermuda high is incorrectly simulated in the GCM2000 simulation. In particular, the GCM simulations do not show a region of low correlation extending northward from the Gulf of Mexico, but instead a region of low correlation is situated well to the north over Kansas. Zhu and Liang (2013) show the Bermuda high is not well simulated in the majority of GCMs. This has important implications for the simulation of the ozone response to temperature over large sections of the country.

In this study, we introduce a new spectral method using multivariate extreme value theory to measure extremal dependence between temperature and ozone in both the observations and model simulations. We find, through the use of this new metric, joint extremes of temperature and ozone occur together up to approximately 35 % of the time in a few regions, although on average their joint occurrence is significantly less. Previous studies have used different methodologies to capture the extremal dependence in measurements between ozone and temperature (Sun et al., 2017; Schnell and Prather, 2017; Zhang et al., 2017). The analysis here uses a somewhat different methodology so it cannot be compared quantitatively with the previous results, but qualitatively the overall patterns are similar to those found previously. The advantage of the spectral method for finding joint extremes of temperature and ozone is that it gives detailed information about the joint extremes and is not restricted to a particular quantile of the distribution. It can be used to forecast joint extremes even out of the range of available samples.

The various model simulations differ in their simulation of φ (measuring the joint spectral extremes of ozone and temperature), again suggesting ozone may respond differently to high temperatures in the different simulations. In the central part of the country where the rhs of the temperature distribution is particularly wide, φ is relatively small (Figs. 9, S7) in all but the REFC1SD CESM1 simulation (Figs. 9, S7). φ is high in the CMAM and the MRI simulations in the northeastern US but not in the CESM1 REFC1SD or GCM2000 simulations; conversely, in the CESM1 REFC1SD, the GCM2000 and the MRI simulations φ is high in the southeastern US but not in the CMAM simulation. φ is high in the majority of the simulations in the northwestern states.

In general (with some exceptions along the US west coast) the geographical pattern of averaged MDA8 ozone conditioned on daily maximum temperature greater than the 90th percentile $\left(\mathrm{mean}\right({\mathrm{O}}_{\mathrm{3}}|T>\mathrm{90}\mathit{\%})-\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right))$ is qualitatively similar to that of φ in the model simulations and the measurements. However, it is important to note that while the measurement sites used here are sufficiently dense in the eastern US to resolve some of the regional features, they are nowhere dense enough to resolve regional features in the western two-thirds of the country. The CASTNET measurements show the northeast region and some sites in the southeast have the largest response of ozone to temperature extremes (Fig. 7). These same regions have the highest spectral dependence between ozone and temperature extremes (Fig. 9) and the highest measured correlations between temperature and ozone conditioned on temperature greater than the 90th percentile level (Fig. 6). Note that the CASTNET measurement sites in the far southeastern part of the country have comparatively small measured variations in relatively extreme temperature $\left(\mathrm{mean}\right(T|T>\mathrm{90}\mathit{\%})-\mathrm{mean}\left(T\right))$. Thus, at some sites over the southeastern US the sensitivity of ozone to changes in temperature is relatively large. In contrast, over the northeastern US both the relative variations in the width of the rhs of the MDA8 ozone distribution ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|{\mathrm{O}}_{\mathrm{3}}>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$) and the rhs of the daily maximum temperature ($\mathrm{mean}\left(T\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}(T)$) are relatively high and thus the ozone sensitivity to temperature change is relatively small.

Comparing the GCM2100 and the GCM2000 simulations, the mean future temperature increases everywhere in relation to the current climate, although in some places the width of the rhs of the future temperature distribution decreases. The future difference between summertime maximum daily temperatures at the 90th percentile minus mean maximum daily temperature increases by up to 20 to 30 % compared to the present day (Ψ(T,T)) over parts of the southern Mississippi basin extending to the northern Midwest and the northwest coast (Fig. 4). The pattern of this increase bears a striking resemblance to those locations where measured summertime inter-annual extreme temperatures increase relative to inter-annual increases in the mean (Huybers et al., 2014), a process linked to drying of the soils (Huybers et al., 2014). In other parts of the country the relative future increase is small or negative. Note that the increase in (Ψ(T,T)) over the lower Mississippi valley occurs in a location strongly impacted by the Bermuda high and thus warrants further investigation. Note also, this increase occurs in those locations where the CESM1 REFC1SD simulation gives a low correlation between ozone and temperature.

To what extent are the relative changes in the future width of the MDA8 ozone distribution determined by the relative future changes in that of the maximum daily temperature distribution? In most locations the future width of the ozone distribution decreases. An exception is in the Midwest where the increase in the future width of the maximum daily temperature distribution is most pronounced. The spatial correlation between Ψ(T,T) and Ψ(O₃,O₃) is significant but weak, with a correlation coefficient of 0.3. This suggests a weak relationship between changes in the rhs of the future MDA8 ozone distribution and the future maximum daily temperature distribution. Overall, future ozone is less responsive to temperature than present-day ozone ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$) but the effect is small averaged over the continental US (−0.26 ppb), suggesting only relatively modest temperature suppression. The ratio of future sensitivities to temperature compared to the present varies regionally, ranging from −50 % to +50 % (Fig. 7). Interestingly, ozone does become more responsive to temperature changes in the lower Mississippi valley (Fig. 7), in precisely the region in which the width of the temperature distribution increases (Fig. 4). However, even in this region the rhs of the future ozone distribution does not become significantly wider than its present-day values (Fig. 4).

There have been different predictions as to whether climate change increases future ozone extremes with respect to the increase in the mean (e.g., see Sun et al., 2017). On average the width of the rhs of the future MDA8 ozone distribution increases slightly by 0.26 ppb in the future simulation ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|{\mathrm{O}}_{\mathrm{3}}>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$) (Fig. 3, but also see Fig. 5 for the relative change in the 20-year return level), but in many locations the relative width decreases (Ψ(O₃,O₃), Fig. 4). Where it increases, the increase is always less than 20 %. Our results generally suggest that the increase in future ozone is primarily due to a shift in the ozone distribution and not due to an increase in ozone at the high end.

6 Conclusions

We investigate high temperature and ozone extremes and their joint occurrence over the United States during the summer months (JJA) in measurements and simulations of the present and future climate. Three simulations using the CESM1 with chemistry were analyzed: the CESM1 CCMI reference experiment using specified dynamics (REFC1SD) between 1992 and 2010, a 25-year present-day simulation branched off of the CCMI REFC2 simulation in the year 2000 (GCM2000) and a 25-year future simulation branched off of the CCMI REFC2 simulation in 2100 (GCM2100). Distinct from the CCMI REFC2 simulations, the emissions and long-lived greenhouse gas distributions (except CO₂) are held constant in the GCM2000 and GCM2100 simulations at values representative of the year 2000. In addition, we analyzed the REFC1SD simulation in four additional models with data available at a sufficiently high temporal frequency: the CHASER, CMAM, MOCAGE and MRI models. All the CCMI REFC1SD simulations (CESM1 REFCS1D, CHASER, CMAM, MOCAGE and MRI) and the GCM2000 simulation have a large bias in maximum daily ozone. Scaled ozone biases are 12 and 21 ppb, respectively, in the CESM1 REFC1SD and GCM2000 simulations. Consistent with many global model simulations, the ozone bias is particularly pronounced over the eastern US. The simulation of daily maximum temperatures shows considerable variability among the various model simulations.

The average global mean daily maximum temperature change between the present-day simulation (GCM2000) and the future simulation (GCM2100) is 2.1 ^∘C, less than the 2.8 ^∘C difference in the parent CCMI REFC2 simulations. The difference between these sets of simulations is most likely attributable to the fact the GCM2100 simulation includes the effect of increased CO₂ forcing in the future but does not account for the impact of projected future aerosol decreases. Thus, in the GCM2100 simulation the relatively large aerosol radiative forcing acts as a buffer against the increased CO₂. Over the continental US ozone increases by approximately 2.1 ppb between the GCM2000 simulation and the GCM2100 simulation.

The main conclusions from this study are as follows.

• Five out of six of the simulations analyzed underestimate the width of the measured tail at the high end of the ozone distribution in the present climate, despite the fact that all simulations overestimate the mean ozone. The 20-year return period of ozone minus its mean is underestimated by more than 9 ppb in the GCM2000 simulation evaluated over all CASTNET sites, while in the REFC1SD CESM1 simulation it is underestimated by somewhat more than 3 ppb. The 20-year return period of temperature minus its mean is generally about 2 ^∘C less than measured in both the GCM2000 and the CESM1 REFC1SD simulations. Despite large biases in mean daily maximum temperature the bias in the width of the rhs of the temperature tails in the CMAM and MRI simulations $\left(\mathrm{mean}\right(T|T>\mathrm{90}\mathit{\%})-\mathrm{mean}\left(T\right))$ is less than 1 ^∘C.

• We propose a new method to measure the joint extremes of temperature and ozone by calculating the spectral density (φ) of the joint extremes of ozone and temperature. This measure of the joint extremes is not restricted to a particular quantile of the distribution but can be used to forecast joint extremes even out of the range of available samples. While in many areas of the country MDA8 ozone and maximum daily temperature are highly correlated, the correlation is reduced significantly at the higher end of the distributions. Measures of spectral density are less than about 0.35 everywhere, so that only about a third of the time do extreme temperatures coincide with extremely high ozone. Observations show that φ is highest in the northeastern US and in the southeastern US. To some extent this response is consistent with the ozone response to extreme temperatures ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|T>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$). To the extent that the measurements are dense enough to define the spectral density geographically, the simulations capture much of the measured pattern.

• In all simulations there is a geographical mismatch between where the rhs of the simulated maximum daily temperature distribution is large and where the rhs of the simulated MDA8 ozone distribution is large. Thus, while ozone concentrations are often correlated with temperature, the regions of high ozone extremes do not necessarily match the regions of high temperature extremes. All things being equal, we might expect the ozone distribution to be a slave to the temperature distribution so that regions of particularly high temperature extremes might also be expected to have particularly high ozone extremes. However, this is not the case. The highest temperature extremes tend to occur in the Midwest while the highest ozone extremes tend to occur in the eastern part of the country. Regions with high ozone precursor emissions are known to increase the ozone–temperature slope, making ozone in regions with high precursor emissions sensitive to smaller temperature variations. Other complicating factors such as the importance of biogenic emissions or regional meteorological differences may also complicate the distributional relation between ozone and temperature.

• The various model simulations show some rather pronounced differences in the ozone–temperature relationship. These differences suggest that ozone will respond rather differently to temperature changes in the various simulations. Differences between the CESM1 REFC1SD and the GCM2000 simulations can be attributed in part to differences in the meteorology. The response of the REFC1SD simulation is qualitatively better than that of the GCM2000 simulation. We hypothesize that the differences in these simulations are meteorologically induced and may, at least in part, be attributed to a poor simulation of the Bermuda high in the GCM2000 simulation.

• In the future climate the ozone and temperature distributions shift to the right. Our results generally suggest that the increase in both future temperature and future ozone is primarily due to a shift in the distributions, not to an increase in the extremes. Overall, the rhs of the temperature distribution increases slightly, with the largest increase, when evaluated at the CASTNET sites, in the Midwest. In some locations the increase in the rhs of the temperature distribution approaches 30 %, in other locations it decreases up to 10 %. The pattern of increase is what might be expected from soil–moisture feedbacks. On average the width of the rhs of the future ozone distribution increases slightly by 0.26 ppb in the future simulation ($\mathrm{mean}\left({\mathrm{O}}_{\mathrm{3}}\right|{\mathrm{O}}_{\mathrm{3}}>\mathrm{90}\mathit{\%})-\mathrm{mean}({\mathrm{O}}_{\mathrm{3}})$), with regional increases of up to 20 % and decreases of up to 10 % over parts of the northeastern US and much of the western two-thirds of the country. At CASTNET sites increases in the 20-year return period minus the mean are only found on average in the Midwest. The correlation between relative changes in the high end of the future temperature distribution (Ψ(T,T)) and the ozone distribution (Ψ(O₃,O₃)) is 0.3, relatively small but still significant. Thus an increase in the rhs of the future temperature distribution may have some impact on the future ozone distribution. However, the correlation is weak, suggesting other complicating factors. It is possible that a stronger relationship would emerge from longer model simulations. In any case the region where the rhs of the temperature distribution increases does not correspond to the region where the width of the rhs of the ozone distribution is large. Elsewhere in the world, perhaps in regions with strong soil–moisture feedbacks and high emissions of ozone precursors, a future amplification in the future temperature distribution would have more dramatic impacts on future ozone extremes.

Appendix A: Univariate regular variation

To understand the basic characteristics of extreme distributions, we should introduce the notion of regular variation. A regularly varying function is a function whose behavior at infinity follows a power-law function. That is, a regularly varying function with an index α can be explained by

$$\begin{array}{}\text{(A1)}& \underset{t\to \mathrm{\infty }}{lim}{\displaystyle \frac{F\left(tx\right)}{F\left(t\right)}}=x^{\mathit{\alpha }}\end{array}$$

for all x>0.

Regularly varying functions are studied in many fields and one of the applications that we will use here is to estimate the tail indices α of extreme ozone and extreme temperature distributions in order to estimate N-year return levels of those variables. Alternatively, we can fit the ozone or temperature distributions to the generalized Pareto distribution and estimate the shape parameters which are equivalent to the reciprocal of the tail indices (shape $={\mathit{\alpha }}^{-\mathrm{1}}$) if the shape parameter is positive. (Phalitnonkiat et al., 2016) suggest a procedure to estimate shape parameters using a combination of Hill estimators and maximum likelihood estimators.

Appendix B: Ranks method

Let us consider two-dimensional random vectors (X,Y). When the tail part of X's distribution and the tail part of Y's distribution are independent, we would expect that X and Y are unlikely to yield extreme values at the same time, and vice versa. This observation suggests that when we plot only extreme points from (X,Y), the points would appear to be around the axes if X and Y are extremely independent, and vice versa. This is actually true in higher dimensions as well. However, the tool described above for measuring the dependence between variables only applies to variables with the same marginal tail indices.

Among the methods suggested by Resnick (2007), we use a transformation that essentially normalizes the tail indices of all components to 1 without calculating or estimating the tail indices α_j for each $j=\mathrm{1},\mathrm{\dots },d$. This method is called the ranks methods. The major benefit from this method is that we can avoid the marginal tail index estimation, which reduces numerical errors; however, the drawback is that the transformation itself destroys the iid (independent and identically distributed) property of the data and makes it more complicated to obtain asymptotic distributions; see Einmahl et al. (2001). The method can be carried out as follows. Let ${\mathit{X}}_i=(X_i^{\left(\mathrm{1}\right)},\mathrm{\dots },X_i^{\left(d\right)}),i=\mathrm{1},\mathrm{\dots },n$ be d-dimensional vectors. Denote the rank of $X_i^{\left(j\right)}$ by

$$\begin{array}{}\text{(B1)}& r_i^{\left(j\right)}:=\sum _{m=\mathrm{1}}^n{\mathrm{1}}_{[X_m^{\left(j\right)}\ge X_i^{\left(j\right)}]}.\end{array}$$

For a fixed k>0 and for each $i=\mathrm{1},\mathrm{\dots },n$ we transform X_i into a rank vector by

$$\begin{array}{}\text{(B2)}& \left(X_i^{\left(\mathrm{1}\right)},\mathrm{\dots },X_i^{\left(d\right)}\right)\mapsto \left({\displaystyle \frac{k}{r_i^{\left(\mathrm{1}\right)}}},\mathrm{\dots },{\displaystyle \frac{k}{r_i^{\left(d\right)}}}\right).\end{array}$$

We consider a point $\left(X_i^{\left(\mathrm{1}\right)},\mathrm{\dots },X_i^{\left(d\right)}\right)$ as jointly extreme if $\left|\left|\left(\frac{k}{r_i^{\left(\mathrm{1}\right)}},\mathrm{\dots },\frac{k}{r_i^{\left(d\right)}}\right)\right|\right|>\mathrm{1}$, where $\left|\right|\cdot \left|\right|$ is a norm in ℝ^d. In this case, we use the L² norm or least squares. We use the transformed vectors to estimate the spectral measure (or angular measure in the case of two-dimensional vectors).

Appendix C: Estimating spectral measure

To estimate the spectral measure from the data in two-dimensional polar coordinates, we measure the angles between the transformed points $\left(\frac{k}{r_i^{\left(\mathrm{1}\right)}},\frac{k}{r_i^{\left(\mathrm{2}\right)}}\right)$ and the x axis. That is, we can estimate the empirical measure $\widehat{S}$ by

$$\begin{array}{}\text{(C1)}& \widehat{S}\left(A\right)={\displaystyle \frac{\text{\# of extreme points with angles in}A}{\text{\# of extreme points}}},\end{array}$$

where A is a set or an interval. Note that this can be extended to higher dimensions in a similar way.

We may notice that the choice of k has a major role in how we categorize extreme points. The higher k is, the more points, and hence the more extreme points, would lie outside the unit circle. We use the procedure from Nguyen and Samorodnitsky (2013) to estimate k.

Since the angular measure is normalized (i.e., the area under curve from 0 to $\frac{\mathit{\pi }}{\mathrm{2}}$ is 1), we can only consider the area of the “middle” part, which we define to be between $\frac{\mathit{\pi }}{\mathrm{8}}$ and $\frac{\mathrm{3}\mathit{\pi }}{\mathrm{8}}$. Denote this amount by φ:

$$\begin{array}{}\text{(C2)}& \mathit{\phi }:=\widehat{S}\left(\left[{\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}},{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{8}}}\right]\right)\approx \mathrm{area}\left[{\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}},{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{8}}}\right],\end{array}$$

where the area is defined in a notion of kernel density estimation from the angular measure.