Evaluation of gas-particle partitioning in a regional air quality model for organic pollutants

Persistent organic pollutants (POPs) are of considerable concern due to their well-recognized toxicity and their potential to bioaccumulate and engage in long-range transport. These compounds are semi-volatile and, therefore, create a partition between vapour and condensed phases in the atmosphere, while both phases can undergo chemical reactions. This work describes the extension of the Community Multiscale Air Quality (CMAQ) modelling system to POPs with a focus on establishing an adaptable framework that accounts for gaseous chemistry, heterogeneous reactions, and gas-particle partitioning (GPP). The effect of GPP is assessed by implementing a set of independent parameterizations within the CMAQ aerosol module, including the Junge– Pankow (JP) adsorption model, the Harner–Bidleman (HB) organic matter (OM) absorption model, and the dual Dachs– Eisenreich (DE) black carbon (BC) adsorption and OM absorption model. Use of these descriptors in a modified version of CMAQ for benzo[a]pyrene (BaP) results in different fate and transport patterns as demonstrated by regional-scale simulations performed for a European domain during 2006. The dual DE model predicted 24.1 % higher average domain concentrations compared to the HB model, which was in turn predicting 119.2 % higher levels compared to the baseline JP model. Evaluation with measurements from the European Monitoring and Evaluation Programme (EMEP) reveals the capability of the more extensive DE model to better capture the ambient levels and seasonal behaviour of BaP. It is found that the heterogeneous reaction of BaP with O3 may decrease its atmospheric lifetime by 25.2 % (domain and annual average) and near-ground concentrations by 18.8 %. Marginally better model performance was found for one of the six EMEP stations (Košetice) when heterogeneous BaP reactivity was included. Further analysis shows that, for the rest of the EMEP locations, the model continues to underestimate BaP levels, an observation that can be attributed to low emission estimates for such remote areas. These findings suggest that, when modelling the fate and transport of organic pollutants on large spatio-temporal scales, the selection and parameterization of GPP can be as important as degradation (reactivity).


Introduction
Polycyclic aromatic hydrocarbons (PAHs) are a group of lipophilic organic compounds with demonstrated carcinogenicity and potential to bioaccumulate (Diamond and Hodge, 2007;Finlayson-Pitts and Pitts, 1999;Pedersen et al., 2004Pedersen et al., , 2005;;WHO, 2003).Typically emitted to the atmosphere as a mixture of semi-to low-volatile congeners with variable carcinogenic potency, their fate and transport in the environment is difficult to assess.Current limitations involve inadequate knowledge of photochemistry (Keyte et al., 2013), air-surface exchange (Galarneau et al., 2014;Keyte et al., 2013;Lammel et al., 2009), and combustion sources (Bieser et al., 2012;Lammel et al., 2013).In addition, aluminium and steel production, along with petrogenic sources, are significant in some environments (Wenborn et al., 1999).This relationship with anthropogenic activity and population density raised awareness of their potential as health stressors, with hotspots identified in urban industrial settings around the world (Šrám et al., 2013;Zhu et al., 2015).On the other hand, observations of PAHs at remote sites and global circulation model studies indicate long-range transport (Keyte et al., 2013;Sehili and Lammel, 2007).Therefore, and due to resistance to degradation, PAHs are classified as persistent organic pollutants (POPs) by the United Nations Economic Commission for Europe (UNECE) Convention on Long-range Transboundary Air Pollution (CLRTAP) and are listed in the Protection of the Marine Environment of the North-East Atlantic (OSPAR).
A number of studies, which aimed at exploring the atmospheric fate and transport of POPs on various scales, have been presented in the last two decades.These include relatively simple mass balance models developed under the Lagrangian framework, some of which have been extended to multimedia applications for studying the accumulation and exchange of POPs between various compartments (Halsall et al., 2001;Lang et al., 2008Lang et al., , 2007;;Liu et al., 2007;Prevedouros et al., 2004Prevedouros et al., , 2008;;Van Jaarsveld et al., 1997).Eulerian grid models add yet another level of complexity regarding process detail that is accompanied by additional computational demands.Global-scale chemical transport models (Friedman and Selin, 2012;Lammel et al., 2009;Sehili and Lammel, 2007), as well as hemispheric models (Gusev et al., 2005;Hansen et al., 2006;Shatalov et al., 2005), have been applied to study the long-range transport of POPs.Moreover, coupled meteorological and air quality modelling frameworks have presented efforts to further investigate the treatment of physical and chemical atmospheric processes related to POPs in higher resolution (Aulinger et al., 2007;Bieser et al., 2012;Cooter and Hutzell, 2002;Galarneau et al., 2014;Inomata et al., 2012;Matthias et al., 2009;Meng et al., 2007;San José et al., 2013;Silibello et al., 2012;Zhang et al., 2009).Table S1 in the Supplement provides a summary of the most relevant processes covered in the respective regional model implementations and subsequent evolution tracked by the authors.While such studies represent the initial efforts to accurately model POPs on regional scale, important processes such as gas-particle partitioning, heterogeneous reactivity, and a multicompartmental approach for the relevant compounds are not always studied.The Weather Research Forecast along with the Community Multiscale Air Quality for BaP (WRF-CMAQ-BaP) model presented in this article is considered to be the first effort to examine the full spectrum of gas-particle partitioning models relevant to POPs in the same regional-scale modelling application.The model attempts to capture POP aerosol chemistry and dynamics within the CMAQ framework and can be extended to account for mass exchange between adjacent compartments for compounds for which such processes are deemed critical (i.e.soil volatilization).

