A global model perturbed parameter ensemble study of secondary organic aerosol formation

A global model perturbed parameter ensemble of 60 simulations was used to explore how combinations of six parameters related to secondary organic aerosol (SOA) formation affect particle number concentrations and organic aerosol mass. The parameters represent the formation of organic compounds with different volatilities from biogenic and anthropogenic sources. The most plausible parameter combinations were determined by comparing the simulations against observations of the number concentration of particles larger than 3 nm diameter (N3), the number concentration of particles larger than 50 nm diameter (N50), and the organic aerosol (OA) mass concentration. The simulations expose a high degree of model equifinality in which the skill of widely different parameter combinations cannot be distinguished against observations. We therefore conclude that, based on the observations we have used, a six-parameter SOA scheme is under-determined. Nevertheless, the model skill in simulating N3 and N50 is clearly determined by the low-volatility and extremely low-volatility compounds that affect new particle formation and growth, and the skill in simulating OA mass is determined by the low-volatility and semi-volatile compounds. The biogenic low-volatility class of compounds that grow nucleated clusters and condense on all particles is found to have the strongest effect on the model skill in simulating N3, N50, and OA. The simulations also expose potential structural deficiencies in the model: we find that parameter combinations that are best for N3 and N50 are worst for OA mass, and the ensemble exaggerates the observed seasonal cycle of particle concentrations – a deficiency that we conclude requires an additional anthropogenic source of either primary or secondary particles.

1 Introduction

About 20 %–50 % of lower-tropospheric fine aerosol mass in continental mid-latitudes and almost 90 % in the tropics is composed of organic material (Kanakidou et al., 2005). The major fraction of organic aerosol has been found to be secondary (Zhang et al., 2007; Jimenez et al., 2009), formed as a result of atmospheric oxidation of volatile organic compounds (VOCs) leading to secondary organic aerosol (SOA). Estimating atmospheric SOA is important for accurately estimating the anthropogenic aerosol radiative forcing (Maria et al., 2004; Tsigaridis and Kanakidou, 2007; Carslaw et al., 2010; Riipinen et al., 2011; Makkonen et al., 2012; Shrivastava et al., 2017; Tsigaridis and Kanakidou, 2018). Despite the importance of SOA for climate, a comprehensive characterization of atmospheric VOCs, their reaction pathways, and their SOA formation potential has not yet been possible. Characterization of VOCs is challenging because of the very large number of compounds involved and their diverse sources; tens of thousands of VOCs have been identified in the atmosphere and yet more still remain to be detected (Goldstein and Galbally, 2007). For the SOA precursor gases that have been identified, questions remain about their emission sources, their chemical conversion to SOA, and the effects of atmospheric chemical composition and oxidants on SOA formation.

From a modelling perspective a further challenge is how to deal with the enormous chemical complexity of the system to adequately parameterize laboratory and observational results and incorporate them in large-scale models. Models traditionally use laboratory measurements of reaction rate constants and product yields to calculate the production of highly oxidized VOCs from the reaction between precursor VOCs and atmospheric oxidants. Only a very small range of natural VOCs are accounted for, such as monoterpenes or isoprene. Uptake of SOA onto particles is then simulated either time-dependently (often called kinetic uptake) or assuming thermodynamic partitioning (Riipinen et al., 2012). Such “bottom-up” approaches to simulating SOA predict a global SOA budget at the lower end of the total uncertainty range (Kanakidou et al., 2005), most likely because simulations do not capture the full range of VOCs and atmospheric oxidative pathways that lead to a range of products of varying volatilities (Hallquist et al., 2009; Tsigaridis et al., 2014). Studies like Heald et al. (2005), Johnson et al. (2006), and Spracklen et al. (2011) show that models have consistently and significantly underestimated SOA concentrations in different parts of the atmosphere. Tsigaridis et al. (2014) show models largely underestimate the amount of organic aerosol present in the atmosphere with the underestimation being strongest in urban regions based on a study involving 31 global models.

It is well-established that atmospheric organic molecules strongly affect the number concentrations of climate-relevant sized particles by condensing on and growing aerosol particles (Riipinen et al., 2011) or by promoting particle formation (O'Dowd et al., 2002; Zhang et al., 2004; Metzger et al., 2010). Recent studies (Ehn et al., 2014; Jokinen et al., 2015; Kirkby et al., 2016; Tröstl et al., 2016; Bianchi et al., 2019) have established the importance of atmospheric highly oxygenated organic molecules in the formation and growth of new aerosol particles. Subsequent model simulations based on these experimental data have shown that new particle formation involving organic molecules could explain the seasonal cycle of particle concentrations (Riccobono et al., 2014) and provide a source of aerosol in clean pre-industrial environments that is important for climate (Gordon et al., 2016). It has been estimated that global cloud condensation nuclei (at 0.2 % supersaturation) would be about one-quarter lower without biogenic VOCs, and about three-quarters of this effect is caused by the role of organic molecules in nucleation and early growth (Gordon et al., 2017).

The importance of SOA for climate means that large-scale models need to simulate the contributions of these highly oxygenated molecules in nucleation or subsequent particle growth. Several modelling studies (such as Farina et al., 2010; Pye and Seinfeld, 2010; Jathar et al., 2011) have implemented the volatility basis set (VBS) framework proposed by Donahue et al. (2012) for the description of organic partitioning and chemical ageing. Other SOA modelling frameworks have been proposed by Odum et al. (1996), Camredon et al. (2007), and Parikh et al. (2011). However increased model complexity inevitably requires more computational resources and more runtime, both of which need to be considered carefully by modellers. Shrivastava et al. (2011) found comparable predictions of observed OA between a nine-species VBS approach and a simplified two-species VBS approach, the latter being a factor of 2 lower in computational cost. Shrivastava et al. (2011) concluded the two-species approach is well-suited to represent the complex evolution of atmospheric organic aerosols. Riipinen et al. (2011) and Scott et al. (2015) propose a simplified SOA formation scheme representing species of two volatilities. Tsigaridis et al. (2014) show global model skill does not increase with model complexity with regard to organic aerosol mass concentrations.

There is the risk that added complexities in models that are not well-constrained by experimental data could increase model uncertainty and thereby introduce more uncertainty in the quantification of the effect of anthropogenic aerosols (Lee et al., 2016). Although model complexity is increased in order to improve the representation of the physical, chemical, or optical properties of SOA (Tsigaridis et al., 2014), a more complex model that matches some observations may not have lower uncertainty in making predictions because of the increased likelihood of compensating errors, often called model equifinality (Beven, 2006). The simulated particle number (and anthropogenic aerosol forcing) is affected by several model parameters and their associated uncertainties, some of which may be compensating for each other. Tuning any one parameter within the model (e.g. the nucleation rate) to improve the model performance against observations of one aspect of the atmospheric aerosol distribution (e.g. the total particle number concentration) may adversely affect the model performance in other outputs (e.g. total aerosol mass). Such tuned observationally constrained models give the impression of low aerosol uncertainty and model robustness but still predict a large range of aerosol forcing in Lee et al. (2016).

In this paper we use a perturbed parameter ensemble of 60 simulations of a global aerosol model to examine the effect of six uncertain parameters in SOA formation on simulated organic aerosol mass and particle number concentrations. We compare three model outputs against observations: the number concentration of particles larger than 3 nm dry diameter (N3) and 50 nm dry diameter (N50) and the mass concentration of organic aerosol (OA). Our primary aim is to understand the sensitivity of N3, N50, and OA to the combinations of model input parameters that control the formation of SOA. We identify parts of parameter space that result in the best agreement with observations of N3, N50, and OA. Further, we identify parameter combinations that produce models that are indistinguishable in terms of their simulations of number and mass concentrations.

2 Methods

2.1 GLOMAP global aerosol model

We use the GLOMAP (Global Model of Aerosol Processes) modal aerosol microphysics model (Mann et al., 2010), which is an extension to the TOMCAT 3-D chemical transport model (Chipperfield, 2006). The model has a horizontal resolution of 2.8^∘ × 2.8^∘, with 31 hybrid σ-pressure levels from the surface to 10 hPa. Large-scale atmospheric transport in the model for 2008 is driven by ERA-Interim reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) at 6-hourly intervals. The aerosol distribution is simulated using four hydrophilic modes (nucleation, Aitken, accumulation, and coarse) and one non-hydrophilic Aitken mode. The aerosol phase has four components: sulfate, sea salt, black carbon, and organic carbon. Where more than one species is contained in a mode, we assume internal mixing (Mann et al., 2010).

