Semivolatile POA and parameterized total combustion SOA in CMAQv5.2: impacts on source strength and partitioning

Abstract

Mounting evidence from field and laboratory observations coupled with atmospheric model analyses shows that primary combustion emissions of organic compounds dynamically partition between the vapor and particulate phases, especially as near-source emissions dilute and cool to ambient conditions. The most recent version of the Community Multiscale Air Quality model version 5.2 (CMAQv5.2) accounts for the semivolatile partitioning and gas-phase aging of these primary organic aerosol (POA) compounds consistent with experimentally derived parameterizations. We also include a new surrogate species, potential secondary organic aerosol from combustion emissions (pcSOA), which provides a representation of the secondary organic aerosol (SOA) from anthropogenic combustion sources that could be missing from current chemical transport model predictions. The reasons for this missing mass likely include the following: (1) unspeciated semivolatile and intermediate volatility organic compound (SVOC and IVOC, respectively) emissions missing from current inventories, (2) multigenerational aging of organic vapor products from known SOA precursors (e.g., toluene, alkanes), (3) underestimation of SOA yields due to vapor wall losses in smog chamber experiments, and (4) reversible organic compounds–water interactions and/or aqueous-phase processing of known organic vapor emissions. CMAQ predicts the spatially averaged contribution of pcSOA to OA surface concentrations in the continental United States to be 38.6 and 23.6 % in the 2011 winter and summer, respectively.

Whereas many past modeling studies focused on a particular measurement campaign, season, location, or model configuration, we endeavor to evaluate the model and important uncertain parameters with a comprehensive set of United States-based model runs using multiple horizontal scales (4 and 12 km), gas-phase chemical mechanisms, and seasons and years. The model with representation of semivolatile POA improves predictions of hourly OA observations over the traditional nonvolatile model at sites during field campaigns in southern California (CalNex, May–June 2010), northern California (CARES, June 2010), the southeast US (SOAS, June 2013; SEARCH, January and July, 2011). Model improvements manifest better correlations (e.g., the correlation coefficient at Pasadena at night increases from 0.38 to 0.62) and reductions in underprediction during the photochemically active afternoon period (e.g., bias at Pasadena from $-$5.62 to$-$2.42 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$). Daily averaged predictions of observations at routine-monitoring networks from simulations over the continental US (CONUS) in 2011 show modest improvement during winter, with mean biases reducing from 1.14 to 0.73 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$, but less change in the summer when the decreases from POA evaporation were similar to the magnitude of added SOA mass. Because the model-performance improvement realized by including the relatively simple pcSOA approach is similar to that of more-complicated parameterizations of OA formation and aging, we recommend caution when applying these more-complicated approaches as they currently rely on numerous uncertain parameters.

The pcSOA parameters optimized for performance at the southern and northern California sites lead to higher OA formation than is observed in the CONUS evaluation. This may be due to any of the following: variations in real pcSOA in different regions or time periods, too-high concentrations of other OA sources in the model that are important over the larger domain, or other model issues such as loss processes. This discrepancy is likely regionally and temporally dependent and driven by interferences from factors like varying emissions and chemical regimes.

Introduction

A substantial fraction of atmospheric particles are made of organic compounds . This is true over the continental United States, where pollutants introduced by human sources often interact in complex ways with compounds from natural sources . In order to fully describe the impacts that human activities have on public health, regional haze, and climate change via airborne particles, we must be able to quantify the capacity for relevant organic compounds to form and maintain airborne particulate mass . To achieve this end, atmospheric models routinely rely on estimates of volatility (i.e., vapor pressure or saturation concentration) and solubility in a predominantly aqueous phase. Although research efforts demonstrate clearly the importance of organic solubility in water on phase partitioning and particle viscosity , improvements to the conceptual model of organic compounds partitioning into an organic-rich particulate phase alone are useful before the entire system is treated holistically.

It is clear from available experimental data that primary organic aerosol (POA), which is operationally defined as the population of organic compounds that are emitted in the particulate phase, is a more dynamic quantity than originally prescribed in chemical transport models . Not only do these studies show that POA evaporates upon dilution of exhaust from multiple sources and that it is often more volatile than secondary organic aerosol (SOA; which had been treated as semivolatile in models), but the semivolatile vapors that are produced are most likely susceptible to continued oxidation in the gas phase. first demonstrated, by applying a chemical transport model to predict air quality over the eastern US (EUS) in July 2001, that treating primary particulate compounds as semivolatile improves model performance for bulk organic aerosol (OA) mass predictions and better predicts the relative importance of POA and SOA, which is formed from the oxidation products of volatile organic compound (VOC) gases. Moreover, predicted that organic compounds with intermediate volatility (IVOCs), that is, volatility higher than that of semivolatile organic compounds (SVOCs) and lower than traditionally defined VOCs, were emitted in significant quantities and that their oxidation products could condense and form OA downwind of combustion sources. An IVOC mass emission scale factor of 1.5 times the POA $+$ SVOC emission rate led to significant contribution of OA derived from IVOCs when applied in their model. However, that scale factor was uncertain, as was the overall efficiency with which those IVOCs formed SOA .

More recently, both experimental and model studies have sought to constrain combustion-related SVOC and IVOC emissions as well as their subsequent SOA yields . Progress has been made characterizing the overall significance of this OA production pathway, but the variability among vehicles and other combustion sources is uncertain, as are the effects of operating (e.g., load, burn temperature) and environmental conditions. Additional limitations exist in our understanding of other important processes for SOA formation. First, multigenerational aging of anthropogenic and biogenic condensable organic vapors could affect OA partitioning if chamber experiments do not proceed for a long enough duration to fully capture the processes, and these processes lead to a net increase or decrease in the volatility of the population . Second, it is likely that previous estimates of SOA yields from well-known VOC precursors (e.g., toluene, xylene, alkanes, isoprene, terpenes) derived from smog chamber experiments were biased by unaccounted loss of semivolatile organic vapors to the chamber walls . The magnitude of this bias likely varies substantially among experimental systems (i.e., precursor and oxidant identity) and experimental conditions e.g., seed concentration, relative humidity, temperature;. Finally, aerosol liquid water may absorb soluble organic vapors and shift the partitioning to enhance OA concentrations as explored by .

Several models have now incorporated the semivolatile POA partitioning phenomena demonstrated by , relying on a framework popularly referred to as the volatility basis set VBS;. Examples of these models include box-scale , regional-scale , and global-scale applications. Each implementation has made assumptions about the emissions and SOA yields of IVOCs, the reaction rates and products of multigenerational aging, and, in some cases, the biases from chamber wall losses of semivolatile vapors. Despite wide application of the VBS framework, large-scale modeling studies, which include many other uncertainties from emissions inventories, scavenging processes, and oxidant concentrations have made slow progress towards reducing the uncertainty of the many free parameters related to OA emission, formation, and processing. However, the details of the implementation in CTMs affect OA source apportionment . Current implementations in global models result in a spread of 1–2 orders of magnitude in OA concentrations across models, especially aloft and in remote areas . Therefore, caution should be exercised when adopting an SOA formation mechanism that appears detailed but is largely unconstrained.

showed that an optimized two-parameter fit for SOA formation from anthropogenic sources (constraining the amount and timescale of SOA formation) could successfully capture trends observed in multiple field studies (DeCarlo et al., 2010). That parameterization reasonably reproduced aerosol mass spectrometer (AMS) observations within and downwind of a megacity. The SIMPLifiEd parameterization of combustion SOA (SIMPLE) demonstrated by involved one VOC surrogate which, when oxidized by hydroxyl radicals, formed one condensable vapor that partitioned irreversibly to the particle phase. The use of nonvolatile SOA as an approximation was based on the observation that real SOA had low volatility and generally much lower than used in prior models . The uncertain parameters, the emission of the VOC and its oxidation reaction rate, were optimized for performance in the regional model CHIMERE. This method was reapplied to analyze urban-scale data at Pasadena . Most recently, applied the SIMPLE method in the Community Multiscale Air Quality model version 5.2 (CMAQv5.2) to Pasadena and showed substantially improved model–measurement agreement, while applied it to the SE US and showed agreement with fossil/non-fossil carbon observations.

We document here the implementation and evaluation of semivolatile partitioning of POA and addition of observed, but unspeciated, SOA in CMAQv5.2. The model approach uses the VBS framework to account for the dynamic partitioning and aging of POA emissions, without modifying the two-product approach for SOA formation from traditional SOA precursor VOCs (e.g., toluene, xylene, alkanes, isoprene, terpenes). We further introduce a new surrogate species, potential SOA from combustion emissions (pcSOA) to account for missing mass from IVOC oxidation, multigenerational aging of (anthropogenic) secondary organic vapors (from IVOC and VOC precursors), biases in SOA yields from vapor wall losses, and enhanced organic partitioning to the condensed aqueous phase. In addition to these sources, pcSOA could account for mass from oxidation of as-yet unidentified sources of SOA precursors. The contribution of this surrogate species to the total OA burden is governed by two parameters, analogous to the simple method of – the emission rate of the SOA precursor and the precursor reaction rate with hydroxyl radicals. In the following sections we present a comprehensive evaluation of the improvement in CMAQ OA predictions at two horizontal resolutions (4 and 12 km) over multiple seasons and years, using multiple gas-phase chemical mechanisms. We explore the sensitivity of CMAQ OA to two parameters directly impacting pcSOA formation and recommend important areas of future work.

