Regional Sources of Airborne Ultrafine Particle Number and Mass Concentrations in California

Regional concentrations and source contributions are calculated for airborne particle number concentration (Nx) and ultrafine particle mass concentration (PM0.1) in the San 15 Francisco Bay Area (SFBA) and the South Coast Air Basin (SoCAB) surrounding Los Angeles with 4 km spatial resolution and daily time resolution for selected months in the years 2012, 2015, and 2016. Performance statistics for daily predictions of N10 concentrations meet the goals typically used for modeling of PM2.5 (MFB< ± 0.5 and MFE < 0.75). The relative ranking and concentration range of source contributions to 20 PM0.1 predicted by regional calculations agree with results from receptor-based studies that use molecular markers for source apportionment at four locations in California. Different sources dominated regional concentrations of N10 and PM0.1 because of the different emitted particle size distributions and different choices for heating fuels. Nucleation (24-57%) made the largest single contribution to N10 concentrations at the ten 25 regional monitoring locations, followed by natural gas combustion (28-45%), aircraft (210%), mobile sources (1-5%), food cooking (1-2%), and wood smoke (0-1%). In contrast, natural gas combustion (22-52%) was the largest source of PM0.1 followed by mobile sources (15-42%), food cooking (4%-14%), wood combustion (1-12%), and aircraft (2-6%). The study region encompassed in this project is home to more than 25M 30 residents, which should provide sufficient power for future epidemiological studies on the


Introduction
Numerous epidemiological studies have identified positive correlations between exposure to ambient particulate matter (PM) and increased risk of respiratory and cardiovascular diseases, premature mortality and hospitalization (Pope et al., 2002;Pope et al., 2004;Pope et al., 2009;Dockery and Stone, 2007;Ostro et al., 2015;Ostro et al., 40 2006; Ostro et al., 2010;Brunekreef and Forsberg, 2005;Fann et al., 2012;Gauderman et al., 2015;Miller et al., 2007). Most of these studies have not fully addressed ultrafine particles (UFPs; Dp<0.1µm) because these particles make a very small contribution to total ambient PM mass (Ogulei et al., 2007). Toxicity studies suggest that UFPs may be especially dangerous to human health since they have higher toxicity per unit mass (Li et 45 al., 2003;Nel et al., 2006;Oberdorster et al., 2002) and can penatrate the lungs and enter the bloodstream and secondary organs . These toxicology results are suggestive but more epidemiological evidence is required before the threat to public health from UFPs can be fully assessed.
Most previous UFP epidemiology studies are based on particle number concentration (Nx 50 -the number of particles with diameter less than X nm) measured at fixed sites using commercially-available instruments. These devices are expensive and they require regular maintence which limits the number of measurement sites that can be deployed.
Translating measured Nx into population exposure estimates is also difficult because UFP concentrations change more rapidly over shorter distances than PM2.5 (Hu et al.,55 2014b; Hu et al., 2015;Hu et al., 2014a). Land use regression (LUR) models could potentially be used to interpolate UFP concentrations between sparse measurement locations, but the atmospheric processes governing Nx concentrations are highly nonlinear and (so far) sufficient training data is not generally available for LUR models to estimate Nx exposure over a large enough population to support a definitive epidemiology 60 study (Montagne et al., 2015). Previous attempts to use regional reactive chemical transport models to predict Nx in highly populated regions have focused on nucleation, yieldeding a wide range of predicted concentrations and only modest agreement with measurements when different nucleation algorithms were used (Elleman and Covert, 2009b;Zhang et al., 2010;Elleman and Covert, 2009a). Obtaining accurate exposure 65 estimates to Nx in highly populated regions therefore remains a major challenge in UFP epidemiological studies.
Recent work has examined UFP mass (PM0.1) as an alternative metric for UFP exposure, and demonstrated that PM0.1 can be predicted with reasonable accuracy over large populations using regional reactive chemical transport models (Hu et al., 2014b;Hu et al., 70 2014a). The PM0.1 exposure fields developed using this technique have been used in multiple epidemiological studies that revealed associations with mortality and pre-term birth (Ostro et al., 2015;Laurent et al., 2016). Despite the success of studies using PM0.1, techniques that estimate Nx exposure are still needed because a large number of ongoing UFP studies are based on Nx and it is possible that PM0.1 and Nx are associated with 75 different types of health effects.
Here we extend the previous work using regional reactive chemical transport models for UFPs to include Nx in the San Francisco Bay Area (SFBA) and the South Coast Air Basin (SoCAB) region around Los Angels which are the two most densely populated major metropolitan location in California. Source contributions to PM0.1 and Nx are tracked 80 using the University of California, Davis / California Institute of Technology (UCD/CIT) regional reactive chemical transport model with 4 km spatial resolution. Predicted concentrations during the year 2012 are compared to measurements available at ten regional monitoring sites. The spatial distribution fields of different particle metrics (Nx, PM0.1, PM2.5) are combined with population distributions to estimate exposure. To the 85 best of our knowledge, this is the first integrated study of both UFP number and mass using a regional reactive chemical transport model in California.