Model development
The model development framework is based on the CMAQ version 4.7.1 modelling system (Byun and Schere, 2006;Byun and Ching, 1999) with the CB05 gas-phase chemistry mechanism (Yarwood et al., 2005) and the AERO4 version of the aerosol module (Binkowski, 2003).In addition, CMAQ contains modules representing advection, eddy diffusion, and in-cloud and below-cloud scavenging with precipitation.Advection and diffusion satisfy mass conservation and include removal by dry deposition, while in-cloud and precipitation processes simulate aqueous chemistry and wet deposition by cloud droplets (Byun, 1999a,b).The aerosol module (Binkowski, 2003) determines concentrations of trimodal size-distributed particulate material with diameters less than 10 µm.According to this model, each mode can be described as an internal mixture of several compounds; the mixing state on the single-particle level is not addressed, but is considered equivalent to internal mixing.Aitken and accumulation modes include particles with diameters ∼ 2.5 µm or less (PM 2.5 ), which are either emitted or produced by gasto-particle conversion processes.A thermodynamic mechanism, based on temperature and relative humidity, determines air concentrations of water, sulfate, ammonium, and nitrate aerosols in the fine modes (Nenes et al., 1998;Zhang et al., 2000).Organic and elemental carbon species are also included in these modes.The third mode represents coarse material having diameters between 2.5 and 10 µm and includes dry inorganic material emitted by natural and anthropogenic sources.Other processes like coagulation, particle growth, aerosol-cloud interactions, and new particle formation are included in the treatment.
CMAQ adaptations are required to fully address the behaviour of semi-volatile organic compounds (SOCs) in the atmosphere and include additional homogeneous gas-phase reactions, a modular algorithm to describe the mass exchange of SOCs between the gas and particulate phases, and a de-scription of deposition for each phase state of the compound of interest.For the case of benzo[a]pyrene (BaP), one gasphase species and three particulate-sorbed species (one for each aerosol mode) are introduced.It is assumed that only the gaseous molecule undergoes chemical degradation through reaction with the hydroxyl (OH) radical, leading to an unspecified product that does not undergo further chemical reaction.A summary of the physico-chemical data used in this study is presented in Table 1.Furthermore, a module that accounts for heterogeneous chemistry was included for compounds for which such reactions represent a significant loss according to experimental studies (e.g.BaP reaction with O 3 ).

Modular implementation of gas-particle partitioning (GPP) in CMAQ
The distribution of POPs between the gas and the particulate phases is widely recognized as the most important parameter in describing their fate and transport in the atmosphere.PAHs like benzo[a]pyrene can be adsorbed to mineral surfaces or elemental carbon, absorbed by organic aerosol, or absorbed by aerosol water (Lohmann and Lammel, 2004).Therefore, the aerosol module of CMAQ was enhanced with an array of models that assume different processes to determine gas-particle partitioning, i.e. an adsorption model (Junge-Pankow, Junge, 1977;Pankow, 1987), an organic matter (OM) absorption model (K oa , Harner and Bidleman, 1998), and a dual OM absorption and elemental or black carbon (EC/BC) adsorption model (Dachs and Eisenreich, 2000;Lohmann and Lammel, 2004).In addition, a method for calculating the absorption in particulate-phase water was included in the module, employed for all aerosol modes in a manner similar to the GPP processes (Aulinger et al., 2007;Cooter and Hutzell, 2002).
The new GPP algorithm for CMAQ brings together existing implementations for pesticides (Cooter and Hutzell, 2002), PAHs (Aulinger et al., 2007;San José et al., 2013), polychlorinated biphenyls (PCBs) (Meng et al., 2007), dibenzo-p-dioxins, and dibenzofurans (PCDD/Fs) (Zhang et al., 2009) in a unified module that allows for flexible GPP model combinations, making it applicable to a wide range of SOCs.The following assumptions were taken into consideration: (1) instantaneous relaxation to phase equilibrium (no mass transport kinetic limitations), (2) the compound does not irreversibly transform in the particulate phase, and (3) sorption processes do not interact.By extending the gasphase and modal aerosol treatments in CMAQ for SOCs, these assumptions allow the use of simple ratios to determine the concentration of each form, following where a i is the particulate concentration in each mode (i = 1 = Aitken, 2 = accumulation, 3 = coarse), g is the gaseous concentration, and φ i is the particulate fraction of the com-pound in mode, with a general formulation (Aulinger et al., 2007;Cooter and Hutzell, 2002).
This way, all partitioning processes, adsorption (φ ad i ), absorption in OM (φ ab i ), and absorption (solution) in aerosol water (φ aq i ), are simultaneously treated.Following the rule of mass conservation for the compound in both phases, the total concentration c tot can be introduced: leading to a system of three linear equations.
which yield the particulate concentrations a i in each mode.Subsequently, using the previous equation, the gaseous concentration g can be calculated.

