Global simulation of semivolatile organic compounds – development and evaluation of the MESSy submodel SVOC (v1.0)

The new submodel SVOC for the Modular Earth Submodel System (MESSy) was developed and applied within the ECHAM5/MESSy Atmospheric Chemistry (EMAC) model to simulate the atmospheric cycling and air–surface exchange processes of semivolatile organic pollutants. Our focus is on four polycyclic aromatic hydrocarbons (PAHs) of largely varying properties. Some new features in input and physics parameterizations of tracers were tested: emission seasonality, the size discretization of particulate-phase tracers, the application of poly-parameter linear free-energy relationships in gas–particle partitioning, and re-volatilization from land and sea surfaces. The results indicate that the predicted global distribution of the 3-ring PAH phenanthrene is sensitive to the seasonality of its emissions, followed by the effects of considering re-volatilization from surfaces. The predicted distributions of the 4-ring PAHs fluoranthene and pyrene and the 5-ring PAH benzo(a)pyrene are found to be sensitive to the combinations of factors with their synergistic effects being stronger than the direct effects of the individual factors. The model was validated against observations of PAH concentrations and aerosol particulate mass fraction. The annual mean concentrations are simulated to the right order of magnitude for most cases and the model well captures the species and regional variations. However, large underestimation is found over the ocean. It is found that the particulate mass fraction of the benzo(a)pyrene is well simulated, whereas those of other species are lower than observed.

1 Introduction

The atmospheric cycling of semivolatile organic compounds (SOCs) is particularly complex because of partitioning across phases and air–surface exchange processes, including multihopping (or “grasshopper effect”; Semeena and Lammel, 2005) and accumulation in ground compartments such as seawater, soil, vegetation, and ice/snow. Many SOCs do resist degradation in environmental compartments, and hence are persistent. In the regulation of chemical substances and in international chemicals legislation (e.g., UNEP, 2017), model-based quantifications of the overall environmental residence times (persistence) and the long-range transport potentials are requested or encouraged to be applied.

Global and regional distribution and transport of SOCs has been studied using multimedia fate (box) models and chemistry transport models (CTMs) (Scheringer and Wania, 2003). The multimedia models describe the whole or part of the globe as a few zones of homogeneous environmental characteristics (Wania and Mackay, 1999; Mackay, 2010). These models are used as tools to assess the influences of environmental parameters and change in pollutant levels in multiple compartments (Dalla Valle et al., 2007; MacLeod et al., 2005; Lamon et al., 2009). On the other hand, CTMs generally imply the application of three-dimensional Eulerian models coupled with surface and chemistry modules (e.g., Ma et al., 2003; Hansen et al., 2004; Malanichev et al., 2004; Gusev et al., 2005; Semeena et al., 2006; Gong et al., 2007; Friedman and Selin, 2012; Galarneau et al., 2014; Shrivastava et al., 2017). The addition of a surface module aims to describe air–surface exchange processes and biogeochemical cycles of contaminants, whereas a chemistry module describes the changes in air concentrations due to phase partitioning and chemical transformations. Compared to the multimedia models, CTMs have better spatial and temporal resolutions but require more computational effort. They are suitable for use to investigate the variability and episodic character of environmental fate and transport. To date, pollutants addressed in model studies were persistent organic pollutants, such as dichlorodiphenyltrichloroethane (DDT), polychlorinated biphenyls (PCBs), hexachlorocyclohexanes (HCHs), polycyclic aromatic hydrocarbons (PAHs), and more recent so-called emerging pollutants (e.g., MacLeod et al., 2011).

The sensitivity of distributions to specific processes of SOC cycling and related input parameters has been the focus of CTM-based studies (Semeena et al., 2006; Sehili and Lammel, 2007; Friedman and Selin, 2012; Galarneau et al., 2014; Thackray et al., 2015). Sehili and Lammel (2007), for instance, suggest that the gas–particle partitioning and particulate-phase oxidation scenarios have significant influences on the long-range atmospheric transport of PAHs. This finding is supported by Friedman and Selin (2012), who furthermore concluded that the effects are higher than those of irreversible partitioning and of increased aerosol concentrations.

This study presents the new multicompartment module (submodel) SVOC for the Modular Earth Submodel System (MESSy; Jöckel et al., 2006, 2010). MESSy provides a modular framework for simulations accounting for various degrees of complexity and to facilitate continuous future submodel improvements. The submodel has been applied using the general circulation model ECHAM5 (Roeckner et al., 2003, 2006) as a base model. In connection with the ECHAM5/MESSy Atmospheric Chemistry (EMAC) model, SVOC encompasses a 3-D atmosphere and 2-D surface compartments (soil, vegetation, snow, and ocean mixed layer), and considers multicompartment fate and exchange processes, such as emission, phase partitioning, wet and dry deposition of gases and particles, degradation, and air–surface gas exchange, including re-volatilization. SVOC is developed and intended to be applied for the study of all potentially re-volatilizing and gas–particle-partitioning (hence semivolatile) compounds. Nevertheless, the focus of this submodel development is the global distribution of four PAH species of largely varying properties. PAHs enter the atmospheric environment as by-products of all technological combustion processes (Shen et al., 2013) and of open fires (Gullett et al., 2008). They are ubiquitous pollutants of particular environmental and health concern (WHO, 2003; Laender et al., 2011; Lammel, 2015) and due to their continuous global emissions. Here we describe the submodel development, compare the results to observations, and assess the significance of four model features to PAH distributions and fate. These features are the temporal resolution of emissions, the size discretization of particulate-phase tracers (bulk or modal), the choice of the gas–particle partitioning scheme, and re-volatilization from surfaces.

2 Methods

2.1 Model descriptions

The global model applied in this study is the ECHAM5/MESSy Atmospheric Chemistry Climate model (EMAC), a three-dimensional Eulerian model for the simulations of meteorological variables, gases, aerosols, clouds, and other climate-related parameters. EMAC combines the general circulation model ECHAM5 (here version 5.3.02) (Roeckner et al., 2003, 2006) with the Modular Earth Submodel System (MESSy version 2.50; Jöckel et al., 2006, 2010). The atmospheric component of ECHAM5 derives four prognostic variables, namely vorticity, divergence, temperature, and the logarithm of surface pressure in truncated series of spherical harmonics, whereas specific humidity, cloud water, and cloud ice are represented in grid point space. MESSy provides a modular framework to define atmospheric dynamics, chemistry, transport, and radiative transfer processes. For a more detailed description of the EMAC model, evaluation and relevant studies, refer to Jöckel et al. (2006, 2010) and http://www.messy-interface.org (last access: 12 July 2019). The list of MESSy process-based modules, hereinafter submodels, applied in this study are summarized in Table 1.

Lauer et al. (2007)Pozzer et al. (2006)Roeckner et al. (2006)Tost et al. (2006b, 2010)Tost (2006)Kerkweg et al. (2006a)Pringle et al. (2010)Landgraf and Crutzen (1998)Tost et al. (2007)Sander et al. (2011)Kerkweg et al. (2006b)Kerkweg et al. (2006b)Roeckner et al. (2006); Jöckel et al. (2006)Tost et al. (2006a)Kerkweg et al. (2006a)

Table 1 Summary of MESSy process submodels used in this study.

The new MESSy submodel SVOC for simulating the fate and cycling of SOCs in the global environment is presented. Processes involved in the submodel include gas–particle partitioning, volatilization from the surface, dry and wet depositions, and chemical and biotic degradations. These processes are connected to other MESSy submodels. For example, deposition of gas-phase SOCs are calculated by the submodels SCAV and DDEP, aerosol microphysics by GMXe, gas-phase chemistry mechanisms by MECCA, and ocean–air flux exchange by AIRSEA. Figure 1 illustrates the SVOC structure within the EMAC system and its interactions with other MESSy submodels. More details on some process parameterizations are given in the following section. A user manual can be found in the Supplement with the list of submodel input and output variables.

Figure 1 Overview of EMAC-SVOC model structure, the cycling processes in SVOC submodel, and its interaction with other MESSy submodels. SMIL (submodel interface layer) and SMCL (submodel core layer) are components of MESSy coding standard, see Jöckel et al. (2006) for further details.

2.2 Parameterizations of cycling processes in multiple compartments

2.2.1 Representation of SOC in particulate phase

The parameterizations of aerosol microphysical processes for SOCs such as gas-to-particle partitioning and dry and wet deposition depend on the way the particulate phase is represented in the model. Here, there are two approaches employed in the submodel to represent the particulate-phase SOC: (1) it is assumed as a bulk species or (2) the particle sizes are resolved into n continuous (modal) distributions. The former will be hereinafter referred to as the bulk scheme and the latter referred to as the modal scheme. In the modal scheme, n is equal to the seven log-normal modes of the GMXe submodel, four with hydrophilic coating (ns: nucleation soluble, ks: Aitken soluble, as: accumulation soluble, cs: coarse soluble) and three hydrophobic (ki: Aitken insoluble, ai: accumulation insoluble, ci: coarse insoluble) (Pringle et al., 2010). Each mode is treated as an individual tracer.

2.2.2 Partitioning between gas phase and particulate phase

Gas–particle partitioning is assumed to take place when SOC is in equilibrium between the gas and particulate phases. The concentration of the species that is bound to particles (C_particle) is calculated with

$$\begin{array}{}\text{(1)}& C_{\mathrm{particle}}=\mathit{\theta }\times (C_{\mathrm{particle}}+C_{\mathrm{gas}}),\end{array}$$

and the particulate mass fraction (θ) is defined as

$$\begin{array}{}\text{(2)}& \mathit{\theta }={\displaystyle \frac{C_{\mathrm{particle}}}{C_{\mathrm{particle}}+C_{\mathrm{gas}}}}={\displaystyle \frac{K_{\mathrm{p}}\times C_{\mathrm{PM}}}{\mathrm{1}+K_{\mathrm{p}}\times C_{\mathrm{PM}}}}={\displaystyle \frac{K_{\mathrm{p}}^{\prime }}{\mathrm{1}+K_{\mathrm{p}}^{\prime }}},\end{array}$$

