Changes of Anthropogenic Precursor Emissions Drive Shifts of Ozone Seasonal Cycle throughout Northern Midlatitude Troposphere

. Simulations by six CMIP6 Earth System Models indicate that the seasonal cycle of baseline tropospheric ozone at northern midlatitudes has been shifting since the mid-20 th Century. Beginning in ~1940 the magnitude of the seasonal cycle 20 increased by ~10 ppb (measured from seasonal minimum to maximum), and the seasonal maximum shifted to later in the year by about 3 weeks. This shift maximized in the mid-1980s, followed by a reversal - the seasonal cycle decreased in amplitude and the maximum shifted back to earlier in the year. Similar changes are seen in measurements collected from the 1970s to the present. The timing of the seasonal cycle changes is generally concurrent with the rise and fall of anthropogenic emissions that followed industrialization and subsequent implementation of air quality emission controls. A quantitative 25 comparison of the temporal changes of the ozone seasonal cycle at sites in both Europe and North America with the temporal changes of ozone precursor emissions across the northern midlatitudes and find a high degree of similarity between these two temporal patterns. We hypothesize that changing precursor emissions are responsible for the shift in the ozone seasonal cycle; this is supported by the absence of such seasonal shifts in southern mid-latitudes where anthropogenic emissions are much smaller. We also suggest a mechanism by which changing emissions drive the changing seasonal cycle: increasing allow summertime photochemical production of ozone to become more important than ozone transported the stratosphere and increasing VOCs lead to progressively greater photochemical ozone production in the summer months, thereby increasing the amplitude of the seasonal ozone cycle. Decreasing emissions of both precursor classes then reverse these changes. The quantitative parameter values that characterize the seasonal shifts provide useful benchmarks for evaluating model simulations, both against observations and between models.

3 mentioning this shift. For example, studies have documented increasing springtime ozone (Lin et al., 2015), often in combination with decreasing summer ozone (Chan et al., 2010;Cooper et al., 2012;Lin et al., 2017), which indicates that the 70 seasonal cycle is shifting towards a springtime maximum.
Some studies discuss another aspect of the varying seasonal ozone cycle: changes in its amplitude. This finding is most evident in measured data collected across the US (Simon et al., 2014), and specifically in the eastern US (Strode et al., 2015). Other studies do not explicitly mention the changing amplitude but still provide evidence for this phenomenon.
According to measurements and models, these papers find ozone decreasing in summer, when it has typically been highest 75 (Hogrefe et al., 2011;Cooper et al., 2012;Parrish et al., 2012;2013 and references therein;Simon et al., 2014;Lin et al., 2017), and concurrently increasing in the winter, when it has typically been lowest (Bloomer et al., 2010;Chan et al., 2010;Cooper et al., 2012;Parrish et al., 2012;Lin et al., 2017). In combination, these changes imply that the amplitude of the seasonal ozone cycle has decreased from its level in the 1990s, a time period characterized by decreasing precursor emission concentrations across northern midlatitudes.