Junge-Pankow adsorption model
The first GPP scheme is based on the model of exchangeable SOC adsorption to aerosols presented by Junge (1977) and later critically reviewed by Pankow (1987).
The equation above relates the mass fraction of chemical adsorbed particles in each mode i (φ ad i ) to the subcooled liquid vapour pressure of the compound (p L • , Pa) and the aerosol surface density (S i , m 2 of particle surface per m 3 of air).The parameter c J (unit, Pa m) depends on the chemical nature of both the adsorbent and the adsorbate and cannot be expected to be a constant, as confirmed by field experiments (Lammel et al., 2010).Nevertheless, a value of c J = 0.172 Pa m has been suggested (Pankow, 1987) and is often used in modelling applications.The surface area parameter, S i , which is a characteristic of the aerosol required for the Junge-Pankow (JP) model, is calculated at each time step for each mode within the aerosol module of CMAQ (Binkowski, 2003).
The empiric JP model, despite assuming adsorption to be the only relevant partitioning process, may eventually not be free from absorptive contributions (to c J ), but these should be very low.The fitting in the empiric parameterization was done for all possible vapours, many of these not soluble in octanol.No c J for BaP derived from ambient measurements has ever been published.However, we recently found in background aerosol of central Europe (Košetice 2012-2013, n = 162 aerosol samples, c J derived from experimental S/V data; Shahpoury et al., 2016) that c J for BaP is very low, at 0.01-5.46Pa cm (median 0.44 Pa cm), supporting the perception that a possible implicit absorptive contribution to c J must be very low.

K oa absorption model
The second GPP scheme follows the work of Pankow (1994Pankow ( , 2007)), with refinements from Harner and Bidleman (1998) that described how the particulate fraction of a compound could be derived from the aerosol-air partitioning coefficient K p .Following this approach for each mode i, the predicted particulate fraction φ ab i is linked to the total mass concentration of suspended particles c TSP i with the following equation: Pankow (1987,1998) demonstrated that the gas-particle partitioning coefficient, K p , and the use of p L • as its descriptor are valid both for adsorption and absorption of SOCs.In the case of absorptive partitioning, Pankow (1998) derived where f OM i is the mass fraction of OM in the particle that can absorb gaseous SOCs, R is the ideal gas constant (= 8.314 Pa m 3 mol −1 K −1 ), T (K) is absolute temperature, M OM (g mol −1 ) is the mean molecular weight of the OM phase, γ OM is the activity coefficient of the selected SOC in the OM on a mole fraction basis, and p L • (Pa) is the vapour pressure of the pure SOC (subcooled liquid in the case of solids).Partitioning coefficients normalized to the OM content were introduced to further support the absorption hypothesis (Finizio et al., 1997).Pankow (2007) and Harner and Bidleman (1998)

replaced p L
• by the octanol-air partitioning coefficient, K oa .
where M oct , M OM (g mol −1 ) are the molecular weights of octanol (= 130 g mol −1 ) and the OM phase, γ oct , γ OM are the activity coefficients of the SOC in octanol and OM, respectively, and ρ oct is the density of octanol (0.82 kg L −1 ).Assuming that octanol imitates organic matter in PM, Harner and Bidleman (1998) suggested that the ratio of γ oct /γ OM and M oct /M OM can be assumed to be 1.However, it was later suggested that M OM could be much higher, particularly in secondary organic aerosols containing polymeric structures (Kalberer et al., 2004); a mean value of 500 g mol −1 was later suggested by Götz et al. (2007), which results in M oct /M OM of 0.26.Under the assumption that γ oct /γ OM and MW oct /MW OM = 1, the equation above can be simplified to the following:

Dual OM absorption and black carbon adsorption model
Finally, the third GPP scheme extends the previous OM absorption model to include adsorption onto particulate soot, according to evidence showing the measured K p to be exceeding the predicted K p based solely on absorption in OM (Dachs and Eisenreich, 2000;Fernandez et al., 2002;Lohmann and Lammel, 2004;Ngabe and Poissant, 2003).
Considering the contribution of those two additive processes, a modified partition coefficient can be employed.
where f BC i is the fraction of black carbon in the particulate matter, a atm−BC and a soot are the available surfaces of atmospheric black carbon (BC) and diesel soot, respectively, K SA is the partitioning coefficient between diesel soot and air.A density value of 1 kg L −1 has often been assumed for BC; however values in the range of 1.7-1.9kg L −1 have been measured and recommended by Bond and Bergstrom (2006).Thus, the equation can be simplified.
This method is subject to uncertainties (Goss, 2004), but is accepted and suitable (Dachs et al., 2004 and references therein).Soot-air partition coefficients (K SA ) were calculated as the ratio of soot-water adsorption constants K SW and the inverse Henry's law constant (H ), with K SW values adopted from Bärring et al. (2002).As a step to address uncertainties contributed by different methods in calculating K SA , a thermodynamic estimation model suggested by van Noort ( 2003) was also tested: The soot-specific surface area was set to 18.21 m 2 g −1 , derived as the geometric mean of surface areas for traffic, wood, coal, and diesel soot (i.e.59.4, 3.6, 8.2, and 62.7 m 2 g −1 , respectively) (Jonker and Koelmans, 2002).Additionally, in order to estimate uncertainties arising from K oa and its temperature dependence, two different parameterizations were tested, the first based on the work of Beyer et al. (2000) and the second following the temperature-dependent parameterization suggested by Odabasi et al. (2006).

