APFoam 1.0: integrated computational fluid dynamics simulation of O₃–NO_x–volatile organic compound chemistry and pollutant dispersion in a typical street canyon

Urban air quality issues are closely related to human health and economic development. In order to investigate street-scale flow and air quality, this study developed the atmospheric photolysis calculation framework (APFoam 1.0), an open-source computational fluid dynamics (CFD) code based on OpenFOAM, which can be used to examine microscale reactive pollutant formation and dispersion in an urban area. The chemistry module of APFoam has been modified by adding five new types of reactions, which can implement the atmospheric photochemical mechanism (full O₃–NO_x–volatile organic compound chemistry) coupled with a CFD model. Additionally, the model, including the photochemical mechanism (CS07A), air flow, and pollutant dispersion, has been validated and shows good agreement with SAPRC modeling and wind tunnel experimental data, indicating that APFoam has sufficient ability to study urban turbulence and pollutant dispersion characteristics. By applying APFoam, O₃–NO_x–volatile organic compound (VOC) formation processes and dispersion of the reactive pollutants were analyzed in an example of a typical street canyon (aspect ratio $H/W=\mathrm{1}$). The comparison of chemistry mechanisms shows that O₃ and NO₂ are underestimated, while NO is overestimated if the VOC reactions are not considered in the simulation. Moreover, model sensitivity cases reveal that 82 %–98 % and 75 %–90 % of NO and NO₂, respectively, are related to the local vehicle emissions, which is verified as the dominant contributor to local reactive pollutant concentration in contrast to background conditions.

In addition, a large amount of NO_x emissions, especially NO, is beneficial to the reduction of O₃ concentrations since NO consumes O₃. Background precursors (NO_x/VOCs) from boundary conditions only contribute 2 %–16 % and 12 %–24 % of NO and NO₂ concentrations and raise O₃ concentrations by 5 %–9 %. Weaker ventilation conditions could lead to the accumulation of NO_x and consequently a higher NO_x concentration but lower O₃ concentration due to the stronger NO titration effect, which would consume O₃. Furthermore, in order to reduce the reactive pollutant concentrations under the odd–even license plate policy (reduce 50 % of the total vehicle emissions), vehicle VOC emissions should be reduced by at least another 30 % to effectively lower O₃, NO, and NO₂ concentrations at the same time. These results indicate that the examination of the precursors (NO_x and VOCs) from both traffic emissions and background boundaries is the key point for understanding O₃–NO_x–VOCs chemistry mechanisms better in street canyons and providing effective guidelines for the control of local street air pollution.

1 Introduction

With the rapid urbanization worldwide, air pollution in cities, such as haze and photochemical smog characterized by high levels of particulate matter and/or surface ozone (O₃), has become one of the most concerning global environmental problems (Lu et al., 2019; Wang et al., 2020). Recently, observational data have shown that PM_2.5, one of the major pollutants in cities, has decreased by 30 %–50 % across China due to strict air quality control measures (Zhai et al., 2019). At the same time, 87 %, 63 %, 93 %, 78 %, and 89 % of the observational stations in China have shown a decreasing trend for CO, NO₂, SO₂, PM₁₀, and PM_2.5 over the last 5 years, respectively (Fan et al., 2020). Various data indicate that air quality in China has been significantly improved. Unlike other pollutants, however, O₃ concentrations have increased in major urban clusters of China (Lu et al., 2018). Severe O₃ pollution episodes still exist and happen frequently (Wang et al., 2017). Therefore, research into reactive pollutants such as O₃, which has adverse effects on human health (Goodman et al., 2015; H. Liu et al., 2018; Sousa et al., 2013), crops (Rai and Agrawal, 2012), building materials (Massey, 1999), and vegetation (Yue et al., 2017), is of great significance to the further improvement of air quality, especially in urban areas.