Organic aerosol module configuration

In developing CMAQv5.2, we have updated the OA module to be more consistent with current understanding of POA emission and to better represent the magnitude of SOA formation from anthropogenic combustion sources.

Properties of semivolatile POA species in CMAQv5.2.

Species name	Molecular weight	$C^{\ast ,\mathrm{1}}$	N${}_{\mathrm{C}}$	O : C	H${}_{\mathrm{vap}}$	POA emission
(Particle/vapor)	$\mathrm{g}$ mol${}^{-\mathrm{1}}$	$\mathrm{\mu }$g m${}^{-\mathrm{3}}$			$\mathrm{kJ}$ mol${}^{-\mathrm{1}}$	Fraction${}^{\mathrm{2}}$
ALVPO1/VLVPO1${}^{\mathrm{3}}$	218	0.1	13.0	0.185	89	0.09
ASVPO1/VSVPO1${}^{\mathrm{4}}$	230	1	14.5	0.123	85	0.09
ASVPO2/VSVPO2	241	10	16.0	0.073	81	0.14
ASVPO3/VSVPO3	253	100	17.5	0.032	77	0.18
AIVPO1/VIVPO1${}^{\mathrm{5}}$	266	1000	19.0	0.00	73	0.5
ALVOO1/VLVOO1${}^{\mathrm{6}}$	136	0.01	5.0	0.886	93	n/a
ALVOO2/VLVOO2	136	0.1	5.5	0.711	89	n/a
ASVOO1/VSVOO1	135	1	6.0	0.567	85	n/a
ASVOO2/VSVOO2	135	10	6.5	0.447	81	n/a
ASVOO3/VSVOO3	134	100	7.0	0.345	77	n/a

${}^{\mathrm{1}}$ $C^{\ast }$ values are defined at reference temperature 298 K. ${}^{\mathrm{2}}$ . ${}^{\mathrm{3}}$ A $=$ aerosol; V $=$ vapor; LV $=$ low volatility; PO $=$ primary organic; ${}^{\mathrm{4}}$ SV $=$ semivolatile; ${}^{\mathrm{5}}$ IV $=$ intermediate volatility; ${}^{\mathrm{6}}$ OO $=$ oxidized organic. n/a: not applicable.

POA semivolatile partitioning and aging

The model accounts explicitly for the gas–particle partitioning of primary organic compounds and their multigenerational oxidation products with a range of volatility, similar to the 1.5D volatility basis set scheme of . Table  lists the new species in the aerosol module as well as important properties. Like , we choose surrogates with volatilities that are relevant for typical atmospheric loadings (i.e., 0.1–1000 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$). related the volatility, expressed in terms of an effective saturation concentration $C^{\ast }$, to both carbon number and O : C using structure activity relationships, vapor pressure estimation methods like SIMPOL , and available vapor pressure measurements, and we apply their methods here (see Supplement). For the directly emitted species (LVPO1, SVPO1, SVPO2, SVPO3, and IVPO1), the carbon number (C17-C19) and O : C (0–0.1) of the new surrogate species are chosen to be consistent with laboratory and field observations of the properties of POA from combustion sources . Previously, used a separate tracer, PNCOM, in CMAQ to track the contribution of non-carbon organic matter to the total POA mixture. During the emission input generation procedure, PNCOM emission factors were informed by the OM : OC (the ratio of total organic mass to organic carbon mass) of individual sources (e.g., mobile, 1.25; cooking, 1.4; biomass burning, 1.7). In the present model we sum the primary organic carbon and non-carbon emission factors in order to preserve the source-aware mass emissions developed by . Because the new primary surrogate species have fixed O : C (Table ), the source resolution in O : C is lost. This limitation can be overcome in future model versions by treating, for example, biomass burning OA independently with additional model species.

For the oxidation products of the directly emitted species (LVOO1, LVOO2, SVOO1, SVOO2, SVOO3), we have chosen higher OM : OC (and thus lower carbon number) values than to bound the typical observed range of OM : OC reported from ambient studies, including and for multiple US and global sites (1.3–2.25) and for routine-monitoring networks in the continental US (0.79–2.15). The molecular weights of these species are calculated using the given carbon number and OM : OC, while assuming each representative species comprises only carbon, oxygen, and hydrogen atoms. The H : C is calculated from O : C using insights from – a common assumption for deriving molecular weights of VBS species in chemical transport models. The species enthalpy of vaporization is assumed to be a linear function of $C^{\ast }$ , and the Henry's law coefficient is prescribed at 2 $\times$ ${\mathrm{10}}^{\mathrm{8}}$ $\mathrm{M}{\mathrm{atm}}^{-\mathrm{1}}$, which is in the range reported for oxidized organic vapors by and .

We employ one volatility distribution, from , for the partitioning calculation of the primary organic emissions from all combustion sources. found some differences in their results when they used separate distributions for biomass burning, gasoline, and diesel sources. The distribution results in more aggressive POA evaporation than the source-resolved multi-distribution parameterization. showed that the parameterization resulted in POA volatility that was similar to thermodenuder observations. We performed a sensitivity simulation using a single modified distribution calculated by summing the source-resolved distributions weighted by their total emissions in the California (CAL) domain. This simulation resulted in maximum hourly POA increases of 10–15 %. Future work will incorporate emerging measurements of primary organic compound volatility information directly into the generation of emission inventories, along with estimating sensitivities to key factors like ambient temperature, operating conditions, and fuel type. We do not apply a scaling factor to the POA input emission factor to, for example, introduce new SVOCs to the model. Our approach assumes that existing inventory POA emission factors are actually representative of the total gas- plus particle-phase mass with saturation concentration below about 3200 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$. It is unclear to what extent SVOCs are missing from the emission factors that inform current inventories as these experiments were conducted at low dilution ratios and enhanced partitioning to the particulate phase. This issue should be resolved by taking into account temperature and organic aerosol loading in past and future laboratory-scale experiments when that data are available.