Dissolution into aerosol water
An assumption of the CMAQ aerosol module is that organics influence neither the water content nor the ionic strength of the system.While this assumption may not be valid for many atmospheric aerosols, sufficient basic data are not available to treat the system in a more complete way.In a similar manner, the dissolution of PAHs in the aqueous electrolyte was calculated using the inverse Henry's law constant (Aulinger et al., 2007;Cooter and Hutzell, 2002;Lohmann and Lammel, 2004).Lohmann and Lammel (2004) derived the partitioning coefficient K aq p,i , although they conclude that as long as a solid-air interface exists, the contribution of dissolved PAHs to the gas-particle partitioning will always be negligible.To account for this, a "wet aerosol" switch was introduced and associated with the relevant partitioning coefficient.
where f wa,i is the mass fraction of aerosol water in each mode i and K wa is the water-air partitioning coefficient, which is equal to H .The aerosol water content is computed in CMAQ based on a variation of the Zdanovskii, Stokes, and Robinson (ZSR) method (Byun and Ching, 1999;Kim et al., 1993;Stokes and Robinson, 1966).As suggested by Aulinger et al. (2007) for CMAQ, only if the ratio of ammonium sulfate to aerosol water drops below the solubility of ammonium sulfate, the aerosol is treated as wet and absorption into aerosol water takes effect, while adsorption to inorganic material is switched off.

Gas-phase and heterogeneous reactions
Reactions with the hydroxyl radical (OH), ozone (O 3 ), and the nitrate radical (NO 3 ) may determine the atmospheric fate and long-range transport potential of SOCs.Gas-phase reactions with OH are believed to be the dominant loss pathway for semi-volatile PAHs (3-4 rings) like anthracene and pyrene (Atkinson and Arey, 1994;Finlayson-Pitts and Pitts, 1999;Keyte et al., 2013).On the other hand, in the case of non-volatile PAHs (5 and more rings, like BaP), loss due to heterogeneous reaction with ozone is believed to dominate (Finlayson-Pitts and Pitts, 1999).Experimental studies of the kinetics of PAHs sorbed to particles indicate strong influences of mixing state (single-particle anisotropy), morphology, and phase state of the particles (Kwamena et al., 2004;Perraudin et al., 2007;Pöschl et al., 2001;Shiraiwa et al., 2009;Zhou et al., 2012).Yet, detailed heterogeneous chemistry models would require corresponding additional parameters and a more sophisticated surface layer model (Kaiser et al., 2011;Shiraiwa et al., 2010;Springmann et al., 2009).However, in the studies by Kahan et al. (2006) and Kwamena et al. (2004Kwamena et al. ( , 2007)), the reactivity could be well described by a Langmuir-Hinshelwood mechanism (i.e. one species [BaP] is adsorbed to the surface while the second species, ozone, is in phase equilibrium).This approach is followed in the CMAQ aerosol module by introducing the degradation rate coefficient k in the following form: where k max is the maximum rate coefficient (0.060 ± 0.018 s −1 ), K O 3 is the ozone gas to surface equilibrium constant (0.028 ± 0.014 × 10 −13 cm 3 ), and c O 3 is the concentration of ozone (mol cm −3 ).As this reaction neglects the substance fraction not accessible for gaseous O 3 molecules (so-called burial effect; Zhou et al., 2012), an upper estimate for BaP reactivity is simulated.

Dry and wet deposition
SOCs and their degradation products are removed from the atmosphere through dry and wet deposition.Dry deposition is modelled in CMAQ as a one-way flux out of the lowest layer of the atmosphere.Particulate forms have dry deposition velocities determined by Brownian diffusion, turbulence, and sedimentation.Mass and diameter of the mode are used to estimate the sedimentation velocity (Binkowski, 2003;Byun and Ching, 1999).The deposition velocity of gaseous molecules is based on an electrical resistance analogue that includes aerodynamic (R a ), quasi-laminar boundary layer (R b ), and surface canopy (R c ) resistance (Pleim et al., 1999).R a is computed using similarity theory with heat flux.R b depends on land-use-specific friction velocity and molecular characteristics of gases.R c consists of several components, including bulk stomatal, dry cuticle, wet cuti-cle, ground, and in-canopy aerodynamic resistance (Pleim et al., 1999).Cuticle and ground resistances are based on ozone and sulfur dioxide observations and neglect absorption into organic plant or soil material.Values for time-dependent R c parameters are determined via land use and mesoscale meteorological models, which provide fractional vegetation cover, leaf area index, fractional leaf area wetness, and leaf stress associated with radiation, root zone soil moisture, temperature, and humidity (Pleim and Xiu, 1995).
Wet deposition of gaseous molecules occurs through incloud and below-cloud scavenging.Both processes are assumed to have the same scavenging efficiency and dependency on Henry's law, water content of the cloud, and precipitation rate (Chang et al., 1986).Wet scavenging parameterization of the particulate-phase molecules follows the method of Roselle and Binkowski (Byun and Ching, 1999) and depends on the aerosol mode (size) and time required to remove all liquid water from a cloud volume.No ice-phase scavenging is included in the fourth version of the aerosol module in CMAQ.

Model implementation and study design
The model domain included the whole of Europe and neighbouring countries as shown in Fig. 1.The domain is within the area covered by the UNECE-CLRTAP (and its POP protocol, covering PAHs) and was specifically chosen to exploit recently prepared high-resolution emission data (Bieser et al., 2011) and make use of the only monitoring network in the region (EMEP; see below).Computational simulations were performed on a 36 × 36 km Lambert Conformal Conic grid using offline hourly meteorological fields, generated by the Weather Research Forecasting (WRF) model version 3.4.1 using GFS meteorological reanalysis data (spatial resolution 0.5 • × 0.5 • ) as input.The parameterization of WRF included the following schemes: Milbrandt-Yau double-moment 7-class microphysics (Milbrandt and Yau, 2005a, b), rapid radiative transfer model (RRTMG) longwave and shortwave scheme (Iacono et al., 2008), asymmetric convective model of the planetary boundary layer (PBL) (Pleim, 2007), Pleim-Xiu surface layer model (Pleim and Xiu, 2003;Xiu and Pleim, 2001), and the improved version of Grell-Devenyi (Grell, 1993;Grell and Devenyi, 2002) ensemble scheme for cumulus parameterization.On the vertical dimension, the domain was resolved in 36 layers, following the eta coordinate system.Subsequently, CMAQready meteorological input fields for 2006 were prepared using the Meteorology-Chemistry Interface Processor (MCIP) (Otte and Pleim, 2010).Computational simulations were performed in a distributed computing infrastructure operated and managed by the Czech National Grid Organization Meta-Centrum NGI.