2.2 SOA scheme

The model described in Mann et al. (2010) includes one SOA species produced from oxidation of monoterpenes only. In this study we produce six SOA species from oxygenated organic compounds derived from the oxidation of monoterpene, isoprene, and anthropogenic sources. Monthly monoterpene and isoprene emissions used in the model are generated from the Community Land Model by Sarah Monks (MEGANv2.1; Guenther et al., 2012) Emission of anthropogenic VOC is parameterized using the approach in Spracklen et al. (2011) – we use CO emissions from the MACCity inventory and assume a $\mathrm{SOA}/\mathrm{OC}$ mass ratio of 1.4. Atmospheric oxidants (OH•, O₃, NO₃•) are taken from 6-hourly monthly mean values calculated offline from a TOMCAT simulation and interpolated to the model chemical time step (Monks et al., 2017).

Figure 1 shows a schematic of the treatment of SOA in this study. Gas-phase oxygenated organic compounds (ox-VOCs) are represented by three classes: extremely low-volatility (ELVOC), low-volatility (LVOC), and semi-volatile (SVOC) organic compounds.

Figure 1 Schematic showing the SOA formation scheme in GLOMAP-mode and the six oxidized VOCs (products of photochemical oxidation of emitted VOCs that eventually produce SOA) whose concentrations are perturbed in this study. The six oxidized VOCs represent three volatility categories – extremely low-volatility organic compounds (ELVOCs, blue), low-volatility organic compounds (LVOCs, green), and semi-volatile organic compounds (SVOCs, red). Prefix “A” indicates precursor gases of anthropogenic origin and “B” indicates biogenic origin. The schematic shows the precursor gases and oxidants that react to produce these ox-VOCs, their relative volatility (ELVOC < LVOC < SVOC), and the mechanism (nucleation for ELVOCs, kinetic condensation for LVOCs, and mass-based partitioning for SVOCs) by which they add to the condensed phase (represented here by the five modes: nucleation soluble, Aitken soluble, accumulation soluble, coarse soluble, and Aitken insoluble modes).

The ELVOCs are assumed to derive only from biogenic sources and nucleate to form new particles that either grow or are lost to pre-existing larger particles (Gordon et al., 2016). The LVOCs are assumed to condense kinetically and irreversibly (i.e. with zero vapour pressure) on all particles, and the SVOCs are assumed to partition into all particles, except those in the nucleation mode, in proportion to the pre-existing organic mass in the mode (Scott et al., 2015). The LVOCs and SVOCs are further divided into biogenic (prefix B) and anthropogenic (prefix A). The biogenic precursors are split into monoterpenes (suffix M) and isoprene (suffix I), with the monoterpenes allowed to form both LVOCs and SVOCs and the isoprene forming only SVOCs. The precursors, oxidants, and role of each ox-VOC in SOA formation are defined in Table 1. The ox-VOCs in this new scheme of SOA formation produced from bimolecular oxidation reactions of VOCs are B_ELVOC, B_LVOC, B_SVOC_M, B_SVOC_I, A_LVOC, and A_SVOC. The six perturbed parameters in this study are scaling factors or reaction yields that control the concentrations of these six ox-VOCs.

Table 1 List of oxidized VOCs implemented in this study, the volatility class they represent, whether produced from biogenic or anthropogenic sources, how they take part in atmospheric SOA formation, parent VOC, and oxidants that react to produce each ox-VOC. MT stands for monoterpene, IP stands for isoprene, and CO stands for carbon monoxide indicating anthropogenically sourced VOC.

2.3 Microphysical processes

SOA formation in the model starts with B_ELVOC and sulfuric acid via nucleation. Nucleation rates in the model (Kirkby et al., 2016; Gordon et al., 2016) determine the formation of clusters of 1.7 nm dry diameter. Thereafter their growth up to 3 nm sizes in the model is determined using the equation of Kerminen and Kulmala (2002), which takes into account the losses during initial growth of clusters. Clusters that reach a dry diameter of 3 nm are added to the nucleation mode as freshly nucleated particles. Thus N3, the number concentration of particles with dry diameter greater than 3 nm (in reciprocal cubic centimetres), represents the total particle number concentration in the model.

Once particles appear in the nucleation mode they may either grow using sulfuric acid and ox-VOCs (as described above) or get scavenged by larger particles acting as a condensation sink. Particles that reach a dry diameter of 50 nm can act as cloud condensation nuclei in the atmosphere. Thus N50, the number concentration of particles with a dry diameter greater than 50 nm (in cm⁻³), represents the number of climate-relevant sized particles in the model.

Aerosol particles are removed through dry deposition, sedimentation, nucleation scavenging, and impact scavenging. Dry deposition accounts for gravitational settling, Brownian motion, impaction interception, and particle rebound and predominantly removes particles smaller than 50 nm. Processes represented under wet deposition are nucleation scavenging and impact scavenging.

2.4 Perturbed parameter ensemble

Simultaneous perturbation of six model parameters within each of their chosen ranges forms a 6-D parameter space within which we explored the competing and compensating effects of these parameters on model outputs. An ensemble of 60 simulations was produced, each with a different combination of the six parameters. Simulations were run for the year 2008. To ensure good coverage of the 6-D parameter uncertainty space, the maximin Latin hypercube sampling technique (McKay et al., 1979) was used to choose the combinations of parameters. The Latin hypercube sampling method samples the entire 6-D parameter space and creates a space-filling design of 60 combinations of parameter values such that the distance between any two points in the 6-D space is maximized. The effectivity of the Latin hypercube sampling method to create space-filling parameter combinations has been discussed in Lee et al. (2011, 2013). The relative variation in the ox-VOCs in each simulation is shown in Fig. 2.

