urbanChemFoam 1.0: Large-Eddy Simulation of Non-Stationary Chemical Transport of Traffic Emissions in an Idealized Street Canyon

This paper describes a large-eddy simulation based chemical transport model, developed under the OpenFOAM framework, implemented to simulate dispersion and chemical transformation of nitrogen oxides from traffic sources in an idealized street canyon. The dynamics of the model, in terms of mean velocity and turbulent fluctuation, are evaluated using available stationary measurements. A transient model run using a photostationary reaction mechanism for nitrogen oxides and ozone subsequently follows, where non-stationary conditions for meteorology, background concentrations, and traffic 5 emissions are applied over a 24-hour period, using regional model data and measurements obtained for the City of Berlin in July, 2014. Diurnal variations of pollutant concentrations indicate dependence on emission levels, background concentrations, and solar state. Comparison of vertical and horizontal profiles with corresponding stationary model runs at select times show that, while there are only slight differences in velocity magnitude, visible changes in primary and secondary flow structures can be observed. In addition, temporal variations in diurnal profile and cumulative species concentration result in significant 10 deviations in computed pollutant concentrations between transient and stationary model runs.

speciated spatio-temporal redistribution of annual NO X emission levels provided by the Berlin City Senate. Urban background concentrations for NO, NO 2 , and O 3 are obtained from the Berlin Air Quality Measurement Network (Berliner Luftgüte-Messnetz, BLUME). In-canyon vertical and horizontal profiles will be compared against corresponding stationary model runs at specific times of the day, in addition to diurnal variations of species concentration at specific locations.

95
The prognostic equations employed in the present model will be formulated using a weakly compressible framework. This generalizes the theoretical foundation, and opens up possibilities of direct coupling with other compressible regional models such as WRF in the future. These equations form the dynamic core of the modelling framework, to which the conservation equations for mass, momentum, energy, and species are to be solved numerically to determine the spatial and temporal distribution and evolution of the prognostic variables, namely, wind velocity components (u i ), pressure (p), internal energy (e), and 100 mass mixing ratios for each gas-phase component species (y q ): where the notation D t (·) ≡ ∂ t (·)+∂ j u j (·) represents the substantial derivative, ρ is the density of the gas mixture, h ≡ e−p/ρ is the mixture enthalpy, and r is the net chemical production rate of each gas-phase species. The Coriolis coefficient is denoted by 2 sin (φ) ω j , as a function of the planetary angular velocity ω, and the latitude of the point of interest φ. The shear stress tensor τ ij in Equation (2) is represented by τ ij ≡ µ [(∂ i u j + ∂ j u i ) − (2/3) (∂ j u i ) δ ij ], where µ is the dynamic viscosity, and δ ij is the Kronecker delta. Enthalpy changes from chemical reactions are deemed negligible and thus are not considered in 110 Equation (3).
Due to the computational effort required to directly resolve all possible turbulent scales in any flow of practical interest, some form of statistical treatment must be adopted. In the LES model approach, this is realized by decomposing the instantaneous scalar variable ϕ into a resolved scale, ϕ , and a residual scale, ϕ ∆ , through spatial convolution with a low-pass filter kernel F (Leonard, 1974): The filter width ∆ is associated with the smallest turbulent scale that can be retained in the resolved scale. In practice, ∆ is derived from the local mesh dimensions, for example, the cube root of the cell volume, such that ∆ ≡ (δx 1 · δx 2 · δx 3 ) 1/3 .
Thus all turbulent scales residing within the resolved scale, that is, those with spatial scales larger than the attenuation threshold defined by F (∆), can be captured directly with the computational grid, while the residual turbulent scales, which are commonly 120 referred to as the subgrid scale (SGS), are modelled.
As with Reynolds averaging for time-averaged scalars, additional terms will result from the filtering operation of equations (1 -4), deriving from the advective term in the substantial derivative D t ϕ : which implies where σ ϕ ij is known as the SGS term, and in the case of the momentum equation (2) this is referred to as the SGS stress tensor, and is denoted simply as σ ij . Further, in a compressible flow framework, where density also varies, additional SGS terms from density filtering will be introduced in Equation (1). This can be avoided, however, by performing the filtering operation on a density averaged scalar, where ϕ = ρϕ / ϕ , through a process known as Favre filtering (Favre, 1965). The corresponding 130 filtered transport equations for equations (1 -4) then become: where the notation for Favre-filtered substatial derivative is written as D t (·) ≡ ∂ t (·) + ∂ j u j (·). In addition, Vreman et al. (1995) show that nonlinearities arising from filtering of diffusion terms can be neglected for weakly compressible flows based on a priori direct numerical simulation (DNS) data. The expanded form for the corresponding SGS terms for equations (7 -140 10) are presented below: 145 Detailed discussions pertaining to the modelling and parameterization of the remaining SGS terms described in equations (12 -14) can be found in Martín et al. (2000).
Numerous approaches are available for modelling the SGS terms (Garnier et al., 2009). A common method, first proposed by Smagorinsky (1963), is to model the SGS stresses by invoking the gradient transport theorem: The closure of subgrid scale viscosity µ ∆ is, in turn, achieved by using representative turbulent scales. Based on the one equation model of Deardorff (1980) for incompressible flows, Yoshizawa (1986) proposed a similar closure scheme for weakly compressible flows using the subgrid scale turbulent kinetic energy (TKE) 155 and the subgrid scale viscosity µ ∆ can then be evaluated as follows: in which both C k and C are dimensionless single valued constants set to 0.094 and 1.05 respectively.