where C_PM is the concentration of particulate matter or PM (µg m⁻³), K_p is the temperature-dependent particle–air partition coefficient (m³ µg⁻¹), and $K_{\mathrm{p}}^{\prime }$ is the dimensionless K_p.

In a model configuration using size-resolved particles (viz. the modal scheme), each SOC tracer is introduced in the model as eight different species: seven aerosol particles in ns, ks, as, cs, ki, ai, ci modes and one in the gas phase. The particulate fraction of the species in mode i (θ_i) is calculated using Eq. (2) with $C_{{\mathrm{PM}}_i}$ and $K_{{\mathrm{p}}_i}$ being the PM mass concentration and aerosol–air partition coefficient in the corresponding mode, respectively. The gaseous concentration C_gas is calculated using the sum of K_p values across modes, as well as total C_particle and C_PM:

$$\begin{array}{}\text{(3)}& C_{\mathrm{gas}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{7}}{C}_{{\mathrm{particle}}_i}/\sum _{i=\mathrm{1}}^{\mathrm{7}}{C}_{{\mathrm{PM}}_i}}{\sum _{i=\mathrm{1}}^{\mathrm{7}}{K}_{{\mathrm{p}}_i}}}.\end{array}$$

It is noted that this approach may not hold the constraint of mass consistency and is thus subject to further corrections. For the current study, the effects from this problem are expected to be minimal, given the fact that PAHs in the particulate phase are mainly distributed in the accumulation mode (Lammel et al., 2010, and references therein).

For K_p calculation four options of gas–particle partitioning schemes are available in SVOC: (1) a parameterization that is based on adsorption onto aerosol surface (Junge, 1977; Pankow, 1987), (2) absorption into organic matter (Finizio et al., 1997), (3) a combination of the two ways of organic matter absorption and black carbon adsorption (Lohmann and Lammel, 2004), and (4) multiple phases of the two-way sorption system (Goss and Schwarzenbach, 2001; Endo and Goss, 2014; Shahpoury et al., 2016). Two schemes used in this study are described below.

• Lohmann–Lammel scheme. The Lohmann–Lammel scheme takes into account an adsorption onto black carbon (BC) surface in addition to absorption into organic matter (OM) (Lohmann and Lammel, 2004). This dual-sorption theory empirically calculates K_p according to the following relation,

  $$\begin{array}{}\text{(4)}& \begin{array}{rl}{K}_{\mathrm{p}}& ={\mathrm{10}}^{-\mathrm{12}}\left(f_{\mathrm{OM}}{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{oct}}{\mathrm{MW}}_{\mathrm{oct}}}{{\mathit{\gamma }}_{\mathrm{OM}}{\mathrm{MW}}_{\mathrm{OM}}{\mathit{\rho }}_{\mathrm{oct}}}}{K}_{\mathrm{oa}}\right.\\ & \left.+f_{\mathrm{BC}}{\displaystyle \frac{a_{\text{atm-BC}}}{a_{\mathrm{soot}}{\mathit{\rho }}_{\mathrm{BC}}}}{K}_{\mathrm{sa}}\right),\end{array}\end{array}$$

  where ρ_BC is the density of BC (assumed as 1 kg L⁻¹), ρ_oct is the density of octanol (0.82 kg L⁻¹ at 20 ^∘C), K_sa is the partition coefficient between diesel soot and air, a_atm-BC is the available surface of atmospheric BC (m² g⁻¹), and a_soot is the specific surface area of diesel soot (m² g⁻¹). The adsorptive properties of diesel soot are selected to represent the atmospheric BC because this material is considered the most significant type of BC in polluted air.

  The K_sa value is calculated as a function of sub-cooled liquid vapor pressure $p_{\mathrm{L}}^{\mathrm{0}}$ using an estimate suggested by van Noort (2003),

  $$\begin{array}{}\text{(5)}& \log K_{\mathrm{sa}}=-\mathrm{0.85}\log p_{\mathrm{L}}^{\mathrm{0}}+\mathrm{8.94}-\log \left({\displaystyle \frac{\mathrm{998}}{a_{\mathrm{soot}}}}\right),\end{array}$$

  where a_soot in the model is set as 18.21 m² g⁻¹.

• Poly-parameter linear free-energy relationships (ppLFER) scheme. The concept of poly-parameter linear free-energy relationships (ppLFER) for the prediction of equilibrium partition coefficients is introduced by Goss and Schwarzenbach (2001), and its application in environmental chemistry has been reviewed by Endo and Goss (2014). This approach can describe a composite of different types of interactions between gas-phase species and aerosols. In contrast, single-parameter LFERs only correlate the partition coefficient to the sub-cooled liquid vapor pressure or the octanol–air partition coefficient of the species, and hence only valid within the group of compounds for which they were developed.

  In the study, the ppLFER scheme is incorporated into SVOC in which it defines K_p as the sum of individual partition coefficients representing surface adsorption and bulk-phase absorption processes to inorganic and organic aerosols. The formulation of K_p is adopted from Shahpoury et al. (2016) and is described as follows

  $$\begin{array}{}\text{(6)}& \begin{array}{rl}{K}_{\mathrm{p}}=& {\displaystyle \frac{K_{\mathrm{p}}^{\prime }}{C_{\mathrm{PM}}}},\\ K_{\mathrm{p}}^{\prime }=& K_{\mathrm{EC}}\times a_{\mathrm{EC}}\times C_{\mathrm{EC}}\times {\mathrm{10}}^{-\mathrm{6}}+\\ & K_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}\times a_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}\times C_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}\times {\mathrm{10}}^{-\mathrm{6}}+\\ & K_{\mathrm{NaCl}}\times a_{\mathrm{NaCl}}\times C_{\mathrm{NaCl}}\times {\mathrm{10}}^{-\mathrm{6}}+\\ & K_{\mathrm{DMSO}}\times {\displaystyle \frac{C_{\mathrm{WSOM}}}{{\mathit{\rho }}_{\mathrm{DMSO}}}}\times {\mathrm{10}}^{-\mathrm{6}}+\\ & K_{\mathrm{PU}}\times \mathrm{0.2}\times C_{\mathrm{WIOM}}\times {\mathrm{10}}^{-\mathrm{12}}+\\ & K_{\mathrm{hexadecane}}\times \mathrm{0.8}\times {\displaystyle \frac{C_{\mathrm{WIOM}}}{{\mathit{\rho }}_{\mathrm{hexadecane}}}}\times {\mathrm{10}}^{-\mathrm{6}},\end{array}\end{array}$$

  where K_EC, $K_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$, and K_NaCl are the substance partition (adsorption) coefficients (${\mathrm{m}}_{\mathrm{air}}^{\mathrm{3}}{\mathrm{m}}_{\mathrm{surface}}^{-\mathrm{2}}$) for elemental carbon or diesel soot, ammonium sulfate, and sodium chloride aerosol surface–air systems, respectively. K_DMSO is the substance partition (absorption) coefficient for the dimethyl sulfoxide–air system (${\mathrm{L}}_{\mathrm{air}}{\mathrm{L}}_{\mathrm{DMSO}}^{-\mathrm{1}}$). K_PU is the substance partition (absorption) coefficient for the polyurethane–air system (${\mathrm{m}}_{\mathrm{air}}^{\mathrm{3}}{\mathrm{kg}}_{\mathrm{PU}}^{-\mathrm{1}}$). K_hexadecane is the substance partition (absorption) coefficient for the hexadecane–air system (${\mathrm{L}}_{\mathrm{air}}{\mathrm{L}}_{\mathrm{hexadecane}}^{-\mathrm{1}}$); a_EC, $a_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$, and a_NaCl are the adsorbent-specific surface areas: 18.21, 0.1, and 0.1 ${\mathrm{m}}_{\mathrm{surface}}^{\mathrm{2}}{\mathrm{g}}_{\mathrm{adsorbent}}^{-\mathrm{1}}$, respectively. ρ_DMSO and ρ_hexadecane are the dimethyl sulfoxide and hexadecane densities of 1.1×10⁶ and 0.77×10⁶ g m⁻³, respectively. C_EC, $C_{({\mathrm{NH}}_{\mathrm{4}}{)}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$, C_NaCl, C_WSOM, C_WIOM are the concentration ($\mathrm{\mu }{\mathrm{g}}_{\mathrm{substance}}{\mathrm{m}}_{\mathrm{air}}^{-\mathrm{3}}$) of elemental carbon (here black carbon), ammonium sulfate, sodium chloride, water-soluble organic matter, and water-insoluble organic matter, respectively.

  The ppLFER scheme calculates the sorptive partition coefficient for every aerosol system, as summarized in Table S1 in the Supplement. Each coefficient requires information on system parameters (e, s, a, b, v, l), and the constant c, as shown in Table S2. The Abraham solute descriptors (E, S, A, B, V, and L) are substance specific; for the species selected in this study, refer to Table S6. All the predicted partition constants are adjusted to environmental temperature using the van 't Hoff equation

  $$\begin{array}{}\text{(7)}& \ln K_{\left(T\right)}=\ln K_{\left(T_{\mathrm{0}}\right)}-{\displaystyle \frac{\mathrm{\Delta }H}{\mathrm{R}}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{0}}}}\right),\end{array}$$

  where ΔH is the enthalpy of solvent–air phase transfer in joules per mole (J mol⁻¹). This variable is system specific and calculated by applying the ppLFER equations given in Table S3 and input parameters given in Table S4. The sequence of K_p calculation from ppLFER analysis in SVOC is illustrated in Fig. S1 in the Supplement.

2.2.3 Volatilization