From the perspective of the cause of urban air pollution, traffic-related emissions are the major part of airborne pollutant sources, including the precursors of O₃, NO_x (= NO + NO₂), and volatile organic compounds (VOCs) (Degraeuwe et al., 2017; Kangasniemi et al., 2019; Keyte et al., 2016; Pu and Yang, 2014; Wild et al., 2017; Wu et al., 2020). It is believed that the production of O₃ comes from NO₂ photolysis. Generally, in a clean atmosphere, the produced O₃ would be consumed by the NO titration effect. However, with the involved VOCs, NO concentrations become lower due to the consumption of RO₂ (the production of VOCs and OH, $\mathrm{VOCs}+\mathrm{OH}\to {\mathrm{RO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$), which weakens the NO titration effect and consequently leads to O₃ accumulation (Seinfeld and Pandis, 2016). In China, previous studies have shown that 22 %–52 % of total CO, 37 %–47 % of total NO_x, and 24 %–41 % of total VOC emissions are contributed by vehicle emissions in urban areas (Li et al., 2017; Zhang et al., 2009; Zheng et al., 2014, 2009).

Numerical simulation using air quality models is considered an effective method to investigate the formation pattern and dispersion of reactive pollutants. Based on length scales, the air flow and air quality modeling in cities are commonly categorized into four groups, i.e., street scale (∼100 m), neighborhood scale (∼1 km), city scale (∼10 km), and regional scale (∼100 km) (Britter and Hanna, 2003). Due to the complex geometry and nonuniformity in building distribution within cities, computational fluid dynamics (CFD) simulation has recently gained popularity in the urban climate research (Toparlar et al., 2017). Different from the typical mesoscale (∼1000 km) and regional-scale (∼100 km) air quality models, CFD has better performance in microscale pollutant dispersion within the urban street canyon (∼100 m) or urban neighborhoods (∼1 km), which are restricted spaces with more complicated turbulent mixing and poorer ventilation conditions than rural areas (Zhong et al., 2015). Besides the shorter physical processes in microscale urban models (∼100 m–1 km, ∼10–100 s), the rather fast chemical processes of NO₂ photolysis and NO titration with the complex chain of VOC reactions also require finer-resolution models (Vardoulakis et al., 2003). For instance, CFD models with fine grids (∼0.1–1 m) and small time steps (∼0.1 s) have been effectively adopted to simulate these high-resolution spatial and temporal variations in urban areas (Sanchez et al., 2016).

Table 1 Overview of the CFD studies with their photochemical mechanisms.

With the rapid growth of the high-performance computing (HPC) platforms, computational power is no longer an obstacle. CFD simulation shows the good application prospects for urban microclimate research (Fernandez et al., 2020; Garcia-Gasulla et al., 2020). Many CFD models coupled with photochemical reaction mechanisms have been developed to investigate the street-scale air quality problem in recent years (see Table 1). More commonly, simple photochemical mechanisms with only three reactions (Leighton, 1961) are adapted in CFD models. This mechanism can simulate the NO_x and O₃ dispersion with a lower computational requirement. Many previous studies have investigated the pivotal factors that affect the reactive pollutant distribution within the street canyon by using a CFD model with a simple photochemical mechanism, such as a street–building aspect ratio (He et al., 2017; Zhang et al., 2020; Zhong et al., 2015), ambient wind conditions (Baker et al., 2004; Merah and Noureddine, 2019), thermal effects (Baik et al., 2007), or emissions from vehicles (C. W. Liu et al., 2018; Y. Zhang et al., 2019). However, due to the simple photochemical mechanism ignoring the effect of other nitrogen oxides and VOCs on the photochemistry, some studies have recently applied the full photochemical mechanism in CFD models to reduce the uncertainty of pollutant simulation. Photochemical mechanisms contain NO_x–O₃–VOC reactions and photochemistry, such as CBM-IV (Garmory et al., 2009; Kwak et al., 2013; Kwak and Baik, 2012, 2014), GEOS-Chem (Kim et al., 2012; Park et al., 2016), RCS (Bright et al., 2013; Zhong et al., 2017), and CCM (Sanchez et al., 2016) and are successfully coupled with CFD models and applied to analyze the street-scale pollutant dispersion.

Currently, most of the simulation studies have been carried out via the application of commercial CFD software. This software is rather simple to operate, which is effective in saving time when setting up the simulation case. However, the commercial codes are usually closed source, which is a “black box” for users (Chatzimichailidis et al., 2019). In this case, adjustments to the equations and parameters or modifications to the model are difficult for some specific simulations. Therefore, an open-source CFD code for atmospheric photolysis calculation, APFoam 1.0, was developed in this study. Open-Source Field Operation and Manipulation (OpenFOAM) was selected as the platform for the APFoam framework, as OpenFOAM has good performance regarding computing scalability and low uncertainty levels, which shows its applicability for large-scale CFD simulations with million-level grid numbers (Robertson et al., 2015). In addition, the solvers in OpenFOAM for specific CFD problems can be developed by using the appropriate packaged functionality, which simplifies the difficulty of programming. Furthermore, OpenFOAM has also been developed with various pre- and post-processing utilities that are convenient for data manipulation (OpenFOAM Foundation, 2018).

This paper is organized as follows. Section 2 presents a full description of the new chemistry module and simulation solver. Model validation for the photochemical mechanism, turbulence simulation, and pollutant dispersion compared with the chemical box model and several wind tunnel experiments is discussed in Sect. 3. In Sect. 4 a series of sensitive cases have been set up to investigate the contribution of the key factors to the reactive pollutants in a typical street canyon (aspect ratio of building height to street width, $H/W=\mathrm{1}$). Finally, the conclusion and future research plans for the APFoam framework are summarized in Sects. 5 and 6.

2 Model description

2.1 General overview

The APFoam framework has been developed based on OpenFOAM, which is an open-source code for CFD simulation. For the numerical solution, APFoam uses finite-volume method (FVM) to discretize the governing equations and adopts arbitrary three-dimensional structured or unstructured meshes. All variables of the same cell are stored at the center of the control volume (CV), and complex geometries can be easily handled with FVM (Chauchat et al., 2017). In APFoam, laminar, Reynolds-averaged Navier–Stokes equation (RANS), and large-eddy simulation (LES) methods are available for turbulence solution. Additionally, APFoam also has complete boundary conditions to choose from for numerical simulations.

Based on OpenFOAM, APFoam has been developed to conduct photochemical simulations within the atmosphere. Different from the general chemical reaction types, there are some new types of gaseous reactions for describing photochemical processing. More details will be introduced in Sect. 2.2.

Figure 1 Flow diagram of the simulation setup in the APFoam framework.

To make it easier to get started with APFoam, the structure of the simulation case folder is consistent with OpenFOAM. Figure 1 shows the flow diagram of the simulation setup in APFoam. The solver of APFoam, APChemFoam for one-dimensional (1D) chemistry solving, is modified from the solver ChemFoam, which will be introduced in detail in Sect. 2.3. Other three-dimensional (3D) solvers are modified from the solver reactingFoam, including APreactingFoam, solving flow field and chemical reactions simultaneously in one time step; APonlyChemReactingFoam, solving only chemical reactions with a certain flow field; and APsteadyReactingFoam, solving flow field and chemical reactions simultaneously in steady state. More details will be presented in Sect. 2.4.

For running the simulation (see Fig. 1), mesh files, configuration files, and initial and boundary condition files should be prepared before the simulation. Mesh files can be made in various ways, such as making a blockMesh application executable by using data from blockMeshDict or fluentMeshToFoam and converting the .msh file to OpenFOAM format. For the APFoam simulation, all required configure files are also listed in Fig. 1. In addition, user-defined functions can be loaded during the run time without recompiling the program via writing-related configuration files in the system folder. As for initial conditions, the initial states of turbulence, the environment (e.g., temperature, pressure), and chemical species are necessary for the simulation. The results of APFoam contain the wind flow and pollutant concentrations, which can be processed by the Paraview in OpenFOAM or any other CFD post-processing tools.

2.2 Chemistry module

For photochemical calculation, there are five types of reactions in the NO_x–O₃–VOC mechanism. In the original version of OpenFOAM, these types of reactions are not included and should be added to the chemistry module prior to simulation. These types of the reactions are described as follows.

1. Arrhenius reactions.

   Arrhenius reactions are the basic reaction in the mechanism, and the rate of the Arrhenius reaction is calculated as

   $$\begin{array}{}\text{(1)}& k=A\cdot {\left({\displaystyle \frac{T}{\mathrm{300}}}\right)}^B\cdot \exp \left(-{\displaystyle \frac{E}{T}}\right),\end{array}$$

   where A, B, and E are the parameters of the reaction rates and T is the temperature of the mixture in degrees kelvin.

2. Photolysis reactions.

   Photolysis reactions are first-order reactions, and the photolysis rate is calculated as

   $$\begin{array}{}\text{(2)}& k_{\mathrm{phot}}=\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}J\left(\mathit{\lambda }\right)\cdot \text{abs}\left(\mathit{\lambda }\right)\cdot \text{QY}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda },\end{array}$$

   where k_phot is the first-order rate for the photolysis reaction; λ₁ and λ₂ are the photolysis wavelength ranges according to the specific species; and J (λ), abs (λ), and QY (λ) are the intensity of the light source, absorption cross section, and the quantum yield for the reaction at wavelength λ, respectively.

   In reality, the photolysis rate could be calculated by other photolysis rate models, such as Fast-J (Wild et al., 2000) or TUV (Madronich and Flocke, 1999), or obtained from a photolysis data set, such as IUPAC (Atkinson et al., 2004). Since the photolysis rate does not depend on temperature, in the current version of APFoam the model does not consider the variation of light intensity, and the photolysis rates are obtained from the literature (Carter, 2010) rather than online calculation in order to improve calculation efficiency.

3. Falloff reactions.

   The rate of falloff reactions is a function of temperature and pressure that is calculated as follows:

   $$\begin{array}{}\text{(3)}& k\left(T,M\right)=\left\{{\displaystyle \frac{k_{\mathrm{0}}\left(T\right)\cdot \left[M\right]}{\mathrm{1}+\frac{k_{\mathrm{0}}\left(T\right)\cdot \left[M\right]}{k_{\inf }\left(T\right)}}}\right\}\cdot F^z,y\end{array}$$

   where $z={\left\{\mathrm{1}+{\left[\frac{{\log }_{\mathrm{10}}\mathit{\left\{}\frac{k_{\mathrm{0}}\left(T\right)\cdot \mathit{\left\{}M\mathit{\right\}}}{k_{\inf }\left(T\right)}\mathit{\right\}}}{N}\right]}^{\mathrm{2}}\right\}}^{-\mathrm{1}}$, [M] is the concentration of third body, which depends on total pressure, and F is the broadening factor. k₀ and k_inf are the rates of the Arrhenius form at the low-pressure limit and high-pressure limit, respectively.

4. Three-k reactions.

   The rate of “three-k” reactions depends on three reaction rates in Arrhenius form. The rate is calculated as

   $$\begin{array}{}\text{(4)}& \begin{array}{rl}k\left(T,M\right)& =k_{\mathrm{0}}\left(T\right)+k_{\mathrm{3}}\left(T\right)\cdot \left[M\right]\\ & \cdot \left(\mathrm{1}+{\displaystyle \frac{k_{\mathrm{3}}\left(T\right)\cdot \left[M\right]}{k_{\mathrm{2}}\left(T\right)}}\right),\end{array}\end{array}$$

   where k₀, k₂, and k₃ are the three reaction rates and [M] is the concentration of third body.

5. Two-k reactions.

   The rate of “two-k” reactions depends on two reaction rates in Arrhenius form. The rate is calculated as

   $$\begin{array}{}\text{(5)}& k\left(T,M\right)=k_{\mathrm{1}}\left(T\right)+k_{\mathrm{2}}\left(T\right)\cdot \left[M\right],\end{array}$$

   where k₁ and k₂ are the two reaction rates and [M] is the concentration of third body.

In the current version of APFoam, two atmospheric photochemical mechanisms are included in the model, SAPRC07 (Carter, 2010) and CB05 (Yarwood et al., 2005). For SAPRC07, two versions of the chemical mechanism are available, which are CS07A and SAPRC07TB. CS07A is one of the condensed versions of the mechanism, which contains 52 species and 173 reactions. SAPRC07TB is a more complicated version and even contains toxic species, with 141 species and 436 reactions. As for CB05, a basic version with 51 species and 156 reactions is optional in the model. In Sect. 3.1, CS07A has been validated, while the other two mechanisms are not verified in this study but are still an available option for users.

2.3 One-dimensional chemical solver (APchemFoam)