Thermophysical and Chemistry Models
In the current study, the photostationary NO -NO 2 -O 3 mechanism (Leighton, 1961) will be considered. The mechanism 160 comprises the photodissociation of NO 2 into NO, the titration of O 3 by NO to form NO 2 , and an intermediate third-body reaction for in situ ozone production through oxygen radicals originated from NO 2 photolysis: Reaction (R1) describes the photodisocciation of NO 2 , whose rate constant k is calculated as a function of the solar zenith angle (χ): Meanwhile, the rate constants for reactions (R2) and (R3) are based on the Arrhenius rate law:

170
where A and T A are the corresponding pre-exponential factor and activation temperature for each reaction. All parameters for the aforementioned reactions are tabulated with reference in Table 1. net production rate r q in equation (14) is determined by summing up the generation and consumption rates of the corresponding species across all relevant reactions. As non-limiting examples, the net production rates for NO, NO 2 , and O 3 for the current chemical mechanism are calculated as follows: Transport and thermophysical properties for these species are available from McBride et al. (1993). Temperature dependence of specific heat for species q at constant pressure (c q p ) is presented as a polynomial function: where the polynomial coefficients a i are obtained using least-square fit, and R q is the specific gas constant for the species of 185 interest. Meanwhile, dynamic viscosity (µ) is represented using Sutherland's law: where A S (in Pa·s·K −1 ) and T S (in K) are species-specific single-valued coefficients. The thermodynamic state of the mixture is determined by the ideal gas relation:

190
where R is the specific gas constant for the gas mixture.

Model Implementation
The open source finite volume computational continuum mechanics framework OpenFOAM (Weller et al., 1998) has been used as the modelling foundation for this work, due to its ability to handle compressible reactive flows in an unstructured grid arrangement. More specifically, the current flow solver and related data processing utilities are derived from OpenFOAM version

Domain Discretization and Decomposition
The solution domain is disretized using the OpenFOAM blockMesh utility, which generates a conformal hexahedral grid. A 220 core region is defined by a 9 m lengthwise section about the domain origin, covering the entire span and height of the street canyon, where grid cells are uniform and have near unity aspect ratio. Grid cell grading is introduced beyond the core region, so that cell size gradually expand at a rate according to a user-specified final expansion ratio relative to the core cell size, which dictates the size of the furthermost cells. For the present study, this expansion ratio is set to 1:2 in the spanwise direction, 1:3.2 and 1:4 in the lengthwise and heightwise directions respectively. Two mesh configurations of core cell sizes of 1 m (coarse 225 mesh) and 0.5 m (fine mesh) have been deployed in the evaluation study, corresponding to a total cell count of 43 500 and 344 000 cells. Following the results from the stationary evaluation (Section 3), the transient study is conducted only at the fine mesh level. Figure 2 shows the mesh cell distribution for the fine mesh arrangement (core cell size 0.5 m), showing core and graded regions.
The numerical schemes used in the discretization of the prognostic equations are second-order accurate in time and space. A 230 collocated grid approach is used, in which pressure-momentum coupling is achieved using a fourth-order interpolation scheme (Rhie and Chow, 1983). For each time step, each discretized prognostic equation is solved interatively using the biconjugate gradient method (van der Vorst, 1992) with a diagonal incomplete lower-upper matrix factorization (ILU) preconditioner, until the global normalized residual falls below a specified convergence criterion, in this case 10 −5 . Mass conservation is enforced at every iteration by solving an iterative pressure correction equation until the aforementioned residual threshold is 235 reached. Typically, the pressure correction equation converges in two to three iterations, and about five to eight iterations for the prognostic equations at each time step for the sizes of grid cells and time step considered in this study. Grid decomposition for parallel computing is performed in three dimensions using a dual recursive bipartitioning algorithm (Chevalier and Pellegrini, 2008).