Additionally, some of the semivolatile primary particle mass that evaporates condenses back to the particle phase after oxidation. Although direct observations of primary combustion-derived SVOC aging are rare, both field results and their relatively large carbon numbers support that they are susceptible to OH oxidation and SOA formation in the atmosphere. We apply an OH oxidation reaction with a rate constant equal to 4.0 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$ cm${}^{\mathrm{3}}$ molec${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ to simulate this aging, consistent with previous modeling studies . The stoichiometric coefficients for this oxidation step are derived using 2D-VBS theory for low-NO${}_x$ systems and are documented in supporting information. The branching ratio for functionalization and fragmentation processes is parameterized as in Eq. (). $${\displaystyle }{\mathit{\beta }}_{\mathrm{frag}}=(\mathrm{O}:\mathrm{C}{)}^{\mathrm{0.4}}$$ The aging configuration tends to produce lower volatility (thus more particulate) mass than that used by and produces slightly less total organic mass (OM) than the full 2D-VBS model of (see Fig. S1a in the Supplement). The bulk O : C enhancement due to aging is more vigorous in this model compared to the full 2D VBS. However, the rate of O : C enhancement is about half that of the POA aging scheme implemented in older CMAQ versions (see Fig. S1b).

SOA from traditional precursors

SOA formed from oxidation of traditional VOC sources (isoprene, monoterpenes, sesquiterpenes, benzene, toluene, xylene, alkanes, and PAHs) are documented in and . The SOA yields of these compounds has been unchanged in the development of CMAQv5.2, although updates to their molecular properties (solubility, molecular weight, OM : OC, etc.) were incorporated to maintain consistency throughout the model processes . CMAQv5.2 also includes aqueous processing of isoprene epoxides, glyoxal, and methylglyoxal, as well as production of oligomeric species from particle-phase reaction of traditional SOA compounds. We do not employ the recent liquid–liquid phase partitioning algorithms implemented by .

Properties of potential SOA from combustion emissions (pcSOA) species and its precursor, pcVOC, for the base case (BASE) CMAQ simulation.

Property	Value
Molecular weight	170 g mol$-\mathrm{1}$
Emissions scale factor (pcVOC/POA)	0.0568 mol g$-\mathrm{1}$
OH oxidation rate constant ($k_{\mathrm{OH}}$)	1.25 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$ ${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$
Condensable product saturation concentration ($C^{\ast }$)	${\mathrm{10}}^{-\mathrm{5}}$ $\mathrm{\mu }$g m${}^{-\mathrm{3}}$

Potential SOA from combustion emissions

Because the individual contributions to the total bias from missing IVOC oxidation, missing multigenerational aging of VOC oxidation products, under-representing SOA yields because of chamber wall losses, and uptake of organics by aerosol liquid water are all uncertain, we introduce one surrogate aerosol species, pcSOA, to address the total missing OA mass. Three essential features are common among the aforementioned uncertain processes: they involve emission of gas-phase compounds, some degree of photooxidation of those compounds, and condensation of the resulting products. Thus, we simulate this formation process analogously to the SIMPLE approach of . A new surrogate VOC species (potential VOC from combustion emissions, pcVOC) is introduced with an emission rate that is scaled to the POA mass emission rate. This species does not partition directly to the particle phase and, in that respect, behaves as a VOC in the model. It is oxidized with OH to form a low-volatility condensable vapor – potential secondary organic gas from combustion emissions, pcSOG (Table ). For the BASE (base emission, base reaction) simulation we assume an OH reaction rate constant equal to the optimal rate constant reported by , and we convert the optimal emission scale factor of , which is scaled from CO emissions to one that can be applied to POA emissions. To convert from ${\mathrm{g}}^{-\mathrm{1}}$ CO to ${\mathrm{g}}^{-\mathrm{1}}$ POA, we multiply by the average ratio of CO emissions to POA emissions in the CONUS11 (Table 3) emissions dataset, 0.82. A promising alternative to scaling emissions to POA emissions could be to apply an SOA yield directly to the oxidation of total non-methane organic gases (NMOG), as demonstrated in and implemented by for vehicle emissions. Unfortunately, NMOG emissions inventory data, as they currently stand, do not always take into account the evaporation of primary vapors nor are they always treated consistently by typical emissions processing practices. We consider the implementation of semivolatile POA and missing combustion SOA in CMAQv5.2 an important step toward an even more rigorous future treatment, whether that involves scaling emissions to CO or NMOG. As such, we have prioritized minimizing the need for reprocessing emissions inputs for users of previous CMAQ versions, but this does leave room for future refinement of source-specific SOA formation from unspeciated vapor emissions.

The emissions of pcVOC are applied to all combustion-related area and point sources of primary organic carbon except for wildfires (classified as point sources), as it has been shown that fire plumes do not exhibit an SOA enhancement as strong as those dominated by urban emissions . Thus, important sources like on-road vehicles, power plants, and commercial cooking are all included. and showed that residential as well as commercial cooking sources of SOA are potentially an important contribution to non-fossil carbon in Pasadena, and residential cooking sources are missing from the inventories used here. The inclusion of both residential wood burning and agricultural burning area sources in the pcSOA parameterization should be investigated further and potentially treated separately in the future by adding more detailed online source resolution to the CMAQv5.2 OA model. The partitioning between pcSOG and pcSOA is governed by a saturation concentration equal to 1 $\times$ ${\mathrm{10}}^{-\mathrm{5}}$ $\mathrm{\mu }$g m${}^{-\mathrm{3}}$. Although heterogeneous reactions are implemented for other SOA types in CMAQ, no further reactions are included for pcSOA. Additionally, photolysis leading to degradation of low-volatility material is not considered in the model . The low volatility of pcSOA makes this OA configuration behave similarly in many ways to the approach of , which attempts to account for semisolid OA behavior and oligomer formation; however, there are important differences. All OA mass (including pcSOA) is treated as an absorbing medium in our approach and pcSOA is created through direct oxidation and condensation of gas-phase species, not through a particle-phase transformation. Future work will address the role of diffusion-limited OA absorptive partitioning and its dependence on relative humidity. The emission scale factor and oxidation rate constant are uncertain parameters, and we present a sensitivity analysis of them in order to both demonstrate the relevance of this uncertainty to total OA concentrations and to recommend acceptably performing parameters for use in operational CMAQv5.2 simulations.

CMAQ model scenarios used for this study.

Name	Spatial domain	Duration	Gas-phase chemical	Reference
			mechanism
CAL	Southwestern US	3 May–30 Jun 2010	CB05e51	;
EUS	Eastern US	1–30 Jun 2013	SAPRC07tic
CONUS11	Continental US	1 Jan–31 Dec 2011	CB05e51
CONUS02	Continental US	1–30 Jan; 1–31 Jul 2002	CB05e51

Model application

We apply CMAQv5.2 with the enhanced OA module to three distinct simulation domains: the continental United States, the eastern United States, and California. The time periods and chemistry options, detailed below, are chosen to leverage existing observations while also evaluating the model under various configurations. The individual simulations are summarized in Table .

Observations

In order to connect the performance of the updated model with past CMAQ versions, we evaluate predictions against organic carbon (OC) observations from routine-monitoring networks including the Interagency Monitoring of Protected Visual Environments (IMPROVE) and Chemical Speciation Network (CSN) datasets. We note that IMPROVE data had a 27 % low bias relative to collocated SEARCH observations during summer 2013 in the SE US, which is thought to be due to evaporation during sampling and transport . That bias was not observed for non-summer samples. Additionally, we focus on three intensive measurement campaigns: the California Research at the Nexus of Air Quality and Climate Change (CalNex), the Carbonaceous Aerosols and Radiative Effects Study (CARES), and the Southern Oxidant and Aerosol Study (SOAS). At these sites, CMAQ OA predictions are evaluated against High-Resolution Time-of-Flight Aerosol Mass Spectrometer (HR-ToF-AMS) data. The Southeastern Aerosol Research Characterization (SEARCH) network sites located in the southeast US provide semi-continuous OC observations also used to evaluate model capability during the SOAS time period. Two urban sites, Jefferson Street, Atlanta, GA, and Birmingham, AL, and one rural site, Yorkville, GA, are selected. During the CARES and SOAS campaigns, thermodenuder measurements allow us to make a quantitative evaluation of the OA volatility distribution predicted by the model.

The CalNex campaign characterized atmospheric composition at two sites in southern California, Pasadena, and Bakersfield, from 15 May to 29 June 2010 . The Pasadena site was located 18 km northeast and generally downwind of downtown Los Angeles, while the Bakersfield site was situated at the southern end of the San Joaquin Valley (SJV). Several studies characterized OA properties and likely sources at Pasadena and the larger LA basin with observations , while others combined measurements with models at multiple scales . On average, organic compounds comprised 41 % of the total PM${}_{\mathrm{1}}$ mass concentrations at Pasadena and exhibited a strong diurnal-averaged peak concentration in the mid-afternoon 15:00 local time;. At Bakersfield, OA compounds comprised 56 % of the total PM${}_{\mathrm{1}}$ mass, and fossil fuel sources were determined to contribute 80–90 % of those components . Average O : C at Pasadena ranged from about 0.38 to 0.48 throughout the campaign with some extreme values extending to 0.8.

The CARES campaign took place in northern California from 2 to 28 June 2010 and conducted measurements at two ground sites, Sacramento and Cool, and from the DOE G-1 aircraft . The Cool site is a suitable location to observe interactions among anthropogenic and biogenic pollutants due to its proximity to outflow from Sacramento and the forested foothills of the Sierra Nevada mountains . Like those at the CalNex stations, observations during CARES indicated that organic compounds dominated ($\sim$ 80 %) the particle-phase mass composition at both sites . and showed that the combination of anthropogenic and biogenic pollutants in the vicinity of Sacramento could enhance SOA formation significantly beyond instances when just anthropogenic or just biogenic sources contributed.

The SOAS campaign also targeted observations of anthropogenic–biogenic interactions but at three sites (Centreville, AL; Birmingham, AL; and Look Rock, TN) in the substantially differing southeastern US environment during June and July 2013 . Anthropogenic emissions, including NO${}_x$ and ${\mathrm{SO}}_{\mathrm{2}}$, have large influence on isoprene and monoterpene SOA in both rural and urban platforms . and found that organic nitrate chemistry likely played a substantial role in the formation of observed less-oxidized oxygenated OA (LO-OOA) compounds at the Centreville site. also found that SOA from isoprene oxidation contributed 18 % of the total OA averaged throughout the campaign.

Model configuration and analysis

We configured CMAQv5.2 to simulate air quality over California at a horizontal resolution of 4 km by 4 km during May and June 2010. The input parameters are documented extensively by , , and and are summarized briefly here. Meteorological inputs are generated using the Weather Research and Forecasting Model (WRF) Advanced Research WRF core version 3.1 . evaluated the meteorological parameters (temperature, wind speed, wind direction, etc.) and found acceptable performance. The boundary conditions were generated from running a coarser resolution CMAQ simulation from December 2009 through June 2010. That coarse simulation used boundaries driven by a GEOS-Chem (v8-03-02) simulation from the same period . Electrical generating unit and other point source emissions were estimated specifically for 2010 using continuous emissions monitoring data. Mobile source emissions were generated using a combined approach of generation by the SMOKE-MOVES integration platform and projection by the California Air Resources Board. Other anthropogenic emissions are based on the 2011 National Emissions Inventory (NEI) version 1 . Day-specific fire emissions are incorporated into the model inputs, although they do not impact concentrations in Pasadena to a large extent . and also found low impact from biomass burning at both CARES sites. Residential meat-cooking emissions are not included in the emissions inventory as pointed out by . That study found that CMAQ underpredicted cooking-influenced OA by a factor of 2–4, depending on the time of day and assumptions about the volatility of meat-cooking emissions. Biogenic vapor emissions are generated by the Biogenic Emission Inventory System (BEIS) v3.14 . A modified version of Carbon Bond 5 (CB05e51) simulated gas-phase chemical kinetics including ozone formation/destruction and OA formation among other processes.

The SOAS campaign was addressed with a second simulation over the eastern United States during June 2013 at 12 km by 12 km horizontal resolution. The WRF core model version 3.6.1 was applied for meteorological inputs, and boundary conditions were obtained from a 36 km by 36 km CMAQ simulation using boundary conditions informed by GEOS-Chem . As explained by and , anthropogenic emissions were based on the Environmental Protection Agency (EPA) National Emissions Inventory 2011 v1, and biogenic emissions were obtained from BEIS v3.6.1. This implementation included enhanced detail in the formation of isoprene SOA and organic nitrate kinetics and partitioning. Fire emissions were based on the latest version of the Satellite Mapping Automated Reanalysis Tool for Fire Incident Reconciliation system (SMARTFIRE, www.airfire.org/smartfire/). The Statewide Air Pollution Research Center (SAPRC07) mechanism was used to predict gas-phase chemistry.

In order to extend results further, we perform a suite of simulations over the continental United States at 12 km by 12 km horizontal resolution for the entire year 2011 (CONUS11), thereby probing seasonal differences in model performance. Year-to-year variability is assessed with additional runs over the continental US during January and July 2002 (CONUS02). Anthropogenic, biogenic, and day-specific fire emissions for the CONUS simulations are documented by . The CONUS02 meteorological inputs were generated by WRF version 3.1, and anthropogenic emissions were derived from the 2002 NEI . The OA boundary conditions are downscaled from GEOS-Chem simulations, and we prescribe an OM : OC for these compounds equal to 2.0. We estimated organic carbon, which was measured at each site, from the model output by taking into account the species-dependent O : C and calculating an OM : OC using the approach of .

Details of base and sensitivity simulation exploring pcSOA formation parameters.

Simulation	pcVOC emission scale	$k_{\mathrm{OH}}$	Model domains
	factor (mol g${}^{-\mathrm{1}}$)	(${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{molec}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$)
Nonvolatile POA (nvPOA)	n/a	n/a	CAL, EUS, CONUS11, CONUS02
Base emission, base reaction (BASE)	0.0568	1.25 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, EUS, CONUS11 (Jan, Jul)
Low emission, base reaction (LEBR)	0.0387	1.25 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, EUS, CONUS11, CONUS02
High emission, base reaction (HEBR)	0.07	1.25 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, CONUS11 (Jan, Jul)
Base emission, low reaction (BELR)	0.0568	1.0 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, CONUS11 (Jan, Jul)
Base emission, high reaction (BEHR)	0.0568	2.0 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, CONUS11 (Jan, Jul)
Low emission, high reaction (LEHR)	0.0387	2.0 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, CONUS11 (Jan, Jul)
High emission, low reaction (HELR)	0.07	1.0 $\times$ ${\mathrm{10}}^{-\mathrm{11}}$	CAL, CONUS11 (Jan, Jul)

n/a: not applicable.

To probe sensitivity to the uncertain pcVOC emission scale factor and oxidation rate constant, we explore a series of perturbation runs varying both of these parameters (Table ). The parameter combinations are chosen to bound acceptably performing estimates of the parameters at Pasadena . We convert the perturbed pcVOC emissions scale factors from ${\mathrm{g}}^{-\mathrm{1}}$ CO to ${\mathrm{g}}^{-\mathrm{1}}$ POA with the same conversion factor used for the BASE simulation. In addition to these sensitivity cases, we run a case that uses the same OA module configuration as CMAQv5.1 (nvPOA). Specifically, the nvPOA case includes nonvolatile POA aging consistent with and SOA formation from traditional VOC precursors. No pcSOA is implemented for the nvPOA case. We apply all of these cases to the CAL domain and to both January and July of the CONUS11 domain. Meanwhile, due to computation time constraints, we only examine the BASE case, the best-performing CONUS11 case (low emission, base reaction; LEBR), and the nvPOA case for the other months of 2011 and for the EUS and CONUS02 domains.

For semi-quantitative evaluations generally within individual or among similar datasets, we utilize standard statistical metrics to characterize model performance. These metrics include the mean bias (MB; $\mathrm{\mu }$g m${}^{-\mathrm{3}}$), mean normalized error (MNE), and correlation coefficient ($r$). For more rigorous quantitative analysis, we employ the normalized mean bias factor (NMBF, Eq. ) and normalized mean absolute error factor (NMEF, Eq. ) developed by . $${\displaystyle }\mathrm{NMBF}=S\left[\exp ,\left(\left|\ln ,{\displaystyle \frac{\sum M_i}{\sum O_i}}\right|\right)\right]-\mathrm{1},$$ $${\displaystyle }\mathrm{NMEF}={\displaystyle \frac{\sum \left|M_i,-,O_i\right|}{{\left(\sum ,O_i\right)}^{[\mathrm{1}+S]/\mathrm{2}}{\left(\sum ,M_i\right)}^{[\mathrm{1}-S]/\mathrm{2}}}},$$ where $M$ are model predictions, $O$ are observations, and $S=(\stackrel{\mathrm{‾}}{M}-\stackrel{\mathrm{‾}}{O})/|\stackrel{\mathrm{‾}}{M}-\stackrel{\mathrm{‾}}{O}|$. These metrics were specifically designed to address biases that arise from comparing datasets with substantially different mean values (e.g., urban and rural locations) and widely varying degrees of over- and underprediction. If NMBF is positive (e.g., 0.6), then the model overestimates the observations by a factor of 1 $+$ NMBF (1.6), on average. If NMBF is negative (e.g., $-$0.6), the model underestimates the observations by a factor of 1 $-$ NMBF (1.6), on average. The NMEF is interpreted as follows: if NMEF $=$ 1.8, then the model gross error is 1.8 times the mean observation for overprediction or 1.8 times the mean model prediction for underprediction, on average.

Results and discussion

Evaluation against continuous OA mass observations

CMAQv5.2 with semivolatile POA and pcSOA (BASE) successfully captured the day–night OA variability substantially better than the model with nonvolatile POA (nvPOA) (Fig. ). Across all sites, improved predictions were due directly to both the volatilization of primary emissions in the morning and night and to the increased role of secondary formation via pcSOA in the afternoon. We expect these two phenomena to play the largest role in urban locations, like Los Angeles, with high combustion emissions and active photochemistry. At Pasadena, the daytime MB improved from $-$5.62 to $-$2.42 from the nvPOA to the BASE case, consistent with , while the correlation improved from 0.67 to 0.79. The remaining underprediction is qualitatively consistent with , who concluded that current estimates of SOA mass from gasoline and diesel emissions cannot produce the SOA observed in Pasadena if unaccounted-for non-vehicular sources contribute heavily to the pool of SOA precursors. Positive matrix factorization (PMF) analysis of the observed OA into surrogate components revealed the strong influence of semivolatile oxygenated organic aerosol (SV-OOA, dominated by compounds emitted locally) in generating the afternoon OA loadings . The pcSOA approach reproduced this buildup of daytime urban SOA much closer than the nvPOA case.

Organic aerosol concentrations ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$) observed (gray) with the HR-ToF-AMS at sites in California in 2010 (Pasadena, Bakersfield, Sacramento, and Cool). Also shown are the model-predicted distributions at each site using the nonvolatile (pink) and base-case semivolatile (green) configurations. The boxes denote the 25th and 75th percentiles of each dataset, while the whiskers extend to the most extreme points. (a) and (b) show data for daytime (08:00–20:00 LT) and nighttime hours, respectively. Model values are projected to PM${}_{\mathrm{1}}$ to correspond roughly to the size cutoff of the AMS.

[Figure omitted. See PDF]

The BASE case performed similarly to the nvPOA case in Pasadena at night. There, the PMF analysis indicated that the low-volatility oxygenated organic aerosol (LV-OOA) surrogate dominated with roughly constant diurnal-averaged concentrations $\sim$ 2 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$;. LV-OOA had a significant contribution from regional (nonurban) biogenic OA. Other OA components such as hydrocarbon-like OA (HOA), cooking-influenced OA (CIOA), and locally sourced OA (LOA) contributed less throughout the day ($\sim$ 0.5–1 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$), although the concentrations of CIOA were elevated at night. The lack of significant observed OA growth at night–early morning and the lack of observed correlation between aerosol liquid water content and the fraction of condensed water-soluble organic compounds indicated a less important role of nighttime chemistry due to nitrate radicals or aqueous phase uptake .

Aggregate daytime predictions at the Bakersfield site improved dramatically as well (Fig. ). Bulk factor analysis of AMS and Fourier-transform infrared spectroscopy (FTIR) observations indicated five principal PM${}_{\mathrm{1}}$ source categories including aromatic SOA (24 %), alkane SOA (41 %), OA formed at night (10 %), SOA from petroleum operations (14 %), and vegetative detritus (10 %). Multiple studies have characterized the importance of organic nitrates for nighttime SOA composition at the site , while daytime OA formation appears to be driven by oxidation of primary vapors followed by condensation . showed that CMAQv5.0.2 with nonvolatile POA underpredicted PM${}_{\mathrm{2.5}}$ organic carbon concentrations at both Pasadena and Bakersfield, with fractional biases (FB) of $-$53 and $-$144 %, respectively, generally consistent with results from CMAQv5.2 with nonvolatile POA (FB of $-$76 and $-$79 %, respectively).

The BASE case underpredicted concentrations at Bakersfield at night, and this was partially due to missing particle-phase organic nitrate compounds that are not formed in the gas phase by the CB05e51 mechanism used in the CAL simulation (see Fig. S4). Amine compounds contributed about 10 % or less of OA at Bakersfield and Pasadena ; these are missing from the particle-phase representation in this version of CMAQ as well. The day–night trends are consistent with the results of WRF-Chem application to the same campaign . Differences that emerge include a stronger OA daytime peak predicted by CMAQ at Pasadena and higher nighttime concentrations predicted by WRF-Chem at Bakersfield. These differences may have resulted from the use of slightly different anthropogenic emission inventories and selection of different modules for planetary boundary layer dynamics.

At Sacramento, the nvPOA case predicted significantly higher total OA concentrations at night (MB $=$ 0.48 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$) than during the day (MB $=$ $-$1.13 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$), while the BASE model agreed better with observations (nighttime MB $=$ $-$0.35; daytime $=$ $-$0.41 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$). The magnitude of the day–night performance for the rural Cool site was encouraging (MB $=$ $-$0.46 and $-$0.05, respectively), as was the correlation ($r=\mathrm{0.74}$ and 0.63, respectively). isolated three factors using PMF analysis at Cool, including hydrocarbon-like organic aerosol (9 % of total OA), a less-oxidized oxygenated OA (50.3 %) that peaked in the early evening, and a more-oxidized oxygenated OA (MO-OOA, 40.7 %) that stayed relatively constant throughout the day with slightly elevated concentrations at night. The authors concluded that the LO-OOA was correlated with urban transport while the MO-OOA was likely influenced by biogenic VOC oxidation. The fact that about 60 % of the total OA mass was related to urban sources helps explain why the introduction of POA partitioning and pcSOA mass in CMAQ would make a difference to predictions at the Cool site, in the absence of regional fire events. As urban plumes move downwind, pcSOA and POA aging in the BASE case replace most of the evaporated POA and may add extra mass to the total OA burden. For example, the additional anthropogenic OA mass improved model performance at the Cool site, which is often directly downwind of Sacramento, but had less impact on the rural Centreville, AL, site (see Fig. S4).

Observed (gray) and modeled (pink, green) organic aerosol ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$) and organic carbon ($\mathrm{\mu }$gC m${}^{-\mathrm{3}}$) concentrations at sites in the southeast US. OA concentrations at the SOAS sites, Centreville and Look Rock, were measured with HR-ToF-AMS, while OC concentrations at the SEARCH sites, Birmingham, Atlanta, and Yorkville were inferred as the difference between total carbon measured by ambient particulate carbon monitors and elemental carbon measured by aethalometers. Also shown are the model-predicted distributions at each site using the nonvolatile (pink) and base-case semivolatile (green) configurations. The boxes denote the 25th and 75th percentiles of each dataset, while the whiskers extend to the most extreme points. (a) and (b) show data for daytime (08:00–20:00) and nighttime hours, respectively. Model values for the SOAS sites are projected to PM${}_{\mathrm{1}}$ to correspond roughly to the size cutoff of the AMS, while for the SEARCH sites the sum of the Aitken and accumulation modes was applied. All model data are produced from the EUS simulation, which uses SAPRC07tic and occurs during June 2013.

[Figure omitted. See PDF]

The BASE case, run for the EUS domain (Fig. ), significantly underpredicted observations at both Centreville, Look Rock, and Yorkville (daytime MB $=$ $-$2.42, $-$3.39, and $-$0.66 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$, respectively) although the OA update does improve predictions slightly compared to the nvPOA case (daytime MB $=$ $-$2.72, $-$3.78, and $-$1.16 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$, respectively). Correlation coefficients improved at all locations during both day and night, with one exception at Centreville during night. Biogenic SOA was the main OA component during SOAS , and CMAQv5.2 potentially under-represents these compounds due to uncertainties in NO${}_x$-dependence, oxidant loadings, and missing organic-water interactions , as they are not yet included in this version of the model. Although previous studies have demonstrated vapor wall losses to have a relatively small effect on SOA formation from $\mathit{\alpha }$-pinene ozonolysis and $\mathit{\beta }$-pinene oxidation by nitrate radicals , it is unclear to what extent vapor losses may influence SOA yields from other monoterpene and sesquiterpene precursors or from reaction with OH. Further studies on vapor wall losses are necessary. The BASE case predicts that the composition of total OA at Centreville was 65 % biogenic SOA, 30 % combustion SOA, and 5 % POA. A prior simulation using GEOS-Chem and a similar implementation of the SIMPLE parameterization for urban and biomass burning SOA identified a contribution of 28 % of the OA from urban sources, which was consistent with the fossil fraction of the carbon at the Centreville site . Predictions from the BASE model are similar, especially since combustion SOA and POA are from both anthropogenic and biomass burning combustion sources. The model predictions for the EUS domain include semi-explicit treatments of isoprene SOA and terpene nitrate SOA formation , following the approach of .

The inclusion of semivolatile POA and pcSOA did not lead to much change in southeastern US urban centers, Birmingham and Atlanta (see Fig. S4), but there was improvement at night in both locations (MB reductions of $\sim$ 50 %). estimated that primary organic carbon was overpredicted by a factor of 1.8 in the southeast US, thereby compensating for the underestimated SOA in that study.

Figure  compares CMAQ OA species to AMS factors derived from PMF analysis at Pasadena, Bakersfield, and Cool. PMF analysis was not available at the Sacramento site, and the Centreville and Look Rock sites showed negligible presence of HOA throughout the SOAS campaign. For the nvPOA case, POA was calculated as the sum of CMAQ species POC and PNCOM and evaluated against observed HOA, even though some fraction of the material is aged. showed that this aging is predicted to add a modest amount of non-carbon mass to POA in California ($\sim$ 20 % or less). For the BASE case, POA was summed from the fresh primary organic emissions (LVPO1, SVPO1, SVPO2, SVPO3, and IVPO1) and evaluated against HOA. This model–observation comparison clearly indicates, as many past studies have, that most OA, even in urban areas is partially or heavily oxygenated and that models using an approach like the nvPOA case replace OA formed by secondary processing with primary emissions. Thus, those models get an approximately correct estimate of OA concentrations for the wrong reasons and often with the wrong diurnal profile see Sect. ;. At Pasadena, the nvPOA case overestimated HOA concentrations by a factor of 2.0. Treating the POA as semivolatile in the model resolved this discrepancy by reducing the predicted HOA factor. The combination of primary vapor aging and introduction of pcSOA resolved the OOA underprediction (MB improved from $-$4.45 to $-$0.81 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$). Figure  shows similar results at Bakersfield and Cool: average HOA predictions were reduced dramatically from the nvPOA to BASE case, while OOA predictions increased but still underpredicted observations by $-$1.3 and $-$0.28 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$. The underprediction disappears if the HOA- and OOA-observed and predicted concentrations are normalized by the CO enhancement ($\mathrm{\Delta }$CO; see Fig. S2). Because normalizing by $\mathrm{\Delta }$CO reduces the influence of errors in transport and dilution, this metric better isolates the performance of the specific chemical and microphysical processes under investigation here. However, confidence in CO model performance from CMAQ is uncertain and should be evaluated regularly and in greater detail in the future. and noted that the CO concentrations predicted by CMAQ during CalNex were systematically 10–30 % lower, and potentially emblematic of a general limitation of CMAQ capabilities for CO prediction.

Comparison of OA factors derived from positive matrix factorization (PMF) of HR-ToF-AMS observations (gray) at sites in California (Pasadena, Bakersfield, and Cool). Also shown are the model-predicted concentrations of each factor at each site using the nonvolatile (pink) and base-case semivolatile (green) configurations. (a) compares estimations for hydrocarbon-like OA (HOA) and (b) for oxygenated OA (OOA). The boxes denote the 25th and 75th percentiles of each dataset, while the whiskers extend to the most extreme points.

[Figure omitted. See PDF]

Both and used the slope of OOA as a function of odd oxygen (O${}_x$ $=$ ${\mathrm{O}}_{\mathrm{3}}$ $+$ ${\mathrm{NO}}_{\mathrm{2}}$) as an indicator of how well the model captured the general magnitude for urban SOA formation and its dependence on available oxidants. found this slope at Pasadena to be 146 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$, which is more consistent with the BASE model prediction (103 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$) than with the nvPOA model (7 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$). While the BASE model does better in capturing urban OOA, it still underpredicts its magnitude of OOA formation by 35 % according to this metric. The underprediction appears more pronounced at the highest O${}_x$ levels, according to Fig. . found similar improvement in this metric when applying a version of CMAQ that explicitly included aging of first-generation anthropogenic OOA (slope equal to 72 $\mathrm{\mu }$ g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$). They noted that CMAQ predictions could match the observed slope well if reactive anthropogenic SOA precursor vapors were added to account for missing SVOCs and IVOCs. It is not surprising that the BASE model performance falls between these two cases since the pcSOA model species was designed to account for these and other uncertain SOA formation processes. Model performance between the nvPOA and BASE cases improves similarly at Bakersfield ($m_{\mathrm{obs}}=\mathrm{80}$ $\mathrm{\mu }$ g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$, $m_{\mathrm{svPOA}}=\mathrm{54}$ $\mathrm{\mu }$ g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$) and Cool ($m_{\mathrm{obs}}=\mathrm{75}$ $\mathrm{\mu }$ g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$, $m_{\mathrm{svPOA}}=\mathrm{64}$ $\mathrm{\mu }$ g m${}^{-\mathrm{3}}$ ppbv${}^{-\mathrm{1}}$). As with the HOA–OOA comparison, normalizing by $\mathrm{\Delta }$CO improves the model performance even further, possibly canceling errors in meteorological phenomena and dispersion processes (see Fig. S3). Because CO performance is uncertain in CMAQ as previously stated, a detailed examination of the sensitivity to the pcVOC emission scale factor and OH reaction rate constant is useful.

OOA concentrations ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$) at California sites as a function of the measured O${}_x$ (O${}_{\mathrm{3}}$ $+$ NO${}_{\mathrm{2}}$) concentration. A background concentration, O${}_x$ (BG), of 13.5 ppbv is assumed, consistent with . The boxes behind each trend indicate the 25th and 75th percentiles of the data. The whiskers extend to the 10th and 90th percentiles. The solid horizontal lines in each box identify the median, and the solid curves indicate the means of each model run. Model values are projected to PM${}_{\mathrm{1}}$ to correspond roughly to the size cutoff of the AMS.

