2.1 GEM-MACH overview

GEM-MACH is an online, 3-D chemical transport model, which is embedded in GEM, ECCC's operational numerical weather prediction model (Côté et al., 1998b, a; Moran et al., 2010) (available online here: https://github.com/mfvalin?tab=repositories, last access: 3 July 2018). Online models such as GEM-MACH improve air-quality chemical prediction performance by reducing interpolation errors between different model coordinate systems and removing the input/output time and disk storage required for the transfer of meteorological input files to their offline CTM counterparts (e.g., Baklanov et al., 2014). The coupling to meteorology is a one-way process in this version, whereby chemistry does not influence the meteorology. More detailed description of the gas-, aqueous-, and particle-phase process representations of GEM-MACH, and an evaluation of its performance for common pollutants such as ozone, particulate matter (PM), and ammonia appears in Moran et al. (2013), Makar et al. (2015a, b), Gong et al. (2015), and Whaley et al. (2018b). Here we will focus on the model changes made to include PAH species and processes.

GEM-MACH is used to provide ECCC's twice daily, 48 h operational public forecasts of criteria air pollutants (ozone, nitrogen oxides, and PM), as well as the Air Quality Health Index (https://ec.gc.ca/cas-aqhi/, last access: 3 July 2018.). To reduce the computational burden for forecasting, the PM size distribution is represented using a simplified sectional treatment consisting of two size bins, a fine-fraction bin for particles with Stokes diameter from 0 to 2.5 µm and a coarse-fraction bin for particles with Stokes diameter from 2.5 to 10 µm (Moran et al., 2010), with sub-binning used for those particle processes requiring a finer particulate size distribution. Here, we utilize the research version of GEM-MACH version 2, revision 2476, with this two-size-bin representation as our starting point for PAH modifications. The model grid used corresponds to a rotated latitude–longitude map projection with 2.5 km horizontal grid spacing and a hybrid vertical coordinate with 80-level vertical discretization spanning the atmosphere from the surface to 0.1 hPa.

2.2 Model modifications for benzene and PAH species

Our modifications to GEM-MACH include adding benzene and 7 gas-phase and 14 particle-phase (7 species × 2 size bins) PAHs to the species arrays. The gas–particle partitioning subroutine described in Galarneau et al. (2014) was also added, but with updated partitioning parameters (see Sect. 2.2.1). Since PAHs have very small concentrations relative to criteria air contaminants, as in Galarneau et al. (2014), we assume they do not have a significant effect on oxidant concentrations (O₃ and OH). Thus, the PAHs in GEM-MACH-PAH make use of the outcomes of the model's gas- and aqueous-phase chemistry in a diagnostic fashion for PAH oxidation. Processes in which the PAHs directly participate include advection, vertical diffusion, plume rise of major point source emissions, aerosol particle microphysics, in- and below-cloud scavenging, and dry and wet deposition of both gas and particle phases. Some of these processes and/or their controlling parameters were updated relative to Galarneau et al. (2014) and are described in the subsections below. Like AURAMS-PAH, the total (gas+particle) PAH emissions were treated as gas-phase emissions in GEM-MACH-PAH, since these quickly repartition between particles and gas phases following emission. The non-PAH and PAH emissions are described further below. The MACH part of the model code is available online (Whaley et al., 2018a; https://doi.org/10.5281/zenodo.116225), and the call sequence of the code is included in Sect. S5 of the Supplement.

2.2.1 Gas–particle partitioning

As PAHs are semi-volatile organic compounds that partition between the particulate and gaseous phases of the atmosphere, their partitioning is a major determinant of their atmospheric fate (Bidleman, 1988). Despite decades of study (Junge, 1977; Yamasaki et al., 1982; Bidleman and Foreman, 1987; Pankow, 1987; Smith and Harrison, 1996; Dachs and Eisenreich, 2000; Lohmann and Lammel, 2004; Keyte et al., 2013), the mechanisms responsible for PAH partitioning and its spatiotemporal variability are not well-understood. AURAMS-PAH included two parametrizations to calculate gas–particle partitioning: Junge–Pankow, JP (Junge, 1977; Pankow, 1987) and Dachs–Eisenreich, DE (Dachs and Eisenreich, 2000), both of which, when applied for partitioning in AURAMS-PAH, assigned too much PAH mass to the gas phase. The two schemes resulted in surprisingly similar gas–particle partitioning (Galarneau et al., 2014). We carried out post-processing and analysis on the AURAMS-PAH model output from both schemes as well as the observations of gas and particle PAHs from the Galarneau et al. (2014) study, in order to determine which scheme to proceed with in GEM-MACH-PAH, and how it could be improved.

Measured PAH partitioning typically takes the linear form of

$$\begin{array}{}\text{(1)}& \log K_{\mathrm{p},k}=m_K\log p_{L,k}^{\circ }+b_K,\end{array}$$

where K_p,k is the partitioning coefficient for each PAH species, k:

$$\begin{array}{}\text{(2)}& \log K_{\mathrm{p},k}=\log \left[{\displaystyle \frac{C_{\mathrm{p}}/C_{\mathrm{TSP}}}{C_{\mathrm{g}}}}\right],\end{array}$$

and C_p, C_TSP, and C_g are the concentrations of the particulate PAH, the total suspended particles, and the gas-phase PAH, respectively. $p_{\mathrm{L},k}^{\circ }$ is the sub-cooled liquid vapour pressure of the kth gas, and m_k and b_k are empirically derived coefficients. This linear relationship (where the logK_p of all PAH species in any given measurement sample fall on a line relative to their logp_L) is common among homologous compound groups such as PAHs, polychlorinated biphenyls (PCBs), and polychlorinated dioxins and furans. However, the JP formulation only allows for $m_K=-\mathrm{1}$ (see Sect. S2 in the Supplement for more detailed information and a derivation). Conversely, observation-based estimates show a wide variety of $\left|m_K\right|$ values that are usually less than 1 (e.g., Fig. S2a in the Supplement), and this could be the reason why the AURAMS-PAH JP model results under-predicted the particulate fraction.