240
For the purpose of evaluating the LES street canyon model, a stationary run is performed for which LDA flow measurement data are available from Li et al. (2008) for a scale model infinite canyon arrangement, also having an height-to-width aspect ratio of two, conducted in a water tunnel. The horizontal velocity of this simulation is prescribed at the freestream boundary, denoted as U ∞ . An initially quiescent velocity field is imposed in the canyon region, while a linear spanwise flow profile -  to adequately capture the profiles. There is a discrepancy at x/W = 0.25, but the salient features are still well represented.
In comparison, the coarse mesh results underpredict the velocity magnitude. In light of the fine mesh results, this observation indicates an issue with grid convergence. In addition, substantially better agreement in turbulent fluctuation (u and w , lower left and right panels, respectively) can be found with the fine mesh results. There are still minor differences between the fine grid LES results and the LDA measurements, but can be attributed to the use of only resolved turbulent fluctuations in deriving the time-averaged results. In addition, as Zhong et al. (2015) also indicated, the circumstances under which the experiment (Li et al., 2008) is conducted, methods of vortex generation for instance, are not modelled, which could attribute, at least in part, 270 to the observed differences.
The mean flow field in the canyon along the lateral plane of symmetry has been obtained by averaging the instantaneous velocity field at interval of 300 s for over the 3600 s sampling period for the fine mesh configuration. It is presented in Figure 4 using a line integral convolution representation (LIC; Cabral and Leedom, 1993). It depicts a typical skimming flow regime occurring in idealized urban canyon arrays (Oke, 1988), where a free shear layer develops over the canyon at z/H 275 = 1 resulting from turbulent interaction between the building roof and the freestream flow upstream. Two vertically stacked primary recirculation zones exist inside the canyon, in line with existing works on deep street canyons (Baik and Kim, 1998;Li et al., 2008). However, additional, minor vortices are also observed at the corners at the street level, which may lead to additional isolation of street level emissions, for instance, from traffic sources.
Although, in the skimming flow regime, there is little mass transfer between the bulk flow and the canyon, as indicated by To further evaluate the resolution of turbulent scales under each mesh configuration, the fraction of TKE as its simulated total -that is, the sum of resolved and SGS TKE -is calculated along the three vertical stations, as presented in Figure 5. The free shear layer on top of the canyon roof region is accompanied by a sudden increase in turbulence kinetic energy upstream in the canyon core, which dissipates downstream (Driver and Seegmiller, 1985). However, a decrease in resolved scale TKE 290 fraction can be observed in this region. In particular, the attenuation in turbulence resolution is significantly higher in the coarse mesh than in the fine mesh configuration, with a minimum of 36.0% and 73.6% for the respective mesh density. This indicates a diminished capability of the coarse mesh to adequately capture turbulent structures explicitly, which contributes to the discrepancies in the vertical profiles observed in Figure 3. Nevertheless, the resolved TKE fraction in the canyon core and the bulk flow regions (that is, outside of the free shear flow region) remains very high for both grid configurations. Based 295 on a 95% confidence interval for regions away from walls and the free shear layer, the resolved fraction is calculated to be 86.0% -87.5% for the coarse mesh, and 92.6% -93.3% for the fine mesh. This suggests that the turbulent resolution in the free shear region could be summarily resolved by local mesh refinement. However, based on Figure 5, the fine mesh configuration provides an adequate resolution of turbulent flow structures. Hence subsequent discussions will be made exclusively with the 0.5 m mesh configuration. In addition, there is a reciprocal relationship between the concentrations of NO and NO 2 , as well as O 3 . This is expected as the O 3 originates from the freestream and is titrated through reaction (R3) by the NO originating from street level traffic.
However, the concentration of in-canyon O 3 is not substantially different from freestream levels, at least in comparision with 365 the NO and NO 2 levels. This is, in part, due to the in situ production of O 3 through photodissocation in reaction (R1). The exception is 20:00, where the hourly sampling and averaging takes place with the absence of the photolysis reaction (R1), as the solar zenith angle (χ) reaches 90 • at 19:25, leaving the titration reaction (R3) the only active reaction in the mechanism.
Combined with a low traffic emission, as indicated in Figure 6(D), this leads to a moderate decrease in O 3 while depleting the NO in the canyon. The stratification of NO X at each horizontal level follow the direction of the vortex at the respective elevation, with higher concentration levels found towards the right wall in the lower vortex (Figures 9 and 10), and the left wall in the upper vortex ( Figure 11). However, at the 2 meter level, the peaks of NO X concentration resulting from diffusion of street level emissions, as depicted in Figure 1(C), is not evidently shown. This is indicative of the dominance of spanwise advection near the street 380 surface in the lower vortex, which results in the aforementioned high concentration region towards the right side of the lower vortex. By comparison, the concentration of O 3 remains more or less stable at all horizontal levels.
Finally, Figure 12 shows the diurnal cycles in the concentrations of NO, NO 2 , and O 3 in the street canyon from the transient model run, at nine different locations corresponding the intersections between the vertical (Figure 8) and horizontal profiles (Figures 9 to 11). Both the NO and O 3 profiles exhibit a typical dirunal cycle, which show strong dependence on solar activity, 385 since only the titration reaction (R3) remains active before sunrise and after sunset (that is, |χ| ≥ 90 • ), as previously pointed out. Thus, in conjunction of emission rates from traffic NO X during these times, a drastic reduction in NO can be seen. Further, the O 3 concentration at all monitoring points follow the same tendency as the background O 3 shown in Figure 6(C), except at the time immediately following the cessation of solar activity (at 19:25), where an acute reduction is observed due to NO X still presently available, mainly from emission sources, for the titration reaction (R3).