[Figure omitted. See PDF]

Diurnal OA concentration profile ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$) observed (black line, gray bars) with the HR-ToF-AMS at sites investigated during the CalNex, CARES, and SOAS campaigns. Also shown are the model-predicted distributions for the nvPOA case, the BASE case, and two sensitivity cases (LEBR and HEBR). The boxes behind each trend indicate the 25th and 75th percentiles of the data. The whiskers extend to the 10th and 90th percentiles. The solid horizontal lines in each box identify the median, and the solid curves indicate the means of each model run. Model values are projected to PM${}_{\mathrm{1}}$ to correspond roughly to the size cutoff of the AMS.

[Figure omitted. See PDF]

Model sensitivity to uncertain pcSOA parameters

Figure  shows diurnal-averaged profiles for the observations and several model cases including nvPOA, BASE, HEBR (high emission, base reaction), and LEBR. The cases with pcSOA captured the variability throughout the day, including the timing of the afternoon OA peak, and the day-to-day variability, at all California sites better than the nvPOA case. At Pasadena, the mean bias improved from $-$4.01 to $-$2.48 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ and the correlation coefficient increased from 0.15 to 0.81. Performance at Sacramento was particularly improved in the BASE case (MNE $=$ 0.24, $r=\mathrm{0.91}$). The pcVOC emission scale factors for HEBR and LEBR were set to the upper and lower bounds of the optimal parameter set reported by . The HEBR case predicted the magnitude of the afternoon OA peak at Bakersfield, Sacramento, and Cool but still underpredicted the peak at Pasadena. However, that underprediction is consistent with the underprediction at night and so may be a result of missing nighttime sources or processes. The effect of perturbing pcVOC emissions on concentrations at Centreville and Look Rock was more subdued. At those sites, pcSOA contributed 24 and 36 % to the total OA mass in the PM${}_{\mathrm{1}}$ range, respectively, in the BASE case. Generally, pcSOA concentrations were reduced by 30 % when the LEBR parameters were employed. Accounting for the pcSOA reduction and the effect of decreased total-OA-absorbing media on other OA surrogates resulted in a 13 % decrease in OA PM${}_{\mathrm{1}}$ concentrations.

