Sensitivity to meteorology of regional contributions to air pollution in eastern Canada: part 1: ozone and NO_X

Abstract

Air pollution is well known to have harmful effects on human health. The mechanisms of production, elimination and transport of air pollutants are sensitive to atmospheric conditions. Therefore, future concentrations of air pollutants, and the transboundary transport of these air pollutants, may vary with climate change. Based on GEOS-Chem chemical transport model output, we created emulators to capture the effects of meteorology on the concentrations of O₃ and NO_X in eastern Canada. We then used the emulators with regional climate model output to project how climate change will affect the concentrations of these pollutants. We project that spring O₃ concentrations in eastern Canada will increase by ∼5 ppb in the RCP8.5 scenario, primarily due to increases in temperature. We project that NO_X concentrations will increase by less than 0.5 ppb, except for the greater toronto area where we project increases of more than 1 ppb. We also created emulators based on additional model simulations in which anthropogenic emissions were set to zero in one of three regions: the province of Quebec, the rest of Canada, and the US. We used these additional emulators to project how the contributions from each of these three regions to O₃ and NO_X concentrations may be altered due to climate change. We project that climate change will increase the US contribution to O₃ concentrations in eastern Canada more than the Canadian and Quebec contributions. Higher temperatures increase the efficiency of O₃ formation from anthropogenic precursors, thus enhancing pre-existing disparities due to the different quantities of O₃ precursor emissions from each region. This result suggests that climate change has the potential to exacerbate the export of air pollution across political boundaries.

---

1 Introduction

Air pollution is well known to have harmful effects on human health.^1–3 However, the lifetimes of many common air pollutants span from days to weeks, allowing these air pollutants to cross national and sub-national borders, negatively affecting human health in other countries and even other continents than those where they were emitted.^4–7 The mechanisms of production, elimination and transport of air pollutants are sensitive to atmospheric conditions.^8–10 Therefore, future concentrations of air pollutants could vary with climate change, and the transboundary transport of these air pollutants may also vary with climate change.

Our analysis focuses on three pollutants: fine particles with a diameter of less than 2.5 µm (PM_2.5) NO_X (nitrogen oxide (NO) + nitrogen dioxide (NO₂)) and O₃ (ozone). We will examine PM_2.5 in a follow-up paper and we discuss NO_X and O₃ in this paper. These three pollutants are known to be associated with adverse effects on human health.^1–3 In addition, their concentrations are currently used to calculate air quality indices in Canada and in other countries. Thus, it is important to better understand the sources of these pollutants as well as their future changes in order to develop policies that enable economic growth while reducing the impacts of poor air quality on affected populations.

Climate change is expected to alter NO_X concentrations through several mechanisms. Natural NO_X emissions from lightning and soil microbial activity are expected to be affected by climate change.^9,11 A decline in lighting over eastern Canada has been observed over recent decades.¹² However, the intergovernmental panel on climate change (IPCC) AR6 report assigns a low level of confidence to the sign of change in lightning NO_X emissions due to fundamental uncertainties in parametrizations of lightning NO_X, and quantitative estimates of climate-induced changes in soil NO_X are not yet represented in Earth system models.⁹ Higher temperatures increase peroxyacetylnitrate (PAN) decomposition rates, leading to local increases in NO₂ in polluted emission regions,^11,13 but higher water vapour mixing ratios also convert NO₂ to nitric acid (HNO₃) which can be quickly removed by wet deposition, thereby decreasing NO_X concentrations.^14,15

As O₃ formation rates are strongly influenced by NO_X concentrations, all of the previous relationships between NO_X and meteorology will have subsequent effects on O₃ concentrations. These effects are complicated in that both the magnitude and sign of this influence varies based on local concentrations of NO_X and volatile organic compounds (VOCs). Previous studies have also noted that there are complex interactions between O₃ concentrations and water vapour concentrations, as hydroxyl radical (OH) concentrations increase with water vapour concentrations, subsequently influencing concentrations of both NO_X and VOCs.^8–10 A warmer climate is predicted to decrease surface O₃ concentrations in regions far from pollution sources, because increased water vapour concentrations are expected to reduce the lifetime of O₃.¹⁴ Biogenic emissions of VOCs are generally expected to increase for small increases in temperature, but this relationship is both nonlinear and complex, depending on plant species and other ambient meteorological conditions.¹⁶ Near anthropogenic pollution sources, O₃ concentrations are expected to increase by several parts per billion (ppb).^8–10,14,17 The sensitivity of O₃ to projected climate change was estimated by an ensemble of four Earth system models as part of CMIP6, as the difference between simulations with sea-surface temperatures evolving due to climate change and a second set of simulations with sea-surface temperatures fixed at present-day values following the AerChemMIP experiment configuration.⁹ The four models were evenly split (two positive, two negative) on the sign of the change in O₃ in the most populated regions of eastern Canada for global average temperature increases of 1 or 1.5 °C. The same ensemble predicted average decreases of about 0.5 to 1 ppb for global average temperature increases of 2 or 2.5 °C.

Correlations between future changes in O₃ and surface air temperature for six models participating in CMIP6 were further examined by Turnock et al.¹⁸. Both emissions and climate varied during the experiment (scenario ssp370). The models disagreed on the sign of the change in O₃ in North America and on the sign of the correlation between O₃ and surface temperature in this region, indicating an unclear signal of the impact of climate change on O₃ concentrations in this region.

Both Kelly et al.¹⁹ and Turnock et al.¹⁷ have examined changes in O₃ due to climate change in the absence of changes in emissions. The former study projected small changes (generally less than 1 ppb) in O₃ concentrations over Canada by 2050, excepting the Great Lakes region where O₃ concentrations were projected to increase by several ppb. The latter study found that annual-mean surface concentrations of O₃ decreased by 13% in 2095 in regions far from major sources of anthropogenic pollution, due to increased quantities of water vapour which increase O₃ destruction. Conversely, they also found that climate change would locally increase O₃ production rates and increase annual-mean surface O₃ concentrations (up to 5 ppb by 2095) in the polluted continental regions of the northeastern United States (US). They note that this is a stronger response than occurs in other CMIP6 models and could be due to the larger temperature variation simulated by the model they used (UKESM1). Eastern Canada lies downwind of the northeastern US, and no statistically significant O₃ concentration changes were found for most of eastern Canada.

The effect of climate change on atmospheric transport is mixed in previous studies. Doherty et al.²⁰ and Murazaki and Hess²¹ have specifically studied changes in transport, but the effects on pollutant concentrations at the Earth's surface in eastern Canada are unclear. Doherty et al.¹⁴ found an increase in transport from the eastern US to eastern Canada, but the effects on concentrations of O₃ were small relative to changes due to differences in other processes, including formation and destruction of O₃. In contrast, Barnes and Fiore²² suggests that changes in the latitude of the jet stream over northeastern North America will have effects on the variability in O₃ concentrations at the Earth's surface and the relationship between O₃ concentrations and temperature.

Few previous studies have examined the effects of climate change on the fate of O₃ formed from anthropogenic precursors from a single source region in isolation, and even fewer studies have considered source regions smaller than the continental scale. This renders it more difficult to distinguish whether O₃ increases in a particular region are primarily due to the increased formation and survival of O₃ from local precursors, or an increase in transport, formation and survival of O₃ formed from precursors that were emitted upwind. To investigate how transboundary transport in eastern Canada may change in the future, we used a chemical transport model (GEOS-Chem) to quantify the concentrations of pollutants originating from the cross-border transport of emissions. GEOS-Chem has previously been evaluated specifically against observations in the Quebec region,^23,24 and has also previously been shown to adequately reproduce the observed response of O₃ to temperature in the US and Europe.¹¹ From the GEOS-Chem results, we developed an emulator that calculates the concentrations of O₃ and NO_X using several meteorological parameters. The emulator described here is essentially a simplified model linking meteorological variables and pollutant concentrations, and it was developed following similar methods to those described in Shen et al.²⁵. In other words, its aim is to capture the effects of changes in meteorology on the concentrations of pollutants in eastern Canada. We created three more emulators based on further GEOS-Chem simulations without anthropogenic emissions from either Quebec, the rest of Canada (excluding Quebec), or the United States to assess the sensitivity of O₃ and NO_X to changes in meteorology in the absence of anthropogenic emissions from each of these three regions. These four emulators are then driven by data for future climate change scenarios to project changes in pollutant concentrations due to climate change. More specifically, we use results from a regional climate model, which can better resolve small-scale features than the global models used in most of the studies referenced above. This allows us to examine how the regional anthropogenic contributions to O₃ and NO_X concentrations may change due to climate change. By examining which meteorological variables are responsible for the change, we can also investigate whether these changes are due to changes in transport or due to changes in formation efficiency and survival.

The sections of this study are organized as follows. Section 2 explains the methodologies used in detail, including the GEOS-Chem chemical transport model (Section 2.1), the creation of the emulator (Section 2.2), and the regional climate model data used (Section 2.3). Then, Section 3 presents our results: GEOS-Chem results in Section 3.1, an evaluation of the emulator in Section 3.2, and projections of the effect of climate change on O₃ and NO_X in Section 3.3. Finally, the conclusions of this report are summarized in Section 4.