390
While the trend of NO 2 concentration in the afternoon and early evening (that is, between 15:00 and 21:00) is in line with the traffic emissions, the overall diurnal variation does not reflect the morning traffic peak, but instead only gradually increases following sunrise at 02:35. As the morning period is accompanied with a comparatively low background O 3 concentration, the photodissociation of NO 2 becomes the dominant reaction in the mechanism, leading to consumption of the emitted NO 2 from the traffic sources. As the O 3 concentration increases in the afternoon, the rate of titration increases correspondingly. Therefore 395 the afternoon NO 2 peak coincides with a decrease in NO concentration.

Concluding Remarks
A weakly compressible, reactive finite volume LES solver, based on the OpenFOAM computational continuum mechanics framework (Weller et al., 1998), has been developed for modelling chemical transport of gas-phase pollutant species in an with time. A deep street canyon with a height-to-width aspect ratio of two (Zhong et al., 2015), has been used throughout this study. The solver has been evaluated using stationary LDA measurements from Li et al. (2008) The results of this study demonstrate the ability for this solver in performing urban-scale chemistry transport modelling at both a stationary and non-stationary capacity in a grid-independent manner. As such, it can be readily adapted to model a larger domain with more realistic geometry, such as a city block. For models covering such domain regions, however, input data such 415 as meteorology and emission profiles need be prescribed at higher spatial and temporal resolutions, which can be accomplished by coupling the urban scale model with outputs from regional models, WRF-Chem for example. Particular emphasis can be placed on the the role of absorptive, reflective, and refractive mechanisms of urban structures in the redistribution of solar radiant energy, and their impact on the overall advective and convective thermal stratification in their immediate vicinity.
Further, the method of temporal-spatial distribution and NO X emissions from annual aggregate, and the subsequent partition 420 into NO and NO 2 contributions, as outlined in Section 4.1, can be extended to cover a road network, where emissions can be mapped as discretized surface fluxes in a similar manner as the simpleEmission boundary condition described in Section 2.3. In addition to NO X , the emissions of PM, carbon monoxide (CO), and volatile organic compounds (VOCs) can be processed in a similar way. While the photostationary NO -NO 2 -O 3 mechanism (Leighton, 1961)  Mechanisms such as CBM-IV (Gery et al., 1989), involving 32 species and 81 reactions, or, to a lesser degree, the SMOG mechanism (Damian et al., 2002), made up of 13 species and 12 reactions, can be considered in subsequent investigations.  Lupaşcu, and J. Quedenau for their valuable inputs and support in the conception and preparation of this work.