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 emissions are applied over a 24 h 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 deviations in computed pollutant concentrations between transient and stationary model runs.

1 Introduction

The study of dispersion and chemical transformation of pollutants at urban scales (less than 200 m; Britter and Hanna, 2003) is not only of immediate scientific and engineering interest, but it is also an indispensable exercise for socioeconomists and policy makers seeking to quantify their decisions. At this scale, the local population is in close proximity with the pollutant sources, such as nitrogen oxides (NO_x) from traffic sources. Based on recent reports from the European Environmental Agency, nitric oxide (NO) and nitrogen dioxide (NO₂) derived from road traffic constitute close to 40 % of total NO_x emissions in the EU-28 nations (EEA, 2019a). In addition to being an irritant and a potential carcinogen, NO₂ acts as precursor to a number of harmful substances, such as nitric acid (HNO₃), ground-level ozone (O₃), and particulate matter (PM). Further, it is estimated that the associated aggregated global economic impact, including diminished productivity, health care expenditure, and reduced agricultural crop yields, could increase from USD 330 billion in 2015 to up to USD 3300 billion in 2060 (OECD, 2016).

In 2017, about 10 % of all air quality monitoring stations among the EU-28 have reported NO_x exceedances of an annual mean limit of 40 µg m⁻³, with the majority observed in urban traffic stations (EEA, 2019b). In Germany, NO_x exceedances were reported in over 40 % of urban traffic stations in 2018 (Minkos et al., 2020). Although these figures already represent a continually downward trend from 2014 levels (EEA, 2019a), additional measures are still required to meet the European air quality targets (EEA, 2019b). In conjunction with rigorous air quality monitoring programs and regional-scale modeling (Kuik et al., 2016; von Schneidemesser et al., 2017), urban-scale modeling can offer additional insights into localized formation and dispersion of pollutants and their impact on air quality at regional and urban scales, thus providing actors involved with environmental policy crucial information for identifying effective mitigation strategies.

As the basic geometrical unit in a typical urban environment, a street canyon consists of two rows of buildings separated by a road section. The aspect ratio between the building height (H) and the road section width (W) of the canyon is an important geometrical characteristic, as it has significant influence on pollutant transport (Liu et al., 2005). The flow in the freestream region above the canyon induces recirculation zones within the canyon, as demonstrated in numerical studies of infinite canyon arrangements (Baik and Kim, 1998; Li et al., 2008), which restricts exchange of air mass to the shear layer between the canyon and freestream regions (Driver and Seegmiller, 1985). The geometrical arrangement of these flow structures could be altered significantly, however, when various parts of the street canyon (that is, the building walls and street surface) are maintained at different temperatures, through solar heating, for instance, which, in combination with the prevailing wind, could increase stratification of pollutant concentrations in the canyon (Kim and Baik, 2001; Madalozzo et al., 2014; Xie et al., 2007). In particular, the formation of convective cells under a full or mixed convective regime introduces mean flow along the street canyon, whose magnitude depends on the canyon geometry and Rayleigh number, which cannot be assumed as negligible (Davies-Jones, 1970).

The Reynolds-averaged Navier–Stokes (RANS) equations are a popular numerical framework for modeling turbulent dynamics due to its simplicity and correspondingly low computational cost. While the RANS approach has been shown to adequately represent mean flow behaviors in some cases (Baik and Kim, 1998; Kwak and Baik, 2014), in the presence of large-scale unsteady mixing, such as vortex shedding from building surfaces, mean dynamics is not always sufficient, as these motions are not universal according to the Kolmogorov hypothesis (Pope, 2000), seen in part by large discrepancies observed in simulated pollutant concentrations for dispersion scenarios in comparison with measurement data (Tominaga and Stathopoulos, 2009). On the other hand, using large-eddy simulation (LES), large-scale turbulent motions can be resolved explicitly, leaving closure modeling to the smaller scales, which exhibit universal character to some extent. Furthermore, unlike RANS, where the prognostic equations are already averaged in time, the resolved scale statistical quantities in LES are evaluated after the computation has taken place. Hence LES can arguably provide a much more accurate and complete representation of the evolution of turbulent flow structures, as well as chemical transport of pollutant species.

Further, a number of studies have considered the coupling of dynamics with reaction kinetics relevant to urban traffic emissions. One such simple mechanism is the photostationary NO–NO₂–O₃ mechanism (Leighton, 1961). In a street canyon, photochemical steady state is well maintained in the individual vortices and ground-level corners where pavements are located, while deviations are observed close to the traffic emission sources and in the roof-level shear layer (Baker et al., 2004). Deviations from the photo steady state indicate heterogeneity due to imperfect mixing of NO_x (Zhong et al., 2015) or presence of NO_x oxidizing agents other than O₃, such as hydroperoxy (HO₂) or organic peroxy (R-O₂) radicals, formed mainly through oxidation of carbon monoxide (CO) and hydrocarbons (R-H) initiated by hydroxyl (OH) radicals (Carpenter et al., 1998). Nevertheless, the photostationary mechanism has been shown to be adequate in modeling the concentrations of NO_x and O₃. In comparison with more detailed mechanisms, such as the Carbon Bond Mechanism (CBM-IV; Gery et al., 1989), Garmory et al. (2009) reported a difference of about 1 ppb for concentrations of NO, NO₂, and O₃. Similarly, similar levels of discrepancy are also observed in Bright et al. (2013) with the Reduced Chemical Scheme (RCS).

The present study seeks to extend the current corpus on idealized street canyon modeling in two aspects. First, previous studies have been conducted in stationary meteorological and pollutant conditions, inasmuch transient studies are carried out under statistically steady wind and constant emission rates (Kwak and Baik, 2014; Kristóf et al., 2020). Thus, the possibility of introducing diurnal variations on all forcings and their effects on the formation and spatial distribution of pollutant species can be contemplated in this investigation. This could, for instance, help conceptualize more comprehensive parametrization schemes, in which effects of transient conditions could be considered. Second, using the idealized street canyon as a starting point, the chemical transport model can be further conceived and developed to simulate more realistic urban environments using detailed chemical kinetics, coupled with meteorological conditions and background concentrations extracted from measurements or regional models.