2 Methods

2.1 Description of the GEOS-Chem simulations

We performed simulations using the GEOS-Chem chemical transport model (version 13.3.4, doi:10.5281/zenodo.5764874).^26,27 GEOS-Chem is driven by assimilated meteorological data from the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) which is produced by NASA's Global Modeling and Assimilation Office (GMAO). A global simulation at a resolution of 2° latitude and 2.5° longitude was performed to generate the boundary conditions which are then used in a second high-resolution simulation. The high resolution simulation is carried out over a domain which extends from 35° N to 65° N and from 90° W to 50° W in order to include the province of Quebec in its entirety, as well as the strong sources in the Great Lakes region and the northeastern US. This second simulation has a resolution of 0.5° latitude and 0.625° longitude (Fig. 1). We chose our simulation periods to sample a wide range of meteorological conditions as follows: The earliest year of MERRA-2 data available is 1980, and therefore the earliest full year we could simulate while allowing for spin-up is 1981. To limit the computational resources required, we simulated only the first five available years of each decade: 1981 to 1985, 1990 to 1994, 2000 to 2004, and 2010 to 2014 inclusive. We allow one month of spin-up time for each five-year simulation to ensure that the model is in a stable state that is not disrupted by the conditions at the start of the simulation. The atmosphere is resolved using 47 vertical layers extending from the surface to an altitude corresponding to an atmospheric pressure of 0.01 hPa. The vertical resolution of the model is about 100 m near the surface, but becomes coarser at higher altitudes.

Fig. 1 Geographical domain of the study region. The pale shaded region indicates the domain of the high-resolution GEOS-Chem simulations and the green shaded region indicates the domain of the emulator. Canadian provinces and territories are labelled by postal abbreviation. NL: Newfoundland and Labrador, PE: Prince Edward Island, NS: Nova Scotia, NB: New Brunswick, QC: Quebec, ON: Ontario, MB: Manitoba, NU: Nunavut.

All anthropogenic emissions are set to those of the year 2010, in order to study only the effect of meteorology on pollutant concentrations without interannual changes in emissions. Anthropogenic emissions for Canada are provided by the Air Pollutant Emissions Inventory (APEI), and anthropogenic emissions for the US are provided by the National Emissions Inventory for 2016,²⁸ multiplied by a species-specific scaling factor to account for changes in emissions between 2010 and 2016. Anthropogenic emissions outside of Canada and the US, as well as shipping emissions globally, are provided by the Community Emissions Data System version 2 (CEDSv2).²⁹ Fossil fuel and biofuel emissions of C₂H₆ are overwritten by the updated inventory of Tzompa-Sosa et al.³⁰. Outside of Canada and the US, emissions of C₃H₈ are overwritten by the inventory of Xiao et al.³¹. Aircraft emissions are provided by the Aviation Emissions Inventory Code (AEIC).^32,33

We allow natural emissions to vary from year to year in the simulations. Natural emissions of several chemical species are calculated online (rather than prescribed) in the GEOS-Chem model and can therefore be influenced by meteorological variability (see Keller et al.³⁴ for more details). We therefore take into account the impacts of meteorology on the concentrations of O₃ and other pollutants through changes in natural emissions. These emissions include NO_X emissions from lightning and ground processes, sea salt emissions, desert dust emissions, and biogenic VOC emissions. Biogenic VOC emissions in particular are important for O₃ formation, and vary in our simulations depending on leaf area index, solar radiation, temperature, and soil moisture.¹⁶ We also include biomass burning emissions, derived from historical global biomass burning emissions for the Coupled Model Intercomparison Project (CMIP)³⁵ to cover the full temporal period of our simulations. In northeastern North America, biomass burning emissions are primarily due to unmanaged wildland fire, as opposed to managed agricultural burning.

To examine the sensitivity of air pollutant concentrations to transboundary transport, we performed a base simulation with no changes in emissions and three sensitivity simulations, each with anthropogenic emissions removed in a different region, listed in Table 1. For the first sensitivity simulation (noQC), anthropogenic emissions are not permitted within the borders of Quebec. For the second simulation (noRoC), anthropogenic emissions are not permitted within the borders of Canada, except emissions coming from Quebec (Rest of Canada, RoC). For the third simulation (noUS), anthropogenic emissions are not permitted within the borders of the contiguous US. Emissions from shipping and aircraft are removed within the boundaries of the specified region, but we have not changed emissions from wildfires or natural sources. We remove the anthropogenic emissions in both the global 2° latitude and 2.5° longitude resolution simulation used to generate boundary conditions and the 0.5° latitude and 0.625° longitude higher-resolution simulations used for our analysis. We note that for the purposes of masking emissions, each cell in the model is considered to be entirely within a province or country. The resolution of the provincial and national emissions masks is the same as the resolution of the model. From the differences in the pollutant concentrations between the base simulation and the sensitivity simulations, we can calculate the fraction of the pollutant concentrations due to sources within Quebec, due to sources within the RoC and due to sources within the contiguous US.

Table 1 GEOS-Chem simulations performed

Name	Description
Base case	All emissions enabled
noQC	No anthropogenic emissions from within Quebec
noRoC	No anthropogenic emissions from within Canada, excluding Quebec
noUS	No anthropogenic emissions from within the contiguous United States

---

2.2 Creation of statistical models

We selected seven meteorological variables to use in creating our statistical models (i.e. the emulator): the velocities of the eastward (west-to-east) and northward (south-to-north) wind components at 850 hPa above the surface (U₈₅₀ and V₈₅₀), near-surface relative humidity (RH), sea level pressure (P₀), near-surface temperature (T), surface precipitation (P_r ) and the downward shortwave radiative flux at the surface (DSR). To capture the variability associated with short-term events, the daily averages of all variables considered in this study were used. We chose these variables based on five criteria:

1. They are available from the MERRA-2 meteorological data used to drive GEOS-Chem;

2. They would be expected to influence the Earth's surface concentrations of at least one of our target pollutants through physical mechanisms represented in GEOS-Chem;

3. They are commonly available as results of regional climate models;

4. They vary over time but are spatially two-dimensional; and

5. The cross-correlation between two selected variables was not high (R² < 0.8) in a preliminary data set.

In instances where multiple variables satisfied criteria 1–4, but had large cross-correlations, we used expert judgment to select the variable that was most likely to be physically relevant or most likely to be consistently defined within the MERRA-2 reanalysis, applied within GEOS-Chem, and defined within the regional climate models. For example, we chose the near-surface temperature instead of the surface skin temperature, as the former is more commonly used in studies of air pollution and climate. Due to the very large correlation between the daily means of these two metrics, our results would likely not have differed if we had chosen the surface skin temperature instead of the near-surface temperature.

We additionally excluded the planetary boundary layer height (PBLH) due to differences in the definition used between the MERRA-2 dataset and the regional climate model data used in this study: In MERRA-2, the PBLH is defined as the lowest model level at which the turbulent (heat) diffusivity falls below a prescribed threshold value.^36,37 In the Canadian Regional Climate Model, the PBLH is defined as the altitude at which the bulk Richardson number first exceeds a critical threshold of 0.25, marking the transition from a turbulent to a stably stratified regime.^38,39 The former therefore employs a diffusivity-based diagnostic, whereas the latter applies a stability-based diagnostic grounded in the Richardson-number criterion. These two approaches are likely to produce substantially different diagnosed values of the PBLH.

The statistical models relating meteorology to air pollutant concentrations were created in a manner similar to that described in Shen et al.:²⁵ for each spatial grid cell, the statistical model is a multiple linear regression of nine terms: the time series of the seven meteorological variables in the same grid cell and two additional terms derived using singular value decomposition (SVD) to reduce the three-dimensional large mesoscale meteorology (two horizontal dimensions and one temporal dimension) to two time-varying scalars. To avoid confounding seasonal or spatial variability with relationships between pollutant concentrations and meteorological variables, we created independent statistical models for each month and for each horizontal grid cell. This also permits the relationships between meteorology and pollutant concentrations to vary spatially and seasonally. We also created independent statistical models for each sensitivity simulation. For this study, we have therefore created 2 (O₃ and NO_X) × 12 months × 2912 spatial grid cells (56 latitudes × 52 longitudes) × 4 sensitivity simulations = 279552 statistical models that comprise the four emulators.

As stated, in addition to local meteorological variables at the central grid cell, the regression includes two predictors (S₁(t) and S₂(t)) which summarize large mesoscale meteorological patterns that influence pollutant concentrations not captured by local variables alone. Including them allows the model to account for both local meteorological effects and the larger-scale context in a low-dimensional way. We used SVD to reduce the three-dimensional large mesoscale meteorology to the two time-varying scalars S₁(t) and S₂(t) through the following process:

We normalized the time series of each of the seven meteorological variables in each grid cell by subtracting the temporal mean and dividing by the temporal standard deviation. We then considered the meteorological variables over a 9 × 9 grid cell area (4.5° latitude by 5.625° longitude) centred on the cell for which we were creating the statistical model. We calculated the correlation coefficient between the time series of each normalized meteorological variable in each grid cell and the time series of the pollutant concentrations in the central grid cell. These correlation coefficients can be considered as a matrix of dimension 7 (the number of meteorological variables) by 81 (the number of grid cells). We performed an SVD of this matrix of correlation coefficients, retaining only the first two SVD components. We excluded the central grid cell, because we later include the weather data from the central grid cell directly in the multilinear regression. The SVD has the form:

F = ULV^T (1)

where F is the matrix to be decomposed, U is a matrix of weights associated with the columns of F, V is a matrix of weights associated with the rows of F, and L is a diagonal matrix with non-negative numbers on the diagonal. The diagonal values in L represent the relative importance of each SVD component. In our study, we choose F to be the matrix of Pearson correlation coefficients, structured with each column corresponding to a spatial grid cell (81 columns, corresponding to the number of grid cells), and each row corresponding to a meteorological variable (7 rows, corresponding to the number of meteorological variables). U therefore represents the spatial weights of the decomposition (81 columns) and V represents the weights of the meteorological variables of the decomposition (7 columns). The number of rows in U, L, and V for the full decomposition is equal to the number of meteorological variables, but we retain only the first two.

As an example, we show in Fig. 2 the first two SVD components for April O₃ in the cell that contains the city of Montreal. The spatial weights of the first SVD component show little variation, and as we will show in Sect. 3.2, the first SVD component variable weights are similar to the correlations between O₃ and the meteorological variables in April. It is therefore likely that the first SVD component represents the autocorrelation between meteorological variables in the cell of interest and the surrounding cells, e.g., high temperatures in Montreal are associated with high temperatures in the area surrounding Montreal. The second component shows a different spatial pattern: more P_r to the northwest and less to the southeast, and lower values of DSR, T, and P₀ to the northwest and higher values to the southeast. This spatial pattern probably represents a storm to the northwest of Montreal in which cyclonic winds transport emissions of O₃ precursors from the west and southwest, particularly from the Great Lakes region.

Fig. 2 The variable weights (V) and the spatial weights (U) of the first two SVD components. Top row: first component. Bottom row: second component.

We then applied the inverse SVD process to the normalized time series of meteorological variables to calculate a new time series of the magnitude of each SVD mode:

S_k(t) = U_k^TM(t)V_k (2)

where U and V are, respectively, the spatial weights and meteorological variable weights matrices derived using eqn (1). The counter k refers to the column of U^T and V and the SVD mode, and has a value of 1 or 2. Additionally, M(t) is the matrix of normalized meteorological variables over the spatial area (dimension 7 by 81) as a function of time. The result of this process is S₁(t) and S₂(t), each of which is a scalar time series. These time series indicate how strongly each SVD-identified pattern is expressed at a given time step and whether the pattern is in a positive mode or a negative mode at that time step. Since the SVD was performed on the correlations of the two-dimensional large mesoscale meteorology with the pollutant concentrations in the grid cell of interest, these S₁(t) and S₂(t) terms will in turn correlate with the pollutant concentrations.

Finally, we performed a multiple linear regression of the concentrations of a given pollutant against the meteorological variables in the same grid cell and against the two SVD modes (S₁ and S₂). Thus, for each pollutant X, the statistical model is written as shown below:

[X](t) = a·U₈₅₀(t) + b·V₈₅₀(t) + c·RH(t) + d·P₀(t) + e·T(t) + f·P_r(t) + g·DSR(t) + h·S₁(t) + i·S₂(t) (3)

where a, b, c, d, e, f, g, h and i are the multiple linear regression coefficients. Since we are using the daily averages of each variable to create these statistical relationships, for a 30 day month, the number of data points for each variable used to create the regression will be 30 days 19 years = 570 data points.

We will refer to the set of statistical models as “the emulator” in the remainder of this article. As the statistical model for each grid cell requires GEOS-Chem training data for the surrounding region included in the SVD, the emulator domain is smaller than the high-resolution GEOS-Chem domain (Fig. 1). We note that the emulator reproduces relationships between meteorological variables and the pollutant concentrations, but does not distinguish the mechanisms behind these relationships. For example, the emulator may reproduce an increase in O₃ associated with an increase in T, but the analysis described here does not allow us to unambiguously distinguish temperature-driven increases in photochemistry, temperature-driven increases in biogenic VOC emissions, or increases in stagnation event frequency correlated with temperature.

2.3 Regional climate model data

To project future changes in pollutant concentrations, we applied the emulator to output from a regional climate model. We used output from version 5 of the canadian regional climate model (CRCM5),^40,41 which was developed by the Center for the Study and Simulation of Climate at Regional Scales (ESCER) at UQAM (Université du Québec à Montréal) with the collaboration of Environment and Climate Change Canada (ECCC). The CRCM5 v3.3.3.1 outputs were generated and provided by Ouranos. The simulation domain covers North America (AMNO22d2: 380 × 340 grid points including a 20-point sponge and halo zone surrounding the domain) with a horizontal grid-size mesh of 0.22° on a rotated latitudevlongitude grid (about 25 km resolution), using 10 minute time steps. This gives a free zone for analysis of 340 × 300 grid points.

We list in Table 2 the CRCM5 simulations used in this study. The regional simulations were performed using the projected evolution of greenhouse gas concentrations as defined by the IPCC RCP4.5 and RCP8.5 scenarios,⁴² as were the global driving models. For historical and climate change projections, each simulation was driven by atmospheric and oceanic 6-hourly fields taken from the first member (simulation r1i1p1) of each of the following four CMIP5 ESMs:

Table 2 The CRCM5 simulations used in this study

Simulation name	Driving model	Scenario	Temporal window
bby	CanESM2	Historical	1950-01 to 2005-12
bbz	CanESM2	RCP8.5	2005-01 to 2100-12
bca	CanESM2	RCP4.5	2006-01 to 2100-12
bcc	CNRM-CM5	Historical	1950-01 to 2005-12
bcd	CNRM-CM5	RCP8.5	2006-01 to 2100-12
bce	CNRM-CM5	RCP4.5	2006-01 to 2100-12
bcg	MPI-ESM-LR	Historical	1950-01 to 2005-12
bch	MPI-ESM-LR	RCP8.5	2006-01 to 2100-12
bdf	MPI-ESM-LR	RCP4.5	2006-01 to 2100-12
bcj	GFDL-ESM2M	Historical	1950-01 to 2005-12
bck	GFDL-ESM2M	RCP8.5	2006-01 to 2100-12
bcr	GFDL-ESM2M	RCP4.5	2006-01 to 2100-12

---

1. The Canadian Earth System Model version 2 (CanESM2/T63 corresponding approximately to 2.81° on a horizontal linear grid^43,44),

2. The Geophysical Fluid Dynamics Laboratory Earth System Model (GFDL-ESM2M/latitude v. longitude grid of 90 × 144 points⁴⁵),

3. The Centre National de Recherches Météorologiques Earth System Model (CNRM-CM5/atmospheric resolution of T127 ∼0.95°, archived on a latitude v. longitude grid of 128 × 256 points⁴⁶), and

4. The Max Planck Institute for Meteorology Earth System Model (MPI-ESM-LR/T063L40, ∼1.875°, latitude v. longitude grid of 96 × 192 points⁴⁷).

This was done by applying a smooth spectral nudging of large scales^48,49 to the horizontal wind component within the CRCM5 domain interior. The spectral nudging configuration consists of large-scale features being defined with a half-response wavelength of 1177 km and a relaxation time of 13.34 h. These large scales are imposed inside the CRCM5 domain and vary along the vertical: the nudging strength is set to zero from the surface to a height of 500 hPa and increases linearly to the top of the model's simulated atmosphere (10 hPa). The CRCM5 output was aggregated to match the MERRA2 grid using the xESMF Python package.⁵⁰

3 Results and discussion

3.1 GEOS-Chem simulations

We begin by comparing the results of the base simulation with the three sensitivity simulations (Table 1) for the historical period. We have already discussed in more detail in a previous publication²³ the effects on O₃ and NO_X of anthropogenic emissions from the US, Quebec and the RoC based on results from a similar set of simulations, so we will only provide an overview here.

As expected, the largest differences in NO_X concentrations are in the emission region set to zero in the sensitivity simulations, and the absolute concentration differences are larger in densely populated southern Canada than in sparsely populated northern Canada (Fig. A1 and A2). In addition, the differences in concentrations are generally greater in January (Fig. A1) than in July (Fig. A2). One factor contributing to this seasonal cycle is the seasonal cycle of boundary layer heights: due to warmer surface temperatures, boundary layer heights are greater in summer than in winter. The same mass of pollution is therefore diluted in a greater vertical layer in summer than in winter. Therefore, in the absence of seasonal differences in chemistry, emissions, or horizontal transport, we would expect higher concentrations of any locally emitted pollutant in winter than in summer.