Model Description
The UCD/CIT chemical transport model used in the current study has been succesfully applied in sevaral previous studies in the San Joaquin Valley (SJV) and the SoCAB (Ying 90 et al., 2008b;Ying et al., 2008a;Hu et al., 2015;Hu et al., 2017;Chen et al., 2010;Held et al., 2004;Held et al., 2005;Hixson et al., 2010;Hixson et al., 2012;Hu et al., 2012;Kleeman and Cass, 2001;Kleeman et al., 2007;Kleeman et al., 1997;Mysliwiec and Kleeman, 2002;Rasmussen et al., 2013;Ying and Kleeman, 2006;Zhang and Ying, 2010). The model includes algorithms for emissions, transport, 95 dry deposition, wet deposition, gas phase chemistry, gas-to-particle conversion, coagulation, and some condensed phase chemical reactions. Nucleation was added to the model for the first time in the current study using the ternary nucleation (TN) mechanism involving H2SO4-H2O-ammonia (NH3) (Napari et al., 2002). As was the case in previous studies using this algorithm, the resulting nucleation rate was adjusted using a tunable 100 nucleation parameter set to 10 -5 for new particle nucleation (Jung et al., 2010). The Kerminen and Kulmala (2002) parameterization was added in order to bridge the gap between the 1 nm particle nuclei and their appearance into the smallest size bin of the UCD/CIT model (~10 nm). The nuclei growth rate (GR) in the Kerminen and Kulmala (2002) parameterization is one of the factors that accounts for the competition between 105 the condensation and nucleation of over-saturated compounds until the nucleated particles grow to the size of the smallest bin in the regional model, at which point this competition is represented explicitly by the model operators. In the current study, the GR for nucleated sulfate particles was calculated using the diffusion-limited condensation rate of sulfuric acid based on the recommendation of Kerminen and Kulmala. Once 110 particles reach ~10nm, the full operators in the model calculations predict growth by condensation of sulfuric acid, nitric acid, ammonia, and secondary organic aerosol (SOA). Perturbation studies were conducted in the current analysis to test the effect of GR with a box model configured to represent a single grid cell using the full set of model operators. Initial conditions in the SAPRC11 gas-phase mechanism were 0.04 ppm O3, 115 0.05 ppm NO, 0.0 ppm NO2, 0.05 ppm HCHO, 0.1 ppm ISOPRENE, 0.1 ppm BENZENE, and 0.01 ppm ALK5. A nucleation event was initiated at 8am by setting H2SO4 concentrations to 10 7 molecules cm -3 and NH3 concentrations to 100 ppt. The nominal GR was multiplied by a factor ranging from 0.5 to 2.0 to test the sensitivity of the model results. Figure 1 illustrates the growth of nucleated particles between 5am and 120 12 noon for conditions representing July in California. The number concentration of nucleated particles increases from zero to values between 2500 -3000 # cm -3 . SOA condenses on the particles causing their size to increase above 100nm. Coagulation and deposition processes remove particles over time. Three separate simulations are illustrated in Figure 1 using the nominal GR along with perturbations of 0.5*GR and 125 2.0*GR. These model perturbations fall almost exactly on top of the basecase simulations, suggesting that results are not overly sensitive to GR during the first few seconds of nuclei growth before calculations are handed off to the regional model algorithms.
130 Figure 1: Simulated particle nucleation event followed by growth due to SOA condensation under conditions representing July in California. Vertical axis displays the mean diameter of the nuclei mode while color represents the particle number concentration.
Several previous modeling studies have been conducted to evaluate the performance of 135 the ternary nucleation mechanism on predicted Nx using global and regional models. Jung et al. (2010) found that a scaled version of the ternary H2SO4-NH3-H2O nucleation theory (Napari et al.,2002 with a supplemental 10 -5 nucleation tuning factor) added to the PMCAMx-UF model produced Nx predictions in reasonable agreement with observations. The study of Westervelt et al. (2013) also showed that the ternary 140 nucleation parameterization (with a supplemental 10 -5 nucleation tuning factor) added to the Goddard Earth Observing System global chemical transport model (GEOS-Chem) produced reasonable Nx predictions on average when compared with measurements at five locations spanning various environments. Jung et al. (2008) ((Napari et al., 2002) with a supplemental 10 -5 nucleation tuning factor) was a suitable nucleation scheme for 3-D chemical transport models. Although there have been numerous significant efforts to 150 incorporate nucleation algorithms into three-dimensional regional and global models (Jung et al., 2008;Jung et al., 2010;Westervelt et al., 2013;Zhang et al., 2010;Yu et al., 2015;Dunne et al., 2016;Fanourgakis et al., 2019), nucleation modeling studies are still in the early stages of development and further efforts are needed to reduce the uncertainty in both the nucleation rate and growth mechanisms.