In the APFoam framework, a one-dimensional chemistry solver (i.e., chemistry box model) called APchemFoam is included in the model. This solver only concerns the chemical concentration and reaction heat variation during simulation, and calculations are started from initial conditions within a single cell mesh. The concentration and energy equation are described as follows (OpenFOAM Foundation, 2018):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \mathit{\rho }{Y}_i}{\partial t}}=k_i\left(Y_i,T\right),\text{(7)}& {\displaystyle }h=u_{\mathrm{0}}+{\displaystyle \frac{p}{\mathit{\rho }}}+\underset{\mathrm{0}}{\overset{t}{\int }}{\displaystyle \frac{\dot{q}}{\mathit{\rho }}}\mathrm{d}\mathit{\tau },\text{(8)}& {\displaystyle }h=\sum Y_i\left(\mathrm{\Delta }{h}_{\mathrm{f},i}^{\mathrm{0}}+\underset{T_{\mathrm{0}}}{\overset{T}{\int }}{C}_{\mathrm{p},i}\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }\right),\text{(9)}& {\displaystyle }p={\displaystyle \frac{\mathit{\rho }RT}{M_{\mathrm{ave}}}}=\sum p_i=\sum {\displaystyle \frac{Y_i}{M_i}}\mathit{\rho }RT,\end{array}$$

where Y_i is the species mass fraction, k_i is the reaction rate, T is the temperature of the mixture, h is the specific enthalpy, u₀ is the initial energy, p is the pressure, ρ is the density of the mixture, $\dot{q}$ is the heat from reaction, $\mathrm{\Delta }{h}_{\mathrm{f},i}^{\mathrm{0}}$ and C_p,i are the enthalpy of formation at reference temperature T₀ and the constant pressure-specific heat (a function of temperature) of species i, R is the gas constant, M_ave is the average molar weight, and p_i and M_i are partial pressure and the molar mass of species i. In addition, either p or ρ should be set as a constant for the simulation according to the needs of research. The other is calculated by Eq. (9).

2.4 Three-dimensional (3D) CFD solver with photochemical reaction

As mentioned above, three 3D solvers for atmospheric photochemical CFD calculation, including APreactingFoam, APonlyChemReactingFoam, and APSteadyReactingFoam, are developed in the APFoam framework.

For APreactingFoam, flow field, chemical reaction, and pollutant dispersion are solved simultaneously in the same time step in this solver. Firstly, the continuity-governing equation in this solver is

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }U\right)=\mathrm{0}.\end{array}$$

Besides, the momentum governing equation is

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathit{\rho }U}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }UU\right)-\mathrm{\nabla }\cdot \mathit{\tau }=-\mathrm{\nabla }p.\end{array}$$

Additionally, the energy governing equation is

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\rho }h}{\partial t}}& +\mathrm{\nabla }\cdot \left(\mathit{\rho }Uh\right)+{\displaystyle \frac{\partial \mathit{\rho }K}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }UK\right)\\ & -\mathrm{\nabla }\cdot \left({\mathit{\alpha }}_{\mathrm{eff}}\mathrm{\nabla }h\right)={\displaystyle \frac{\partial p}{\partial t}}+\dot{q},\end{array}\end{array}$$

where U is the velocity vector of the air flow, τ is the viscous stress tensor, K is the specific kinetic energy, and α_eff is the effective thermal diffusivity coefficient.

Pressure–velocity coupling schemes for solving the flow field use the PIMPLE algorithm, a merged PISO–SIMPLE algorithm in OpenFOAM toolkit. This algorithm uses a steady-state solution (SIMPLE algorithm) for the flow field within the time step. When the defined tolerance criterion is reached, this algorithm uses the PISO algorithm in the outer correction loop and moves on in time (Holzmann, 2017). The PIMPLE algorithm allows for larger Courant numbers (Co>1) so that the time step can be increased to reduce the computation time. Even so, the time step (Δt) generally follows the Courant–Friedrichs–Lewy (CFL) condition to maintain numerical stability, which is as follows:

$$\begin{array}{}\text{(13)}& \mathit{Co}={\displaystyle \frac{U\mathrm{\Delta }t}{\mathrm{\Delta }x}}\le \mathrm{1},\end{array}$$

where Δx is the grid size.

Besides, the governing equation for the reactive species transportation is

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \mathit{\rho }{Y}_i}{\partial t}}=-\mathrm{\nabla }\cdot \left(\mathit{\rho }U{Y}_i\right)+\mathrm{\nabla }\cdot \left(\mathit{\mu }\mathrm{\nabla }{Y}_i\right)+{\left[\mathrm{\Delta }{Y}_i\right]}_{\mathrm{chem}}+E_i,\end{array}$$

where μ is the kinematic viscosity and ${\left[\mathrm{\Delta }{Y}_i\right]}_{\mathrm{chem}}$ is the concentration change of species i from the chemical reaction. As mentioned above, ${\left[\mathrm{\Delta }{Y}_i\right]}_{\mathrm{chem}}$ is calculated following Eqs. (6)–(9). E_i is the emission source of species i. The chemistry is solved by the ordinary deferential equation (ODE) solvers in OpenFOAM library, in which the chemical reactions can be integrated by automatically dividing the flow time step into several time sub-steps. The APreactingFoam is only designed to solve compressible fluids because the simulation results are more likely to be unstable and divergent when the chemistry and flow field are solved simultaneously under the incompressible fluids. This is a limitation of OpenFOAM code that is widely known by its users.

APonlyChemReactingFoam is only capable of solving the chemical reaction and species dispersion in the same time step under a certain flow field. The solution of turbulent fluids governing the equation is switched off. The purpose of developing this solver is to save computation time and reduce repetitive simulations. In general, the atmospheric chemical reactions have negligible effect on the flow field. Therefore, when the example cases to be studied do not involve the flow field change, this solver is suitable for this kind of simulation. The governing equation for the reactive species transportation is consistent with APreactingFoam (Eq. 14).

APSteadyReactingFoam is developed for solving the chemical reaction and species dispersion under the steady-state flow field. This solver is only designed to solve compressible fluids for the same reason as APreactingFoam. In this solver, a pressure–velocity coupling scheme switches to a SIMPLE algorithm for steady-state solution. The continuity governing equation in this solver is

$$\begin{array}{}\text{(15)}& \mathrm{\nabla }\cdot \mathit{\rho }U=\mathrm{0}.\end{array}$$

Besides, the momentum governing equation is

$$\begin{array}{}\text{(16)}& \mathrm{\nabla }\cdot \left(\mathit{\rho }UU\right)=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left(\mathit{\mu }\mathrm{\nabla }U\right).\end{array}$$

In addition, the energy governing equation is

$$\begin{array}{}\text{(17)}& \mathrm{\nabla }\cdot \left(\mathit{\rho }Uh\right)+\mathrm{\nabla }\cdot \left(\mathit{\rho }UK\right)-\mathrm{\nabla }\cdot \left({\mathit{\alpha }}_{\mathrm{eff}}\mathrm{\nabla }h\right)=\dot{q}.\end{array}$$

As for reactive species, the governing equation for the transportation still applies Eq. (14) as well in order to ensure the stability of chemical reaction.

3 Model validation

3.1 Photochemical reaction mechanism

To verify the accuracy of the chemical reaction solution and species concentration calculation, APFoam results are compared with the results from SAPRC box modeling software (Carter, 2010). For the chemical mechanism, CS07A is selected for validation in this study, and the simulation time is set to 24 h without diurnal variation (i.e., the chemical reaction rate is constant during the simulation), allowing the reactants to fully react and verifying the stability of the model.

Figure 2 The comparison of concentration results between APFoam and the SAPRC box model.

Figure 2 shows the concentrations of 52 species from two models at 24 h, which is the last time step of the simulation. In general, APFoam results have a good agreement with the SAPRC box model. The simulation results for other species from two models are basically consistent, except that some species have large errors when the magnitude is very small (Fig. 2f–h).

