Error induced by neglecting subgrid chemical segregation due to inefﬁcient 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 inefﬁcient 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 conﬁguration with an A-B-C second-order chemical scheme, urban surface emission 5 ﬂuxes 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 is increasingly overestimated by the imposed rate with an increased horizontal scale of emission heterogeneity. Coarse-grid 10 models with resolutions commensurable to regional models give reduced yet still signiﬁcant 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. Code and data availability. Source ﬁles of the DNS model TLab and further documentation can be found at https://github.com/turbulencia/tlab (accessed 07 April 2020). Primary data and scripts used in the analysis and other supplementary information that may be useful in repro- ducing the author’s work are archived by the Max Planck Institute for Meteorology and can be obtained at http://hdl.handle.net/21.11116/ 0000-0006-11C3-A.


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 stabiliy, 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 mode, 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 5 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 10 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 15 be of similar importance as the dynamical term in the flux equations for the NO−NO 2 −O 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 (LES), which can resolve the most energetic eddies in the boundary layer. 20 Many of these LES studies focus on the convective boundary layer (CBL), in which the imbalance between updraft and downdraft transport produces a large segregation of the reactants (Wyngaard and Brost, 1984;Chatfield and A. 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 25 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 ∼ 1, while k eff ∼ k when Da ∼ 0.1. Recently, the effect of segregation due to inefficient turbulent mixing on chemical 30 reaction has been considered as the cause of the miscalculation of large-scale models in a number of scenarios. One of the most discussed issues is the under-prediction 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 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 ∼ −50% 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-to very high-NO X surface emissions. Another possible factor is related to the emission heterogeneity 5 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 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 10 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 emission 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 urbanlike 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 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 DNS simulations with homogeneous emissions with an idealised second-order A-B-C chemistry scheme (A + B → 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 25 chemistry scheme can generally represent any second-order chemical reactions commonly seen in an urban environment.
The results from our DNS simulations are then degraded into 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 30 our simulation domain. However, this comparison is not suitable in application to urban environments, as the regional models simulating urban air quality typically have a much finer resolution with mesh size of the order of 1 km. Therefore, we degrade our DNS results to coarse-grid models with mesh size of 1 km 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 simulations 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. 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 10 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 15 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 3 − 10 4 (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 8 . Fortunately, relevant turbulence statistics such as variances, covariances and entrainment rates tend towards Reynolds-number independence once the Reynolds number increases beyond 10 3 − 10 4 (Dimotakis, 2000;Mellado 20 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 25 resolutions. Nonetheless, DNS is gaining traction as it approaches a degree of Reynolds number similarity that allows certain extrapolation to atmospheric conditions (Monin and Mahesh, 1998;Jonker et al., 2012Jonker et al., , 2013Waggy 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 30 resolutions (Sullivan and Patton, 2011). In addition, with strong emission fluxes, subgrid turbulent parametrisation in LES simulations may induce significant errors near the surface where the tracer is emitted (Vinuesa andPorté-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 Sections 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 model TLab to perform our computational experiments of turbulent mixing of reactive species in the convective boundary layer. The dynamical part of the simulations performed in this work follows similar settings in Garcia and Mellado (2014) and Van Heerwaarden et al. (2014). In the model, the Navier-Stokes equations 5 of incompressible fluids in the Boussinesq approximation are solved to calculate the buoyancy and the velocity fields: where u(x, t) is the velocity vector with components (u, v, w) in the directions x = (x, y, z) at time t, respectively. p is the 10 modified pressure divided by the constant reference density, and b(x, t) is the buoyancy, expressed as background buoyancy profile is set as b 0 (z) = N 2 z, where N 2 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.  The size of the computational grid is 720 × 720 × 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 120L 0 × 120L 0 × 34.4L 0 , where L 0 is the reference length scale. Together with the reference time and velocity time scales T 0 and U 0 , L 0 is used to non-dimensionalise the Navier-Stokes equations (see Appendix A). With typical atmospheric 25 parameters of L 0 ∼ 100 m and U 0 ∼ 1 m s −1 , 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 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 versus height in the first 5 km is shown in Figure 7 (black line). 30 We use a sixth-order compact scheme to calculate the spatial derivatives, and a forth-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 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), , Garcia and Mellado (2014) and .
The simulations are terminated after a total simulation time equivalent to 4.5 hours, during which the boundary layer grows 5 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 [b b ] 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 time scales 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 (Deardorff, 1970). This turbulent timescale will be used to calculate the Damköhler number Da in later sections.

Homogeneous emissions
An archetypical entrainment-emission configuration, as typically used by other LES past studies such as Krol et al. (2000) and 15 Vinuesa and Vilà-Guerau de Arellano (2005), is adopted in our DNS simulation (see left panel of Figure 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 0 ) is constantly and homogeneously fixed at 2 ppb. Tracer C is formed in a second-order chemical reaction between Tracer A and B: 20 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 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 25 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.  (2005)), mflux, sflux and ssflux respectively. As a reference, in the ssflux run, the char- 30 acteristic mixing ratio scale for Tracer A ( A 0 ) is equal to that of Tracer B ( B 0 ). These characteristic scales are used in the non-dimensionalisation of Equation 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 4.75 × 10 −4 and 4.75 × 10 −3 ppb −1 s −1 . In the real atmosphere, the rate of chemical reaction between NO and O 3 is 4.75 × 10 −4 ppb −1 s −1 , corresponding to the cases with slow chemistry. In the meanwhile, the reaction rate between isoprene and OH is 0.1 ppb −1 s −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 5 in Table 1.

Heterogeneous emissions
For the simulations with heterogeneous emissions, Tracers A and B are emitted alternately from patches on the surface with widths of 1 km, 2 km and 6 km respectively. The width of these emission patches are referred in this study as the length of 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 2 ) from the oxidation of volatile organic compounds (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 20 instance, the reaction rate between NO and acetone peroxy radical (CH 3 COCH 2 O 2 ), which is a common VOC species found in urban environments (Brasseur and Jacob, 2017;Li, 2019), is 0.503 ppb −1 s −1 (Manion et al., 2008), giving an example for the fast-chemistry cases in our simulations with heterogenous emissions.

Quantifying the chemical-turbulence interaction
Two numbers are mainly employed to quantify the effect of turbulent mixing on chemical reaction, namely the Damköhler 25 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 A + B → C, the Damköhler number of Tracer A is given by 1 Note the length of heterogeneity referred in this study is equivalent to half of the length denoted in Ouwersloot et al. (2011). 2 Note that the adopted emission fluxes are doubled from the values in the simulations with homogeneous emissions in order to conserve the total fluxes.
Similarly the Damköhler number of Tracer B can be written as 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 5 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 10 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 typical 10 2 − 10 3 s, to which the chemical lifetime of NO X and O 3 in the boundary layer is comparable (10 2 − 10 5 ). Therefore, the turbulent motions in the CBL potentially affect any chemical reaction occurring at a rate higher 15 than that between NO and O 3 .
To quantify the chemical-turbulence interaction at the molecular diffusion spatiotemporal 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.: 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 can be written as where the angle-bracketed and dashed terms are the means and the deviations from the means respectively, of the Reynoldsdecomposed concentrations.
The error induced by neglecting such chemical segregation in a complete-mixing model, in which both Tracer A and B are assumed to be evenly distributed throughout the boundary layer, can be written as where the segregaton coefficient 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 If Tracers A and B are emitted simultaneously and at the same location so that their concentrations are correlated, then 5 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 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. 10 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: 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 15 layer at z is then where the height-dependent segregation coefficient is 3 Results of the DNS simulations 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 as a pseudo first-order reaction, and the corresponding production term kAB ≈ k A (plotted in 25 the top panel of Figure 2) is dependent only on the concentration of Tracer A, where k = 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 Figure 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 updrifting air. 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 5 consumed, more Tracer B can flow down the boundary layer without being consumed, and 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 (Da A,f = 5.71 > 1), the mixing between Tracer A and B hence becomes less efficient, resulting in a larger segregation between the two tracers.
3.1.2 Comparison with previous LES studies 10 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) 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 Figure 2. The chemical system now transits to the so-called "diffusion-limited" reaction (Sykes et al., 1994), where the reaction between Tracers A and Tracer B depends on the availability of Tracer B, and hence is Tracer The shift of the reaction from Tracer A-limiting to Tracer B-limiting can also been seen from the profiles of the vertical fluxes of Tracer C on Figure 3 (sflux in red and ssflux in cyan). The fluxes of Tracer C now become negative, indicating that Tracer C is downdraft-correlated, as is Tracer B. As the emission flux of Tracer A increases, the maximum of the production 5 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  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 20 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 3 . The vertical profiles of the horizontally-averaged effective chemical reaction rate in Figure 4 show that the maximum segregation between Tracer 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. Sometimes the segregation caused by an increase of 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 6km-sflux runs (check Table 2 for their 10 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-time increase in k. This results in similar amounts of Tracer C produced in the two 6km-sflux runs, with the mixing ratio of Tracer C equal to 1.58 ppm 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 Figure 5, 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 Figure 6. In general, log(k eff /k) follows a square law with log(Da f ), i. e. log(k eff /k) = 30 −a fit (log Da f + 1) 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 , or k eff /k = exp(−a fit (log Da f + 1) 2 ). The x-intercept is chosen to be -1 because k eff ∼ k when Da f 0.1 (or log Da f ∼ −1). As dx increases, a fit also increases (see the values in Figure 6). Their relation is however unlikely to be linear.
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 in order of a few kilometres with multiple vertical levels within the boundary layer. Such models are often employed when modelling urban areas. To evaluate 5 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-km and 3-km, with 12 and 4 nodes in the horizontal direction respectively) and two vertical resolutions (32-lev 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  Table 2) to obtain the boundary layer-averaged errors from the coarse-grid models induced by neglecting the subgrid chemical segregation (E) (see Equation 8). These errors are listed in Table 3. The second column is the error from a complete-mixing model (E cm ), calculated from Equation 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 25 to the complete-mixing model, usually overestimate the actual chemical reaction rate.

Homogeneous emissions
We take a closer look at the results of the simulations with homogeneous emissions as described in Section 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 30 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-lev models, and further to 15 -51% in the 64-lev models. However, this also means that even with the highest resolution the model errors are still noticeably significant. Also, for 5 out of 7 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 is that Tracers A and B are segregated along the vertical direction.
5 Figure 9 shows the vertical profiles of the horizontally-averaged error ([E](z), see Equation 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. 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 coarsegrid models, as the model grids are slightly offset from the emission patches. In these runs, the coarse-grid models mix the 20 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. 25 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.

Implications for regional modelling of urban environments
Although the DNS simulations in this work are only conducted in idealised conditions, they can provide insights 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 (Section 3.1) show significant errors in the low-resolution model under  (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 25 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. (2016Li et al. ( , 2017) in which multiple-reaction chemical systems are employed, our work mostly focuses on an idealised second-order chemical reaction of 30 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 3 , 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-and very high-NO X cases of Kim et al. (2016) respectively. Li (2019) also performed simulations similar to this study with the NO−NO 2 −O 3 triad, and found less significant errors than with the A-B-C chemical system. Also, this error does not increase with Da lim 5 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 Section 4 gives an estimation of the error in a regional chemical-transport model when neglecting subgrid 10 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 first few vertical layers (e. g. Freitas et al. (2011)) and employing boundary layer schemes to parametrise eddy transport in the lower troposphere (e. g. Hong et al. (2006)). These measures may diminish the excessive overestimation in 15 the coarse-grid models reported in Section 4.2 in cases of 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 parametrise 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 20 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 3 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 on resolving chemical segregation and providing accurate chemical 25 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 model induced by neglecting subgrid chemical segregation. While our work suggests a correction of k eff with dependencies on Da lim and dx, other work suggest that such correction should also include other variables such as updraft/downdraft fluxes (Petersen and Holtslag, 1999), variance of the reacting species, entrainment/emission fluxes (Petersen and Holtslag, 1999;Vinuesa and 30 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  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 parametrisation 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-dimension 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 parametrisations in areas 5 with intense emission and large source heterogeneity.

Other factors to be considered when studying chemical segregation in urban environments
Due to computational limitations of the DNS simulations, 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 that the boundary layer height gradually increases. Hence, the simulation can only approximate the time 10 when the boundary layer grows from sunrise to mid-afternoon. This also raises the issue of the time required for statistical equilibrium to attain. 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 15 other scenarios with different weather conditions (such as cloud-top boundary layer, e. g. Li et al. (2016Li et al. ( , 2017) or when the atmosphere is stably stratified (typical for nighttime, e. g. Galmarini et al. (1997); Geyer and Stutz (2004) An important source of surface forcings in an urban boundary layer is undoubtedly from the urban structures (buildings 25 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 on the vertical distribution of chemical species (Huszar et al., 2020). Therefore, our DNS simulations are more suitable in 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. 30 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 to further studies with urban canopies as the vegetation are mixed within the urban canopies in many cities.

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 estimate the resultant error in the regional chemical-transport models by conducting direct numerical 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 25 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 30 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.
Code and data availability. Source files of the DNS model TLab and further documentation can be found at https://github.com/turbulencia/tlab (accessed 07 April 2020). Primary data and scripts used in the analysis and other supplementary information that may be useful in repro- . With the Prandtl number (P r = ν/κ) set to 1, the system only depends on the reference buoyancy Reynolds number 15 which is set to 15000 in all simulations. Re 0 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 7 − 10 8 ) is computationally impossible. However, with Reynolds number similarity, one can perform DNS with moderate values, which still allows 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 20 of turbulent flows become practically independent of Re (Tennekes et al., 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 Section 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 is therefore defined as A 0 = F A U −1 0 . For Tracer 25 B, the characteristic scale is set as its background concentration in the free troposphere ( B 0 ). For Tracer C, the characteristic scale is defined as C 0 = L 0 U −1 0 k A 0 B 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 A = A/ A 0 , c B = B/ B 0 and c C = C/ C 0 . With these definitions, the rates of the three tracers can be expressed as the following equations: where the dimensionless rate constants are Competing interests. The authors declare that they have no conflict of interests.

10
Acknowledgements. This work has been financially supported by Max Planck Institute for Meteorology as part of the doctoral thesis of CWYL (Li, 2019). The computation of the simulations presented in this work was completed on the Mistral supercomputer of DKRZ.
DKRZ also contributed to the storage of the data presented in this work.