We then varied both of the pcSOA parameters in the BASE model configuration as described in Table  to quantify the sensitivity of model predictions holistically. Figure  shows the NMBF and NMEF metrics calculated for each combination of input parameters investigated as well as those of the nvPOA run. These factors were computed for model–observation pairs grouped by each station individually and then averaged together. The conclusions of this analysis did not change when these metrics were calculated for all the pairs without grouping into sites first. When comparing against continuous data collected during CalNex and CARES, the nvPOA case performed worse than every sensitivity case for both metrics. For NMBF, the optimal emission factor likely lies between the base and high values, as long as the appropriate OH reaction rate is chosen to constrain the total mass formed. It is not clear whether the parameter space explored constrained the NMEF values, although they may possibly plateau at 0.48 as one moves up and to the right on Fig. . These findings are consistent with those of , although that analysis parameterized the VOC precursor emissions to the emission of CO rather than POA. The emission factor used here avoids the potentially problematic CO model predictions by relying on the more robust POA emission rate as an indicator. An area of future work will be to develop the capability to use total VOC emissions from all combustion sources as the basis for pcVOC emissions since they should be even more correlated.

Evaluation of the effect of uncertain parameters for pcSOA formation, emission factor, and OH reaction rate constant for the CAL domain. The normalized mean bias factor (NMBF) and normalized mean error factor (NMEF) symmetrically represent both under- and overprediction.