For further investigation, relative error (RE, %) for each species at each time step and mean relative error (MRE, %) are calculated for the selected species with large bias. These statistics are calculated as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\text{RE}}_{i,t}={\displaystyle \frac{C_{\mathrm{APFoam},t}-C_{\mathrm{SAPRC},i,t}}{C_{\mathrm{SAPRC},i,t}}}\times \mathrm{100}\mathit{\%},\text{(19)}& {\displaystyle }{\text{MRE}}_i={\displaystyle \frac{\sum _{t=\mathrm{1}}^n\left|\frac{C_{\mathrm{APFoam},i,t}-C_{\mathrm{SAPRC},i,t}}{C_{\mathrm{SAPRC},i,t}}\right|}{n}}\times \mathrm{100}\mathit{\%},\end{array}$$

where RE_i,t is the relative error of species i at time step t, $C_{\mathrm{APFoam},i,t}$ is the concentrations of species i at time step t from APFoam, $C_{\mathrm{SAPRC},i,t}$ is the concentrations of species i at time step t from SAPRC box model, MRE_i is the mean relative error of species i, and n is the total number of the time step.

Figure 3 The time series of (a) relative error (%) and (b) concentrations (ppmv) of the six species with largest bias.

Overall, most of the RE_i,t are less than 1 % in the concentration range between 0 to 10⁻²⁰ ppmv (i.e., the concentrations under realistic conditions, ∼0 to 10⁻¹⁷ ppbv), indicating that the simulation error of APFoam is less than 1 % during the whole simulation period. However, there are six species with RE and MRE greater than 1 %: TERP, ISOPRENE, OLE1, OLE2, IPRD, and ARO2. The MRE of these six species are 44.0 %, 40.7 %, 7.74 %, 38.5 %, 7.71 %, and 1.20 %, respectively. Additionally, Fig. 3 shows time series of RE and concentrations for these high RE and MRE species. In Fig. 3a, RE values in the early stage of the simulation (t=0–180 min) are less than 1 % for these species. However, the REs of TERP, ISOPRENE, and OLE2 increase dramatically after t=180 min. The REs can even reach up to 190.4 %, 297.0 %, and 867.4 %, respectively, in the following simulation. It should be noted that at the later time the REs of these three species have no values because they are consumed during the chemical reaction and their concentrations from SARPC box model become zero. The significant increase in the RE values of OLE1 and IPRD begins at t=1020 min, with maximum RE values of 60.1 % and 60.9 %, respectively. Relatively, the RE of ARO2 is smaller, with a value of 5.0 %.

Figure 3b illustrates the concentration variation of these six species with the worst agreement from two models. It can be found that the concentrations of these six species keep dropping during the whole simulation period. Combined with the result of Fig. 3a, the dramatic increase of RE is due to the significant concentration decrease of these six species. In this study case, these six species are continuously consumed without supplement, which results in the concentrations of these species tending towards 0. For extremely small numbers, the processing of different model is diverse. Thus, when the reduction of magnitude exceeds 10⁻⁵ to 10⁻⁶ ppbv, the RE would become much larger between the two models. In the realistic situations, the concentrations of the species would not be completely consumed with continuous emission sources and boundary conditions. Therefore, the photochemical reaction simulation results of APFoam could be reliable, and the overall errors might be less than 1 %.

3.2 Numerical settings and validation studies in urban flow modeling

It is well known that large-eddy simulations (LES) perform more accurately when simulating urban turbulent characteristics than the Reynolds-averaged Navier–Stokes (RANS) simulations. However, RANS models (e.g., k-ε models) are still more widely utilized because of the disadvantages of LES, such as their much higher computational time and resource requirements, the difficulties in setting appropriate wall boundary conditions and defining the time-dependent domain inlet, and the challenges in developing advanced sub-grid-scale models. Among the RANS turbulence models, in contrast to the modified k-ε models (e.g., realizable and re-normalization group (RNG) k-ε models), although the standard k-ε model performs worse in predicting turbulence in the strong wind region of urban districts (e.g., separate flows near building corners), the prediction accuracy is better when simulating the low-wind-speed region (e.g., weak wind in a 2D street canyon sheltered by buildings at both sides) (Tominaga and Stathopoulos, 2013; Yoshie et al., 2007). Hence, as one of the widely adopted RANS methods, the standard k-ε model is selected to solve the incompressible steady-state turbulent flows in a 2D street canyon.

Figure 4 (a) Wind tunnel experiment in a 2D street canyon. (b) The single street canyon CFD domain setups in the scaled model, (c) the inlet profile (Line E), and measurement profiles (Line F) in the street canyon.

To further evaluate the numerical accuracy of the turbulence flow simulation, a scaled CFD case is performed under the estimation of wind tunnel data. In the wind tunnel experiments (Fig. 4a), in total 25 rows of building models are set along the wind direction with the working section that is 11 m long, 3 m wide, and 1.5 m tall. For each row, building height (H), building width (B), and street width (W) are 12, 5, and 5 cm (i.e., aspect ratio $H/W=\mathrm{2.4}$), respectively. The span-wise (or lateral) length is $L=\mathrm{1.25}\mathrm{m}>\mathrm{10}H$, which is sufficiently long to ensure the 2D flow characteristics in the street canyon (Hang et al., 2020; Oke, 1988; K. Zhang et al., 2019); i.e., the flow in the targeted street region is determined by the external flow above it but includes small impacts from the lateral boundaries. Free-flow wind speed in the wind tunnel experiment is 13 m s⁻¹.

Figure 4b and c show the schematic diagrams of the CFD simulation domain setting for a single full-scale street canyon simulation. H and W of the street canyon are set as 24 and 10 m, with spatial scale ratio of 200:1 compared to the wind tunnel experiment. The corresponding Reynolds number ($\mathit{Re}=\frac{U_{\mathrm{ref}}H}{\mathit{\upsilon }}$) in the full-scale flow CFD validation ($H/W=\mathrm{2.4}$, H=24 m) is about 2.14×10⁷ and that in the wind-tunnel-scale experiments ($H/W=\mathrm{2.4}$, H=0.12 m) is 1.9×10⁵, which satisfies the requirement of Reynolds number independence (the critical level is about 8.7×10⁴ with the $H/W$ of 2) (Chew et al., 2018; Yang et al., 2020, 2021). The normalized wind profiles with two scales can be compared for validation purposes. Such validation techniques have been adopted in the literature (Hang et al., 2020; Yang et al., 2021). Besides, the building width B and L_y in the CFD simulation are 10 and 3.2 m (2$H/\mathrm{15}$), respectively, assuming that only a section ($L_y=\mathrm{2}H/\mathrm{15}$) of long street canyons adopted with symmetry conditions is applied at two lateral boundaries. The minimum grid size in this case is 0.2 m, with an expansion ratio of 1.2 from the wall surface toward the around it, which refers to the grid independence tests from our previous research (K. Zhang et al., 2019). The upstream domain inlet profiles along Line E and comparison of profiles along Line F (Fig. 4c) are measured by the laser Doppler anemometry (LDA) system in wind tunnel tests. Additionally, CFD inlet profiles of stream-wise velocity (u) and turbulent kinetic energy (TKE) are fitted following the profiles in experimental data (Fig. 5).

Figure 5 The inlet profile of (a) stream-wise velocity and (b) turbulent kinetic energy in a single street canyon case.

All governing equations for the flow and turbulent quantities are discretized by FVM, and the SIMPLE scheme is used for the pressure and velocity coupling. The under-relaxation factors for the pressure term, momentum term, k term, and ε term are 0.3, 0.7, 0.8, and 0.8, respectively. CFD simulations do not stop until all residuals become constant. Typical residuals at convergence are $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$, $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$, and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ for U_x, U_y, and U_z, respectively; $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ for continuity; $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ for k; and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ for ε.

Figure 6 The stream-wise velocity profiles along Line F.

Figure 6 shows the stream-wise velocity profiles of simulation results and experimental data along the centerline (Line F) of the street canyon. The predicted wind profile agrees well with the wind tunnel data. One main vortex structure is formed in the street canyon. The center of the main velocity (i.e., where stream-wise velocity is 0) also matches well between simulation and experiment.

Furthermore, some statistical parameters, including normalized mean-square error (NMSE), fractional bias (FB), and correlation coefficient (R) are calculated by the following equations:

$$\begin{array}{}\text{(20)}& {\displaystyle }\text{NMSE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(O_i-P_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{O}_i{P}_i}},\text{(21)}& {\displaystyle }\text{FB}={\displaystyle \frac{\mathrm{2}\left(\overline{O}-\overline{P}\right)}{\overline{O}+\overline{P}}},\text{(22)}& {\displaystyle }R={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left[\left(O_i-\overline{O}\right)\left(P_i-\overline{P}\right)\right]}{{\left[{\sum }_{i=\mathrm{1}}^n{\left(O_i-\overline{O}\right)}^{\mathrm{2}}\right]}^{\mathrm{0.5}}{\left[{\sum }_{i=\mathrm{1}}^n{\left(P_i-\overline{P}\right)}^{\mathrm{2}}\right]}^{\mathrm{0.5}}}},\end{array}$$

where n is the total number of measurement points, O_i is the experimental data at measurement point i, P_i is the CFD result at measurement point i, $\overline{O}$ is the mean value of experimental data at all points, and $\overline{P}$ is the mean value of CFD results at all points. According to the previous studies (Chang and Hanna, 2005; Sanchez et al., 2016), the model acceptance criteria for an urban configuration are NMSE < 1.5, $-\mathrm{0.3}<\text{FB}<\mathrm{0.3}$, and R>0.8. In this simulation case, the respective NMSE, FB, and R are 0.01, −0.04, and 0.99, respectively (Table 2), which shows the good performance of APFoam in flow field simulation.

Table 2 Static values of the turbulence flow simulation.

Figure 7 The schematic diagram of the 2D pollutant dispersion simulation setting.

3.3 Pollutant dispersion in a 2D street canyon

Currently, there are rarely wind tunnel experiments with chemical reactions. Thus, the pollutant dispersion accuracy in a 2D street canyon is validated by wind tunnel experimental data with tracer gas (Meroney et al., 1996), following previous studies (He et al., 2017; Zhang et al., 2020). The wind tunnel and the CFD domain configuration are presented in Fig. 7. A total of 28 rows of the wooden bar with 27 street canyons are set from upstream toward downstream along the inflow, and the street axis is perpendicular to the wind direction. Both the height (H) and width (B) of the bar are 0.06 m, and the street canyon width (W) is also 0.6 m, i.e., the aspect ratio ($H/W$) is 1 in this study case. A pollutant line source of ethane (C₂H₆) is set to emit the pollutant in the targeted street canyon. Following the wind tunnel configuration, there are 20 bars upstream and 8 bars downstream of the targeted street canyon. Eight measurement points are set in the targeted street canyon, with four (P4, P5, P6, P7) of them on the leeward side and the other four (P11, P12, P13, P14) on the windward side. The positions of the measurement points are demonstrated in Fig. 7. Pollutant concentrations at each measurement point are normalized with respect to that of the P7 (C_i/C₇) within the street canyon (Sanchez et al., 2016; Santiago and Martín, 2008). For the CFD simulation, the APreactingFoam solver with the standard k-ε model is applied to solve the compressible unsteady-state turbulent flow field and pollutant dispersion. In order to be consistent with the wind tunnel experiment setting, the photochemical mechanism is not used in the simulation. The minimum grid size in this case is 0.5 mm, with an expansion ratio of 1.1 from the wall surface toward its surroundings, and the inlet velocity is constant at 3 m s⁻¹ in the simulation. The time step of the simulation is set as $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ s in this validation case.

Figure 8 Normalized concentrations of CFD and experimental data at each measurement point in the 2D dispersion case.

As a result of the comparison, Fig. 8 shows the normalized concentrations between the CFD simulation and experimental data. In general, the model slightly overestimates the C₂H₆ concentrations on the windward side. However, at P4 the model concentrations for the top of the leeward side are lower than the experimental data, and the simulation results at P5 overestimate the concentrations of the pollutants. In this simulation case, the respective values of NMSE, FB, and R are 0.06, −0.13, and 0.95 (Table 3), which shows the good performance of APFoam in 2D pollutant dispersion simulation.

Table 3 Static values of the 2D pollutant dispersion simulation.

Figure 9 (a) The simulation domain of 3D pollutant dispersion and (b) the measurement points' locations in the street canyon.

3.4 Pollutant dispersion in a 3D street canyon

As mentioned in Sect. 3.3, 3D pollutant dispersion validation with tracer gas is conducted in this study, following the previous study (Y. Zhang et al., 2019). Simulation results are also compared with the wind tunnel experimental data (Chang and Meroney, 2001). The CFD domain configuration is presented in Fig. 9a. In this case, six buildings are set in the domain. Building height (H) and street canyon width (W) are both 0.08 m with $H/W=\mathrm{1}$. Building length (L_x) and building width (L_y) are 0.276 and 0.184 m, respectively. The distance between the buildings and the domain inlet, side boundary, top boundary, and domain outlet is 5H, 5H, 10H and 15H, respectively, for simulating realistic results (Tominaga et al., 2008). Within the target street canyon, there are also eight measurement points (four are on the leeward side, and four are on the windward side) for measuring the concentrations (Fig. 9b). In addition, six more measurement points are also set on the top of the downstream building. Pollutant concentrations at each measurement point in this simulation case are normalized with respect to the P5 (C_i/C₅) within the street canyon. The source of the C₂H₆ is set as an inlet at the bottom of the target street canyon. The size of the source is 0.005 m in width and 0.092 m in length, and it is set in the middle of canyon. The release velocity is 0.01 m s⁻¹ toward the top boundary and the mass fraction of the C₂H₆ is 1 (pure gas of C₂H₆). For the 3D pollutant dispersion simulation, an APreactingFoam solver with a standard k-ε model is applied to solver-compressible unsteady-state turbulent flow and pollutant dispersion. The photochemical mechanism is not used in the simulation. The minimum grid size in this case is 0.0005 m with an expansion ratio of 1.1 from the wall surface toward the surrounding area. The time step of the simulation is set as $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ s in this validation case as well. Meanwhile, the inlet velocity and TKE profile are also retrieved from and fitted by the experimental data (Fig. 10).

Figure 10 The inlet profile of (a) stream-wise velocity and (b) turbulent kinetic energy in the 3D dispersion case.

Figure 11 Normalized concentrations of CFD and experimental data at each measurement point in the 3D dispersion case.

Figure 11 shows the comparison results between the CFD simulation and experimental data. Overall, the CFD simulation in the 3D dispersion case slightly overestimates the concentrations in the street canyon. As for P23 and P24, the simulated results also overestimate the concentrations, as they are affected by the higher concentrations predicted within the street canyon. Similarly, statistical variables such as NMSE, FB, and R are calculated to evaluate the performance of the model. As shown in Table 4, the value of NMSE, FB, and R is 0.16, −0.21, and 0.93 in the 3D dispersion case, respectively, which agrees with the acceptance criteria. In general, APFoam also shows the good performance of the 3D pollutant dispersion simulation.

Figure 12 Schematic diagram of (a) the CFD simulation domain and (b) the probe point locations.

Table 4 Static values of the 3D pollutant dispersion simulation.

4 Numerical results in the case study

4.1 Simulation configuration and CFD setting

In this study, APFoam with a CS07A photochemical mechanism is applied for the street air quality simulation. As shown in Fig. 12a, the street aspect ratio ($H/W$) is set with building height (H=24 m), street width (W=24 m), and span-wise street length (L=30 m). A telescoping multigrid approach is adopted in the simulation with a minimum grid size of 0.2 m and an expansion ratio of 1.2 from the building walls to the surrounding area. The total grid number is about 87 300 for the whole CFD domain. The top and two lateral boundaries of the domain are set up as the symmetry boundary conditions. The emissions area is set up at the bottom of street canyon with a pollutant source size of 18 m (width, W_E) × 30 m (length, L_E) × 0.3 m (height, H_E), representing traffic emissions near street level, and emissions data are obtained from our previous work (Wu et al., 2020). In this study, the emissions of NO_x, VOCs, and CO are $\mathrm{4.37}\times {\mathrm{10}}^{-\mathrm{8}}$, $\mathrm{2.34}\times {\mathrm{10}}^{-\mathrm{8}}$, and $\mathrm{2.03}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m⁻³ s⁻¹ (i.e., ∼35, ∼200, and ∼170 ppbv s⁻¹), respectively. The NO and NO₂ are separated from NO_x by a ratio of 9:1, which is similar to the previous study (Baik et al., 2007). VOCs are speciated following the SAPRC mechanism, and the emission fraction of the species is obtained from the literature (Carter, 2015).