155
In the current study, emission, transport, deposition, and coagulation of UFPs were simulated using operators developed for the UCD/CIT model framework, leading to modification of the particle size distribution and the subsequent Nx concentrations.
Dynamic condensation / evaporation is considered for all particle size bins with predicted UFP growth rates of 2-3 nm hr -1 or higher under favorable conditions. The regional 160 model operators are not well suited for the most extreme changes to the particle size distribution that occur within the first few seconds or minutes after emissions to the atmosphere (such as within 300 m of roadways). Dedicated simulations can predict the dynamic condensation/evaporation of particles at distances of 10's of meters downwind of the roadway (Zhang et al., 2005;Zhang et al., 2004) mostly due to the partitioning of 165 SOA (Anttila and Kerminen, 2003;Trostl et al., 2016), but these calculations are too expensive for domains spanning thousands of km. Regional calculations such as those illustrated in the current study rely on emissions characterization measurements that include a few minutes of aging to capture the "near-field" emissions of particle size and composition that can then be used as the starting point for regional model calculations. In 170 some cases, evaporation of UFPs in the first few seconds after release to the atmosphere is therefore represented by reducing the primary emissions of nano-particles based on measurements conducted at high dilution factors (Xue et al., 2018a) or using measurements of particle volatility to estimate the evaporation at high dilution factors (May et al., 2013a;May et al., 2013b;Kuwayama et al., 2015). All of the results presented 175 in the current analysis focus on regional UFP concentrations with 4km resolution.
The model domains used in the study are shown in Figure 2. The parent domain with 24 km horizontal resolution covered the entire state of California (referred to as CA_24 km) and the two nested domains with 4 km horizontal resolution covered the SFBA + SJV + South Sacramento Valley air basins (referred as SJV_4 km) and the SoCAB surrounding 180 Los Angeles (referred as SoCAB_4 km). The UCD/CIT model was configured with 16 vertical layers up to a height of 5 km above ground level, with 10 layers in the first 1 km.
Previous studies have shown that this vertical configuration captures the air pollution system above California (Hu et al., 2014a;Hu et al., 2014b;Hu et al., 2015). Particulate number, mass, and composition are represented in 15 size bins, with particle diameters 185 being centered within equally spaced logarithmic size interval spanning the diameter range from 0.01 to 10µm. Nucleated particles were initialized in a 16 th size bin with initial diameter of 0.007 µm.

Meteorological Fields
Hourly meteorological fields during the modeling period were generated by the Weather Research and Forecasting (WRF) model version 3.4 with three nested domains that had horizontal resolutions of 36 km, 12 km and 4 km, respectively. In the present simulations, the WRF model was configured with 50 vertical layers (up to 100 hpa) and four-dimensional data assimilation (FDDA) nudging was utilized to improve the agreement between model predictions and observed meteorological patterns (Otte, 2008b, a). WRF 205 predictions for wind speed, temperature, and relative humidity were compared to measurements for seven counties in the SFBA and two counties in SoCAB (see Table   S2). Temperature has mean bias (MB) within ~0.2 • C and root-mean-square errors (RMSE) between 4-5 • C. Wind speed has mean fraction bias (MFB) within ±0.20 and RMSE generally <2.0 m/s. This level of performance is consistent with performance of 210 WRF in previous studies conducted in California (Zhao et al., 2011;Hu et al., 2015).