Anthropogenic emissions reduce the O₃ concentrations in regions with high emissions in January due to titration by NO_X (Fig. 3), but they increase O₃ in July (Fig. 4). This indicates that the production of O₃ is VOC-limited in wintertime in the large urban areas of northeastern North America, and local emissions decrease O₃ by up to 5 ppb in Toronto and Montreal. It should be noted, however, that O₃ concentrations are generally much higher in summer than in winter and exceedances of the 8 hour O₃ standard are generally not observed in winter. In contrast, July O₃ appears to be NO_X-limited and reducing emissions to zero leads to a decrease in O₃.

Fig. 3 Differences in January O₃ surface concentrations between the sensitivity simulations and the base simulation. Top row: absolute differences. Bottom row: relative differences. The three columns show the results of the noQC (left), noRoC (centre) and noUS (right) sensitivity simulations. Negative values indicate that concentrations were lower in the base simulation than in the simulation with emissions removed.

Fig. 4 Differences in July O₃ surface concentrations between the sensitivity simulations and the base simulation. Top row: absolute differences. Bottom row: relative differences. The three columns show the results of the noQC (left), noRoC (centre) and noUS (right) sensitivity simulations.

Before examining the emulator, we wish to investigate how the chosen meteorological variables relate to the concentrations of O₃ and NO_X in the GEOS-Chem results. For this purpose, we show in Fig. 5 the annual-mean correlations between O₃ concentrations and the seven meteorological variables. We note that O₃ concentrations are correlated with higher temperatures for nearly all of the domain. As we described in Sect. 1, increases in O₃ concentrations with higher temperatures over polluted regions are consistent with other studies.^8–10,14,17 We also note positive correlations with V₈₅₀ and U₈₅₀ for the majority of the domain, which is associated with the transport of O₃ precursors from the southwest. The correlations with the other variables are weaker.

Fig. 5 Annual means of correlations between O₃ concentrations and the meteorological variables.

There is a strong seasonal dependence of the correlations between O₃ concentrations and meteorological variables, as revealed by histograms of the correlations between O₃ concentrations and meteorological variables for each month (Fig. A3), and the contrast between maps of the correlations between O₃ concentrations and meteorological variables for July (Fig. A4) and January (Fig. A5). The July correlations are more similar to those for the rest of the year than the January correlations. In winter, the correlations with most variables change sign around urbanized areas near the Great Lakes and the St. Lawrence River. For example, the correlation with T becomes negative while remaining positive in less-urbanized areas, and the correlation with P₀ becomes positive while remaining negative in many less-urbanized areas. For U₈₅₀ and V₈₅₀, it is likely that weaker winds lead to higher concentrations of NO_X which contributes to O₃ formation in summer and O₃ titration in winter. It is also likely that increases in RH and T increase the rate of destruction of O₃ in these urbanized locations in winter.

The correlations between the meteorological variables and NO_X do not depend on season so strongly as O₃ (Fig. A6). Therefore, we will only discuss the annual averages of the correlations between NO_X concentrations and the meteorological variables, shown in Fig. 6. The correlations between NO_X and V₈₅₀ or U₈₅₀ are likely due to the transport of NO_X and its precursors from urban source regions. NO_X concentrations are more strongly correlated with T, RH, and V₈₅₀, but these correlations are all negative south of ∼41° N. Racherla and Adams¹⁵ found previously that the response of simulated NO_X concentrations to climate change depended on existing O₃ concentrations, with NO_X increasing due to climate change under low-O₃ conditions and decreasing under high-O₃ conditions. They explained the increases as due to an increase in PAN decomposition, and the decreases due to increased loss of NO_X to HNO₃, which is in turn due to an increase in the NO₂ : NO ratio from increased O₃ chemical formation. We note that the region of our domain with negative correlations between NO_X concentrations and T, RH, and V₈₅₀ corresponds roughly with the region where temporal-mean O₃ concentrations are greater than 35 ppb. It is therefore likely that the correlations with T and RH are primarily due to decreases in PAN lifetimes in the northern part of the domain and to increased loss of NO_X to HNO₃ in the southern part of the domain. It is also likely that the geographical distinction between these two regimes is sensitive to changes in anthropogenic emissions, which were held fixed in our study at 2010 levels.

Fig. 6 Annual means of correlations between NO_X concentrations and the meteorological variables.

3.2 Validation of the emulator and uncertainty analysis

We note that the emulator is not intended to exactly replicate pollutant concentrations, as we have intentionally not included any factors in the emulator for changes in natural or biogenic emissions over time. Insofar as these emissions are well-correlated with the meteorological parameters used in the emulator, the effects of changes in meteorology on the emissions will implicitly be included in the emulator, although our experimental design does not allow us to distinguish the effects of these changes in emissions from other effects of the meteorology on the pollutant concentrations. Additionally, even where meteorology has a direct effect on emissions, the effect of transport will be to further weaken correlations between emissions and pollutant concentrations, as there will be a delay in time and a displacement in space between the meteorology that drives the emissions and the meteorology used to derive the emulator. Some emissions (e.g. wildfires) will not be well correlated with meteorological factors over the simulated period, and this variability cannot be captured by the emulator. Nevertheless, we find that the emulator captures much of the variability in the O₃ and NO_X concentrations.

As mentioned in Sect. 2.2, we create individual emulators for each month. Due to the difficulties in performing calculations on the entire data set at the same time due to the large amount of data, we evaluate each monthly emulator independently and summarize the goodness-of-fit metrics using the averages of the monthly evaluations. We use Pearson's correlation coefficient (R), mean error (ME), and normalized mean error (NME) as goodness-of-fit metrics. We show the spatial and temporal averages of these metrics in Table 3, and we show the annual means of these metrics for O₃ in Fig. 7 and for NO_X in Fig. A7.

Table 3 Spatial and temporal averages of goodness-of-fit metrics

Species	Mean R	Mean ME [ppb]	Mean NME [%]
O3	0.62	2.95	11
NOX	0.60	0.31	47

---

Fig. 7 Annual means of goodness-of-fit metrics between the emulator and the training data for O₃. Left: correlation coefficient. Centre: mean error. Right: normalized mean error.

For O₃, the annual averages of R are >0.6 for the majority of eastern Canada, and the annual averages of NME are <15% for nearly all of eastern Canada. Correlations decrease further away from source regions. In these regions, O₃ concentrations are more strongly influenced by long-range transport, and the emulator cannot fully represent the variability in long-range transport.

For NO_X, in general, ME values are higher at higher concentrations, and NME values are lower at higher concentrations. This pattern is expected because small relative errors in high concentrations yield large absolute errors, and small absolute errors in low concentrations yield large relative errors. Correlations between the NO_X emulator results and the GEOS-Chem results are highest around the Great Lakes and north of the St. Lawrence River.

In order to evaluate the uncertainties due to the choice of training data for the emulator, we recreated the emulator according to the process described above 20 times, each time excluding one year of training data. Since the goal of the emulator is to make projections of changes in pollutant concentrations due to future climate change, we focus on the variability of relationships between the emulator output and the input weather variables. For a simple multilinear regression, this would normally be represented by the variability of the slopes (the coefficients a, b, c, … in the eqn (3)). However, due to the spatial autocorrelation in meteorology that is frequently captured by the first SVD component, as shown in Fig. 2, the value of g is not independent of other slopes of the eqn (3), and a naive interpretation of the variability of the slopes would overestimate the uncertainty of the emulator.

Instead, we will examine the partial derivative of the emulator results with respect to spatially homogeneous changes in weather variables. For example, the partial derivative of the eqn (3) with respect to a spatially homogeneous change in RH would be:

where U₁ and U₂ are the spatial weights of the first and second SVD components, and V_1RH and V_2RH are the variable weights corresponding to RH of the first and second SVD components. We will evaluate the uncertainty in the emulator as the variability of these partial derivatives when we recreate the emulator excluding a single year from the input data. This approach will properly capture the variability in the sensitivity of the emulator's predicted pollutant concentration to meteorology in the same grid cell. However, if the spatial weights of one of the SVD components frequently total zero or near zero, the variability of that component will not be captured by these indicators.

We show in Table 4 the spatial and temporal averages of the standard deviation of the partial derivatives of O₃ with respect to the seven meteorological variables. We show analogous statistics for NO_X in Table 4. For comparison, we also include the spatial and temporal averages of the absolute values of the partial derivatives of O₃ or NO_X with respect to the seven meteorological variables, when calculated using the full set of training data. The means of the standard deviations are at most 17% of the means of the absolute values, and they are an order of magnitude smaller in most cases. For another perspective, we calculate another metric similar to the coefficient of variation:

Table 4 Statistics describing the partial derivatives of O₃ and NO_X: mean of absolute values, mean of standard deviation, and median of normalized standard deviation