[Figure omitted. See PDF]

Evaluations of NMBF and NMEF for the CONUS11 simulations (Fig. ) result in different optimal pcSOA parameters. Daily averaged model predictions are paired with observations at sites from both the IMPROVE and CSN networks for January and July 2011. For this application, the nvPOA case results in a lower NMEF (0.83) than half of the sensitivity cases, which are those with emission factors or reaction rates higher than those of the BASE model (NMEF $=$ 0.82). Further, the nvPOA case has a lower NMBF (0.39) than all cases but the LEBR case (0.27). The parameter space is unbounded here, and even lower emission factor–reaction rate combinations could potentially reduce the NMBF and NMEF, but those solutions are not supported by the analysis for the CAL domain. The need for higher emissions or reaction rates to support higher OA concentrations at the California sites is not corroborated by the CONUS11 data, indicating that differences in the sources or production pathways that drive OA concentrations at these two scales are not completely captured with the uniform application of pcSOA. In the eastern US, for example, relative humidity or acidity may play a larger role for OA formation and deposition. Meanwhile in the wintertime cases, large wood burning area emissions may not result in substantial net OA formation downwind. Although the simulation results shown here take this possible feature into account for wildfire sources, residential and other smaller-scale wood combustion sources are assumed to produce pcSOA consistent with fossil fuel sources. If these sources are significantly overpredicted, then lower pcSOA production rates would yield better agreement at the CONUS scale for the wrong reasons. This possibility will be explored by applying a source-resolved emission inventory to determine how urban SOA formation from vehicles and cooking sources might be better simulated in the future.

Evaluation of the effect of uncertain parameters for pcSOA formation, emission factor, and OH reaction rate constant against OC measurements from IMPROVE and CSN networks. Data from January and July only of the CONUS11 simulations were used for this analysis.

[Figure omitted. See PDF]

Continent-wide organic carbon (OC) evaluation of the nvPOA case and the best-performing sensitivity case (LEBR) against routine-monitoring data from IMPROVE and CSN sites for an annual simulation during 2011. (a) Distribution of OC observations from both networks individually and combined. (b) Distributions of OC bias for each network and combined. (c) Histogram of error changes from the nvPOA to the svPOA (LEBR) case at specific stations aggregated throughout the year. (d) Regional distribution of OC observations and predictions throughout the annual simulation.

[Figure omitted. See PDF]

Evaluation against routine-monitoring data

We extended the analysis from specific measurement campaigns to a thorough evaluation of a year-long dataset at the continental US scale, CONUS11. For this application, we used daily averaged results from the LEBR model case, which performed best in the January and July CONUS11 scenarios (discussed in Sect. ). Figure  summarizes the model performance across both the IMPROVE and CSN networks for the entire year. The MB across all CSN model–observation pairs decreased substantially from the nvPOA (0.88) to the LEBR (0.39) case, while the MB at IMPROVE sites stayed roughly the same. The aggregate improvement across all sites resulted in about a 50 % reduction in MB, a 12 % reduction in mean error, and very little change to the correlation coefficient. As shown in Fig. c, error decreases occurred at a majority of sites, with large error decreases at some sites (up to and exceeding 3 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$). A minority of sites showed small increases in error, but none showed increases above 0.5 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$. Segregating the data regionally (Fig. d) indicated similar changes to the model OC distribution across the US. The new model yielded slightly less variable OC values (i.e., shorter whiskers) and only small changes to the median OC predictions. Changes to the explanatory power of the model data were mixed, with some regions indicating better correlations (northeast, southeast, Plains), some worse correlations (Midwest, southwest) and some staying the same (northwest). The consistency between the performance of the nvPOA and LEBR cases for the CONUS11 domain was likely due to the coordination of several confounding factors including spatiotemporal variability, uncertainty in emissions and reactivity, and meteorological errors (e.g., wet deposition, boundary layer dynamics). The evaluations presented in Sect.  demonstrated the extensive influence of temporal averaging: at many sites throughout the CONUS domain, the nvPOA case likely overestimated OC concentrations at night and underestimated them during the day, and these compensating errors averaged out. This stresses the critical need to evaluate models like CMAQ at multiple spatial and temporal scales where data are available.