Therefore, we proceeded with the Dachs–Eisenreich formulation in GEM-MACH, but with improved parameters. The Dachs–Eisenreich (DE) partitioning formulation was adapted from work examining water–sediment partitioning (Dachs and Eisenreich, 2000). The DE expression for K_p (Eq. S2.2.1) is related to the octanol–air and soot–air (K_SA,k) partitioning coefficients, the latter depending on the soot–water (K_SW,k) and air–water partitioning coefficients. The soot–water partitioning coefficients are highly uncertain. Their values in the literature span two orders of magnitude for the same compound (Dachs and Eisenreich, 2000; Bucheli and Gustafsson, 2000; Jonker and Koelmans, 2002; Xu et al., 2012). K_SW,k from Jonker and Koelmans (2002), was used in AURAMS-PAH. However, using the 2002 measurement data and their average m_K, we have determined new K_SW,k values based on ambient observations that improve the DE particulate fraction representation (see Sect. S2.3 in the Supplement for this process, and Table 1 for the original and new values). The purple boxes in Fig. 2 represent the results of the AURAMS-PAH partitioning module, making use of our new K_SW,k values instead of the originals.

Figure 2All-site ensemble of modelled/measured ratios of logK_p (a) and particulate fraction (b) for each PAH. Shown are the “adjusted model” (purple, from Eq. S2.2.4) and the original AURAMS-PAH model (green). Box and whiskers: thick line is the median, boxes extend to the 25th and 75th percentiles, and whiskers extend to the minimum and maximum.

Download

Dachs and Eisenreich (2000)Jonker and Koelmans (2002)Xu et al. (2012)Bucheli and Gustafsson (2000)

Table 1Original (from Galarneau et al., 2014, based on Jonker and Koelmans, 2002) and adjusted (based on AURAMS-PAH model–measurement analysis for North America) K_SW values, along with other K_SW cited in the literature. All values are × 10⁵, and unitless.

Download Print Version | Download XLSX

While the adjusted K_SW values in Table 1 are significantly different from those in the original model (based on Jonker and Koelmans, 2002 as adjusted by the relative contributions of PM mass to the domain total in the inventory of Galarneau et al., 2007), particularly for lower molecular-weight species, they fall within the range of values found in the literature (e.g., Dachs and Eisenreich, 2000; Bucheli and Gustafsson, 2000; Jonker and Koelmans, 2002; Xu et al., 2012). The large correction and large range of published K_SW,k values (which lack a temperature dependence) lends further support for the need to measure K_SA,k directly on soots of atmospheric relevance so that they need not be estimated based on the highly uncertain/variable K_SW,k.

2.2.4 Dry and wet removal of PAHs and benzene

Gas-phase dry deposition follows a multiple resistance approach and single-layer “big leaf” approach (Wesely, 1989; Zhang et al., 2002; Makar et al., 2018) with a temperature dependency for Henry's law constants and for water solubility (Sander, 1999; Ma et al., 2010). Dry deposition of benzene and PAHs is also output by the model; however, measurements of the deposition flux of these species were unavailable during the study period.

PAH particle-phase dry deposition is treated following Zhang et al. (2001), resulting in size-dependent particle deposition velocities.

Gas-phase benzene and PAHs undergo cloud and rain scavenging via Henry's law. Henry's law partition coefficients (K_AW) for the seven PAHs are proportional to inverse temperature. The mass of benzene and PAHs in the gas-phase (as opposed to the aqueous phase in cloud droplets and raindrops) is derived from

$$\begin{array}{}\text{(5)}& K_{\mathrm{AW},k}={\displaystyle \frac{m_{\mathrm{gas}}/V_{\mathrm{air}}}{m_{\mathrm{aq}}/V_{\mathrm{h}\mathrm{2}\mathrm{o}}}}=m_{\mathrm{AW}}/T+b_{\mathrm{AW}},\end{array}$$

and solving for m_gas:

$$\begin{array}{}\text{(6)}& m_{\mathrm{gas}}={\displaystyle \frac{\frac{V_{\mathrm{air}}}{V_{\mathrm{h}\mathrm{2}\mathrm{o}}}{K}_{\mathrm{AW},k}{m}_i^{\mathrm{gas}}}{(\mathrm{1}+\frac{V_{\mathrm{air}}}{V_{\mathrm{h}\mathrm{2}\mathrm{o}}}{K}_{\mathrm{AW},k})}},\end{array}$$

where, $m_i^{\mathrm{gas}}$ is the initial mass of the gas-phase PAH before the Henry's law partitioning, and the remaining PAH mass is scavenged to the liquid rain or cloud phases (m_aq). Once the gaseous PAHs are scavenged to cloud-water (which is recalculated at each time step), they are subject to rain-out (cloud-to-rain conversion) process. At the end of each chemistry step, the fraction of the dissolved tracers not removed by the rain-out process will be returned to gas phase. For the fraction that do go into rain water, there is a parametrization in the rain scavenging code for re-evaporation of the rain before it reaches the ground. That fraction of PAHs would also return to the atmosphere. The remaining fraction that is not re-evaporated is counted towards the wet deposition in the model, which is irreversible (i.e., no re-emission of PAHs from the ground after they are deposited).

Note that Okochi et al. (2004) reported that assuming Henry's law equilibrium for benzene under-predicts the extent of wet-deposition. In the absence of a suitable alternative parameterization, we used Henry's law partitioning and obtained a conservative estimate of wet deposition for benzene. Note that benzene wet deposition is not evaluated in this paper as no measurements are available.