Emissions
Capturing the spatio-temporal pattern of emissions is crucial to determine the fate and transport of PAHs.Several attempts have been carried out during the past decades, with varying spatial and temporal resolutions in historical, current, and future projection terms for Europe (Bieser et al., 2012;Pacyna et al., 2003).Currently, four database sources with gridded BaP emissions exist for Europe and Russia: the TNO database (Denier van der Gon et al., 2006), the EMEP database (Mantseva et al., 2004), a dataset for POPCYCLING-Baltic project (Pacyna et al., 2003), and a global emission inventory for the year 2008.The SMOKE for Europe (SMOKE-EU) model (Bieser et al., 2011) output used in this study follows a methodology that is based on the TNO database and has been evaluated in a number of subsequent CMAQ applications (Aulinger et al., 2011;Bewersdorff et al., 2009;Matthias et al., 2008).In order to spatially disaggregate BaP emissions in a more realistic way, SMOKE-EU improves on using a linear dependency on the population density as surrogate factor by introducing different relevant relationships based on source classification and meteorological parameters.For example, the relationship to wood availability and the level of urban agglomeration have been introduced in an effort to include the effect of burning biomass.The resulting spatial distribution and annual average emission flux of BaP for our domain during the year 2006 is shown in Fig. 1.In addition, BaP emissions are shaped temporally using diurnal and weekly functions.As residential combustion is associated with the majority of BaP emissions, a strong dependence on the season and spatial variability is expected.To account for this dependence, heating degree days have been introduced as a proxy (Aulinger et al., 2011).A major strength of the SMOKE-EU inventory is the provision of CMAQ-ready emissions for all species involved in the desired chemical mechanism (e.g.CB5) and aerosol module (e.g.aero4).Following the trimodal representation of aerosol in CMAQ, 0.1 % of the emissions of BaP was allocated to the Aitken mode, while the rest was emitted in the accumulation mode, as this is the default treatment with all primary organic aerosol species in the model.

Evaluation data and sensitivity simulations
The spatio-temporal coverage of BaP in air and deposition in Europe and Russia is limited, in many regions sporadic or not available at all.For 2006, concentration in air is available from eight stations (n = 173), while deposition information is available from six stations of the EMEP monitoring network (Aas and Breivik, 2012).An overview of the data, including site codes, geographical information, and the sampling strategies, is provided in Table S2.It is evident that during 2006, only 4 of the total 10 sites provide information on concentration and deposition simultaneously.Additionally, the sampling strategy varies widely among the stations, leading to a non-homogeneous pool of data.For example, at the Košetice (CZ0003R) and Niembro (ES0008R) stations, BaP is measured for 24 h once a week, while in Aspvreten (SE0012R) and in Pallas (FI0036R) samples cover one week and are taken once every month.On the other hand, the Latvian sites (LV0010R, LV0016R) are sampled on a monthly basis, while for the High Muffles location (GB0014R) a 3month sample is obtained four times in 2006.Finally, the sites located in Germany (DE0001R, DE0009R) provided information only in terms of wet deposition and on a monthly basis during 2006.
Based on the set of GPP mechanisms described before, the effect of each parameterization was evaluated with the EMEP measurements by incrementally testing the models involved, while keeping the same meteorological and emission inputs.Table 2 summarizes the five annual scenarios that were developed for 2006: (1) Junge-Pankow, (2) scenario 1 with water dissolution, (3) scenario 2 with the Harner-Bidleman (HB) K oa absorption model, (4) scenario 2 with dual OM absorption and black carbon adsorption model Dachs-Eisenreich (DE) scheme, and (5) scenario 4 supplemented by a heterogeneous reactivity module.A spin-up time of 7 days was allowed for each simulation in order to maintain a low influence of the initial conditions.Comparisons with measurements of BaP in ambient air were produced using the statistical computing framework of the R language, and the model performance was evaluated, building on a wide range of statistics calculated using the Models-3 and openair packages (Carslaw and Ropkins, 2012).
of large cities (Moscow, Paris, Vienna, Madrid, and Istanbul, among others).In addition, annual mean BaP concentrations according to the fully expanded scenario 5 (KW) are depicted in Fig. 2b.Comparing the output from the aforementioned simulations revealed higher BaP levels with scenario 5, which can be associated with the geographic distribution of emissions.This is the first indication that the chosen partitioning approach has a significant impact on the model results for total BaP concentrations.However, it is also evident that in southern Europe and the Mediterranean region, the BaP losses associated with the addition of an upper estimate of heterogeneous reactivity with O 3 (scenario 5, Fig. 2b) may have a limiting effect on BaP long-range transport (LRT) potential.This can be supported by the prominent BaP gradients calculated over the Mediterranean Sea, acting as a sink that is driven by low BaP emissions and elevated O 3 levels.The overall model patterns show a good agreement with recent regression methods mapping BaP levels that provide with a partial coverage of our EU domain.
With respect to previous studies involving simulations of BaP over Europe, we find lower BaP concentrations, but within the same order of magnitude as those estimated by initial CMAQ applications (Aulinger et al., 2007;Bewersdorff et al., 2009;Matthias et al., 2009).This supports the perception that BaP emissions over Europe have decreased substantially during the recent decades, a feature improved in subsequent evolution of the emissions model which led to better agreement with the more recent CMAQ applications (Bieser et al., 2012).A closer look over Italy and the Iberian Peninsula reveals a good agreement with previous studies depicting BaP distributions for these regions (San José et al., 2013;Silibello et al., 2012).Finally, discrepancies over eastern Europe indicate that a more elaborate emissions inventory (i.e.domestic heating adjustments, fire activities) would be beneficial to address the elevated BaP levels predicted by recent hybrid modelling approaches (Guerreiro et al., 2016a).