We perturb the concentrations of SOA-producing ox-VOCs by perturbing the yields of bimolecular oxidation reactions for B_LVOC, B_SVOC_M, B_SVOC_I, A_LVOC, and A_SVOC. The baseline molar yields for each of these ox-VOCs before perturbation were 13 % for B_LVOC, 13 % for B_SVOC_M, and 3 % for B_SVOC_I, producing approximately 40 Tg yr⁻¹ of SOA from biogenic sources. The total anthropogenic ox-VOCs are split equally between A_LVOC and A_SVOC (effectively a 50 % yield of the total anthropogenic each), together producing approximately 63 Tg yr⁻¹ of SOA. Within the ensemble the yield of each of the above is perturbed from 0 to 20 times the baseline for biogenic ox-VOCs (B_LVOC, B_SVOC_M, B_SVOC_I) and from 0 to 5 times the baseline for anthropogenic ox-VOCs (A_LVOC and A_SVOC). We perturb the concentrations of ox-VOCs taking part in the model nucleation scheme, B_ELVOC, using a scaling factor between 0 and 10, where 0 is equivalent to sulfuric-acid-only nucleation and 10 is equivalent to ELVOC concentrations about a factor of 10 higher than reported in Kirkby et al. (2016). The ranges (Table 2) were chosen to encompass a wide range of uncertainties and simplifications in the model. These include VOC emission uncertainties, structural uncertainty such as neglected precursor gases (e.g. sesquiterpenes), uncertainty in yields of the ox-VOCs from oxidation reactions, simplifications to the oxidation pathways (in GLOMAP-mode only single-stage oxidation products are represented), and uncertainty in SOA due to neglecting the volatility distribution and re-evaporation of SOA (Donahue et al., 2011, 2012). Changing the ox-VOC concentrations by changing the yields of chemical reactions has the same effect in the model as perturbing the emissions of the parent VOCs. A yield above 100 % in the bimolecular reactions should therefore be interpreted as an increase in the total concentration of reactants. The advantage of varying the yield of ox-VOCs rather than the emissions of VOCs is that a perturbation applied to one ox-VOC does not affect the production or loss of the other ox-VOCs when they have the same parent VOC (i.e. they have uncorrelated effects across the 6-D parameter space).

Table 2 Perturbed parameters and ranges. Also listed are the precursor VOC gases that produce the ox-VOCs, the value of ox-VOC yields in the default model version (unperturbed), how the perturbation for each parameter is implemented in the model (“absolute” for replacing the default yield value by the perturbation, “scaled” for scaling the default yield value by the perturbation), and the amount of SOA produced from each parameter in the default model version.

Figure 2 The relative variation in the six perturbed parameters (%) for each ensemble member (numbered 1 to 60). Each hexagon (grey dashed area) represents the 6-D parameter space and the positions of the black dots show the position of each parameter within its range for the specific ensemble member. The dots are joined and shaded green for easy identification of explored parameter space in each ensemble member. Anticlockwise from top, the black dots represent parameter settings for B_ELVOC, B_LVOC, B_SVOC_M, B_SVOC_I, A_LVOC, and A_SVOC respectively. Example interpretation: in simulation 19 (fourth row, first hexagon) B_SVOC_I and A_SVOC concentrations are towards the lower ends of the respective ranges being explored for each of them while concentrations of A_LVOC are towards the high end of the A_LVOC range.

2.5 Observations

Figure 3 shows a map of ground-based observation stations used to compare the surface-level number concentrations of particles with dry diameter greater than 3 nm (N3 in reciprocal cubic centimetres), particles with dry diameter greater than 50 nm (N50 in reciprocal cubic centimetres) and the mass concentration of organic aerosol (OA in microgrammes per cubic metre) predicted by the ensemble.

N3 observations cover 34 ground stations worldwide as used in Spracklen et al. (2010). Measurements of N3 were made between 1994 and 2009 using either condensation particle counters (CPCs), scanning mobility particle sizers (SMPSs), differential mobility particle sizers (DMPSs), or diffusion aerosol spectroscopes (DASs). The N3 data set is fully described in Spracklen et al. (2010) and has been used in previous studies such as Riccobono et al. (2014), Gordon et al. (2016), and Dunne et al. (2016). N50 observations are from 31 ground stations worldwide including sites in Europe as described in Asmi et al. (2011) and at African (Vakkari et al., 2013), Indian (Hyvärinen et al., 2010), Canadian (Jeong et al., 2010; Takahama et al., 2011), and polar (Asmi et al., 2016; Hansen et al., 2009) sites obtained from individual projects and online data portals. Measurements of N50 were made between 2007 and 2015 using DMPS or SMPS instruments. OA observations cover 41 ground stations worldwide from the Global Aerosol Synthesis and Science Project (GASSP) database (Reddington et al., 2017). OA measurements were made between 1990 and 2015 using the aerosol mass spectrometer (AMS) and the associated aerosol chemical speciation monitor (ACSM) that characterize the mass and chemical composition of particulate matter (Canagaratna et al., 2007; Ng et al., 2011).

All station data were averaged to create monthly mean values. Station heights were matched to model pressure levels for each month using barometric altitude. Stations cover a wide range of atmospheric conditions such as continental boundary layer (CBL, e.g. Hyytiälä, Harwell, Botsalano), marine boundary layer (MBL, e.g. Mace Head, Trinidad Head, Sable Island), and free tropospheric (FT, e.g. Nepal, Jungfraujoch, Pico Espejo, Mauna Loa) sites. Errors in measurements are estimated to be around 30 % on average, depending strongly on the spatial heterogeneity of sources (Reddington et al., 2017).

Figure 3 Locations of ground-based sites where model–observation match is compared for N3 (34 locations, symbols in light), N50 (31 locations, symbols in blue), and OA (41 locations, symbols in red).

To examine the performance of the ensemble members, we use statistical metrics including correlation coefficient, normalized mean bias factor (NMBF), and Taylor skill score (TSS) (Taylor, 2001), in the following sections. NMBF is an unbiased and symmetric metric with a range from −∞ to +∞, with 0 corresponding to exact agreement. It is calculated using the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{NMBF}=\left\{\begin{array}{l}\frac{{\sum }_{i=\mathrm{1}}^N(M_i-O_i)}{{\sum }_{i=\mathrm{1}}^N{O}_i},\text{ if }\overline{M}\ge \overline{O}\\ \frac{{\sum }_{i=\mathrm{1}}^N(M_i-O_i)}{{\sum }_{i=\mathrm{1}}^N{M}_i},\text{ if }\overline{M}<\overline{O},\end{array}\right.\end{array}$$

where for ensemble members $i=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots }N\mathit{\}}$, M_i and O_i are the modelled and observed variables, the $\overline{M}$ represents the mean across all ensemble members, and $\overline{O}$ represents the mean across all observations. NMBF $=\mathrm{1.5}(-\mathrm{1.5})$ means the model is biased towards overestimating (underestimating) observations by a factor of 2.5 (NMBF+1).