Figure 12b shows the probe point locations for the numerical case, wherein temporal variations of reactive pollutant concentrations are monitored, which include three points at a pedestrian height of z=1 m (near the street bottom at the leeward side (LB), center (CB), and windward side (WB)) and two other probe points near the street top (ST, z=23 m) and street center (SC, z=12 m = 0.5H).

A power-law velocity vertical profile is adopted for the inflow boundary condition, which is described as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }{U}_{\mathrm{in}}\left(z\right)=U_{\mathrm{ref}}\times {\left({\displaystyle \frac{z-H}{z_{\mathrm{ref}}}}\right)}^{\mathit{\alpha }},\text{(24)}& {\displaystyle }{k}_{\mathrm{in}}\left(z\right)={\left(U_{\mathrm{in}}\left(z\right)\times I_{\mathrm{in}}\right)}^{\mathrm{2}},\text{(25)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{in}}\left(z\right)={\displaystyle \frac{C_{\mathit{\mu }}^{\frac{\mathrm{3}}{\mathrm{4}}}k{\left(z\right)}_{\mathrm{in}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathit{\kappa }z}}.\end{array}$$

Here the reference velocity U_ref is 3 m s⁻¹, the reference height z_ref is 24 m, the turbulence intensity I_in is 0.1, the power-law exponent α is 0.22 (He et al., 2017; K. Zhang et al., 2019, 2020), the Von Kármán constant κ is 0.41, and the C_μ is 0.09.

In addition, the initial and inlet background concentrations for O₃, NO, NO₂, VOCs, and CO are 60, 5, 15, 40, and 400 ppbv, respectively, which are obtained from an observation campaign (Liu et al., 2008). For meteorological conditions, the temperature is 300 K and the operating pressure is 1013.25 hPa.

In all simulation cases, the steady-state turbulence field is first solved in advance. The result of turbulent flow drives the chemistry solution from t=0. During t=0–30 min (1800 s), the emission and chemistry solution are turned on under the statistically steady turbulent flow, reaching a quasi-dynamic and photostationary steady state. Data from the next 60 min (t=30–90 min, 1800–5400 s) are used for analysis. The time step of the chemistry solution is set as 0.1 s in all numerical cases.

Table 5 Description of all simulation cases.

A description of all simulation cases is given in Table 5. All CFD simulations are finished with the Tianhe II supercomputer and supported by National Supercomputer Center in Guangzhou. To investigate the effect of the chemical mechanism, the background condition of the precursors (BC), emissions (Emis), and wind conditions (U_ref) of the reactive pollutant concentrations in the street canyon, the cases of BC_zero_out, Emis_zero_out, and Uref0.5 are set up in numerical simulations. In Case_BC_zero and Case_Emis_zero, the precursors of O₃ (i.e., NO_x and VOCs) are removed from the domain inlet (background boundary conditions) and pollutant source emissions, respectively, and then we compare the results with the base case. In Case_Uref50%, the U_ref is reduced by 50 % to investigate the contribution of wind conditions to the chemical reaction. However, in Case_simple_mech, only three photochemical reactions (Leighton, 1961) are considered in the simulation.

$$\begin{array}{c}{\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+hv\to \mathrm{NO}+\mathrm{O}{(}^{\mathrm{3}}\mathrm{P})\\ {\displaystyle }\mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right)+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{O}}_{\mathrm{3}}+\left[M\right]\\ {\displaystyle }{\mathrm{O}}_{\mathrm{3}}+\mathrm{NO}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\end{array}$$

In order to improve the air quality within the urban area, some cities have tried to implement traffic control policies to reduce the pollutants from vehicle emission sources. Thus, four emission control scenarios are carried out to investigate the effect of emission reduction. Case_Emis_Ctrl50% is the scenario where 50 % of the traffic volume is reduced by applying the odd–even license plate policy (i.e., reducing 50 % of the total vehicle emissions). Case_Emis_Ctrl_VOC20%, Case_Emis_Ctrl_VOC30%, and Case_Emis_Ctrl_VOC40% are the scenarios that apply the stricter VOC control measures (corresponding to 20 %, 30 %, and 40 % more VOC emission reduction, which is a 60 %, 65 %, and 70 % reduction of total VOC emission, respectively) on the vehicles using traffic control policies.

Additionally, the change rate is used to reveal the effect of different factors on pollutant concentrations in the street canyon. For each pollutant, the change rate (CR_p) for different cases is defined as

$$\begin{array}{}\text{(26)}& {\text{CR}}_{\mathrm{p}}\left(\mathit{\%}\right)={\displaystyle \frac{C_{\mathrm{case}}-C_{\mathrm{base}}}{C_{\mathrm{base}}}}\times \mathrm{100}\mathit{\%},\end{array}$$

where C_case and C_base are the concentrations regarded as the condition change case and base case, respectively.

4.2 The comparison of pollutant distribution among the 3D CFD solvers

To investigate the difference between the APonlyChemReactingFoam, APreactingFoam, and APSteadyReactingFoam results, comparisons of O₃, NO, NO₂, and CO distribution are conducted in a $H/W=\mathrm{1}$ street canyon in this study. For APonlyChemReactingFoam, the flow field is treated as the incompressible steady-state flow and pre-solved using the SIMPLE method. The under-relaxation factors and residual threshold for convergence are the same as the setting in Sect. 3.2. Chemical reaction and pollutant dispersion are solved under the steady-state flow for 90 min. For the APreactingFoam and APSteadyReactingFoam cases, turbulence flow, chemical reaction, and pollutant dispersion are solved simultaneously for 90 min. The results in Fig. 13 and all subsequent figures are the pollutant dispersion at 90 min.

Figure 13 The comparison of (a–c) wind speed, (d–f) O₃, (g–i) NO, (j–l) NO₂, and (m–o) CO between APonlyChemReactingFoam, APreactingFoam, and APSteadyReactingFoam.

As depicted in Fig. 13, the wind speed in the APonlyChemReactingFoam case (Fig. 13a) is lower than that in the APreactingFoam (Fig. 13b) and the APSteadyReactingFoam (Fig. 13c) cases. The reason for the difference is most likely due to the different turbulence flow algorithm, where the turbulence is treated as incompressible steady flow, compressible unsteady flow, and compressible steady flow in APonlyChemReactingFoam, APreactingFoam, and APSteadyReactingFoam, respectively. Because of the slight difference in wind speed, the concentrations of APonlyChemReactingFoam (Fig. 13d, g, j, m) for pollutants are higher (due to the lower wind speed) than those in the APreactingFoam (Fig. 13e, h, k, n) and APSteadyReactingFoam (Fig. 13f, i, l, o) cases.

Table 6 The elapsed time of the three solvers.

Table 6 shows the elapsed time of these three simulations in same $H/W=\mathrm{1}$ street canyon for the 90 min simulation. In total, the elapsed time of the APonlyChemReactingFoam case (226 min) is slightly longer than that of the APreactingFoam (214 min) and APSteadyReactingFoam (217 min) cases when employing 192 CPU cores (16× Intel^® Xeon^® E5-2692) for the simulation. However, if the flow field has been determined and there is no need to recalculate in the simulation case, the APonlyChemReactingFoam only takes 191 min to solve the chemical reaction and pollutant dispersion, which is 11 % less time than APreactingFoam.

Many previous studies have treated the urban air turbulence as incompressible steady-state flow and investigate the pollutant dispersion successfully (He et al., 2017; Ng and Chau, 2014; K. Zhang et al., 2019, 2020; Y. Zhang et al., 2019). The APonlyChemReactingFoam is applied in the study over a shorter period of time to analyze the photochemical reaction process in the street canyon.

4.3 Pollutant concentration distribution with a full chemistry mechanism vs. simple chemistry

As shown in Fig. 13a, one main clockwise vortex is formed in the street canyon with $H/W=\mathrm{1}$. The wind speed (WS) is small near the vortex center; i.e., the minimum wind speed is approximated at 0.03 m s⁻¹, which is only 1 % of the speed at the domain inlet. The distributions of pollutants in the street canyon, such as O₃, NO, and NO₂, are also swirling (Fig. 13a, d, g, j).