Species	Variable	Mean of absolute values	Mean of standard deviation	Median of normalized standard deviation
O3	V850	0.120 ppb/m s^−1	0.0109 ppb/m s^−1	12%
	U850	0.123 ppb/m s^−1	0.0106 ppb/m s^−1	11%
	RH	0.0992 ppb/%	0.00868 ppb/%	10%
	P0	0.0582 ppb/hPa	0.00881 ppb/hPa	21%
	T	0.636 ppb/K	0.0322 ppb/K	5%
	Pr	0.0812 ppb/mm day^−1	0.0130 ppb/mm day^−1	18%
	DSR	0.0190 ppb/W m^−2	0.00227 ppb/W m^−2	18%
NOX	V850	0.0157 ppb/m s^−1	0.00115 ppb/m s^−1	7%
	U850	0.0164 ppb/m s^−1	0.00106 ppb/m s^−1	6%
	RH	0.0153 ppb/%	0.000756 ppb/%	6%
	P0	0.00486 ppb/hPa	0.000830 ppb/hPa	18%
	T	0.0204 ppb/K	0.00196 ppb/K	9%
	Pr	0.0118 ppb/mm day^−1	0.00128 ppb/mm day^−1	17%
	DSR	0.00171 ppb/W m^−2	0.000210 ppb/W m^−2	14%

---

1. First, we normalize the standard deviations by the absolute values of the partial derivatives calculated using the full set of training data for each grid cell and each month.

2. Next, we calculate the medians. We choose medians, because there are partial derivative values close to zero, which give very large values to the normalized standard deviations, and this makes the means of the normalized standard deviations misleading. We note that the largest median of the normalized standard deviations for O₃ has a value of 21%. The value for T, the variable that has the most significant changes due to climate change and has the strongest correlations with O₃, is only 5%. The results for NO_X are similar, with the medians of the normalized standard deviations ranging from 6% for U₈₅₀ and RH to 18% for P₀.

3.3 Future projections based on CRCM5 data

As described in Sect. 2.3, we apply the emulator to the results of four versions of CRCM5, each having been used to simulate two GHG (greenhouse gas) scenarios: RCP8.5 and RCP4.5. In addition to the two future scenarios, the emulator is applied to simulations of the historical period. We list these simulations in Table 2. This section describes our projections of changes in pollutant concentrations caused by climate change between the periods of 1986 to 2006 and 2080 to 2100. We choose periods of 20 years to reduce the effect of interannual variability. We report the differences in concentrations between these two periods rather than the absolute concentrations provided by the emulator, since the differences are a more robust result. Only the means of the four CRCM5 simulations are shown in this section, but the agreement between the four simulations is indicated using the following notation: Vertical hatching indicates that the standard deviation among the CRCM5 simulations is larger than twice the mean of the change obtained from the CRCM5 ensemble (i.e. the four simulations have low agreement on the magnitude of the change), and diagonal hatching indicates that the simulations are evenly distributed according to the sign of the change (two positive, two negative). Additionally, we use horizontal hatching to indicate that the change is insignificant (p > 0.05) in more than one simulation according to Welch's t-test.

In both the RCP8.5 (Fig. 8) and the RCP4.5 (Fig. A8) scenarios, changes due to climate in V₈₅₀ and U₈₅₀ are not significant for the majority of our domain according to the four simulations used. There is a robust increase in temperature which is greater in the north for both scenarios, about 5 K for RCP8.5 and about 3 K for RCP4.5. There are increases in P_r that are significant in most of eastern Canada of 0.5 to 1.0 mm day⁻¹ in the RCP8.5 scenario. We project significant increases over the Atlantic Ocean in the southeast of our domain and significant decreases over Hudson's Bay and northern Quebec in the northeast of our domain for P₀ (+0.5 hPa and −2 hPa in RCP8.5) and DSR (+5 W m⁻² and -15 W m⁻² in RCP8.5). The changes in DSR are associated with opposing changes in cloud cover. We project increases in RH of up to 2% in RCP8.5 over the Atlantic Ocean and portions of Canada and the US east of the Great Lakes. Three of the four CRCM5 simulations project decreases in RH over northeastern Canada, including Newfoundland and Labrador and much of Ontario and Quebec, but the simulations do not agree on the magnitude of these decreases. For all of the meteorological variables, the spatial patterns are similar in RCP8.5 and RCP4.5, but the magnitudes of the changes are roughly a factor of two smaller in the RCP4.5 scenario.

Fig. 8 CRCM5 ensemble-mean projected changes in meteorological variables for the RCP8.5 scenario. Hatching is described in Sect. 3.3.

By applying the emulator to the CRCM5 results, we project how these changes in meteorology will affect spring (March, April, and May) O₃ concentrations in the RCP8.5 (Fig. 9) and RCP4.5 (Fig. 10) scenarios. We focus on the spring, because O₃ shows maximum monthly concentrations in spring. Spring O₃ concentrations are projected to increase by about 4 ppb in Quebec in the RCP8.5 scenario, but only by about 2 ppb in the RCP4.5 scenario. We also show in Fig. 9 and 10 the changes in O₃ due to each meteorological variable in isolation. As O₃ is best correlated with T, and T has the most significant projected changes of all the meteorological variables, it is not surprising that the increases in O₃ concentrations are largely explained by the increases in T. Nevertheless, we note that there are also small mitigating contributions from the decrease in DSR and the increase in P_r, and small contributions correlated with changes in RH.

Fig. 9 CRCM5 ensemble-mean projected changes in Spring (March, April, May) O₃ concentrations for the RCP8.5 scenario due to each meteorological variable in isolation. Total projected change is shown in the lower right. Hatching is described in Sect. 3.3.

Fig. 10 CRCM5 ensemble-mean projected changes in Spring (March, April, May) O₃ concentrations for the RCP4.5 scenario due to each meteorological variable in isolation. Total projected change is shown in the lower right. Hatching is described in Sect. 3.3.

We also apply the emulator to the CRCM5 results in project changes in annual-mean NO_X concentrations in the RCP8.5 (Fig. 11) and RCP4.5 (Fig. A9) scenarios. Smaller changes in NO_X concentrations are projected compared to the increases in O₃ concentrations. We project greater changes in NO_X concentrations in urbanized environments. We project that NO_X concentrations will increase in most of eastern Canada by less than 0.5 ppb, with increases in the Greater Toronto Area of more than 1 ppb. We also show in Fig. 11 and A9 the changes in NO_X due to each meteorological variable in isolation. Consistent with the correlations shown in Fig. 6, NO_X concentrations in the southern part of our domain decrease associated with increases in T, likely due to an increase in loss to HNO₃. This is partially offset by increases in concentrations associated with increases in RH. Therefore, we project decreases in NO_X concentrations in the southern part of our domain.

Fig. 11 CRCM5 ensemble-mean projected changes in annual mean NO_X concentrations for the RCP8.5 scenario due to each meteorological variable in isolation. Total projected change is shown in the lower right. Hatching is described in Sect. 3.3.

We also created emulators based on the results of the three GEOS-Chem sensitivity simulations without anthropogenic emissions from the three regions: noQC, noRoC, and noUS. We applied these emulators to the CRCM5 ensemble output. To assess future changes in transboundary transport, we use the emulators for each scenario to predict spring O₃ and annual-mean NO_X for each scenario for each of the years 1986 to 2006 and 2080 to 2100. We calculate the change between these two periods for the full emulator domain (36.75 N to 62.75 N, 87.5 W to 52.5 W, shown in Fig. 1) and the parts of the domain covering the province of Quebec, the rest of Canada, and the United States. As listed in Table 5, spring O₃ concentrations within Quebec are projected to increase by 3.36 ppb in the base-case RCP8.5 scenario due to climate change. Without emissions from either Quebec or the RoC, this increase is slightly smaller at 3.23 or 3.17 ppb. The increase is much smaller without US emissions, at 2.42 ppb. We therefore project that US emissions will contribute 0.94 ppb more to Quebec-mean O₃ concentrations in 2080 to 2100 under an RCP8.5 scenario due solely to changes in climate. Similarly to O₃, projected increases in annual-mean NO_X concentrations averaged across Quebec are smallest for the noUS scenario, indicating that US emissions are projected to contribute 0.0413 ppb more to NO_X concentrations due to climate change, while these increases are 0.0063 and 0.0112 ppb for emissions from Quebec and the RoC, respectively. In all cases, the projected increases for the RCP4.5 climate scenario for each emissions scenario are roughly half of those for the RCP8.5 scenario.

Table 5 Projected changes in Spring (March, April, May) O₃ and annual-mean NO_X concentrations in ppb between 1986 to 2006 and 2080 to 2100 for each climate scenario and emissions sensitivity case. Results are averaged across the full emulator domain (36.75 N to 62.75 N, 87.5 W to 52.5 W, shown in Fig. 1) and the parts of the domain covering the province of Quebec, the rest of Canada, and the United States

