Aerosol – computational fluid dynamics modeling of ultrafine and black carbon particle emission , dilution , and growth near roadways

Many studies have shown that on-road vehicle emissions are the dominant source of ultrafine particles (UFPs; diameter < 100 nm) in urban areas and nearroadway environments. In order to advance our knowledge on the complex interactions and competition among atmospheric dilution, dispersion, and dynamics of UFPs, an aerosol dynamics–computational fluid dynamics (CFD) coupled model is developed and validated against field measurements. A unique approach of applying periodic boundary conditions is proposed to model pollutant dispersion and dynamics in one unified domain from the tailpipe level to the ambient near-road environment. This approach significantly reduces the size of the computational domain, and therefore allows fast simulation of multiple scenarios. The model is validated against measured turbulent kinetic energy (TKE) and horizontal gradient of pollution concentrations perpendicular to a major highway. Through a model sensitivity analysis, the relative importance of individual aerosol dynamical processes on the total particle number concentration (N) and particle number–size distribution (PSD) near a highway is investigated. The results demonstrate that (1) coagulation has a negligible effect on N and particle growth, (2) binary homogeneous nucleation (BHN) of H2SO4–H2O is likely responsible for elevated N closest to the road, and (3) N and particle growth are very sensitive to the condensation of semi-volatile organics (SVOCs), particle dry deposition, and the interaction between these processes. The results also indicate that, without the proper treatment of the atmospheric boundary layer (i.e., its wind profile and turbulence quantities), the nucleation rate would be underestimated by a factor of 5 in the vehicle wake region due to overestimated dilution. Therefore, introducing atmospheric boundary layer (ABL) conditions to activity-based emission models may potentially improve their performance in estimating UFP traffic emissions.