Figure  illustrates the role of seasonal variability. The nvPOA case overpredicted daily averaged OC concentrations in the winter months and underpredicted them in the summer months with correlation coefficients in the range of 0.3–0.64. This bias trend is a well-documented feature of models that employ the nvPOA in general . The LEBR case reduced almost all of the bias in October–December 2011 and some of the bias in January and February, possibly a result of the more accurate temperature dependence in the semivolatile POA model. POA plays a larger role in the winter months than in the summer because of favorable temperature-dependent partitioning. However, without the pcSOA species, the LEBR case would likely have slightly underpredicted the observations in the winter time (see Figs. S5–S7). Meanwhile, performance during the summer months improved slightly for the LEBR case, with significant reduction in bias for July. The transition season months showed smaller bias reductions, with 0.22, 0.14, 0.12, and 0.14 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ bias reductions for March, April, September, and October, respectively. During these months, the opposing changes induced by POA volatilization and pcSOA formation had highly variable impacts across the US, with dependence on temperature, boundary layer height, and oxidant loadings. The correlation coefficient did not change appreciably from the nvPOA case to the LEBR case, in contrast to the dramatic changes in correlation seen when comparing against California hourly data. A regional evaluation of the NMBF and NMEF metrics broken down by season illustrates that performance in almost every region of the US improved throughout the year from the nvPOA to the LEBR case (Fig. ). The nvPOA case had the worst-performing predictions during summer, and the LEBR updates moved all but the northeast region into the statistical area defined as acceptable by . Performance during the summer improved for all regions and improved substantially for the northwest. The NMBF was reduced for all regions in the fall scenarios although the NMEF changed only slightly. The only degraded performance occurred in the northwest in spring, where the NMBF increased by about a factor of 2 and the NMEF by a factor of 3. This is likely driven by errors in the emissions and aging of POA from fire events and should be investigated in the future as more is learned about the identity and fate of the compounds emitted by this group of sources.

Monthly model performance for the nvPOA (blue) and LEBR (orange) cases throughout the 2011 simulation. The top panel indicates the mean OC bias, while the bottom panel indicates the correlation coefficient of each model run with the CSN and IMPROVE observations.

[Figure omitted. See PDF]

OA spatiotemporal composition in CMAQ

The total OA surface concentrations predicted by the model vary considerably from winter to summer for many regions of the nation (top row of Fig. ). In the winter, urban locations emerge as the sources mainly driving the regional OA distribution, whereas, in the summer, these urban centers compete with heavily forested regions of the country like the southeast and the Sierra Nevadas. Although the abundance of biogenic OA observed for 6 weeks at the Centreville site is qualitatively reproduced by the model, significant uncertainties exist in our knowledge of the air–surface exchange of important biogenic precursors, their SOA yields, and the partitioning properties of their oxidation products. The annually averaged contribution of POA (unreacted primary compounds in the particulate phase) predicted by the model is in the range of 5–10 %, increasing to 30 % in some cities, especially at high elevations with colder temperatures, which is in the same range as the observations . These maxima are driven by wintertime episodes when the atmospheric lifetime of POA against evaporation and gas-phase oxidation increases substantially. The LEBR model predicted that in the summer there is little contribution of these fresh emissions to particulate loadings as most of the mass was predicted to be oxidized quickly in the gas phase. CMAQv5.2 with the LEBR configuration should be evaluated against AMS-inferred HOA concentrations at cities characterized by higher wintertime POA fractions (e.g., ).

Regional and seasonal model performance at routine-monitoring stations for the CONUS11 simulation as a function of the normalized mean bias factor (NMBF) and normalized mean error factor (NMEF) for OC predictions. Plotted regions include the northeast (NE), Midwest (MW), southeast (SE), Plains (PL), northwest (NW), southwest (SW), and total (T) and are defined as visualized in Fig. .

[Figure omitted. See PDF]

Spatial distribution of products form the LEBR CONUS11 simulation. Products include the OA concentration (row 1), fraction of POA (row 2), pcSOA concentration (row 3), and change in model-predicted POA from the nonvolatile POA model to the semivolatile POA model (row 4). Maps illustrate the median of all annual surface data (a, d, g, j), winter months (December, January, February; b, e, h, k), and summer months (June, July, August; c, f, i, l).

[Figure omitted. See PDF]

Organic carbon (OC) distributions for routine-monitoring stations and model predictions at the same locations and times. The meaning of the boxes and whiskers is explained in Fig. .

[Figure omitted. See PDF]

The pcSOA concentrations are driven both by variability in primary pcVOC emissions and by efficient photochemical oxidation, as illustrated by Fig.  (third row). The summertime concentration field was similar to that of the total OA concentrations and was relatively dispersed in the eastern US like other regional secondary pollutants (e.g., ozone, particulate nitrate). Interestingly, the wintertime concentration field in the east was more dispersed than in the summer because of the longer lifetime of pcVOC against oxidation to pcSOG. The increased role of residential wood combustion emissions led to elevated pcSOA concentrations in the northeast US compared to those in the summer. The pcVOC emission sources in the western US were visually evident from the pcSOA concentration field in both the summer and winter simulations, mostly due to relatively less population density in western states compared to the east. An interesting dynamic associated with pcSOA appeared in California. In the summer, enhanced photochemistry caused high pcSOA concentrations in southern California (4 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$). Meanwhile, concentrations were elevated but not as much in the San Joaquin Valley. The relationship was flipped in the wintertime where concentrations were low in southern California and much higher in the valley. The wintertime San Joaquin Valley enhancement was likely due to lower boundary layer heights enhancing the oxidative capacity of the airshed and increasing the formation rate of pcSOA. Future work will investigate episodes like this to determine if more can be learned about how the contributions of various SOA pathways vary seasonally. The percent change in POA concentrations from the nvPOA to the LEBR cases was highest during the summer and far from urban sources, as expected. For this study, the POA boundary conditions were renamed SOA and highly affect the percent change in POA close to the model boundaries. In the summer time, the change in POA concentrations in the domain interior varied from $-$100 to $-$70 % change in California. The winter simulations showed more variability ($-$100 to $-$55 % change), indicating that even when treating POA as semivolatile, it is still an important component of urban and suburban particulate loadings and important to consider in exposure assessments.

Multiyear organic carbon trends

For many policy and environmental applications, it is important for the model to capture trends in ambient concentrations in response to emissions changes. In order to probe CMAQv5.2 sensitivity to changes in real-world emissions from 2002 to 2011, we applied the model (both the nvPOA and LEBR cases) to the continental US domain during January and July 2002 (CONUS02) and compared the modeled OC predictions to available CSN and IMPROVE observations (Fig. ). The OC trend for IMPROVE sites between January 2002 and January 2011 were relatively flat and the model reproduced the behavior for both cases. The January CSN site trends were more complex. While the observations indicated a decrease of 0.92 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$ yr${}^{-\mathrm{1}}$, both the nvPOA and LEBR model cases predicted increases of 0.87 and 0.65 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$ yr${}^{-\mathrm{1}}$, respectively. Because this deviation occurs in the context of wintertime scenarios, uncertainties from emissions or misrepresented boundary layer dynamics likely play a more important role. With more realistic emissions from this source, both CMAQv5.2 model configurations would predict greater OC, potentially overpredicting the observations by a similar magnitude as the 2011 datasets.

Both CSN and IMPROVE sites during July show a similar trend, although it is exacerbated for CSN sites. Here, the observations indicate dramatic OC decreases ($-$5.72 % ${\mathrm{yr}}^{-\mathrm{1}}$), while the nvPOA and LEBR models predict increases by 4.28 % and 3.65 % ${\mathrm{yr}}^{-\mathrm{1}}$, respectively. Clearly the model is missing an important source or pathway for OC formation in the 2002 simulation. The 2011 predictions overlap reasonably with the observations, so this pathway either diminishes substantially when moving to the 2011 scenario or there is a coincident offset by other pathways that may be overestimated by the model in 2011. The 2002 CSN observations are potentially biased by the application of the National Institute for Occupational Safety and Health (NIOSH) method; the IMPROVE (2002 and 2011) and the CSN 2011 observations were gathered using the IMPROVE method, which is consistent with the characterization of the POA emission factors used to build the emission inventory. Because the differences between the NIOSH and IMPROVE methods affects the observed ratio of OC to elemental carbon (EC), the potential magnitude of OC bias can be estimated by inspecting the bias in EC. For the CSN 2002 dataset, the MB for EC is 0.2 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$, indicating that most of the underprediction at CSN sites in summer 2002 is not related to discrepancies between measurement methods.