1. Soil
   

   For soil volatilization, two parameterization schemes are implemented in the SVOC submodel: the Jury scheme (Jury et al., 1983, 1990) and the Smit scheme (Smit et al., 1997). The latter was applied in the study which is based on the volatilization of pesticides from the surface of fallow soils (Smit et al., 1997). The volatilization occurs upon partitioning over three soil phases (solid, gas, and liquid). The concentration of the chemical in the soil system (kg m⁻³) is formulated as

   $$\begin{array}{}\text{(8)}& C_{\mathrm{soil}}=Q\times C_{\mathrm{vapor}},\end{array}$$

   and the capacity factor Q is given by

   $$\begin{array}{}\text{(9)}& Q=\mathit{\psi }+\mathit{\phi }{K}_{\mathrm{wa}}+{\mathit{\rho }}_{\mathrm{soil}}{K}_{\mathrm{wa}}{K}_{\mathrm{sl}},\end{array}$$

   where ψ and φ are the volume fractions of air and moisture, respectively; ρ_soil is the soil density; K_wa is the water–air partition coefficient where $K_{\mathrm{wa}}=\mathrm{1}/K_{\mathrm{aw}}$; and K_sl is the solid–liquid partition coefficient. K_aw is calculated based on the Henry's law constant:

   $$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{aw}}=\mathrm{1}/\left(H\mathrm{R}T\right)\mathrm{and}\text{(11)}& {\displaystyle }& {\displaystyle }H=k_{\mathrm{H}}^{\ominus }\times \exp \left[{\displaystyle \frac{-\mathrm{\Delta }{H}_{\mathrm{soln}}}{\mathrm{R}}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{0}}}}\right)\right],\end{array}$$

   where H is the temperature-adjusted Henry coefficient (M atm⁻¹), R is the dry air gas constant (8.314 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), T is the environment temperature (K), ΔH_soln is the enthalpy of dissolution (J mol⁻¹).

   K_sl can be set equal to the sorption coefficient to soil organic matter K_om times the fraction of organic carbon in soil $f_{{\mathrm{OM}}_{\mathrm{s}}}$. Since K_om data are not available, the coefficient for sorption to soil organic carbon K_oc was used to estimate K_sl:

   $$\begin{array}{}\text{(12)}& K_{\mathrm{sl}}=\mathrm{0.56}{K}_{\mathrm{oc}}{f}_{{\mathrm{OM}}_{\mathrm{s}}}.\end{array}$$

   Mackay and Boethling (2000) has suggested a reasonably good regression relationship between K_oc and octanol–water partition coefficient K_ow for PAHs:

   $$\begin{array}{}\text{(13)}& \log (K_{\mathrm{oc}}/\mathrm{1000})=\mathrm{0.823}\log (K_{\mathrm{ow}}/\mathrm{1000})-\mathrm{0.727}.\end{array}$$

   where the factor of 1000 is needed because K_oc and K_ow are expressed in cubic meter per kilogram (m³ kg⁻¹), whereas in the original regression they used milliliter per gram (mL g⁻¹).

   Once Q is computed, the dimensionless fraction of the chemical in the gas phase F_gas is then calculated as

   $$\begin{array}{}\text{(14)}& F_{\mathrm{gas}}={\displaystyle \frac{\mathit{\psi }}{Q}}.\end{array}$$

   In the Smit scheme, an empirical relation was established between F_gas and cumulative volatilization (CV in percent of substance deposit). CV was determined based on field and greenhouse experiments with numerous pesticides at 21 d after application. For normal to moist field conditions, CV is expressed as

   $$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathrm{CV}=& \mathrm{71.9}+\mathrm{11.6}\log \left(\mathrm{100}{F}_{\mathrm{gas}}\right);\\ & \mathrm{6.33}\times {\mathrm{10}}^{-\mathrm{9}}<F_{\mathrm{gas}}\le \mathrm{1},\end{array}\end{array}$$

   and for dry field conditions

   $$\begin{array}{}\text{(16)}& \begin{array}{rl}\mathrm{CV}=& \mathrm{42.3}+\mathrm{9.0}\log \left(\mathrm{100}{F}_{\mathrm{gas}}\right);\\ & \mathrm{0.2}\times {\mathrm{10}}^{-\mathrm{6}}<F_{\mathrm{gas}}\le \mathrm{1}.\end{array}\end{array}$$

2. Vegetation
   

   Smit et al. (1998) derived an equation for the cumulative volatilization, CV, from plants against vapor pressure P_v (mPa) at 7 d after application based on field and climate chamber experiments of pesticide volatilization (Eq. 17).

   $$\begin{array}{}\text{(17)}& \mathrm{CV}={\mathrm{10}}^{\mathrm{1.528}+\mathrm{0.466}\log P_{\mathrm{v}}};P_{\mathrm{v}}\le \mathrm{10.3}.\end{array}$$

   For compounds with P_v above 10.3 mPa, CV is set at 100 % of deposit. Temperature adjustments were made for P_v using the Clausius–Clapeyron equation:

   $$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}\left(\ln P_{\mathrm{v}}\right)}{\mathrm{d}T}}=-{\displaystyle \frac{\mathrm{\Delta }{H}_{\mathrm{vap}}}{\mathrm{R}{T}^{\mathrm{2}}}}.\end{array}$$

3. Snow and glaciers
   

   The parameterization of substance loss by volatilization from snow pack follows Wania (1997), whereby the process is calculated using a consecutive cycle of an equilibrium partitioning among four phases followed by a contaminant loss. The four phases considered are liquid water, organic matter contained in the snowpack, snow pores (air), and an ice–air interface. Fugacity capacity factors for these phases are expressed with the following relations:

   $$\begin{array}{}\text{(19a)}& {\displaystyle }\text{air}\left(\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{Pa}}^{-\mathrm{1}}\right),Z_{\mathrm{a}}=\mathrm{1}/\mathrm{R}T;\text{(19b)}& {\displaystyle }\begin{array}{rl}& \mathrm{water}\left(\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{Pa}}^{-\mathrm{1}}\right),\\ & Z_{\mathrm{l}}=K_{\mathrm{wa}}/\mathrm{R}T=K_{\mathrm{wa}}{Z}_{\mathrm{a}}\end{array};\text{(19c)}& {\displaystyle }\begin{array}{rl}& \text{organic carbon}\left(\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{Pa}}^{-\mathrm{1}}\right),\\ & Z_{\mathrm{o}}=Z_{\mathrm{l}}\mathrm{0.41}{K}_{\mathrm{ow}}^{\prime };\end{array}\text{(19d)}& {\displaystyle }\begin{array}{rl}& \text{ice–air interface}\left(\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{Pa}}^{-\mathrm{1}}\right),\\ & z_i=K_{\mathrm{ia}}/\mathrm{R}T=K_{\mathrm{ia}}{Z}_{\mathrm{a}};\end{array}\end{array}$$

   where R is the dry air gas constant (8.312 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), T is the air temperature (K), K_wa is the water–air partition coefficient (unitless), $K_{\mathrm{ow}}^{\prime }$ is the dimensionless octanol–water partition coefficient, and K_ia is the ice surface–air partition coefficient (m). K_ia at 20 ^∘C is estimated using K_wa and water solubility $C_{\mathrm{w}}^{\mathrm{s}}$ (mol m⁻³),

   $$\begin{array}{}\text{(20)}& \log K_{\mathrm{ia}}\left(\mathrm{20}{}^{\circ }\mathrm{C}\right)=-\mathrm{0.769}\log C_{\mathrm{w}}^{\mathrm{s}}-\mathrm{5.97}+\log K_{\mathrm{wa}},\end{array}$$

   and further extrapolated to other temperatures using enthalpy of condensation of solid (ΔH_subl in J mol⁻¹),

   $$\begin{array}{}\text{(21)}& \begin{array}{rl}\log K_{\mathrm{ia}}\left(T\right)=& \log K_{\mathrm{ia}}\left(\mathrm{20}{}^{\circ }\mathrm{C}\right)\\ & +{\displaystyle \frac{\mathrm{0.878}\mathrm{\Delta }{H}_{\mathrm{subl}}}{\mathrm{2.303}\mathrm{R}}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{293}}}\right).\end{array}\end{array}$$

   An equilibrium fugacity f_s is thereby determined by

   $$\begin{array}{}\text{(22)}& f_{\mathrm{s}}={\displaystyle \frac{M_{\mathrm{sp}}}{Z_{\mathrm{a}}{v}_{\mathrm{a}}+Z_{\mathrm{l}}{v}_{\mathrm{l}}+z_i{A}_{\mathrm{snow}}{\mathit{\rho }}_{\mathrm{mw}}+Z_{\mathrm{o}}{v}_{\mathrm{o}}}},\end{array}$$

   where M_sp is the amount of chemical contained in snowpack (the model here applies snow burden of the chemical in kg m⁻²) and v_a, v_l, and v_o are the volume fractions of air, liquid water, and organic matter in snowpack (m³ m⁻³), respectively. For this study, v_a and v_l values are set to 0.3 and 0.1, respectively, whereas v_o is zero assuming no polluted snow. A_snow is the specific snow surface area (m² g⁻¹). In Daly and Wania (2004), a value of 0.1 m² g⁻¹ for A_snow was used for snow accumulation period and a linear decrease from 0.1 to 0.01 m² g⁻¹ was used during the snowmelt period. In SVOC submodel, a value of 0.025 m² g⁻¹ is adopted for A_snow to represent a fairly aged snowpack. ρ_mw is the density of snowmelt water and here is taken as 7×10⁵ g m⁻³.

   Volatilization rate ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is calculated by applying

   $$\begin{array}{}\text{(23)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{sp}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{1}}{U_{\mathrm{7}}{Z}_{\mathrm{a}}}+\frac{\mathrm{1}}{U_{\mathrm{5}}{Z}_{\mathrm{l}}+U_{\mathrm{6}}{Z}_{\mathrm{a}}}}}\times {\displaystyle \frac{f_{\mathrm{s}}}{h_{\mathrm{s}}}},\end{array}$$

   where h_s is the snow depth (m), U₅ is the snow–water phase diffusion mass transfer coefficient (m h⁻¹), U₆ is the snow–air phase diffusion mass transfer coefficient (m h⁻¹), and U₇ is the snow–air boundary layer mass transfer coefficient (m h⁻¹). U₅ and U₆ are calculated from molecular diffusivities in air and water (Eqs. 24 and 25), whereas a typical value of 5 m h⁻¹ is adopted for U₇.

   $$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{5}}=B_{\mathrm{w}}{\displaystyle \frac{v_{\mathrm{l}}^{\mathrm{10}/\mathrm{3}}/(v_{\mathrm{a}}+v_{\mathrm{l}}{)}^{\mathrm{2}}}{\ln \mathrm{2}{h}_{\mathrm{s}}}},\text{(25)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{6}}=B_{\mathrm{a}}{\displaystyle \frac{v_{\mathrm{a}}^{\mathrm{10}/\mathrm{3}}/(v_{\mathrm{a}}+v_{\mathrm{l}}{)}^{\mathrm{2}}}{\ln \mathrm{2}{h}_{\mathrm{s}}}},\end{array}$$

   where B_w and B_a are the molecular diffusivities (m² h⁻¹) in water and air, respectively. In the model, B_a is derived from the molecular weight (MW) as $B_{\mathrm{a}}=\frac{\mathrm{1.55}}{{\mathrm{MW}}^{\mathrm{0.65}}}$ cm² s⁻¹, whereas B_w is set as 1×10⁴ less than B_a, following Schwarzenbach et al. (2005).