Species	Climate	Emissions	Domain mean	Within QC	Within RoC	Within US
MAM O3	RCP8.5	Base	3.60	3.36	3.64	4.75
		noQC	3.55	3.23	3.57	4.80
		noCA	3.53	3.17	3.48	4.88
		noUS	2.20	2.42	2.59	1.26
	RCP4.5	Base	1.74	1.50	1.70	2.42
		noQC	1.72	1.45	1.67	2.44
		noCA	1.72	1.42	1.63	2.48
		noUS	1.05	1.07	1.20	0.66
Annual-mean NOX	RCP8.5	Base	0.0507	0.0763	0.0804	0.0524
		noQC	0.0503	0.0700	0.0808	0.0566
		noCA	0.0456	0.0650	0.0666	0.0511
		noUS	0.0252	0.0349	0.0392	0.0437
	RCP4.5	Base	0.0243	0.0385	0.0390	0.0233
		noQC	0.0238	0.0348	0.0392	0.0248
		noCA	0.0219	0.0329	0.0326	0.0220
		noUS	0.0116	0.0174	0.0185	0.0207

---

To investigate these future changes in transboundary transport in more detail, we calculate the differences between the results of the emulators for the base scenario and each of the noQC, noRoC, and noUS scenarios in order to estimate the contribution of each region to spring O₃ concentrations. This calculation is done for each of the years 1986 to 2006 and 2080 to 2100. Then, we calculate the changes in these contributions between the two periods. In this way, we evaluate future changes in the contributions of each region to air pollution in eastern Canada. We will discuss first the results for the RCP8.5 GHG scenario, and then the results for the RCP4.5 scenario.

We begin with the CRCM5 ensemble-mean changes in the absolute contributions of the different regions (Quebec, RoC, and US) to the spring O₃ concentrations (Fig. 12). We project that future increases in O₃ concentrations in eastern Canada are primarily driven by U.S. emissions, due to the larger anthropogenic emissions of O₃ precursors from the U.S. region. Changes in the amount of O₃ due to emissions from Quebec and RoC are less than 1 ppb, and we even project decreases in the contributions from Quebec and the RoC to O₃ to the south and southeast of these emission regions. As the CRCM5 results do not show a significant change in U₈₅₀ or V₈₅₀ (Fig. 8), the changes in contributions (including the decreases) do not appear to be due to a change in transport patterns. To investigate this further, we have calculated the effect of changing each meteorological variable independently on the contributions to O₃ concentrations from each region. As shown in Fig. A10, the changes in the contributions are due almost entirely to the changes in T. As T increases throughout our domain, this implies that in the southern part of our domain, the increase in O₃ due to T in the base case emulator is less than in the noQC and noROC emulators. In other words, the relationship between O₃ and T in this region is less positive in the base-case GEOS-Chem simulation than in the noQC or noRoC simulations. One possible cause for this could be that the same transport patterns that bring Canadian emissions southward also typically bring colder air from the north, decreasing T. While we have attempted to control for cross-correlation between transport direction and T by including U₈₅₀ and V₈₅₀ in our emulator, it is still possible that the correlation between T and transport direction causes the emulator to overestimate the relationship between O₃ and T in the absence of Quebec or RoC emissions. If so, then this same effect may also be contributing to an overestimate of the effect of emissions from each of the three regions on O₃ concentrations to the north of the emissions regions, including the influence of US emissions on eastern Canada.

Fig. 12 Spring (March, April, May) means of CRCM5 ensemble-mean changes in the contributions of different regions to the O₃ concentrations. The GHG scenario is RCP8.5. Hatching is described in Sect. 3.3. Top: changes in absolute contributions. Bottom: changes in relative contributions. Left: contribution from Quebec (base – noQC). Centre: contribution from the RoC (base – noRoC). Right: contribution from the US (base – noUS).

We use the same process as above to investigate the changes in the relative contributions from each region (e.g. (base – noQC)/base) and calculate the changes due to climate change in the proportion of pollutant concentrations due to emissions from each region. In Fig. 12, we show the springtime means of the changes in the relative contributions to O₃ concentrations. Changes in relative contributions are mostly consistent with changes in absolute contributions, but negative changes extend further north for the noQC and noRoC cases. This result can be explained by the following reasoning: If Quebec or RoC emissions are not responsible for more O₃, and US emissions are responsible for more O₃, it follows that emissions from Quebec or the RoC are responsible for a smaller fraction of the O₃. We also note that all maps show positive values in northern Ontario, Quebec, and Labrador. Here, the proportion of O₃ due to anthropogenic emissions from each of the three regions increases, as the proportion due to natural sources or due to transport from regions outside North America decreases.

We apply the same analysis to annual-mean NO_X concentrations (Fig. 13). Increases in NO_X concentrations due to anthropogenic emissions from Quebec are <0.1 ppb everywhere in our domain. Increases due to RoC emissions are <0.1 ppb outside of the Greater Toronto Area. We project that the contribution of US emissions to NO_X concentrations in eastern Canada will increase by up to 0.2 ppb. As with O₃, the proportion of NO_X concentrations due to US emissions increases throughout most of eastern Canada, and the proportion due to Quebec and RoC emissions frequently decreases as the absolute US contribution increases more than the absolute Quebec or RoC contributions. We also project that anthropogenic US emissions will contribute less to NO_X concentrations south of ∼40° N, likely due to increases in removal via HNO₃.

Fig. 13 Annual means of CRCM5 ensemble-mean changes in the contributions of different regions to the NO_X concentrations for the RCP8.5 scenario. Hatching is described in Sect. 3.3. Top: changes in absolute contributions. Bottom: changes in relative contributions. Left: contribution from Quebec (base – noQC). Centre: contribution from the RoC (base – noRoC). Right: contribution from the US (base – noUS).

Lastly, we apply this analysis of the changes in regional contributions to the RCP4.5 scenario. The RCP4.5 results have the same spatial patterns as for the RCP8.5 scenario, but the magnitude of the changes in the contributions is about half of the magnitude in the RCP8.5 scenario, as shown for spring O₃ concentrations in Fig. A12 and for annual-mean NO_X concentrations in Fig. A13.

4 Conclusions

The increases in NO_X near emitting regions that are projected here are consistent with previous studies using different methods.¹³ Increases in O₃ near emitting regions are also consistent with other studies.^8,9,19 We do not project decreases in O₃, even far from the transmitting regions. We note, however, that our projections are sensitive to the anthropogenic emissions that we used in our training data, fixed at 2010 levels. Other recent studies do not project unambiguous increases in O₃ concentrations throughout eastern Canada.^17,18,51 An explanation for this difference is that the emissions scenarios used in these other studies project a decrease in anthropogenic emissions of O₃ precursors in the future, and lower concentrations may weaken the relationship between O₃ and temperature, or even reverse it. As the relationship between NO_X and temperature is also dependent on O₃ concentrations, we also expect the projections of future NO_X changes shown here to be sensitive to the assumptions regarding anthropogenic emissions.

According to our projections, the origin of most of the anthropogenic precursors giving rise to our projected changes in NO_X concentrations will be local emissions, while increases in O₃ concentrations will be due primarily due to precursors from the US. This is not due to a change in transport direction, as we do not find a clear signal in changes in U₈₅₀ or V₈₅₀ due to climate change. Instead, the increases in temperature increase the efficiency of O₃ formation from anthropogenic precursors, thus enhance pre-existing disparities due to the different quantities of O₃ precursor emissions from each region. If this result is robust to other regions globally, it would suggest that climate change has the potential to exacerbate the export of pollution across political boundaries. However, we project decreases in the contribution of the RoC and QC to O₃ concentrations south of the Canada-US border that are difficult to mechanistically explain as being caused by an increase in temperature. While we have attempted to control for cross-correlation between transport direction and temperature by including wind velocity variables in our emulator, it is still possible that the correlation between temperature and transport direction (northerly winds bringing both colder air and Canadian emissions into the US) is contributing to the emulator projecting less influence of Canadian emissions on US O₃ concentrations when temperatures increase. If so, then this same effect may also be contributing to an overestimate of the effect of emissions from each of the three regions on O₃ concentrations to the north of the emissions regions, including the influence of US emissions on eastern Canada. Due to the nonlinear relationships between concentrations and precursor emissions, we expect these results to be sensitive to assumptions regarding the anthropogenic emissions. We therefore encourage further work investigating the interactions between climate change, the export of air pollution across political boundaries, and assumptions regarding the future evolution of anthropogenic emissions.

The O₃ climate penalty means that mitigation policy has to work harder to achieve the same result. In addition, for Canada, international cooperation will be necessary to address the problem if most of the projected increase in O₃ concentrations is associated with emissions occurring in the United States. The Canadian Ambient Air Quality Standards (CAAQS) are set by the Canadian Council of Ministers of the Environment, and the current standard for O₃ is 60 ppb (calculated from the 3 year average of the annual 4th highest of the daily maximum 8 hour average O₃ concentration).⁵² Observations from monitoring stations in eastern Canada, including its two largest cities Toronto and Montréal⁵³ show that ambient levels are often close to or exceed this standard. For example, the City of Montréal reported for the metropolitan region an observed concentration of 58 pbb for 2021–2023,⁵⁴ calculated as described above. It is therefore expected that the number of O₃ exceedances will increase due to the projected changes in meteorological parameters under the RCP8.5 and RCP4.5 scenarios.