80
Ozone precursor emission changes have been hypothesized as the cause of shifts in the seasonal ozone cycle in some studies, which have reached a consensus about how emissions affect the seasonal cycle. In the absence of large anthropogenic precursor emissions, the seasonal ozone maximum occurs in the spring (Logan et al., 1985; 2014 and references therein). Higher emissions correspond to a seasonal cycle of larger magnitude with seasonal maximum ozone occurring later in the year (in summer); likewise, lower emissions correspond to a seasonal cycle of smaller amplitude 85 with an earlier-occurring spring maximum (Parrish et al., 2013;Clifton et al., 2014;Cooper et al., 2014;Strode et al., 2015).
Most of these studies investigate areas affected only by baseline conditions or areas where local emissions have been controlled in recent years, and thus capture the decrease in amplitude and shift towards an earlier maximum. However, there is some analysis of areas with increasing emissions, and the seasonal cycle grew in amplitude with a progressively later maximum. For example, NOX emissions roughly tripled since 1990 in parts of China, and summertime ozone has increased 90 at many polluted sites (e.g., sites directly downwind of these increasing emissions) by up to 2 ppb/year (Li et al., 2017); thus on local to regional scales in China, increasing emissions correspond to a growing seasonal cycle with a shift towards summer. Most of the European and North American studies were largely based on observations and simulations from the late 20 th Century and early 21 st Century, when emissions were generally decreasing. One exemption is the study of Marenco et al. (1994) that noted the preindustrial 19 th Century seasonal maximum occurred in the spring at a remote European site, but 95 that maximum had shifted towards the summer by the 1980s. Taken together, these results suggest that the increase of anthropogenic precursor emissions during industrial development shifts the ozone seasonal cycle toward the summer, and reductions in those emissions allow the seasonal cycle to shift back toward the preindustrial condition.
Other studies identify correlations between precursor emissions and a changing seasonal cycle at sites separated geographically (instead of at the same site studied across a period of time). For example, sites in eastern Canada are subject 100 to less pollution than sites in the eastern US and subsequently show smaller summertime and larger wintertime ozone concentrations, evidence that the amplitude of the seasonal cycle is smaller in the absence of precursor emissions (Chan et 4 al., 2010). Across the same sites, a springtime ozone maximum is observed for more pristine Canadian sites, while the more polluted eastern US displays a summertime maximum (Chan et al., 2009). Similarly, Clifton et al. (2014) note that the polluted northeast region of the US displays a summertime seasonal ozone maximum, while the more pristine intermountain 105 west region displays a springtime seasonal maximum.
A quantitative understanding of the link between precursor emissions and the seasonal ozone cycle will benefit air quality policy development. Rieder et al. (2018) note that changes to the seasonal ozone cycle may influence the timing and number of days of ozone exceedance above the US National Ambient Air Quality Standard (NAAQS). As such, understanding the seasonal cycle -including how it changes in response to changing emissions -may usefully inform air quality control 110 managers across the world in setting future ozone standards in efforts to reduce the harmful impacts of surface ozone (Lin et al., 2017). A changing future climate will bring further uncertainty, including the possibility of an ozone-climate penalty in which rising temperatures can offset tropospheric ozone reductions from precursor emission control (Rasmussen et al., 2013;Clifton et al., 2014;Jaidan et al., 2018). Moreover, Schnell et al. (2016) report that a warming climate may cause changes to the seasonal cycle of tropospheric ozone independent from seasonal cycle effects of changing emissions. Jaidan et al. (2018) 115 and Rieder et al. (2018) describe the possibility of a similar "methane penalty" in which rising methane concentrations offset the reduction of other ozone precursors. Under some modeled scenarios, a future reversal of the present-day seasonal ozone cycle may be possible due to ozone precursor emission reductions in combination with the ozone-climate penalty (Clifton et al., 2014). We can reduce some of this uncertainty by understanding the interactions between emissions and atmospheric impacts (e.g., the seasonal ozone cycle) as fully as possible.

120
Despite the extensively documented record of shifts in the seasonal ozone cycle, no previous study has quantitatively analyzed measured data and model simulation results from across the northern midlatitude region, examined shifts in the amplitude and phase of the seasonal ozone cycle, quantitatively analyzed changing precursor emissions alongside seasonal cycle shifts, and proposed the mechanisms by which changing emissions affect the seasonal cycle; this paper aims to accomplish these tasks. We examine sites representative of baseline conditions in both Western Europe and North America.

125
Given the zonal similarity of ozone at northern midlatitudes, our analysis is expected to be representative of the baseline troposphere throughout northern midlatitudes. We investigate seasonal ozone cycle changes that began ~75 years ago, before reliable ozone measurements are available; thus, we rely on historical simulations from Coupled Model Intercomparison Project Phase 6 (CMIP6) Earth System Models (ESMs) as our primary basis for seasonal ozone cycle analysis. We compare these simulation results to available observations. A previous study of seasonal ozone cycle found that the previous 130 generation of Earth system models poorly simulated the seasonal cycle, including changes to it (Parrish et al., 2013).
However, Griffiths et al. (2020) find that CMIP6 ESMs capture the general shape of the observed seasonal ozone cycle averaged between 30° and 90° N, despite a positive bias of 3-4 ppb in overall ozone concentrations. Thus, CMIP6 ESMs may be more reliable for ozone seasonal cycle analysis than previous models.
In this work we investigate two quantities that define the seasonal cycle of tropospheric ozone: the amplitude (the 135 difference between the annual average and the annual maximum or minimum ozone concentrations) and the phase (the 5 timing of annual maximum ozone concentrations). Model simulations indicate that both of these quantities have changed over past decades; we compare their shifts with temporal changes in ozone precursor emissions that are prescribed in the models. In the following sections, we describe our analytical methods, present the analysis results, and discuss those results within a broader context. The overall goal of the paper is to provide a quantitative analysis of the shifting seasonal cycle of 140 tropospheric ozone at northern midlatitudes. Since the seasonal cycle reflects the sources and loss of ozone, quantifying it provides the opportunity for comparison of the simulated seasonal cycle between different models, and for comparison of model simulations with the limited record of observations. Comparing models and observations is an important way to gain insight into the performance of the models.