Emissions
The emission inventories used in the SFBA were developed by the BAAQMD for the year 2012 based on the regulatory inventory provided by the California Air Resources Board for that same year. The SFBA inventory was processed using the Sparse Matrix Measurements conducted in parallel with the current study found that particles emitted from natural gas combustion in home appliances were semi-volatile when diluted by a factor of 25 in clean air, but particles emitted from reciprocating engines did not evaporate under the same conditions (Xue et al., 2018a). Near-field emissions from all 225 natural gas sources combustion sources other than reciprocating engines were therefore set to 30% of their nominal levels. A map of the natural gas emissions distribution is shown in Supporting Information ( Figure S3). SMOKE results were transformed into size-resolved emissions of particle number, mass, 230 and composition using measured source profiles through an updated version of the emissions model described by Kleeman and Cass (1998). The PM profiles used for each source type were specified as weighted averages from each of the detailed sources within each broad category as summarized in Table S1. Detailed PM source profiles for major sources of ultrafine particulate matter are based on measurements conducted during 235 source tests (Li and Hopke, 1993;Kleeman et al., 1999Kleeman et al., , 2000Robert et al., 2007a;Robert et al., 2007b;Mazaheri et al., 2009). In most cases, these emissions size distributions strongly influence the size distributions of particles in the ambient atmosphere (see Figures S1 and S4). A more detailed discussion of the emissions processing has been presented in a previous study .

Statistical Evaluation
According to Taylor's Hypothesis (Shet et al., 2017), it is expected that the spatial distribution of model results is more important than the temporal distribution when evaluating performance. In the current study model performance evaluations are limited 245 to the locations where measurements were made. Therefore, the temporal distribution is also considered by comparing predicted vs. measured daily average Nx, PM2.5 and individual PM2.5 species mass concentrations.
The evaluation data set was compiled from several measurement networks including the sites operated by staff at the Bay Area Air Quality Management District (BAAQMD), the 250 IMPROVE sites, the MATES IV sites and FRM sites. In order to account for the uncertainty in predicted wind fields and spatial surrogates used to place emissions, "bestfit" model results were created by identifying the closest match within 3 grid cells of each measurement location. "Best-fit" model performance for PM2.5 at routine monitoring sites ( Figure 2) meets the performance criteria suggest by Boylan and Russell (Boylan 255 and Russell, 2006) (mean fractional error (MFE) ≤ +0.75 and mean fractional bias (MFB) ≤ ±0.5) (Table S4). Table S5 shows the MFB and MFE values of gaseous species of O3, NO, NO2, CO and SO2 using daily averages across all measurement sites during the entire simulated period. Gaseous species of O3, CO, NO, NO2 and SO2 have MFBs within ± 0.3 and MFE less than 0.5, indicating consistent behavior between predictions and 260 measurement for these species. The ability of UCD/CIT predictions for key gas species, mass and chemical component concentrations in the PM0.1 and PM2.5 size fractions was also evaluated in previous studies (Ying and Kleeman, 2006;Ying et al., 2008a;Ying et al., 2008b;Hu et al., 2012;Chen et al., 2010;Held et al., 2005;Hu et al., 2015;Hu et al., 2017;Venecek et al., 2018). The performance of the UCD/CIT air quality model in these 265 studies generally meets standard model performance criteria. Of greatest interest in the current study, predicted "best-fit" N10 values were compared to measured N7 values at four sites in the SFBA (Santa Rosa, San Pablo, Redwood City and Livermore) and six sites in SoCAB (Anaheim, Central Los Angeles, Compton, Huntington, Inland-valley and Rubidoux). N7 measurements in the SFBA were made using an Environmental Particle Previous studies conducted at Fresno, California, suggest that N7-10 accounts for approximately 8% of N7 (Watson et al., 2011), and so some amount of negative bias is 275 expected when comparing predicted N10 to measured N7. The evaluation results for "best-fit" N10 summarized in Table 1 follow this expected trend but mean fractional bias (MFB) and mean fractional error (MFE) at each comparison site still meet the PM2.5 performance criteria suggested by Boylan and Russell (2006). This level of performance is comparable to the results from a previous UFP number simulation conducted in 280 Northern California using a modified version of the WRF-Chem model (Lupascu et al., 2015). The level of agreement between predicted "best-fit" and measured PM2.5, individual PM2.5 species, key gas species and N10 builds confidence in the model skill for UFP predictions in the current study.  Table 2 below summarizes the predicted correlations between daily-average particle 290 number concentrations and PM2.5 along with the measured correlations for these metrics.