The Canada-US Air Quality Agreement is the current environmental treaty between the two countries that addresses the transboundary transport of air pollutants, specifically of SO₂, NO_X, and VOCs.⁵⁵ The most recent review and assessment of the treaty in 2023 concluded that Canada and the US have met their commitments with respect to ground-level O₃, but also that it continues to have significant impacts on public health and agricultural production.⁵⁶ The results of our current study suggest that the O₃ climate penalty will make meeting the CAAQS standard increasingly difficult in the future, and that the Canada-U.S. Air Quality Agreement may need to be updated with more stringent specific objectives concerning ground-level O₃ precursors in order to ensure the continued success of the agreement.

A further question beyond the scope of this work is how the effects of a given feasible policy intervention, such as a 10% decrease in NO_X emissions from a particular region, would vary with climate change. We chose to remove all anthropogenic emissions of all species in each region in order to attribute the full contribution of each region to pollutant concentrations, and we note that due to the nonlinear relationships between pollutant concentrations and precursor emissions, we do not expect that e.g. the effect of reducing the emissions by 10% will be 10% of the effect of reducing the emissions by 100%. This is true for NO_X, but particularly problematic for O₃, where concentrations do not necessarily increase for a given increase in NO_X or VOC emissions. The results shown here would suggest that the effects on O₃ concentrations of a given policy intervention will be sensitive to climate change, but this should be investigated further in future work.

Author contributions

RS performed the analysis, created the visualizations, and wrote the original manuscript. PH was responsible for conceptualization and funding acquisition. RS, PH, HC, and BM all contributed to methodology development. TS procured and aggregated the CRCM5 output. All authors contributed to review and editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The GEOS-Chem output used in this study and the fit parameters that comprise the emulators and are available on the Borealis dataverse repository at https://doi.org/10.5683/SP3/X8ONIC The code for GEOS-Chem can be found at https://zenodo.org/records/5764874 with doi:https://doi.org/10.5281/zenodo.5764874. The version of the code employed for this study is version 13.3.4. The CRCM5 output is available upon request by contacting simulations_ouranos@ouranos.ca.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ea00112a.

Acknowledgements

The CRCM5 computations by Ouranos and the GEOS-Chem calculations were made on the supercomputers beluga and narval managed by Calcul Québec and the Digital Research Alliance of Canada (alliancecan.ca). The operation of this supercomputer received financial support from Innovation, Science and Economic Development Canada and the Ministère de l’Économie et de l’Innovation du Québec. This work was supported by the Ministère de l’Environnement, de la Lutte contre les changements climatiques, de la Faune et des Parcs du Québec. We also thank NSERC and CFCAS for the funding of the development of the CRCM5. We thank Charmaine Holloway for assistance with creating Fig. 1. We thank Dr Cynthia H. Whaley for helpful comments on a draft of this manuscript.

Notes and references