145
In this work, we seek to quantify the ozone seasonal cycle based on a small set of parameter values that reflect the amplitude and phase of that cycle. To accomplish this quantification, we analyze monthly mean ozone concentrations from ESM simulations as well as observations, to the extent they are available. Monthly means have sufficient temporal resolution to capture seasonal changes, while effectively averaging over most variability driven by diurnal and meteorological changes.
Our goal is to investigate the long-term changes in the seasonal cycle over the period included in the CMIP6 historical 150 simulations (Eyring et al., 2016), which were designed to extend form 1850 to the near present, which was set as 2014.
Ozone varies systematically on decadal scales, and also has temporal variability on interannual scales (i.e., on the scale of a few years) driven by changes in large scale transport patterns in the troposphere. For our purposes, this sub-decadal variability is "noise"; we minimize the obscuring effects of this variability by selecting analysis techniques that effectively average over this variability.  Table S2 gives references 6 for descriptions of these ESMs and their model output. Results of CMIP6 model simulations are archived at the Earth System Grid Federation (https://esgf-node.llnl.gov/projects/cmip6/) and are freely available to download. We obtained 165 monthly mean ozone concentrations for all six ESMs at the model levels that correspond to the selected comparison locations. Where available, a mean of multi-ensemble members was calculated for each model from the CMIP6 historical simulations over the period 1850-2014.

Fit Equations
We fit time series of monthly mean ozone concentrations with the following equation, which has separate functions for the 170 average long-term change (LTC) and the superimposed seasonal cycle (SC): where t is time in years. We find representing LTC(t) by a 5-term power series (i.e., a 5-term polynomial): LTC(t) = intercept + slope*t + curve*t 2 + d*t 3 + e*t 4 , captures the long-term changes and allows detrending of the monthly means for quantification of the seasonal cycle. The 175 choice of five polynomial terms in the power series is arbitrary: a range of polynomial terms (2 to 12) successfully detrends the monthly means without affecting seasonal cycle parameter values; Section S1 of the Supplement gives further details.
Quantification of the seasonal cycle is complicated by significant shifts in both the amplitude and phase of the seasonal cycle over the last ~75 years of the time series, as illustrated in Figure 1. Annual maximum ozone moved from primarily March and April (pink, purple, or blue) before 1900 to primarily June (blue-green and light green) by the 1980s, while the 180 vertical spread of the time series increased from ~12 ppb before 1900 to ~20 ppb by 1980. SC(t) must capture both the preindustrial seasonal cycle and the seasonal cycle shifts beginning around 1940. A 2-term Fourier series quantifies the preindustrial seasonal cycle (PISC) in the detrended monthly means: where A1 and φ1 are the amplitude and phase, respectively, of the fundamental (one sine cycle year -1 ), and A2 and φ2 are the 185 amplitude and phase, respectively, of the second harmonic (two sine cycles year -1 ) of the Fourier Series. As discussed in Section S2 of the Supplement, the fundamental is generally larger in magnitude than the second harmonic; together these two harmonics capture nearly all variance associated with the seasonal cycle, so higher order harmonics are not quantified. harmonic to consistently shift across locations and models. Note that in Equation 4, the A1 and φ1 parameters characterize the amplitude and phase of the preindustrial seasonal cycle, since the Gaussian maxima contribute negligibly before ~1900.

Inclusion of two
Thus, the sums (A1 + r) and (φ1 + rφ) derived here are most appropriate to compare to the A1 and φ1 values derived from a similar equation without the Gaussian terms (e.g., in Parrish et al., 2019), since that earlier work analyzed observations collected only in the last few decades.
In the following analysis ozone time series are fit to this equation, which consistently captures more than 95% of the variance in the time series of monthly mean ozone. In these fits, A linear function: quantifies the early long-term change in average precursor emissions. The EG term of Equation 6 is given by a Gaussian function that is consistent with an increase and then decrease in emissions primarily driven by anthropogenic activity: 215 EG(t) = rem*exp(-((t-mem)/sem) 2 ) This term is analogous to the Gaussian functions included in Equations 4 and 5; rem represents the maximum, mem represents the year of that maximum, and sem represents its width. Substitution of equations 7 and 8 into Equation 6 gives Equation 9: which captures more than 98% of the variance in the precursor emission time series.

220
Fitting of ozone and precursor emissions time series with these similar equations allows for quantitative comparison between the ozone seasonal cycle shift and the growth and decrease in emissions. Specifically, comparison of the ozone seasonal cycle Gaussian parameters in Equations 4 and 5 with the precursor emission parameters in Equations 8 and 9 is the basis for examining the correlation between changes in the seasonal cycle and changing ozone precursor emissions.
The multivariate regression fits of Equations 5 and 9 to time series of monthly mean ozone concentrations and annual 225 emissions, respectively, quantify confidence limits for all derived parameter values. In this work 95% confidence limits are tabulated and discussed throughout. However, these confidence limits only reflect the variability of the time series about the functional form fit to that time series, and this approach assumes that each member of the time series is an independent 8 variable with no autocorrelation within the time series; hence it must be recognized that these confidence limits are lower estimates of the actual uncertainties of the derived parameter values.

Selected CMIP6 Simulation Locations
The CMIP6 ESMs provide monthly mean ozone concentrations on global grids. To focus our investigation on northern midlatitudes, model-simulated monthly mean ozone time series are taken from model cells at three locations in western Europe and three in western North America. European locations are surface sites at Hohenpeissenberg, Germany and 235 Jungfraujoch, Switzerland, and in the FT above Jungfraujoch at altitudes between 5 and 6 km. North American locations are located in California -one surface site in the US at Lassen Volcanic N.P., and in the FT above Trinidad Head at two different altitudes: between 0.9 and 1.2 km in different models (we refer to this site as Trinidad Head at 1 km for simplicity), and between 5 and 6 km. Table 1a  America. The sites on each continent are all within ~500km; given the pronounced zonal similarity of ozone concentrations at midlatitudes (Parrish et al., 2020), the geographic separation between these sites has negligible impact on ozone concentrations, so the two sets of three sites are representative of their respective continents at different altitudes or elevations. In summary, the locations are selected to provide an altitude-dependent contrast between the western regions of the two continents.
9 a Although sonde measurements were not conducted above Jungfraujoch, we analyzed model simulation results above Jungfraujoch for comparison with the European FT measured data set (see Table 1b) b Offshore location selected for comparison with onshore and is only considered in the Supplement. Details included in Section S4.

Ozone Observations
Although model simulations are our main basis for analysis, observational data are also considered. Shifts in the amplitude and phase of the seasonal ozone cycle are generally apparent in observational records that span long enough time periods. impetus behind averaging measurements from three surface sites and three sonde data sets is to reduce the impact of ozone variability in any one data set, and to thereby obtain a more precise quantification of ozone over Western Europe. These same data sets have been considered in previous studies of western European baseline ozone concentrations 2020). The Hohenpeissenberg data discussed by Parrish et al. (2014) are here extended through 2016, and the European alpine and sonde data sets are the same as those analyzed by Parrish et al. (2020). The North American data sets 275 include one spanning 30 years, 1998-2017, at Lassen Volcanic N.P., and one spanning ~ 21 years, late 1997-early 2018, from sondes launched from Trinidad Head. We consider average sonde measurements between 5 and 6 km, and between 0.5 and 1.0 km. These North American data sets are also the same as those analyzed by Parrish et al. (2020). Table 1b summarizes location details of observational data considered in the analysis. 280 10