Introduction
Many studies have shown that vehicle emissions are the dominant source of ultrafine particles (UFPs; diameter < 100 nm) in urban areas and near-roadway environments.For example, about 95 % of UFPs (diameter = 50 ∼ 100 nm) observed near a US freeway were apportioned to fresh vehicular emissions (Toner et al., 2008).In the more confined environment of a street canyon, over 99 % of particles in number were found to be below 300 nm and number concentrations of particles in this size range were found to be linearly correlated with the traffic volume (Kumar et al., 2008).Due to their small size and abundance in number, recent toxicological and epidemiological studies suggest a strong correlation between adverse health effects and personal exposure to UFPs (e.g., Brugge et al., 2007;Ruckerl et al., 2007;Valavanidis et al., 2008).A recent review study by Schlesinger (2007) pointed out that the health impacts of chemical constituents, such as sulfate, seem to be inconsistent across all epidemiological studies.Comparing epidemiological studies of heart rate variability in humans, Grahame (2009) suggests that differences in accuracy of exposure information for health-relevant emissions may explain conflicting study results.This has led to an urgent need to study the temporal and spatial variations of local traffic emission in the vicinity of roadways.
With the growing concern of adverse health effects from exposure to UFPs, the gradients of vehicle-emitted pollutants (such as CO, NO x , and UFPs) have been measured in the ambient atmosphere near roadways (e.g., Beckerman et al., 2008;Reponen et al., 2003;Pirjola et al., 2006;Zhu et al., 2009).For example, Zhu et al. (2009) found that elevated particle numbers decay exponentially on the downwind side of three different types of roadways with increasing distance and reach background levels within a few hundred meters.Karner et al. (2010) synthesized field measurements of nearroadway pollutants from over 40 monitoring studies and investigated the concentration-distance relationship.The variation of UFP concentrations near roadways among studies is likely affected by factors including meteorological conditions (wind speed, ambient temperature, relative humidity, and atmospheric stability), traffic characteristics (volume and fleet composition), the geometry of roadways, and aerosol transformation processes (nucleation, coagulation, condensation/evaporation, and dry deposition).However, field measurements alone are often associated with such limitations as low spatial or temporal resolution in sampling, conclusions restricted by local meteorology, and difficulties in separating the effects of interactive processes.
Therefore, numerical modeling of UFPs has been conducted to address these limitations.Due to the challenge of resolving processes with very different scales, a two-stage dilution modeling strategy, including tailpipe-to-road and roadto-ambient dilutions, has been proposed (Zhang and Wexler, 2004;Zhang et al., 2004).In the first stage (i.e., tailpipeto-road), strong vehicle-induced turbulence (VIT) results in fast and strong dilution (dilution ratio ∼ 1000 in 1 s) and triggers nucleation and condensation/evaporation.While in the road-to-ambient stage, atmospheric boundary layer turbulence (ABLT) continues to dilute exhaust particles with ambient air accompanied with particle size changes due to condensation/evaporation.A review study by Carpentieri et al. (2011) has shown that, with recent advances in numerical modeling, computational fluid dynamics (CFD) models can be valuable tools for nanoparticle dispersion in the first stage of dilution.In addition to the limited spatial scale of the dispersion investigated, other limitations in these most recent modeling studies include particle number-size distribution (PSD) and chemical composition not being explicitly resolved (Chan et al., 2010).Recent modeling studies of UFP dispersion on street level, on the other hand, have crudely simplified treatment of vehicular emission, VIT, and aerosol dynamics (Gidhagen et al., 2003(Gidhagen et al., , 2004a;;Kumar et al., 2009).
Most recently, Wang et al. (2013) proposed a two-stage simulation approach to integrate the tailpipe-to-road dispersion into the road-to-ambient dispersion stage for the first time.As the authors noted, however, the proposed approach remains computationally demanding, especially when both particle size and chemical composition need to be resolved.To effectively model UFP dynamics and dispersion near roadways in a single unified tailpipe-to-ambient domain, a unique approach of applying periodic boundary conditions to the computational domain is proposed in this paper.Compared to a road-to-ambient dispersion modeling approach, the advantage of a unified domain is that the uncertainty due to a simplified or non-existent treatment of VIT can be greatly reduced by explicitly modeling VIT.With VIT being explicitly modeled, aerosol dynamics (such as nucleation, condensation, and evaporation) triggered by the rapid firststage dilution can be properly incorporated into dispersion models to study their effects on roadside air quality.From a modeling perspective, such a unified model provides a tool to link individual tailpipe emissions (controlled laboratory measurements) to roadside air quality (ambient field measurements), which is a noted challenging task (Keskinen and Ronkko, 2010).As the main focus of this paper, we present the development and validation of a multi-component sectional aerosol dynamics-CFD coupled model to account for the complex dilution, dispersion, and dynamics of UFPs immediately after tailpipe emission to ambient background.

Multi-phase approach to the mixture of atmospheric gas and aerosol
Following previous studies (e.g., Li et al., 2006b;Uhrner et al., 2007;Wang et al., 2013;Albriet et al., 2010), the commercial CFD code ANSYS Fluent is used to model turbulent flow around realistically shaped vehicles.An Euler-Euler approach to multi-phase atmospheric air flow is coupled with new particle formation, transformation, and dry deposition processes.Among the general multi-phase models available in Fluent, the mixture multi-phase model is chosen for the present work due to its superior numerical efficiency compared to the Eulerian multi-phase model (ANSYS, 2009a).The concept of phase in the Fluent multi-phase model is defined in a broad sense as an identifiable class of material that has a particular inertial response to and interaction with the flow (ANSYS, 2009b).The gas phase and the particulate phase in the real atmosphere are represented in the model by the primary phase and a number of secondary phases, respectively.The number of secondary phases is determined by the number of discrete size bins used to resolve the particle number-size distribution.The flow field of the mixture phase is obtained by numerically solving Reynoldsaveraged Navier-Stokes equations (RANS) including con-servation equations of mass, momentum and energy with the standard k-ε turbulence model.For the transport of gas and particulate species, Fluent predicts the local mass fraction of each species by solving the advection-diffusion equations.The volume fraction of each secondary phase in a control volume (equivalent to the number concentration of particles of the same size) is obtained by numerically solving the continuity equation for the secondary phase with a specified source term due to aerosol dynamic processes.The diffusive mass flux in Fluent is modeled as the sum of two components: molecular and turbulent diffusion (e.g., Eq. 4 in Di Sabatino et al., 2007).Turbulent diffusion due to VIT and ABLT are the main dilution mechanisms for pollutants in the near-road environment (Zhang and Wexler, 2004).The key parameter governing modeled turbulent diffusion of pollutants using the RANS approach is the turbulent Schmidt number (Sc t ), which is defined as the ratio of the turbulent momentum diffusivity and the turbulent mass diffusivity.Analyzing a widely distributed range of Sc t (0.2-1.3) in literature versus the commonly used values (0.7-0.9), Tominaga and Stathopoulos (2007) found that, for plume dispersion in open country for example, a smaller value of Sc t might be used to compensate the underestimated turbulent momentum diffusion.They further suggested adopting its standard value in simulations without this type of underestimation.Therefore, given the successful model validation on turbulent kinetic energy (TKE) (discussed in Sect.4.1.1),the standard value of 0.7 for Sc t is used in our study.Thus, the advection, the turbulent mixing, and the diffusion of gases and particles are inherently treated by Fluent through the continuity equation for each phase.Aerosol dynamic processes, which change the chemical components in particle and gas phases, are integrated through the source terms in continuity equations, and incorporated into Fluent through user-defined functions (UDF).

Aerosol dynamics
Each secondary phase is a particulate phase composed of mixed chemical components within a specified size range.The density of particles in a given size bin is dynamically computed by Fluent based on the volume-weighted mixing law.From the continuity equation for each secondary phase p, the volume fraction of the secondary phase (α p ) is obtained by solving where u is the velocity field of the primary (gas) phase, u dr,p is the drift velocity of secondary phase, S p = M i=1 S p,i is the rate of mass transfer for phase p, S p,i is the rate of mass transfer for species i in phase p, and M is the total number of chemical species in the model.For phase p, the local mass fraction of each species (Y p,i ) is predicted by solving a convection-diffusion equation for the ith species, given S p,i due to the aerosol dynamical processes described in Sect.2.2.1-2.2.4.The total mass of the gas and particulate phases is conserved, while the particle number is diagnosed from the predicted mass.The number concentration of particles in size bin p (N p , in particles/cm 3 ) is computed from the ratio of the phase volume fraction solved by Fluent to the particle volume of a certain size: where D p is the diameter (in m) for particles in size bin p.Additionally, the local mass concentration of chemical component i from particles in size bin p is calculated from the phase volume fraction (α p ) and the local mass fraction (Y p,i ) as The underlying implementation of aerosol dynamics is a multi-component, size-resolved, sectional aerosol model, as described as follows.

Nucleation
Immediately after tailpipe emissions, new particles form by homogeneous nucleation with initial particle size around 1.5-2.0nm in the first few milliseconds of exhaust cooling and dilution (Kulmala et al., 2007).A qualitative investigation by Zhang and Wexler (2004) found that sulfuric acid-induced nucleation could be the dominant new particle production process.The experimental study conducted by Arnold et al. (2006) observed a positive correlation between gaseous sulfuric acid and particle number in the exhaust of a passenger diesel car burning ultra-low sulfur fuel, indicting an important role of sulfuric acid-induced nucleation.The sulfuric acid gas emission rate is estimated based on fuel sulfur content following Uhrner et al. (2007).The parameterization of binary homogeneous nucleation (BHN) of H 2 SO 4 -H 2 O (Vehkamaki et al., 2003) developed specifically for engine exhaust dilution conditions is implemented in this study.This parameterization has already been successfully used in a number of different aerosol-CFD applications (e.g., Uhrner et al., 2007Uhrner et al., , 2011;;Albriet et al., 2010;Wang and Zhang, 2012).

Coagulation
Particles in the exhaust plume collide due to random (Brownian) motion and turbulent mixing to form larger particles, which is called coagulation.The coagulation process reduces N (mainly in the smaller size range) while preserving the aerosol total mass.However, it modifies the particle numbersize distribution, and internally mixes particles of different chemical composition over the population.Coagulation may be driven by Brownian motion, turbulent flow conditions, gravitational collection, inertial motion, and turbulent shear.Individual coagulation rate coefficients (or coagulation kernels) due to the above driving forces are calculated in this work based upon Jacobson (2005), with consideration of particle flow regimes and convective Brownian diffusion enhancement.The overall coagulation rate coefficient is the summation of individual coefficients.

Condensation and evaporation
A complex mixture of condensable gases, including water vapor, sulfuric acid, and semi-volatile organics (SVOCs), is emitted from the tailpipe after fuel combustion in the engine.During the strong dilution and cooling stage of up to a few seconds after emission (Zhang and Wexler, 2004), supersaturation of these condensable gases occurs and favors the diffusion-limited mass transfer process from the gas phase to the pre-existing particle phase.Following primary emission and nucleation, condensation of SVOCs was suggested by a number of studies (e.g., Clements et al., 2009;Wang and Zhang, 2012;Albriet et al., 2010;Uhrner et al., 2011;Mathis et al., 2004) to be responsible for the rapid growth of nanoparticles in the exhaust plume.As the reverse of condensation, evaporation occurs due to further dilution of condensable gases to subsaturation level in the air surrounding exhaust particles.It was suggested by field measurements of freeway emissions from predominantly gasoline vehicles that lower ambient temperature may favor the condensation of organic species to the particle phase (Kuhn et al., 2005).In this study, the net mass transfer rate of a condensable gas from/to an existing particle with multiple components is driven by the difference between the bulk partial pressure and the saturation vapor pressure above the particle surface (Jacobson, 2005).The calculation of species mass transfer rate implements corrections to the diffusion coefficient, the thermal conductivity of air, and the saturation vapor pressure over curved particle surfaces to reflect its dependence on particle size and chemical composition.

Dry deposition
Driven by mechanisms such as Brownian diffusion, turbulent diffusion, sedimentation, and advection, dry deposition removes particles at the air-surface interface when they contact and remain on the surface (Jacobson, 2005).Brownian diffusion is more effective in removing smaller particles due to their larger diffusion coefficient, while sedimentation is more important for larger particles whose fall speeds are much higher.In the current study, parameterization of particle dry deposition follows the size-resolved dry deposition scheme developed by Zhang et al. (2001).The effect of turbulent mixing on particle dry deposition is taken into account by the locally calculated friction velocity.This parameterization has been successfully validated and implemented in a number of air quality and climate studies (e.g., Gong et al., 2003;Pye and Seinfeld, 2010), and it has recently been improved and extended (Petroff and Zhang, 2010).The recent development accounts for more detailed characteristics of the surface canopy, and suggests possible overestimation of dry deposition velocity for particles in the fine mode.Thus, our current study is likely biased to overestimate the removal of UFPs by dry deposition.

Modeling turbulence
For turbulence modeling, although the large eddy simulation (LES) approach has been reported to be a more promising solution, the standard k-ε turbulence model is implemented in this work for several reasons.First, the high computational demands of the LES approach prevent its application for modeling the dispersion and transformation of multiple pollutants with complex geometry (i.e., gas and particle emissions and aerosol dynamics from multiple vehicles in this study).Compared to RANS closures, the LES approach is at least 1 order of magnitude more computationally expensive (Rodi, 1997).Second, a proper treatment of the atmospheric boundary layer (ABL) has proven to be crucial to dispersion modeling studies (Blocken et al., 2007b;Zhang, 1994).Recent advances made by Balogh et al. (2012) and Parente et al. (2011a, b) permit a general and practical means to include the ABL using the standard k-ε model.To achieve this in LES simulations, on the other hand, inflow conditions would have to be carefully generated with additional, significant computational overhead (Xie and Castro, 2008;Li et al., 2006b).Finally, RANS models agree reasonably well with experimental data in predicting mean flow and pollutant concentrations (e.g., Labovsky and Jelemensky, 2011;Sklavounos and Rigas, 2004).Kim et al. (2001) successfully modeled the dispersion of a truck exhaust plume in a wind tunnel using k-ε turbulent closure focusing on rapid dilution and turbulent mixing of exhaust CO 2 .Although CFD codes have been widely adopted in pollutant dispersion modeling, the accuracy of such simulations can be seriously compromised when wall functions based on experimental data for sand-grain roughened pipes and channels are applied at the bottom of the computational domain (Blocken et al., 2007b;Riddle et al., 2004).Attempts have been made to better predict ABL flow by changing turbulent model constants, tuning boundary profiles, and modifying wall functions and turbulent transport equations (Blocken et al., 2007a;Pontiggia et al., 2009;Li et al., 2006a;Alinot and Masson, 2005;Hargreaves and Wright, 2007;Labovsky and Jelemensky, 2011;Balogh et al., 2012;Parente et al., 2011a, b).Among the aforementioned studies, recent advances made by Balogh et al. (2012) and Parente et al. (2011a, b) are implemented in this work, which permits a general and practical means to include ABL using the standard k-ε model in the CFD code, Fluent.The ABL profiles of mean velocity, TKE, and dissipation rate for atmospheric flow under neutral stratification conditions (Richards A modified wall function for turbulent mean velocity following Parente et al. ( 2011b) is implemented through UDF and applied to wall adjacent cells, where E = υ z 0 u * and z + = (z+z 0 )u * υ .To keep the default constant value of σ ε in the standard k-ε model, a source term is added to the dissipation rate equation as follows: Furthermore, we adopt the approach by Parente et al. (2011a) to allow a gradual transition from Eqs. (4-6) (i.e the undisturbed ABL) to the wake region simulated by the standard k-ε model (Supplement, Fig. S4).

FEVER field study
The Fast Evolution of Vehicle Emissions from Roadway (FEVER) study was conducted to monitor pollutant gradients perpendicular to a major highway north of Toronto, Canada (Highway 400, Hwy 400; 43.994 • N, 79.583 • W).The model developed and tested in this paper was designed to simulate the FEVER observations.A complete description of the monitoring strategies of the FEVER project was documented by Gordon et al. (2012a, b), the BC emission rate for gasoline vehicles was estimated by Liggio et al. (2012), and the rapid organic aerosol production under intense solar radiation was investigated by Stroud et al. (2014).
The site under investigation is a six-lane (25 m across from the lane edges) highway, mainly surrounded by flat agricultural fields and some trees lining the side roads, with negligible local pollution sources other than vehicular emissions.To validate modeled VIT, the on-road TKE data measured by the Canadian Regional and Urban Investigation System for Environmental Research (CRUISER) mobile laboratory were compared with modeled TKE.The on-road TKE data were measured by two 3-D sonic anemometers during passenger vehicle chasing experiments on 6 days between 20 August and 15 September 2010.To validate modeled nearroad dispersion, a case study period of 14 and 15 September 2010 between 05:00 and 08:00 LT was chosen for comparison.The near-road TKE data were measured by a 3-D sonic anemometer at a 3 m tower located 22 m east of the road center.Wind speed and direction data were measured by an Air-Pointer system (Recordum GmbH), averaged every minute, 34 m east of the road center.As shown in Table 1, the predominant wind direction was approximately perpendicular to the highway and the median Monin-Obukhov length indicates near-neutral stability conditions.The CRUISER mobile lab housed instrumentation to measure BC, CO 2 , and UFPs while driving transects perpendicular to the highway.Following a previous study (Gordon et al., 2012a), data were filtered for winds within 45 • of the highway normal, which results in removing less than 5 % of the data.In addition, particle size distributions between 14.6 and 661.2 nm were measured at two fixed sites with scanning mobility particle sizer spectrometers (SMPS) every 3 min and averaged for 05:00-06:00 and 06:00-08:00 LT of 14 and 15 September 2010 for model validation.

Computational domain and flow boundary conditions
The sizes of computational domain used for near-road dispersion modeling in this study are summarized in Table 2.The computation domain for the base case simulation, for example, is shown in Fig. 1a.The top of the domain is set to 50 m above the ground so that the turbulent flow near the surface is not affected by the top boundary (the x-y plane in purple mesh).The horizontal dimension of 375 m perpendicular to the highway (x axis) is determined by considering the availability of measurements and the extent of pollutant dispersion (Gordon et al., 2012a).Both dimensions of the domain are in compliance with the recommendations for CFD simulation of flows in the urban environment (Franke, 2007).To reduce computational overhead, first, the actual two-way, six-lane traffic fleet is represented by a one-way, three-lane traffic fleet (as shown in Fig. 1b) while conserving the total traffic volume; furthermore, translational periodic boundary conditions are applied to the x-z planes (in cyan mesh) to account for the effects of the stable, continuous traffic fleet on the highway by numerically repeating the computational domain in the direction of y axis.Based  on the measured traffic volume of about 104.3 passenger vehicles min −1 , traveling speed of approximately 120 km h −1 (or 33.3 m s −1 ), and the assumed average vehicle length of 4.5 m, the average y axis distance (bumper to bumper) between two vehicles traveling in adjacent lanes is calculated as 11.5 m assuming all three lanes are evenly occupied.Thus, the horizontal dimension of 48 m along the highway (y axis) is calculated based on the measured traffic volume for weekday early morning rush hours between 06:00 and 08:00 LT.The whole domain is meshed into 871 065 unstructured hexahedral cells with the finest ones concentrated around the moving vehicles, tailpipes, and their wake regions and immaterially above the ground.
In our simulation, the vehicles are set to be stationary, nonslip walls with a roughness length of 0.0015 m (Wang and Zhang, 2009), while the blowing air has two velocity components: the first component towards the vehicles (or the negative y axis direction) of 33.3 m s −1 in magnitude, and the second component perpendicular to the highway (or the negative x axis direction) in a form given by Eq. ( 4).The first component of wind velocity accounts for the relative movement between the moving vehicles and still air, and the second component describes the observed wind speed according to a fully developed ABL wind profile under neutral stratification.Thus, the upwind side boundary parallel to the road (the x-y plane in purple mesh) and the top boundary are set as velocity inlets.The ground surface is set to have the same velocity magnitude as the running vehicles but in the opposite direction.The modified wall function (Eq.7) and the additional source term to dissipation rate (Eq.8) as described in Sect.2.3 are applied to the ground surface to account for a fully developed ABL turbulent flow.Translational periodic boundary conditions are set to the x-z planes of the domain, and a pressure-outlet boundary is applied to the boundary (in red mesh) at the far-end side to the highway.Each tailpipe, 52 mm in diameter, is specified as mass flow inlet with a mass flow rate of 0.055 kg s −1 and an exhaust temperature of 480 K (Uhrner et al., 2007).An O-grid composed of seven inflation layers, whose thickness gradually increased from 0.003 m, around vehicles is used to allow the standard wall functions to apply to the fully turbulent layer around moving vehicles.Crucial but not mentioned in previous CFD models, a high spatial resolution applied here results in a dimensionless wall distance (y + = ρµ t y/µ) of about 90 at the vehicle surface, which is well within the suggested range of 30 to 300 for the standard wall functions to apply (ANSYS, 2009b).Simulation results of turbulence and pollutant concentrations at site B (33 m east of the highway center) show no significant dependence on a further refined grid (Supplement, Fig. S2).

Chemical boundary conditions: background concentrations and traffic emissions
In addition to meteorological and traffic data, chemical data of gases and particles are required as part of CFD boundary conditions.According to the source type, the required chemical data are divided into two categories: background concentrations and traffic emission rates.The mass concentrations of background gaseous and particulate species from the FEVER field measurements are listed in Table 3, with their corresponding values used as model input.The background gas phase includes dry air (O 2 and N 2 ), water vapor, and CO 2 .The CO 2 volume fraction and relative humidity (RH) value are converted to mass fractions to specify the species input values for the velocity-inlet boundaries.For the particulate phase, a bi-modal lognormal particle size distribution is assumed (as summarized in Table 5).The parameters of background particles are obtained from the FEVER measurements at about 100 m upwind of Hwy 400.Given the total background N and PSD, the volume fractions of individual size bins are obtained, and the mass of black carbon (BC) and organic aerosol (OA) is distributed into each size bin according to the ratio of their background mass concentrations listed in Table 3. Vehicles driving on the highway continuously emit a complex mixture of gases and particles.It is not possible to include a complete set of gaseous and particulate species in the model, which also is not numerically practical.In this study, the tailpipe emission rates of the gaseous and particulate species are summarized in Table 4. Currently, the exhaust gas is composed of CO 2 , H 2 O, H 2 SO 4 , SVOCs, and N 2 , which are key species to the aerosol dynamics and dispersion.The treatment of H 2 SO 4 as direct emission rather than a mixture of SO 2 and SO 3 followed by hydrolysis has been explained in Sect.2.2.1.It was suggested by modeling single exhaust plumes (e.g., Albriet et al., 2010;Uhrner et al., 2007Uhrner et al., , 2011) ) that SVOCs are likely to be responsible for the rapid growth in particle size when they condense on UFPs.Following Albriet et al. (2010), pyrene (C 16 H 10 ), N nonadecane (C 19 H 40 ), and N pentacosane (C 25 H 52 ) are introduced to represent the polycyclic semi-volatile organic compounds, the semi-volatile alkanes between C 14 and C 22 , and semi-volatile alkanes between C 23 and C 29 , respectively.The mass fractions of the above three groups of SVOCs are based on Albriet et al. (2010), and the total mass emission rate of SVOCs is set as 0.0186 g km −1 (Pye and Seinfeld, 2010).All SVOCs initially from the tailpipe are assumed to exist only in the gas phase, but are subject to interactions with the particle phase through condensation/evaporation upon immediate dilution with the surrounding air.To reduce the number of species considered in the model, the non-volatile fraction of primary organic aerosol (POA) from tailpipes is assumed to share the properties of the background OA, i.e., with an average molecular mass of 300 g mol −1 and an average density of 1.5 g cm −3 .This assumption is not likely to affect our results because the amount of the non-volatile fraction of POA from tailpipes is very small compared to the background OA.N and PSD for tailpipe emissions are based on a recent study by Nikolova et al. (2011b), which provides an emission rate according to traffic volume and type.As pointed out by Nikolova et al. (2011a), however, the parameterization they originally proposed implicitly accounts for a fast nucleation process.As indicated by laboratory measurements (Ronkko et al., 2007;Kirchner et al., 2009), the nucleation mode particles have a nonvolatile core in the exhaust of a heavy-duty diesel vehicle; however, they are completely volatile under 280 • C in the exhaust of a diesel passenger car.Thus, we assume in this study that N of nucleation mode particles from all passenger cars are from BHN, while those from heavy-duty vehicles (HDVs) have a solid core of BC and nonvolatile POA.Given the mass flow rate of the tailpipe exhaust, the mass fraction of each individual species can be estimated from its mass emission rate listed in Table 4.These mass fractions are used to specify chemical boundary conditions for tailpipes.

Results and discussions
Turbulent mixing of tailpipe emissions with the ambient air largely determines the initial dilution and the threedimensional distribution of the traffic pollutants downwind of roadways.Thus, the modeled TKE is first compared against on-road and near-road TKE measurements reported by Gordon et al. (2012b).For model validation results in Sect.4.1, two scenarios with different traffic conditions (base case: 06:00-08:00 LT and half traffic case: 05:00-06:00 LT) are considered.The modeled CO 2 and BC concentrations and PSDs are compared with the FEVER field measurements.Finally, the impacts of individual aerosol dynamical processes on UFPs and model sensitivity to the treatment of ABL are investigated in Sect.4.2.A total of five sensitivity runs are performed based upon the base case (06:00-08:00 LT).Four sensitivity runs are conducted by turning off a single aerosol dynamical process for each run, and the results are compared with the base case in Table 6.An additional sensitivity run is conducted without maintaining modeled ABL profiles through Eqs.(7-8).

Turbulent kinetic energy
Both theoretical (e.g., Zhang and Wexler, 2004;Ketzel and Berkowicz, 2004) and monitoring (e.g., Zhou and Levy, 2007) studies have concluded that dilution is the dominant mechanism governing N of UFPs, and turbulent kinetic energy (TKE) measures the strength of mixing and dilution.Figure 2 compares the modeled TKE with the measurements from the FEVER chasing experiments of passenger cars on highways.The x axis in Fig. 2     passenger cars and the y axis is the 10 s average TKE.The modeled TKE values in Fig. 2 are calculated for individual vehicle wakes super-imposed on an estimated background on-road TKE of 2.4 m 2 /s 2 (Gordon et al., 2012b).The average traveling speed during the FEVER chasing experiments was about 20 m s −1 , and 84.5 % of the measurements were taken at a chasing speed between 15 and 25 m s −1 .Therefore, model simulations are conducted for passenger cars traveling at 15, 20, and 25 m s −1 .As shown in Fig. 2, the modeled TKE in the wake of a vehicle traveling at a speed of 15-25 m s −1 agrees well within the 25th and 75th percentile of the measurements; furthermore, the variations among modeled TKE in Fig. 2 show the sensitivity of on-road TKE to vehicle type and traveling speed.With these scenarios agreeing within the 25th and 75th percentile of the measurements, it is clear that the turbulent mixing within individual vehicle wakes on the highway can be reasonably well modeled.As a measure of the turbulent dilution under perpendicular wind conditions, near-road TKE is also modeled and compared to the measurements at a tower located 22 m east of the road center.Two time periods (05:00-06:00 and 06:00-08:00 LT) with distinctly different traffic volumes are considered.For both time periods, the atmospheric boundary layer was neutrally stratified.The average traffic volume, however, increased from 54.9 to 104.3 vehicles per minute, as shown in Fig. 2 of Gordon et al. (2012a).Although the tower is stationary, the distance and time along the wind trajectory from the highway center vary with changing wind direction and speed.Therefore, the evolution of turbulence with distance can be   investigated based on measurements at a fixed location in a Lagrangian sense.There are 120 and 240 measurements of TKE taken for the periods of 05:00-06:00 and 06:00-08:00 LT, respectively, and they are binned by distance as shown in Fig. 3.The obtained wind trajectory distances vary between 20 and 80 m with about 93 % of them concentrating on the first bin (20-40 m).For the period of 05:00-06:00 LT, the measured (modeled) TKE at a distance of 20-40 m from the highway center is in the range of 0.46-0.80(0.58-0.73) m 2 s −2 .Similarly, for the period of 06:00-08:00 LT, the measured (modeled) TKE lies in the range of 0.55-0.90(0.65-0.95) m 2 s −2 .Although the observed TKE decay is limited in spatial resolution for both time periods, the comparison in Fig. 3 shows an adequate agreement with the field measurements and suggests turbulent mixing in a roadside environment can be successfully modeled even with varying traffic volumes.

Near-road concentration gradients: CO 2 and BC
As a chemically passive gas species in vehicular emissions, CO 2 is an ideal indicator of atmospheric mixing of tailpipe exhaust with ambient air.In a previous study (Kim et al., 2001), CO 2 was experimentally measured inside a single turbulent plume of heavy-duty truck exhaust and successfully modeled with the standard k-ε model in the CFD code Fluent.The focus of this study, however, is the horizontal concentration gradient on the downwind side of a highway.Figure 4a shows the concentration of CO 2 (ppmv) as a function of downwind distance from the center of Hwy 400 for the morning period of 06:00-08:00 LT.FEVER measurements were first corrected to wind trajectory distance, grouped into 20 m bins between 50 and 350 m, and then plotted in median concentrations and 25th and 75th percentiles.Modeled CO 2 concentrations closely follow the decreasing trend of the median values of the FEVER measurements, and agree well within the 25th and 75th percentiles.However, the model tends to underestimate CO 2 concentrations by about 6 ppmv in the first 50 m (i.e., 50-100 m) of downwind distance and overestimate by about 8 ppmv in the last 50 m (i.e., 250-300 m).Similarly, the concentration-distance relationship for particulate BC is shown in Fig. 4b.Modeled BC concentrations are also within the 25th and 75th percentiles exhibiting a trend with distance similar to the median of the measured values.Similar to CO 2 , minor underestimations (15 %) between 50 and 100 m from Hwy 400 and slight overestimations (20 %) after 100 m were observed for particulate BC.
The behavior of the model suggests slightly overestimated pollutant concentrations between 0 and 3 m above the ground within 100 m distance to the road, possibly due to underesti-  mated vertical mixing by the model.Beyond 100 m from the road, vertical diffusion of the near-surface pollutants results in the overestimations.There are two factors that may explain the underestimated vertical mixing closest to the road.First, the modeled road structure is missing a 1 m high barrier at the highway center, which could potentially lift nearsurface pollutants under cross-wind conditions (Ning et al., 2010;Hagler et al., 2011).Second, midsize and heavy-duty trucks are neglected due to their small fractions in total traffic, which emit pollutants at a greater height (up to 4 m) than (0.5 m) passenger cars (Gordon et al., 2012a).

Particle size distribution
The fate of atmospheric particles depends strongly on PSD, which is the result of the complex influences of mobile emissions, atmospheric dilution, and transformation processes.In Fig. 5, the predicted PSDs are compared to SMPS measurements for two selected periods of the early morning rush hours (i.e., 05:00-06:00 and 06:00-08:00 LT) at two fixed locations: sites B and C, located 34 and 300 m from Hwy 400, respectively.The early morning rush hours are subdivided into the above two periods based on hourly averaged traffic flow.According to Gordon et al. (2012a), the traffic flow of 06:00-08:00 LT (105 veh min −1 ) almost doubled the average traffic of 05:00-06:00 LT (55 veh min −1 ), while the ambient conditions (such as atmospheric stability class, in-coming solar radiation, and wind velocity) remained approximately constant.
Comparing the measured PSDs at two sites, it was found that, for all measured particle sizes, the number concentrations decreased significantly when particles were transported from 34 to 300 m downwind of the highway.The observed total particle number concentrations decreased by a factor of about 2.5 and 5.7 between these two locations for the periods of 05:00-06:00 and 06:00-08:00 LT, respectively.Similar to many previous roadside monitoring studies reviewed by Pant and Harrison (2013), the measured PSDs showed distinct multi-modal size regimes (i.e., nucleation, soot, and accumulation modes) in the measured size range of 15-700 nm.Tri-modal lognormal curve fitting of the observed PSDs revealed a nucleation mode at 20-25 nm, a soot mode at 65-75 nm, and an accumulation mode at 160-380 nm.In agreement with Zhu et al. (2002a, b), the nucleation mode particles dominated N and decreased much faster (by a factor of 8 in this study) than that of the accumulation particles (only a factor of 2 in this study).It was also found that the geometric mean diameter of the nucleation mode increased from 20.0 (at site B) to 24.7 nm (at site C), which may be attributed to the condensation of SVOCs (Clements et al., 2009) and the coagulation of nucleation mode particles (Zhu et al., 2002b).
The comparison in Fig. 5 demonstrates an adequate agreement between the modeled and the observed PSD at both distances under different traffic conditions.For the peak traffic hours during 06:00-08:00 LT, the model estimated total particle number concentrations are 6.25 × 10 4 and 1.66 × 10 4 particles cm −3 at sites B and C, respectively, with approximately 10 % underestimations compared to the observations.Second, the dominant nucleation mode was properly captured by the current model, as well as its decreasing trend with increasing distance from the highway.Furthermore, the nucleation mode particles were modeled to grow from 19.8 to 24.3 nm in geometric mean diameter with increasing distance away from the highway.This agrees exceptionally well with the observations.Similar conclusions can be drawn from the 05:00-06:00 LT comparison.
However, the model clearly underestimated the number concentrations of particles of 100-730 nm in diameter.This discrepancy can be attributed, at least partially, to the missing non-tailpipe emissions in the current model, such as brake wear, road-tire interaction, and re-suspension of road dust as reviewed by Kumar et al. (2013).Although road dust particles formed mechanically by frictional contact between road surface and tire or between break system components are assumed to be primarily coarse particles, both laboratory experiments (Dahl et al., 2006;Gustafsson et al., 2008) and real-world measurements (Mathissen et al., 2011) recently observed a significant portion of particles to be in the range of 6-700 nm in diameter.On the other hand, the estimated emission factors for sub-micrometer particles generated by the road-tire interaction under steady driving condition based on these studies vary significantly, indicating that the emis-sion strength tends to be very site specific.Thus, the underestimated particles larger than 100 nm might be a result of missing estimates of non-tailpipe emissions for the underlying site.

Role of aerosol dynamical processes
Along with dilution, aerosol dynamical processes (i.e., condensation/evaporation, coagulation, nucleation, and dry deposition) may interact with one another and modify N and PSD in near-road environments.In this section, the relative importance of the above aerosol dynamical processes is investigated by conducting simulations with individual processes removed and comparing against the base case simulation, in which all dynamical processes are considered by the model.The obtained N and the geometric mean diameter of nucleation mode particles from this sensitivity analysis are summarized in Table 6.The base case simulation demonstrates that in moving from site B to site C, the nucleation mode particles decrease by approximately a factor of 3, and the geometric mean diameter increases by 4.5 nm.The modeled soot mode and accumulation mode particles are excluded from the analysis due to significant underestimations compared to the measurements, as discussed in the previous section.The results of excluding particle dry deposition process are investigated first because its impact on N can interact with particle nucleation and condensation processes, as discussed later in this section.
When particle dry deposition is deactivated in the model, nucleation mode particle numbers increase significantly (∼ 23 and 53 % at site B and C, respectively), resulting in 1-5 nm smaller geometric mean diameters compared to the base case as listed in Table 6.The modeled particle dry deposition velocity is up to 0.2 m s −1 for the smallest particles of 3-5 nm in diameter due to strong Brownian diffusion.Our results show that particle dry deposition plays a significant role in governing N in the vicinity of roadways between 30 and 300 m.Gidhagen et al. (2004b) estimated that dry deposition removes only about 12 % of total particles near a Swedish highway, in contrast to our estimation of 15-35 %.This discrepancy may be due to the different treatment of atmospheric boundary layer turbulence in both studies.Gidhagen et al. (2004b) introduced an artificial source of turbulence into their model to mimic the observed atmospheric dilution of NO x near the road, while the theoretically based method by Parente et al. (2011b) combined with the measured ABLT was implemented in this study.The modeled VIT and ABLT have been validated against the measurements in Sect.4.1.1.
Compared with the base case simulation, the predicted mean diameters at both sites without condensation remain nearly unchanged from the tailpipe emission of 15 nm.At the same time, the predicted particle number concentrations without condensation are about 1 order of magnitude lower than the base case, and are the lowest among all scenarios.The implication of this is twofold.In agreement with previous modeling studies (i.e., Wang and Zhang, 2012;Uhrner et al., 2007;Uhrner et al., 2011;Albriet et al., 2010), it strongly suggests that the condensation of SVOCs is responsible for the growth of nucleation mode particles during their atmospheric transport.It also reflects the strong interactions between particle growth and removal processes in governing the simulated particle number.Without the condensational growth of nucleation mode particles, new particles formed due to the BHN mechanism remain in the smallest size bin of 3-5 nm in diameter.Immediately after formation, these particles are subject to efficient removal by dry deposition due to their small particle sizes, resulting in the lowest N among scenarios.This result implies that controlling tailpipe SVOC emissions may indirectly help reduce UFP number concentrations in the vicinity of roadways.
For the scenario without BHN of H 2 SO 4 -H 2 O, the geometric mean diameters at both sites are similar to the base case with slightly more condensational growth in size.However, the particle number concentrations are underestimated by 36 and 10 % at site B and C, respectively, compared to the base case.This implies that over 60 and 90 % of the nucleation mode particles at site B and C are attributed to HDV emissions with non-volatile cores.The result also shows that the BHN of H 2 SO 4 -H 2 O has the greatest impact on the particle population closest to the road.This is because the particles formed through BHN are much smaller in size than those directly emitted with non-volatile cores around 15 nm in diameter.Thus, particles of BHN origins are subject to faster dry deposition removal, and contribute less to N at greater distances in the near-road environment.However, it should not be ignored in air quality modeling studies of mobile emissions, especially within the first 100 m of the roadways.
The scenario excluding particle coagulation alone results in the least impact on both N and geometric mean diameter of nucleation mode particles near the road.The results strongly agree with both timescale analysis (Zhang and Wexler, 2004) and previous CFD modeling studies (Wang and Zhang, 2012;Albriet et al., 2010;Gidhagen et al., 2004b).However, the coagulation process was suggested to be important under mild to weak atmospheric dilution conditions, such as street canyons (Gidhagen et al., 2004a) and road tunnels (Gidhagen et al., 2003).

Role of atmospheric boundary layer
Previous studies have shown that accurate CFD simulation of the ABL (including its wind profile and turbulence quantities) is essential for atmospheric dispersion of inert pollutants.For example, Gorle et al. (2009) investigated the effect of atmospheric TKE on the dispersion of particles of 1 µm in diameter, and concluded the impact was significant.In their study, however, aerosol dynamical processes were not con- sidered, nor were their interactions with the ABLT.Here, a sensitivity analysis on the ABL is performed to investigate the impact of the ABL on UFP formation and dispersion in the near-roadway environment.Specifically, the base case simulation is compared with a test simulation where only the wall-function modifications of Eqs. ( 7)-( 8) are not applied.
Figure 6a shows the model predicted UFP concentrations of the base case and the test simulations, along with the integration of SMPS data at two fixed locations for the 06:00-08:00 LT morning rush hours.The predicted concentrations of UFPs from the test simulation are lower by about 1 × 10 4 particles/cm 3 (or 18 % of the background corrected peak concentration) near the center of the highway compared to the base case.However, the concentration difference for CO 2 (as shown in Fig. 6b) between the two simulations is only slightly different (∼ 10 % of its background corrected peak value).The underestimated concentrations of both pollutants are due to unrealistic acceleration of the surface wind and changes in the TKE profile in the upstream region of the computational domain, as discussed in Blocken et al. (2007b).Thus, the concentration underestimation for CO 2 is the result of the overestimated dilution near the surface, where vehicular exhaust occurs.
In addition to the overestimated dilution effect on particle dispersion, the impact of the ABL on UFP number concentrations is enhanced by the reduced nucleation rate due to the underestimation of gaseous precursors (i.e., H 2 SO 4 and H 2 O) in the vehicle wake regions.The maximum nucleation rate from both simulations is around 2.9 × 10 16 particles m −3 s −1 , which are in the range of 1-6.3 × 10 16 particles/m 3 /s from single vehicle exhaust plume simulations (Wang and Zhang, 2012;Uhrner et al., 2007).Although the maximum nucleation rate is not sensitive to ABL profiles, the area-averaged nucleation rate of the cross section in the exhaust pipe plane behind the vehicle is underestimated by a factor of 5 due to the overestimated dilution behind vehicles in the sensitivity run.This comparison strongly suggests that the concentration of UFPs from mobile sources may be even more sensitive to the ABL conditions than inert gaseous species.It also implies that introducing ABL conditions to activity-based emission models (such as Nikolova et al., 2011b) may potentially improve their performance in estimating UFP traffic emissions.

Conclusions
In this study, an aerosol dynamics-CFD coupled model is applied to a single unified computational domain to investigate the dynamics and dispersion of UFPs from tailpipe exhaust to the near-road environment.The interactions among individual exhaust plumes are explicitly modeled within the tailpipe-to-ambient computational domain.The unique application of translational periodic boundary conditions effectively reduces the size of the computational domain and allows fast multiple-scenario simulations of size-resolved and chemical-component-resolved aerosol dynamics.This paper has demonstrated that, together with field measurements, the model is an effective tool which can be used to advance our knowledge on the formation and dispersion of UFPs in the near-road environment.This information is needed to help develop parameterizations of sub-grid processes ultimately to improve air quality model simulations over urban areas.
The model was successfully validated with FEVER field study measurements of both on-road and near-road TKE.The results indicate that the strength of turbulent mixing of pollutants due to VIT and the ABLT is properly captured by the model, leading to good agreement between modeled and measured concentrations for CO 2 and BC.For UFPs, the modeled PSDs demonstrated adequate agreement with measurements at two fixed locations near a major highway, under different traffic conditions.Sensitivity analysis indicated that the modeled N and PSD of UFPs are sensitive to H 2 SO 4 -H 2 O binary homogeneous nucleation, condensation/evaporation of SVOCs, and particle dry deposition.However, for such an unconfined near-road environment as in this study, coagulation appears to have a negligible effect on UFPs.Results also suggest that UFPs from mobile sources may be even more sensitive to ABL conditions than inert species because the average nucleation rate in vehicle wakes is very sensitive to the dilution of H 2 SO 4 .Therefore, introducing ABL conditions to activity-based emission models may potentially improve their performance in estimating UFP traffic emissions.
The next step of our study is to conduct multiple-scenario simulations and to provide parameterizations on mobile emissions to large-scale air quality models.Ultimately, a unified model of UFP dynamics and dispersion can help predict exposures in the vicinity of roadways to support risk assessments and health effect studies.Better treatment of these processes in air quality models is also expected to allow for more accurate predictions of the impacts specific vehicle emission control strategies may have on near-road air quality and subsequent health and environmental benefits.

Table 2 .
Fig. 1.Computational domain (a) and running vehicles and ground mesh (b).Purple mesh indicates velocity-inlet boundaries (left and top); Red mesh indicates pressure-outlet boundary (right); Black mesh indicates wall boundaries (bottom and cars); and cyan mesh indicates translational periodic boundaries (lateral).

Figure 1 .
Figure 1.Computational domain (a) and running vehicles and ground mesh (b).Purple mesh indicates velocity-inlet boundaries (left and top); red mesh indicates pressure-outlet boundary (right); black mesh indicates wall boundaries (bottom and cars); and cyan mesh indicates translational periodic boundaries (lateral).

Fig. 2 .
Fig. 2. Comparison of the on-road TKE from the passenger vehicle chasing experiments of the FEVER project (black line) and model simulations (red, blue and purple lines).The error bars represent the 25 th and 75 th percentiles of the measured on-road TKE.PC stands for passenger vehicle and SUV stands for sport utility vehicle.V is the average travelling speed (m/s) of PC in chasing experiments or the vehicle speed used in model simulation.

Figure 2 .
Figure 2. Comparison of the on-road TKE from the passenger vehicle chasing experiments of the FEVER project (black line) and model simulations (red, blue, and purple lines).The error bars represent the 25th and 75th percentiles of the measured on-road TKE.PC stands for passenger vehicle and SUV stands for sport utility vehicle.V is the average traveling speed (m s −1 ) of PC in chasing experiments or the vehicle speed used in model simulation.

Fig. 3 .
Fig. 3. Comparison of the TKE from the FEVER observations at a roadside tower and model simulations for morning rush hours: 05:00-06:00 (red) and 06:00-08:00 a.m.(blue).Measurement data are plotted in solid lines and model simulation results are plotted in dashed lines.The error bars on top of measurement data represent 25 th to 75 th percentiles.

Figure 3 .
Figure 3.Comparison of the TKE from the FEVER observations at a roadside tower and model simulations for morning rush hour: 05:00-06:00 (red) and 06:00-08:00 LT (blue) (measurements in image referred to in a.m.).Measurement data are plotted as solid lines and model simulation results are plotted as dashed lines.The error bars represent 25th to 75th percentiles.

Fig. 4 .
Fig. 4. Comparison of modelled and measured near-road concentrations of CO2 (ppmv) and BC (µg/m 3 ) on the downwind side of Hwy-400.Median concentrations from FEVER measurements are plotted in black solid lines, and modelled concentrations are plotted in red dashed lines.The grey areas represent measurements within 25 th and 75 th percentiles.

Figure 4 .
Figure 4. Comparison of modeled and measured near-road concentrations of CO 2 (parts per million volume, ppmv) and BC (µg m −3 ) on the downwind side of Hwy 400.Median concentrations from FEVER measurements are plotted in black solid lines, and modeled concentrations are plotted in red dashed lines.The gray areas represent measurements within 25th and 75th percentiles.

Fig. 6 .Figure 6 .
Fig. 6.Predicted UFP number (a) and CO2 (b) concentrations as a fu 4 the centre of Hwy-400. 5 Figure 6.Predicted UFP number (a) and CO 2 (b) concentrations as a function of distance to the center of Hwy 400.

Table 1 .
Median values of the measured meteorological data and model input.

Table 3 .
Background concentrations of particulate and gaseous species considered in the model.

2014 12638 L. Huang et al.: Aerosol-computational fluid dynamics modeling of ultrafine and black carbon particle emissionTable 4 .
Tailpipe gaseous and particulate species mass emission rates (g km −1 driven) for gasoline engines.

Table 5 .
Particle number-size distribution parameters assumed by the model.

Table 6 .
Number concentration and mean diameter of the nucleation mode particles predicted by the model under different scenarios.Normalized bias (in percentage) was calculated as (base case -scenario)/base case × 100 %.