Figure 3NAPS, NATTS, IADN, and Hamilton measurement sites in the model domain.

Download

Where temperatures are < 0 ^∘C below-cloud, or $<-$15 ^∘C in-cloud, scavenging of gas-phase benzene and PAHs by snow and cloud-ice is done via surface adsorption following the formulation used in Franz and Eisenreich (1998), which was also used by Wania et al. (1999), Lei and Wania (2004) and Friedman and Selin (2012):

$$\begin{array}{}\text{(7)}& W_{\mathrm{g}}=K_{\mathrm{ia},k}\left(\mathrm{SA}\right)\mathit{\rho }.\end{array}$$

In Eq. (7), W_g is the gas scavenging ratio (equal to the concentration of PAH in snow over the concentration of PAH in air – both in moles m⁻³), K_ia,k is the interfacial adsorption coefficient (equal to the mass adsorbed per surface area of snow to the atmospheric vapour phase concentration – both in ng m⁻³), SA is the specific surface area of the snow crystal, for which we use a constant 1 m² g⁻¹ based on literature values for fresh snow precipitation, which are highly variable, and for which no clear relationship with temperature or wind speed has been found (Hoff et al., 1998; Hanot and Dominé, 1999; Domine et al., 2007; Hachikubo et al., 2014), ρ is the density of ice (0.917 g cm⁻³), and K_ia,k is calculated from the following (Franz and Eisenreich, 1998):

$$\begin{array}{}\text{(8)}& \log \left(K_{\mathrm{ia},k}\right)=-\mathrm{1.2}\log p_{\mathrm{L},k}^{\circ }-\mathrm{5.82}.\end{array}$$

Eq. (7) is used to determine the fraction of PAH mass in the gas and snow/cloud-ice phases.

Particle-phase PAHs are treated as passive tracers that undergo wet removal along with the modelled aerosol particles (Gong et al., 2006; Wang et al., 2010). The cloud and precipitation processes above are applied sequentially in the model using operator splitting, and the amount of PAH deposited from wet deposition is output by the model. Table S4 in the Supplement provides all of the constants used in the model.

2.3 Model setup for two three-month simulations

GEM-MACH-PAH (rev2488) was run from 8 May to 13 August 2009 and from 18 October 2009 to 5 January 2010, where the first week from each period is treated as a spin-up period (for chemical concentrations to stabilize and for the initial condition effects to be negligible: e.g., Samaali et al., 2009), and were not used in our evaluation. The time periods were chosen to coincide with as many PAH concentration and deposition measurements as possible, while limiting simulation duration to reduce computational expenses.

The chemical initial and boundary conditions for the outer nest North American domain were taken from a 1-year MOZART simulation for all pollutants (Emmons et al., 2010; Pendlebury et al., 2017), except for benzene and PAHs. Initial and boundary conditions for PAHs and benzene were set to zero for the North American domain as its boundaries are generally away from PAH and benzene sources (e.g., over the ocean), and are also very distant from the Pan Am domain. The simulations for the nested Pan Am region were run using the chemical initial and boundary conditions from the 10 km North American model run.

Figure 4(a) Map of 2-week summertime average fluoranthene concentrations in Hamilton, Ontario: (left) from measurements and GEM-MACH-PAH model, and (right) their differences and ratios. Note the grey dots are measurement sites that had missing data during this time period. (b) Spatial variability in the Hamilton data in (left) summer, and (right) winter.

Download

The model simulation was carried out in sequence of 27 h staggered simulations starting at 00:00 UTC, in order to re-initialize meteorology with the analysis at 10 km resolution. The first 3 h in the 2.5 km domain were discarded as spin-up to reduce the dependency on the 10 km resolution meteorological initial conditions. Each 27 h simulation used the chemical concentrations from the end of the previous simulation as initial conditions for the next 27 h, and this sequence continued until each 3-month period was complete.

GEM-MACH-PAH was run in the two-size-bin mode to represent the PM size distribution, which means that particles fall in either fine mode (PM_2.5 – diameter 2.5 µm or less) or coarse mode (PM₁₀–PM_2.5 – diameter 2.5–10 µm).

3.1 Hamilton measurement campaign

Ambient measurements of PM_2.5 and 16 PAH species were collected from a dense network of measurement sites in Hamilton, Ontario during June/July 2009 and December 2009. These measurements are described in Anastasopoulos et al. (2012), where they found a high level of intra-urban variability for the PAHs; 3–4 times more variable than PM_2.5 concentrations.

There were 43 measurement sites operating during the summer period (see Fig. 4a), and 46 sites during the winter period. All measurements are from 2-week integrated time frames (24 June to 8 July and 2 to 16 December) taken with URG personal pesticide samplers, which collected gas and particle-phase PAHs less than 2.5 µm in diameter in 40 m³ of sample air. The PM_2.5 measurements were made at the same sites using a three-stage Harvard cascade impactor. The particle-phase PAHs (up to 2.5 µm in diameter) were collected on a Teflon filter, gas-phase PAHs were collected in polyurethane foam (PUF), and total (gas+particle) PAHs concentrations were reported in ng m⁻³ as determined by gas chromatography/mass selective detection. PM sample filter masses were determined by gravimetric analysis (Anastasopoulos et al., 2012).