Precursor Emissions
Annual mean ozone precursor emissions were derived from ESM emission inventories integrated over the northern midlatitude region between 30° and 60° N for the 1850-2014 simulation period. The primary analysis examines emissions of NOX and VOCs from anthropogenic (Hoesly et al., 2018) and biomass burning sources (van Marle et al., 2017) that were provided as a common emission inventory to be used by all models (including the six in this study) in CMIP6 simulations.

290
As discussed in further detail in Section S5 of the Supplement, the anthropogenic emissions dominate this inventory.
Although there are small seasonal cycles in these emissions, these seasonal cycles are either approximately constant over the entire time interval, or their relative magnitudes are small compared to that of the seasonal cycle of ozone; further discussion is included in Section S5.
Even though the six ESMs used the same prescribed anthropogenic and biomass burning emissions, Figure 1 of Griffiths 295 et al. (2021) shows that subtle differences remain in NOx emissions and even greater differences in CO and biogenic VOC emissions between models. Differences in the VOC emissions arise because the speciated VOC emissions that were provided had to be mapped onto the chemical mechanisms in the individual models, and this mapping may not fully account for the total VOC emissions prescribed. The emissions that are the focus of our analysis have been taken from the prescribed emission inventory; we have not further diagnosed the exact northern midlatitude emissions actually used in each individual 300 model. Some of the models were able to provide quantifications of emissions from biogenic and other natural sources for evaluation. These emissions varied between models based on model-specific chemistry and parameterizations, and included biogenic VOC emissions (specifically, isoprene) and dimethyl sulfide (DMS) from oceans, and NOX emissions from soil and lightning. Methane is considered independently due to its very long lifetime compared to other VOCs; all of the ESMs use 305 prescribed global annual mean values of CH4 concentrations as input at the surface throughout the whole historical period (Meinshausen et al., 2017). Further details of model-specific natural emissions and CH4 are given in Section S6 of the Supplement.