One possible formation pathway of particle-phase OC mass that could have been higher in 2002 involves perturbations of aerosol liquid water concentrations. , , and argued that higher aerosol water concentrations, possibly resulting from increased sulfate concentrations in the past, could have enhanced the partitioning of semivolatile organic compounds, depending on their solubility, or they could have enhanced reactive uptake rates of VOC oxidation products. This is consistent with , who showed that 88 % of OA were water soluble at the Centreville site during SOAS and that only 25–40 % of the total water-soluble organic compounds were in the particle phase throughout the day. Figure  indicates the trend toward higher OC biases as sulfate increases. This trend occurred despite the fact that sulfate and elemental carbon aerosol predictions did not underestimate observations as did OC predictions. CMAQv5.2 includes the aerosol–water dependent pathway of isoprene OA formation via IEPOX (isoprene-epoxydiol) processing. A similar interaction among products of anthropogenic VOC oxidation may be responsible for the discrepancy in the urban-influenced CSN sites. Another potential OC pathway, organic nitrate formation, could have played a larger role in 2002 and diminished over the last decade due to decreased NO${}_x$ emissions (e.g., ). Those authors showed that a 25 % reduction in NO${}_x$ led to a 25 % reduction in organic nitrate SOA. Since 2002, NO${}_x$ concentrations in the US have decreased due to the implementation of the NO${}_x$ Budget Trading Program (or NO${}_x$ SIP Call and the Clean Air Interstate Rule). For the combined 2002 and 2011 CSN dataset, the OC underprediction improved from almost $-$5 $\mathrm{\mu }$gC m${}^{-\mathrm{3}}$ at 50–60 pbbV NO${}_x$ to unbiased at 10 ppbv NO${}_x$ and less. Although these relationships are not mechanistic proof of a relationship between OC aerosol formation and ${\mathrm{SO}}_{\mathrm{2}}$ or NO${}_x$ emissions, they do encourage further scrutiny of potential co-benefits that reducing the latter two pollutants has on OA burden and thus public health. Finally, it is possible that the ratio of POA to intermediate and high-volatility unspeciated organic combustion emissions changed from 2002 to 2011 as new emissions controls and fuel formulations were adopted . Thus, the analysis presented here may speak to the robustness of using the pcSOA approach for regulatory applications. Fully quantifying the impacts of the pcSOA assumption would require comparing to simulations that include both speciated, source-specific, and year-specific NMOG estimates in the emissions input generation process and accurate speciated SOA yields in the CTM.

Organic carbon bias for CSN summer time data during 2002 and 2011 (combined) as a function of observed sulfate concentrations (a) and observed NO${}_x$ concentrations (b). The meaning of boxes and whiskers is explained in Fig. . Blue (2011) and black (2002) solid lines quantify the number of data points for each year as functions of the observed pollutants.

[Figure omitted. See PDF]

Conclusions

We have shown that CMAQv5.2 with a semivolatile representation of POA compounds and an introduction of the lumped pcSOA model species predicts organic aerosol concentrations acceptably against high-time-resolution field campaign data and daily averaged routine-monitoring network products. In almost all cases, the model improves diurnal and seasonal trends compared to the model with nonvolatile POA – a consequence of better representation of both temperature sensitivity and photochemical oxidation dependence. Because the BASE model is able to achieve similar performance to other more-complicated approaches to OA formation, we advise the OA fate and transport modeling community to consider adding complexity only when there is enough experimental data to provide an associated reduction in uncertainty.

The discrepancy between the optimal parameters inferred from the CAL and CONUS11 simulations highlights the importance of not overstating the meaningfulness of pcSOA or the value of the optimal parameters presented here. Because these parameters are fit to available observations, they lose some explanatory power as one extrapolates to other domains or time periods. However, the essential features of modeled pcSOA formation, the oxidation of a VOC precursor followed by condensation to the particle phase, broadly represent the transformations of most of the likely pathways to OA formation needed by the model. This approach parallels the more-detailed, but still very uncertain, approaches used in similar CTMs simulating OA from primary sources while maintaining consistency with existing regulatory inventories developed over decades. For these reasons, we feel confident that CMAQv5.2 with semivolatile POA and pcSOA is a substantial step toward a more accurate representation of OA formation, even though future work will focus on replacing the pcSOA species with specific mechanisms constrained by direct measurements.

Based on the limitations discussed in this work, we recommend the chemical transport modeling community address the following concerns:

• The entire volatility spectrum of primary organic emissions, including IVOCs, should be incorporated directly into emission inventories and applied to model input as a function of relevant parameters like combustion technology and fuel type.

• Because of their unique properties and variable OA formation potential, organic compound emissions from biomass burning sources will likely need to be treated distinctly from fossil fuel combustion sources in CTMs. For wildfires the net added mass due to SOA formation, relative to the initial POA mass, is 1–2 orders of magnitude lower than for urban sources .

• Aerosol–water interactions are an important complexity to incorporate consistently in CTMs. These interactions not only impact OA loading but have important implications for deposition rates, diurnal profiles, multiyear trends, and model response to the simulated control strategy.

• Organic nitrates were detected at multiple sites during the CalNex, CARES, and SOAS campaigns. A consistent treatment of these compounds and the effect nitrate group addition has on volatility and solubility is needed.

• In general, greater efforts are needed to improve the conceptual link between the gas- and aerosol-phase chemical components in CTMs. An example of this approach is described for isoprene epoxide degradation by and monoterpene nitrate formation by .

We further identify the following potential measurement efforts as useful for addressing critical gaps in our current understanding of OA formation:

• Winter-time ambient observations of urban OA in multiple cities with varying meteorology. Observations are needed at high time resolution (at or greater than 1 h) in order to discern important sensitivities to photooxidation, temperature, and atmospheric water content. Important properties to constrain are the POA fraction, volatility/solubility, and aging timescales.

• Detailed transport model comparisons to existing and upcoming measurements of ambient biomass burning plumes including characterization of potential downwind SOA enhancement as well as aging effects on volatility/solubility (e.g., ).

• Updated quantitative estimates of SOA yields as a function of volatility/solubility that take into account losses of low-volatility vapors to chamber walls.

We conclude that the important OA mass formation processes are most likely regional in scale and that OA concentrations tend to be regionally distributed in the US, especially in summer. However, we also emphasize that, despite the diminished contribution direct POA emissions have in the model predictions as a result of evaporation, they are still a critical component to consider, especially for exposure assessments near primary sources, urban centers, and in the winter. Moreover, the evaporated POA compounds likely react in the vapor phase and may contribute to particle concentrations again quickly or they may contribute to the formation of oligomer species, about which little is known with certainty.

The multiyear trends in Figs.  and  highlight the difficulty in integrating OA formation and properties into models to be used for regulatory applications. If interactions with regulated pollutants (e.g., ${\mathrm{SO}}_{\mathrm{2}}$, NO${}_x$) are important contributors to OA formation, then the applicability of a model configuration to a specific scenario may change over time as the concentration of those pollutants decreases. Pathways that used to be important may become less so and vice versa, thus demonstrating the benefit of evaluating the model for a wide array of atmospheric conditions. Characterizing model results with this dynamic view helps improve the robustness of the model over time and guide future policy decisions and research.

The CMAQv5.2 code for both the semivolatile and nonvolatile POA options is available via the GitHub repository (https://github.com/usepa/cmaq). Data included in figures are available upon request to the contact author. CalNex measurement data are available through the Earth System Research Laboratory data portal (https://www.esrl.noaa.gov/csd/groups/csd7/measurements/2010calnex/, CalNex measurement data, 2012), CARES data are available through the DOE Atmospheric Radiation Measurement program portal (http://www.archive.arm.gov/discovery/#v/results/s/fiop::aaf2009carbonaerosol, CARES measurement data, 2010), and SOAS data are available through the Southeast Atmosphere Study website (https://esrl.noaa.gov/csd/groups/csd7/measurements/2013senex/, SOAS measurement data, 2014).

The Supplement related to this article is available online at https://doi.org/10.5194/acp-17-11107-2017-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

The US EPA, through its Office of Research and Development, supported the research described here. The research has been subjected to agency administrative review and approved for publication but may not necessarily reflect official agency policy. The authors gratefully acknowledge Heather Simon (Office of Air Quality Planning and Standards) and Brian Eder (National Exposure Research Laboratory) from the US EPA for their helpful comments during the agency review process. Jose L. Jimenez was supported by EPA STAR 83587701-0 and National Science Foundation (NSF) AGS-1360834. Shang Liu and Lynn M. Russell were supported by the California Air Resources Board (contract 09–328). Funding for data collection during the CARES campaign was provided by the Atmospheric Radiation Measurement (ARM) Program sponsored by the US Department of Energy (DOE), Office of Biological and Environmental Research (OBER). Rahul A. Zaveri was supported by the US DOE's Atmospheric Systems Research (ASR) program under contract DE-AC06-76RLO 1830 at the Pacific Northwest National Laboratory. Qi Zhang and Ari Setyan were supported by DOE ASR DE-FG02-11ER65293. Lu Xu and Nga L. Ng acknowledge NSF grant 1242258 and US EPA STAR grant RD-83540301. Edited by: Manabu Shiraiwa Reviewed by: two anonymous referees