Measurement-driven Large Eddy Simulation of dispersion in street canyons of variable building height Journal of Wind Engineering & Industrial Aerodynamics

The ventilation ef ﬁ ciency of three periodic building patterns of equal total volume and packing density is investigated: street canyons bounded by H uniform and 0.5 H /1.5 H variable height buildings, via wind tunnel measurements and Large Eddy Simulation. The numerical model utilizes the Transient Wind Forcing method to take the effect of eddies larger than the domain size into account, and a special Lagrangian dispersion model, which allows for the calculation of particle trajectories exceeding the periodic boundaries of the LES domain. The spatial and temporal characteristics of the concentration responses of both pulse-like and steady point sources located at the street surfaces are analyzed. It is shown that the variable building height has a favorable effect on urban ventilation in densely built areas: the average near-ground concentration can be reduced by up to 70% for the variable height buildings in a staggered arrangement. In terms of velocity, turbulence and concentration distributions, the model results are consistent with the experiments, verifying the applicability of the model for comparative air quality studies. Since the presented dispersion model is capable of handling dynamic changes in wind direction and magnitude, the accuracy of the TWF model can potentially overcome the limitations of wind tunnel tests.


Introduction
As urbanization advances, access to healthy air is becoming increasingly important in the lives of city dwellers. The individuals' exposure to traffic-related air pollutants and airborne pathogens is increasing due to the declining dilution associated with the concentration of pollutant sources. There are some feasible technical measures to limit transmission and promote ventilation, the effectiveness of which can be assessed by experimental or numerical simulation methods. The results of field measurements, wind tunnel tests, and microscale meteorological models are strongly related to the examined geometry as well as to the flow characteristics of the atmospheric boundary layer, which are also geometry dependent. The results are, therefore, difficult to generalize.
Periodically repeated patterns, such as building blocks or junctions, are key features of the urban canopy. By exploiting the periodicity, the number of model parameters can be significantly reduced, facilitating the generalization of the model results. Consequently, repetitive geometric patterns are preferred subjects of field experiments (Biltolft, 2001;Kanda, 2009;Chen et al., 2020), wind tunnel measurements (Leitl et al., 2007;Yee and Blitolft, 2004) and numerical models (Eichhorn and Balcz o, 2008;Santiago et al., 2010;Dejoan et al., 2010;R akai and Krist of, 2010). For the latter, the repetitive surface patterns can be easily constructed by applying periodic boundary conditions at the lateral sides of the simulation domain, thus taking into account an infinitely large repetitive structure in a small modeling box. A further advantage of this method is that it does not require the exact definition of inlet boundary conditions, which is particularly beneficial for CFD models utilizing Large Eddy Simulation (LES)a turbulence model that can produce more accurate results in urban flow and dispersion compared to the conventional RANS approaches (Tominaga and Stathopoulos, 2011;Tominaga et al., 2013). The inlet boundary conditions for LES may be defined according to the method proposed by Kataoka and Mizuno (2002); alternatively, the method proposed by Xie and Castro (2008) and Xie (2011) using velocity time series from field observations supplemented with artificially generated turbulence. Nevertheless, the inlet boundary conditions introduce additional model parameters and modeling uncertainties into the analysis. Instead of prescribing the inlet velocity and turbulence profiles, in the periodic models, the flow can be created by a spatially distributed driving force, the volume integral of which balances the drag force of the surface objects, such as buildings and vegetation.
Most periodic models are driven by a constant intensity volumetric driving force (Ciofalo, 1996;Su et al., 1998;Cheng et al., 2003;Baker et al., 2004;Cai et al., 2008;Letzel et al., 2008;Kikumoto and Ooka, 2012;Krist of and Füle, 2017) which requires the specification of only two model parameters: the horizontal components of the pressure gradient. The disadvantage of the static pressure gradient-driven LES models is the fact that the size of the resolved turbulent structures is limited by the size of the computational domain, which also limits the accuracy of dispersion results since the larger turbulent structures have a substantial contribution to the transport processes. The literature review by Tominaga and Stathopoulos (2016) also encourages the investigation of the influence of the atmospheric turbulence scales on dispersion.
The above-described problem can be solved by using a transient driving force (Krist of et al., 2020). Measurement data or the results of larger-scale meteorological models can be used to control the driving force in a way that the horizontal components of the calculated velocity follow the prescribed wind speed at a designated observation point.
An essential element of the Transient Wind Forcing model presented in our previous study (Krist of et al., 2020) is that the turbulent spectrum is divided into microscopic, mesoscopic, and macroscopic parts based on the size of the smallest and largest turbulent structures that can be resolved in the model. The microscopic turbulence is considered using a Subgrid-Scale Stress (SGS) model; the LES model resolves the mesoscopic turbulence, so it appears in the results in the form of velocity fluctuations over time. Furthermore, the effect of macroscopic turbulence is taken into account by the transient driving force, so the frequency of the propulsion must match the lowest frequency of the mesoscopic turbulence. The relaxation time of the propulsion control (Δt 0 ) must be chosen in proportion to the ratio of the domain length and the average wind speed measured at the reference point. The temporal resolution of the measured data set must be finer than the relaxation time of the force control.
At the time of publishing the TWF model, the authors were not aware of the Real-Time Boundary wind Condition (RTBC) model previously published by Zhang et al. (2011), which was, to our present knowledge, the earliest measurement-driven periodic LES model. Comparing the two methods, a significant difference is that the RTBC model does not fit the time parameter of the propulsion to the turbulence time scale resolved in the model; moreover, the resolution of the measurement data set is well below the control's cut-off frequency, which resulted in a series of steady flows. The most recent version of the RTBC model  is more elaborate in this regard, due to fitting the time parameter of the method to the frequency corresponding to the peak of the pre-multiplied power-spectrum of the measured turbulence.
There is also a crucial difference in the choice of the location of the gauging point: in the case of the RTBC model, the observation point is located near the roof of the buildings (z/H ¼ 1.1), so the gauge intermittently fell in the wake of the building; consequently, the driving wind speed does not adequately represent the wind speed above the urban canopy. Presumably, this caused the simulation to produce windless episodes often. The spatial distribution of the driving force is homogeneous horizontally for both the RTBC and the TWF models, and it follows a prescribed vertical profile. Using a driving force with a Gaussian profile, the TWF model produced average velocity and turbulence distributions in street canyons with a height-to-width ratio of H/W ¼ 2 and 3 that showed good agreement with nearly full-scale measurement data. The RTBC model uses a power-law distribution for the driving force in the vertical direction, which assumption is not yet validated.
In the time-dependent flow field created by the TWF model, the dispersion process is investigated by the Lagrange model, i.e., by determining the trajectories of particles. We track the position of massless particles within the domain and record the number of their streamwise and lateral periodic jumps, which can be used to calculate the position of the particles in a horizontally unlimited space. This method makes it possible to generalize the results of periodic dispersion models to an aperiodic dispersion process, such as the transmission of contaminants from a source canyon to subsequent downstream street canyons. Our previous study (Krist of et al., 2020) has proven that including the large-scale velocity fluctuations in the simulationinstead of working only with the mean wind characteristicscan significantly influence the resultant pollutant distribution, namely, a wider plume and a significant temporal variation of the pollutant concentration in the source canyon was observed. Similar results were shown previously by Xie (2011): when including 30 s resolution mesoscale wind data to the street-scale LES simulations, the predicted dispersion was in significantly better agreement with the field measurements than when steady wind conditions were applied.
Previous street canyon models usually employ lengthwise line sources of constant intensity (Meroney et al., 1996;Chan et al., 2002;Simo€ ens and Wallace, 2008;Gromke et al., 2016), which represent the steady emission at the traffic lanes. The emission intensity distribution of a spatially inhomogeneous line source, however, can be composed as the sum of the intensity of its sub-surfaces represented by point sources. A similar summation can be performed in time: the changing source intensity can be represented by the sum of intermittently active sources of various strengths.
To generalize the results, we analyze the effect of point sources, thus allowing for the investigation of spatially and temporally inhomogeneous emissions, such as exhaust gases leaving a multi-story car park, or even aerosolized pathogens exiting the human respiratory system. The advantage of the Lagrange method is that any number of point sources can be modeled simultaneously because the emission point of each particle is known, thus the effect of different source locations can be distinguished. Also, in wind tunnel models with heterogeneous roof height, it is necessary to test several point sources located along the street, in accordance with the geometrical changes.
The transmission between a source and an observation point is characterized by the local dose distribution, determined by the total residence time of the particles per unit volume. The dose distribution does not dependin terms of the ensemble meanon the duration of the emission, only on the total mass of the pollutants emitted, provided that the observation time is long enough to include the entire concentration wave (rise, plateau, and decay; see Fig. 1) caused by the emission. The resulting dose distributions can be used to evaluate the effect of either steady or pulse-like sources. For the latter, a representative concentration response can be determined by averaging the concentration waves generated by several gas-puffs.
For the validation of the numerical model, wind tunnel experiments were performed in the framework of this research, investigating the flow and dispersion processes in various parallel street canyon series subjected to perpendicular wind relative to their length axis. The recorded timeresolved two-component LDA data was also used for driving the TWF model. Based on the experimental results, we quantitatively evaluate the accuracy of the numerical model in terms of mean velocity, turbulence quantities, and average concentration distribution of tracer gas emitted from steady point sources of different locations.
The heterogeneity of the building height strongly influences the ventilation of street canyons by modifying the roof-level flow structures and mixing processes. The wind tunnel measurements of Pascheke et al. (2008) and the LES simulations of Boppana et al. (2010) showed that around staggered roughness elements of varying height (0.28H … 1.72H in multiple steps) the turbulence intensity was higher, and canopy average concentration (up to H) was lower than for a uniform block height. Based on field experiments, Kanda and Moriizumi (2009) concluded that for not too sparsely located buildings of two block heights (H 1 and H 2 ), the mass transfer coefficient characterizing the ventilation of the building arrangement can be approximated by the mass transfer coefficient of an array of uniform homogeneous building heights equal to the difference between the aforementioned block heights (H 1 -H 2 ). Nosek et al. (2016Nosek et al. ( , 2017 investigated the flow and dispersion processes in a network of street canyons representative for European inner cities with both uniform and heterogeneous roof heights. The wind tunnel experiments revealed that the building non-uniformity can improve the pollutant removal capability of complex building arrangements as well, in the case of numerous wind directions. The favorable effects of building height variability were also observed for typical metropolitan building configurations as well: strictly arranged high-rise buildings of different heights showed enhanced pollutant removing capabilities, see Hang et al. (2012), Chen et al. (2017) and the references therein.
As the building height heterogeneity may be a practical planning concept to improve urban ventilation since it does not adversely affect the economic parameters, the building volume per unit area, the building  packing density and the installation raster area are kept constant in our investigations (following Krist of and Füle, 2017). In this paper, we analyze the dose responses to street level pollutant emissions in long, parallel street canyons of constant and variable building height, using numerical models and wind tunnel experiments.

Wind tunnel measurements
Wind tunnel experiments were performed to investigate the flow and dispersion characteristics of three periodic building arrangements using the closed-circuit horizontal (G€ ottingen-type) wind tunnel of the K arm an Wind Tunnel Laboratory of the Department of Fluid Mechanics at the Budapest University of Technology and Economics. The wind tunnel has a circular cross-section of 2.6 m diameter at the open test section of 3.8 m length.

Wind tunnel setup and building patterns
For the wind tunnel measurements, three periodically repeated building configurations were constructed from styrofoam. The baseline case consisted of 23 small-scale (M ¼ 1:200), 1/1 (height-to-width, H/W) aspect ratio uniform street canyons, positioned perpendicular to the wind direction. The canyons were constructed using buildings of constant B Â H ¼ 100 mm Â 100 mm (breadth by height) rectangular cross-section, placed at W ¼ 100 mm distance one after another; hence, the offset between two adjacent streets is S ¼ 200 mm. The 22 canyonsthe length of which exceeded the width of the measurement sectionare numbered from À10 to þ11 in the streamwise direction, with tracer gas sources located in the 0th canyon.
In order to assess the effect of the heterogeneity of the buildings on air quality, two more building arrangements were designed. Towers were constructed by elevating and lowering the rooftops in the middle section of 1.25 m width, in 6 segments by 0.5H. Hence, in each row, three taller and three shorter segments were created, each of B Â T ¼ 100 mm Â 208.33 mm plan area, see Fig. 2 and Fig. 3. In the matrix arrangement, the tall towers are located behind another tall tower in streamwise (x) direction; therefore, the roof height is varying only in the lateral (y) direction. Finally, in the staggered arrangement, the tall towers are shifted laterally relative to the previous and following rows, i.e., the heterogeneity of the roof height is present in both x and y directions. Note that all three building configurations have identical total volume.
The building models were subjected to a uniform approach flow. 4.75H upstream of the test section, a 0.6H tall baffle was placed (see Fig. 2), in order to reduce the size of the separation bubble forming over the first few rows of buildings. It was found that after the so-called 0th canyon (physically the 11th in streamwise direction), the flow can be considered periodic. Two identical tracer gas sources of 28 mm diameter were flush-mounted in the base plate of the building models. As shown in Fig. 3, the sources are located in the middle of the 0th canyon in streamwise direction (x ¼ 0) and with a lateral offset of half a tower width (y/T ¼ AE0.5).

Velocity measurements
The distribution of the horizontal velocity components for each building pattern was measured using two-component Laser Doppler Anemometry in the 0th and 4th canyons. The complete measurement instrumentation is listed in Table 1. The optics of the LDA device was mounted on the probe traversing system of the wind tunnel to enable precise and repeatable positioning. The system incorporates Bragg-cells for both velocity components, so the sign of the velocity components could be identified. The pass-through time of each particle signal was recorded and further used as weights when calculating the average velocity (gate-time weighting).
The velocity measurement in each gauging point took 150 s with a sampling frequency range of 100 … 1000þ Hz (depending on how efficiently the seeding is able to access the gauging points), resulting in time  statistics of the horizontal velocity components with less than 5% statistical uncertainty in the streamwise direction, except for the points with extremely large turbulence intensity values. In the reference points, the representativeness uncertainty values are 0.57% for the uniform canyons, 0.56% for the matrix towers, and 1.06% for the staggered towers. The representativeness was calculated from the turbulent time scale, according to Kaimal and Finnigan (1994). The absolute measurement uncertainty is 0.1 m/s based on reproducibility measurements. The ambient temperature and pressure, as well as the flow velocity outside the boundary layer, were monitored continuously during the experiment, and the LDA results were later scaled according to the minor variation of the wind tunnel bulk velocity. The average relative deviation of the mean streamwise velocities originating in the different measurement locations (either in the 0th or in the 4th canyons, calculated as the average of ju 0th À u 4th j=u 0 ) was 4.4% in the 0 … 2H range (see Fig. 4), which is the focal area of the present study. The Reynolds number based on the reference mean velocity at 2H and the roof height H was calculated around Re H ¼ 31 000, which exceeds the threshold of 12 000, above which the flow can be considered Reynoldsinsensitive for H/W ¼ 1 street canyons. Moreover, it is located in the range between Re H ¼ 31 000 … 58 000, for which the flow can also be considered Re-independent for H/W ¼ 2 street canyons (Chew et al., 2018). The Reynolds-independence study was performed at four different heights for each building pattern in the range of Re H ¼ 12 000 … 45 000; the results are shown in Fig. 4. It can be concluded that the flow field can be considered Reynolds number independent above Re H ¼ 30 000.
Apart from logging the time statistics, full-length time series were also recorded for driving the numerical model, at z 0 /H ¼ 2 height over the left source. The measured velocity time series were later resampled to the simulation time steps for the propulsion of the CFD models (for the details, see Section 2.2.2).

Concentration measurements
The pollutant distribution of the building arrangements in the wind tunnel was mapped using Fast Flame Ionization Detection (Fast FID). The applied tracer gas representing the traffic-induced air pollutants was pure methane (100% CH 4 ), which was continuously emitted from one of the two sources at a time. To control the volume flow rate, pre-calibrated smart mass flow meters were applied (see the details in Table 1).
To take the characteristic source locations into accountin other words, whether the emission is placed after a short or a tall building, see Fig. 3 -, the measurements were performed for both scenarios in the matrix and staggered tower arrangements as well, and only once for the uniform canyons, with only the left source being active. The resultant dispersion fields were sampled for 150 s, in each measurement point, with a response time of 2.05 ms, resulting in 12% relative and (0.28 ⋅ 1 c*) absolute uncertainty, based on 4 reproduced measurements (Fig. 5). The relative uncertainty is calculated as the 90th percentile of the coefficient of variation (relative standard deviation), and the absolute uncertainty is obtained as the mean of the coefficient of variation in the far-field. The FID device was calibrated every 30 min during the measurements using two gas mixtures of known concentration; moreover, the results were compensated for the minor temporal changes in the bulk velocity and in the background concentration of the wind tunnel as well.
The measured concentrations were normalized according to  (11)). When comparing single source between experiment and simulation, A was kept constant, equal to the plan area of the large simulation domain (A ¼ 4ST). Note that in the wind tunnel experiment, the sources were in operation long enough before the concentration sampling began for the dispersion field to reach a steady-state, and they remained active during the entire length of the measurement; hence, for the present experiments, t s ¼ t avg .

CFD simulations
The Transient Wind Forcing model, which serves as the basis of the present CFD simulations, was developed by Krist of et al. (2020). In the framework of this study, it is implemented in the Finite Volume Method-based multi-purpose CFD software, ANSYS Fluent 19.3, using a three-dimensional transient modeling approach, similar to our previous investigations.

Numerical meshes and boundary conditions
For all of the three periodic building patterns (uniform canyons, matrix towers, staggered towers), the elementary repeated geometrical unit was simulated, as shown in Figs. 3 and 6. The length (X) and width (Y) of the domain is defined according to the size of the repeated elementary geometry, and the height of the domain is determined based on the maximum roof height and the length of the domain in the dominant wind direction: An evenly spaced hexahedral mesh was constructed to discretize the computational domain. The spatial resolution was Δx ¼ Δy ¼ Δz ¼ W/32, i.e., the street width was divided into 32 cells in all three cases. For the assessment of the mesh convergence, two additional meshes were created for the uniform canyons case of W/22 and W/48 resolution. The wall boundary layer is not directly resolved, since previous studies such as Xie and Castro (2006) have shown, that the mass exchange is principally governed by the large eddies, and the near-wall mesh refinement cannot significantly improve the results. Hence, even more recent studies such as Li et al. (2008), Liu et al. (2012), Michioka et al. (2014), Llaguno-Munitxa et al. (2017) and Castro et al. (2017) prefer to employ evenly spaced computational grids. Xie and Castro (2006) also conclude that resolving the characteristic size by 16 elements can be sufficient for capturing the building-scale flow structures and the surface drag using LES.
No-slip walls were applied on the ground and the surface of the buildings. Periodic boundary conditions were assumed at the horizontal boundaries of the domain, both in the x and in the y directions. At the top of the domain, a symmetry boundary condition was defined. The details of the simulation cases are compiled in Table 2.

Measurement-based driving force and turbulence modeling
The main driving force of the flow is the turbulent momentum exchange with the higher atmospheric layers, which is taken into account in the present numerical model as a time-dependent volume source of momentum, with constant intensity in both horizontal directions, i.e., it is only dependent on the vertical coordinate and on time. The driving force represents the macroscopic turbulence; in other words, large flow structures characterized by a minimum size, which was assumed equal to the length of the domain in flow direction. It was concluded based on the spectral energy density distributions of the wind tunnel measurements, that around 70% of the turbulent kinetic energy is contained in the macroscopic fluctuationsi.e. below the frequency defined by 1/Δt 0at the center of the propulsion.
The time scale of the macroscopic turbulence, which is proportional to the average flow-through time, is calculated according to in which X [m] is the length of the domain in the dominant flow direction, and C t [-] is the time parameter of the TWF method. In the present study, C t ¼ 0.5 was used.
To impose large-scale turbulence in the TWF model, the time series of the horizontal velocity components were recorded in the wind tunnel during the measurements for each building pattern. The reference location is at 2H height over the left pollutant source (x 0 ¼ 0, y 0 /T ¼ 0.5, z 0 / H ¼ 2) in all cases.
The time series are employed for the propulsion of the model using volume sources of momentum. The driving force is controlled in a way that the wind velocity in the numerical model follows the reference velocity time series with Δt 0 relaxation time. The propulsion, therefore, shows a low pass filter-like spectral behavior, with a À3 dB amplitude drop corresponding to the cut-off frequency f c ¼ 1/Δt 0 . C t < 1 was chosen to compensate for this phenomenon.
Imposing the reference time series at z 0 /H ¼ 2 is effectively an immersed boundary condition. Reproducing the velocity field above the reference height is out of the scope of the TWF approach; however, to allow free eddy motion, the simulation domain is kept taller. Between z/ H ¼ {2.2, 4.2} the flow in the CFD model is 16.3% slower than in the wind tunnel. Note that regarding the 0th … 4th canyons, only 0.2%, 10.9% and 10.8% of the particles were transported above the reference heighttherefore affected by the lower velocityin the case of the uniform canyons, matrix towers and staggered towers, respectively.
When moving away vertically from the reference location, the intensity of the driving force may decrease. The vertical distribution of the intensity of the wind forcing should be chosen according to the coherence of the flow. In the present study, the vertical propulsion distribution is described by a Gaussian profile with a constant radius of where C L is the length parameter of the TWF method. In the present study, C L ¼ 0.5 was used, similarly to the value of C t . The spatial and temporal variation of the propulsion in the TWF modelwhich is representing the macroscopic flow structures exceeding the domain sizetherefore can be expressed as the product of the timedependent velocity control a(t), and the spatial distribution of the driving force G(z), illustrated in Fig. 6.    If the relaxation time τ is too large, the velocities produced by the CFD model cannot follow the reference time series properly. On the other hand, the application of an appropriately long relaxation time allows for the applied Large Eddy Simulation turbulence model to create mesoscopic turbulence around the reference point, preventing the wind forcing from acting as a direct velocity constraint. An exception to this latter condition is the initial period of Δt 0 duration, in which a shorter relaxation time is used to speed up the initialization of the flow field. The relaxation time is increased exponentially within the startup transient, which leads to the following formula: in which Δt [s] is the time step size of the simulation; moreover, t [s] is the physical time, and t start [s] is the starting time of the simulation. The time step size of the simulation was constant in each case, and it was chosen to satisfy the Courant-Friedrichs-Lewy stability condition, i.e., the Courant number (C ¼ u ⋅ Δt/Δx) is less than one in the entire domain at all times.
The mesoscopic turbulence, i.e., the temporal velocity fluctuations were modeled using Large Eddy Simulation, and the effects of microscopic turbulence were taken into account using the Smagorinsky-Lilly sub-grid-scale stress model (Smagorinsky, 1963;Lilly, 1992), using C s ¼ 0.1, recommended by Shah (1998) for a flow past a blunt obstacle.

Pollutant transport model
The dispersion of the traffic-induced air pollutants was investigated using a Lagrangian discrete phase model in the ANSYS Fluent system. Massless, inert particles were emitted at the location of the two circular tracer gas sources placed in the wind tunnel model (at x ¼ 0, y/T ¼ AE0.5, z/H ¼ 0.05). The number of times each particle had already crossed the periodic boundaries was recorded in two user-defined variables per particle. This way, based on the number of previous jumps and the current position of the particles within the domain, the full non-periodic 3D trajectories can be reconstructed in the horizontally unbounded space at any time instance.
The tracer gas emissions were modeled using two approaches. Firstly, simple pointwise particle injections were used, with a constant emission rate of one particle per time step, which is a discrete representation of the continuous, spatially distributed sources featured in the wind tunnel experiments. Secondly, to investigate the finite source area and the temporal discretization effects, each of the pollutant sources was also modeled using 1000 particles distributed evenly on a circular plate, with their center at the above-defined locations, and their diameter identical to the tracer gas emission of the wind tunnel (28 mm).
The spatially distributed particle clouds were emitted at 2, 3, and 4 s of flow time, at the beginning and the end of the simulation time step. It was revealed that the large-scale dispersion of the pollutants, i.e., the temporal variation of the particles in the subsequent street canyons, is practically independent of the diameter of the source (tested for 2, 10, and 28 mm). The first 2 s of the simulations were used for the initialization of the turbulent flow field, and the continuous particle emission was enabled only after this period. In addition, at the end of the simulations another 8 s and 2 s was run with the sources switched off (for the uniform canyons and for the matrix and staggered towers, respectively), for the vast majority of the particles to leave the investigated area (see the results later in Fig. 13a).
The transmission relationship of an arbitrarily chosen locationfor example, a cell of the numerical meshand one particle source can be specified by the local dose, given by in the above equation, d [s/m 3 ] is the dose, calculated for each cell c and particle p of the simulation domain in each particle time step s. N s [-] is the number of particle time steps in the simulation up to the current physical time, N p [-] is the number of particles, Δt p,s [s] is the particle time step size corresponding to the pth particle in the sth particle time step, and V c [m 3 ] is the volume of the cell. The Inside(c,p,s) function takes the value of one if particle p is present in cell c at the end of the time step s; otherwise, it yields zero. It is an important aspect of the definition of the dose thatas a consequence of the temporal summationthe particle emission rate does not need to be constant in time.
It can be seen that the calculation of the dose can be realized as the incrementation of a user-defined scalar per cell in each particle time step s by p ⋅ Δt p,s /V c . The dose distribution was calculated using several User Defined Memory (UDM) fields covering the periodic building patterns. The number of UDMs required for the separate representation of the dose fields corresponding to different sources can be calculated by multiplying the number of investigated rasters by the number of particle sources.
The average normalized concentration distribution can be calculated based on the dose field following Eq. (8).
Note that since one particle is emitted from the steady sources per time step, N p ¼ t s /Δt for a single point source. Comparing the above equation with the first definition of the dimensionless concentration given in Eq. 1, it can be shown that the time-averaged concentration can be obtained from the dose as Because of the limited number of particles, the concentration field computed from the dose results of the Lagrange model is noisy; therefore, it needs to be smoothened in space (Fig. 7). Therefore, a Gaussian diffusion filtering procedure was applied to the scalar fields. After the periodic boundary conditions are changed to symmetries, and the solution of the flow equations are switched off, the generic transport equation with zero velocity is solved on the user scalar fields representing the concentration distribution, in each of the rasters separately. The spatial filter size was chosen based on the average mean particle count under roof height. The average distance between the particles is in which d avg [s/m 3 ] is the average dose within the canopy (z H max ). The diffusion coefficient of the smoothing process (D [Pa•s]) was chosen according to Fig. 7. Dose fields in one raster of the staggered towers case before (left) and after (right) the spatial filtering. Fig. 8. Comparison of the pointwise velocity and turbulence results between the experiment and the CFD simulations for the three building patterns: uniform canyons (left), matrix towers (middle), and staggered towers (right).
in which t diff [s] is the temporal length of the filtering, which was kept constant (10 s) in all cases.
Based on the near-ground or canopy average concentration, the mass Stanton number characterizing the ventilation efficiency of the different building patterns can be obtained as the reciprocal of the mean normalized concentration:

Solution methods and convergence criteria
For the numerical solution, the Bounded Central Differencing Scheme flux formulation and the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm were used, along with the Bounded Second Order Implicit scheme for temporal discretization, which are the default options for LES simulations in ANSYS Fluent 19.3 (ANSYS, 2019). The iterative solution was carried on in each time step until the residuals of the flow equations were reduced by at least three orders of magnitude compared to the initial state but in a maximum of 20 iterations. For the unsteady particle tracking in Fluent, the trapezoidal method was utilized for numerical integration. To model the microscale turbulent diffusion, the Discrete Random Walk model was enabled for the particles from all sources, with the time scale constant uniformly set to 0.15 (default value). In the scalar transport equationused for the spatial filtering of the concentration field -, the second order flux formulation was used.

Mesh convergence and the impact of the domain size
The mesh convergence analysis was executed on the model of uniform canyons using the small domain (see Fig. 6 and Table 2 for the details). As Large Eddy Simulation is applied, the spatial resolution affectsbeyond the discretization errorthe width of the resolved turbulent spectrum, i.e. the boundary between the meso-and microscopic turbulence.
Lower velocities were found within the canyon in the case of the finest mesh (W/Δx ¼ 48), which showed a closer agreement with the measurement data compared to the coarser meshes (W/Δx ¼ 32 and 22). This can be a consequence of the fact that the relatively thin shear layer at roof height (z/H ¼ 1) is sensitive of the resolution. However, it must be noted that the variation of the numerical results obtained on different meshes falls into the range of the difference between the experimental profiles obtained in different canyons (see Figs. 8 and 4, respectively). The turbulence intensity results did not improve substantially on the finer meshes, which points out that most of the inaccuracies can be attributed to other model parameters.
In terms of the mean concentration distributions, the application of the larger simulation domains (wider and taller domains with the corresponding propulsion parameters, see Fig. 6 and Table 2) substantially improve the agreement with the experimental data, especially in the source canyon of the matrix towers arrangement.

Validation for velocity and turbulence
The time average and the fluctuation of the streamwise and lateral velocity components were compared (wind tunnel experiments vs. CFD simulations) for all three building patterns in the 0th canyon. For the matrix and staggered tower arrangements, the comparison was carried out behind both a tall and a short building (x ¼ 0; y/T ¼ AE0.5). The propulsion side profiles (y/T ¼ 0.5) are indicated with darker, while the ones on the passive side are denoted by lighter colors in Fig. 8. The profiles obtained in the smaller domains are shown in greyscale. For each profile, 20 gauging points were used, spaced evenly in the vertical direction (z/H ¼ 0.1, 0.2, …, 1.9, 2). Note that the horizontal location of the vertical profiles is identical to that of the sources shown in Fig. 3.
It is shown in Fig. 8 that the normalized mean velocity profiles show good agreement between the experimental and simulation results, characterized by a 95% correlation. The largest differences can be found in the streamwise mean velocity component over roof height in the case of the towers (z/H > 1.5), which is a consequence of the horizontally homogeneous driving force being relatively close to the top of the buildings.
The turbulence results of the CFD model display a good qualitative agreement with the wind tunnel data, which can also be seen in Fig. 8. The correlation between the computed and measured turbulence data is 83%. It is shown by both methods that the turbulence intensity is significantly increased below roof height by constructing towers, which suggests that the variable roof height promotes canopy ventilation. There is a remarkable turbulence deficit at the reference location (z/H ¼ 2) which could be mitigated by reducing the time parameter of the model (C t ).
Another source for inaccuracy is that the periodic LES model is only able to generate mesoscopic turbulence caused by the wind shear above the roofs. However, in the wind tunnel, a strong shear layer is formed at the leading edge of the building array (upstream from the À10th canyon), which generates turbulence that can be observed even in the gauging points at the 0th canyon. The existence of such an inner boundary layer was also observed in other wind tunnel experiments, see Castro et al. (2017). The turbulence deficit is less pronounced for the matrix and staggered tower configurations, since more kinetic energy (TKE) is produced by the canopy, suppressing the turbulence originating from the leading edge of the building array.

Validation for concentration
The lateral concentration distributions are shown in Fig. 9, and the vertical profiles are shown in Fig. 10. Generally, the model captures the important tendencies of the measured mean concentration distribution. The correlation coefficient between the measured and computed concentration data is 0.77.
The best match can be found within the 0th canyons, except for the close vicinity of the sources, which can be explained by the fact that the wind tunnel measurements feature a 28 mm diameter flush-mounted sources, while the particles in the CFD model are injected at a single point. It was found in the framework of the present study that the temporal variation of the number of particles in separate street canyons is practically independent of the diameter of the source (tested for 2, 10, and 28 mm).
The mean concentration field is generally overestimated within the downstream canyons; furthermore, the concentration boundary layer over the buildings is thinner in the CFD simulations compared to that observed in the wind tunnel measurement, which is the most pronounced in the case of the uniform canyons. The differences might be attributed to the following reasons: 1. In some cases, the turbulent kinetic energy is underestimated, which affects the dispersion as well. minor crosswise flow within the canyons, which is not present in the CFD model (see Fig. 8), causing the additional dilution of the tracer gas. 4. In the low concentration zones, the relative measurement error is higher.

Comparison of the building patterns
3.4.1. Dispersion fields of different sources In each street canyon, uniform or not, as a consequence of the canyon vortex (Fig. 11a) or similar rotating flow structures, the locally emitted pollutants are swept toward the leeward wall, resulting in an imbalance between the pedestrian exposure at the leeward and windward sides. This asymmetry in concentration is not present in the canyons downstream from the sources.
The introduction of building height variability generates lateral velocity as well as increased turbulence under roof height, therefore promoting the dilution of the pollutants. For both tower structures, there is a downwash at the windward side of the tall buildings, which is deflected in both crosswise (y) directions as it reaches ground level. Therefore, diverging streamlines can be observed in front of the tall buildings, thus the pollutants emitted at these locations are spread at street level (Fig. 11cd). On the other hand, as a consequence of the converging streamlines upstream from the short buildings, the pollutants are locally accumulated (Fig. 11be). Regardless of the position of the source; however, in the case of variable roof height, the exposure of the pedestrians and the residents of the nearby buildings are significantly smaller than that of the uniform canyons, and negligible in the subsequent streets according to the model results.

Characteristic vertical profiles
In this section, the periodic building patterns are globally characterized by the vertical profiles of velocity, turbulence and concentration, obtained by horizontal averaging in one cell thick layers (Δz ¼ H/32).
In Fig. 12 similar near-ground velocity ratios (u/u 0 , v/u 0 ) can be observed for all periodic building arrangements, but there is significantly more turbulence both in streamwise and in lateral directions when towers are present. The variable roof height results in rougher canopies, which cause smaller velocity gradients at z/H ¼ 1.5 (compared to z/H for the uniform canyons), as well as substantially higher turbulence in the boundary layer. Rough building patterns are prone to extract more energy from the atmospheric boundary layer, which leads to a decrease in wind speed in the long run. This effect is not addressed in the present study, in contrast to our previous investigations (Krist of and Füle, 2017).
The Lagrangian approach applied in the present model allows for analyzing the entire dose field of each source, as the particles are tracked in the infinite space. It is advantageous that the individual vertical concentration profiles of each investigated raster can be summed to obtain the total footprint of a source, which characterizes the canopy ventilation efficiency. In the present study, for the matrix and staggered towers, two injections are placed in one raster in accordance with the variable roof height; therefore, their combined result was divided by two.
The vertical profiles of mean normalized concentration shown in Fig. 12 confirm the superiority of the matrix and the staggered tower arrangements over the uniform street canyons. By constructing buildings of different height while keeping the building volume constant, both the near-ground and canopy average concentrations can be significantly decreased, resulting in substantially more efficient ventilation. The dilution coefficients of the building patterns with variable roof height are

Time-dependent particle results
The location of the particles emitted from each source was recorded in every 0.05 s, which served as a basis for calculating the time-dependent particle counts in the first few canyons. The investigated volumes for all three building patterns can be defined by x/S 2 [i -0.25, i þ 0.25] and z/ H 2 [0, 1]identical to the below-roof part of the uniform canyonswith i being the number of the canyon of interest.
As discussed in Section 2.2.3, there are two types of sources placed in the domains: continuous point-like emissions, injecting one particle per time step; and spatially distributed sources of 28 mm diameter, injecting 1000 particles simultaneously, at six different times (2, 3, and 4 s of flow time, at the beginning and the end of the time step). Fig. 13 shows the particle count relative to the total number of particles emitted from either the steady or the pulse-like sources during the entire simulation (N all,steady and N all,pulse , respectively). Each of the curves of Fig. 13b shows the average of the six particle clouds, with the flow time relative to the time of emission (t 0 ).
It can be seen in Fig. 13a, that the average number of particles accumulating in the source canyons are substantially lower for the matrix (left source: À43%, right source: À55%) and the staggered towers (left: À62%, right: À80%) compared to the uniform canyons. Similar tendencies can be observed in the downstream canyons; however, it must be noted that the mean particle count in the subsequent streets is one order of magnitude smaller. In this regard, it does not seem crucial that the extreme values corresponding to the staggered, and especially, the matrix towers can exceed the extremities of the uniform canyons. The curves also reveal that far from the source canyon, the particle clouds emitted from the left and right sources show increasing correlation; i.e., the point of origin has a decreasing influenceas expected based on the observations of Xie and Castro (2009). Fig. 13b shows how particles emitted from pulse-like sources (such as pollutants from a bus stop, or aerosols from a sneeze) escape the canopy. It can be concluded that it takes less time for the particles in the matrix and especially in the staggered arrangement to clear the 0th canyon compared to the uniform canyons. Similar tendencies can be observed in the downstream streets, and the rates of the pollutants escaping the 0th canyon after switching off the continuous sources also show an analogous behavior in Fig. 13a. The particle reflux into the subsequent canyons is lower for the uniform canyons relative to the towers, resulting in lower maxima of the corresponding curves in Fig. 13b. But, again, the pedestrian exposure in the downstream streets is at least one order of magnitude lower, thus this effect is less important.
The particle count in a control volume can be interpreted as concentration, and the integral of the curves give the corresponding dose values. The area below the corresponding curves between Figs. 13a and b is practically identical for the entire length of the simulation, supporting the approach that the dispersion field of spatially and temporally inhomogeneous emissions can be constructed via the convolution of the concentration responses of the individual sources.

Conclusions and outlook
The technical objective of the present study is the assessment of the ventilation efficiency of three periodic building patterns of equal total volume: uniform H/W ¼ 1 street canyons and a combination of 0.5H and 1.5H high towers, both in matrix, and in staggered arrangement. The investigation was carried out via wind tunnel experiments and numerical modeling, namely, the Transient Wind Forcing model introduced by the previous work of the authors (Krist of et al., 2020). The TWF model utilizes periodic boundary conditions to simulate the flow field around the periodically repeated buildings. It is an important feature of the model that due to the time-dependent measurement-based propulsion of the flow, the macroscale turbulencein other words, the flow structures exceeding the domain sizecan also be taken into account. It was shown that at the center of the propulsion around 70% of the TKE is contained in these large-scale eddies. The TWF model utilizes Lagrangian particle tracking to simulate the dispersion of air pollutants emitted from near-ground sources. The number of periodic jumps in the x and in the y direction is recorded for each particle; therefore, their path can be followed outside the borders of the original periodic domain. Based on the local residence time of the particles, the pollutant dose distribution for both steady and pulse-like emissions can be calculated in the horizontally unlimited space. The present investigation is also aimed at the validation of the TWF-based aperiodic dispersion model.
Note that the accuracy of TWF model results is also affected by the correct choice of the model parameters (C t , C L , z 0 , X, Y, Z). The present research did not cover the full optimization of these model parameters due to the high computational demand of LES calculations.
During the validation, it was shown that the velocity and turbulence intensity profiles characterizing the flow field are reproduced with good accuracy (correlation coefficients: 95% and 83%, respectively). Based on the periodic flow fields, the aperiodic dose distributioncaused by pointlike sourceswas calculated. The mean normalized concentration distribution also showed good correspondence with the measurement data, especially in the source canyons, which are the most critical locations in terms of pedestrian exposure to traffic-induced air pollutants. Reasonable agreement was found between the experimental and simulation results, characterized by a 77% correlation. Furthermore, the time-dependent particle counts within the street canyons was analyzed, which indicated the superiority of the matrix and staggered towers compared to the uniform canyons, both for steady and for pulse-like sources. It was observed for both the matrix and staggered tower structures that the particles are able to escape the source canyon significantly faster compared to the uniform canyons; therefore, the timeaverage of the normalized particle count below average roof height is reasonably lower in the source canyon (matrix towers: 50%; staggered towers: 70%). Since the cumulative dispersion field of spatially and temporally inhomogeneous emissions can be constructed via the convolution of the concentration responses of point sources, similar to the ones presented in this paper, they can be crucial input data for operational forecast systems.
It was shown by normalized concentration results of the CFD model that the variation of the roof height can significantly improve the ventilation of the periodic building configurations. It was concluded, the near-ground and the canopy average concentrations show a 48 … 70% decrease for the matrix tower arrangement, and a 60 … 74% decrease for the staggered tower arrangement compared to the uniform canyons, while the total building volume is kept constant between the three cases.
The presented dispersion model is capable of taking the dynamic changes in wind direction and magnitude into account; therefore, the accuracy of the TWF model can potentially exceed the limitations of wind tunnel tests.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.