Effect of GPP parameterizations and heterogeneous reactivity with O 3
Figure 3 and Table 3 illustrate the differences resulting from the incrementally tested scenarios (Table 2), expressed as annual mean BaP concentrations over the entire domain of interest and the individual 28 EU member states.As expected, the effect of dissolution into aerosol water (Fig. 3a) is limited throughout most of the domain, with coastal urban areas being the most influenced (i.e.Istanbul, North Sea).According to scenario 2, BaP dissolution in aerosol liquid water led to a reduction of less than 0.36 % of ground-level BaP in Europe compared to the baseline JP scheme (Table 3).Introducing the HB absorption scheme (scenario 3) resulted in ∼ 119.2 % higher average BaP concentrations (Fig. 3b) with extremes observed in areas associated with emission sources in central Europe (i.e.Po Valley) and large cities (i.e.Moscow, Paris, Vienna, and Istanbul).Expressed in terms of domainwide atmospheric lifetimes (defined as τ = b tot /F em , where b tot denotes the total atmospheric burden and F em the total emissions), scenario 3 leads to a 210 % increase for BaP over scenario 2. Similarly, shifting to the DE scheme (scenario 4) appears to have an additive effect on near-ground BaP concentrations with total lifetimes reaching a domain-average increase of 13 %.Compared to scenario 3, scenario 4 affects a larger area of the domain (Fig. 3c), even though the magnitude of the effect is less pronounced due to its association with the availability of black carbon.The domain mean of modelled elemental carbon is 0.0212 µg m −3 and the spatial distribution matches its emission pattern (Fig. S3).These levels are significantly lower than in the field studies in the urban air of south-western Europe (EEA, 2013;Milford et al., 2016) and may subsequently lead to a hindered performance of the DE scheme.Finally, the introduction of heterogeneous chemistry (scenario 5) had, as expected, a substantially negative effect on the BaP concentrations within our domain, culminating in major conurbations (i.e.Po Valley, Moscow, and Istanbul; see Fig. 3d).Heterogeneous degradation due to the reaction with O 3 accounts for an approximate 18.8 % reduction of the mean ground-level BaP and a 25 % reduction of its atmospheric lifetime in the European domain (Table 3).As mentioned in Sect.4.1, the BaP concentration gradients reveal a strong spatial link with emission sources in the south, as well as high O 3 concentrations typically found along the coasts of the Mediterranean Sea (Fig. S4).CMAQ calculates ozone patterns that are in accordance to the findings of recent analyses, exhibiting a similar performance to measurements at the European level (de Smet et al., 2009).

Seasonality and relevant parameters affecting BaP fate and transport
Seasonally disaggregated distributions of BaP, as simulated in scenario 4, are presented in Fig. 4. Strong emissions during wintertime and autumn have an apparent effect on the average ground-level concentrations.In addition, a normalized effect of long-range transport for BaP can be observed in the vicinity of Istanbul and other coastal metropolitan areas, over parts of the large water bodies of the Mediterranean and the Black Sea.This scenario, however, does not include the treatment of the heterogeneous reaction with ozone and, therefore, provides an upper-bound estimate of the LRT potential of BaP.On the other hand, Fig. 5 demonstrates the seasonal effect of the heterogeneous reaction of BaP with ground-level ozone on BaP concentrations.Distributions of BaP concentration differences reveal a stronger reduction during winter that is influenced by central and southern European emission sources and midlatitude ozone availability.Interestingly, the differences between winter and summertime appear to be comparable in terms of magnitude, despite the substantially weaker sources during summer that shift the process to the south.
An important factor for the performance of the BaP model is the mass concentration and composition of the aerosols in the three different modes and their relationship with the temperature of the ambient air. Figure 6 illustrates the aerosol distributions as simulated at model cells representing three different settings (urban: the city of Vienna, suburban: CZ0003R, remote: FI0036R) during a typical winter and summer day.As we move further away from the vicinity of the sources, we notice not only a decrease of the over-all BaP levels but also a change in the distribution influenced by each individual simulated scenario.While the basic Junge-Pankow implementation results in higher coarsemode BaP concentrations, the HB and DE schemes accentuate the strength of the accumulation mode where most of the BaP aerosol mass is emitted.The disproportionately higher coarse-BaP fractions calculated at the remote grid cells during the summertime can be explained by conditions favouring long-range transport.
Regarding the parameterizations of the HB and DE schemes, the calculation of the parameter K oa has a direct effect on the relevant model results (scenarios 3 and 4).Use of the regression parameters for the semilogarithmic temperature-dependent form proposed by Odabasi et al., (2006) results in a mixed effect over its predecessor scheme, which is seasonally disaggregated and presented in Fig. S5.These figures suggest that the partitioning based on each K oa approach is sufficiently similar for most of the domain during the winter.However, differences between the two expressions with respect to simulating summertime BaP concentrations over central Europe were noted.In addition, simulations that address the uncertainties contributed by the method for estimating K SA were performed for the winter period of 2006 using the DE scheme.Figure S6 reveals similar results, with BaP levels slightly higher in cells with strong emissions sources.However, with the exception of Moscow, the average difference BaP concentration in urban cells did not exceed 0.02 ng m 3 when K SA was calculated based on the thermodynamic estimation model suggested by van Noort (2003).Similarly to K oa , results suggest that the partitioning based on each K SA approach is sufficiently similar.