310
Ozone in the troposphere varies on a wide spectrum of temporal scales, which makes it difficult to quantify a particular contribution to that variability. To isolate the seasonal cycle, we examine time series of monthly mean ozone concentrations.
Monthly means integrate over the short-term variability driven by diurnal cycles and short-term meteorological changes, which effectively removes their influence. The time series considered here span a maximum of 165 years, which allows significant influence from "longer-term" (i.e., on the scale of decades to centuries) variations driven by ozone precursor 315 emission changes and climate variations. We isolate the seasonal cycle from these longer-term changes by detrending the monthly mean concentrations, which we accomplish by subtracting LTC(t) in Equation 1 from the time series of monthly means. A regression fit of the five-term polynomial given in Equation 2 to the time series of monthly means gives values for the five polynomial coefficients; Section S1 of the Supplement discusses more details of the determination of LTC(t). Figure   2 shows the example time series of Figure 1 with the fit to Equation 2 indicated by the black curve, which quantifies the 320 longer-term temporal change that underlies the time series of monthly means. examine are described sufficiently by the sum of the first two harmonics: the fundamental (one sine cycle yr -1 ) and the second harmonic (two sine cycles yr -1 ) as indicated in Equation 3. The fundamental is generally larger in magnitude than the second harmonic, except for the two lower-elevation North American sites, for which the two harmonics were approximately equal in magnitude during the preindustrial period. In combination, the fundamental and second harmonic capture almost all the variance associated with the seasonal cycle. A quantitative Fourier analysis that provides the basis for this harmonic 335 analysis and the inclusion of only the first two harmonics is detailed in Section S2 of the Supplement. The detrended seasonal cycles, including their evolution over the course of the 1850-2014 period, are analyzed for the six ESM simulations at the six selected northern midlatitude locations; this is the primary basis of our analysis, which is discussed in the next three subsections.

340
To understand the magnitude and timing of changes to the seasonal ozone cycle that began near the middle of the 20th Century, it is important to quantify the seasonal cycle before those changes began, i.e., the preindustrial seasonal cycle. Only  13 separated from ground-based sources of emissions than are surface sites. We assume these locations are representative of the FT baseline seasonal ozone cycle with little influence from local or regional emissions; thus, they are appropriate for our 360 initial analysis. At both FT locations, the preindustrial seasonal cycle is similar in character; it is determined largely by the fundamental, which generally reaches its seasonal maximum in May or June. Figures S4 and S5 give plots similar to Figure   3 for the four lower-elevation locations, with discussion in Section S3 of the Supplement.
Between different models, there are important qualitative similarities in the simulated preindustrial seasonal cycle at most locations. First, the fundamental is larger in magnitude than higher-order frequencies, except for Trinidad Head at 1 km and 365 Lassen NP. Second, the maximum of the fundamental occurs in the late spring or early summer, which drives an overall seasonal maximum that also occurs in the late spring or early summer. Marenco et al. (1994)

Seasonal cycle shifts across northern midlatitudes
Across all models and all locations, shifts in both the amplitude and phase of the seasonal cycle are ubiquitous. Importantly, the presence of a seasonal cycle shift in the FT indicates it is a hemisphere-wide phenomenon, rather than limited to a localized environment. Figure 4 compares preindustrial and modern-day seasonal cycle simulations from one example model 390 with the observed modern-day seasonal cycle in the FT over Europe. This is the same example time series shown in Figures   1 through 3. The modern-day seasonal cycle is larger in amplitude with a later maximum compared to the preindustrial seasonal cycle. These changes are primarily driven by the changing fundamental, rather than the second harmonic, which makes only a small contribution in the FT and is not statistically different between the preindustrial and modern-day simulations. The modern-day simulated seasonal cycle approximates, but does not exactly match observations from the past 395 two decades; the simulated seasonal cycle is smaller in amplitude with an earlier maximum than the measured seasonal cycle.
The noted that the linear fit to the phase shift is closely related to an earlier analysis approach (Parrish et al., 2013) that also quantified the phase shifts from observations at some of these same sites.

Connection between ozone precursor emissions and the seasonal cycle
All six CMIP6 ESM simulations incorporated the same ozone precursor emission inventory for anthropogenic and biomass burning sources. Figure 7a illustrates the temporal evolution of these nonmethane VOC and NOX emissions integrated annually and over the entire northern midlatitude region (30° to 60° N). The curves in Figure 7b are fits of 440 Equation 9 to those emissions; these fits (with the underlying linear increases) capture more than 98% of the variance in the time series of annual emissions. Equation 9 is designed to provide Gaussian function fits to the emissions, so that the derived parameters can be directly compared to the Gaussian parameters that quantify the shift of the ozone seasonal cycle. Figures S10 and S11 compare the Gaussian parameters from the emission fits with those derived from fits to the model simulated ozone at individual sites.