Figure 14 Pollutant distribution of (a) $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$, (b) ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$, and (c) $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ in the base case. (d) Schematic diagram of the formation mechanism of photochemical reactions (Tang et al., 2006).

Leeward-side O₃ concentration in the base case is less than the windward side, while NO and NO₂ are the opposite. At the corner of leeward side, the minimum value of O₃ and maximum value of NO_x appear, with less than 20 ppbv for O₃, more than 200 ppbv for NO and 140 ppbv for NO₂ (Fig. 13a, d, g, j), respectively. Meanwhile, due to the higher NO emissions, the ratios of NO and NO₂ are higher at the bottom of the street canyon (Fig. 14a). The larger $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ values indicate that the titration effect from NO (${\mathrm{O}}_{\mathrm{3}}+\mathrm{NO}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}$) and ozone depletion would be stronger, leading to the lower O₃ concentrations in this area.

However, on the windward side, NO_x concentrations are less than that on the leeward side. This is because the NO_x from the emission source first affects the leeward side, which leads to the high concentrations in this area. As the wind flows, the concentrations of NO_x gradually decrease due to the wind diffusion and the dilution effect. With the comparison of background, the windward NO and NO₂ concentrations increase by approximately 35 and 55 ppbv, respectively. On the one hand, pollutants from emissions are transported along the flow, which increases the concentrations. On the other hand, VOCs in street canyons react with OH via the chemical reactions and generate HO₂ and RO₂. Theis RO₂, HO₂, and O₃ would react with NO and generate NO₂, leading a higher increment for NO₂ (Fig. 14b–c). It should be noted that the O₃ concentrations could be higher due to the lower depletion reaction with NO compared to that of the leeward side.

Figure 15 Change rate of (a) O₃, (b) NO, (c) NO₂, and (d) $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ between the simple chemistry and full chemistry mechanisms.