The Taylor skill score (TSS) is calculated as

$$\begin{array}{}\text{(2)}& \mathrm{TSS}={\displaystyle \frac{\mathrm{4}(\mathrm{1}+R{)}^{\mathrm{4}}}{(\widehat{{\mathit{\sigma }}_f}+\frac{\mathrm{1}}{\widehat{{\mathit{\sigma }}_f}}{)}^{\mathrm{2}}(\mathrm{1}+R_{\mathrm{0}}{)}^{\mathrm{4}}}},\end{array}$$

where $\widehat{\mathit{\sigma }}$ is the normalized standard deviation (${\mathit{\sigma }}_{\mathrm{model}}/{\mathit{\sigma }}_{\mathrm{observation}}$), R is the Pearson correlation coefficient, and R₀ is the maximum correlation attainable by the model, assumed to be 1. As the model variance approaches the variance in the observations and R approaches R₀, i.e. the model is most skilful, TSS approaches unity. As the model variance approaches zero or as the correlation coefficient between model and observation becomes more negative, TSS approaches zero. TSS thus takes into account both how well the model simulates the observed pattern (correlation coefficient) and how close model observation agreement is (variance). The full statistics (TSS, NMBF, and R) calculated for each simulation within the ensemble is presented in Table A1.

3 Results

3.1 Global and regional aerosol mass and number

The global mass of SOA produced in the model simulations ranges from 220 to 850 Tg yr⁻¹ across the 6-D parameter range. Our global SOA range covers the upper end of the 50–380 Tg yr⁻¹ global SOA found in Spracklen et al. (2011) after constraint of an earlier version of the GLOMAP model to global AMS OA observations. Although individual ox-VOC yields were varied between 0 and 20 times their baseline yields, the total global production of SOA varies by only a factor of 4, and the lowest value is 220 Tg yr⁻¹. Global mean N50 varies by a factor of 2.39 within the ensemble (from 377 to 903 cm⁻³). Modelled N3 is more sensitive to changes in ox-VOCs than N50; global mean values vary by a factor of 3.5 across the 60 ensemble members (from 531 to 1889 cm⁻³).

Figure 4 shows that the ensemble members produce large regional variations in OA, even when they predict similar global values. For example, simulations 9 (subplot 4.3) and 36 (subplot 5.3) have similar global mean OA concentrations (6.47 and 6.49 µg m⁻³ respectively) but they simulate very different OA over the highly polluted regions of South Asia. Such regional variations are dependent on the parameter settings of the ox-VOCs. In this case simulations 9 and 36 particularly differ in the contribution from B_SVOC_M and A_LVOC. Modelling efforts need to focus on capturing the competing and compensating effects of ox-VOCs contributing to different stages of particle formation and growth, rather than detailed representation of any one contribution. This also emphasizes the need to compare model outputs with regional as well as global metrics and observations to determine whether the model is performing well.

Simulations 45 (subplot 2.9) and 49 (subplot 5.10) demonstrate the role of anthropogenic VOCs on the simulated aerosol size distribution and OA mass. Global mean N3, N50, and OA in simulations 45 and 49 are all in the upper quartile of the ensemble's output distribution for these quantities. Both of these simulations have high concentrations of B_ELVOC, which promotes nucleation, and moderate to high B_LVOC, which promotes the survival of nucleated particles, contributing to the high global mean particle number and mass concentrations. In simulation 45 both biogenic and anthropogenic ox-VOCs contribute significantly, with B_ELVOC, B_SVOC_I, and A_LVOC concentrations being predominant (Fig. 2). In contrast contributions to SOA are predominantly biogenic in simulation 49 (Fig. 2). This difference in parameter design is reflected in the global distributions. Simulation 49 with low anthropogenic contributions predicts lower OA concentrations in the highly polluted SE Asian region. Anthropogenic ox-VOCs in the model favour the loss of smaller particles (which are more susceptible to condensation sinks) and the growth of larger particles (mass-based partitioning of A_SVOC). Therefore in the simulation with low anthropogenic SOA more smaller particles survive but fewer particles hold substantial mass, leading to lower OA mass concentrations. This is also supported by the considerably higher N3 but only slightly higher N50 number concentrations in the SE Asian region in simulation 49 compared to simulation 45. Anthropogenic sources of SOA affect the pre-industrial and present-day atmospheres differently in global models (Carslaw et al., 2013). Under-represented anthropogenic sources in global models therefore have greater climatic implications than under-represented biogenic SOA sources.

Figure 4 Global annual mean anomaly of organic aerosol mass at the surface (OA in microgrammes per cubic metre) produced within the ensemble. Each subplot shows the anomaly of an ensemble member (numbered between 1 and 60) from the ensemble mean OA. The global mean OA is given on each subplot. The subplots are arranged in order of increasing global mean OA.

Figure 5 Global annual mean anomaly of N3 number concentration at the surface (cm⁻³) produced within the ensemble. Each subplot shows the anomaly of an ensemble member (numbered between 1 and 60) from the ensemble mean N3. The global mean N3 is given on each subplot. The subplots are arranged in order of increasing global mean N3.

Figure 6 Global annual mean anomaly of N50 number concentration at the surface (cm⁻³) produced within the ensemble. Each subplot shows the anomaly of an ensemble member (numbered between 1 and 60) from the ensemble mean N50. The global mean N50 is given on each subplot. The subplots are arranged in order of increasing global mean N3.

Figures 5 and 6 show the global distribution of N50 and N3. The panels in each figure are ordered according to increasing global mean N3. There is a general increase in N50 with increasing N3, but 11 of the simulations clearly show low N50 concentrations despite high N3 (simulations 14, 39, 47, 56, 50, 38, 9, 15, 59, 51 in Fig. 6). All of these simulations have one aspect in common – very low yields of B_LVOC (Fig. 2). High concentrations of nucleating B_ELVOC in these simulations initially facilitate the formation of particles, but low yields of B_LVOC, especially when coupled with one or more of high B_SVOC_I, A_LVOC, and A_SVOC, suppress particle growth up to 50 nm. B_SVOC_I, A_LVOC, and A_SVOC are not spatially co-located with the nucleated clusters produced from biogenic ox-VOCs and facilitate the growth of larger particles, which then increases the condensation sink for nucleated clusters, thereby effectively suppressing the growth of nucleated clusters to N50-relevant sizes in the above simulations with low B_LVOC yields.