Figure 1 Schematic of the solution domain representing the infinite street canyon geometry from Zhong et al. (2015): (a) basic dimensions and locations of cyclic and freestream boundaries, (b) isometric view of the domain showing lateral dimension and corresponding cyclic boundaries, and (c) top view of ground surface indicating size and position of simulated traffic lanes.

2 Model description

As illustrated in Fig. 1, the solution domain representing an infinite deep street canyon will be used for this study throughout. To facilitate discussion, x will be referred to as the spanwise direction, y the lengthwise direction, and z the vertical direction. The origin of the domain is located on the street level, at the intersection between the spanwise (x−z) and lengthwise (y−z) planes of symmetry. The span of the street section (W) is 18 m, and the building height (H) is 36 m, resulting in an height-to-span aspect ratio of two. The freestream boundary is located at 76 m above the canyon roof, or 112 m above street level, from which meteorological conditions are specified. The canyon roof extends 9 m on each side of the spanwise direction, as shown in Fig. 1a. The street section length is set to 40 m, as shown in Fig. 1b, centered about the domain origin. Traffic emissions are introduced as two 1.5 m wide sections in equidistant arrangement, as shown in Fig. 1c, and are assumed to be continuous and spatially uniform along both sections. Cyclic conditions are imposed on the spanwise and lengthwise boundary pairs to ascertain periodicity in the corresponding directions, as well as to allow turbulent flow structures to freely develop across the solution domain. To prevent accumulation of pollutants through the spanwise cyclic boundaries while maintaining background concentrations, chemical species will be introduced from the upstream spanwise boundary at specific, spatially uniform mixing ratios, while a zero gradient boundary condition will be imposed on the downstream spanwise boundary, allowing excess pollutants to leave the solution domain.

Initially, the dynamics of this model framework will be evaluated using the stationary laser Doppler anemometry (LDA) measurements of Li et al. (2008). This is followed by a transient simulation run with reaction kinetics activated. Boundary conditions representative of a working day in Berlin in July 2014 will be used. Diurnal variations for ambient wind and thermodynamics are extracted on an hourly basis from the regional model results of Kuik et al. (2018) using the Weather Research and Forecasting (WRF) Model (Skamarock et al., 2008). Meanwhile, hourly traffic emissions are derived from speciated spatiotemporal redistribution of annual NO_x emission levels provided by the Berlin City Senate. Urban background concentrations for NO, NO₂, and O₃ 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.

2.1 Prognostic equations

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 modeling 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 mass mixing ratios for each gas-phase component species (y_q):

$$\begin{array}{}\text{(1)}& {\displaystyle }{D}_t\mathit{\rho }=\mathrm{0},\text{(2)}& {\displaystyle }{D}_t\left(\mathit{\rho }{u}_i\right)=-{\partial }_{i}p+{\partial }_j{\mathit{\tau }}_{ij}-{\mathit{ϵ}}_{ijk}\left[\mathrm{2}\sin \right(\mathit{\varphi }\left){\mathit{\omega }}_j\right]\left(\mathit{\rho }{u}_k\right),\text{(3)}& {\displaystyle }{D}_t\left[\mathit{\rho }h+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathit{\rho }{u}_i{u}_i\right)\right]-{\partial }_{t}p={\partial }_j[{\partial }_{j}e+{\mathit{\tau }}_{ij}{u}_i],\text{(4)}& {\displaystyle }{D}_t\left(\mathit{\rho }{y}_q\right)={\partial }_j\left(D_{pq}{\partial }_j{y}_q\right)+\mathit{\rho }{r}_q,\end{array}$$

where the notation $D_t(\cdot )\equiv {\partial }_t(\cdot )+{\partial }_j{u}_j(\cdot )$ represents the substantial derivative, ρ is the density of the gas mixture, $h\equiv e-p/\mathit{\rho }$ is the mixture enthalpy, and r is the net chemical production rate of each gas-phase species. The Coriolis coefficient is denoted by 2sin (ϕ)ω_j, as a function of the planetary angular velocity ω, and the latitude of the point of interest ϕ. The shear stress tensor τ_ij in Eq. (2) is represented by ${\mathit{\tau }}_{ij}\equiv \mathit{\mu }[{\partial }_i{u}_j+{\partial }_j{u}_i-(\mathrm{2}/\mathrm{3}\left)\right({\partial }_i{u}_i{\mathit{\delta }}_{ij}\left)\right]$, 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 Eq. (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):

$$\begin{array}{}\text{(5)}& \langle \mathit{\phi }\rangle =\underset{\forall }{\int }\mathit{\phi }\left(\mathit{\xi }\right)F(x-\mathit{\xi }|\mathrm{\Delta })\mathrm{d}\mathit{\xi }.\end{array}$$

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 $\mathrm{\Delta }\equiv (\mathit{\delta }{x}_{\mathrm{1}}\cdot \mathit{\delta }{x}_{\mathrm{2}}\cdot \mathit{\delta }{x}_{\mathrm{3}}{)}^{\mathrm{1}/\mathrm{3}}$. Thus all turbulent scales residing within the resolved scale, i.e., 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 referred to as the subgrid scale (SGS), are modeled.

As with Reynolds averaging for time-averaged scalars, additional terms will result from the filtering operation of Eqs. (1)–(4), deriving from the advective term in the substantial derivative D_t〈φ〉:

$$\begin{array}{rl}{\partial }_j\left(u_j\mathit{\phi }\right)={\partial }_j[& \langle u_j\rangle \langle \mathit{\phi }\rangle +\\ & \underset{{\mathit{\sigma }}_{ij}^{\mathit{\phi }}}{\underbrace{\langle u_j\rangle {\mathit{\phi }}_{\mathrm{\Delta }}+(u_j{)}_{\mathrm{\Delta }}\langle \mathit{\phi }\rangle +(u_j{)}_{\mathrm{\Delta }}(\mathit{\phi }{)}_{\mathrm{\Delta }}}}],\end{array}$$