4. Ocean
   

   In the study, the ocean is represented as a surface mixed layer of a depth varying spatially and in time without lateral transports. The mixed layer depths were obtained from (de Boyer Montégut et al., 2004). The SOC volatilization flux from the sea surface is parameterized based on the two-film model of Liss and Slater (1974) and is calculated within the AIRSEA submodel (Pozzer et al., 2006). Note that no ocean and sea-ice dynamics were included in the simulations.

   Sorption of SOCs in water to suspended particulate matter (colloidal or sinking detritus) is neglected. Therefore, SOC concentration in surface seawater, and hence volatilization from sea surface, is overestimated, in particular for very lipophilic (log iK_ow>6) substances. This bias is negligible for the substances studied here (PAHs), which are less lipophilic or volatilization is limited by vapor pressure (e.g., benzo(a)pyrene). Forces from strong winds and dissolved or particulate organics in seawater are transferred to air via sea spray, which adds to particulate OM in air over the ocean (O'Dowd et al., 2008; Qureshi et al., 2009). This process is neglected in the model.

2.2.4 Dry deposition

Dry deposition is simulated using deposition velocities. For gas-phase SOCs, the velocities are calculated by the DDEP submodel (Kerkweg et al., 2006a), whereas particulate-bound SOCs are assumed to deposit at similar rates to other aerosols whose velocities are also computed by DDEP. If the modal scheme is selected (see Sect. 2.2.1), the particle deposition velocity $v_{\mathrm{d}}^{\mathrm{SOC}}$ at mode i is equal to the aerosol deposition velocity $v_{\mathrm{d}}^{\mathrm{aer}}$ at the respective mode. On the other hand, for the bulk scheme, $v_{\mathrm{d}}^{\mathrm{SOC}}$, is computed as a weighted average of $v_{\mathrm{d}}^{\mathrm{aer}}$ from the four BC modes (ki, ks, as, and cs) where the weight is the surface area of BC. This approach is most relevant for PAHs as they are assumed to be predominantly transported by sorption to BC. The above relations are formulated as

$$\begin{array}{}\text{(26a)}& {\displaystyle }& {\displaystyle }\text{modal scheme: }& {\displaystyle }{v}_{\mathrm{d},i}^{\mathrm{SOC}}=& {\displaystyle }{v}_{\mathrm{d},i}^{\mathrm{aer}};& {\displaystyle }& {\displaystyle }\text{(26b)}& {\displaystyle }& {\displaystyle }\text{bulk scheme: }& {\displaystyle }{v}_{\mathrm{d},\mathrm{bulk}}^{\mathrm{SOC}}=& {\displaystyle }{\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{4}}{S}_{{\mathrm{BC}}_i}\times v_{\mathrm{d},i}^{\mathrm{aer}}}{\sum _{i=\mathrm{1}}^{\mathrm{4}}{S}_{{\mathrm{BC}}_i}}};& {\displaystyle }& {\displaystyle }\end{array}$$

and the BC surface area per unit volume, S_BC (cm² cm⁻³), is given by

$$\begin{array}{}\text{(27)}& S_{{\mathrm{BC}}_i}=\mathrm{4}\mathit{\pi }{\left[r_i\exp \left({\ln }^{\mathrm{2}}{\mathit{\sigma }}_{{\mathrm{g}}_i}\right)\right]}^{\mathrm{2}}{N}_i\times {\displaystyle \frac{C_{{\mathrm{BC}}_i}}{C_{{\mathrm{aer}}_i}}},\end{array}$$

where N_i is the number concentration for mode i (cm⁻³), r_i is the number radius (cm), ${\mathit{\sigma }}_{{\mathrm{g}}_i}$ is the geometric standard deviation, $C_{{\mathrm{BC}}_i}$ is the BC concentration (µg m⁻³) in mode i, and $C_{{\mathrm{aer}}_i}$ is the sum of aerosol concentrations in the same mode (µg m⁻³).

2.2.5 Wet deposition

Wet deposition is applied to both gas and particulate SOCs. The gaseous fraction is scavenged into cloud and rain droplets according to diffusion limitation, Henry's law equilibrium, and accommodation coefficient, and this process is parameterized and solved empirically in the SCAV submodel (Tost et al., 2006a). Particulate-phase SOCs are scavenged in convective updrafts, rainout and washout, and cloud evaporation, with the rate being proportional to BC wet scavenging; hence the change in SOC concentration is described as

$$\begin{array}{}\text{(28a)}& {\displaystyle }& {\displaystyle }\text{modal scheme: }{\displaystyle \frac{\mathrm{\Delta }{C}_{{\mathrm{SOC}}_i}}{\mathrm{\Delta }t}}={\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{SOC}}_i}}{{\mathit{\mu }}_{{\mathrm{BC}}_i}}}\times {\displaystyle \frac{\mathrm{\Delta }{C}_{{\mathrm{BC}}_i}}{\mathrm{\Delta }t}};\text{(28b)}& {\displaystyle }& {\displaystyle }\text{bulk scheme: }{\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{SOC}}}{\mathrm{\Delta }t}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{SOC}}}{\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\mu }}_{{\mathrm{BC}}_i}}}\times \sum _{i=\mathrm{1}}^{\mathrm{4}}{\displaystyle \frac{\mathrm{\Delta }{C}_{{\mathrm{BC}}_i}}{\mathrm{\Delta }t}};\end{array}$$

where µ is the particle volume mixing ratio (${\mathrm{mol}}_{\mathrm{SOC}/\mathrm{BC}}{\mathrm{mol}}_{\mathrm{air}}^{-\mathrm{1}}$) and Δt is the model time step (s). Note that Eq. (28a) imposes a restrictive prerequisite, namely BC and the particle-bound SOC have similar size distributions. When this condition is not met, there is a high possibility of an artificial mass being produced, usually to the largest aerosol mode. To solve this problem, a correction factor is applied and defined as a function of the ratio of positive fluxes to negative fluxes integrated across levels and modes.

2.2.6 Atmospheric degradation

The atmospheric degradation of SOCs in the gas phase as well as within aerosol particles are explicitly treated in SVOC. The gas-phase chemical mechanism is calculated within the MECCA submodel (Sander et al., 2011). SOC gaseous degradation is from photochemical reactions with OH, NO₃, and O₃ radicals, which follow a second-order transformation, with the rate constants k⁽²⁾ obtained from laboratory studies. The $k_{\mathrm{OH}}^{\left(\mathrm{2}\right)}$ value is typically higher, suggesting that oxidation with OH radical is the dominant loss pathway.

Most models do not consider oxidation rate of particulate-phase SOCs as experimental aerosols studied in laboratory cover only a small part of relevant atmospheric aerosols. For PAHs, such as benzo(a)pyrene, which stays mostly in the particulate phase, the degradation is more efficient by surface reactions with O₃ (Shiraiwa et al., 2009, and references therein) with the rate depending on the substrate. The SVOC submodel includes the degradation process of PAHs on aerosol particles from the O₃ reaction with one assumption; that is, the heterogeneous reaction does not lead to a change in the oxidant concentration. Due to a limited number of kinetic studies of heterogeneous reactions of PAHs, only two species are considered (phenanthrene and benzo(a)pyrene). Nevertheless, the submodel structure provides a relatively straightforward approach to allow more species in the future.