In contrast, simulations 13, 24, 35, and 46 have low concentrations of B_ELVOC but relatively high B_LVOC. Despite the low B_ELVOC concentrations, which produces fewer nucleated particles, the relatively high concentration of B_LVOC ensures that more of the nucleated particles reach 50 nm diameter in these simulations. Consequently for these simulations the simulated global mean N3 concentrations are in the lower quartile within the ensemble but global mean N50 concentrations are in the interquartile range within the ensemble. B_LVOC can compensate for B_ELVOC to some extent and clearly stands out as the most important controller of climate-relevant particle number concentrations.

Including new, more accurate nucleation pathways into models is unlikely to improve the model performance with respect to N50 (a highly relevant model output for estimation of climate-relevant aerosol–cloud interactions) unless the models also include adequate representation of B_LVOCs. Several studies have investigated nucleation capability and nucleation pathways of atmospheric molecules (Kulmala et al., 1998, 2004; Kirkby et al., 2011; Kurtén et al., 2008; Almeida et al., 2013; Riccobono et al., 2014; Kirkby et al., 2016), whereas the contribution of organic molecules to sub-3 nm cluster growth is relatively recent knowledge, and the molecules involved are largely unidentified (Tröstl et al., 2016). The significance of B_LVOC is explored and established further in later sections.

Figure 7 shows how the global mean N50 and OA concentrations depend on the six ox-VOC parameters. The shading, blue for OA and red for N50, indicate whether the global mean values within the ensemble fall in the upper quartile (Q4 indicated by dark shade), inter-quartile range (IQR indicated by medium shade), or lower quartile (Q1 indicated by light shade). The colour coding clearly shows the multi-variate relationship between simulated N50, OA, and ox-VOC parameters. High values of OA may be associated with various combinations of ox-VOC parameters and with both high or low particle number concentrations. The figure is consistent with the challenge faced by state-of-the-art global models – despite simulating the particle number concentrations (N3 and N50) reasonably well models consistently under-predict OA mass concentrations (Kanakidou et al., 2005; Tsigaridis et al., 2014). The challenge in predicting both particle number concentrations and OA mass is re-visited in later sections.

N50 concentrations in Fig. 7 have a strong relationship with B_ELVOC up to about 5 times the ELVOC yield of 3.2 % (Kirkby et al., 2016) used in the model, above which there is more scatter in N50 caused by the other model processes and parameters. N50 is also related to B_LVOC, but with more scatter than for B_ELVOC. These relationships show that N50 concentrations are strongly controlled in part by the production of B_ELVOC, which causes nucleation, and by the production of B_LVOC, which grows the nucleated clusters via kinetic condensation. There is no clear relationship between N50 and A_LVOC production. The likely reason for this is that anthropogenic VOCs are not spatially co-located with the biogenically produced B_ELVOC, so there are fewer nuclei in polluted regions and hence much less effect of the A_LVOC on the growth of nuclei to larger sizes.

OA concentrations in Fig. 7 are found to be unrelated to B_ELVOC concentration, showing that new particle formation has little effect on simulated OA mass in our model. Increases in all the other ox-VOC parameters generally increase OA, although there is a lot of scatter, particularly with the anthropogenic parameters, indicating a strong multi-variate relationship for simulated OA. Global mean OA shows the strongest dependence on B_LVOC and B_SVOC_M, and the highest global mean OA (darkest red in Fig. 7) is simulated when B_LVOC and B_SVOC_M yields are more than 7.5 times the baseline yield of 13 % (above 100 % in Fig. 7; producing over 113 Tg yr⁻¹ each).

We make two more observations from Fig. 7. The simulated OA mass concentrations seem to have a stronger relationship with SVOCs than LVOCs, because SVOCs partition to larger particles which already hold substantial mass, thereby having a greater impact on OA mass. Secondly OA concentrations appear to have a steeper increase with increases in the biogenic ox-VOCs than their anthropogenic counterparts. This is because in the current SOA scheme A_LVOCs or A_SVOCs grow fewer particles than their biogenic counterparts which have the same spatial distribution as the nucleated particles they produce. As a result, changes in the concentrations of anthropogenic ox-VOCs have a lesser impact on simulated OA mass. The involvement of anthropogenic precursors in particle formation and cluster growth is likely to change this picture (Molteni et al., 2018).

Overall we find that when particle number concentrations are low, the main difference between simulations that produce a high amount of OA (simulations 16, 23, 35, and 46 in Fig. 4) and those that do not (see simulations 34, 55, and 58 in Fig. 4) is the relative concentrations of B_LVOC, which grows freshly nucleated clusters before they can be scavenged by coagulation (see Fig. 2 for parameter combinations). When particle number concentrations are high (due to high B_ELVOC or B_LVOC or both), the mass of OA produced is determined by the combined effects of all other ox-VOCs. Parameter combinations such as in simulations 45 and 49 produce some of the highest global mean N3, N50, and OA within the ensemble. In parameter combinations in which B_ELVOC and/or B_LVOC dominate significantly over the rest of the ox-VOCs (such as simulations 8, 21, 41, 54, and 57), the increased competition between small particles for growth cause the available high-volatility ox-VOCs to distribute on more particles, causing a smaller increase in the particle mass.

Figure 7 Global annual mean OA (top panel) and N50 (bottom panel) against the perturbed range of the six ox-VOCs for the 60 ensemble members. The three shades of blue in the top panel divide the ensemble members into three categories according to the mean N50 concentrations: Q1_N50 for N50 values in the lower quartile, IQR_N50 for the inter-quartile range, and Q4_N50 for N50 values in the the upper quartile. The three shades of red in the bottom panel depict the same as above but for global mean OA concentrations in the ensemble.

3.2 Ensemble comparison with observations – structural deficiencies in the model

Figure 8 shows the monthly-mean time series of N3 at 34 ground sites. The 60-member ensemble is able to encompass the observations in all months at only 3 out of 34 sites. In most locations the annual mean model bias in each of the 60 ensemble members ranges between a factor of 3 underestimation to a factor of 2 overestimation, with underestimation being more prevalent. Riccobono et al. (2014), Dunne et al. (2016), and Gordon et al. (2016) have previously reported underestimation of wintertime particle number concentrations by the model. Here we find the wintertime underestimation continues even when B_ELVOC or B_LVOC ox-VOC parameters are at their highest settings. In addition to underestimation of modelled concentrations in winter, some ensemble members overestimate particle concentration in the summer (for example Aspvreten in Fig. 8). These combined biases mean that the model overestimates the strength of the seasonal cycle compared to observations (see Appendix Figs. A1 and A2).

Figure 9 shows the monthly-mean time series of N50 at 31 ground sites. The 60-member ensemble is able to encompass the observations in all months at only 2 out of 31 sites. Like N3, the model bias for N50 in each of the 60 ensemble members ranges between a factor of 3 underestimation and a factor of 2 overestimation, with underestimation being more prevalent. The correlation coefficients (calculated for each ensemble member at each location using monthly mean simulated and observed N50 concentrations) in 22 out of 31 sites are higher than 0.5 (figure not shown). The best correlation coefficients are observed in non-urban sites and the maximum underestimation and poorest correlation coefficients for N50 are observed at the polluted sites of Ispra and Marikana (Asmi et al., 2011; Vakkari et al., 2013). In contrast the ensemble performs significantly better at Hohenpeissenberg and Zugspitze, both high-altitude sites free from nearby anthropogenic influence only about 458 km from Ispra, and at Botsalano, representing a semi-clean environment, which is about 150 km from Marikana.

