Error induced by neglecting subgrid chemical segregation due to inefficient turbulent mixing in regional chemical-transport models in urban environments

We employed direct numerical simulations to estimate the error on chemical calculation in simulations with regional chemical-transport models induced by neglecting subgrid chemical segregation due to inefficient turbulent mixing in an urban boundary layer with strong and heterogeneously distributed surface emissions. In simulations of initially segregated reactive species with an entrainment-emission configuration with an A–B–C second-order chemical scheme, urban surface emission fluxes of the homogeneously emitted tracer A result in a very large segregation between the tracers and hence a very large overestimation of the effective chemical reaction rate in a complete-mixing model. This large effect can be indicated by a large Damköhler number (Da) of the limiting reactant. With heterogeneous surface emissions of the two reactants, the resultant normalised boundary-layer-averaged effective chemical reaction rate is found to be in a Gaussian function of Da, and it is increasingly overestimated by the imposed rate with an increased horizontal scale of emission heterogeneity. Coarse-grid models with resolutions commensurable to regional models give reduced yet still significant errors for all simulations with homogeneous emissions. Such model improvement is more sensitive to the increased vertical resolution. However, such improvement cannot be seen for simulations with heterogeneous emissions when the horizontal resolution of the model cannot resolve emission heterogeneity. This work highlights particular conditions in which the ability to resolve chemical segregation is especially important when modelling urban environments.

1 Introduction

Turbulence mixes initially segregated reactive species in the boundary layer and allows chemical reactions to occur. However, for fast chemical reactions with the chemical timescale shorter than the turbulent timescale, turbulent motions mix the reactants so slowly that they remain segregated rather than reacting. This segregation can be a result of the inefficient mixing due to the state of turbulence and its driver, such as thermal stability, canopy–atmosphere interaction and cloud processes, and/or a result of the heterogeneity of surface emissions (Vilà-Guerau de Arellano, 2003). In an Eulerian chemical-transport model, these processes and characteristics often occur at a length scale smaller than the size of a model grid box that the resultant chemical segregation cannot be represented, as the chemical species are assumed to be instantaneously and homogeneously mixed within a grid box. The negligence of such subgrid segregation induces potential errors to the calculation of the chemical transformation within a model grid of a large-scale model.

Efforts have been made to examine and quantify this error under different turbulent and chemical regimes in a range of atmospheric environments. Earliest studies can be dated back to Damköhler (1940) and Danckwerts (1952), which focused on quantifying the effect of turbulence on combustion processes, also applicable to other chemical transformations such as chemical reactions, by introducing quantitative definitions of scaling parameters including the Damköhler number (Da, the ratio between the turbulent and chemical timescales) and the segregation coefficient (I_S, the ratio between the covariance and the product of the mean concentrations). Donaldson and Hilst (1972) adopted the discussion in the context of atmospheric reactions. This line of research was then continued by a number of investigations that used first-order and second-order closure methods to study the profiles and budgets of the fluxes of chemical reactive species in the boundary layer. For example, Fitzjarrald and Lenschow (1983) and Vilà-Guerau de Arellano and Duynkerke (1993) reported that the chemical terms can be of similar importance as the dynamical term in the flux equations for the $\mathrm{NO}-{\mathrm{NO}}_{\mathrm{2}}-{\mathrm{O}}_{\mathrm{3}}$ chemical system in the surface layer. Verver et al. (1997) further stated that ignoring the higher-order chemical terms in the flux equations would result in deterioration of model performance. As these closure methods failed to resolve vertical turbulent mixing and horizontal fluctuations, the investigation of the topic was extended with the use of large-eddy simulations (LESs), which can resolve the most energetic eddies in the boundary layer.

Many of these LES studies focus on the convective boundary layer (CBL), in which the imbalance between updraught and downdraught transport produces a large segregation of the reactants (Wyngaard and Brost, 1984; Chatfield and Brost, 1987). Such LES studies were often performed with idealised cases with a bottom-up tracer emitted from the surface and top-down tracer entrained from the free troposphere with a simple second-order chemistry scheme. For instance, Schumann (1989) pointed out that the segregation between the two tracers depends on the Damköhler number (Da), the concentration ratio of the two species and the specified initial conditions, pinning the use of Da as an indicator to estimate whether turbulent motions are significantly affecting chemical reactions in the flow. Vinuesa and Vilà-Guerau de Arellano (2005) introduced the concept of an effective chemical reaction rate (k_eff) to quantify the actual boundary-layer-averaged reaction rate that accounts for the effect of the chemical segregation. They also reported a drop of k_eff up to 20 % from the imposed chemical reaction rate k when Da is ∼1, while k_eff≈k when Da is ∼0.1. Recently, the effect of segregation due to inefficient turbulent mixing on chemical reaction has been considered the cause of the miscalculation of large-scale models in a number of scenarios. One of the most discussed issues is the underprediction of the concentration of hydroxyl radical (OH) in global models. A number of studies employed LES models over forestal areas with sophisticated chemical mechanisms involving biogenic volatile organic compounds (VOCs) to simulate the resultant segregation between isoprene and OH (e.g. Patton et al., 2001; Brosse et al., 2017; Dlugi et al., 2019). However, the magnitude of the segregation coefficient I_S shown in these studies is in general less than 20 %, which is too small to explain the observed discrepancy, where I_S of $\sim -\mathrm{50}\mathrm{\%}$ is necessary (Butler et al., 2008).

There are a number of factors that may increase chemical segregation. One possible factor is to increase the surface emission of the bottom-up tracer. For instance, Kim et al. (2016) indicated the segregation between isoprene and OH increases when switching from very low NO_X to very high NO_X surface emissions. Another possible factor is related to the emission heterogeneity of the bottom-up tracer. Krol et al. (2000) found that, by imposing non-homogeneous emission of a generic hydrocarbon (RH) at the surface, the segregation coefficient between RH and OH can reach up to −30 %, compared to only −3 % when RH is emitted homogeneously. Ouwersloot et al. (2011) implemented in their LES simulations patches of surface isoprene emissions with different fluxes on the African savannah, and they showed that the segregation increases with the width of the patches. Kaser et al. (2015) further pointed out that with the consideration of surface heterogeneity, chemical segregation can locally slow down the isoprene chemistry by up to 30 %. Most of the previous LES studies focus on scenarios under rural conditions, but these two factors should be even more relevant for urban environments with intense emissions of pollutants and high spatial heterogeneity.

The aim of the present work is to investigate the effect of inefficient turbulent mixing on chemical reactions in an urban-like boundary layer with strong and heterogeneously distributed surface emissions and to account for the errors induced by neglecting the resultant subgrid chemical segregation in relatively coarse regional models. While previous studies focused on agricultural and rural conditions where the emission fluxes are relatively low (∼O (0.01) ppb m s⁻¹), our simulations address cases of strong emission fluxes in typical urban values (∼O (0.1–1.0) ppb m s⁻¹). In two rare related studies in urban air conditions, Baker et al. (2004) reported from their LES model of an urban street canyon significant deviations of the concentrations of the $\mathrm{NO}-{\mathrm{NO}}_{\mathrm{2}}-{\mathrm{O}}_{\mathrm{3}}$ triad from the equilibrated values to the photostationary state depending on the turbulent structure in the canyon, while Auger and Legras (2007) obtained high values of instantaneous segregation under certain emission configurations in urban areas with a chemical system of 44 species. On top of their conclusions, our study aims to explain why this strong segregation occurs under urban conditions and on which parameters the errors induced by neglecting the segregation in large-scale models depend. To achieve this aim, we perform direct numerical simulation (DNS) with homogeneous emissions with an idealised second-order A–B–C chemistry scheme ($\mathrm{A}+\mathrm{B}\to \mathrm{C}$) with emission fluxes extended to urban values, in addition to a set of simulations with heterogeneous emissions. With varied reaction rates, the idealised second-order A–B–C chemistry scheme can generally represent any second-order chemical reactions commonly seen in an urban environment.