PM0.1 and N10 Source Apportionment in California
The UCD/CIT model uses a moving sectional approach to conserve particle number and mass while letting particle radius change due to condensation and evaporation (Kleeman 315 et al., 1997). The method to calculate source contributions to number concentration is performed for each moving section individually. Number is explicitly conserved and correctly apportioned to sources in this algorithm. Each particle source type / moving size bin includes an artificial tracer equal to 1% of the primary particle mass. The mass of this tracer is related to the number of particles by the equation 320 tracer_source_i * 100 = N_source_i * 3.14159/6 * Dp_bin * ρ_i (eq1) where ρ_i is the density of primary particles emitted from source i. This equation can be easily rearranged to solve for N_source_i as a function of tracer_source_i in each size bin. Condensation/evaporation changes the particle diameter as semi-volatile components move on and off the particle but this does not change tracer_source_i or 325 N_source_i. As a result, the moving sectional approach greatly simplifies the source apportionment of particle number compared to other models that use fixed particle size bins with condensation / evaporation transferring material between bins.
Coagulation complicates source apportionment calculations for particle number because coagulation events conserve particle mass but destroy particle number. The model 330 calculations treat the most frequently occurring coagulation events between very small particles and very large particles in a manner analogous to condensation. When two particles coagulate, the mass of the smaller particle is added to the mass of the larger particle. The number concentration of the smaller particle is discarded while the number concentration of the larger particle stays constant. This slightly reduces the accuracy of 335 source apportionment calculations for particle number in the larger size bins because the tracer_source mass in the larger size bin is no longer proportional to the number concentration from that source. This issue is relatively minor since size bins larger than 1µm that act as the dominant sink during particle coagulation events typically account for less than 5% of the total number concentration.