We also find that as the normalized mean bias factor (calculated between each simulation and observed values) decreases, the calculated correlation coefficient weakens. This implies that the model has a structural deficiency that cannot be resolved by perturbing the model parameters. The strong link between B_ELVOC and N3 indicates that the weakening of the correlation coefficient with the improvement of normalized mean bias factor is related to nucleation: higher ELVOC production rates increase annual mean N3 concentrations, but the summer concentrations are affected much more than the winter concentrations, which weakens the correlation. We suggest that a missing particle source such as anthropogenic pollutants – which are at a maximum in the winter due to low boundary layer height and increased local emissions from sources such as domestic heating – will rectify the model bias significantly. Alternatively, the poor model performance may be improved by exploring uncertainties in other parts of the model unrelated to SOA (Lee et al., 2013; Yoshioka et al., 2019).

Figure 8 Seasonal cycle of simulated (solid coloured lines) and observed (black stars) monthly mean surface-level N3 concentrations at 34 ground-based sites. Each coloured line is an ensemble member.

Figure 9 Seasonal cycle of simulated (solid coloured lines) and observed (black stars) monthly mean surface-level N50 concentrations at 31 ground-based sites. Each coloured line is an ensemble member.

3.3 Model skill across the 6-D parameter space

We now explore how the model skill in simulating observed N3, N50, and OA varies across the 6-D ox-VOC parameter space (Figs. 10, 11, and 12). Each scatter plot in Figs. 10, 11, and 12 shows the relationship between two ox-VOCs for each of the 60 ensemble members. Note the values for B_ELVOC yields shown in Figs. 10, 11, and 12 and Figs. A3, A4, and A5 are scaling factors which have to be multiplied to the baseline model yields (Table 2) to get the B_ELVOC yields for the ensemble members. For the rest of the ox-VOCs the values shown are yield values in percent which can be converted to teragrammes per year using Table 2. Because it is a 6-D space, it is also important to note that the other four parameters are varying randomly across each plane.

We use Figs. 10, 11, and 12 to identify patterns of dependencies of the Taylor skill scores for N3, N50, and OA on the ox-VOC yields within the 6-D parameter space. A weak dependency between an ox-VOC and model skill does not imply that the contribution of the ox-VOC to OA and particle number concentration is unimportant. Rather, it implies that within the current modelling framework its contribution can be compensated for by changes in other ox-VOCs.

To identify the plausible and implausible parts of the parameter space using the patterns of dependencies, the ensemble simulations (denoted by triangles in Figs. 10, 11, and 12) in the subplots are shaded blue to red. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent high/good Taylor skill score within the ensemble. We note the relative rank of the simulations in Taylor skill score and their relative positions in each 2-D subplot and use this information to identify clusters of blue or red triangles in the parameter space. For absolute values of Taylor skill scores of each simulation, see Table A1 and Figs. A3, A4, and A5 (which are Figs. 10, 11, and 12 labelled with simulation number).

Some clear patterns in model skill are apparent across the six dimensions. N3 skill depends strongly on B_ELVOC, B_LVOC, and A_LVOC (Fig. 10). The skill is generally lower (clusters of blue triangles) for B_ELVOC yields less than 6.4 % from ozonolysis of α-pinene (i.e. twice the baseline yield; see Fig. 10 first column, fifth row), irrespective of the value of other parameters (left column in Fig. 10). The N3 skill is also generally low for values of A_LVOC greater than about 95 Tg yr⁻¹ (yield corresponding to 150 % in Fig. 10), although there are a few simulations that have reasonable skill (triangles in lighter shades of red amidst mostly blue triangles) above this value (Fig. 10 second row). We note two additional regions in the 6-D space dominated by low model skill in N3 – where B_ELVOC yield is greater than 19.8 % and B_LVOC is less than 113 Tg yr⁻¹ (Fig. 10 first column, fifth row, cluster of blue triangles in the bottom right corner of the subplot corresponding to 6 times the baseline yield on the x axis and a yield of 100 % on the y axis) and where the sum of anthropogenic LVOC and SVOC is greater than 127 Tg yr⁻¹ (Fig. 10 fifth column, first row, cluster of blue triangles in the top right corner of the subplot corresponding to 200 % yield on both x and y axes). There is also a general increase in skill for high values of B_LVOC (second column).

N50 skill has the strongest dependency on B_LVOC. The model is most skilful (clusters of red triangles) for B_LVOC greater than about 113 Tg yr⁻¹ (corresponding to a yield of 100 %, Fig. 11 second column) with a general increase in skill for higher values of B_LVOC. N50 skill is generally lower (blue triangles) for B_ELVOC yields less than twice the baseline yield or 6.4 %, although other parameter values and particularly high B_LVOC improve model skill in some cases (Fig. 11 first column, fifth row; note the cluster of red triangles in the top left corner of the subplot compared to the same area in the subplots above). The dependence of N50 skill on A_LVOC is much weaker than for N3, with high and low model skills spread across the entire parameter range. In contrast the N50 skill tends to be low for A_SVOCs greater than about 95 Tg yr⁻¹ (150 % yield, Fig. 11 first row).

Figures 10 and 11 show that model simulation of N3 and N50 is most skilful when B_ELVOC production is a factor of 2 to 8 higher than the baseline model B_ELVOC yields of 3.2 % and 1.2 % from O₃ and OH• oxidation reactions respectively (Kirkby et al., 2016). With less than 6.4 % B_ELVOC yields from O₃ and OH• oxidation reactions, the model–observation match is consistently poor (shades of blue in Figs. 10 or 13). The best estimate of B_ELVOC yield to obtain reasonable agreement with observed N3, N50, and OA in our model (denoted by shades of red in Fig. 13) is about a factor of 4 higher than the baseline B_ELVOC yield from ozonolysis of α-pinene (Kirkby et al., 2016).

OA skill has the strongest joint dependency on B_LVOC and B_SVOC_M (Fig. 12). The joint distribution suggests that the skill is poor if the sum of these SOA production rates exceeds about 226 Tg yr⁻¹ (up to 200 % yield of both; Fig. 12 second column, fourth row). However, for all other parameters there are skilful and unskillful simulations right across the 6-D parameter space.

Together, the variations in skill for N3, N50, and OA concentrations across the six dimensions show that the parameter space for a high N3 or N50 skill score does not overlap with the parameter space for a high OA skill score (Figs. 10, 11, and 12). This gives an insight as to why models that are fine-tuned to simulate particle number concentrations underestimate OA mass concentrations in the atmosphere – also identified as a persistent challenge for state-of-the-art global models (Kanakidou et al., 2005; Spracklen et al., 2011).