475
Interpretation of the confidence limits quoted for the derived parameters is difficult. The multivariate regressions utilized to fit the model simulations, observations and emissions return parameter values with 95% confidence limits, which are plotted in Figures 3, 4 and 8; many are not visible because they are smaller than the plotted symbols. These confidence limits are underestimated (see Section 2.2) due to autocorrelation in the time series of monthly mean ozone concentrations. An 480 independent estimate of the confidence limit of each overall average parameter value can be obtained from the variance of the individual parameter values included in the average. If one assumes that each seasonal shift parameter must be identical a b c d e f at all six locations, and that each model simulation at each site provides an independent determination of that parameter value, then the confidence limit of the average can be estimated from the square root of the variance divided by the number of independent model determinations (36 if the fits to each of the six model simulations returns a parameter value at each of 485 six locations). Such upper limits are included in parentheses in the bottom line of Table 3 and are included as the blue error bars on the overall averages in Figure 8; they are larger by factors of 2.5 to 17 than those derived from the weighted averages of the parameters from the regression fits. In quantitative comparisons of the parameters from observations and emissions with those simulated, this issue with the confidence limits must be considered. Table 3. Gaussian parameters that define changing emissions and seasonal ozone cycle shifts over northern midlatitudes. First two 490 rows give fit parameters for total anthropogenic and biomass burning ozone precursor emissions integrated across the entire northern midlatitude region (30° to 60° N); second two rows give parameters for fits to observed and model simulated seasonal ozone cycles. Positive r values for the phase shift indicate a seasonal cycle shifting towards an earlier annual maximum. The seasonal cycle Gaussian parameters are averaged over the six locations considered in the analysis.

Discussion and Conclusions
We analyze the seasonal cycle of tropospheric ozone over the historical period, as simulated by six CMIP6 Earth System Models and deduced from available observations at six northern midlatitude locations in western Europe and western North 500 America. Over the time period of the model simulations (1850-2014), the seasonal cycles shifted significantly in both phase and amplitude at all locations, including within the free troposphere. The seasonal cycles simulated by the models remained generally constant from 1850 until well into the 20 th century; this preindustrial seasonal cycle is shown in Figure 3 for two FT locations and in Figures S4 and S5 for four lower-elevation locations. In the period from approximately 1920-1940 the seasonal cycle amplitude began to increase, and the seasonal maximum began to shift to later in the year. These changes 505 reached their maximum extent late in the 20 th Century, after which they began to reverse -the seasonal cycle decreased in amplitude and the annual maximum shifted back to earlier in the year. Gaussian functional fits quantify these shifts.
Observations are available for at most only the last 44 years of the model simulations; within their large uncertainties (see error bars in Figure 8) the available measurements indicate seasonal cycle shifts similar to those simulated. Figure 4 illustrates these shifts as simulated by one model at one location; it shows that the fundamental harmonic is the primary 510 contributor to both the seasonal cycle and its shifts. Figure 4 also compares the simulated modern-day seasonal cycle with 20 that derived from observations. Figures 5 and 6 show comparisons of the shifting amplitude and phase of the fundamental harmonic among all models and with available observations at the six locations considered. The Introduction discussed extensive literature reports of modelled and observed changes in the seasonal ozone cycle throughout northern midlatitudes over the most recent three to four decades; the seasonal cycle shifts examined here are generally consistent with those 515 reports.
Throughout northern midlatitudes, on average (blue symbols in Figure 8) the simulated shifts in both the amplitude and phase of the fundamental of the seasonal cycle maximize at similar times (~1985; Figures 8a, b) with the amplitude shift having a somewhat smaller width (~30 years; Figure 8c) than the phase shift (~40 years; Figure 8d). At the maxima, the fundamental amplitude (Figure 8e) had increased by ~5.5 ppb (i.e., a ~11 ppb increase in the difference between the seasonal 520 minimum and maximum), and the seasonal maximum (Figure 8f) had shifted to ~3 weeks later in the year. For comparison, the average simulated preindustrial seasonal cycle in the free troposphere had an amplitude of ~7 ppb and a seasonal peak near June 1 (Figures 5 and 6). The sparse measurement record from the European alpine sites and Hohenpeissenberg (red and green points in Figure 8 and entries in Table 3) agrees well for the timing of the maximum shifts, but suggests somewhat smaller seasonal cycle changes in the widths and magnitudes of those shifts; however, the large uncertainty in the 525 observational determinations should be noted.
The model simulations exhibit large variability, both among models and locations (compare points on right side of graphs in Figures 8, S10 and S11); however, it is difficult to judge if this variability is statistically significant. Here we identify some aspects of this variability that appear to be robust. First, the relative spread among the model averages in the phase shift is greater than that in the amplitude for all three parameters (year of maximum, and the width and magnitude of the 530 Gaussians quantifying the shifts). Second, both the amplitude and phase shifts appear larger and more varied at lower elevations compared to the FT (compare lower graphs in Figures 5 and 6 with the FT results in the upper graphs, and the low and high elevation averages in Figures 8e,f and S12); since the anthropogenic emissions are located at the surface, this behavior may reflect the greater influence from local and regional emissions at the surface sites compared to the more isolated locations in the FT. Third, Hohenpeissenberg (located at a relatively low elevation in central Western Europe) 535 generally shows the largest amplitude shifts in the model simulations as well as in the measurements, although the measurement results are highly uncertain. At Hohenpeissenberg (Figure 5c) all six models simulated the timing of the maximum amplitude shift (i.e., the m parameter) within the uncertainty of that derived from the measurements (1987 ± 6 years). This temporal agreement occurs despite disagreement (by a factor of ~2) in the maximum fundamental amplitude (peaks of Gaussian curves in the Figure 5c) and disagreement (up to a factor of ~3) by 2 of the 6 models in the amplitude of 540 the preindustrial fundamental (horizontal portion of the curves at the left of figure 5c). There is poorer agreement regarding the phase shift at Hohenpeissenberg, with the simulated maxima occurring between 1984 and 2000 in the six model simulations; the timing of the maximum phase shift derived from the measurements is not precisely defined, but its confidence limits include (nearly) all of the model results. Fourth, the greatest variability of the simulated phase shifts is seen at Trinidad Head at 1 km (Figure 6f), which is the lowest elevation North American location considered here; there the 21 maximum of the fundamental is found to occur in nearly every month of the year over the simulation period in at least one of the model simulations, although the seasonal cycle amplitude is relatively small at this location, which makes determination of that maximum difficult. A possible explanation for these divergent model results is that this location is on the edge of two transitions -the MBL to FT and the marine to continental environment -which may be a particularly difficult situation for the models, which have coarse horizontal resolution, to simulate.