1. D. M. Stieb, S. Judek and R. T. Burnett, J. Air Waste Manage. Assoc., 2002, 52, 470–484.
2. C. A. Pope and D. W. Dockery, J. Air Waste Manag. Assoc., 2006, 56, 709–742.
3. F. Laden, J. Schwartz, F. E. Speizer and D. W. Dockery, Am. J. Respir. Crit. Care Med., 2006, 173, 667–672.
4. D. R. Reidmiller, A. M. Fiore, D. A. Jaffe, D. Bergmann, C. Cuvelier, F. J. Dentener, B. N. Duncan, G. Folberth, M. Gauss, S. Gong, P. Hess, J. E. Jonson, T. Keating, A. Lupu, E. Marmer, R. Park, M. G. Schultz, D. T. Shindell, S. Szopa, M. G. Vivanco, O. Wild and A. Zuber, Atmos. Chem. Phys., 2009, 9, 5027–5042.
5. S. C. Anenberg, J. J. West, H. Yu, M. Chin, M. Schulz, D. Bergmann, I. Bey, H. Bian, T. Diehl, A. Fiore, P. Hess, E. Marmer, V. Montanaro, R. Park, D. Shindell, T. Takemura and F. Dentener, Air Qual. Atmos. Health, 2014, 7, 369–379.
6. C.-K. Liang, J. J. West, R. A. Silva, H. Bian, M. Chin, Y. Davila, F. J. Dentener, L. Emmons, J. Flemming, G. Folberth, D. Henze, U. Im, J. E. Jonson, T. J. Keating, T. Kucsera, A. Lenzen, M. Lin, M. T. Lund, X. Pan, R. J. Park, R. B. Pierce, T. Sekiya, K. Sudo and T. Takemura, Atmos. Chem. Phys., 2018, 18, 10497–10520.
7. Hemispheric transport of air pollution 2010: prepared by the Task Force on Hemispheric Transport of Air Pollution acting within the framework of the Convention on Long-range Transboundary Air Pollution, ed. F. Dentener, T. Keating, H. Akimoto, N. Pirrone, S. Dutchak, A. Zuber, C. On Long-Range Transboundary Air Pollution, U. Nations and U. T. F. on Emission Inventories and Projections, United Nations, New York ; Geneva, 2010.
8. E. von Schneidemesser, P. S. Monks, J. D. Allan, L. Bruhwiler, P. Forster, D. Fowler, A. Lauer, W. T. Morgan, P. Paasonen, M. Righi, K. Sindelarova and M. A. Sutton, Chem. Rev., 2015, 115, 3856–3897.
9. S. Szopa, V. Naik, B. Adhikary, P. Artaxo, T. Berntsen, W. D. Collins, S. Fuzzi, L. Gallardo, A. Kiendler-Scharr, Z. Klimont, H. Liao, N. Unger and P. Zanis, in Climate Change 2021: the Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, ed. V. Masson-Delmotte, P. Zhai, A. Pirani, S. L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M. I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J. B. R. Matthews, T. K. Maycock, T. Waterfield, O. Yelekçi, R. Yu and B. Zhou, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2021, pp. 817–922.
10. X. Lu, L. Zhang and L. Shen, Curr. Pollut. Rep., 2019, 5, 238–260.
11. W. C. Porter and C. L. Heald, Atmos. Chem. Phys., 2019, 19, 13367–13381.
12. W. R. Burrows, B. Kochtubajda and G. Fricska, Atmos.-Ocean, 2025, 63, 125–146.
13. R. M. Doherty, M. R. Heal and F. M. O'Connor, Environ. Health, 2017, 16, 118.
14. R. M. Doherty, O. Wild, D. T. Shindell, G. Zeng, I. A. MacKenzie, W. J. Collins, A. M. Fiore, D. S. Stevenson, F. J. Dentener, M. G. Schultz, P. Hess, R. G. Derwent and T. J. Keating, J. Geophys. Res. Atmos., 2013, 118, 3744–3763.
15. P. N. Racherla and P. J. Adams, Atmos. Chem. Phys., 2008, 8, 871–885.
16. A. B. Guenther, X. Jiang, C. L. Heald, T. Sakulyanontvittaya, T. Duhl, L. K. Emmons and X. Wang, Geosci. Model Dev. (GMD), 2012, 5, 1471–1492.
17. S. T. Turnock, R. Allen, A. T. Archibald, M. Dalvi, G. Folberth, P. T. Griffiths, J. Keeble, E. Robertson and F. M. O'Connor, Earths Future, 2022, 10, e2022EF002687.
18. S. T. Turnock, R. J. Allen, M. Andrews, S. E. Bauer, M. Deushi, L. Emmons, P. Good, L. Horowitz, J. G. John, M. Michou, P. Nabat, V. Naik, D. Neubauer, F. M. O'Connor, D. Olivié, N. Oshima, M. Schulz, A. Sellar, S. Shim, T. Takemura, S. Tilmes, K. Tsigaridis, T. Wu and J. Zhang, Atmos. Chem. Phys., 2020, 20, 14547–14579.
19. J. Kelly, P. A. Makar and D. A. Plummer, Atmos. Chem. Phys., 2012, 12, 5367–5390.
20. R. M. Doherty, C. Orbe, G. Zeng, D. A. Plummer, M. J. Prather, O. Wild, M. Lin, D. T. Shindell and I. A. Mackenzie, Atmos. Chem. Phys., 2017, 17, 14219–14237.
21. K. Murazaki and P. Hess, J. Geophys. Res. Atmos., 2006, 111, 2005JD005873.
22. E. A. Barnes and A. M. Fiore, Geophys. Res. Lett., 2013, 40, 2839–2844.
23. R. Stevens, C. Poterlot, N. Trieu, H. A. Rodriguez and P. L. Hayes, Environ. Sci. Adv., 2024, 3, 448–469.
24. N. Fakhri, R. Stevens, A. Downey, K. Oikonomou, J. Sciare, C. Afif and P. L. Hayes, Atmos. Chem. Phys., 2024, 24, 1193–1212.
25. L. Shen, L. J. Mickley and L. T. Murray, Atmos. Chem. Phys., 2017, 17, 4355–4367.
26. I. Bey, D. J. Jacob, R. M. Yantosca, J. A. Logan, B. D. Field, A. M. Fiore, Q. Li, H. Y. Liu, L. J. Mickley and M. G. Schultz, J. Geophys. Res. Atmos., 2001, 106, 23073–23095.
27. R. J. Park, D. J. Jacob, B. D. Field, R. M. Yantosca and M. Chin, J. Geophys. Res. Atmos., 2004, 109, D15204.
28. National Emissions Inventory Collaborative, 2016v1 Platform, US EPA, 2019, https://www.epa.gov/air-emissions-modeling/2016v1-platform.
29. R. M. Hoesly, S. J. Smith, L. Feng, Z. Klimont, G. Janssens-Maenhout, T. Pitkanen, J. J. Seibert, L. Vu, R. J. Andres, R. M. Bolt, T. C. Bond, L. Dawidowski, N. Kholod, J.-i. Kurokawa, M. Li, L. Liu, Z. Lu, M. C. P. Moura, P. R. O'Rourke and Q. Zhang, Geosci. Model Dev. (GMD), 2018, 11, 369–408.
30. Z. A. Tzompa-Sosa, E. Mahieu, B. Franco, C. A. Keller, A. J. Turner, D. Helmig, A. Fried, D. Richter, P. Weibring, J. Walega, T. I. Yacovitch, S. C. Herndon, D. R. Blake, F. Hase, J. W. Hannigan, S. Conway, K. Strong, M. Schneider and E. V. Fischer, J. Geophys. Res. Atmos., 2017, 122, 2493–2512.
31. Y. Xiao, J. A. Logan, D. J. Jacob, R. C. Hudman, R. Yantosca and D. R. Blake, J. Geophys. Res. Atmos., 2008, 113, 2007JD009415.
32. M. E. J. Stettler, S. Eastham and S. R. H. Barrett, Atmos. Environ., 2011, 45, 5415–5424.
33. N. W. Simone, M. E. J. Stettler and S. R. H. Barrett, Transport. Res. Transport Environ., 2013, 25, 33–41.
34. C. A. Keller, M. S. Long, R. M. Yantosca, A. M. Da Silva, S. Pawson and D. J. Jacob, Geosci. Model Dev. (GMD), 2014, 7, 1409–1417.
35. M. J. E. van Marle, S. Kloster, B. I. Magi, J. R. Marlon, A.-L. Daniau, R. D. Field, A. Arneth, M. Forrest, S. Hantson, N. M. Kehrwald, W. Knorr, G. Lasslop, F. Li, S. Mangeon, C. Yue, J. W. Kaiser and G. R. van der Werf, Geosci. Model Dev. (GMD), 2017, 10, 3329–3357.
36. E. L. McGrath-Spangler and A. S. Denning, J. Geophys. Res. Atmos., 2012, 117, 2012JD017615.
37. J. Guo, J. Zhang, K. Yang, H. Liao, S. Zhang, K. Huang, Y. Lv, J. Shao, T. Yu, B. Tong, J. Li, T. Su, S. H. L. Yim, A. Stoffelen, P. Zhai and X. Xu, Atmos. Chem. Phys., 2021, 21, 17079–17097.
38. D. J. Seidel, Y. Zhang, A. Beljaars, J. Golaz, A. R. Jacobson and B. Medeiros, J. Geophys. Res. Atmos., 2012, 117, 2012JD018143.
39. R. McTaggart-Cowan, P. A. Vaillancourt, A. Zadra, S. Chamberland, M. Charron, S. Corvec, J. A. Milbrandt, D. Paquin-Ricard, A. Patoine, M. Roch, L. Separovic and J. Yang, J. Adv. Model. Earth Syst., 2019, 11, 3593–3635.
40. A. Martynov, R. Laprise, L. Sushama, K. Winger, L. Šeparović and B. Dugas, Clim. Dyn., 2013, 41, 2973–3005.
41. L. Šeparović, A. Alexandru, R. Laprise, A. Martynov, L. Sushama, K. Winger, K. Tete and M. Valin, Clim. Dyn., 2013, 41, 3167–3201.
42. M. Meinshausen, S. J. Smith, K. Calvin, J. S. Daniel, M. L. T. Kainuma, J.-F. Lamarque, K. Matsumoto, S. A. Montzka, S. C. B. Raper, K. Riahi, A. Thomson, G. J. M. Velders and D. P. Van Vuuren, Clim. Change, 2011, 109, 213–241.
43. V. K. Arora, J. F. Scinocca, G. J. Boer, J. R. Christian, K. L. Denman, G. M. Flato, V. V. Kharin, W. G. Lee and W. J. Merryfield, Geophys. Res. Lett., 2011, 38, n/a.
44. K. Von Salzen, J. F. Scinocca, N. A. McFarlane, J. Li, J. N. S. Cole, D. Plummer, D. Verseghy, M. C. Reader, X. Ma, M. Lazare and L. Solheim, Atmos.-Ocean, 2013, 51, 104–125.
45. J. P. Dunne, J. G. John, A. J. Adcroft, S. M. Griffies, R. W. Hallberg, E. Shevliakova, R. J. Stouffer, W. Cooke, K. A. Dunne, M. J. Harrison, J. P. Krasting, S. L. Malyshev, P. C. D. Milly, P. J. Phillipps, L. T. Sentman, B. L. Samuels, M. J. Spelman, M. Winton, A. T. Wittenberg and N. Zadeh, J. Clim., 2012, 25, 6646–6665.
46. A. Voldoire, E. Sanchez-Gomez, D. Salas Y Mélia, B. Decharme, C. Cassou, S. Sénési, S. Valcke, I. Beau, A. Alias, M. Chevallier, M. Déqué, J. Deshayes, H. Douville, E. Fernandez, G. Madec, E. Maisonnave, M.-P. Moine, S. Planton, D. Saint-Martin, S. Szopa, S. Tyteca, R. Alkama, S. Belamari, A. Braun, L. Coquart and F. Chauvin, Clim. Dyn., 2013, 40, 2091–2121.
47. B. Stevens, M. Giorgetta, M. Esch, T. Mauritsen, T. Crueger, S. Rast, M. Salzmann, H. Schmidt, J. Bader, K. Block, R. Brokopf, I. Fast, S. Kinne, L. Kornblueh, U. Lohmann, R. Pincus, T. Reichler and E. Roeckner, J. Adv. Model. Earth Syst., 2013, 5, 146–172.
48. S. Riette and D. Caya, in Research Activities in Atmospheric and Earth System Modelling, ed. H. Ritchie, WMO, Geneva, 2002, pp. 7.39–7.40.
49. L. Separovic, R. De Elía and R. Laprise, Clim. Dyn., 2012, 38, 1325–1343.
50. J. Zhuang, R. Dussin, D. Huard, P. Bourgault, A. Banihirwe, S. Raynaud, B. Malevich, M. Schupfner, F. Fernandes, S. Levang, C. Gauthier, A. Jüling, M. Almansi, R. Scott, S. Rasp, T. J. Smith, J. Stachelek, M. Plough, P. Manchon, R. Bell, R. Caneill and L. Xianxiang, pangeo-data/xESMF: v0.8.2, 2023, https://zenodo.org/record/4294774.
51. P. Zanis, D. Akritidis, S. Turnock, V. Naik, S. Szopa, A. K. Georgoulias, S. E. Bauer, M. Deushi, L. W. Horowitz, J. Keeble, P. Le Sager, F. M. O'Connor, N. Oshima, K. Tsigaridis and T. Van Noije, Environ. Res. Lett., 2022, 17, 024014.
52. C. C. of Ministers of the Environment, Le Conseil canadien des ministes de l’environment, Canadian Ambient Air Quality Standards | Normes canadiennes de qualité de l’air ambiant, https://www.ccme.ca/en/air-quality-reportslide-7.
53. K. P. for Ontario, National Air Quality Management System | Air Quality in Ontario 2021 Report, 2025, https://www.ontario.ca/document/air-quality-ontario-2021-report/national-air-quality-management-systemsection-7.
54. S. de l’environnement, Environmental Assessment Report 2023: Air Quality In Montréal, 2024, https://donnees.montreal.ca/dataset/24d31955-5590-47b9-a6f4-ffdac36429ba/resource/10aeb25f-7c89-41d2-9948-e40cfae9234a/download/en-bilan-rsqa-2023.pdf.
55. Agreement between the Government of Canada and the Government of the United States of America on Air Quality, 2000, https://www.epa.gov/system/files/documents/2024-08/original-aqa-text-with-ozone-annex.pdf.
56. Environment Climate Change Canada United States Environmental Protection AgencyReview and assessment of the Canada-United States Air Quality Agreement (AQA), Environment and Climate Change Canada = Environnement et changement climatique Canada, Gatineau QC, [Cat. No.: En4-651/2024E-PDF]. edn, 2024.

---