The results from our DNS runs are then degraded to lower resolution to mimic the calculations from regional models. Previous studies often compared their LES results with a mixed-layer or complete-mixing model, which assumes the whole simulation domain to be in the same model grid. This assumption may still be reasonable for forested areas, as the global and mesoscale models in use typically have a mesh size of the order of 10 km, which is comparable to the size of our simulation domain. However, this comparison is not suitable in applications to urban environments, as the regional models simulating urban air quality typically have a much finer resolution with mesh sizes of the order of 1 km. Therefore, we degrade our DNS results to coarse-grid models with mesh sizes of 1 and 3 km with multiple vertical resolutions to estimate the errors from regional models.

The structure of this paper is as follows. The next section introduces the DNS model adopted in this work and the settings of the simulations. Subsequently, the results of our DNS runs are first presented with cases of homogeneous emissions and then of heterogeneous emissions. This is followed by the comparison between the results from the degraded coarse-grid models and those from the DNS model to account for the errors from regional models. The implication of our results for regional models applied to urban environments is then discussed, and conclusions are provided at the end.

2 Model description

2.1 Dynamical settings

The relation between turbulence and chemistry in a convective boundary layer is investigated in this work by means of direct numerical simulation. There are two common approaches to numerically simulate turbulent flows, namely large-eddy simulation (LES) and direct numerical simulation (DNS). The former applies a low-pass filter and models the subgrid-scale effect on the filtered variables. The later solves the original Navier–Stokes equations but often with a molecular viscosity that is larger than the value in the real application represented (Monin and Mahesh, 1998; Pope, 2004; Mellado et al., 2018). Hence, both approaches are restricted to low-to-moderate Reynolds numbers. The Reynolds number is defined in terms of the subgrid-scale viscosity and numerical diffusivity in LES and in terms of the molecular viscosity in DNS. It can be interpreted as the scale separation between the largest and the smallest resolved scales (i.e. the Kolmogorov and Batchelor scales in DNS), which is in turn related to the number of grid points in the simulation and hence the computational resources needed for the simulation. Current computational resources allow the Reynolds numbers in simulations to reach the order of 10³–10⁴ (the range within which the Reynolds number of our simulations lies; see Appendix A), which are substantially smaller than the typical atmospheric values of 10⁸. Fortunately, relevant turbulence statistics such as variances, covariances and entrainment rates tend towards Reynolds-number independence once the Reynolds number increases beyond 10³–10⁴ (Dimotakis, 2000; Mellado et al., 2018). When studying the properties of the smallest resolved scales and their substantial effect on large-scale dynamics, DNS eliminates the uncertainty introduced by subgrid-scale models and numerical artefacts, which complicates the interpretation and the comparison of results from different LES models (Mellado et al., 2018). DNS provides results that depend only on the achievable Reynolds number in the simulation, and one can perform sensitivity analyses to estimate such dependence. The disadvantage of DNS is the higher computational cost, since DNS typically utilises higher-order numerical schemes and higher resolutions. Nonetheless, DNS is gaining traction as it approaches a degree of Reynolds number similarity that allows for certain extrapolation to atmospheric conditions (Monin and Mahesh, 1998; Jonker et al., 2012, 2013; Waggy et al., 2013; Mellado et al., 2018; van Hooft et al., 2018; Mellado, 2019).

In this work, we opt for DNS because we are studying the variance and covariances of various fields in which entrainment may play an important role. It is known that the smallest resolved scales might become important for these variables in typical resolutions (Sullivan and Patton, 2011). In addition, with strong emission fluxes, subgrid turbulent parameterisation in LES simulations may induce significant errors near the surface where the tracer is emitted (Vinuesa and Porté-Agel, 2005, 2008). Therefore, here DNS provides an alternative for studying the topic and also a tool for an intercomparison study of the LES results. In Sect. 3.1.1 and 3.1.2, we will intercompare the results from our DNS model with the LES results by adopting the same initial conditions as in Vinuesa and Vilà-Guerau de Arellano (2005) for some of our simulations.

We employed the direct numerical simulation tool Turbulence Laboratory (TLab) to perform our computational experiments of turbulent mixing of reactive species in the convective boundary layer. Source files with the implementation of TLab and further documentation can be found at https://github.com/turbulencia/tlab (last access: 13 January 2021). The dynamical part of the simulations performed in this work follows similar settings as in Garcia and Mellado (2014) and Van Heerwaarden et al. (2014). In the model, the Navier–Stokes equations of incompressible fluids in the Boussinesq approximation are solved to calculate the buoyancy and the velocity fields:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\\ {\displaystyle }& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathrm{\nabla }\cdot (\mathit{u}\otimes \mathit{u})=-\mathrm{\nabla }p+\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{u}+b\mathit{k},\\ {\displaystyle }& {\displaystyle \frac{\partial b}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{u}b\right)=\mathit{\kappa }{\mathrm{\nabla }}^{\mathrm{2}}b,\end{array}$$

where u(x,t) is the velocity vector with components ($u,v,w$) in the directions $\mathit{x}=(x,y,z)$ at time t, respectively. p is the modified pressure divided by the constant reference density, and b(x,t) is the buoyancy, expressed as $b=g({\mathit{\theta }}_v-{\mathit{\theta }}_{v,\mathrm{0}})/{\mathit{\theta }}_v$. The background buoyancy profile is set as b₀(z)=N²z, where N² is the Brunt–Väisälä frequency. The surface buoyancy flux is set to be F_b (see Appendix A). The parameters ν and κ are the kinematic viscosity and the molecular diffusivity, respectively. k is the unit vector in the z direction. The system is statistically homogeneous in the horizontal direction, such that its statistics depend on z and t only.

A no-penetration, no-slip boundary condition is imposed at the surface, and a no-penetration, free-slip boundary condition is imposed at the top boundary. Neumann boundary conditions are imposed for the buoyancy and velocity fields at both the top and the surface to maintain constant fluxes. The velocity and buoyancy fields are relaxed towards zero and N²z, respectively, at the top of the computational domain. The height of the top boundary is adjusted so that the turbulent region is far enough from the top to avoid significant interaction. Periodic boundary conditions are implemented in lateral directions.

The size of the computational grid is $\mathrm{720}\times \mathrm{720}\times \mathrm{512}$ for all simulations (the number of vertical layers is 512). Stretching is applied vertically to increase the resolution near the surface in order to resolve the surface layer and to stretch the grid size in the upper portion of the domain end to further separate the top boundary from the turbulent region. The total simulation domain size is $\mathrm{120}{L}_{\mathrm{0}}\times \mathrm{120}{L}_{\mathrm{0}}\times \mathrm{34.4}{L}_{\mathrm{0}}$, where L₀ is the reference length scale. Together with the reference time and velocity timescales T₀ and U₀, L₀ is used to non-dimensionalise the Navier–Stokes equations (see Appendix A). With typical atmospheric parameters of L₀∼100 m and U₀∼1 m s⁻¹, the horizontal resolution of the model in use is equivalent to 15 m×15 m with a total domain size of 12 km×12 km. The resolution is considerably higher than those adopted in the previous LES studies, and it is fine enough to resolve convective turbulent motion in the convective boundary layer. In the vertical direction, the grid spacing increases from Δz=6.72 m at the surface to Δz=266.28 m at the top of the domain at 16 km. The vertical grid spacing vs. height in the first 5 km is shown in Fig. 7 (black line).

Figure 1 Schematic diagrams of the configurations of the DNS runs with homogeneous emissions (a) and heterogeneous emissions (b).

Table 1 The names and simulation parameters of the DNS runs, including (a) the imposed chemical reaction rate k and the surface emission flux of tracer A, F_A, for the simulations with homogeneous emissions and (b) plus the length of heterogeneity (dx) for the simulations with heterogeneous emissions.

We use a sixth-order compact scheme to calculate the spatial derivatives and a fourth-order Runge–Kutta scheme to advance the equations in time. The pressure Poisson equation is solved by applying a Fourier decomposition along the horizontal directions and solving the resultant set of finite-difference equations in the vertical direction to machine accuracy (Mellado and Ansorge, 2012). The sensitivity to grid resolution has been tested up to the triple-velocity correlation in the transport term of the evolution equation for the turbulence kinetic energy, observing an error of less than 5 % at the surface when doubling the spatial resolution. Further details of the numerical algorithm and the grid-sensitivity studies can be found in Mellado (2010, 2012), Garcia and Mellado (2014), and Mellado et al. (2017).