which implies

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{ij}^{\mathit{\phi }}=\langle u_j\mathit{\phi }\rangle -\langle u_j\rangle \langle \mathit{\phi }\rangle ,\end{array}$$

where ${\mathit{\sigma }}_{ij}^{\mathit{\phi }}$ is known as the SGS term, and in the case of the momentum Eq. (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 Eq. (1). This can be avoided, however, by performing the filtering operation on a density averaged scalar, where $\stackrel{\mathrm{̃}}{\mathit{\phi }}\equiv \langle \mathit{\rho }\mathit{\phi }\rangle /\langle \mathit{\rho }\rangle$ through a process known as Favre filtering (Favre, 1965). The corresponding filtered transport equations for Eqs. (1)–(4) then become

$$\begin{array}{}\text{(7)}& {\displaystyle }{\stackrel{\mathrm{̃}}{D}}_t\langle \mathit{\rho }\rangle =\mathrm{0},\text{(8)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{D}}_t(\langle \mathit{\rho }\rangle {\stackrel{\mathrm{̃}}{u}}_i)& =-{\partial }_i\langle p\rangle +{\partial }_j({\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{ij}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij})\\ & -{\mathit{ϵ}}_{ijk}\left[\mathrm{2}\sin \right(\mathit{\varphi }\left){\mathit{\omega }}_j\right](\langle \mathit{\rho }\rangle {\stackrel{\mathrm{̃}}{u}}_k),\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{D}}_t[\langle \mathit{\rho }\rangle \stackrel{\mathrm{̃}}{h}& +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\langle \mathit{\rho }\rangle {\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_i)]-{\partial }_t\langle p\rangle \\ & ={\partial }_j\left[\right({\partial }_j\stackrel{\mathrm{̃}}{e}+{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{ij}{\stackrel{\mathrm{̃}}{u}}_i)-({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^h+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^{u_i}\left)\right],\end{array}\text{(10)}& {\displaystyle }{\stackrel{\mathrm{̃}}{D}}_t(\langle \mathit{\rho }\rangle {\stackrel{\mathrm{̃}}{y}}_q)={\partial }_j(D_{pq}{\partial }_j{\stackrel{\mathrm{̃}}{y}}_q-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^{y_q})+\langle \mathit{\rho }\rangle {\stackrel{\mathrm{̃}}{r}}_q,\end{array}$$

where the notation for Favre-filtered substatial derivative is written as ${\stackrel{\mathrm{̃}}{D}}_t(\cdot )\equiv {\partial }_t(\cdot )+{\partial }_j\left[{\stackrel{\mathrm{̃}}{u}}_j\right(\cdot \left)\right]$. 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 Eqs. (7)–(10) are presented below:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}=\langle \mathit{\rho }\rangle (\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j),\text{(12)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^h=\langle \mathit{\rho }\rangle (\stackrel{\mathrm{̃}}{u_{j}e}-{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{e})-(\stackrel{\mathrm{̃}}{\langle p\rangle {\mathit{\tau }}_{ij}}-\langle p\rangle {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{ij}){\mathit{\delta }}_{ij},\text{(13)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^{u_i}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle \mathit{\rho }\rangle (\stackrel{\mathrm{̃}}{u_i{u}_i{u}_j}-\stackrel{\mathrm{̃}}{u_i{u}_i}{\stackrel{\mathrm{̃}}{u}}_j),\text{(14)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}^{y_q}=\langle \mathit{\rho }\rangle (\stackrel{\mathrm{̃}}{u_j{y}_q}-{\stackrel{\mathrm{̃}}{u}}_j{\stackrel{\mathrm{̃}}{y}}_q).\end{array}$$

Detailed discussions pertaining to the modeling and parameterization of the remaining SGS terms described in Eqs. (12)–(14) can be found in Martín et al. (2000).

Numerous approaches are available for modeling 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:

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}\approx -\mathrm{2}{\mathit{\mu }}_{\mathrm{\Delta }}\left({\stackrel{\mathrm{̃}}{S}}_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\stackrel{\mathrm{̃}}{S}}_{ij}{\mathit{\delta }}_{ij}\right),\end{array}$$

where ${\stackrel{\mathrm{̃}}{S}}_{ij}\equiv (\mathrm{1}/\mathrm{2})({\partial }_i{\stackrel{\mathrm{̃}}{u}}_j+{\partial }_j{\stackrel{\mathrm{̃}}{u}}_i)$. 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 $k_{\mathrm{\Delta }}\equiv (\mathrm{1}/\mathrm{2})(\stackrel{\mathrm{̃}}{u_i{u}_i}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_i)$:

$$\begin{array}{}\text{(16)}& {\stackrel{\mathrm{̃}}{D}}_t(\langle \mathit{\rho }\rangle k_{\mathrm{\Delta }})={\partial }_j(\mathit{\mu }+{\mathit{\mu }}_{\mathrm{\Delta }}){\partial }_j{k}_{\mathrm{\Delta }}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{ij}{\partial }_j{\stackrel{\mathrm{̃}}{u}}_i-{\displaystyle \frac{C_{\mathit{ϵ}}\langle \mathit{\rho }\rangle k_{\mathrm{\Delta }}^{\mathrm{3}/\mathrm{2}}}{\mathrm{\Delta }}},\end{array}$$

and the subgrid scale viscosity μ_Δ can then be evaluated as follows:

$$\begin{array}{}\text{(17)}& {\mathit{\mu }}_{\mathrm{\Delta }}=C_k\langle \mathit{\rho }\rangle \mathrm{\Delta }{k}_{\mathrm{\Delta }}^{\mathrm{1}/\mathrm{2}},\end{array}$$

in which both C_k and C_ϵ are dimensionless single-valued constants set to 0.094 and 1.05, respectively.

2.2 Thermophysical and chemistry models

In the current study, the photostationary NO–NO₂–O₃ mechanism (Leighton, 1961) will be considered. The mechanism comprises the photodissociation of NO₂ into NO, the titration of O₃ by NO to form NO₂, and an intermediate third-body reaction for in situ ozone production through oxygen radicals originated from NO₂ photolysis:

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\text{NO}}_{\text{2}}+h\mathit{\nu }\to \text{NO}+{\text{O(}}^{\text{3}}\text{P)}\text{;}\mathit{\lambda }\le \mathrm{420}\text{nm},\text{(R2)}& {\displaystyle }{\text{O(}}^{\text{3}}\text{P)}+{\text{O}}_{\text{2}}+\text{M}\to {\text{O}}_{\text{3}}+\text{M},\text{(R3)}& {\displaystyle }\text{NO}+{\text{O}}_{\text{3}}\to {\text{NO}}_{\text{2}}+{\text{O}}_{\text{2}}\text{.}\end{array}$$

Reaction (R1) describes the photodissociation of NO₂, whose rate constant k is calculated as a function of the solar zenith angle (χ):

$$\begin{array}{}\text{(18)}& j_{\mathrm{R}\mathrm{1}}=\left\{\begin{array}{l}{J}_{\mathrm{0}}\cos {}^m\left(\mathit{\chi }\right)\exp \left[{\displaystyle \frac{-n}{\cos \left(\mathit{\chi }\right)}}\right]\left|\mathit{\chi }\right|<\mathrm{90}{}^{\circ }\\ \mathrm{0},\text{otherwise.}\end{array}\right.\end{array}$$

Meanwhile, the rate constants for Reactions (R2) and (R3) are based on the Arrhenius rate law:

$$\begin{array}{}\text{(19)}& k_{\text{R}i}=A_i\exp [-(T_{\mathrm{A}}{)}_{\text{R}i}/T],\end{array}$$

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.

Saunders et al. (2003)Huie et al. (1972)Burkholder et al. (2015)

Table 1 Reaction rate parameters for the photostationary NO-NO₂-O₃ mechanism. Units are kmol, m³, s, and K.

In the current study, air is modeled as a mixture of its primary constituents, namely N₂ and O₂, but the effects of H₂O and CO₂ are not considered at this stage. Together with the additional gaseous species introduced in Reactions (R1) to (R3), i.e., NO, NO₂, O(³P), and O₃, Eq. (14) consists of six species transport equations that need be solved. The individual species net production rate r_q in Eq. (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₂, and O₃ for the current chemical mechanism are calculated as follows:

$$\begin{array}{}\text{(20)}& {\displaystyle }{r}_{\text{NO}}=j_{\mathrm{R}\mathrm{1}}\left[{\text{NO}}_{\text{2}}\right]-k_{\mathrm{R}\mathrm{3}}\left[\text{NO}\right]\left[{\text{O}}_{\text{3}}\right],\text{(21)}& {\displaystyle }{r}_{{\text{NO}}_{\text{2}}}=-j_{\mathrm{R}\mathrm{1}}\left[{\text{NO}}_{\text{2}}\right]+k_{\mathrm{R}\mathrm{3}}\left[\text{NO}\right]\left[{\text{O}}_{\text{3}}\right],\text{(22)}& {\displaystyle }{r}_{{\text{O}}_{\text{3}}}=k_{\mathrm{R}\mathrm{2}}\left[{\text{O(}}^{\text{3}}\text{P)}\right]\left[{\text{O}}_{\text{2}}\right]\left[M\right]-k_{\mathrm{R}\mathrm{3}}\left[\text{NO}\right]\left[{\text{O}}_{\text{3}}\right].\end{array}$$

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_p^q$) is presented as a polynomial function:

$$\begin{array}{}\text{(23)}& c_p^q=R_q\sum _{i=\mathrm{0}}^{\mathrm{4}}{a}_i{T}^i,\end{array}$$

where the polynomial coefficients a_i are obtained using a least-squares fit and R_q is the specific gas constant for the species of interest. Meanwhile, dynamic viscosity (μ) is represented using Sutherland's law:

$$\begin{array}{}\text{(24)}& \mathit{\mu }={\displaystyle \frac{A_{\mathrm{S}}{T}^{\mathrm{1}/\mathrm{2}}}{\mathrm{1}+{\displaystyle \frac{T_{\mathrm{S}}}{T}}}},\end{array}$$

where A_S (in Pa s K⁻¹) and T_S (in K) are species-specific single-valued coefficients. The thermodynamic state of the mixture is determined by the ideal gas relation:

$$\begin{array}{}\text{(25)}& p=\mathit{\rho }RT,\end{array}$$

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

2.3 Model implementation

The open-source finite-volume computational continuum mechanics framework OpenFOAM (Weller et al., 1998) has been used as the modeling 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 7, maintained by CFD Direct Ltd. All pertinent modifications and augmentations at the source level are performed in the C++ programming language. The OpenFOAM framework contains core classes for LES and gas-phase thermophysics, chemistry, and reaction kinetics. To adopt OpenFOAM for urban-scale chemical transport modeling, a custom solver, urbanChemFoam, derived from the existing rhoReactingFoam, has been developed to enable geodetic and solar calculations using a Julian-Day-based date and time module. The solar state, represented by the solar zenith angle (χ) can either be calculated at each time step using the algorithm of Meeus (1998), via a user-defined reference time and location, or set to a fixed value by the user throughout the model run.

Additional components have also been developed to operate in conjunction with urbanChemFoam. As existing chemical reactions in OpenFOAM are functions only of the thermodynamic state, a new reaction class, called irreversiblephotolysisReaction, has been developed to model photodissociation as a function of the cosine of the solar zenith angle. Further, surface emissions, for instance from traffic sources, are provided in units of mass flux (i.e., kg s⁻¹). These are typically introduced in chemical transport models as volumetric source terms in the interior cell nodes of the corresponding surface boundaries. On the other hand, in OpenFOAM, emissions must be expressed in terms of fluxes of mass mixing ratios (in kg s⁻¹ per kilogram of mixture). Since the mass of the mixture, and therefore the corresponding mixing ratio flux, is not known at the boundary surface, this necessitates the implementation of a new boundary condition, simpleEmission, in which the mixture mass at the boundary is calculated using the adjacent interior cells, allowing surface emissions to be prescribed as mass flux boundaries. To prevent spurious accumulation of pollutant species through the spanwise cyclic boundaries, the default cyclic boundary patches in OpenFOAM have been augmented to behave as a Neumann boundary condition (zeroGradient) on the downstream boundary of the flow, while maintaining background concentrations on the upstream boundary, similar to a Dirichlet boundary condition (fixedValue). These modified boundary patches are referred to as cyclicZeroGradient and cyclicFixedValue, respectively. Finally, a preprocessing utility, initCanyon, has been developed to provide an initial velocity perturbation in the solution domain, allowing turbulence flow structures to freely develop between the spanwise and lengthwise cyclic boundaries.

Figure 2 Grid distribution of the canyon geometry in the (a) x−z plane and (b) y−z plane. Cell size in the canyon core region (dashed line) is $\mathrm{0.5}\times \mathrm{0.5}\times \mathrm{0.5}$ m³.

2.4 Domain discretization and decomposition

The solution domain is discretized using the OpenFOAM blockMesh utility, which generates a conformal hexahedral grid. A 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 an aspect ratio near unity. Grid cell grading is introduced beyond the core region so that cell size gradually expands 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 in the lengthwise directionm and 1:4 in the heightwise direction. Two mesh configurations of core cell sizes of 1 m (coarse 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 (Sect. 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 colocated 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 iteratively 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⁻⁵. Mass conservation is enforced at every iteration by solving an iterative pressure correction equation until the aforementioned residual threshold is 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).

3 Stationary evaluation

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 with a 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 – i.e., $u\left(z\right)=U_{\mathrm{\infty }}(z-H)/(z_{\mathrm{\infty }}-H)$, where z is the elevation above the canyon and z_∞ and U_∞ are the elevation and velocity at the domain freestream boundary, respectively – is assigned above the canyon. A value of 3.52 m s⁻¹ has been used for the spanwise velocity at the freestream boundary. Turbulence in the flow domain is initiated using a superposition of randomized velocity perturbation, with a magnitude of up to 10 % of the freestream velocity, into the initial flow field. Simulations are conducted using the coarse (1 m) and fine (0.5 m) mesh configurations, both at a uniform time step of 0.01 s. Each simulation is given a 6 h (21 600 s) spin-up period before sampling begins for a period of 3600 s at a frequency of 100 Hz, from which turbulent statistics are constructed.

To coincide with the LDA measurements, sampling points for the LES simulations are located along vertical stations at $x/W=-\mathrm{0.25}$, 0, and +0.25, denoting upstream, central, and downstream positions with respect to the direction of freestream flow, at the lateral plane of symmetry (y=0). The time-averaged velocity components and subgrid scale turbulent kinetic energy (SGS TKE) at each sampling point are calculated as the arithmetic mean of the respective time series. The turbulent velocity fluctuation in each direction is, in turn, determined using the L2 norm of the central difference of the instantaneous velocity components. The resolved scale TKE can then be calculated using its definition $k=(\mathrm{1}/\mathrm{2})\left({\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_i\right)$. Note that the LDA mean and fluctuation measurements represent a Reynolds average of the flow field, while the LES model data in the present study attempt to duplicate these results using only data contained in the resolved scale.

Figure 3 Normalized vertical profiles of the of resolved mean horizontal (top left) and vertical (top right) velocity components, as well as horizontal (bottom left) and vertical (bottom right) velocity fluctuations, for stations located at $x/W=-\mathrm{0.5}$ (left), 0 (middle), and +0.5 (right) in the idealized street canyon geometry. Solid lines indicate the present study using fine mesh; dashed lines indicate the present study using coarse mesh; and symbols represent the LDA measurements from Li et al. (2008).

Figure 3 shows the comparison between the present LES results on both mesh arrangements with the LDA measurements at the aforementioned vertical stations. In this figure, both the fine and coarse mesh results are in excellent agreement with respect to horizontal time-averaged velocity ($\stackrel{\mathrm{‾}}{u}$, upper-left panel), with the exception of a slight overprediction above the canyon at $x/W=\pm \mathrm{0.25}$. For the time-averaged vertical velocity component ($\stackrel{\mathrm{‾}}{w}$, upper-right panel), the fine mesh solution is able to adequately capture the profiles. There is a discrepancy at $x/W=\mathrm{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 lower-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 they 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 modeled, which could attribute, at least in part, to the observed differences.

Figure 4 Line integral convolution representation of the velocity vector along the x−z plane. Inset indicates direction of freestream flow.

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 Fig. 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=\mathrm{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 there is little mass transfer between the bulk flow and the canyon in the skimming flow regime, as indicated by the RANS vertical velocity component ($\stackrel{\mathrm{‾}}{w}$) near the canyon roof in Fig. 3. Based on the mean flow structures indicated in Fig. 4, weak vertical interaction still exists downstream in the canyon between the primary recirculation zones. It can be seen that pollutants from the street level can be carried upwards by the bulk flow, first through the lower recirculation zone in an counterclockwise direction downstream. Further exchange then takes place between the lower and upper recirculation zones, eventually reaching the canyon roof clockwise around the upper recirculation zone, where it merges with the bulk flow, again downstream in the canyon.

Figure 5 Fraction of resolved scale TKE of total (resolved and SGS) TKE for the three stations in the street canyon. Solid lines indicate fine mesh, and dashed lines indicate coarse mesh.

To further evaluate the resolution of turbulent scales under each mesh configuration, the fraction of TKE as its simulated total – i.e., the sum of resolved and SGS TKE – is calculated along the three vertical stations, as presented in Fig. 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 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 Fig. 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 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 Fig. 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.

4 Transient study with photostationary NO–NO₂–O₃ mechanism

The transient study consists of two sets of model runs. The first model set is a fully transient run, where boundary conditions are prescribed using the diurnal profiles for meteorological conditions, solar state, background concentrations, and traffic emissions. For the second set, a series of stationary runs are conducted with boundary conditions fixed at 08:00, 12:00, 16:00, and 20:00 UTC. Each run begins 6 h prior to the start of the data sampling period and is conducted using the fine mesh arrangement (core cell size 0.5 m) and a constant time step of 0.05 s. Vertical profile data are extracted at $x/W=-\mathrm{0.25}$, 0 and +0.25, and in-canyon horizontal profile data are extracted at z=2 m and $z/H=\mathrm{0.5}$ and 0.75.

Figure 6 Prescribed conditions for the transient simulation: (a) horizontal freestream wind velocity and solar zenith angle; (b) ambient temperature and density; (c) urban background NO, NO₂, and O₃ concentrations; and (d) workday traffic NO and NO₂ emissions per unit street section length.

4.1 Model configuration

Figure 6 shows the diurnal profiles for meteorology, background concentrations, traffic emissions, and solar state applied to the transient model run, representative of a typical work day for the city of Berlin in July 2014. The freestream horizontal (u) wind, as well as the solar zenith angle (χ), is presented in Fig. 6a, while the air temperature and pressure are illustrated in Fig. 6b. The profiles for wind, pressure, and temperature have been obtained from the hourly WRF model output of Kuik et al. (2018) by first interpolating the pressure coordinates to the freestream height (112 m) of the canyon and averaging them at each hour across every day in the period of July 2014. A neutral boundary layer is assumed in order to suppress lengthwise convective cells formed under thermal stratification, which introduces grid dependence in the numerical solution in a cyclic finite-length street canyon (Davies-Jones, 1970). The solar zenith angle is calculated on line at each time step. Sunrise and sunset are defined when the solar zenith angle reaches the 90^∘ crossing. For the 24 h simulation period, sunrise occurs at 02:55:06.83 UTC, and sunset takes place at 19:25:26.13 UTC. In the mean time, the solar zenith angle reaches a minimum of 29.42^∘ at 11:10:24.67 UTC.

The profiles for background concentrations of NO, NO₂, and O₃, shown in Fig. 6c, are obtained from hourly averages on all full working days (Tuesdays and Wednesdays; Mondays and Thursdays have been excluded to minimize end-effects carried over from adjacent non-working days) from observation data of urban background stations located on Amrumer Straße (DEBE010) and Nansenstraße (DEBE034) in Berlin. Days with missing measurement data are discarded entirely, leaving 17 d with intact hourly measurements. The traffic NO and NO₂ emissions, as presented in Fig. 6d, are estimated from annual NO_x aggregate for the city of Berlin in 2014 using hourly traffic data and vehicle fleet composition provided by the Berlin City Senate. Temporal redistributions of annual traffic emissions are inferred by traffic density for major street segments, from which a representative flow speed and level of service are determined. Emission factors for NO_x can then be calculated assuming a standard fleet composition for the street segment of interest. The NO_x emission is then partitioned into NO and NO₂ according to the Gridding Nomenclature for Reporting (GNFR) for the road transport source sector categories (Cranier et al., 2019), with a molar fraction of 0.95 for NO and 0.05 for NO₂ for gasoline exhaust (F1) and 0.7 (NO) and 0.3 (NO₂) for diesel exhaust (F2).

All runs are conducted with four cluster nodes, each comprising two Intel Xeon Platinum 9242 processors and 384 Gb of physical memory. Based on OpenFOAM run statistics, each simulation hour of model run requires on average 8623±148 s in processor time or 8788±145 s in wall clock time to complete, with a 95 % confidence interval over all transient and stationary runs. All output data are written to file every 30 s of simulation time, and corresponding hourly means are calculated at the top of each hour using a central window average spanning 3600 s.

Figure 7 Comparison of LIC representation of hourly mean velocity vector along the x−z plane for transient (top) and stationary (bottom) model runs at various points of time in UTC.

4.2 Results and discussion

Figure 7 shows the qualitative comparison based on the LIC representation of the hourly average velocity field between the fully transient run and the stationary runs at each specific time. As with Fig. 4, two vertically stacked primarily recirculation regions are visible in all cases, as well as secondary vortices located at the street-level corners. Differences in fluid motions can be observed between the transient and stationary model runs at the same time, as indicated by the relative displacement in the primary vortex locations. This is a result of the differences in temporal variations in the meteorological conditions, which are held constant at specific times in the stationary runs, as indicated in the difference in hourly vertical profile of horizontal velocity component in the freestream region ($z/H>\mathrm{1}$) between the transient and stationary runs, as shown in Fig. 8. The momentum exchange in and out of the canyon is restricted by the shear region, defined by a strong mean shear (${\partial }_z\overline{u}$) at the canyon roof ($z/H=\mathrm{1}$). The in-canyon ($z/H\le \mathrm{1}$) velocity component profiles follow a general feature proportional to the magnitude of the freestream velocity. Nevertheless, while the diurnal variation in freestream velocity is approximately 2 m s⁻¹, as illustrated in Fig. 6a, the corresponding change in the in-canyon velocity magnitude is in only in the order of a few centimeters.

Figure 8 Comparison of vertical station profiles of mean velocity components and species concentrations for transient and stationary model runs at various points of time in UTC. Solid lines indicate transient runs, and dashed lines indicate stationary runs.

However, deviations of up to about 20 ppb in pollutant concentrations can be observed in the vertical profiles in Fig. 8, even though the rate of hourly traffic NO and NO₂ emissions are identical for the indicated times between transient and stationary simulations. The discrepancy is particularly large at 08:00 and 20:00 UTC, This is due to a two factors contributing to the difference between the two types of model runs. First, in addition to meteorological conditions, the rates of reactions, background concentrations, and emissions are also held constant in stationary simulations, while the transient run follows a diurnal variation, which results in a difference in equilibrium composition in the street canyon and freestream regions. This also brings about the second factor; i.e., the difference in boundary conditions also leads to a difference in cumulative pollutant concentration at the beginning of the sampling period.

Two observations can be made from this comparison. First, to the extent of the meteorological conditions considered in this study, there is little change in the in-canyon flow field between stationary and non-stationary runs. This is not surprising, considering that the shear layer in the roof region attenuates the dynamics inside the canyon from those above the buildings. Accordingly, the in-canyon wind should be less sensitive to the temporal variations in the freestream flow. On the other hand, significant changes can be observed in pollutant concentrations under stationary and non-stationary treatment. As such, temporal variations in emissions and background conditions, which are not considered in stationary runs, are significant and cannot be neglected. In light of the above observations, where wind and flow distributions are of primary concern, stationary model runs would likely suffice. However, if the focus of the study is on pollutant transport and transformation, temporal variations of background concentrations and emissions will have an appreciable influence on the model concentrations, and thus a non-stationary approach should be considered.

Figure 9 Comparison of horizontal station profiles of mean velocity components and species concentrations at an elevation of 2 m for transient and stationary model runs at various points of time in UTC. Solid lines indicate transient runs, and dashed lines indicate stationary runs.

A cascade of concentrations for all species considered can be identified, corresponding to the individual primary vortices seen in Fig. 4, concurring with the numerical results of Bright et al. (2011) for a street canyon of unity aspect ratio. The concentrations of NO and NO₂ generally increase along each level of cascade (and hence each primary vortex), due to the higher flow velocity in the outer regions of each vortex, causing pollutants to be transported first to the top of the vortex. A portion of NO and NO₂ is then transferred at the top of the lower vortex into the bottom of the upper vortex through shear motion. The bulk motion in the upper vortex leads to an increase in vertical concentrations along the vortex, albeit in lower concentrations, before leaving the canyon through the shear layer at the roof level.

In addition, there is a reciprocal relationship between the concentrations of NO and NO₂, as well as O₃. This is expected as the O₃ 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₃ is not substantially different from freestream levels, at least in comparison with the NO and NO₂ levels. This is, in part, due to the in situ production of O₃ through photodissociation in Reaction (R1). The exception is 20:00 UTC, 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 UTC, leaving the titration Reaction (R3) the only active reaction in the mechanism. Combined with a low traffic emission, as indicated in Fig. 6d, this leads to a moderate decrease in O₃ while depleting the NO in the canyon.

Figure 10 Comparison of horizontal station profiles of mean velocity components and species concentrations at a normalized elevation of $z/H=\mathrm{0.5}$ for transient and stationary model runs at various points of time in UTC. Solid lines indicate transient runs, and dashed lines indicate stationary runs.

Figure 11 Comparison of horizontal station profiles of mean velocity components and species concentrations at a normalized elevation of $z/H=\mathrm{0.75}$ for transient and stationary model runs at various points of time in UTC. Solid lines indicate transient runs, and dashed lines indicate stationary runs.

Horizontal profiles of velocity and concentrations are presented in Figs. 9, 10, and 11, evaluated at 2 m, 18 m ($z/H=\mathrm{0.5}$), and 27 m ($z/H=\mathrm{0.75}$) above the canyon street level respectively, representing street level, between the two vertical primary vortices, and the upper primary vortex. The direction of each vortex is indicated by the mean horizontal and vertical velocity components, with the lower vortex traveling counterclockwise, and clockwise for the upper vortex. There are slight differences ($\sim {\mathrm{10}}^{-\mathrm{2}}$ m s⁻¹) in velocity distributions between the transient and corresponding stationary runs, although both sets of runs follow a similar velocity profile. The velocity magnitude also increases with elevation as it approaches freestream flow.

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 (Figs. 9 and 10) and the left wall in the upper vortex (Fig. 11). However, at the 2 m level, the peaks of NO_x concentration resulting from diffusion of street-level emissions, as depicted in Fig. 1c, is not evidently shown. This is indicative of the dominance of spanwise advection near the street 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₃ remains more or less stable at all horizontal levels.

Figure 12 Hourly averaged time series of NO, NO₂, and O₃ concentrations at (a–c) $z/H=\mathrm{0.75}$, (d–f) $z/H=\mathrm{0.5}$, and (g–i) z=2 m at different horizontal stations over a 24 h period. The solar zenith angle (χ) crosses the 90^∘ threshold at 02:35 and 19:25 UTC.

Finally, Fig. 12 shows the diurnal cycles in the concentrations of NO, NO₂, and O₃ in the street canyon from the transient model run at nine different locations corresponding the intersections between the vertical (Fig. 8) and horizontal profiles (Figs. 9 to 11). Both the NO and O₃ profiles exhibit a typical diurnal cycle, which show strong dependence on solar activity, since only the titration Reaction (R3) remains active before sunrise and after sunset (i.e., $\left|\mathit{\chi }\right|\ge {\mathrm{90}}^{\circ }$), as previously pointed out. Thus, in conjunction with emission rates from traffic NO_x during these times, a drastic reduction in NO can be seen. Further, the O₃ concentration at all monitoring points follows the same tendency as the background O₃ shown in Fig. 6c, except at the time immediately following the cessation of solar activity (at 19:25 UTC), where an acute reduction is observed due to NO_x still presently available, mainly from emission sources, for the titration Reaction (R3).

While the trend of NO₂ concentration in the afternoon and early evening (that is, between 15:00 and 21:00 UTC) 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 by a comparatively low background O₃ concentration, the photodissociation of NO₂ becomes the dominant reaction in the mechanism, leading to consumption of the emitted NO₂ from the traffic sources. As the O₃ concentration increases in the afternoon, the rate of titration increases correspondingly. Therefore the afternoon NO₂ peak coincides with a decrease in NO concentration.

5 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 modeling chemical transport of gas-phase pollutant species in an idealized street canyon under conditions where meteorology, background pollutant concentrations, and traffic emissions vary 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), showing excellent agreement in mean flow and turbulent velocity fluctuations. Whereas previous studies on street canyons have been conducted under stationary conditions, a subsequent transient model run with chemical transport, time-varying boundary conditions representing a typical working day for the city of Berlin in July 2014 covering a 24 h period have been prepared. Comparisons have also been made with stationary runs corresponding to certain points of time during the diurnal cycle. While vertical and horizontal profiles for velocity components indicate minor differences between the stationary runs and the transient runs at corresponding times, a large deviation in pollutant species concentration can be observed. This discrepancy is caused by time variations in ambient conditions, which also leads to differences in cumulative concentrations prior to solution sampling, both factors not considered in stationary runs. Diurnal variations in pollutant concentrations can also be inferred through the transient model run at specific locations in the canyon, showing typical profiles of NO, NO₂, and O₃ that correspond to traffic emissions, background concentrations, and solar zenith angle.

The results of this study demonstrate the ability of this solver in performing urban-scale chemistry transport modeling in 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 as meteorology and emission profiles need to be prescribed at higher spatial and temporal resolutions, which can be accomplished by coupling the urban scale model with outputs from regional models, e.g., WRF-Chem. Particular emphasis can be placed on 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 spatiotemporal distribution and NO_x emissions from annual aggregate and the subsequent partition into NO and NO₂ contributions, as outlined in Sect. 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 Sect. 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₂-O₃ mechanism (Leighton, 1961) serves as the integral component of the OpenFOAM chemistry model for the current study, modeling transformation processes of CO and VOCs will certainly require chemical kinetics that include reactions of these species and their intermediate derivatives, such as OH, R-H, and R-O₂. 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.

Appendix A: Nomenclature

A1 Roman symbols

ai	Polynomial coefficient for the ith power term
A	Arrhenius pre-exponential factor, in units of [m, kmol, s, K]
AS	Sutherland law reference viscosity factor [kg m−1 s−1 K${}^{-\mathrm{1}/\mathrm{2}}$]
Ck	Parametric constant for SGS viscosity []
Cϵ	Parametric constant for SGS TKE dissipation rate []
cp	Specific heat at constant pressure [J K−1 kg−1]
D	Binary diffusion coefficient [kg m−1 s−1]
e	Internal energy [J]
g	Gravitational acceleration [m s−2]
h	Enthalpy [J]
h	Planck constant (in Reaction R1) [J s]
i, j, k, q	Generic indices
J0	Photolysis rate at χ=0 [s−1]
jR	Photolysis rate [s−1]
k	Turbulent kinetic energy [m2 s−2]
kR	Reaction rate constant [m3 kmol−1 s−1]
m, n	Photolysis rate parameters []
p	Pressure [Pa]
R	Specific gas constant [J K−1 kg−1]
r	Net species production rate [m3 s−1]
Sij	Shear strain rate tensor [s−1]
T	Temperature [K]
TA	Arrhenius activation temperature [K]
TS	Sutherland law reference temperature [K]
U	Reference velocity [m⋅s-1]
u, v, w	Velocity components [m s−1]
u′, v′, w′	Velocity fluctuation [m s−1]
x, y, z	Cartesian coordinate [m]
yq	Mass mixing ratio for species q []

A2 Greek symbols

Δ	Subgrid filter length scale [m]
δx	Local grid size [m]
δij	Kronecker delta []
ϵijk	Permutation tensor []
λ	(Photon) wavelength [m]
μ	Dynamic viscosity [Pa s]
ν	(Photon) vibrational frequency [s−1]
ξ	Cartesian coordinate (alternative to x) [m]
ρ	Density [kg m−3]
${\mathit{\sigma }}_{ij}^{\mathit{\phi }}$	SGS term for scalar φ
τij	Shear stress tensor [Pa]
ϕ	Latitude []
φ	Generic scalar variable
χ	Solar zenith angle []
ω	Angular velocity [s−1]

A3 Mathematical operators

dx(⋅)	Direct derivative: $\mathrm{d}(\cdot )/\mathrm{d}x$
${\partial }_x(\cdot )$	Partial derivative: $\partial (\cdot )/\partial x$
${\partial }_i(\cdot )$	Partial derivative (shorthand): $\partial (\cdot )/\partial x_i$
Dt(⋅)	Substantial derivative: $d_t(\cdot )+{\partial }_i\left[u_i\right(\cdot \left)\right]$
F(x|Δ)	Spatial low-pass filter kernel with cut-off length Δ
[φ]	Volumetric concentration
$\stackrel{\mathrm{‾}}{\mathit{\phi }}$	Temporal (RANS) average
〈φ〉	Spatial average
$\stackrel{\mathrm{̃}}{\mathit{\phi }}$	Density-weighted (Favre) average: $\langle \mathit{\rho }\mathit{\phi }\rangle /\langle \mathit{\rho }\rangle$

A4 Acronyms and abbreviations

BLUME	Berlin Air Quality Measurement Network (Berliner Luftgüte-Messnetz)
CBM	Carbon bond mechanism
DNS	Direct numerical simulation
GNFR	Gridding nomenclature for reporting
ILU	Incomplete lower–upper (matrix factorization)
LDA	Laser Doppler anemometry
LES	Large-eddy simulation
PM	Particulate matter
RANS	Reynolds-average Navier–Stokes (equations)
RCS	Reduced chemical scheme
SGS	Subgrid scale
TKE	Turbulent kinetic energy
VOC	Volatile organic compound
WRF	Weather Research and Forecasting (Model)