Our results reveal the problem of model equifinality, highlighted for the whole aerosol model by Lee et al. (2016). Equifinality means that there are multiple ways (i.e. parameter combinations) of achieving the same model skill against observations, which makes it difficult to identify the best model (Beven, 2006). An important consequence of equifinality is that fine-tuning any one aspect of the model (for example, the nucleation mechanism) to achieve best model–observation agreement for any one variable (e.g. the particle number concentration) can be achieved with a wide range of settings of other parameters (e.g. parameters controlling overall OA mass production). While this may not affect the overall model skill in the particular evaluation, the various parts of equally plausible parameter space may result in very different model behaviour in, say, climate projections.

Figure 10 Taylor skill score for model simulations against N3 observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (%) over which each of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation, and the color of the triangle indicates its Taylor skill score for N3. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent high/good Taylor skill score. Figure A3 shows the same plot with ensemble members numbered. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Figure 11 Taylor skill score for model simulations against N50 observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (%) over which each of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation and the color of the triangle indicates its Taylor skill score for N50. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent high/good Taylor skill score. Figure A4 shows the same plot with ensemble members numbered. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Figure 12 Taylor skill score for model simulations against OA observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (%) over which each of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation and the color of the triangle indicates its Taylor skill score for OA. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent a high/good Taylor skill score. Figure A5 shows the same plot with ensemble members numbered. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Figure 13 summarizes Figs. 10, 11, and 12 showing the model skill score in all three model outputs across the entire parameter space for all six ox-VOCs in 1-D. We find five simulations that are shaded red for all three model outputs across the six parameters: simulations 8, 17, 19, 21, and 41. The parameter combinations and the resulting Taylor skill scores are listed in Tables 3 and A2.

Table 3 Yield of B_ELVOC and SOA (Tg yr⁻¹) for other ox-VOCs with the corresponding Taylor skill scores for five ensemble members that are shaded red for all three outputs N3, N50, and OA in Fig. 13. The parameter combinations used to produce the above ox-VOCs in the ensemble are listed in Table A2.

One ensemble member, simulation 41, scores reasonably well (Taylor skill scores of 0.28, 0.11, and 0.14 for N3, N50, and OA) in simulating observed mass and number of particles. It is the only simulation for which Taylor skill score in each of the three outputs is among the top 10 highest scores within the ensemble (see Table A1). The score is highest for N3, second highest for N50 (0.12 being the highest score), and sixth highest for OA (0.19 being the highest score). For this simulation B_ELVOC is 4.1 times the baseline yield (i.e. about 13 % B_ELVOC yield from O₃ and 5 % yield from OH•), with about 157 Tg yr⁻¹ of SOA from B_LVOC, 33 Tg yr⁻¹ from B_SVOC (monoterpenes + isoprene), 38 Tg yr⁻¹ from A_LVOC, and 13 Tg yr⁻¹ from A_SVOC. In this simulation the SOA production pathways are characterized by (i) high concentrations of both B_ELVOC and B_LVOC that ensure particle production via nucleation and subsequent growth of nucleated clusters, (ii) relatively high A_LVOC concentrations that further help to sustain the growth of small particles in the nucleation mode, and (iii) modest yields of SVOCs that ensure effective growth of particles up to sizes relevant to cloud condensation nuclei (CCN), improving skill scores for N50 and OA, while at the same time restraining the condensation sink and the loss of too many growing particles by coagulation scavenging. With Taylor skill scores of 0.28, 0.11, and 0.14 for N3, N50, and OA, simulation 41 has much scope for improvement. Nevertheless it exemplifies the characteristics required to improve the simulation of SOA in the model.

Figure 13 The 1-D representation of the 6-D space-filling design of experiments shaded by model performance. The 18 subplots (labelled in plot) correspond to the design of experiments for six ox-VOCs, each shaded according to the Taylor skill score (TSS) in three model outputs (N3, N50, OA – in that order for each ox-VOC). Each subplot shows the parameter range perturbed for the yield (%) of the corresponding ox-VOC* (x axis). Each point on a subplot represents a simulation, and the colour of the point indicates the performance of the simulation (TSS) against observed N3, N50, or OA. Darker shades of blue indicate low/poor TSS and darker shades of red represent high/good TSS. The plot identifies plausible and implausible parameter space for each ox-VOC. Simulations 8, 17, 19, 21, and 41 that are shaded red for all three model outputs across the parameter space of all six ox-VOCs are labelled. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

4 Conclusions

We have used a perturbed parameter ensemble of 60 model simulations to explore how uncertainty in six biogenic and anthropogenic precursors affects organic aerosol mass and particle number concentrations. The ranges for each parameter were chosen to encompass maximum uncertainty associated with organic compounds that affect three different stages of SOA formation − nucleation, cluster growth, and particle growth. Simultaneous perturbations of the six parameters using a Latin hypercube sampling technique allow the effects of parameter combinations rather than just individual parameters on model outputs to be explored. Three model outputs, the number concentration of particles larger than 3 nm diameter (N3), the number concentration of particles larger than 50 nm diameter (N50), and the organic aerosol (OA) mass concentration, were compared against observations, and the model skill score was then used to determine the skilful parts of parameter space.

The results expose a high degree of equifinality in the SOA model in which there are multiple ways of generating similar outputs (particle concentrations and OA mass). This is to be expected in a system with six free parameters and only three output variables of interest – that is, our six-component SOA model is underdetermined. Equifinality, or compensating parameter effects, limit the extent to which the best set of parameters can be identified by comparing the simulations against observations. Our results suggest that the effects of three categories of volatile organic compounds can be detected by comparing the 60 simulations against observations: ELVOC, LVOC, and SVOC. B_ELVOC is crucial for the formation of particles via nucleation. Thereafter contributions from LVOCs and SVOCs contribute to the growth of freshly-nucleated particles to produce a realistic N50 concentration and SOA mass.

B_ELVOC strongly influences model skill scores in N3 and to a lesser extent in N50 Fig. 13). When B_ELVOC is low (< twice the baseline yield of 3.2 %), the ensemble consistently underestimates N3 and N50 number concentrations, irrespective of the availability of other ox-VOCs. We find the best model skill scores in N3, N50, and OA are achieved when the ELVOC yield from precursor VOCs is between 6–26 %, with the most plausible ELVOC yield estimate being around 12.8 %. Previously reported ELVOC yields from α-pinene ozonolysis are at the lower end of this constrained range (3.2 % with an uncertainty range of +100 %$/-$60 % reported by Kirkby et al. (2016), 7 ± 4 % reported by Ehn et al. (2014), and 3.4 ± 1.7 % reported by Jokinen et al. (2015)). Ehn et al. (2014) and Jokinen et al. (2015) report higher yields of ELVOCs from other common biogenic VOCs, and nucleating ELVOCs from anthropogenic sources have been reported by Molteni et al. (2018). Our range defines the plausible boundary for modelling nucleating ELVOCs by eliminating parameter space that results in low model skill scores. The best estimate for nucleating ELVOCs for the best model skill score within our plausible range is determined by the parameter combinations of other VOCs.