The reaction rate coefficient for particulate-phase phenanthrene with O₃ at aerosol surfaces was derived from laboratory experiments using chemically unspecific model aerosol (silica) with PAH surface coverage of less than a monolayer (Perraudin et al., 2007). To this end, the second-order rate coefficient, k⁽²⁾, in ${\mathrm{cm}}^{\mathrm{4}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ was derived from the reported phenanthrene (PHE) decay kinetics, $k_{{\mathrm{O}}_{\mathrm{3}},\mathrm{het}}^{\left(\mathrm{2}\right)}$ (${\mathrm{cm}}^{\mathrm{3}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$), as $k^{\left(\mathrm{2}\right)}=k_{{\mathrm{O}}_{\mathrm{3}},\mathrm{het}}^{\left(\mathrm{2}\right)}/\left(\frac{S}{V}\right)=\left(\mathrm{6.2}\pm \mathrm{4.8}\right)\times {\mathrm{10}}^{-\mathrm{17}}$ ${\mathrm{cm}}^{\mathrm{4}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, with $\frac{S}{V}$ (cm⁻¹) being the experimental aerosol surface concentration (0.56±0.43 cm⁻¹ in Perraudin et al., 2007). In the submodel, $k_{{\mathrm{O}}_{\mathrm{3}},\mathrm{het}}^{\left(\mathrm{2}\right)}$ (${\mathrm{cm}}^{\mathrm{3}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) is calculated using the ambient aerosol surface concentration. As for benzo(a)pyrene, the pseudo-first-order rate coefficient, $k_{{\mathrm{O}}_{\mathrm{3}},\mathrm{het}}^{\left(\mathrm{1}\right)}$ in reciprocal seconds (s⁻¹), was derived from surface-adsorbed benzo(a)pyrene reaction with O₃ on solid organic and salt aerosols following the Langmuir–Hinshelwood mechanism (Kwamena et al., 2004).

2.2.7 Biotic and abiotic degradations

Biotic and abiotic processes in surface compartments contribute to the degradation of chemicals and are strongly dependent on local environmental conditions, e.g., nutrient contents, water, temperature, PH, and light. In SVOC, these factors are not explicitly quantified. The degradation is alternatively described as following a first-order decay law (Eq. 29). The 10 K temperature warming is assumed to double the rate of degradation (Eq. 30), following recommendations in chemical risk assessment (European Commission, 2000) and consistent with findings, such as a two-time increase in the growth of hydrocarbon-degrading microbes found in soils (Thibault and Elliott, 1979).

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle \frac{\partial C_{{\mathrm{SOC}}_{\mathrm{s}}}}{\partial t}}=-k_{\mathrm{sfc}}\times C_{{\mathrm{SOC}}_{\mathrm{s}}},\text{(30)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{sfc}\left(T\right)}=k_{\mathrm{sfc}\left(T_{\mathrm{ref}}\right)}\times {\mathrm{2}}^{\frac{T-T_{\mathrm{ref}}}{\mathrm{10}}},\end{array}$$

where $C_{{\mathrm{SOC}}_{\mathrm{s}}}$ is the substance concentration (kg m⁻³) in surface compartments (i.e., soil, vegetation, or ocean) and k_sfc is the first-order decay rate (s⁻¹). T_ref is the reference temperature, i.e., 298 K for soil and 273 K for ocean. Note that the degradation in vegetation is calculated assuming the same k_sfc for the soil compartment.

2.3 Input data

2.3.1 Kinetic and physicochemical properties

The model simulations were performed for four PAH species: phenanthrene (PHE), pyrene (PYR), fluoranthene (FLT), and benzo(a)pyrene (BaP). To simulate the fate and environmental distribution of these species, the model requires some physicochemical properties as summarized in Table S5 of the Supplement. These include equilibrium partition coefficients and their related energies of phase transfer. The characteristics from PHE to BaP are indicated by decreasing volatility (as molar mass increases), increasing K_oa and K_ow, and decreasing water solubility (as $C_{\mathrm{w}}^{\mathrm{s}}$ and Henry's coefficients decrease). The properties also include the second-order rate coefficients for homogeneous oxidation with OH, O₃, and NO₃ except for BaP where the gaseous reaction is switched off. Heterogeneous oxidation by O₃ is simulated only for PHE and BaP. Furthermore, the model also requires compound solute parameters for simulations using the ppLFER gas–particle partitioning scheme (Table S6).

2.3.2 Emissions and other model input

As model input, several emission datasets were employed in the study. Emission estimates for PAHs were obtained from the annual mean inventory of Shen et al. (2013) for the year 2008. They applied regression and technology split methods to construct country-level emissions for six categories (coal, petroleum, natural gas, solid wastes, biomass, and an industrial process category) or six sectors (energy production, industry, transportation, commercial/residential sources, agriculture, and deforestation/wildfire) before further regridding the emissions to a ${\mathrm{0.1}}^{\circ }\times {\mathrm{0.1}}^{\circ }$ grid.

Emissions of aerosol species such as organic carbon (OC), black carbon (BC), mineral dust (DU), and sea salt (SS) were included. For BC and OC, the Representative Concentration Pathway (RCP) 6.0 emission scenario of the IPCC (Intergovernmental Panel on Climate Change) (van Vuuren et al., 2011) was used and accessible via ftp://ftp-ipcc.fz-juelich.de/pub/emissions/gridded_netcdf (last access: 15 February 2018). Emissions are calculated for anthropogenic, biomass burning, ship, and aircraft. The RCP database provides a seasonality only for the biomass burning and ship emissions. In this study, seasonal-scale factors were applied to the anthropogenic emission whereby the seasonality was based upon the monthly variation of the Hemispheric Transport of Air Pollutants (HTAP) v2.2 anthropogenic emission inventory (Janssens-Maenhout et al., 2015). BC emissions from all sectors were assumed to be hydrophobic. For OC, it was assumed to be 65 % hydrophilic and 35 % hydrophobic upon biomass burning emissions and to be 100 % hydrophobic upon anthropogenic and ship emissions. Both OC and BC were emitted at Aitken mode, which spans the size range from about 5 to 50 nm in diameter. A factor of 1.724 was used to scale the OC emissions to primary organic matter (OM). It is noteworthy that the formation of secondary organic aerosols (SOAs) from atmospheric oxidation and condensation of volatile organic compounds (VOCs) were not treated in the simulations. In the model, particulate organic matter is emitted and transported as a bulk aerosol species (OM).

DU and SS emissions were computed online by the ONEMIS submodel (Kerkweg et al., 2006b). DU emission flux is calculated based on wind speed at 10 m altitude and soil parameters (Schulz et al., 1998). The emission of SS particles by bubble bursting is described as wind-speed-dependent particle mass and number fluxes at accumulation (50–500 nm) and coarse (>500 nm) modes. In ONEMIS, the fluxes are determined from precalculated lookup tables following the Guelle et al. (2001) parameterization. SVOC submodel accounts for OM fraction in the SS mass fluxes (J_SS), and the fraction is estimated using a 10 m wind (v₁₀)-dependent empirical relationship derived from Fig. 2a in Gantt et al. (2011). Equation (31) below is used to calculate the OM mass fluxes, J_OM, in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$.

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{J}_{\mathrm{OM}}=& J_{\mathrm{SS}}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{0.78}}{\mathrm{1}+\mathrm{0.03}\exp \left(\mathrm{0.48}{v}_{\mathrm{10}}\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{0.24}}{\mathrm{1}+\mathrm{0.05}\exp \left(\mathrm{0.38}{v}_{\mathrm{10}}\right)}}\right).\end{array}\end{array}$$

Emissions of other gases including volatile organic species (SO₂, CO, NH₃, NO, CH₄, and non-methane hydrocarbon (NMHC)) were prescribed using the IPCC RCP6.0 dataset (van Vuuren et al., 2011). Global estimates for the soil properties dry bulk density and organic matter fraction were obtained from Dunne and Willmott (1996) and Batjes (1996), respectively.

Figure 2 Locations of monitoring stations used in the study. The initial letter of each station ID refers to the individual monitoring network (E: EMEP and AMAP, D: DEFRA, I: IADN, M: MONET-Africa)

2.4 Observational data

The observation data used for model performance evaluation were collected from several surface monitoring networks: the European Monitoring and Evaluation Programme (EMEP) (Tørseth et al., 2012), the Arctic Monitoring and Assessment Program (AMAP) (Hung et al., 2005), the Great Lakes Integrated Atmospheric Deposition Network (IADN) (IADN, 2014), the Department for Environment, Food and Rural Affairs (DEFRA) UK-AIR (Air Information Resource) program (DEFRA, 2010), and the MONitoring NETwork in the African continent (MONET-Africa) (Klánová et al., 2008). These data were screened and quality controlled according to the description in the Supplement Sect. SIII. Final stations with reliable monthly data are depicted in Fig. 2, with 3 stations in the Arctic, 19 in the northern midlatitudes, and 6 in the tropics. The availability of data differs by station, species, and variable of interest; see the site-specific information in Table S10 for total concentration and Table S11 for particulate mass fraction (θ).

The study also compared simulated concentrations in the marine atmosphere to two ship cruise measurement campaigns: (1) on a west to east transect across the tropical Atlantic Ocean (Lohmann et al., 2013) and (2) along the Asian marginal seas, the Indian Ocean, and the Pacific Ocean (Liu et al., 2014). The monthly mean modeled values were compared to daily measurements at each sampling point.

2.5 Experiment design

2.5.1 Model configuration

The model was run on a spectral T42 grid in the horizontal (approximately 2.8^∘ in a lat–long grid) and 19 unevenly distributed layers in the vertical with the top level at 10 hPa. The vertical layers are discretized using a hybrid coordinate (the lowest level follows the terrain and become a surface of constant pressure in the stratosphere). All simulations were run for a 3-year period (i.e., 2007–2009), with a 1-year spinup (i.e., 2006), and nudged toward the European Centre for Medium-range Weather Forecasts (ECMWF) reanalysis data (Dee et al., 2011). Note that the simulation period was selected based on the representative year of PAH emissions (i.e., 2008) and the availability of reliable observation data (see the Supplement Sect. SIII).

2.5.2 Sensitivity to the temporal resolution of emissions and process parameterizations

Factor separation analysis (Stein and Alpert, 1993) was used to quantitatively evaluate the contributions to changes in a particular output variable that result from changing components of model input and physics parameterizations. The model sensitivity to four model components (hereinafter, “factors”) was tested. The four factors were the following:

1. Temporal resolution of emissions (hereinafter, fac1). The PAH emission inventory of Shen et al. (2013) was based on 2008 annual emission totals from all sectors (see Sect. 2.3.2). The emissions were divided over the year using monthly factors derived from BC anthropogenic emissions. Two sets of simulations were carried out to test the sensitivity of model output to the seasonal profile of emission. The first set used constant emissions throughout the simulations, whereas the second set used a monthly emission interval.

2. The size-discretization of particulate-phase PAHs (hereinafter, fac2). The two options for this factor were tested: bulk versus modal (Sect. 2.2.1). Note that with regards to BaP, 95 % of the emissions were assumed to be in particulate phase and for the modal-scheme scenario, all of the emitted particles are treated as the hydrophobic Aitken (ki) tracers.

3. The choice of gas–particle partitioning scheme (hereinafter, fac3). The present study focuses on the comparison between the Lohmann–Lammel and ppLFER schemes for gas–particle partitioning.

4. The influence of re-volatilization (hereinafter, fac4). Model runs with the volatilization process switched off are compared to those runs which have volatilization switched on.

The factor separation technique is described in the Supplement Sect. SIV including the equations used to compute the model sensitivity to four factors. A total of 16 (or 2ⁿ) experiments summarized in Fig. 3 are necessary to supply the complete solution for the factor analysis. The ABLN (annual + bulk + Lohmann–Lammel + no volatilization) experiment is designed to be the base simulation (f₀), in which annual emission (A), the bulk scheme (B), the Lohmann–Lammel scheme (L), and no re-volatilization (N) were applied. SMPW (i.e., Seasonal emission + Modal scheme + PpLFER scheme + With re-volatilization) is referred to as the target simulation (f₁₂₃₄) in which the more sophisticated choice of the four features (factors) was tested. The total (gas + particle) concentration at the lowest model level was selected for its higher relevance with all the factors (compared to, for example, atmospheric burden) and to facilitate direct comparison with observations.

Figure 3 List of experiments performed for the factor separation analysis to study sensitivity to temporal variation in emission and process parameterizations (particulate-phase representation, gas–particle partitioning scheme, and volatilization); L+L: Lohmann–Lammel; PpLFER: poly-parameter linear free-energy relationships.

3 Results and discussion

3.1 Sensitivity tests

The analysis of the factor separation results is given below. For each factor, the analysis includes the assessment of direct effects (${\widehat{f}}_i$) and total interaction effects ($\mathrm{\Sigma }{\widehat{f}}_{ij}+\mathrm{\Sigma }{\widehat{f}}_{ijk}+{\widehat{f}}_{\mathrm{1234}}$) on near-surface PAH concentrations in two seasons, i.e., December–January–February (DJF) and June–July–August (JJA). Figure  shows the respective effects for all factors as relative to the seasonal means of the base experiment (f₀). A positive value indicates a concentration increase with respect to f₀, whereas a negative indicates a decrease. The spatial distributions of f₀ and f₁₂₃₄ seasonal mean concentrations for the four species are shown in the Supplement Sect. SVI Figs. S4 and S5, respectively.

Figure 4 Direct and interaction effects on seasonal-mean near-surface PAH concentrations of (a, b) monthly emissions (i=1), (c, d) the modal scheme (i=2), (e, f) the ppLFER scheme (i=3), and (g, h) volatilization (i=4). The direct effects (a, c, e, g) are expressed as the difference between two distributions (${\widehat{f}}_i=f_i-f_{\mathrm{0}}$), whereas the interaction effects (b, d, f, h) are expressed as the sum of two ($\mathrm{\Sigma }{\widehat{f}}_{ij}$, i≠j), three ($\mathrm{\Sigma }{\widehat{f}}_{ijk}$, $i\ne j\ne k$), and all (${\widehat{f}}_{\mathrm{1234}}$) factor interactions. They are presented as relative to concentrations from the base (f₀) simulation. The figures display the median, 25th and 75th percentiles of the relative effects over each of five main climatic regions. Note that an inverse hyperbolic sine function has been used in scaling the y axes.

We studied the relative effects in five climate zones (Arctic, northern midlatitudes, the tropics, southern midlatitudes, and Antarctica) The global distributions of the relative effects are presented in Figs. S6–S13, whereas Figs. S14–S21 present the relative interaction effects from the individual combination of factors. In the following, we do not look to interpret concentration responses to each interaction term. The reasons for this are that (1) accounting for all such interactions is complicated given the number of factors and (2) higher-order interactions (combinations of more than two factors) are hard to physically interpret.

We further investigate the factor effects on model performance by comparing the predicted seasonal mean near-surface concentrations from 16 experiments against observation data in the Arctic and northern midlatitudes (Supplement Sect. SVII).

3.1.1 Effects of seasonality of emissions

Figure a shows that using monthly emissions increases PAH concentrations in DJF and decreases the concentrations in JJA over the areas from the middle to high latitudes of the Northern Hemisphere (NH). This result is expected and is attributed to emissions during the northern winter (summer) being higher (lower) than annual means and photochemistry being less (more) active. Over the Arctic, the relative changes (${\widehat{f}}_{\mathrm{1}}/f_{\mathrm{0}}$) in DJF show a median increase of approx. 30 % for PHE, PYR, and FLT, and 7 % for BaP, whereas ${\widehat{f}}_{\mathrm{1}}/f_{\mathrm{0}}$ in JJA is weaker in magnitude for PHE and PYR (−16 %) but comparable for FLT (−28 %) and BaP (−5 %). Accounting for seasonality leads to generally lower bias for PHE, and this effect is more pronounced in the middle latitudes rather than in high latitudes (Supplement Sect. SVII, Fig. S22).

In general, ${\widehat{f}}_{\mathrm{1}}/f_{\mathrm{0}}$ becomes smaller over the northern midlatitudes by around half. The upper (lower) quartile of ${\widehat{f}}_{\mathrm{1}}/f_{\mathrm{0}}$ in DJF (JJA) indicates that about one-quarter of the areas of the temperate and polar regions experience an at least 40 % increase (decrease), most were located in northeastern Eurasia (see the left panels of Figs. S6 and S7). Note ${\widehat{f}}_{\mathrm{1}}/f_{\mathrm{0}}$ over the tropics are small to negligible (±1 %) mainly due to little variation in emissions from anthropogenic sectors. PAH concentrations may be higher in dry season due to increased amounts of biomass burning, but they are poorly represented in the current inventory. In southern middle and high latitudes, the direct effects of emission change are substantially opposite in sign to the effects seen in the northern latitudes, being negative in DJF (median ranges from −4 % to −32 %) and positive in JJA (7 %–25 %).

The total interactions between fac1 and other factors generally produce opposite signals to ${\widehat{f}}_{\mathrm{1}}$ over middle and high latitudes in the two seasons. This result indicates that the changes in other factors tend to buffer the influence of monthly emission on increasing or decreasing f₀ concentrations. Some exceptions are seen over parts of East Asia in DJF for all species (Fig. S6, right panels) and over the Southern Ocean in JJA for BaP (Fig. S7, right panels) where the interactions work to reinforce the direct effects. In DJF, the degree of interactions is smaller or comparable to the size of ${\widehat{f}}_{\mathrm{1}}$ for the Arctic and northern midlatitudes but becomes stronger by at least double for Antarctica. The opposite tendency is seen in JJA but only applies to PYR and FLT. In agreement to ${\widehat{f}}_{\mathrm{1}}$, the interaction effects are less apparent over the tropics. Note that the positive effects in ${\widehat{f}}_{\mathrm{14}}$ during local summer tend to be more dominant than the effects in other combinations for PHE, PYR, and FLT (Figs. S18–S20). In the simulation, the presence of re-volatilization in summer tends to suppress ${\widehat{f}}_{\mathrm{1}}$ by promoting more gases available for long-range transport, and thus implies a negative feedback.

3.1.2 Effects of size discretization of particulate-phase tracer

The direct effects of the modal scheme (${\widehat{f}}_{\mathrm{2}}$) vary among species (Fig. c). ${\widehat{f}}_{\mathrm{2}}$ is almost absent for PHE as the species resides almost completely in the gas phase. For PYR and FLT, ${\widehat{f}}_{\mathrm{2}}$ is negative during DJF over northern midlatitudes (${\widehat{f}}_{\mathrm{2}}/f_{\mathrm{0}}$ quartiles range from −5 % to −30 %) and the Arctic (−50 % to −75 %), whereas it is hardly visible in JJA or over other regions. Further analysis reveals stronger particle deposition results when the aerosol phase is discretized into different modes (not shown). In long-range transport under modal aerosol representation, the aerosols are more associated with larger particles, hence particle deposition becomes more effective. The choice of size discretization has only minor effects for atmospheric levels, except for BaP, especially during DJF, for which overestimates are significantly compensated for (Fig. S25). Actually, for BaP, the modal scheme generally decreases the concentrations in the Arctic (as median, −35 % in DJF and −15 % during JJA) and increases (approx. 5 %) those over the mid- and low-latitude landmass (Figs. S8d and S9d, left panels).

As is the case for the direct effects, the interaction contributions are peculiar to individual species (Fig. d). For PHE, the interaction effects in DJF are reflected in negative concentration responses over the Arctic (18 %) and positive over Antarctica (7 %), in contrast to relatively mild influences over other regions or in JJA. For PYR and FLT, the effects are negative over the Arctic both in DJF (quartiles vary from −20 % to −75 %) and in JJA (−6 % to −120 %). It is interesting to note that the interaction between the modal and ppLFER schemes has a major influence on the negative signal (Figs. S15, S16, S19, S20), suggesting that the decrease in simulated concentration associated with the change from bulk to modal could be intensified when the ppLFER scheme is used. In the remaining areas, the interaction effects vary in sign spatially as illustrated in the right panels of Figs. S8b and c to S9b and c. Nevertheless, it shows for both species that maximum influences occur over the Southern Ocean in DJF (where the effects may reach 2 orders of magnitude) and midlatitude landmass in JJA (more than a factor of 5). As for BaP, the median effects are negative (−7 % to −30 %) in both seasons, although some positive signals are apparent in parts of high latitudes while the tropical oceans bear small synergistic effects. Similar to other species, the degree of interactions is stronger than ${\widehat{f}}_{\mathrm{2}}$ by more than a factor of 3 for the majority of grid cells (Figs. S8d–S9d, right panels). The large fractions of the effects are dominated by 2-factor and 3-factor combinations related to the interaction with the ppLFER scheme and/or re-volatilization (Figs. S17, S21).

3.1.3 Effects of the choice of gas–particle partitioning scheme

Figure e shows that the direct effects of the ppLFER scheme (${\widehat{f}}_{\mathrm{3}}$) show little spatial heterogeneities in both seasons and for all species. The effects are barely important for PHE due to a low gas–aerosol partition constant (K_p); ${\widehat{f}}_{\mathrm{3}}$ is positive for PYR and FLT over polar regions and northern midlatitudes especially in winter when low temperature favors partitioning to aerosols (higher K_p). The median of ${\widehat{f}}_{\mathrm{3}}/f_{\mathrm{0}}$ varies from 1 % to 25 % with some parts of Antarctica showing an increase larger than 50 %. For BaP, the effects are overall negative (by at least −5 %) with ${\widehat{f}}_{\mathrm{3}}/f_{\mathrm{0}}$ reflecting a positive north–south gradient (increasing from the Arctic to Antarctica), associated in part with stronger signals over oceans (Figs. S10d and S11d, left panels). In particular under the modal size discretization, the choice of gas–particle partitioning scheme has only minor effects for atmospheric levels, except for BaP for which model overestimates are compensated by the choice of the ppLFER scheme (Fig. S25). Under the bulk size discretization, the ppLFER scheme tends to enhance some of the overestimates in the Arctic summer (FLT, PYR; Figs. S23–S24). The application of ppLFER increases K_p as this module is calculated from not only interaction with BC and OM (as in Lohmann–Lammel scheme) but also with some other aerosol matrices. Higher K_p indicates higher particle mass fraction. For PYR and FLT, this leads to an increase in total atmospheric lifetime as the aerosol phase is not degraded, and can therefore be transported over a larger distance. For BaP, the additional particles are subject to depositions and heterogeneous oxidation by ozone, particularly in regions away from sources. The factor influence is notably too small for PHE as oxidations occur in both phases.

The effects from fac3 interactions vary by region and are relatively stronger than ${\widehat{f}}_{\mathrm{3}}$ (Fig. f). This finding is common to all species and seasons. The degree of effects is weaker for PHE compared to that for other species. However, the interactions increase polar concentrations in local summer by 20 % to a factor of 5, mainly associated with the coupled effect of ppLFER and volatilization (${\widehat{f}}_{\mathrm{34}}$, Figs. S14 and S18). For PYR and FLT, there is a high spatial variability over extratropical regions in local summer, as indicated by the interquartile range (distance between the third and first quartiles). With regard to synergistic terms, ppLFER interactions with the modal scheme and re-volatilization, in 2- or 3-factor combinations, are more important than other contributions (Figs. S15–S16 and S19–S20). For BaP, the interaction effects show negative signals similar to ${\widehat{f}}_{\mathrm{3}}$, suggesting a positive feedback. The interactions exert a stronger influence on the concentrations of the oceans than on that of land, except in the tropics (Figs. S10d and S11d, right panels). The median of relative effects ranges from −1 % to a factor of −10, minimum (maximum) in the northern (southern) extratropics. Two second-order interactions likely make major contributions, i.e., ${\widehat{f}}_{\mathrm{34}}$ which dominates the response over oceans and ${\widehat{f}}_{\mathrm{23}}$ which dominates over land (Figs. S17 and S21).

3.1.4 Effects of re-volatilization

The direct effects of re-volatilization (${\widehat{f}}_{\mathrm{4}}$) are illustrated in Fig. g. ${\widehat{f}}_{\mathrm{4}}$ is positive in the tropics in both seasons, with the median ${\widehat{f}}_{\mathrm{4}}/f_{\mathrm{0}}$ ranging from 5 % to 50 %. Intensive re-volatilization in this region would increase net surface fluxes, thereby increasing concentrations. For PHE, positive ${\widehat{f}}_{\mathrm{4}}$ values are more localized over the tropical landmass, whereas negative ${\widehat{f}}_{\mathrm{4}}$ values are predicted over the tropical ocean (Figs. S12a and S13a; left panels). The positive (negative) effects over land (ocean) areas are also apparent at higher latitudes during most of the year. This reflects the fact that the negative effects on concentrations over ocean act contrary to the positive effects on net surface fluxes, mainly caused by the nonlinear relationships of air–sea gas exchange (deposition and volatilization), air and surface burden, atmospheric oxidation, and emissions. Accounting for re-volatilization compensates for a significant part of underestimates of PHE in the Arctic during summer, but adds to overestimates in midlatitudes (Fig. S22).

For the studied species of mid semivolatility, PYR and FLT, a positive signal is apparent over the high and middle latitudes during local summer in contrast to a negative signal during local winter (Fig. g). Similar to PHE, the negative signal is confined over oceans (Figs. S12b–c and S13b–c; left panels). The summer increases are stronger (20 % to a factor of 10) than the winter decreases (−10 % to −60 %) and the magnitudes are higher in FLT than in PYR. The near-ground concentrations of PYR and FLT are estimated to be ≈30 %–80 % in midlatitudes of which ≈30 % are attributable to re-volatilization. In the Arctic, re-volatilization compensates for ≈60 % of PYR underestimation (Fig. S23) and explains most of, ≈60 %–80 %, the FLT overestimation (Fig. S24). For BaP, ${\widehat{f}}_{\mathrm{4}}$ is positive consistently across regions and seasons (${\widehat{f}}_{\mathrm{4}}/f_{\mathrm{0}}$ ranges from 20 % to a factor of 10), with substantial effects occurring over oceans (Figs. S12d–S13d; left panels). Accounting for re-volatilization creates some overestimates in the Arctic during summer (Fig. S25). It should be noted that the parameterization adopted here to describe volatilization from soils (the Smit scheme) is derived from an experimental study on mid-polar-to-polar pesticides and there is a need to validate and eventually sophisticate the parameterization to apolar substances.

The interactions generally point toward positive effects for the high-to-medium volatility species (Fig. h), despite some negative effects present over parts of the southern (northern) oceans in DJF (JJA) (Figs. S12–S13; right panels). As for BaP, the effects are uniformly negative, inferring the interactions work in opposition to ${\widehat{f}}_{\mathrm{4}}$. The negative response is almost entirely caused by the negative ${\widehat{f}}_{\mathrm{34}}$, i.e., the 2-factor interaction between re-volatilization and the ppLFER scheme (Figs. S17, S21). Compared to ${\widehat{f}}_{\mathrm{4}}$, the degree of interactions is weaker for PHE, except in polar regions during local summer where the interactions could amplify ${\widehat{f}}_{\mathrm{4}}$. The above implies that ${\widehat{f}}_{\mathrm{4}}$ may point in the right direction regardless of the influences from other factor changes. In contrast, the degree of interactions is overall comparable to ${\widehat{f}}_{\mathrm{4}}$ for the other species.

3.2 Model evaluation

Model performance using the sophisticated realization of the four features (factors), i.e., Seasonal emission + Modal scheme + PpLFER scheme + With re-volatilization (SMPW), is presented below. Two predicted variables are evaluated, i.e., total (gas+particle) concentrations and aerosol particulate mass fraction at the lowest model level. The metrics applied are listed in the Supplement Sect. SV.

Figure 5 Seasonal mean total (gas+particle) concentrations of PAHs (ng m⁻³) from observations (solid lines) and simulations (dashed lines) averaged over all stations in the (a) Arctic, (b) northern midlatitudes, and (c) the tropics. Note that a logarithmic scale has been used for BaP concentrations in the Arctic.

Table 2 Statistics comparison of model simulation and observations of total (gas+particle) concentrations of PAHs from stations in the Arctic, northern midlatitudes, and the tropics. N: Number of observed-simulated monthly data pairs; $\stackrel{\mathrm{‾}}{x}$: mean; Q2_x: median; SD_x: standard deviation; GM_x: geometric mean; x: simulated (M) or observed (O) data; MB: mean bias; RMSE: root mean square error; NMB: normalized mean bias; NMBF: normalized mean bias factor; FAC2: factor of 2; FAC10: factor of 10; r: correlation coefficient.

3.2.1 Near-surface air concentration

Firstly, the comparison to land monitoring stations is as follows.

• Central tendency. Table 2 shows statistical indices for near-surface concentrations of atmospheric PAHs from observations and simulations and their comparisons, averaged across stations in the Arctic, northern midlatitudes, and the tropics. We can see that observed mean concentrations are higher for PHE and smaller for BaP over all regions. Furthermore, the Arctic concentrations are lower than those in the northern midlatitudes by a factor of around 20 and those in the tropics by approx. 2 orders of magnitude. The model captures well these species and regional variations, but the magnitudes are both underestimated and overestimated. In the Arctic, it underestimates PHE (MB $=-\mathrm{0.060}$ ng m⁻³) and BaP (MB $=-\mathrm{0.006}$ ng m⁻³) concentrations but slightly overestimates PYR (MB =0.001 ng m⁻³) and FLT (MB =0.04 ng m⁻³). In the NH midlatitudes, the model overestimates the three species predominantly occurring in the gas phase (MB =0.077–0.867 ng m⁻³) but underestimates BaP (MB $=-\mathrm{0.58}$ ng m⁻³). Negative bias is seen in the tropics for three PAHs (MB $=-\mathrm{3.443}$ to −6.851 ng m⁻³). Nevertheless, the comparison of model and observations at individual monitoring stations can be different from the regional mean statistics, as described in the Supplement Sect. SVIII. Comparing all four PAHs, a larger degree of bias is found for BaP, which increases from the northern midlatitudes (NMB $=-\mathrm{0.58}$, NMBF $=-\mathrm{1.40}$, FAC2 =0.31, FAC10 =0.79) to the Arctic (NMB $=-\mathrm{0.92}$, NMBF $=-\mathrm{12.17}$, FAC2 =0.17, FAC10 =0.33).

• Dispersion of monthly concentrations. In the following, the coefficient of variation (CoV) is used to compare the dispersion of concentrations among species of different ranges. CoV was calculated by dividing the standard deviation (SD) of all data points by its mean value ($\stackrel{\mathrm{‾}}{x}$). The observations show high variability (CoV >1) with CoV ranging between 1.12 and 2.14. The simulated concentrations appear to be less dispersed than the observations (CoV =0.78–1.93) except for the Arctic PHE and PYR concentrations. The degree of underestimation is larger in the tropics with CoV being 30 %–50 % smaller than the observations. Furthermore, correlations between predicted and observed concentrations are weaker than those in other regions where r varies between 0.29 and 0.63 (the model reproduces 8 %–40 % of the variance in observed concentrations). Comparing the four species, the simulated BaP shows greater underpredictions of the variability where CoV values are less than half of those observed and correlations are less than 0.2 (accounting not over 4 % of observed variance). Higher variability in BaP measurements (than in model results) can be influenced by strongly varying emissions in source regions that are not reflected in the emission inventory (Matthias et al., 2009).

• Seasonal variation. Figure  compares simulated and observed seasonal cycle of average concentrations for different species and regions. The observed mean concentrations are largest in winter and lowest during summer because of less emission and the strong presence of OH for oxidation. The winter maximum to summer minimum ratio (amplitude) is more pronounced (by more than a factor of 2) in the Arctic than that in the NH midlatitudes. The seasonality between model and observations is in qualitative agreement, particularly over the Arctic (except in summer) and midlatitudes. In the Arctic, the model overestimates the seasonal amplitude of PHE and BaP and underestimates their mean concentrations. The contrast is seen for PYR and FLT. FLT concentration is overestimated by up to a factor of 3 in summer while PYR is quite well predicted. In the NH midlatitudes, the model underestimates the amplitude but overestimates the concentrations of PHE, PYR, and FLT (by typically a factor of 2), whereas a systematic negative bias is found for BaP. In the tropics, both the amplitude and magnitude are too low in the model (for magnitude, by a factor of 2–5).

Additional findings are discussed in the Supplement Sect. SIX related to the comparison between EMAC model results and those from other global PAH modeling studies.

Secondly, the comparison to ship cruise measurements is as follows. Measurements of PHE, PYR, and FLT concentrations over the Atlantic Ocean were taken during a cruise in July 2009 (Lohmann et al., 2013). Figure 5 shows the ship sample concentrations overlaying the simulated PAH concentrations. Sample arithmetic (geometric) means during the whole cruise transect are 322 (209), 95 (88), and 128 (111) pg m⁻³ for PHE, PYR, and FLT, respectively. The model poorly reproduces the remote marine environments and overall underestimates the observations, except at 3 locations along the North American coast. The simulated means across sampling positions are 23 (7), 20 (3), and 39 (2) pg m⁻³, respectively, and the underestimation ranges from a factor of 2 to 1000. The degree of bias is most apparent over the tropical South Atlantic at latitude band 5–15^∘ S.

As reported in Liu et al. (2014), the measured concentrations of BaP over the Asian marginal seas, the Indian Ocean, the South Pacific Ocean, and North Pacific Ocean are 131 (45), 14 (3), 9 (2), and 8 (3) pg m⁻³, respectively, for the arithmetic (geometric) means of all samples. Similar to other species, the model also underestimates the BaP concentrations with mean values being 75 (15), 4 (0.05), 0.09 (0.03), and 0.2 (0.06) pg m⁻³, respectively. The discrepancy appears relatively smaller over the Asian marginal seas as compared to other locations (Fig. 6). A substantial degree of bias is seen over the Indian Ocean covering approximately the area bounded by 70–90^∘ E and 10–30^∘ S, with simulated values being more than 2 orders of magnitude smaller than the observed.

The model tendency to underestimate the marine air concentrations may likely be due to several factors. (a) The grid resolution is not sufficient to reproduce fine-scale processes at the grid points close to shipping tracks; (b) high uncertainties associated with the air–sea gas exchange parameterizations still exist, most notably in the estimation of gas transfer velocity; (c) the global inventory (Shen et al., 2013) may significantly underestimate emissions from ocean shipping and does not well characterize the spatial and temporal variability of biomass burning plumes as another potential point of origin of pollutants in the marine air (Nizzetto et al., 2008); and (d) PAH concentration over remote oceans is controlled by atmospheric components (e.g., temperature, wind speed, boundary layer height, photochemical degradation) and the dynamical and biogeochemical components of the ocean. However, the ocean components have not been covered in the simulation; (e) the particulate-bound PAHs may undergo too fast heterogeneous oxidation (most relevant for BaP), leading to short atmospheric lifetimes and weaker long-range transport. BaP, mostly stays in the particulate phase, presumably also in seawater, and therefore may be somewhat underestimated due to the neglect of sea-spray-driven aerosol suspension.

Figure 6 Simulated concentrations of PHE, PYR, and FLT (pg m⁻³) over the Atlantic ocean overlaid with concentrations from a ship cruise measurement campaign during July 2009 (triangles). Land grid cells are depicted in gray shades.

Figure 7 Simulated BaP concentrations (pg m⁻³) over the four ocean margins overlaid with concentrations from a ship cruise measurement campaign (triangles). Land grid cells are depicted in gray shades.

3.2.2 Particulate mass fraction

Measurements of particulate mass fraction (θ) were available only from the E3 station in Europe and IADN stations (I1–I7) in North America (see Table S11). Table 3 presents summary statistics on monthly mean θ from observations and simulations including some performance metrics. The observed mean θ is smaller for PHE (0.051±0.035) and higher for BaP (0.949±0.067). This result is expected as volatility decreases (hence θ increases) from (lighter) PHE to (heavier) BaP. The θ values for PYR and FLT are larger by over 5 times than those for PHE and lower by around one-third than those for BaP. The model reproduces well the distinct differences among species but underestimates the observed θ for PHE, PYR, and FLT. The degree of negative bias is relatively large in PHE (NMB $=-\mathrm{0.910}$ and NMBF $=-\mathrm{10.145}$), whereas for the isomer pair of PYR and FLT, the model exhibits a similar performance with a slight improvement in PYR (NMB $=-\mathrm{0.410}$ and NMBF $=-\mathrm{0.694}$). With regard to BaP, there is a satisfactorily small bias (MB =0.015, RMSE =0.074, NMB =0.016, and NMBF =0.016) although the observed and simulated values have a very weak correlation (r=0.03).

Table 3 Statistics comparison of model simulation and observations of particulate mass fraction (θ) from a subset of surface stations, as listed in Table S11. N: Number of observed-simulated monthly data pairs; $\stackrel{\mathrm{‾}}{x}$: mean; SD_x: standard deviation; x: simulated (M) or observed (O) data; MB: mean bias; RMSE: root mean square error; NMB: normalized mean bias; NMBF: normalized mean bias factor; FAC2: factor of 2; FAC10: factor of 10; r: correlation coefficient.

Figure 7 shows the seasonal mean θ averaged over 3 years for all PAHs. Observations show that θ for BaP varies less than those for 3–4-ring PAHs. Although the model adequately reproduces this feature as well as seasonal variation of individual species, the simulated θ of PHE, PYR, and FLT is generally lower than the observations (except for PYR in winter). For BaP, differences between model and observations are less than 10 % in all months. The SVOC submodel describes the gas–particle partitioning of atmospheric SOCs as a function of temperature and aerosol phase composition. The underestimation might be related to the fact that the submodel assumes the particle to be fully in equilibrium with the gas phase at all times. It neglects kinetic limitations of molecular diffusivity that could lead to the trapping of particles inside viscous (or semisolid) organic aerosol coatings. This shielding effect increases equilibration times of the particles, thereby reducing part of θ from the mass available for gas–particle partitioning. Deviations from measurements can also be partly attributed to the locations of some stations that are within, or close to, residential and industrial area (namely I4, I6, and I7) where the scale and gradient in anthropogenic emissions are not resolved by the model grid resolution nor represented by the emission inventory.

Figure 8 Seasonal mean particulate mass fraction (θ; unitless) from observations (solid lines) and simulations (dashed lines).

4 Conclusions

The submodel SVOC has been developed and operated within the EMAC model for the application to global distribution and environmental fate of SOCs. In this first development, the focus was set on the predictions of four PAH species: phenanthrene (PHE), pyrene (PYR), fluoranthene (FLT), and benzo(a)pyrene (BaP). Multicompartmental fate and air–surface exchange processes were included in SVOC. Some novel features in PAH modeling were tested, including seasonality in emissions, the modal scheme for particulate-phase tracer representation, the ppLFER scheme for gas–particle partitioning, and re-volatilization from surfaces. The results indicate that using seasonal emission compensates for model biases in the predictions of more volatile species (PHE), whereas the effects of the modal and ppLFER schemes are of less significance. Re-volatilization increases the near-ground concentrations in air, which is found most significant for species of mid semivolatility (PYR and FLT). The attribution of model response to individual features (factors) is blurred by the nonlinear interactions between two or more factors. The effects of these interactions are found to both reinforce (positive feedback) and suppress (hence negative feedback) the effects of the individual factors.

For near-surface concentrations, model bias varies by region and/or species, being negative (positive) in the Arctic within typically a factor of 2–13 (6 % to a factor of 2) for PHE and BaP (PYR and FLT); positive in the northern midlatitudes for PHE, PYR, and FLT by up to a factor of 3; negative in the tropics (by a factor of 2–3); and largely over ocean up to 3 orders of magnitude. The model adequately reproduces the seasonal variation of the particulate mass fraction (θ), but underestimates θ for high-to-medium-volatility PAHs. This might be related to a systematic underestimation of OC by the model, which neglects secondary organic aerosols (SOA). The latter may cause significant underestimation of the overall atmospheric aerosol burden and θ of SOCs, in particular over ocean. Since recently a MESSy submodel, ORACLE, dedicated to the simulations of SOA (Tsimpidi et al., 2014) based on lumping organic species in volatility bins is available. It should be included in future SOC simulations using EMAC.

Moreover, the implicit assumption of instantaneous gas–particle equilibrium for SVOC may cause both over- and underestimates of θ, as interphase mass transfer my be kinetically limited to gaseous sources (hence overestimate of θ) or within the particle bulk (hence underestimate θ), as the PAHs may become trapped within particles during transport (Friedman et al., 2014; Zelenyuk et al., 2012; Mu et al., 2018). For multidecadal studies, the coupling of a 3-D ocean model (coupled with a marine biogeochemistry module) would be needed since the present model application does not allow for horizontal and vertical transports in the deep ocean. For the same reason, contaminant remobilization within deep soil layers should also be introduced. To this end, a multilayer (3-D) soil compartment would be needed to replace the 2-D soil compartment used here.