The simulations are terminated after a total simulation time equivalent to 4.5 h, during which the boundary layer grows convectively. The simulation time represents the hours from sunrise to midday, when the boundary layer is well developed and statistical equilibrium is attained. The boundary layer height z_i is identified as the height where the buoyancy variance $\left[b^{\prime }{b}^{\prime }\right]$ is maximum away from the surface layer (Garcia and Mellado, 2014). At the end of our simulations, the boundary layer height is around 2.2 km. As the boundary layer grows in depth, the corresponding convective velocity and timescales evolve in time. In this study, the convective timescale is adopted as the turbulent timescale t_turb, which is related to the convective velocity w_c by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{t}_{\text{turb}}={\displaystyle \frac{z_i}{w_{\text{c}}}}={\displaystyle \frac{{z_i}^{\frac{\mathrm{2}}{\mathrm{3}}}}{F_{\mathrm{b}}}}\end{array}$$

(Deardorff, 1970). This turbulent timescale will be used to calculate the Damköhler number Da in later sections.

2.2 Chemical settings

2.2.1 Homogeneous emissions

An archetypical entrainment-emission configuration, as typically used by past LES studies such as Krol et al. (2000) and Vinuesa and Vilà-Guerau de Arellano (2005), is adopted in our DNS model (see left panel of Fig. 1). The simulations are conducted with two theoretical tracers (A and B). Tracer A is emitted from the surface at a constant flux F_A at a range of values. Tracer B is entrained from the free troposphere, in which the mixing ratio of tracer B (〈B〉₀) is constantly and homogeneously fixed at 2 ppb. Tracer C is formed in a second-order chemical reaction between tracers A and B:

$$\begin{array}{}\text{(R1)}& k:\mathrm{A}+\mathrm{B}\to \text{C},\end{array}$$

where k is the imposed chemical reaction rate. The mixing ratios of the three species, denoted as A, B and C, respectively, vary with time with rates equal to

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=kAB.\end{array}$$

Neumann boundary conditions are imposed at the boundaries of the computational domain on the mixing ratios of the tracers so that the surface flux of tracer A is constant at F_A, and the surface fluxes of B and C are zero. In the free troposphere, the initial mixing ratios of A and C are zero. The initial mixing ratio of tracer B is decreased linearly down the boundary layer to zero on the surface.