B_LVOC has the strongest influence on model skills in N3, N50, and OA (second column, Figs. 10, 11, and 12). With insufficient B_LVOC, nucleated particles are lost before they can reach climate-relevant sizes in the model (Fig. 6). High B_LVOCs compensate to some extent for low B_ELVOCs, depending on the value of other ox-VOCs (see Sect. 3.1). Both N50 and OA skill scores show the strongest relationship with B_LVOC, but the parts of parameter space that favour high N50 skill and high OA skill contradict each other (second column, Figs. 11 and 12). Most skilful predictions of N50 tend to require over 113 Tg yr⁻¹ of B_LVOC, while the most skilful predictions of OA mass are favoured by B_LVOC lower than 113 Tg yr⁻¹. Our results show this inconsistency may be reconciled when B_LVOC greater than 113 Tg yr⁻¹ (favouring good N3 and N50 skill score) is accompanied by low B_SVOC_M such that the sum of B_LVOC + B_SVOC_M does not exceed about 226 Tg yr⁻¹ (favouring a good OA skill score, due to the joint dependency of OA skill score on these two ox-VOCs − see Fig. 12, cluster of red triangles towards the bottom right corner of subplot in second column, fourth row).

SVOCs (from monoterpene, isoprene, and anthropogenic sources) are important to simulate realistic number concentrations of climate-relevant sized particles such as N50 and OA mass concentrations. SVOCs generated from α-pinene in the model are spatially co-located with nucleated particles and hence have a strong effect on their growth to larger sizes. OA skills in particular show a strong dependence on B_SVOC_M as discussed above.

We cannot determine the plausible or implausible parameter space for those ox-VOC parameters which show a weak relationship with the model skill score. Model skill in N3, N50, and OA have the weakest relationship with B_SVOC_I. This is because in the current model setup the role of B_SVOC_I in growing particles by mass-based partitioning can be easily compensated for by other oxidized VOCs. The relationship between anthropogenic oxidized VOCs and model skill score in OA is also weak, although A_LVOC shows a strong relationship with N3 skill score and A_SVOC with N50 skill score. The weak relationship between anthropogenic SVOCs with simulated OA mass is because the clusters produced by biogenic nucleation are not spatially co-located with the anthropogenic LVOCs or SVOCs. We expect the relationship between anthropogenic VOCs and model skill scores in N3, N50, and OA to strengthen when the role of anthropogenic oxidized VOCs in nucleation and cluster growth (Molteni et al., 2018) are represented. Our results point to a structural deficiency in the model. The perturbed parameter ensemble tends to exaggerate the observed seasonal cycle of particle concentrations, overestimating summer and underestimating wintertime particle concentrations. We suggest this is due to SOA sources in models being predominantly biogenic. Increases in the concentration of nucleating biogenic aerosols or in the complexities of nucleation mechanisms will not improve the model's ability to replicate the observed seasonal cycle. We expect the inclusion of anthropogenic sources in nucleation and in initial cluster growth for SOA particles in models would significantly improve the simulated aerosol seasonal cycle. It is important to explore the anthropogenic contribution to cluster-growth because such pathways could considerably impact model estimates of anthropogenic aerosol forcing by differentially affecting the present-day and pre-industrial atmospheres (Carslaw et al., 2013).

Appendix A

In Fig. 8 several simulations overestimate summertime N3 concentrations, showing the issue is not so much as that the underprediction in winter but the simulation of a largely exaggerated seasonality in total aerosol number concentrations. To validate the above, two simulations from the ensemble (Fig. 8) – simulation 55 with the lowest global mean N3 and simulation 54 with the highest global mean N3 – are shown in Fig. A1. Figure A1 clearly shows with increasing B_ELVOCs (which will have a strong effect on N3) the model increasingly overestimates particle concentrations in summer with little change in particle concentrations in winter, thereby reducing the model–observation (negative) bias but worsening the model–observation correlation coefficient Fig. A2).

Figure A1 Annual cycle of simulated (solid lines) and observed (black stars) monthly mean surface-level N3 concentrations at 34 ground-based sites. The ensemble members shown are those with the lowest and highest number concentrations of global mean surface-level N3 from Fig. 8.

Figure A2 Correlation coefficient (diamonds in blue for winter, yellow for spring, red for summer, and green for autumn) and and normalised mean bias factor (squares in blue for winter, yellow for spring, red for summer, and green for autumn) shown against each perturbed parameter ensemble member (PPEM). The simulations are arranged in order of increasing global mean surface-level N3 (as in Fig. 5) with every sixth simulation labelled on the x axis. The order of simulations on the x axis: 55, 53, 23, 58, 16, 42, 44, 33, 40, 46, 34, 52, 35, 7, 13, 24, 2, 12, 17, 60, 10, 37, 14, 4, 1, 26, 3, 28, 39, 5, 25, 31, 19, 18, 30, 46, 29, 47, 11, 56, 50, 38, 22, 43, 9, 49, 45, 41, 32, 36, 6, 21, 8, 20, 15, 57, 59, 27, 51, 54.

Figure A3 Taylor skill score for model simulations against N3 observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (%) over which of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation (labelled 1–60), and the color of the triangle indicates its Taylor skill score for N3. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent high/good Taylor skill score. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Figure A4 Taylor skill score for model simulations against N50 observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (%) over which of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation (labelled 1–60) and the color of the triangle indicates its Taylor skill score for N50. Darker shades of blue indicate low/poor Taylor skill score and darker shades of red represent high/good Taylor skill score. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Figure A5 Taylor skill score for model simulations against OA observations across the 6-D parameter space. The x and y axes for a subplot show the total range of reaction yields (in percent) over which of the two parameters (as indicated by the plot labels at the top and right for each subplot respectively) is perturbed in the ensemble. Each triangle in a subplot represents a simulation (labelled 1–60) and the color of the triangle indicates its Taylor skill score for OA. Darker shades of blue indicate low/poor Taylor skill score, and darker shades of red represent high/good Taylor skill score. Note the axis for B_ELVOC shows the scaling factor for B_ELVOC yields. Axes for the rest show the corresponding ox-VOC yields.

Table A1 Summary statistics – Taylor skill score (TSS), normalized mean bias factor (NMBF), and Pearson correlation coefficient (R) – for each ensemble member, based on model comparison against observations of N3, N50, and OA. The 10 best and 10 worst simulations in each category are highlighted in bold and designated with (+) and (−) respectively.

Table A2 Scaling factor for B_ELVOC and yields of other ox-VOCs with the corresponding Taylor skill scores for five ensemble members that are shaded red for all three outputs N3, N50, and OA in Fig. 13.