O₃ measurements that were used to correct the Hamilton BaP measurements came from three monitoring sites in the Hamilton region (“Downtown”, “Mountain”, and “West”) from the Ontario Ministry of Environment and Climate Change (MOECC) website for historical air quality data (http://www.airqualityontario.com/history/, last access: 3 July 2018).

3.2 National Air Pollution Surveillance program

NAPS is a Canadian program to provide accurate and long-term air quality data of a uniform standard across the country. NAPS is managed under a cooperative agreement between ECCC and the provinces, territories, and some municipal governments. There are currently 286 NAPS measurement sites in 203 communities located in every province and territory (https://open.canada.ca/data/en/dataset/1b36a356-defd-4813-acea-47bc3abd859b, last access: 3 July 2018).

Under this program, PAH samples were collected over 24 h, beginning and ending at midnight (local time; LT), typically every 6 days, with a sample volume range of 600–800 m³ (Environment Canada, 2013). Benzene samples were collected in 6 L stainless steel canisters over 24 h, starting at midnight, every 3 days (Galarneau et al., 2016).

Within the Pan Am domain, total PAHs (gas+particle-phases combined) and benzene were measured at eight NAPS sites (listed in Table S1 in the Supplement; Fig. 3), and their 2009 data were downloaded from the following url: http://maps-cartes.ec.gc.ca/rnspa-naps/data.aspx (last access: 3 July 2018).

For BaP measurement corrections, the NAPS network also measures hourly O₃ at four of these eight PAH/benzene sites (Windsor, Hamilton, Simcoe, and Egbert). Two of the missing sites (Toronto and Etobicoke) had nearby O₃ measurements from MOECC, but the last two (Burnt Island, and Point Petre), which are rural sites, had no O₃ measurements nearby. Therefore, BaP could only be corrected at six of the eight NAPS sites in the Pan Am domain.

3.3 National Air Toxics Trends Stations network

NATTS is a US program to monitor toxic air pollutants in accordance with the US Government Performance Results Act, which requires the US EPA to reduce the risk of cancer and other serious health effects associated with hazardous air pollutants (HAPS) by achieving a 75 % reduction in toxic chemical emissions in air, based on 1993 levels (US EPA, 2009). Regulated under the Clean Air Act are 188 HAPS species including benzene and the 7 PAHs in this study.

Every 6 days, 24 h ambient air samples are collected starting at midnight (00:00 LT). Analysis of the samples is done by high-resolution gas chromatography/mass spectrometry (GCMS) selective ion monitoring (SIM) mode to get total (gas+particle) PAH concentrations, and benzene concentrations (Eastern Research Group, Inc., 2009).

There are 115 NATTS sites within the model domain (Table S1 in the Supplement, Fig. 3), but only 21 sites measured PAHs, while 113 sites measured benzene. Those data were downloaded from the following url: https://www3.epa.gov/ttnamti1/toxdat.html#data (last access: 3 July 2018). Measurement methods in NATTS are very similar to those of NAPS.

Since O₃ was not measured at the NATTS sites, NATTS BaP was corrected with the nearest O₃ monitor data found at the US EPA and CASTNET websites: https://www.epa.gov/outdoor-air-quality-data/download-daily-data (last access: 3 July 2018) and https://java.epa.gov/castnet/reportPage.do (last access: 3 July 2018), respectively.

3.4 Integrated Atmospheric Deposition Network

IADN was mandated by the 1987 Great Lakes Water Quality Agreement, and was initiated in 1990 to measure atmospheric concentrations of persistent toxic pollutants in the Great Lakes basin. There are nine IADN sites total within our Pan Am model domain, and they are listed in Table S1 of the Supplement (see also Fig. 3). Six of the nine IADN sites report gas- and particle-phase PAH atmospheric concentrations separately (labelled “PAHs” in Table S1), and a different set of six sites report wet deposition of PAHs (labelled “PAH wet dep” in Table S1) using sampled precipitation concentrations (Blanchard et al., 2005). Thus, these data can be used to evaluate the model's gas–particle partitioning and deposition, respectively. Benzene was not measured by IADN, nor was C_TSP or PM₁₀ in 2009, with the unfortunate result that K_p can not be directly calculated from the IADN measurements. O₃ was also not measured by IADN, necessitating the use of the nearest O₃ monitors in order to carry out the BaP oxidation correction. This latter step was only possible for observation stations at Cleveland and Chicago. The other four IADN air sites were rural/background locations, and did not have any O₃ measurements nearby.

PAHs were collected by a high-volume sampler for periods of 24 h beginning at 08:00 EST, every 12 days. At Canadian IADN sites, glass fiber filters and PUF sorbent collected the particulate and gaseous fractions, whereas the US stations collected PAHs with quartz fiber filters and XAD resin (Blanchard et al., 2005). Sample volume for the US method is about 800 m³, but is 400 m³ for the Canadian method to minimize the breakthrough of volatile species during warm summer months (Blanchard et al., 2005).

Wet deposition of PAHs are measured with MIC-B precipitation collectors. The US stations used XAD-2 resin column cartridges for accumulating the organics on a 28-day cumulative basis, while the Canadian stations use a dichloromethane solvent extraction system, also on a 28-day cumulative basis. Both countries collect samples on a monthly basis. Note that one of the six wet deposition sites, St Clair, Ontario (STC), only had valid measurements during February 2009, which was not a time period simulated here. Therefore, only five IADN sites appear in our wet deposition analysis in Sect. 4.4 below.

Note that Point Petre and Burnt Island are NAPS stations co-located with IADN. IADN data were downloaded from the following url: http://donnees.ec.gc.ca/data/air/monitor/monitoring-of-combined-atmospheric-gases-and-particles/anthropogenic-organic-pollutants-in-air/ (last access: 3 July 2018).

4.1 PAH concentrations in the Hamilton region

Figure 5(a) GEM-MACH-PAH model vs. measurement scatter plot of 2-week summertime fluoranthene concentrations at 40+ sites in Hamilton. (b) Frequency distributions of GEM-MACH-PAH model/measurement ratios of PAH concentrations for the Hamilton model–measurement pairs for all sites from both summer and winter. Darker-coloured boxes are results from O₃-corrected BaP measurements.

Download

GEM-MACH-PAH output for gas+fine-PM PAH were compared to measurements of same from the 2009 Hamilton campaign (Anastasopoulos et al., 2012). Figure 4a shows a map of measured and modelled fluoranthene concentrations (14-day average) in the summer time period, as well as their differences and ratios. Here we see that GEM-MACH-PAH has captured intra-city variability, and that the differences between observations and simulated values are, at a maximum, 2.8 times too high. The model is biased low in the upwind/background areas of the city, and a high in the eastern areas of the city (Fig. 4a); this pattern is seen across all seven PAHs. The spatial pattern in the PAH bias is less apparent when PAH ∕ PM_2.5 ratios are plotted (in ng µg⁻¹ – shown in Fig. S6.1 in the Supplement) – removing the dependency on modelling PM correctly (since fractions of the PAH are particulate). Therefore, the spatial pattern in the PAH bias is mainly due to the pattern in the model PM bias, which is shown in Fig. S6.2 in the Supplement.

When the spatial variability is represented by the standard deviation over the mean (σ/mean), the model achieves very similar spatial variability to the measurements (Fig. 4b). The scatter plot of model vs measurements for summertime fluoranthene concentrations (Fig. 5a, FLRT selected as a good example) has a correlation coefficient (R²) of 0.57, and the slope of the best-fit line is 0.96. The other PAH species had similar results, where, except for FLRT, the slopes and R² values were better in the winter than in the summer.

The model bias (given as a model/measurement ratio) for all PAHs is shown as box and whiskers in Fig. 5b. Here we see that wintertime biases are smaller than those in the summertime for all PAHs except for ANTH. The four lightest PAHs (left side of Fig. 5b) have model/measurement ratios near one (except summertime PHEN), but the three heaviest PAHs, are biased high (except for wintertime CHRY). We will see this same pattern for the model bias (small for lighter PAHs, high for BaA and BaP) in the next sections as well.

BaA and BaP are the most reactive of the heavier species; thus, the lack of O₃ correction to the BaA measurements may be partially responsible for the model–measurement differences. However, mean O₃ for the three measurement stations in the Hamilton region was only 20–26 ppbv day⁻¹, which meant that the O₃-corrected BaP was approximately 20 % greater than the reported BaP concentrations. The median BaP bias was brought down to 6.2 from 7.6 in the summer, and 5.5 from 6.3 in the winter – these are shown as the purple boxes in Fig. 5b. Additional reactions with BaA and BaP (e.g., with NO₃) are noted in the literature (Keyte et al., 2013), but were not included in GEM-MACH-PAH at this time, given larger uncertainties in those reactions. However, our model biases appear to indicate that those missing reactions may need to be considered for further model improvement. Indeed, a back-of-the-envelope calculation using a reaction rate reported in Liu et al. (2012) implied that the lifetime of particulate BaA with respect to a heterogeneous reaction with NO₃ (using mean nighttime concentrations of NO₃ from the model) was about 1 h. More research should be undertaken to determine the range of uncertainty for this reaction and whether it should go into PAH models. Additionally, Mu et al. (2018) suggest that the heterogeneous BaP–O₃ reaction should be temperature-, humidity-, and organic aerosol phase state-dependent (none of which were taken into account in the Kwamena scheme used in our work). However, it has been shown that the Kwamena scheme and the Mu scheme produce similar results in mid-latitudes (where our study is located) (Mu et al., 2018). Spring/summertime BaP would be minimally affected, as outdoor temperatures at that time of year resemble the room temperature laboratory conditions that the Kwamena scheme was based on. Additionally, our positive model bias would likely increase in the fall–wintertime, when low temperatures and humidity would increase BaP lifetime in the Mu scheme.

In order to remove the impact of the model's PM predictions on the PAH comparison, we also plotted the PAH ∕ PM_2.5 model-over-measurement ratios (shown in Fig. S6.3). There we see all of the ratios reduced – which improves results for the heavier PAHs, but increases the low bias for the lighter PAHs. The reason the bias decreases for all seven PAHs is that the model PM_2.5 is overestimated by a factor of 2 in the summertime, and a factor of 1.4 in the wintertime (average across all sites in the Hamilton region).

4.2 PAH and benzene concentrations from the NAPS, NATTS, and IADN networks

Modelled 24 h-average total (gas+particle) PAH and benzene concentrations can be evaluated with measurements (taken every sixth day) from the NATTS, NAPS, and IADN surface measurement networks, which sample much of the model domain well (Fig. 3). As with the Hamilton evaluation, the model has very good agreement for seasonal averages at the monitoring network sites for benzene, phenanthrene, anthracene, fluoranthene, and pyrene, which all have model/measurement ratios (red and blue boxes) close to one, and their concentrations (green and orange boxes) overlapping in Fig. 6a. The 24 h average model (daily) and measurements (every sixth day) have been averaged over each of the 3-month time periods. BaA and CHRY are within a factor of 5 of the measurements in the summertime, but worse in the wintertime. BaP is overestimated by the model in both summer and winter by about a factor of 10, although the measurement-corrected BaP has a slightly reduced bias. We have not shown the O₃-corrected BaP measurements in the plots because the changes are small, similar to the Hamilton plot (Fig. 5a).

Figure 6(a) Frequency distributions of GEM-MACH-PAH (green) and measured (orange) benzene (gas) and PAH (gas+particle) seasonal average concentrations at all IADN, NAPS, and NATTS sites. Modelled/measured concentration ratios also shown for summer (red) and winter (blue), with grey lines indicating agreement within an order of magnitude. (b) Modelled/measured concentrations for each daily model–measurement pair, separated by site (FLRT given as example); (c) same as (b) but separated by month.

Download

When the model biases are examined more closely we find a few patterns to determine the cause(s). The following list summarizes some observations from our evaluation of each model–measurement pair (24 h averages, not seasonal averages):

• By site – overestimates: all PAHs are significantly overestimated at the Kennedy Township, Pennsylvania (northwest of Pittsburgh) NATTS site (Fig. 6b). There appears to be a major emissions point source near that station that is emitting too much PAH in our model compared to reality. There are in fact hundreds of point sources in the emissions inventory that are within 20 km of Kennedy Township, but one in particular emits a relatively large amount of VOCs, and is associated with the “Secondary Metal Production; Aluminum; Raw Material Charging” source category, which has very large PAH-to-TOG EFs in Galarneau et al. (2007) because aluminum smelter emissions are largely particulate; therefore, expressing EFs as a large fraction of TOG was somewhat artificial. However, our results indicate that the PAH EFs for that PAH speciation profile (1036b) should be reduced substantially compared to Galarneau et al. (2007). In order to ensure that this facility did not begin operation after 2009 – which is a risk when using a 2011 inventory to model the year 2009 and would result in a large overestimation – we have further verified that the facility existed and was emitting similar VOC amounts in the NEI2008 inventory as well.

  PHEN (Fig. 7c) and ANTH (not shown) are also greatly overestimated in New York City, although none of the other PAHs are biased particularly high there. However, we note that the measurements for New York City appear erroneously low, as the reported PHEN concentrations there are around the same magnitude as those in Underhill, VT (Fig. 7c), which is a background site near a national park.

  Most PAHs are also overestimated at the Gary, Indiana site (Fig. 6b), which may also have a nearby major point emissions source that is too high compared to reality. Furthermore. the heavier PAHs (BaA, CHRY, and BaP) are overestimated at the Toronto Gage Institute NAPS site, but are only slightly higher there than the average model/measurement ratio for those species (not shown).

• By site – underestimates: all PAHs are markedly underestimated at the Liberty, Pennsylvania site (e.g., Figs. 6b, and 7c), implying that there may be industry emissions of PAHs here that are missing, misallocated, or misplaced in the NEI2011 inventory, or an improper PAH speciation profile applied. Similarly, PAHs are underestimated in Buffalo, New York, and at Franklin Furnace, Ohio. As these are not large cities, there may be industrial emissions that are not reported in the NEI2011 emissions inventory (or are reported at too low levels) – perhaps because those facilities shut down in 2010 (or installed emission control technology), which would mean the problem simply lies in using a 2011 inventory to model 2009.

  However, when we further investigated stacks near Buffalo, NY, we found that the facility with the largest CO and VOC emissions had zero PAH emissions. This facility is associated with the generic process of “Primary metal production; By-product Coke Manufacturing”, which did not have an associated PAH-to-TOG profile in Galarneau et al. (2007), because the source category codes that follow it (such as flushing liquor circulation tank, excess-ammonia liquor tank, tar dehydrator, tar interceding sump, tar storage, and so on) are not expected to emit PAHs to air. However, our results imply that the PAH speciation profile for “By Product Coke Oven Stack Gas” (0011b) would have been more appropriate for this facility and its use might eliminate the model bias near Buffalo in future studies.

• By month: all PAH species have lower mod/meas ratios in the summer than in the winter (shown in Fig. 6a by season for all PAHs and in Fig. 6c for FLRT by month) – implying that modelled hydroxyl radical (OH) and/or PM biases (which have strong seasonal cycles) are impacting modelled PAHs. For example, if model OH is too high in the summer, or too low in the winter, this would cause the U-shaped pattern that we see when plotting model/measurement ratio vs. month (Fig. 6c) and it would be particularly pronounced for the lighter, gas-phase PAHs, which it is. Another possibility is seasonal bias in the representation of atmospheric vertical stability: if the modelled stability is too low in the summer and too high in the winter, then winter emissions will tend to be trapped in inversions more than observed, and summer emissions will be diluted by excessive vertical mixing. However, evidence in Makar et al. (2010) and Stroud et al. (2012) suggest that model stability is too high (not too low) for the summer time period in those studies.

• By season: for the four lightest PAHs, the model/measurement ratios are < 1 in the summer, and > 1 in the winter (Fig. 6a). As mentioned above, this is likely due to modelled OH being too high in the summer and too low in the winter. BaA and CHRY follow a similar seasonal difference but do not straddle the ratio = 1 mark.

  BaP, in comparison, has a model/measurement ratio that is slightly higher in the summer than in the winter (Fig. 6a). For BaP, the OH bias could be offset by an opposite O₃ bias in the model. Indeed, it has been shown (Makar et al., 2010; Stroud et al., 2012) that the processes in GEM-MACH cause urban surface O₃ to be too low in the summertime, due to insufficient vertical mixing and excessive titration from NO_x, and surface PM tends to be too high in the wintertime due to an overestimation of wintertime atmospheric stability (e.g., lack of an urban heat island parameterization in the driving meteorology). These factors, together with the BaP measurement bias due to on-filter reaction with O₃, may explain the high model BaP bias.

Figure 7Phenanthrene time series for the summer 2009 period for the IADN (binational), NAPS (Canada), and NATTS (US) networks. Orange represents the measurements, dark green represents the base GEM-MACH-PAH model, light green represents GEM-MACH-PAH with 0.5× mobile emissions, and cyan represents GEM-MACH-PAH with 2.0× mobile emissions.

Download

Thus, the generally high bias of modelled BaP may to be due to additive OH, O₃, and PM model biases (plus the missing O₃ denuder technology in the measurements), impacting BaP more than the other species because BaP has the highest O₃ reactivity and the highest particulate fraction of the seven PAHs examined here.

When the five measurement sites mentioned in “By Site”, above (Kennedy Township, PA; Gary, IN; Liberty PA; Buffalo, NY; and Franklin Furnace, OH) are removed from the NATTS analysis (because errors in their nearby emissions were identified), model–measurement correlation (R) and slopes improve. For example, the model vs. measurement best-fit line slope for PHEN doubles from 0.3 to 0.6 when those sites are removed, and its R increases from 0.16 to 0.32. The slopes and R values of the four heaviest PAHs all move from negative to positive. PYR has the largest improvement, going from a slope of −0.049 and R of −0.028 to a slope of 0.26 and R of 0.35. To the extent that the model prediction errors at the other sites may reflect emissions inaccuracy, having an accurate major point emission inventory for the time period modelled, along with proper PAH speciation profiles are extremely important requirements for modelling PAHs well at high resolution. The cases with large discrepancies mentioned above highlight the need to be as specific as possible when assigning source category codes to facility processes (which is difficult given that there are tens of thousands of point sources in the inventories).

In addition to emissions errors (for which fine temporal information is lacking), transport errors can cause poorer agreement between model and measurements at high spatial- and time-scales. For example, PAH concentrations have sharp gradients and large dynamic range. Thus, if plume transport is off by a few kilometres, the result will be a very large difference between modelled and measured concentrations downwind of point sources.

That said, when doing a paired t test on all data in the model domain to examine whether the summertime and wintertime modelled averages are the same as the measured averages, we found that the model was indistinguishable from the measurements for all PAH species (t<1 and p>0.05), except for wintertime BaP (which has t>1, and p<0.05, even with the O₃-corrected measurements). At finer time scales (e.g., daily model–measurement pairs) only modelled ANTH was statistically indistinguishable from measurements. Therefore, GEM-MACH-PAH can accurately model benzene and PAHs seasonally, but not daily.

4.2.1 Sensitivity of model to mobile emission factors

As discussed in the previous section, ensuring the accuracy of major point source emissions is important for model–measurement agreement near industrial locations. However, those major point source emissions tend to be located far from large population centres where human exposure is concentrated. In our inventory, mobile emissions make up 44, 45, 19, 32, 14, 21, and 30 % of total PAH emissions, for PHEN, ANTH, FLRT, PYR, BaA, CHRY, and BaP, respectively, in our continental model domain, and studies have shown that the bulk of emissions within population centres is likely to originate from on-road mobile sector emissions (Dunbar et al., 2001; Pachón et al., 2013; Kuoppamäki et al., 2014; Miao et al., 2015). Thus, in order to accurately model ambient PAHs in urban centres, the uncertainty in emission factors from on-road vehicles may play a more significant role than major point sources.

Therefore, we carried out 2-week sensitivity simulations (9–24 May and 18 October–2 November 2009) wherein the mobile emissions of PAHs were scaled by factors of 0.5 (halved) and 2.0 (doubled). This is approximately equivalent to the 25th and 75th percentiles in the range of emission factors found in the recent literature.

In Fig. 7, we show the surface PHEN time series from the measurements, base model run, and 0.5 and 2.0 scaled model runs at the IADN, NAPS, and NATTS sites. It is clear that – while a relatively small PAH source overall – changes to mobile emissions cause a large change in ambient PAH concentrations at certain urban locations, such as Philadephia (PA), New York (NY), and Burlington, Etobicoke, and Windsor (ON).

On average, there is about a 20–30 % increase in PAH concentrations when mobile emissions are doubled, and a 5–10 % decrease in PAH concentrations when mobile emissions are halved (Fig. 8, PHEN and BaP shown as examples) – with a larger sensitivity in the summer than the winter, and slightly larger sensitivity at NATTS (US) sites than at NAPS (Canada) sites. The predicted ambient concentrations generally follow the increase or decrease in on-road mobile emissions monotonically.

Figure 8Average percent change in surface PHEN and BaP concentrations by season when PAH on-road mobile emissions are scaled up or down by factors of 2.0 and 0.5, respectively.

Download

4.3 Gas–particle partitioning of PAHs

The IADN network also allows us to assess model predictions of gas–particle partitioning of PAHs at six sites (24 h averages, every 6 days). Figure 9 shows a time series of pyrene particulate fraction (ϕ_k). Both model and measurements show higher ϕ_k in the wintertime, when there are higher PM concentrations for PAH adsorption, and lower temperatures. Generally, the model seems to underestimate ϕ_k at background sites (e.g., Burnt Island), and overestimate ϕ_k at urban sites (e.g, Chicago), and this is true for all PAH species. Thus, in Fig. 10, which shows the results for all PAHs at all sites, the model (green) has a larger range of ϕ_k than the measurements (orange). This is caused by the model overestimating and underestimating PM concentrations at urban and rural sites, respectively. For example, wind-blown dust is not included in the model; however, it is known to be a potentially significant contributor to total PM in rural areas. Also, due to an underestimate of vertical mixing in the model, PM tends to be biased high in urban areas, near emissions, due to a lack of a parameterization for urban heat islands (Stroud et al., 2012).

Figure 9Time series of pyrene particulate fraction at six IADN network sites (BNT represents Burnt Island, CLV represents Cleveland, IIT represents Chicago, PPT represents Point Petre, SBD represents Sleeping Bear Dunes, and STP represents Sturgeon Point). GEM-MACH-PAH values are denoted by (green dots) and IADN measurements by orange dots.

Download

Figure 10GEM-MACH-PAH modelled (green) and measured (orange) particulate fraction (ϕ) of all PAHs at all IADN sites, and their model/measurement ratios (blue). The blue line indicates the one-to-one line, and the grey lines are for ratios of 10 and 0.1.

Download

Comparing Fig. 10 (blue boxes) to Fig. 2b (green boxes), we see significant improvement over the original AURAMS-PAH partitioning due to the improved K_SW parameters described in Sect. 2.2.1.

Generally, the gas–particle partitioning scheme in the model results in model/measurement ratios well within an order of magnitude, given by the gray lines in Fig. 10. ϕ_k for BaA and CHRY are still underestimated, but this may be related to modelled PM errors as noted earlier. We note that the addition of C_TSP or even PM₁₀ measurements at IADN sites (which existed in 2002, but not in 2009) would allow for the calculation of measured partitioning coefficients (logK_p, Eq. 2), which could be used to validate the modelled logK_p in future work. Since K_p takes total suspended particle into account, it removes the dependency on modelled PM; thus, it would increase confidence that the modelled partitioning is working properly, despite model errors in PM.

The fact that the GEM-MACH-PAH model partitioning of BaA and CHRY (and BaP to a lesser extent) puts too much concentration in the gas phase, may help explain why these species in particular are overestimated in the model. While in the gas phase, these species are less likely to be removed from the atmosphere, so their concentrations would erroneously build up in the model.

4.4 Wet deposition of PAHs

Table 2Mean and median GEM-MACH-PAH model/measurement ratios for PAH wet deposition.

Download Print Version | Download XLSX

When compared to the IADN 1-month wet deposition measurements, the model generally overestimates wet deposition for all PAHs, as is shown in Table 2 and Fig. 11. In Fig. 11, the blue lines show the ideal 1:1 model/measurement ratio, and most of the data lie well above these lines. By site (Fig. 11a), the modelled wet deposition was slightly better at urban locations (Toronto and Cleveland) than suburban and background sites (Burlington, Sturgeon Point, and Point Petre). By month (Fig. 11b), the wet deposition from the model is best represented in June and July, whereas, wet deposition is greatly overestimated in the winter, implying that the current snow adsorption parametrization may be too effective at removing PAHs in the model. This may be due to our simplification of a constant snow surface area, which may be set too high, or due to inaccuracies in the modelled or measured precipitation.

Figure 11GEM-MACH-PAH model/measurement wet deposition ratios for all PAHs (a) for five sites (all months) and (b) for four months (all IADN sites).

Download

Figure 12One-month (June 2009) wet deposition of pyrene from the (a) measurements and GEM-MACH-PAH model and (b) their differences and ratios.

Download

The IADN measurements, which report the concentration of PAHs in the collected precipitation (in pg L⁻¹), were converted to pgPAH m⁻² in order to compare to the wet deposition output of the model. However, this conversion assumes that the volume of precipitation reported by IADN was the total precipitation in the container's cross-sectional area, and in fact, it is not. The IADN wet deposition collectors are actually known to not sample all of the precipitation because the samplers are not in the correct configuration to get an accurate precipitation measurement (Dryfhout-Clark, personal communication, 2017). When we compared the actual precipitation amounts (from separate meteorological rain gauge data) to IADN precipitation volumes at the Point Petre location in January 2009, we found that only 68 % of the total precipitation was captured by the wet deposition sampler. Therefore, if that correction factor were applied to all IADN wet deposition measurements, they would increase by a factor of approximately 1.5, which would improve our comparison, but not eliminate the total bias. If IADN sites added separate, accurate rain gauges, then we could apply a “precipitation correction” to the IADN wet deposition measurements in a thorough, consistent way in future work.

Aside from the measurement bias, the modelled PAH wet deposition bias will also be dependent on the model's overall ability to predict accurate precipitation. We compared the modelled daily accumulated precipitation to the precipitation measured with the accurate gauges at Burnt Island and Point Petre, and found that, while the model's median precipitation bias was only about 0.2 mm, there was a large standard deviation, and there were some incidences where the model greatly over-predicted high precipitation events. Those incidences would result in greater modelled wet deposition of PAHs than was measured, and because we sum over a month, there is a significant likelihood of an over-prediction occurring in that long time frame. Indeed, the median ratios in Table 2, which are less sensitive to high outliers than the mean is, are substantially lower than the mean for most species.

Therefore, the model bias in wet deposition would appear to be caused by three additive factors: (1) measurements themselves having a negative bias relative to reality, due to insufficient capture of the net fluxes of precipitation; (2) modelled precipitation being biased high; and (3) a positive model bias in atmospheric PAH concentrations (which was highest for BaA and BaP in particular). Other PAH models have reported similar overestimates of PAH deposition (Matthias et al., 2009; Friedman and Selin, 2012).

The reverse reasoning can be applied, whereby we can see if high atmospheric concentrations of BaA and BaP were caused by too little wet deposition. In this case, since both of these species have correspondingly high wet deposition in the model (Fig. 11), it would appear that underestimation of wet deposition is probably not one of the causes.

Figure 12 shows results from a sample month (June 2009) for pyrene (PYR). The spatial distribution of wet deposition was not captured, with the model predicting lower PYR deposition in Toronto and Sturgeon Point than the measurements, higher at Point Petre, and about equal at Burlington. This spatial pattern is not the same for all PAHs, and even differs by month for the same PAH (e.g., PYR deposition in the next month, July, is low at Point Petre, high at Toronto, and highest in Burlington). The lack of spatial or temporal pattern in the sign and/or magnitude of wet deposition biases indicates that there is no major error with the PAH scavenging scheme itself. Given that wet deposition of PAHs relies on the correct simulation of many model factors (meteorology, scavenging parameters, atmospheric concentrations, etc.), our work suggests that more process studies aimed at quantifying wet deposition are needed. In fact, other PAH models also overestimate PAH deposition (Matthias et al., 2009; Friedman and Selin, 2012).