Model evaluation with EMEP measurements
The EMEP network database and model output from all five scenarios were imported in the R language framework to create measurement-model pairs for further analysis.Figure 7 illustrates the resulting time series plots for selected EMEP sites against CMAQ output extracted from scenario 4.This figure also demonstrates differences in terms of sampling protocols among individual sites, with the extreme of the High Muffles station (4 samples year −1 -Table S2).Aggregate temporal profiles of the modelled and measured BaP concentrations across all sites depicted in Fig. 8 reveal a strong seasonal pattern that is captured well by the model.However, BaP concentrations appear to be underestimated during the colder months.This effect can be attributed to the strength of emissions, which is probably related to inadequate coverage of residential heating in current inventories (wood/biofuel burning).Figure S7 illustrates this by comparing annual average BaP emission fluxes between model cells enclosing EMEP monitoring site locations and selected metropolitan areas of Europe during 2006.While the deficiency of PAH inventories is known, it may limit the usefulness of background measurements (such as those in the EMEP network) in assessing the overall model performance under different GPP scenarios.
It is evident from the metrics obtained across all sites (Table 4 and Taylor diagram summary Fig. S8) that the increased complexity in GPP formulation results in better agreement with the EMEP measurements.For reasons mentioned in Sects.2.2-2.3 and 4.1.1,deviations from this observation are the simulated scenarios that introduce the dissolution to aerosol water (scenario 2) and degradation due to the heterogeneous reaction with ozone (scenario 5).Overall, the dual model (DE) employed in scenario 4 showed the best agreement with the measurements.The DE parameterization also performed best when compared with the JP and HB schemes  on global-scale applications (Lammel et al., 2009).However, when looking at site-specific performance metrics (Table S9), the Košetice location (CZ0003R) reveals a better agreement with the fully expanded scenario 5, which treats degradation by ozone.The Košetice site is the closest to conurbations among the EMEP sites available during 2006.This could indicate that the heterogeneous chemistry, which is implemented as an upper estimate in this model, tends to underestimate BaP levels further away from the sources.The best model performance was observed at the Pallas station (FI0036R), where scenario 4 was capable of matching 50 % of the measurements within a factor of 2. For the remaining remote sites in Finland, the simulation under scenario 5 calculated underestimated BaP levels.Based on the index of agreement (IOA) metric, a ranking of the overall model performance similar to the Pallas site was observed for the rest of the monitoring sites.Having compared sitespecific performance to previous modelling efforts over Europe (Aulinger et al., 2007;Matthias et al., 2009), we reach the same conclusion as pointed out with the spatial distributions; i.e. we attribute the lower-performance metrics presented here to the substantially weaker emission sources in recent years, hence the better agreement with the more recent SMOKE-EU/CMAQ application (Bieser et al., 2012).
With respect to other annual simulations, we observed lower average BaP levels over the Piedmont region compared to the FARM model, which were higher than the EMEP (MSC-E) output used in the same study (Silibello et al., 2012).Finally, a qualitative analysis of the deposition measurements for the few stations that measured BaP showed that CMAQ calculated plausible patterns of dry and wet deposition (see Fig. S10).