Four emission fluxes of tracer A (F_A), at 0.05, 0.25, 0.5 and 1.41 ppb m s⁻¹, are adopted to test the effect of strong emissions in an urban environment. The runs are named VV05 (from its similar initial conditions as Vinuesa and Vilà-Guerau de Arellano, 2005), mflux, sflux and ssflux, respectively. As a reference, in the ssflux run, the characteristic mixing ratio scale for tracer A (〈A〉₀) is equal to that of tracer B (〈B〉₀). These characteristic scales are used in the non-dimensionalisation of Eq. (2) (see Appendix A). Two chemistry cases, namely slow and fast chemistry, are considered for each of the imposed emission fluxes, with respective chemical reaction rates of $\mathrm{4.75}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{4.75}\times {\mathrm{10}}^{-\mathrm{3}}$ ${\mathrm{ppb}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$. In the real atmosphere, the rate of chemical reaction between NO and O₃ is $\mathrm{4.75}\times {\mathrm{10}}^{-\mathrm{4}}$ ${\mathrm{ppb}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, corresponding to the cases with slow chemistry. Meanwhile, the reaction rate between isoprene and OH is 0.1 ${\mathrm{ppb}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ (Karl et al., 2004), giving an example illustrating the cases with fast chemistry. The names and simulation parameters of the eight simulations are listed in Table 1.

2.2.2 Heterogeneous emissions

For the simulations with heterogeneous emissions, tracers A and B are emitted alternately from patches on the surface with widths of 1, 2 and 6 km. The widths of these emission patches are referred to in this study as the length of heterogeneity dx, similarly as in Ouwersloot et al. (2011).¹ The A–B–C chemical system ($\mathrm{A}+\mathrm{B}\to \mathrm{C}$) is implemented with slow and fast chemistry, as described in Sect. 2.2.1. Both tracers A and B are emitted at the same imposed flux, with the two values of emission fluxes at 0.05 ppb m s⁻¹ (VV05) and 0.5 ppb m s⁻¹ (sflux) adopted.² The right panel of Fig. 1 shows a schematic diagram describing the simulation configuration. The names and simulation parameters of the twelves simulations are listed in Table 1. The simulations with heterogeneous emissions are run for 8 h, representing the hours from sunrise to the afternoon. These simulations are run for a longer time than those with homogeneous emissions because the statistical equilibrium takes a longer time to be attained with heterogeneous emissions.

In the context of an urban environment, the second-order chemical reaction imposed in our simulations can be considered analogous to the reaction between NO and peroxyl radical derivatives (RO₂) from the oxidation of VOCs, which is the limiting reaction of ozone production in the VOC-limited regime. The VOCs are often emitted from sources segregated from NO at very high fluxes in urban areas. The involved chemical reactions are often relatively fast. For instance, the reaction rate between NO and acetone peroxy radical (CH₃COCH₂O₂), which is a common VOC species found in urban environments (Brasseur and Jacob, 2017; Li, 2019), is 0.503 ${\mathrm{ppb}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ (Manion et al., 2008), giving an example for the fast-chemistry cases in our simulations with heterogeneous emissions.

Table 2 The initial Damköhler numbers Da_A,i and Da_B,i, and final Damköhler numbers Da_A,f and Da_B,f with the resultant normalised effective chemical reaction rate k_eff∕k in the DNS runs with (a) homogeneous emissions and (b) heterogeneous emissions.

2.3 Quantifying the chemistry–turbulence interaction

Two numbers are mainly employed to quantify the effect of turbulent mixing on chemical reaction, namely the Damköhler number Da and the effective chemical reaction rate k_eff. The first number, the Damköhler number (Da), is defined as the ratio between the turbulent timescale (t_turb) and the chemical timescales (t_chem) (Damköhler, 1940). For a second-order chemical reaction $\mathrm{A}+\mathrm{B}\to \mathrm{C}$, the Damköhler number of tracer A is given by

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }D{a}_{\mathrm{A}}={\displaystyle \frac{t_{\text{turb}}}{t_{\text{chem,A}}}}={\displaystyle \frac{z_i/w_{\text{c}}}{\mathrm{1}/k\langle B\rangle }}=k\langle B\rangle {\left({\displaystyle \frac{z_i^{\mathrm{2}}}{F_b}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

Similarly the Damköhler number of tracer B can be written as

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }D{a}_{\mathrm{B}}=k\langle A\rangle {\left({\displaystyle \frac{z_i^{\mathrm{2}}}{F_b}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

The Damköhler numbers of tracers A and B at the start (Da_A,i and Da_B,i) and at the end (Da_A,f and Da_B,f) of each of our simulations are listed in Table 2. For a slow chemical reaction of two initially segregated reactants with the respective Da≪1, the chemical lifetime of the reactants is long enough for turbulence to mix them in the boundary layer at a rate fast enough that they are homogeneously distributed in a confined volume (e.g. a model grid) in a short time compared to the chemical lifetime. In that case chemical reaction can occur everywhere in that confined volume and the corresponding volumetric-mean production rate of tracer C can be approximated as the product of k and the volumetric means of the reactants (k〈A〉〈B〉). On the other hand, for fast reactions with the respective Da>1, the reactants are not well mixed by turbulence during the theoretical chemical lifetime of the reactants, such that they can react only in a fraction of the confined volume where turbulent motion brings the two species together. In that case the volumetric-mean of the production rate of tracer C also depends on the turbulent timescale and would be overestimated by k〈A〉〈B〉. Hence, the Damköhler number can indicate whether the effect of turbulent mixing on a chemical reaction is significant in the flow. For instance, the eddy turnover timescale for large-scale eddies in a CBL is typically 10²–10³ s, to which the chemical lifetime of NO_X and O₃ in the boundary layer is comparable (10²–10⁵). Therefore, the turbulent motions in the CBL potentially affect any chemical reaction occurring at a rate higher than that between NO and O₃.

To quantify the chemistry–turbulence interaction at the molecular diffusion spatio-temporal scale, the Kolmogorov Damköhler numbers are also calculated in some of simulations. The definition of Kolmogorov Damköhler number is adopted from Vilà-Guerau de Arellano et al. (2004), i.e.

$$D{a}_{\text{k,A}}={\displaystyle \frac{t_{\mathrm{k}}}{t_{\text{chem,A}}}}$$

for tracer A and similarly for tracer B (Da_k,B) by replacing the denominator with t_chem,B. The Kolmogorov timescale t_k is around 10 s for all the simulations (Garcia and Mellado, 2014).

The second number, the effective chemical reaction rate (k_eff), can quantitatively measure the actual chemical reaction rate averaged in a confined volume under the consideration of chemical segregation due to inefficient turbulent mixing (Vinuesa and Vilà-Guerau de Arellano, 2005), and it can be written as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{k}_{\text{eff}}=\left(\mathrm{1}+{\displaystyle \frac{\langle A^{\prime }{B}^{\prime }\rangle }{\langle A\rangle \langle B\rangle }}\right)k,\end{array}$$

where the angle-bracketed and primed terms are the means and the deviations from the means, respectively, of the Reynolds-decomposed concentrations.

The error induced by neglecting such chemical segregation in a complete-mixing model, in which both tracers A and B are assumed to be evenly distributed throughout the boundary layer, can be written as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }E=\mathrm{1}-{\displaystyle \frac{k_{\text{eff}}}{k}}=-I_{\text{S}},\end{array}$$

where the segregation coefficient

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{I}_{\text{S}}={\displaystyle \frac{\langle A^{\prime }{B}^{\prime }\rangle }{\langle A\rangle \langle B\rangle }}\end{array}$$

represents the state of mixing of tracers A and B in the boundary layer (e.g. Danckwerts, 1952; Vinuesa and Vilà-Guerau de Arellano, 2003). If tracers A and B are completely mixed, then I_S=0. If tracers A and B are completely segregated, then $I_{\text{S}}=-\mathrm{1}$. If tracers A and B are emitted simultaneously and at the same location so that their concentrations are correlated, then I_S>1. With the DNS model fully resolving the chemical segregation in the boundary layer, the error from a complete-mixing or coarse-grid model can be calculated by

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }E=k_{\text{eff, cm}}-k_{\text{eff, DNS}},\end{array}$$

where k_{eff, cm} is the calculated boundary-layer-averaged effective chemical reaction rate in a complete-mixing or coarse-grid model, and k_{eff, DNS} is the calculated value in the corresponding DNS model.

When one considers a confined volume within a horizontal layer at a specific height z, the corresponding effective chemical reaction rate also varies with height. The height-dependent effective chemical reaction rate can be written as a function of z:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\left[k_{\text{eff}}\right]\left(z\right)=\left(\mathrm{1}+{\displaystyle \frac{\left[A^{\prime }{B}^{\prime }\right]}{\left[A\right]\left[B\right]}}\right)k=(\mathrm{1}+\left[I_{\text{S}}\right](z\left)\right)k.\end{array}$$

The square-bracketed terms denote the horizontal means of the quantities. The corresponding height-dependent error induced by neglecting chemical segregation by a model that assumes tracers A and B are well-mixed horizontally within a horizontal layer at z is then

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\left[E\right]\left(z\right)=\mathrm{1}-{\displaystyle \frac{\left[k_{\text{eff}}\right]\left(z\right)}{k}}=-\left[I_{\text{S}}\right]\left(z\right),\end{array}$$

where the height-dependent segregation coefficient is

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\left[I_{\text{S}}\right]\left(z\right)={\displaystyle \frac{\left[A^{\prime }{B}^{\prime }\right]}{\left[A\right]\left[B\right]}}.\end{array}$$

3 Results of the DNS runs

3.1 Homogeneous emissions

3.1.1 Reference cases with rural emission flux

For both VV05 runs, all the tracers are relatively well mixed in the mixed layer. Tracer A is largely consumed in the mixed layer, while tracer B is in excess so that its concentration remains essentially constant. The reaction between tracers A and B can then be considered a pseudo first-order reaction, and the corresponding production term $kAB\approx k^{\prime }A$ (plotted in the top panel of Fig. 2) is dependent only on the concentration of tracer A, where $k^{\prime }=kB$ is the pseudo first-order rate coefficient (Hobbs, 2000; Jacobson and Jacobson, 2005). In this situation, tracer A is said to be the limiting reactant, and the chemical reaction between tracers A and B is tracer-A limiting (Zumdahl, 1992). As the production of tracer C is also tracer-A limiting, the vertical fluxes of tracer C (the magenta lines in Fig. 3) in these two cases are both positive, similar to those of tracer A (not shown), indicating that tracers A and C are both correlated with the updraught. As the chemistry shifts from slow to fast, the height with the maximum flux of tracer C transits from near the top of the boundary layer to near the surface. This is due to the increasing lifetime of tracer B with an increasing chemical reaction rate. As more tracer A is consumed, more tracer B can flow down the boundary layer without being consumed; hence, more tracer C is produced at a lower altitude in the boundary layer where tracer A is available. As indicated by a large Damköhler number of tracer A ($D{a}_{\mathrm{A},f}=\mathrm{5.71}>\mathrm{1}$), the mixing between tracers A and B hence becomes less efficient, resulting in a larger segregation between the two tracers.

Figure 2 Colour maps of the distribution of the production term (kAB) at the end of the simulations with homogeneous emissions for (a) the cases slow-VV05 and (b) slow-mflux. (c) Same colour map for the same case as in (a) but with time averaged for the last 30 min of the simulation. (d) Same colour map for the same case as in panel (a) but in the 1 km-64 lev coarse-grid model.

Figure 3 Vertical profiles of the horizontally averaged vertical flux of tracer C $\left[C^{\prime }{w}^{\prime }\right]$ with homogeneous emissions for varied cases time averaged for the last 30 min of the simulations.

3.1.2 Comparison with previous LES studies

We use the results of the two VV05 runs to compare with those presented in the LES study of Vinuesa and Vilà-Guerau de Arellano (2005). The horizontal and vertical resolutions of their LES model are 50 and 25 m, respectively, roughly 4 times coarser than the resolutions of our DNS model. From our DNS model, the resultant boundary-layer-averaged effective chemical reaction rates k_eff are 96.5 % and 86.9 % of the imposed rate k for our slow-VV05 and fast-VV05 runs, respectively (see last column of Table 2). These rates are also similar to the values of 96 % and 85 % reported in Vinuesa and Vilà-Guerau de Arellano (2005). Therefore, our DNS model gives similar boundary-layer-averaged values of k_eff as the LES model adopted in Vinuesa and Vilà-Guerau de Arellano (2005). When comparing our vertical profiles of the horizontally averaged effective chemical reaction rate (Fig. 4) with the vertical profiles of the segregation coefficient presented in Vinuesa and Vilà-Guerau de Arellano (2005), the profiles of the two studies in general share similar shapes, with the largest magnitude of segregation all occurring just below the top of the boundary layer. For instance, the maximum values of [k_eff] are 83 % and 62 % of k in our slow-VV05 and fast-VV05 runs, respectively.³ However, our profiles show more prominent minima near the top of the surface layer, which shows the effect of the higher vertical resolution in the DNS model.

Figure 4 Vertical profiles of the horizontally averaged normalised effective chemical reaction rate [k_eff]∕k with homogeneous emissions for varied cases time averaged for the last 30 min of the simulations. Note that the vertical profiles are not plotted at altitudes where the concentration of tracer B is equal to zero.

3.1.3 The effect of strong emission fluxes

When the emission fluxes of tracer A increase to urban values beyond 0.25 ppb m s⁻¹ in the mflux, sflux and ssflux runs, tracer A accumulates and is in excess in most of the boundary layer. On the other hand, tracer B is almost completely depleted in the lower part of the boundary layer. As the emission flux of tracer A increases and/or the chemistry becomes faster, the chemical timescale of tracer B shortens, so the depth where tracer B can be replenished from the free troposphere becomes narrower. Tracer C is mostly produced near the top of the boundary layer, as indicated in the colour map of the production term in the middle panel of Fig. 2. The chemical system now transits to the so-called “diffusion-limited” reaction (Sykes et al., 1994), where the reaction between tracers A and B depends on the availability of tracer B and hence is tracer-B limiting.

The shift of the reaction from tracer-A limiting to tracer-B limiting can also be seen from the profiles of the vertical fluxes of tracer C in Fig. 3 (sflux in red and ssflux in cyan). The fluxes of tracer C now become negative, indicating that tracer C is downdraught correlated, as is tracer B. As the emission flux of tracer A increases, the maximum of the production term shifts from the top of the surface layer to near the top of the boundary layer closer to the source of tracer B. These cases are also characterised by very large values of Da_B,f instead of a large Da_A,f as in the fast-VV05 run (see Table 2). Since the chemical transformation of the reactant that is relatively less abundant than the other reactant is influenced more by turbulent mixing (Vilà-Guerau de Arellano et al., 2004), Da_B,f is now a better indicator for the role of turbulent mixing on chemical reactions than Da_A,f. Summarising the simulations with homogeneous emissions, one can observe from Fig. 8 (black circles) that the deviation of k_eff from the imposed rate k (or the error from the complete-mixing model) increases with the increased Damköhler number of the limiting reactant (Da_lim), where $D{a}_{\text{lim}}=D{a}_{\mathrm{A},f}$ when the reaction is tracer-A limiting and $D{a}_{\text{lim}}=D{a}_{\mathrm{B},f}$ when the reaction is tracer-B limiting. This transition occurs when Da_lim reaches the order of 10. A similar concept can also be applied to the corresponding Kolmogorov Damköhler number, where the limiting Kolmogorov Damköhler number (Da_k,lim) is introduced (refer to the upper x axis of Fig. 8). The errors induced by neglecting the subgrid chemical segregation become significant when Da_k,lim reaches the order of around 10⁻¹.

For the runs with urban emission fluxes, the boundary-layer-averaged effective chemical reaction rate k_eff is very low even with slow chemistry. With slow chemistry, k_eff is only 23.5 %, 12.3 % and 8.2 % of the imposed rate k in the mflux, sflux and ssflux runs, respectively. With fast chemistry, k_eff further drops to 4.6 %, 3.0 % and 2.3 % of k, respectively. These low values of k_eff imply that if the surface emission of a pollutant is so strong that the reaction is limited by the availability of the entrained reactant, the resultant overestimation of the actual chemical reaction from an assumption of complete mixing can be substantial (in our cases up to 98 %) for any reaction with reaction rate equivalent to or faster than that between NO and O₃. The vertical profiles of the horizontally averaged effective chemical reaction rate in Fig. 4 show that the maximum segregation between tracers A and B still occurs near the top of the boundary layer. With fast chemistry, the effective chemical reaction rate at that altitude can even drop to around 10 % of the imposed rate.

3.2 Heterogeneous emissions

With the emissions of tracers A and B heterogeneously distributed on the surface, the segregation is maximum at the surface with its magnitude decreasing with increasing altitude, as seen from the vertical profiles of the horizontally averaged segregation coefficient ([I_S](z)) in Fig. 5. It is because under this setting mixing is most difficult near the surface where the tracers are emitted separately. This description of the profiles agrees with the results of the simulations with heterogeneous emissions presented in Krol et al. (2000) and Vinuesa and Vilà-Guerau de Arellano (2005), despite that they implemented a different emission configuration with a Gaussian function imposed on the surface emission flux. To examine how the segregation between tracers A and B changes with the emission flux and the imposed chemical reaction rate, the results of the set of simulations with the length of heterogeneity (dx) of 2 km (all plotted in red) are first compared. In the slow-2 km-VV05 run (dashed red line), the magnitude of [I_S] decreases with increasing altitude in the boundary layer, so at a height above 0.8z_i, [I_S] becomes larger than zero, meaning that tracers A and B are not only well mixed but also that their flows are correlated. This results in [k_eff]>k, where the assumption of complete mixing in the model grid will underestimate the reaction rate. This phenomenon is also reported in Ouwersloot et al. (2011), showing [I_S]>0 at higher altitudes in their LES simulations with heterogeneous isoprene emission fluxes. On the other hand, as the surface emission flux increases to 0.5 ppb m s⁻¹ in the slow-2 km-sflux run (red solid line), [I_S]>0 only at the height above 0.98z_i.

Figure 5 Vertical profiles of the horizontally averaged normalised effective chemical reaction rate [k_eff]∕k (and the horizontally averaged segregation coefficient [I_S] (upper x axis)) with heterogeneous emissions for varied cases time averaged for the last 30 min of the simulations.

When shifting from slow to fast chemistry (dotted red line), the magnitude of the segregation further increases, resulting in a boundary-layer-averaged value of −0.95. The segregation at the top of the boundary layer remains high, with $\left[I_{\mathrm{S}}\right]=-\mathrm{0.8}$. Sometimes the segregation caused by an increase in imposed reaction rate k can reduce k_eff so much that merely increasing k can no longer increase the production of tracer C. This happens between the two 6 km-sflux runs (check Table 2 for their resultant k_eff). When shifting from slow chemistry to fast chemistry, k_eff drops 9.18 times, such that it in turn cancels the effect of the 10 times increase in k. This results in similar amounts of tracer C produced in the two 6 km-sflux runs, with the mixing ratio of tracer C equal to 1.58 and 1.60 ppm at the end of the simulations for the slow- and fast-chemistry cases, respectively.

To discuss the effect of the length of heterogeneity (dx) on segregation, the three solid lines in the left panel of Fig. 5, which indicate the slow-sflux runs at different dx (dx=1 km plotted in blue, dx=2 km in red and dx=6 km in green), are compared. In general, the segregation increases with increasing dx. This agrees with the conclusion from Ouwersloot et al. (2011). For the cases with dx=1 km and dx=2 km, one can still observe a decrease in the magnitude of the segregation with increasing altitude, with the corresponding boundary-layer-averaged I_S equal to −0.49 and −0.69, respectively. However, for dx=6 km, the segregation remains almost constant throughout the mixed layer and starts to decrease only below the top of the boundary layer. This also results in a very large segregation throughout the boundary layer, with the boundary-layer-averaged I_S equal to −0.97. This shows a similar phenomenon, as described in Ouwersloot et al. (2011), that as the length of heterogeneity exceeds a few boundary layer heights,⁴ i.e. ∼3 km at the end of the simulations, the boundary layer between the two patches that emit tracers A and B, respectively, barely mix, and the system behaves like two individual boundary layers.

Similar to the cases with homogeneous emissions, the increased segregation with increasing surface emission flux and imposed chemical reaction rate can be indicated by the relatively large values of the final Damköhler numbers. Since the boundary-layer-averaged Damköhler numbers of tracers A and B are statistically the same in the simulations with heterogeneous emissions, here we take the mean of the two numbers to get the averaged final Damköhler number Da_f (equal to $\mathrm{1}/\mathrm{2}(D{a}_{\mathrm{A},f}+D{a}_{\mathrm{B},f})$). We then fit Da_f with the normalised boundary-layer-averaged k_eff∕k for the simulations with each dx (check Table 2 for the values). The results are plotted in Fig. 6. In general, log (k_eff∕k) follows a square law with log (Da_f), i.e. $\log (k_{\text{eff}}/k)=-a_{\text{fit}}(\log D{a}_f+\mathrm{1}{)}^{\mathrm{2}}$, where a_fit is the fitting coefficient varying with dx. Or in another way, k_eff∕k follows a Gaussian function of Da_f, i.e. $k_{\text{eff}}/k=\exp (-a_{\text{fit}}(\log D{a}_f+\mathrm{1}{)}^{\mathrm{2}})$. The x intercept is chosen to be −1 because k_eff≈k when Da_f≲0.1 (or $\log D{a}_f\sim -\mathrm{1}$). As dx increases, a_fit also increases (see the values in Fig. 6). However, their relation is unlikely to be linear.

Figure 6 Fitted curves of the normalised boundary-layer-averaged effective chemical reaction rate k_eff∕k against the final Damköhler numbers Da_f for the simulations with heterogeneous emissions with different lengths of heterogeneity (dx). a_fit refers to the fitted slope for the square law between log (k_eff∕k) and log (Da_f) corresponding to different dx values.

Figure 7 Plot of the vertical grid spacing (dz) (in equivalent to vertical resolution) vs. height (z) in the DNS model (black), the 64 lev coarse-grid models (red) and the 32 lev coarse-grid models (blue) in the first 5 km.

4 Estimating errors from regional models

In the previous sections, the boundary-layer-averaged effective chemical reaction rate is compared with the imposed rate in a complete-mixing model that assumes the tracers to be completely well mixed in the whole boundary layer. However, the horizontal resolutions in regional chemical-transport models are comparatively high (of the order of a few kilometres) with multiple vertical levels within the boundary layer. Such models are often employed when modelling urban areas. To evaluate the importance of the subgrid chemistry–turbulence interaction in these regional chemical-transport models, we degrade our DNS model to coarse-grid models with two horizontal resolutions (1 and 3 km, with 12 and 4 nodes in the horizontal direction, respectively) and two vertical resolutions (labelled 32 and 64 lev, with 32 and 64 nodes in the vertical direction in the whole domain, and with 10 and 20 levels within the boundary layer, respectively, when the boundary layer is fully grown), resulting in four coarse-grid models, named 1 km-64 lev, 1 km-32 lev, 3 km-64 lev and 3 km-32 lev. Figure 7 shows the respective vertical resolutions of the DNS model and the 64- and 32-level coarse-grid models vs. height. These resolutions are selected according to common grid sizes in regional chemical-transport models (e.g. Bouarar et al., 2017; Marécal et al., 2015).

Figure 8 Plot of the model errors E in percentage from the complete-mixing model (black circles) and the four coarse-grid models (+ and × in red and blue) against the corresponding Damköhler number of the limiting reactant (Da_lim) (and the corresponding Kolmogorov Damköhler number of the limiting reactant (Da_k,lim) (upper x axis)) for the simulations with homogeneous emissions.

The tracer concentration fields obtained from the DNS runs are interpolated from the high-resolution DNS model grids to the lower-resolution coarse-grid model grids. The volumetric averages of tracer concentrations in each model grid are then calculated. The statistics of these resolution-degraded concentration fields are then calculated as in Sect. 3. The chemical production terms from a coarse-grid model and the corresponding DNS run are plotted in the top and bottom panels of Fig. 2, respectively, to illustrate the effect of resolution degrading. In contrast to the DNS model, the coarse-grid model can no longer resolve the chemical segregation within each of its model grids. The boundary-layer-averaged normalised effective chemical reaction rates k_eff∕k of the four coarse-grid models averaged over the last 1.5 h of the corresponding simulations are subtracted from the values of k_eff∕k from our DNS model (last column of Table 2) to obtain the boundary-layer-averaged errors from the coarse-grid models induced by neglecting the subgrid chemical segregation (E) (see Eq. 8). These errors are listed in Table 3. The second column is the error from a complete-mixing model (E_cm), calculated from Eq. (6). Note that the errors are listed in percentages. Positive values indicate overestimation of the model from the actual chemical reaction rate, and negative values indicate underestimation. The dominance of positive values shows that the coarse-grid models, similar to the complete-mixing model, usually overestimate the actual chemical reaction rate.

Figure 9 Vertical profiles of the horizontally averaged normalised effective chemical reaction rate [k_eff]∕k from the four coarse-grid models with homogeneous emissions for the slow-mflux (in dark green), fast-mflux (in light green) and slow-sflux (in red) cases time averaged for the last 30 min of the simulations. Note that the vertical profiles are not plotted at altitudes where the concentration of tracer B is equal to zero.

Figure 10 Plot of model errors in percentage, categorised with different lengths of heterogeneity (dx), from the four coarse-grid models (E) against the corresponding errors from the complete-mixing model (E_cm) for the simulations with heterogeneous emissions.

4.1 Homogeneous emissions

We take a closer look at the results of the simulations with homogeneous emissions as described in Sect. 3.1. Figure 8 shows the model errors from the complete-mixing model (black circles) and the four coarse-grid models (+ and × marks) against the corresponding final Damköhler number of the limiting reactant (Da_lim; see the values in Table 2). One can see that in all runs the errors from all coarse-grid models are reduced in comparison to the complete-mixing model. The model errors from the coarse-grid models are largely reduced when Da_lim>10, which corresponds to the simulations with urban surface emission fluxes (mflux, sflux and ssflux). In those runs, the model errors drop from 77 %–98 % in the complete-mixing model to 36 %–69 % in the 32-level models and further to 15 %–51 % in the 64-level models. However, this also means that even with the highest resolution the model errors are still noticeably significant. Also, for five out of seven runs (the slow-VV05 run is excluded here due to the insignificant errors), a larger reduction of the model error is observed with an increased vertical resolution (from blue + to blue ×) than with an increased horizontal resolution (from blue + to red +). The reason behind this is that tracers A and B are segregated along the vertical direction.

Figure 9 shows the vertical profiles of the horizontally averaged error ([E](z); see Eq. 10) of the four coarse-grid models for the slow-mflux and slow-sflux runs. All coarse-grid models fail to show the local minimum at the top of the surface layer, as the surface layer is too thin for these models to resolve. They also show a less prominent minimum at the top of the boundary layer. The distinct underestimation of [E] in the surface layer and entrainment zone indicates that high vertical resolution is especially important to resolve the chemical segregation around these two zones.

4.2 Heterogeneous emissions

The errors from the four coarse-grid models (E) in the simulations with heterogeneous emissions are then plotted against the errors from the corresponding complete-mixing model (E_cm) in Fig. 10. These data points are categorised with different lengths of heterogeneity (dx) implemented in the corresponding simulations. The dashed line indicates the value where E=E_cm. The points below the dashed line show the cases for which the corresponding coarse-grid model improves from the complete-mixing model, while those above the dashed line show a deterioration. Contrary to the simulations with homogeneous emissions the coarse-grid models do not always perform better than the complete-mixing model, especially when the horizontal model grid size is larger than the length of heterogeneity. It is because under this condition, the model grid fails to resolve the heterogeneity of the emission sources. This also applies to the simulations with dx=1 km in the 1 km resolution coarse-grid models, as the model grids are slightly offset from the emission patches. In these runs, the coarse-grid models mix the originally segregated tracers in the lower layers at a higher concentration than the complete-mixing model, in which the latter artificially mixes the tracers all over the whole boundary layer and dilutes the tracer concentrations. That is why the coarse-grid models give larger overestimations. The model with the higher vertical resolution overestimates the reaction rate even more, as its first layers are thinner and hence contain even higher concentrations of the artificially mixed tracers in the lower-level grids.

For the simulations with dx=6 km, the coarse-grid models always perform better than the complete-mixing model. In these simulations, the horizontal grids in the coarse-grid models can always resolve the emission heterogeneity. Note that the reduction of the errors from the coarse-grid models is larger for the increased horizontal than vertical resolution, as tracers A and B are now segregated in the horizontal direction.

Table 3 Resultant model errors in percentage from the complete-mixing model (E_cm) and the four coarse-grid models (E) for the simulations with (a) homogeneous emissions and (b) heterogeneous emissions.

5 Discussions

5.1 Implications for regional modelling of urban environments

Although the DNS runs in this work are only conducted in idealised conditions, they can provide insights into and estimations of the errors induced by neglecting the subgrid chemical segregation due to inefficient turbulent mixing in larger-scale models. The simulations with homogeneous emissions (Sect. 3.1) show significant errors in the low-resolution model under strong emission flux conditions, which are indicated by a large (>1) Damköhler number of the limiting reactant (Da_lim). This suggests that Da_lim>1 can be taken as a condition in the model under which the chemical calculation should be corrected by taking the subgrid chemical segregation into account. One way to implement such correction is to substitute the imposed chemical reaction rate k by the effective rate k_eff which, as shown in Fig. 8, can be expressed as a function of Da_lim. The results from the simulations with heterogeneous emissions (Fig. 6) suggest that k_eff additionally depends on the length of heterogeneity (dx). When considering increasing model resolution to reduce the errors induced by neglecting subgrid chemical segregation, we learn from Sect. 4 that whether it is more effective to increase the horizontal or vertical resolution depends on the direction of the initial segregation of the reactant sources. High resolution is in particular necessary around the surface layer and entrainment zone. When the reactant sources are horizontally segregated, it is essential that the horizontal resolution of the model is fine enough to resolve the emission heterogeneity. Otherwise, increasing the vertical resolution can induce even larger model errors.

There are additional points to note in our estimations. When arriving at our conclusion of the dependency of k_eff on Da_lim, we increase Da_lim by shortening the corresponding chemical timescale with increased imposed k and emission fluxes. One should not neglect the dependency of Da_lim on the buoyancy fluxes, which determine the turbulent timescale. However, we expect to see a similar k_eff dependency on Da_lim when one repeats the simulations with increasing buoyancy fluxes instead.

By increasing the emission fluxes to urban values in Sect. 3.1.3, it is shown that the effective chemical reaction rate can be much smaller than the values reported in previous studies under rural conditions (e.g. Vinuesa and Vilà-Guerau de Arellano, 2003, 2005). However, it should be also noted that the chemical segregation is particularly large in these simulations, because one of the reactants (in this case tracer B) is depleting. Despite the large deviation of k_eff from k, the production term kAB is anyway small due to the low concentration of tracer B. Therefore, the actual effect on the calculated concentration of tracer B may not be as large as one may expect from the overestimation of k_eff. A similar reasoning also applies to tracer A. When the emission flux of tracer A is large, the main process affecting the concentration of tracer A is emission instead of chemistry. However, the effect of chemical segregation on the concentrations of the product tracer C is still worth noticing.

Unlike other studies (e.g. Ouwersloot et al., 2011; Dlugi et al., 2019; Kim et al., 2016; Li et al., 2016, 2017) in which multiple-reaction chemical systems are employed, our work mostly focuses on an idealised second-order chemical reaction of two non-specific reactive species. This approach allows us to interpret our work for any second-order chemical reactions. For a chemical species involved in a multiple-reaction chemical system, like O₃, one can still calculate the net impact of chemical segregation by considering the errors of all reactions in which the species is involved. However, it is also important to notice that the net impact of chemical segregation on such a species would in general be reduced with the increasing complexity of the chemical system, because cycling reactions tend to lengthen the chemical timescale and reduce the corresponding Damköhler number and hence reduce the chemical segregation. For example, with emission fluxes of NO_X comparable to our sflux and ssflux cases, the segregation coefficient between isoprene and OH is only −0.05 and −0.17 in the high-NO_X and very-high-NO_X cases of Kim et al. (2016), respectively. Li (2019) also performed simulations similar to this study with the $\mathrm{NO}-{\mathrm{NO}}_{\mathrm{2}}-{\mathrm{O}}_{\mathrm{3}}$ triad, and they found less significant errors than with the A–B–C chemical system. Also, this error does not increase with Da_lim when Da_lim reached the order of 10, contrary to the increasing trend observed in the A–B–C chemical system. Therefore, we expect the magnitude of chemical segregation to be smaller than the values reported in this study when a multiple-reaction chemical system is included. Future studies can also consider implementing a backward reaction in addition to the idealised second-order chemical reaction to approximate the effect of a multiple-reaction chemical system.

The work in Sect. 4 gives an estimation of the error in a regional chemical-transport model when neglecting subgrid chemical segregation, where the model possesses relatively high horizontal and vertical resolution. The method we employed is to average the concentration fields within the volume of each model grid. Current regional chemical-transport models often include additional measures to prevent accumulation of reactants close to the emission sources, such as distributing the emitted species in the first few vertical layers (e.g. Freitas et al., 2011) and employing boundary layer schemes to parameterise eddy transport in the lower troposphere (e.g. Hong et al., 2006). These measures may diminish the excessive overestimation in the coarse-grid models reported in Sect. 4.2 in cases in which the model fails to resolve emission heterogeneity. However, these measures may introduce additional artificial mixing in other cases. Also, there is an inherited difference in modelling framework that makes the comparison between regional chemical-transport models and turbulence-resolving models like DNS and LES difficult. The deficiency of regional models to only parameterise the vertical mixing in the boundary layer causes its inaccuracy in representing the boundary layer dynamics. For instance, by comparing the model results from WRF-Chem and from the NCAR LES model, Li et al. (2019) reported a weaker vertical mixing in the WRF-Chem simulations than in the LES simulations, causing the upward transports of surface-emitted chemicals and hence the chemical production of OH and O₃ to be undermined aloft in the mixed layer. Our averaging method employed in the coarse-grid models may not be able to account for these undermined upward transports in regional chemical-transport models. Nevertheless, our estimation shows that the model resolution clearly plays an important role in resolving chemical segregation and providing accurate chemical calculations in the models.

One important aim of the study of chemistry–turbulence interaction is to provide a correction to the error in large-scale models induced by neglecting subgrid chemical segregation. While our work suggests a correction of k_eff with dependencies on Da_lim and dx, other works suggest that such correction should also include other variables such as updraught and/or downdraught fluxes (Petersen and Holtslag, 1999), variance of the reacting species, entrainment and/or emission fluxes (Petersen and Holtslag, 1999; Vinuesa and Vilà-Guerau de Arellano, 2003), turbulent fluxes (Dlugi et al., 2014), variance of the emission (Galmarini et al., 1997), magnitude and direction of mean horizontal wind (Ouwersloot et al., 2011), and distance from the sources (Karamchandani et al., 2000). The correction to the vertical exchange coefficient (i.e. the ratio between the turbulent flux and the mean gradient) can also serve as an alternative way to address the change in the transport terms from the effect of chemical reactions (Vilà-Guerau de Arellano and Duynkerke, 1993; Geyer and Stutz, 2004). There have been attempts to implement a parameterisation of subgrid chemistry–turbulence interaction to large-scale models. For example, Molemaker and Vilà-Guerau de Arellano (1998) suggested the use of lookup tables to indicate how chemical reaction rates should be corrected under different physical or chemical scenarios. Lenschow et al. (2016) developed a one-dimensional second-order closure model to account for the vertical turbulent mixing of chemical species ready to be incorporated into large-scale models. Given that the effect of subgrid chemical segregation is non-negligible under urban conditions, modellers should consider applying similar parameterisations in areas with intense emission and large source heterogeneity.

5.2 Other factors to be considered when studying chemical segregation in urban environments

Due to computational limitations of the DNS runs, some other conditions that may be important for turbulent mixing in the urban boundary layer were neglected in our work. First of all, the growth of the boundary layer is driven by a constant buoyancy flux, so the boundary layer height gradually increases. Hence, the simulation can only approximate the time when the boundary layer grows from sunrise to mid-afternoon. This also raises the issue of the time required for statistical equilibrium to be attained. In some cases, this time may be longer than the duration of daylight, after which the atmosphere is no longer convective. This poses a greater problem to cases with strong emission fluxes, as the time required to attain statistical equilibrium is even longer. It may not be practical to run longer than the duration of daylight even though statistical equilibrium is not reached. Second, our simulations only address a convective boundary layer under clear-sky conditions. We do not take other scenarios with different weather conditions (such as cloud-top boundary layer, e.g. Li et al., 2016, 2017) or when the atmosphere is stably stratified (typical for nighttime, e.g. Galmarini et al., 1997; Geyer and Stutz, 2004) into account. For instance, the perturbation in radiation with the existence of clouds closely couples to the turbulent effect (Mellado, 2017) and also modifies chemical reactions (Barth et al., 2002; Vilà-Guerau de Arellano et al., 2005). Third, it may be useful to consider other boundary conditions of the entrained tracers, such as varied entrainment fluxes (e.g. Krol et al., 2000; Albrecht et al., 2016) and varied concentrations in the free troposphere due to long-range transport (e.g. Zyryanov et al., 2012). Fourth, our simulations also do not include any mean horizontal winds (e.g. Ouwersloot et al., 2011) or forcings from the surface other than the constant buoyancy flux, such as varied surface heat fluxes (e.g. Ouwersloot et al., 2011; Van Heerwaarden et al., 2014) and surface roughness.

An important source of surface forcings in an urban boundary layer is undoubtedly from the urban structures (buildings and streets in the urban canopy). The structure of turbulent flow can be significantly altered in the street canyons due to the perturbation in radiation and the exchange of heat and momentum with the urban structures (e.g. Oke, 1997). These urban canopy meteorological forcings also affect vertical turbulent transport and hence the vertical distribution of chemical species (Huszar et al., 2020). Therefore, our DNS runs are more suitable for addressing the vertical mixing caused by the turbulent motions in the mixed layer relatively far away from the surface features that may induce other additional forcings. But in the mixed layer, the urban canopy still potentially affects the chemistry and dynamics in the boundary layer by means of surface roughness and emission heterogeneity. While in this work we have addressed the effect of emission heterogeneity, the effect of surface roughness and other configurations of emission patterns (e.g. Auger and Legras, 2007) can be further studied in the future. Previous work on natural canopies can also be combined into further studies with urban canopies, as the vegetation is mixed within the urban canopies in many cities.

6 Conclusions

We explore in this work the effect of chemical segregation due to inefficient turbulent mixing on chemical reactions in an urban boundary layer, and we estimate the resultant error in the regional chemical-transport models by conducting direct numerical simulation (DNS). As past studies mainly examined scenarios in forestal areas, we focus on urban conditions, specifically with strong emission fluxes and heterogeneous emissions, as both factors can potentially increase chemical segregation.

With homogeneous emissions, our simulations give similar results as past studies using large-eddy simulations when the emission flux of the surface-emitted tracer A is of rural value, in spite of the increase in resolution of our DNS model. On the other hand, increasing the emission flux of tracer A to urban values depletes the entrained tracer B, so its availability limits the reaction. In this situation, the segregation between tracers A and B becomes very large, which results in significant overestimation of the effective chemical reaction rate k_eff in a complete-mixing model. Such cases are always characterised by a very large Damköhler number of tracer B, the limiting reactant, suggesting that the Damköhler number of the limiting reactant Da_lim can be used as an indicator of the effect of turbulent mixing on chemical reaction. When Da_lim>1, chemical calculations in low-resolution models should be corrected by taking chemical segregation into account, and k_eff can be used for such a correction, which can be expressed as a function of Da_lim. From our simulations with heterogeneous emissions, in which tracers A and B are emitted from the surface in alternate patches of different widths, k_eff normalised by the imposed rate k is further found to be in a Gaussian function of the final Damköhler number of the reactants, with the slope of the fitted Gaussian law increasing with the length of heterogeneity. The large overestimation of k_eff reported from some cases of our simulations shows that the segregation between heterogeneously distributed surface-emitting pollutants is potentially important and should not be neglected especially when the corresponding emission flux and/or the length of heterogeneity are large.

To evaluate the errors induced by neglecting the chemical segregation in regional chemical-transport models with higher resolution than a complete-mixing model, we degraded our DNS model to coarse-grid models with two horizontal and vertical resolutions commensurable to regional models. With homogeneous emissions, all the coarse-grid models give smaller errors than the complete-mixing model. Yet the errors from the coarse-grid models remain high for simulations with urban emission fluxes. The improvement is more significant for the increased vertical resolution instead of horizontal resolution, as the initial segregation between tracers A and B is in the vertical direction. All coarse-grid models give the largest overestimations of the height-dependent effective chemical reaction rate near the top of the surface layer and the entrainment zone, indicating that high resolution is most important in these areas. With heterogeneous emissions, the coarse-grid models perform worse than the complete-mixing model when the coarse horizontal model grid cannot resolve the emission heterogeneity. With higher vertical resolution, the respective coarse-grid model gives an even larger error. This illustrates that increasing the model resolution may not improve the model performance when the enhanced model resolution still fails to resolve the emission heterogeneity. For the coarse-grid models which can resolve the emission heterogeneity, the model improvement of the coarse-grid model is more significant for increased horizontal than vertical resolution, as in these cases the initial segregation is in the horizontal direction. This suggests that whether the model improvement is more sensitive to the increase in the horizontal or vertical resolution depends on the direction of initial segregation between the reactants.

Appendix A: Non-dimensionalisation of the model

Our results from the DNS runs are based on the data in the fully developed turbulent regime in the simulation that is established after the initial transient phase. In that regime, the initial conditions have been sufficiently forgotten, and the parameters $\left\{\mathit{\nu },\mathit{\kappa },F_b,N\right\}$ define the system completely. N and F_b are chosen to non-dimensionalise the equations (Fedorovich et al., 2004; Garcia and Mellado, 2014), such that the system yields a reference timescale T₀ as N⁻¹, a reference length scale L₀ as (F_b∕N³)^½ and a reference velocity scale U₀ as (L₀F_b)^⅓. With the Prandtl number ($Pr=\mathit{\nu }/\mathit{\kappa }$) set to 1, the system only depends on the reference buoyancy Reynolds number:

$$R{e}_{\mathrm{0}}={\displaystyle \frac{F_b}{\mathit{\nu }{N}^{\mathrm{2}}}},$$

which is set to 15 000 in all simulations. Re₀ refers to the Reynolds number in the entrainment zone (Garcia and Mellado, 2014). In this configuration, the Reynolds number is ∼1000 in the mixed layer. Conducting DNS with Re in typical atmospheric condition (∼10⁷–10⁸) is computationally impossible. However, with Reynolds number similarity, one can perform DNS with moderate values, which still allows for certain extrapolation of the DNS results to corresponding atmospheric conditions (Dimotakis, 2000; Monin and Mahesh, 1998). Reynolds number similarity is an experimental observation of which properties of turbulent flows become practically independent of Re (Tennekes and Lumley, 1972; Mellado et al., 2018). With the values in our simulations, the Reynolds number in the whole boundary layer lies within the range where the Reynolds number similarity applies. One can refer to Sect. 2.1 for the more detailed explanation of the Reynolds number in DNS models.

The rate equations of the species are non-dimensionalised by introducing the characteristic scales for the mixing ratios. The characteristic scale for tracer A is controlled by its emission flux F_A, and it is therefore defined as $\langle A{\rangle }_{\mathrm{0}}=F_{\text{A}}{U}_{\mathrm{0}}^{-\mathrm{1}}$. For tracer B, the characteristic scale is set as its background concentration in the free troposphere (〈B〉₀). For tracer C, the characteristic scale is defined as $\langle C{\rangle }_{\mathrm{0}}=L_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}k\langle A{\rangle }_{\mathrm{0}}\langle B{\rangle }_{\mathrm{0}}$ to non-dimensionalise the concentration of tracer C with its dimensionless production term. The mixing ratios of the tracers are normalised with these scales so that the normalised mixing ratios are defined as $c_{\text{A}}=A/\langle A{\rangle }_{\mathrm{0}}$, $c_{\text{B}}=B/\langle B{\rangle }_{\mathrm{0}}$ and $c_{\text{C}}=C/\langle C{\rangle }_{\mathrm{0}}$. With these definitions, the rates of the three tracers can be expressed as the following equations:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{L}_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{d}{c}_{\text{A}}}{\mathrm{d}t}}=-K_{\text{A}}{c}_{\text{A}}{c}_{\text{B}},\\ {\displaystyle }& {\displaystyle }{L}_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{d}{c}_{\text{B}}}{\mathrm{d}t}}=-K_{\text{B}}{c}_{\text{A}}{c}_{\text{B}},\\ {\displaystyle }& {\displaystyle }{L}_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{d}{c}_{\text{C}}}{\mathrm{d}t}}=-c_{\text{A}}{c}_{\text{B}},\end{array}$$

where the dimensionless rate constants are $K_{\text{A}}=L_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}k\langle B{\rangle }_{\mathrm{0}}$ and $K_{\text{B}}=L_{\mathrm{0}}{U}_{\mathrm{0}}^{-\mathrm{1}}k\langle A{\rangle }_{\mathrm{0}}$.