Extremal Dependence between Temperature and Ozone over the Continental U.S.

. The co-occurrence of heat waves and pollution events and the resulting high mortality rates emphasizes 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 speciﬁed dynamics (REFC1SD) 5 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 latter 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 U.S. do not occur in regions of relatively high temperature extremes. A new 10 metric, the spectral density, is developed to measure the joint

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 5 the future temperature probability distribution simply shifts to higher temperatures but does not change in shape consistent with measured trends (McKinnon et al. (2016)). On the other hand, 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 U.S. in 2050, attributing this to an increase in stagnation episodes with soil-moisture feedbacks impacting 10 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 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 This article is organized as follows: in Section 2, we describe the datasets and model simulations used; in Section 3, we introduce the statistical procedures used to quantify the relationship between ozone and temperature. In Section 4, we present 20 the results then discuss these results in Section 5. Section 6 gives the conclusions.

Data and Model Descriptions
In this study we examine the simulated and measured relationship between ozone and temperature extremes over the U.S.
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-2010. This allows a robust evaluation of simulated ozone and 25 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 are free-running simulations in that the meteorology, sea surface temperatures (SST) 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.
ing comparisons of ozone return periods. The horizontal grid resolution in the CESM1 simulations analyzed here is 1.9 • ×2.5 • ; 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-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 Heg- 5 glin and Lamarque (2015): 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-hour 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 is available. Details on the additional model simulations is given in Morgenstern et al. (2017). 10 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 CO2 concentration is specified at 369 ppm, representative of the year 2000; in the future simulation (the GCM2100 simulation) the CO2 concentration is specified at 669 ppm, representative of the 2100 concentration of CO2 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 15 are also fixed and are representative of the year 2000. Both the GCM2000 and GCM2100 simulations are 25-year simulations branched off 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 20 is likely due to the fact that the emissions of GHGs 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-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 (Beufort NC, and Lassen Volcanic 25 CA) where the first 2 years or 3 years of data were missing, respectively. See Figure 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 provides a good estimate of 10-meter ozone concen-30 trations as measured by CASTNET (Brown-Steiner et al. (2015)). The maximum daily 2-meter temperature is used in both the CASTNET measurements and the simulations.
To render the data approximately stationary on both the 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 June 1 st and d = 91 refers to August 30 th .
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 m y,a : 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,(1) ≤ x y,(2) ≤ ... ≤ x y,(D) . Then, we calculate a daily ozone deviation: 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 averagex G y,d for each day d over all years Y = 20 (or Y = 19 for REFC1SD). That is, for each day d, we calculate: M d is then smoothed by 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 bŷ where is the standard deviation of day d. Later in the text, we refer to (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 rescalex DS y,d back to its original scale (x res y,d ) by using 15 the formula: where Y = 20 (or Y = 19 for REFC1SD).

Methodology
In this study besides 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 coef-20 ficients 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 0, π 2 . If the mass of the pdf is concentrated near 0 or near π 2 , extremes of the two variables are unlikely to be significant at the same time, which 5 points to an independence of the extremes. On the other hand, if the mass of the pdf is concentrated away from the endpoints 0, π 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 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 π 8 and 3π 8 , to represent the extreme dependence between two variables. Denote this amount by ϕ: See the detailed explanation in Appendix C. Note that the range of ϕ is [0, 1], where ϕ = 1 refers to extreme dependence and extreme independence implies ϕ = 0. 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 ( Figure 1c). In contrast, Figure 1b gives an example of data with low correlation but highly dependent extremes ( Figure 1d).

Results
In this section we compare measured and simulated temperature and ozone records separately (4.1) and then their joint dependence are analyzed in section (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 MRIs 20 models are given in the supplement. In the main body of the paper all the measurements shown are from CASTNET. Addition 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 equation (5).

25
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 (Figures 2a, c, e; Figure S1). 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 [1, −0.9; −0.9, 1] follows a bivariate normal, where µ is the mean vector and Σ is the covariance matrix, and Z ∼ N (µ = 10, σ 2 = 9). The sample size is also n = 10000. The plots show the existence of tail dependence by having high angular density near π 4 (ϕ = 0.75); however, low correlation (true ρ = 0; estimated ρ = −0.05).
simulations have large positive biases in mean temperature ( Figure S1). The spatial correlation between measured and simulated rescaled temperature in the CESM1 REFC1SD and GCM2000 simulation 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 (Figures 2c, e and Table 2), where the regions of high temperatures in the GCM2100 simulation expand prominently with a tongue of the highest temperatures extending throughout the Midwest. 5 In all simulations the width of the high end maximum daily temperature distribution, calculated as the difference between the 90 th percentile and the average maximum daily temperatures (i.e., mean(T |T > 90%)−mean(T )) maximizes in the northern part of the domain, but with a tongue of high temperatures differences extending southwards through the Midwest ( Figure   2b 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 20 and the GCM2000 simulation is 0.24 and 0.23 respectively (with 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 ( Figure 3e and Table 2).
Despite the simulated positive bias in average ozone, the simulated difference between the 90 th percentile and average    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 Ψ.
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 5 temperature distribution by up to 30% in the Southeast US extending northwards through the Eastern Mid-west (see Figure 4b).
In contrast, Ψ(O 3 , O 3 ) is less than 1 over much of the domain ( Figure 4a). Note, however, that the region where is slightly greater than 1, extending from the southeast US northwestward to the Midwest corresponds quite close 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 10 locations in a future climate does the 90 th percentile ozone concentration increase by at least 10% over the increase in the median.   (7)).
As an alternative way of viewing the data we also present the 20-year return levels to describe the marginal extremes ( Figure   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 both the CESM1 REFC1SD and GCM2000 simulations, although more dramatically so in GCM2000 and particularly in the ozone are biased low by approximately 3.3 ppb in the CESM1 REFC1SD simulation and 9.4 ppb in the GCM2000 simulation ( Figure 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 (the CHASER, CMAM, MOCAGE, MRI simulations) 5 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 twenty-year return period MDA8 ozone concentration and the difference between the twenty-year return MDA8 concentration and the mean concentration between the GCM2100 and GCM2000 simulations are: 2.6 and -0.28 10 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 we find an increase at the high end is the Midwest (Table 3).
Simulated differences between 20-year maximum daily return temperatures and mean temperature ( Figure 5; Figure S4) are largest in the northern part of the domain and extend southwards through the Midwest consistent with Figure 2 and Figure S1. 15 The GCM2000 simulation generally captures the measured twenty-year return maximum temperature level while the CESM1 REFC1SD is biased low by almost 3.5 o C ( Figure 5, 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 ( Figure S4).
The GCM2100 maximum daily temperatures ( Figure 2) and 20-year maximum daily temperature return levels increase (Table 3), with a relatively small increase in the temperature difference between the twenty-year return value and the mean 20 temperature. In the Midwest this difference increases by about 3.5 o .

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: 25 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 90 th percentile to measure the relationship of ozone and temperature at higher temperatures.
-The quantity Ψ(O 3 , T ) measures the relative response (against the mean response) of MDA8 ozone at the 90 th percentile 30 level to daily maximum temperature at the 90 th percentile level in the future versus present climate (Equation (7)).
-The quantity ϕ gives an explicit relationship between ozone and temperature extremes (Equation (6)).  The various simulations between the MDA8 temperature and maximum daily ozone correlation show some similarities, but also distinct geographical differences (Figures 6a, c, f; Figure 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  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 (Figures 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 northeast U.S. along the eastern seaboard into the southeastern U.S. The southeastern U.S. is also a region where the right hand side of the temperature distribution is rather narrow (see Figures 2b,5b;Figures S1c,S4c). This suggest that 25 ozone in the southeast U.S. is particularly sensitive to comparatively small changes in temperature. Note, however, the results do seem somewhat at odds with Shen et al. (2016) who find that temperature in the Southeast does not improve their statistical model of ozone exceedances.
The MRI simulation (see Figure S7b) captures this measured pattern the best of all the model simulations with a modelmeasurement correlation coefficient of 0.83. The CESM1 simulations miss the high response over the northeast US: the largest 30 simulated response in the CESM1 simulations extends off the Eastern seaboard into the southeastern US (Figures 7a, b, c), but does not extend to the northeastern U.S. itself. In contrast, the CMAM simulation (see Figure S7a) shows a high sensitivity of ozone to temperature extremes over the northeast U.S., but misses the extension of the response along the eastern seaboard of the U.S. 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 does not.   (7)).
Averaged over the continental US, MDA8 ozone conditioned on maximum daily temperature greater than the 90 th percentile minus average MDA8 ozone (Figure 7c) decreases modestly by 0.26 ppb between GCM2000 and GCM2100. Given the 2.07 ppb future increase in mean ozone [see Figure 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 Measured and simulated scatter plots of deseasonalized MDA8 ozone versus maximum daily temperature are shown for three CASTNET sites in Figure 8. In each plot, the extreme points after the normalization by the ranks method (see Appendix B for 10 the detailed procedure) are shown in red. As described above (Section 3) the extremal dependence between the two variables is characterized by ϕ, where ϕ gives the proportion of the extreme points where 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 (Figure 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 (Figure 8b), the model simulations overestimate 5 the extreme dependence, where approximately 9% of the measured data have simultaneous ozone and temperature extremes; at Beaufort North Carolina (Figure 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 (Figure 9) in the northeastern US and in the southeastern US are related to those sites where ozone shows the most response to high temperatures (Figure 7a, b). Previous 10 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) finds that the conditional probability of a high ozone day given a high temperature day is approximately 50% in the northeast U.S. and 30% in the southeast and mid-Atlantic regions. Schnell and Prather (2017) also find that the co-occurrence of temperature and ozone extremes maximize over the northeast US (occurring 50% or more of the time in their analysis) but decrease towards the Midwest (where joint 15 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 Figure S7f  does not capture the maximum in the southeastern U.S. There are also discrepancies between the simulations in the midwestern U.S. Student's t-test suggests that the GCM2000 fails to capture the extreme dependence between temperature and ozone at the 25 95% level; however, in the CESM1 REFC1SD simulation we cannot reject with a 95% confidence interval the null hypothesis that the simulated and measured ϕ are the same. Consistent with measurements, all simulations using the 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 increase in the future (Figure 9d).

30
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 right hand side (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 ( Figure 5; Figure S4) or the difference between the 90 th percentile temperature and the mean temperature (Figure 2; Figure 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 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 its greater continentality. However, ozone is most sensitive 5 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 ( Figure 5, 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 Figure 6 and Figure S6). Important differences include the location of the region of low correlations in the south central part of the U.S., the relative strength of the correlation in the northeast and southeast U.S. and the pattern of correlations in the western U.S. Thus we might expect the ozone response to 15 temperature extremes will differ in the different simulations. The CASTNET and AQS measurements (Figure 9 and Figure   S6) generally support high temperature-ozone correlations in the Northeast U.S., low correlations along the Gulf 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 northeast US, although it does not fully capture their strong measured correlation (e.g., see Figure 6). In the Southeast the measured response appears to be generally well simulated 5 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.
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 lated with high temperatures and low ozone over the much of the Southeast 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( Figure 6; Figure S6) extending northward from the gulf. We note this region of low correlations in the CMAM and MRI REFC1SD simulations is to the west of that simulated in the 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 5 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, 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 10 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 15 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 pattern of 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. 20 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 right-hand side of the temperature distribution is particularly wide ϕ is relatively small (Figure 9; Figure   S7) in all but the REFC1SD CESM1 simulation (Figure 9; Figure S7). ϕ is high in the CMAM and the MRI simulations in the northeast U.S., but not in the CESM1 REFC1SD or the GCM2000 simulations; on the other hand in the CESM1 REFC1SD, 25 the GCM2000 and the MRI simulations ϕ is high in the southeast U.S., 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 U.S. west coast) the geographical pattern of averaged MDA8 ozone conditioned on daily maximum temperature greater than the 90th percentile (mean(O 3 |T > 90%) − mean(O 3 )) is qualitatively similar to that of ϕ in the model simulations and the measurements. However, it is important to note that while the measurement sites 30 used here are sufficiently dense in the eastern U.S. 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 (Figure 7). These same regions have the highest spectral dependence between ozone and temperature extremes (Figure 9) and the highest measured correlations between temperature and ozone conditioned on temperature greater than the 90 th percentile level ( Figure 6). Note, that the 35 relative extreme temperature (mean(T |T > 90%)−mean(T )). 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 (mean(O 3 |O 3 > 90%)−mean(O 3 )) and the rhs of the daily maximum temperature (mean(T |T > 90%) − mean(T )) are relatively high and thus the ozone sensitivity to temperature change is relatively small. 5 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 90 th 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 (Figure 4). The pattern of this increase bears a striking resemblance to those 10 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. 15 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 3 , O 3 ) is significant, but weak, with a correlation coefficient of 0.3. This suggests a weak relationship between changes in the right-hand side of the future MDA8 ozone distribution 20 and the future maximum daily temperature distribution. Overall, future ozone is less responsive to temperature than present day ozone (mean(O 3 |T > 90%) − mean(O 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 present varies regionally ranging from -50% to +50% (Figure 7). Interestingly, ozone does become more responsive to temperature changes in the lower Mississippi valley (Figure 7), in precisely the region that the width of the temperature distribution increases ( Figure   25 4). However, even in this region the rhs of the future ozone distribution does not become significantly wider than its present-day values (Figure 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 (mean(O 3 |O 3 > 90%) − mean(O 3 )) ( Figure 2, but also see Figure 5 30 for the relative change in the 20-year return level), but in many locations the relative width decreases (Ψ(O 3 , O 3 ), Figure 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.
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-2010, a 25-year present-day simulation branched off the CCMI REFC2 simulation in the year 2000 (GCM2000) and a 25-year future 5 simulation branched off the CCMI REFC2 simulation in 2100 (GCM2100). Distinct from the CCMI REFC2 simulations the emissions and long-lived greenhouse gas distributions (except CO2) 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 sufficiently high temporal frequency: the CHASER, CMAM, MOCAGE and MRI models. All the CCMI REFC1SD simulations (the CESM1 REFCS1D, CHASER, CMAM, MOCAGE and MRI) and the GCM2000 simulation The average global mean daily maximum temperature change between the present-day simulation (GCM2000) and the 15 future simulation (GCM2100) simulation is 2.1 • C, less than the 2.8 • C difference in the parent CCMI REFC2 simulations.
The difference between these set of simulations is most likely attributable to the fact the GCM2100 simulation includes the effect of increased CO2 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 CO2. Over the continental US ozone increases by approximately 2.1 ppb between the GCM2000 simulation and the GCM2100 simulation. 20 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 25 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 (mean(T |T > 90%) − mean(T )) is less than -We propose a new method to measure the joint extremes of temperature and ozone by calculating the spectral density 30 (ϕ) 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 everywhere less than about 0.35, so that only about a third of the time at most do extreme temperatures coincide with extremely high ozone. Observations show that ϕ is highest in the northeast US and in the southeast US. To some extent this response is consistent with the ozone response to extreme temperatures (mean(O 3 |T > 90%) − mean(O 3 )). To the extent 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 5 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 10 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 15 differences suggest that ozone will respond rather differently to temperature changes in the various simulations. Differences between the CESM1 REFC1SD and the GCM2000 simulation can be attributed in part due 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. 20 -In the future climate the ozone and temperatures 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 25 feedbacks. On average the width of the rhs of the future ozone distribution increases slightly by 0.26 ppb in the future simulation (mean(O 3 |O 3 > 90%) − mean(O 3 )), with regional increases up to 20% and decreases up to 10

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 30 with an index α can be explained by 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 (GPD) and 5 estimate the shape parameters which are equivalent to the reciprocal of the tail indices (shape = α −1 ) if the shape parameter is positive. Phalitnonkiat et al. (2016) suggests a procedure to estimate shape parameters using a combination of Hill estimators and Maximum Likelihood Estimators (MLE).

Appendix B: Ranks Method
Let us consider 2-dimensional random vectors (X, Y ). When the tail part of X's distribution and the tail part of Y 's distribution 10 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 extreme 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 15 all components to 1 without calculating or estimating the tail indices α j for each j = 1, ..., 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 done as follows.
For a fixed k > 0 and for each i = 1, ..., n we transform X i into a rank vector by We consider a point X > 1, where || · || is a norm in R d . In this case, we use the L 2 -norm or least squares. We use the transformed vectors to estimate the spectral measure (or angular measure in the case of 2-dimensional vectors).

Appendix C: Estimating Spectral Measure
To estimate the spectral measure from the data in 2-dimensional polar coordinates, we measure the angles between the transformed points k r (1) i , k r (2) i and the x-axis. That is, we can estimate the empirical measureŜ bŷ S(A) = # of extreme points with angles in A # of extreme points , 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 on how we categorize extreme points. The higher k is, the more points 5 would lie outside the unit circle, and hence, the more extreme points. 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 π 2 is 1), we can only consider the area of the 'middle' part, which we define to be between π 8 and 3π 8 . Denote this amount by ϕ: where the area is defined in a notion of kernel density estimation from the angular measure. Other GHG as in GCM2000.
[Rescaled Data] MDA8 ozone averages (ppb) and daily maximum temperature averages ( • C) in different regions over the U.S. 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 deviations (sd) are calculated between 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)