Figure 15 shows the change rate of pollutant concentrations (Fig. 15a–c) and the NO to NO₂ ratio change rate of Case_simple_mech compared with the base case (Fig. 15d). Without the consideration of the VOC-related reactions in the mechanism, there is no consumption of NO by RO₂ in the mechanism ($\mathrm{VOCs}+\mathrm{OH}\to {\mathrm{RO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$, ${\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}\to \mathrm{RO}+{\mathrm{NO}}_{\mathrm{2}}$), and the NO titration effect (${\mathrm{O}}_{\mathrm{3}}+\mathrm{NO}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}$) would be stronger in this case. For O₃, the concentrations are 36 %–58 % lower than that in base case within the street canyon; NO₂ concentrations are also 15 %–40 % lower, and NO could be up to 90 % higher than that in the simple chemistry case. Thus, the NO to NO₂ ratio would be 60 %–150 % higher in the simple chemistry case.

Figure 16 Change rate of (a) O₃, (b) NO, (c) NO₂, (d) $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$, (e) ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$, and (f) $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ at t=5400 s .(g) Time series of the reaction rate of RO₂ ($\frac{\mathrm{d}{\mathrm{RO}}_{\mathrm{2}}}{\mathrm{d}t}$) of Case_BC_zero.

4.4 Influence of background precursors of O₃ on reactive pollutant concentrations

In Case_BC_zero, all background precursors of O₃ (i.e., NO_x and VOCs) from the upstream domain inlet are removed. As depicted in Fig. 16a, the change rates of O₃ are negative, confirming that O₃ concentration becomes lower without the background NO_x and VOCs. In addition, the O₃ reduction rate on the windward side (−5 % to −8 %) is smaller than that on the leeward side (−9 %). The influencing mechanisms of O₃ reduction are complicated and will be explained later.

On the one hand, by analyzing Fig. 16b and c, such reduction rates on the windward side for NO (−12 % to −16 %) and NO₂ (−20 % to −24 %) are greater than those on the leeward side (−2 % to −6 % for NO and −12 % to −15 % for NO₂). Therefore, the ratio of $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ increases from about 9 % to 14 % in the street canyon (Fig. 16d). Overall, this increment of $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ enhances O₃ depletion because the main source of ozone is the photolysis reaction of NO₂. Meanwhile, the main sink is the titration effect of O₃ and NO.

On the other hand, the RO₂ from the oxidation of VOCs with OH will consume the NO ($\mathrm{VOCs}+\mathrm{OH}\to {\mathrm{RO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$, ${\mathrm{RO}}_{\mathrm{2}}+\mathrm{NO}\to \mathrm{RO}+{\mathrm{NO}}_{\mathrm{2}}$), which would affect the NO_x–O₃ circulation. Figure 16e shows that the reduction of ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$ on the windward side is more than that on the leeward side, which indicates that the background VOCs and OH reaction on the windward side are more active. However, due to the slower reaction rate of RO₂ and NO compared to that of HO₂ and O₃ with NO, the conversion of RO₂ to RO by reacting with NO would require more time. The reduction rate of $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ (Fig. 16f) on the leeward side (−18 % to −23 %) is slightly greater than that on the windward side (−16 % to −17 %). Therefore, the influence of RO₂ on NO_x–O₃ circulation would gradually appear on the leeward side, in addition to the flow transportation.

Additionally, Fig. 16g shows the reaction rate of RO₂ ($\frac{\mathrm{d}{\mathrm{RO}}_{\mathrm{2}}}{\mathrm{d}t}$) at the bottom, center, and top point on the centerline of the street canyon with (base) and without (Case_BC_zero) background conditions. At the bottom point (CB), the reduction rate of RO₂ is lower in the base case. This is because the background VOCs and OH reaction consume a portion of NO on the windward side, which leads to a lower consumption of RO₂ with NO. As the simulation continues, NO concentrations would increase due to the continuous release of large amounts of NO from source emissions, and the reduction rate of RO₂ becomes lower in Case_BC_zero due to the lack of background VOCs.

At the top (ST) and center point (SC), however, as the reaction goes on and pollutants mix upwards, the NO concentration could become higher. Therefore, the RO₂ consumption rate is lower without the background RO₂. This reduction indicates that the background conditions mainly consume NO in the street canyon, which leads to the increase of O₃ due to the weakening titration effect.

Figure 17 Change rate of (a) O₃, (b) NO, (c) NO₂, (d) $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$, (e) ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$, and (f) $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ at t=5400 s. (g) Time series of the reaction rate of RO₂ ($\frac{\mathrm{d}{\mathrm{RO}}_{\mathrm{2}}}{\mathrm{d}t}$) of Case_Emis_zero.

4.5 Effects of vehicular source emissions on reactive pollutants

In Case_Emis_zero, the precursors of O₃ (NO_x and VOCs) from near-ground emissions are removed. As shown in Fig. 17a, O₃ concentrations increase by over 30 %–120 % in the whole street canyon compared to the base case. In particular, the O₃ increment on the leeward side is from 80 % to 250 % (not shown here), which is much higher than that on the windward side (30 % to 40 %). However, NO_x concentrations decrease significantly, i.e., the reduction rates are −84 % to −98 % for NO and −76 % to −90 % for NO₂, showing that the NO_x from source emissions is the dominant part of NO_x in the street canyon (Fig. 17b and c). The large reduction of NO_x concentrations induces the increase of O₃ concentrations with weaker titration effect of O₃. Specifically, both the maximum increase for O₃ and minimum reduction of NO and NO₂ appear at the near-ground corner of the leeward side, which is the downwind area of the pollutant source. In addition, due to a larger amount of NO emissions than NO₂ (emission ratio of NO to NO₂ is 9:1), the concentration ratio of $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ considerably decreases with the reduction rates of −30 % to −70 % (Fig. 17d) if vehicular pollutant sources are removed.

Additionally, due to the large reduction of NO and NO₂ concentrations, more OH would react with VOCs instead of NO_x, which increases the RO₂ concentration (Fig. 17e and f). Meanwhile, with the reduction of NO concentration, the consumption of RO₂ significantly decreases, which leads to the dramatic increase of RO₂ concentration. Thus, the ratio of ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$ rises by 115 %–205 % and the ratio of $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ decreases by −60 % to −88 %.

In Fig. 17g, the reaction rate of RO₂ at three points in Case_Emis_zero is positive, which means that the RO₂ keeps being generated but is not consumed among these three points. As mentioned above, RO₂ (the production of VOCs and OH) will consume the NO and weaken the O₃ titration effect with NO. In the base case (Fig. 17g), the reaction rate of RO₂ is negative, which means that RO₂ consumes the NO. However, in Case_Emis_zero, the reaction rate of RO₂ is positive during the whole simulation period, which means that there is not enough NO to react with RO₂ or even O₃ without the vehicular source. Therefore, the source emissions provide a large amount of NO, which enhances the O₃ depletion in the street canyon.

Figure 18 Change rate of (a) O₃, (b) NO, (c) NO₂, (d) $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$, (e) ${\mathrm{RO}}_{\mathrm{2}}/\mathrm{OH}$, and (f) $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ at t=5400 s.(g) Time series of the reaction rate of RO₂ ($\frac{\mathrm{d}{\mathrm{RO}}_{\mathrm{2}}}{\mathrm{d}t}$) of Case_Uref50%.

4.6 Influence of wind velocity reduction on reactive pollutants

Figure 18 shows the change rates of O₃, NO_x, and ratios when the background wind speed decreases from U_ref=3 m s⁻¹ (base) to U_ref=1.5 m s⁻¹ (Case_Uref50%). In Fig. 18a, there is no significant change of O₃ concentration at the center of street canyon. However, in the downwind area of the near-ground pollutant source, O₃ concentration decreases by 5 % to 30 % compared with that of the base case. Interestingly, at the bottom of the leeward side, O₃ has an increase up to 6 % under the half inlet wind speed condition.

Due to the weaker capacity of pollutant dilution caused by the smaller wind speed, the concentrations of NO and NO₂ almost double (i.e., rising by 80 %–98 % in Fig. 18b and c), but the $\mathrm{NO}/{\mathrm{NO}}_{\mathrm{2}}$ change rate has no significant change (−1 % to −3 % in Fig. 18d). Besides, a higher NO_x concentration would react with more OH, which consequently weakens the RO₂ production from VOCs (−8 % to −20 % in Fig. 18e). Meanwhile, the increase of NO concentration consumes more RO₂ to RO, which leads to an 180 % to 340 % increase of $\mathrm{RO}/{\mathrm{RO}}_{\mathrm{2}}$ (Fig. 18f).

Additionally, Fig. 18g illustrates the RO₂ reduction in a street canyon. Because of the higher concentration of NO in Case_Uref50%, the RO₂ reduction rates at three monitoring points are higher than that in the base case, particularly at the bottom of the street canyon (CB). In the early stages of the reaction, the reduction rate of RO₂ at the top point (ST) is slightly lower in Case_Uref50%. This is because the RO₂ concentration at ST is first affected by the background. As the NO concentrations increase in the whole street canyon, the RO₂ consumptions become higher than that in the base case.

Figure 19 The (a) O₃, (b) NO, and (c) NO₂ concentrations at 90 min in different emission control scenarios.

4.7 Emission control strategy on reactive pollutant concentrations

Figure 19 shows the concentrations of O₃, NO, and NO₂ in different NO_x and VOC emission control scenarios at 90 min. In Case_Emis_ctrl50% (the emission of NO_x and VOCs reduces 50 %, i.e., 50 % reduction of traffic volume), the O₃ concentration increases from 19–47 to 29–54 ppbv (Fig. 19a). On the contrary, this control measure for NO and NO₂ is very effective (Fig. 19b and c), and NO and NO₂ concentrations reduce from 47 % to 54 % and 37 % to 40 % in Case_Emis_ctrl50%, respectively.

This indicates that the simple traffic control measures cannot effectively reduce O₃ concentration. This is because most of the urban areas are in VOC-sensitive regions (Ye et al., 2016). When the total number of vehicles decreases under the traffic control measures, the reduction of NO_x is higher than that of the VOCs (due to the larger NO_x emission from vehicles), which leads to a higher VOC-to-NO_x ratio, consequently resulting in a higher O₃ concentration in the street canyon (Sillman and He, 2002). Thus, in order to reduce the concentrations of O₃, the stricter VOC control measures on vehicles should be conducted. Based on the results shown from three other emission control scenarios (Case_Emis_ctrl_VOCs20%, Case_Emis_ctrl_VOCs30% and Case_Emis_ctrl_VOCs40%) in Fig. 19a, the emission of VOCs needs be reduced by another 30 % under traffic control (Case_Emis_ctrl50%) to bring the O₃ concentrations back to the level that they were at when no traffic control measures had been taken (base case).

As for NO and NO₂, when the additional VOC control measures are carried out, the concentrations are higher than those in Case_Emis_ctrl50%. Even so, their concentrations do not still exceed the concentration level before the traffic control (base case), which means that such an emission control scenario is still effective for NO and NO₂. In summary, the control policies of reactive pollutants require the comprehensive consideration of the relationship between precursors and pollutants so that the goal of improving air quality can be achieved.

5 Conclusions

A detailed description of the atmospheric photolysis calculation framework APFoam 1.0 is presented in this paper, and this CFD model is coupled with multiple full atmospheric photochemical mechanisms, including SAPRC07 (CS07A and SAPRC07TB) and CB05. In order to simulate the photochemical process of reactive pollutants, five new types of the reactions, i.e., the new form of the (1) Arrhenius reactions, (2) photolysis reactions, (3) falloff reactions, (4) “three-k” reactions, and (5) “two-k” reactions, have been modified and added into APFoam. Additionally, to verify the model performance, several validations, including a photochemical mechanism (CS07A) with SAPRC box modeling, a flow field, and 2D and 3D pollutant dispersions with wind tunnel experimental data have been conducted in this study. The model results show a good agreement with the SAPRC box modeling and wind tunnel experimental data, indicating that APFoam can be applied in the analysis of microscale urban pollutant dispersion.

Key factors of chemical processes are investigated by applying APFoam with a CS07A mechanism in the simulation of reactive pollutants in a typical street canyon ($H/W=\mathrm{1}$) with a VOC to NO_x emission ratio of ∼5.7 ppbv s⁻¹. In the comparison of chemical mechanisms, O₃ and NO₂ are underestimated by 36 %–58 % and 15 %–40 %, respectively, while NO is overestimated by 30 %–90 % without the consideration of the VOC reactions. Other numerical sensitivity cases (Case_BC_zero, Case_Emis_zero, and Case_Uref50%) reveal that vehicle emissions are the main source of NO and NO₂, with a contribution of 82 %–98 % and 75 %–90 %, respectively. The resident part of the NO_x in the street canyon is contributed by the background concentration. However, vehicle emissions with a large amount of emitted NO_x, especially NO, are the main reason for the decrease of O₃ due to the stronger NO titration effect within the street canyon. In contrast, 5 %–9 % of the O₃ is contributed by the boundary conditions. Ventilation conditions are another reason for the NO_x concentration increments, and the increase of NO_x can be up to 98 % when the wind speed is reduced by half. If there are no chemical reactions, NO_x concentration should rise by 100 % when the wind velocity decreases by 50 % (i.e., ventilation capacity reduces by 50 %) because the Re independence requirement is satisfied. However, O₃ is reduced downwind of the emissions due to the increase of NO concentrations. In order to control and improve the air quality in the street canyon, traffic control policies are effective for NO_x. However, our results indicate that at least another 30 % reduction in vehicle VOC emissions could reduce O₃ concentrations below those of the odd–even license plate policy, with 24 %–32 %, 25 %–28 %, and −6 %–2 % reduction rates of NO, NO₂, and O₃, respectively. Overall, APFoam 1.0, a fully coupled CFD model, can be employed to investigate atmospheric photolysis calculation in urban areas and provide reliable and useful suggestions for the improvement of urban air quality.

6 Future plans

However, in the current version of APFoam, aerosol chemistry is not included in the model, and thus it is necessary to couple it with aerosol processes, such as MOSAIC (Zaveri et al., 2008) or ISORROPIA (Fountoukis and Nenes, 2007; Nenes et al., 1998), in future work. In addition, the photolysis rates in the current model have been fixed without diurnal variation, which means that the model is not suitable for a long-term simulation and needs to be updated in subsequent versions. Moreover, the interaction between radiation and chemical reactions will be investigated by APFoam and validated by the scaled outdoor experiment (Chen et al., 2020a, b) in the future.