Conclusions
This study presented the development of a new modelling system for investigating the dynamics of transport and gasparticle partitioning schemes for POPs at the regional scale over Europe.The implementation of this model is based on the WRF-CMAQ model framework with additional algorithms for GPP schemes and a module that accounts for heterogeneous reactivity of BaP with ozone.Predictions from the WRF-CMAQ-BaP model were compared to measurements obtained from the EMEP monitoring network.Model evaluation has revealed satisfactory agreement to the mea- surements and performance metrics similar to those of previous studies with significantly higher measurement availability.The results presented in this work suggest that the new expanded model is able to simulate the ambient levels of BaP fairly well.It is found that BaP distributions in Europe are very sensitive to the choice of GPP models.While this statement has been also indicated at the global scale (Lammel et al., 2009), GPP has not been adequately explored by previous regional-scale studies (see Table S1 and references therein).We conclude that the dual OM absorption and black carbon adsorption (DE) parameterization (scenario 4-5) offer a better performance, as BaP predictions tend to be closer to the measurements.Introducing an independent upper estimate of heterogeneous BaP reactivity is found to significantly reduce BaP levels throughout the domain, with a re-duction pattern that follows the spatial distribution of ozone and emission sources.Despite the limited BaP monitoring network, a better agreement was observed when complementing WRF-CMAQ-BaP with the upper estimate for heterogeneous reactivity (scenario 5) at certain EMEP sites.Uncertainties and limitations of the current emission inventory estimates of BaP, OM, and BC, along with insufficient laboratory data, are of major concern, hampering efforts to fully evaluate the long-range transport potential, the effect of photochemistry, and interactions of such compounds.In addition, certain simplifications in the treatment of coarse PM and OM within the aerosol module of CMAQ (Binkowski, 2003;Mathur et al., 2008) are posing limitations on fully exploiting the advanced empirical gas-particle partitioning models (K oa , Dachs-Eisenreich) in this specific mode.An implicit overlap of the empiric adsorption model (JP) with the K oa (or another absorption) model (see Sect. 2.1.1)may artificially enhance partitioning, something that has not been addressed by previous modelling efforts.However, such discrepancy is expected to be negligible for BaP, especially when compared to existing uncertainties of output aerosol parameters propagating from emission estimates and from the current design and state of development of the CMAQ aerosol module (AERO).The modelling approach presented here allows a simultaneous estimation of organic pollutants, including semivolatiles (such as most POPs) within the CMAQ framework, and can be used as a supplementary component of a population exposure modelling system.Although similar systems are currently explored, they rely heavily on observation data while using a very basic air quality model structure (Guerreiro et al., 2016a, b).At this time, annual exposure estimates are more useful due to the low accuracy accompanying highly time-resolved models and observations.As previ-ous efforts to evaluate the use of mesoscale models for BaP dispersion studies also conclude, changes in emission patterns and their strengths are of particular importance.While there is undeniable effort to shift from fossil fuels to renewables which has led to updated emission inventories, domestic heating across Europe has also a complex relationship to economic factors and human activities that varies tremendously between countries and fuel types (Lohmann et al., 2006;Saffari et al., 2013).Besides the need for updated emission inventories for the eastern Europe, extended simulations are suggested for years with higher measurement availability and density.Following modelling, studies should focus on quantifying the long-range transport potential and examining the hypothesis that secondary organic aerosol (SOA) may facilitate PAH transport by utilizing more sophisticated aerosol and heterogeneous chemistry parameterizations or submodels (i.e.accounting also for the burial effect; see Sect.2.2).More specifically, recent evidence suggests that, in cold and dry air, accessibility of PAHs in OM is reduced (due to low diffusivity), which might explain the apparent inconsistency of high LRT potential (BaP levels in the Arctic) on one hand and the relative high heterogeneous reactivity measured in the laboratory on the other (Friedman et al., 2014;Zelenyuk et al., 2012;Zhou et al., 2012Zhou et al., , 2013)).To understand the atmospheric lifetime of PAH, interaction between PAH and SOA, and progression of the chemical composition of ambient OM, a new need for additional studies quantifying GPP and LRT potential under a wide range of atmospheric conditions is emerging.In view of such findings, the WRF-CMAQ-BaP modelling system should be extended to study a wide range of additional organic pollutants and processes (e.g.multicompartmental cycling, biodegradation, heterogeneous chemistry).

Data availability
The WRF-CMAQ-BaP output data used in this study are archived in MetaCentrum and are available upon request from the authors (efstathiou@recetox.muni.cz).
The Supplement related to this article is available online at doi:10.5194/acp-16-15327-2016-supplement.

Figure 1 .
Figure 1.Study domain and 2006 annual grid cell mean BaP emission flux (g s −1 ) at the surface.

Figure 2 .
Figure 2. Annual average surface-level BaP concentrations [log c BaP (ng m −3 )] during 2006 plotted in logarithmic scale for (a) the Junge-Pankow (JP) scenario 1 and (b) for the fully expanded GPP scenario 5, which includes a heterogeneous reaction with ozone.

Figure 5 .
Figure 5. Average difference of surface-level BaP concentrations (ng m −3 ) between scenario 4 and the fully expanded GPP scenario 5, which includes the heterogeneous reaction with ozone during the winter and summer months of 2006 (scenario 4-scenario 5).

Figure 6 .
Figure 6.Distribution of BaP mass [log c BaP (ng m −3 )] in aerosols of three different modes in three different cells [urban (Vienna), suburban (CZ0003R), and remote (FI0036R)] for each scenario as calculated during a winter (20 January) and a summer (20 July) day.

Figure 7 .
Figure 7. Model predicted and observed concentration (ng m −3 ) time series at selected EMEP monitoring sites during 2006 for the dual GPP model scenario 4.

Figure 8 .
Figure 8. Model predicted and observed monthly BaP concentrations (ng m −3 ) across all EMEP sites in 2006 for the dual GPP model scenario 4.

Table 1 .
Physicochemical parameters and reaction rate coefficients used in the WRF-CMAQ-BaP model.

Table 2 .
Summary of the simulated scenarios designed to test gas-particle partitioning and heterogeneous reactivity of BaP.

Table 3 .
Average ground-level BaP concentrations, c BaP (pg m −3 ) over the entire domain and in individual 28 EU member states for all simulated scenarios and concentration changes c BaP (%) between subsequent scenarios.

Table 4 .
Comparison between modelled and measured BaP concentrations for the entire 2006 EMEP dataset (n = 173) along with performance metrics and ranking across all sites.
FAC: fraction of predictions within a factor of two, MB: mean bias, MGE: mean gross error, NMB: normalized mean bias, NMGE: normalized mean gross error, RMSE: root mean squared error, r: correlation coefficient, COE: coefficient of efficiency, IOA: index of agreement.