355
Many of the spatial patterns measured for airborne particle number concentrations in past studies have focused on the gradients around roads (see for example (Zhu et al., 2002a;Zhu et al., 2002b;Zhang et al., 2005;Zhang et al., 2004;Sowlat et al., 2016). These gradients are impossible to resolve using a regional model with 4km resolution. A limited set of additional simulations were conducted using the WRF/Chem model 360 configured with Large Eddy Simulation (LES) around Oakland California so that spatial scales down to 250m could be examined. Maps of the predicted ultrafine particle mass concentrations for gasoline, diesel, food cooking, wood combustion, and natural gas combustion particles are shown in Figure 4 below. At 250m resolution, ultrafine particles from diesel engines peak on major transportation corridors while ultrafine 365 particles from gasoline vehicles are more diffuse reflecting their increased activity on adjacent surface streets. Ultrafine particles from natural gas combustion are even more diffuse reflecting contributions from area sources across the region. As the spatial resolution decreases to 1km and then 4km, the fine details around roadways are artificially diluted in the larger grid cells. This process shifts the dominant source of 370 ultrafine particles over roadways from diesel engines at 250m resolution to natural gas combustion at 4km resolution. These simulation results are consistent with measurements of particle number in the proximity of roadways which show that the traffic contribution to particle number concentration decays to background levels within 300 m (Zhu et al., 2002a;Zhu et al., 2002b). The measurements made by Zhu et al. 375 indicate that the traffic contribution to regional number concentration cannot be distinguished from other sources on a regional scale using 4km grid cells which is the focus of this study.

UCD/CIT PM0.1 source contributions compared to Chemical Mass Balance (CMB) results
A recently completed study measured the composition of PM0.1 at four sites in California and calculated source contributions using molecular markers (Xue et al., 2018b). Figures   5 and 6 compare the source contributions to PM0.1 OC concentrations predicted by the UCD/CIT model and "measured" using the molecular marker technique at San Pablo,

390
East Oakland, downtown Los Angeles and Fresno during a summer month (August 2015) and a winter month (February 2016). The "others" category in the molecular marker calculation represents unresolved sources, while in the UCD/CIT model "others" represents the sum of non-residential natural gas source combustion, aircraft emissions, and the sources that were not tagged in the current study. In general, the ranking and

PM0.1 and N10-1000 Source contributions in California
Figures 7-9 and 10-12 show the seasonal variation of major source contributions to primary N10 and PM0.1, respectively. The black circles in Figure 7-9 represent the measured N7-1000 at four BAAQMD sites in SFBA and six MATES sites in Los Angeles and Riverside counties. Predicted "best-fit" N10 follows the same trends as measured 420 seasonal variations of N7 at Livermore, Redwood City, Santa Rosa, Huntington Park, Inland Valley, and Rubidoux. The model over predicts N7 at Anaheim, central Los Angeles, and Compton but overall model performance statistics for N7 are within the target range for PM2.5 applications (see Table 1). Nucleation contributes to N10 at all sites but makes negligible contributions to PM0.1 concentrations. Traffic sources including 425 gasoline-and diesel-powered vehicles make significant contributions to PM0.1 concentrations at each measurement site depending on proximity to major freeways.
Near-roadway effects on ultrafine particle concentrations are not apparent since these locations were chosen to be regional monitors and so they are more than 300 m from the nearest freeway. Predicted contributions from traffic sources are consistent with the 430 molecular marker results illustrated in Figures 5-6. Traffic contributions to regional N10 concentrations more than 300 m away from roadways are even smaller than PM0.1 contributions because the size distribution of particles emitted from motor vehicles peaks at 100 -200 nm (Robert et al., 2007a;Robert et al., 2007b). Wood smoke makes strong contributions to regional PM0.1 concentrations in central California during winter but 435 much smaller contributions in the SoCAB because wood burning is not typically used for home heating in this region. Wood burning contributions to N10 are less dominant in central California because the size distribution of particles emitted from wood combustion peaks at 100-300 nm . The largest primary source of N10 in central California and N10+PM0.1 in the SoCAB is natural gas combustion.

440
Industrial processes and power generation that use natural gas do not follow strong seasonal cycles and so the strength of the natural gas source contributions is somewhat constant across seasons subject to variability caused by meteorological conditions.  contributions to "best-fit" N10 at Anaheim, Central LA, and Compton, respectively.
Results within each month have daily time resolution. 450 Figure 9: Seasonal variation of measured N7 (black circles) and major source contributions to "best-fit" N10 at Huntington, Inland-Valley, and Rubidoux, respectively.

455
Results within each month have daily time resolution. natural gas combustion makes the largest predicted primary contribution to N10 at all the sites that were evaluated. Traditional sources that were tracked including meat cooking, wood smoke, and mobile (gasoline + diesel) accounted for approximately 5-15% of the predicted N10 at the sites selected for study. "Other" sources that were not tagged explicitly in the current study accounted for 5-31% of N10 across these sites. Nucleation 475 is a significant source for of N10 for both BAAQMD sites and MATES sites where sulfur emissions were highest, with contributions ranging from 24-57%.
The strong N10 contribution from natural gas combustion reflects the emitted particle size distribution combined with the ubiquitous use of this fuel in the SFBA and SoCAB regions. The chemical composition and size distribution information for non-residential 480 natural gas combustion emissions used in this study was measured by Hildemann (1991) and Li and Hopke (1993), respectively. Size distributions and volatility were further confirmed during on-going field studies conducted by the current authors (Xue et al., 2018a). The estimated non-residential natural gas combustion particle number and mass size distributions are shown in Figure S1 (left column). Clearly, the majority of particles 485 from non-residential natural gas combustion are typically found in diameters <0.05 µm, while particles emitted from other sources such as wood combustion tend to have slightly larger particle diameter (with lower number concentration per unit of emitted mass).
These natural gas particles grow through the condensation of SOA once in the atmosphere, but they still contribute strongly to N10 concentrations. PM0.1 at Santa Rosa. The different rankings of source contributions to N10 and PM0.1 can be explained by the comparison of particle number-size distribution and particle masssize distribution for the non-residential natural gas and wood burning sources at the four evaluated sites ( Figure S1). Particles emitted from non-residential natural gas combustion and wood burning have number distributions that peak at particle diameters of 0.016-500 0.025 µm and 0.025-0.04 µm, respectively. Non-residential natural gas combustion and wood burning mass distributions, however, peak at particle diameters of 0.025-0.04 µm and 0.10-0.16 µm, respectively. Figure 15-17 show diurnal variations of measured N7-1000 and predicted "best-fit" N10 averaged over days in August and December 2012. These months span the temperature 505 range typically experienced across the year in California. Measured N7-1000 diurnal patterns in August generally peak in the afternoon hours between 12-3pm with an optional morning peak around 6am. The main afternoon peak appears to be related to nucleation events while the smaller early-morning peak appears to be related to early morning human activity including natural gas combustion. The predicted "best-fit" N10 510 diurnal variations in August followed the same trends as measurements at six out of ten sites (Livermore, Anaheim, Compton, Huntington Park, Inland Valley, and Rubidoux).
The model failed to capture the mid-day nucleation event at Redwood City and Santa Rosa possibly due to missing SO2 sources in the emissions inventory upwind from these sites. The model overestimated mid-day peak values at Anaheim and central Los

515
Angeles. In December, the measured N7-1000 diurnal pattern were more distinctly bimodal with the first peak around 7:00-8:00am and the second peak in the evening at around 8pm. This pattern reflects both the emissions activity and the mixing status of the atmosphere throughout the day. The predicted "best-fit" N10 concentration follows this same pattern. Nucleation continues to play a role during winter but does not dominate to 520 the point that it produces a midday peak in N10 concentrations . Non-residential natural gas combustion is predicted to be the largest source of N10 during morning and evening peaks. The diurnal profiles of non-residential natural gas emissions are included in supplemental information ( Figure S2) along with the regional distribution of those emissions ( Figure S3). These diurnal variations of the natural gas combustion emissions 525 were obtained directly from the emissions inventory specified by the California Air Resources Board. Industrial natural gas combustion emissions peak during the daytime with lower values at night. Emissions from electricity generation powered by natural gas peak in the morning and evening. Commercial natural gas combustion emissions may either peak in the morning and evening or they may follow a uniform diurnal profile 530 depending on the specific source and location.      Units are kcount cm -3 .
The dominant factors resolved by these studies have been traffic, urban background, secondary aerosol, wood burning and nucleation (Sowlat et al., 2016;Gu et al., 2011;Ogulei et al., 2007;Kasumba et al., 2009;Wang et al., 2013;Yue et 635 al., 2008;Friend et al., 2013). Particles from natural gas combustion were not separately identified by PMF because they do not contain a unique chemical tracer. It is very likely that natural gas combustion particles are artificially lumped into another source (e.g. traffic) or part of the "urban background" signal identified in previous studies. Natural gas combustion is used extensively in California for electric power, industrial, 640 commercial and residential use (Table S6), and so it seems plausible that this source contributes to ambient UFP concentrations.
The current UFP predictions rely on source profile measurements for wood burning, food cooking, mobile sources, and non-residential natural gas combustion (Cooper, 1989;Harley et al., 1992;Hildemann et al., 1991a;Hildemann et al., 1991b;Houck and L. 645 C., 1989;Kleeman et al., 2008a;Kleeman et al., 2000;Robert et al., 2007b;Robert et al., 2007a;Schauer et al., 1999bSchauer et al., , a, 2001Schauer et al., , 2002aTaback, 1979). All of these size distributions were measured using appropriate instruments and methods by knowledgeable researchers, but some of these past studies were conducted more than a decade ago. Size distribution information for vehicles, natural gas, etc. have been added 650 to the supplemental information ( Figure S4). Changes in fuel composition and emissions control technology in the interim years may have altered the emitted size distributions.
New measurements of particle size distributions emitted from natural gas and biomethane combustion were made in parallel with the current project to confirm the source profile measurements from past studies (Xue et al., 2018a). The results of these measurements 655 are consistent with previous size distribution results (Li and Hopke, 1993).
California has tighter air pollution standards than many other regions in the United States due to the severe air quality problems that have historically occurred in the state.
California therefore has a unique combination of fuels and emissions control technology that may affect the mixture of sources that contribute to atmospheric ultrafine particle 660 concentrations. Venecek et al. (2018) recently used the UCD/CIT air quality model with the 2011 National Emissions inventory to calculate source contributions to PM0.1 in 39 major cities across the United States during peak summer photochemical smog episodes in the year 2010. The findings from this study show that natural gas combustion is a major source of ultrafine particles in the regional atmosphere over urban areas across the 665 United States. The public health questions associated with ultrafine particles emitted by natural gas combustion have wide-ranging implications. Similar levels of ultrafine particle concentrations will likely occur in other regions across the world that extensively use natural gas as a fuel source, although other sources of ultrafine particles may also make strong contributions depending on the total mix of fuels in each region.

670
Recent theories suggest that primary particulate matter composed of semi-volatile organic compounds may evaporate after release to the atmosphere, which may reduce ambient Nx. Measurements conducted in parallel with the current study confirmed that particles emitted from natural gas combustion in home appliances partially evaporated when diluted by a factor of 25 in clean air, but particles emitted from reciprocating engines did 675 not evaporate under the same conditions (Xue et al., 2018a). Future work should verify the accuracy of the size and composition distributions for all natural gas combustion sources given their apparent importance for predicted Nx.
Evidence from both toxicology and epidemiology will be required to assess the effect of UFPs on public health. It is essential to identify and quantify UFP sources based on both 680 mass (PM0.1) and Nx during this process (Friend et al., 2013). An accurate comparison of both PM0.1 and Nx exposure could lay the groundwork for specific assessment of health effects of UFPs and potentially more efficient control strategies for PM emission from major sources (Yue et al., 2008). Ideally, spatial exposure patterns for Nx, PM0.1, and PM2.5 will be sufficiently unique to separate their individual effects in epidemiological 685 studies. Regression statistics for different metrics were calculated by using all grid cells in the model domain of the current study. The correlations between the various particle metrics were: R 2 (PM2.5 vs. N10)=0.35, R 2 (PM2.5 vs. PM0.1)=0.63, R 2 (PM0.1 vs. N10)=0.75.
It seems likely that future epidemiological studies will be able to differentiate between the effects of PM2.5 and Nx based on the low R 2 value. The potential for comparisons 690 between PM2.5 and PM0.1 is less clear cut, but previous work helps understand what may be possible. Ostro et al. (2015) compared the associations between IHD mortality and PM2.5 vs. PM0.1 in the California Teachers Study (CTS) cohort. Associations between IHD mortality and the sum of PM2.5 mass (p-value=0.001) were stronger than associations between IHD mortality and the sum of PM0.1 mass (p-value=0.01) but 695 individual components of mass (EC, OC, Cu, etc) all had stronger associations with IHD mortality in the PM0.1 size fraction than the PM2.5 size fraction.
The current study focuses on outdoor exposure to UFPs that may be useful in future epidemiological studies. Indoor or in-vehicle exposure to UFPs can also be significant (Wallace and Ott, 2011;Rim et al., 2010;Bhangar et al., 2011;Weichenthal et al., 700 2015;Fruin et al., 2008) but characterizing these micro-environments is beyond the scope of the current manuscript.

Conclusions
The UCD/CIT regional chemical transport model has been updated with a nucleation algorithm and combined with the existing size-resolved source profiles of particlualte 705 matter emissions to predict regional source contributions to airborne particle number concentration (N10) and airborne particulate ultrafine mass (PM0.1). Predicted 24-hour average N10 follows the same trend as measured N7 at ten sites across California in summer (Aug) and winter (Dec). Predicted diurnal variation of N10 follows the same trend as measured concentrations at the majority of the evaluation sites in August and

710
December, but the results suggest that further refinement is needed for both primary emissions and nucleation algorithms. Predicted PM0.1 source contributions follow the same trends as PM0.1 source contributions measured in a molecular marker study at four sites across California in summer (August) and winter (December) months. Natural gas combustion is the largest primary source of regional N10 at all locations outside of the 715 immediate vicinity of other major combustion sources. Nucleation contributed strongly to particle number during both the summer and winter months. Traffic sources contributed to N10 but did not dominate over regions more than 300 m away from freeways. Combustion sources such as wood burning, food cooking, and mobile sources made stronger contributions to PM0.1 at heavily urbanized locations. Wood burning for 720 home heating had strong seasonal patterns with peak concentrations in winter while other sources contributed more consistently throughout the seasons. Nucleation made a negligible contribution to PM0.1 over the urban areas at the focus of the current study.
The current study identifies natural gas combustion as an important source of ultrafine particle number and mass concentrations in urban regions throughout California. The (1) makes any warranty, express or implied, with respect to the use of any information, apparatus, method, or process disclosed in this report, or (2) assumes any liabilities with respect to use, or 740 damages resulting from the use or inability to use, any information, apparatus, method, or process disclosed in this report.