550
We also quantify the temporal changes in total northern midlatitude ozone precursor emissions from anthropogenic and biomass burning sources (Figure 7) that are incorporated into the emission inventories assumed by all of the ESMs. Between 1850 and 1940 emissions increased only slowly, with more rapid increases beginning in the mid-20 th Century as the result of rapid industrialization in Europe and North America. By the late 20 th Century, emissions began to decrease as the result of air quality control efforts in more developed countries. Notably, the changing emissions are driven by anthropogenic 555 activity; Section S5 of the Supplement compares the temporal changes of anthropogenic and biomass burning emissions, and shows that it is only the anthropogenic emissions that rise and fall over time, while the biomass burning emissions remain approximately constant.
On average, the parameters that quantify the shifts in the seasonal cycle correlate strongly with those that quantify the emissions. Figures 8 a-d show that the overall model simulation averages of the four parameters that quantify the timing of 560 the shift of the amplitude and phase of the fundamental harmonic closely correspond to the parameters that quantify the temporal evolution of the emissions. There is significant variability in the results from the different models (open circles on the right in the graphs in Figure 8), but that variability is reduced in four regional averages (open circles on the left). There is no consistent, strong difference between the European and North American continents. Both the amplitude and phase shifts apparently maximized earlier in North America than Europe, but there is a great deal of variability among the individual 565 determinations (Figures S10 and S11) so the statistical significance of this apparent difference is uncertain. There also may be significant differences in the shapes of the Gaussian describing the phase shift between the lower-elevation surface sites (i.e., earlier years of maximum shift and greater widths) compared to the higher-elevation sites representative of the FT (Figures 8b,d); and the phase shift at high elevations appears to have maximized later with a smaller width. The maxima of the amplitude and phase shifts (Figures 8e,f and S11) are apparently larger at the low elevation sites, which may reflect more 570 direct impact by anthropogenic emissions.
Based on the temporal correlation between the emission changes and the seasonal cycle shifts shown in Figure 8, we hypothesize that the changing ozone precursor emissions is the cause of the shifts in the seasonal ozone cycle throughout northern midlatitudes. During industrial development, ozone production driven by rising anthropogenic precursor emissions progressively becomes the predominant source of ozone, which shifts the ozone seasonal maximum into the summer, when 575 photochemical ozone production is more important (compared to, e.g., ozone input from stratospheric intrusions, which peaks in the spring). Ozone production driven by anthropogenic activity also increases the amplitude of the seasonal cycle by boosting summertime concentrations while wintertime concentrations are less affected. As emissions decrease, those changes reverse, with the seasonal cycle returning toward the pre-industrial cycle. Although ozone precursor emissions from 22 all sources influence ozone production and the ozone seasonal cycle, it is anthropogenic activity that drives the seasonal 580 cycle changes; more discussion of natural and anthropogenic emissions is given Sections S5 and S6 and Figure S7-S9 of the Supplement. The temporal correlation between the changes in emissions and the ozone seasonal cycle does not necessarily prove our hypothesis, but a comparison of ozone seasonal cycles between southern and northern mid-latitudes does support a causal relationship. As discussed and illustrated in Section S7 of the Supplement, we find no evidence of a significant shift in either the phase or magnitude of the ozone seasonal cycle at southern midlatitudes. The presence of a shift in the ozone 585 seasonal cycle throughout northern midlatitudes and its absence at southern midlatitudes is as expected from our hypothesis, due to the much smaller anthropogenic ozone precursor emissions in the southern hemisphere. For reference, Figure 8 of Crippa et al. (2020) illustrates the dramatic difference in emissions from fossil fuel combustion between hemispheres.
An interesting aspect of the correlation between precursor emissions and the ozone seasonal cycle shifts is the temporal offset in the evolution of the emissions. The Gaussian functions fit to the non-methane VOC emissions in Figure 7 (see 590 Table S3a for parameter values) peaked in ~1983 with a full width at half maximum (FWHM, which is a factor of 1.67 larger than the Gaussian s parameter) of ~47 years, while the fit to the NOX emissions peaked in ~1995 with a FWHM of ~63 years. The shifts in the amplitude and the phase of the average simulated ozone seasonal cycle both reached peaks in ~1985, closely corresponding to the VOC emission peak. The FWHM of the ozone seasonal cycle amplitude shift (~48 years) also closely matches the FWHM of the VOC emissions. In contrast, the FWHM of the ozone seasonal cycle phase 595 shift (~65 years) corresponds more closely to the FWHM of the NOX emissions. A simple hypothesis can provide a qualitative explanation for this correspondence. The VOC emissions provide fuel for photochemical production of ozone; thus these emissions exert primary control of the seasonal cycle amplitude driven by summertime production. The NOX emissions provide the catalyst that determines whether photochemistry produces or destroys ozone -once the NOX emissions are large enough that photochemical production dominates the seasonal cycle and moves the seasonal maximum 600 into the summer, the phase shift ends, since the maximum cannot continue shifting into the autumn, and the seasonal maximum will not shift back until NOX emissions decrease to levels low enough that photochemical production no longer dominates the ozone budget. In summary, we are suggesting that the NOx emissions largely control the timing of seasonal maximum in ozone, while the VOC emissions control the seasonal cycle amplitude. If this hypothesis is correct, consideration of the role of biogenic VOCs could help to explain some of the diversity in the seasonal cycles and shifts seen 605 among the model simulations; as can be seen in Figure 1 of Griffiths et al. (2021) the temporal variation of the biogenic VOCs emissions are significantly different across the models.
Assuming that the above hypotheses are correct, the ozone seasonal cycle shift derived from observations must reflect the time evolution of emissions, and thereby provide tests of the emission estimates upon which the model simulations are based. The measurement records (maximum of 44 years) are so short that the precision of the parameters of the seasonal 610 cycle shift that can be derived from the measurements (see Table 3) is limited, as indicated by the relatively large confidence limits for those parameters included in Figure 8. However, two points can be noted. First, the average year of the maximum shift in the amplitude of the observed ozone seasonal cycle (1990 ± 3 years) is later than the maximum of the VOC emissions (1983 ± 1 year); since we expect these two maxima to be the same, this disagreement may indicate that anthropogenic VOC emissions actually peaked a few years later than indicated in the emission inventory. The uncertainty in 615 the year of the maximum phase shift determined from observations (1985 ± 8 years) prevents a precise comparison between the emission maxima and the phase shift maxima. Second, the widths of the amplitude and phase shifts of the observed seasonal cycle (15 ± 9 and 17 ± 17 years, respectively) appear to be smaller than the widths of the NOX and VOC emissions (38 ± 2 and 28 ± 2 years, respectively).
The seasonal cycle of ozone reflects the annual variability of the sources and sinks of ozone; thus its accurate simulation is 620 expected to present a stringent test for models. Given the paucity of the observational ozone record, both spatially but more importantly temporally, improved confidence in our understanding of changes in the seasonal ozone cycle must primarily come from improved agreement between different model simulations. Our analysis has focused on changes in anthropogenic and biomass burning emissions, which were prescribed from the same source to the extent possible for all models; however, there were differences in implementing the prescribed emissions into the models, mainly from VOCs due to individual 625 requirements of the chemistry scheme within each model. In addition, the representation of natural emissions (e.g. biogenic VOCs emitted from vegetation) differed between individual models, giving variation in the natural to anthropogenic emission ratios between models. Thus, remaining differences in emissions between models may cause some of the intermodel differences. More generally, Griffiths et al. (2020) suggest that differences in the simulation of ozone from CMIP6 models could be due to inter-model variations in the treatment of chemical and physical processes including dynamic 630 transport, stratosphere-troposphere exchange, photolysis, deposition, convection and boundary-layer schemes. There is a need to go beyond direct model-observation comparison studies; for example, multi-model perturbed parameter ensembles can be used to intercompare the sensitivity of models to different input parameters and/or parameterizations (Wild et al., 2020). Notably, in this work we document relatively large seasonal cycle shifts that are common to the entire northern midlatitude baseline troposphere; given the magnitude of these shifts, which we attribute to changing precursor emissions, it 635 may be difficult to independently determine the effects of other factors, e.g. changing climate (Fowler et al., 2008;Clifton et al., 2014), on the northern midlatitude ozone seasonal cycle.
Data Availability: All of the data utilized in this paper are available from public archives; most are fully discussed and referenced in the Acknowledgments section of Parrish et al. (2020). One additional data set is included in this paper, the 640 surface observations at Hohenpeissenberg, which are available from the Tropospheric Ozone Assessment Report data base (Schultz, M.G. et al. (2017).