1 Introduction

More than half of the world’s population lives in urban settlements and the process of urbanization is estimated to progress further in the upcoming decades (United Nations 2018). Air quality in a city context has a critical impact on human health: 7.6% of total global deaths are attributable to air pollution and the ambient PM\(_{2.5}\) produced by combustion processes (including by motor vehicles) was the fifth-ranking mortality risk factor in 2015 (Cohen et al. 2017).

The complex topic of air quality has been extensively studied in recent decades. Vardoulakis et al. (2003) examined models for pollutant dispersion and concentration, while Britter and Hanna (2003) reviewed studies on airflow above and through urban areas with complex structures. Fernando et al. (2010) reported the advancement in research of city flows and air quality improvement. They highlighted that high-resolved numerical simulations should be exploited for improving the mesoscale and multi-scale models. Di Sabatino et al. (2018) discussed the recent advancement of pollution distribution in cities, while Buccolieri and Hang (2019) proposed a collection of studies on ventilation and pollutant dispersion. It is worth noticing that a large part of the literature on this topic studies the effects of the presence of trees or green infrastructure in an urban context (Yazid et al. 2014), but very few focus on building roof morphology and how roof infrastructure can increase natural city ventilation. The present contribution extends the investigation in this direction by analyzing in detail the paradigmatic case of the urban canyon.

Cities show a large range of possible configurations; nevertheless, urban canyons represent a typical basic unit of the city fabric. In this context, the major pollutant source is automobile exhaust gas at the street level, transported within the canyon by one or multiple principal vortices. The main cleaning mechanism is the removal at the building roof level by turbulent mixing and dispersion by atmospheric ventilation. The landmark work of Oke (1988) classified flow regimes with respect to the canyon aspect ratio (the ratio between canyon height and width), describing the overall structures driving the flow dynamics. Subsequent studies addressed the problem of air circulation in city canyons both via measurements and numerical simulations.

Reduced-scale experiments are widely used because they allow investigation of complex building geometries and relative high-velocity flow. Baik et al. (2000) studied the influence of aspect ratio in the development of the turbulent field and the primary canyon vortex. Several configurations were analyzed also considering buildings with different relative height. Brown et al. (2000) presented a measurement dataset of flow past a two-dimensional square building array. The development of the flow along multiple canyons was analyzed: a stable mean flow configuration was reached at the 4th canyon in the streamwise direction, while stable turbulent statistics were recovered at the 7th canyon. Louka et al. (2000) conducted a field study to investigate Kelvin–Helmholtz instabilities developing atop buildings. They found that such instabilities caused an unsteady between-building circulation, which induced episodic pollutant removal from the street canyon. A parametrization of such a phenomenon for dispersion models was given for low wind speeds. Louka et al. (2000) and Li et al. (2007) investigated a series of two-dimensional canyons with different aspect ratios in a water channel. To obtain reliable results, they found that the minimum channel height should be double the building height and that a short canyon length can generate three-dimensional motions by injecting flow in the street-parallel direction. Salizzoni et al. (2011) studied the interaction of two-dimensional canyons and surrounding ambient flow: the momentum transfer was mainly due to coupling between the turbulent structures in the outer flow and those in the shear-layer interface at the roof level. The latter was significantly influenced by turbulence in the external flow. Di Bernardino et al. (2015) undertook a water-channel experiment on two-dimensional canyons with aspect ratio 1 and 2. They showed that the in-canyon circulation strongly depends on the aspect ratio, in contrast with the outer flow, which is less sensitive to this. In subsequent work, Di Bernardino et al. (2018) analyzed the oscillation of the shear-layer, identifying two main frequencies: the largest matched the main vortex shedding frequency and the lower was related to a characteristic wall time estimated utilizing the friction velocity. This shear-layer flapping modulated the pollutant exchange rate. Chew et al. (2018) analyzed the Reynolds number (Re) independence criterion, which states that above a critical value \({ Re}_c = 11{,}000\) the flow non-dimensional statistics are invariant. They proved that this value is not valid for aspect ratios greater than 1.5, and additional threshold values were given for higher ratios.

The development of numerical simulation techniques improved the knowledge of the basic mechanisms underlying pollutant removal, allowing a detailed analysis of the fluid dynamics. Early studies employed Reynolds-averaged Navier–Stokes (RANS) simulations taking advantage of robust algorithms and low computational cost. Subsequently, the large-eddy simulation (LES) technique is preferred because of the higher accuracy and reliability, mainly due to the reproduction of unsteady fluctuations and transitional flows (Tominaga and Stathopoulos 2011). Large-eddy simulation directly resolves the larger scales of motion (with respect to computational cell size) and model the subgrid scales, leading to three-dimensional and time-varying flow simulations. The subgrid scales mainly dissipate turbulence kinetic energy and exhibit universal characteristics and thus can be modelled more effectively (Rodi et al. 2013). When the size of the computational cells becomes of the order of the Kolmogorov’s scale length, the LES reduces to direct numerical simulations. Because a larger part of the physical processes is reproduced, the resulting simulations are more realistic compared to RANS simulations and allow an accurate analysis of transient phenomena, in particular turbulent mechanisms. In this context, it is worth noting that a proper application of LES requires a careful definition of spatial resolution: the convergence of LES results is a complex issue because subgrid-scale (SGS) modelling and numerical errors are connected in a complicated way (Klein 2005; Celik et al. 2005; Klein et al. 2008). Concerning the present topic, LES studies have mostly dealt with a simplified geometry of two-dimensional square canyon: Liu and Barth (2002) simulated passive scalar transport using the Smagorinsky SGS model and wall functions (Smagorinsky 1963). They found that the \(97\%\) of pollutant is retained in the canyon during the simulated period, and that turbulent diffusion was the prevailing mechanism for scalar removal. The location of the pollution source affected the concentration in the canyon; a source shifted near the downwind wall produced lower concentration at the upwind wall. Walton and Cheng (2002) used a more accurate dynamic SGS model that produced a higher turbulence intensity in the core region of the canyon with respect to the RANS simulation with the \(k-\epsilon \) model, leading to a more homogeneous distribution of pollutant. The criterion to properly simulate the coherent structures that develop above a canyon was investigated by Kanda et al. (2004), who suggested a streamwise extension of the computational domain of at least ten times the canyon height. Cui et al. (2004) adopted a Smagorinsky model where the Smagorinsky constant \(c_s\) assumes different values inside and outside the canyon to account for the different turbulence content of the internal/external canyon flow. They noticed that weak ejection events were more frequent at the roof level, but a few strong sweep events contributed much to the momentum transfer. Li et al. (2008) performed LES with a one-equation SGS model and investigated the flow field and mechanisms of pollutant removal in city canyons with different aspect ratios. A similar SGS model was used by Letzel et al. (2008) to analyze the shear-layer intermittency. This was interpreted as a superimposition of a weak plane wake on a dominant plane mixing layer. Cheng and Liu (2011) also used a one-equation closure model and wall function to reproduce three street canyons. Most of the pollutant removal by external flow took place at the downwind wall, where the outflow impinged upon the building, but it also caused the pollution reinjection by entrainment. Michioka et al. (2011) employed a dynamic Smagorinsky model and compared the simulation results with wind-tunnel experiments. A wide range of aspect ratios was discussed. Pollutant removal events were enhanced by the ejection of low-momentum fluid when small-scale coherent structures appeared atop the canyon. Instead, the large-scale coherent structures developing above the canyon did not influence removal. In subsequent work, Michioka et al. (2014) conducted a numerical experiment in a variety of block configurations reproducing an urban-like scenario. Overall, they found that at least \(75\%\) of the total emissions from three-dimensional street canyons were accountable to turbulent motions.

In previous studies on urban canyons, attention was focused on the analysis of removal mechanisms when the building roof is perfectly smooth, in the absence of obstacles or roughness which can greatly affect the canyon–atmosphere exchange. It has been shown that the larger part of removal from within the canyon is due to turbulent events, particularly ejections at the roof level (Michioka et al. 2014). However, the higher turbulence intensity (i.e., higher turbulence kinetic energy) is localized near the downwind building, while turbulent motions are weaker in proximity to the upwind building roof. The basic idea of the present work is to use solid infrastructure to increase the natural turbulence on the upwind side of the canyon and, thus, produce a more effective pollutant removal. To this end, solid rectangular obstacles are placed at the rooftop level. The obstacles generate turbulence motions that destabilize the sharp shear layer atop the canyon, eventually increasing the turbulent mixing between the inner and outer flow. To the best of the authors’ knowledge, this is the first time that such a configuration has been investigated. The canyon fluid dynamics and the pollutant distribution are investigated numerically through highly resolved LES, provided with a dynamic Lagrangian turbulence model. Simulations reproduce a periodic and infinitely long urban canyon, with and without solid obstacles on the upwind building roof. The study of such an archetypal configuration aids in understanding the fundamental physical processes that rule the system and analyzing in detail the effects of the obstacles on the overall fluid dynamics. The results support the evaluation of the impact of roof roughness and roof irregularities, which appear in many practical applications.

The paper is organized as follow: in Sect. 2 the cases under examination are described; in Sect. 3 the simulation methodology and settings are presented; in Sect. 4 the numerical results are validated against laboratory experiment and numerical datasets; in Sect. 5 the effects of the roof obstacles are analyzed and the enticement in pollutant removal is estimated, and in Sect. 6 the final comments are given.

2 Problem Description

As an exploratory study, we numerically investigate the archetypal case of an infinite array of two-dimensional, on average, urban canyons with unity aspect ratio. The ambient airflow is perpendicular to the canyons, while a constant flux of passive scalar is released at the street level. Two geometries are analyzed: a smooth-roof case where the building roofs are flat, and an obstacle-roof case where obstacles are placed at the right corner of the building roofs.

2.1 Geometry

Fig. 1
figure 1

Computational mesh case geometry. The computational domain is bounded by dotted lines. The red line marks the area of pollutant release

Full size image

The case geometries are sketched in Fig. 1b together with the computational domain that is delimited by the dashed lines.

In the smooth-roof case, two square buildings create a canyon that has height \(H=1.0\) m and width \(W=1.0\) m (aspect ratio \(H/W=1\)). The computational domain is \(2H \times 2H \times 3H\) long in the streamwise (x), spanwise (y), and vertical (z) directions, respectively. The origin of the system of reference is the intersection between the horizontal line at the street level and the vertical line passing through the centre of the first building in the x-direction. The ambient flow is directed along the x-direction, perpendicular to the canyon, and is characterized by a freestream velocity \(U_{\mathrm{Ref}}\). This is defined as the average of the streamwise velocity component over the top face of the computational domain.

In the obstacle-roof case, 10 solid obstacles are placed at the right corner of the building roofs. They are cuboidal parallelepipeds equally distributed along the spanwise direction, with dimensions \(0.05H \times 0.05H \times 0.1H\) in x-, y-, z-directions respectively. The index of linear obstacle density, defined as the ratio of the distance between two consecutive obstacles and the obstacle width, is \(\lambda _{\mathrm{obs}} = 3.0\).

The pollutant, of concentration c, is released at the ground level from a band of width 0.7H centred in the street (see red line in Fig. 1b). The pollutant flux per unit area \(q=\lambda \frac{\partial c}{\partial n}\) (\(\lambda \) is the molecular diffusivity, n is the surface-normal direction) is constant. The pollutant source does not extend the entire canyon width to model the presence of a footpath. The building walls are modelled as smooth surfaces.

2.2 Non-dimensional Numbers and Parameters

The non-dimensional number characterizing the urban canyon dynamics is the building Reynolds number, defined as

$$\begin{aligned} { Re} = \frac{U_{\mathrm{Ref}} H}{\nu }, \end{aligned}$$
(1)

where \(\nu \) is the kinematic viscosity of the fluid. An additional parameter is the non-dimensional pollutant concentration,

$$\begin{aligned} C^* = \frac{c}{C_{\mathrm{Ref}}}, \end{aligned}$$
(2)

where \(C_{\mathrm{Ref}}= q/U_{\mathrm{Ref}}\) is the reference concentration.

In the following, lengths are made non-dimensional by means of buildings height H, and velocities are normalized by the freestream velocity \(U_{\mathrm{Ref}}\). Two characteristic times are defined: an outer-flow time, \(T = H/U_{\mathrm{Ref}}\), which represents the time scale of the ambient wind, and the internal-flow time, \(\tau = H/0.1U_{\mathrm{Ref}}\), which is the time scale of the main vortex appearing within the canyon. The latter is found by empirically estimating that the characteristic velocity within the canyon is one-tenth of the ambient characteristic velocity.

3 Simulation Methodology

Air is modelled as an incompressible fluid, and the pollutant concentration is considered a passive scalar, i.e., just transported and diffused by the fluid flow. The equations of motion read

$$\begin{aligned} \frac{\partial u_i}{\partial x_i}&= 0 , \end{aligned}$$
(3)
$$\begin{aligned} \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j}&= -\frac{1}{\rho _0} \frac{\partial p}{\partial x_i} + \nu \frac{\partial ^2 u_i }{\partial x_j \partial x_j}, \end{aligned}$$
(4)

where \(u_i\) are the velocity components, p is the pressure term, \(\rho _0\) is the reference density, and the summation convention over repeated indexes is adopted. The pollutant concentration c is governed by

$$\begin{aligned} \frac{\partial c}{\partial t} + u_i \frac{\partial c}{\partial x_i} = \lambda \frac{\partial ^2 c}{\partial x_i \partial x_i} ; \end{aligned}$$
(5)

being a passive scalar, it is just advected by the flow and diffused by the fluid.

3.1 Large-Eddy Simulation and Subgrid-Scale Model

In the LES approach, the computational grid acts on the governing equations as an implicit spatial filter (e.g., Sagaut 2000 and Rodi et al. 2013 among others). The filtering of the momentum (4) and concentration (5) equations generates additional terms that account for the unresolved subgrid scales of motion. They read, respectively,

$$\begin{aligned} \frac{\partial \tau _{ij}}{\partial x_i}&= \frac{\partial }{\partial x_i} \left( \overline{u_i u_j} - \overline{u}_i \overline{u}_j \right) , \end{aligned}$$
(6)

which are the SGS kinematic momentum fluxes, and

$$\begin{aligned} \frac{\partial h_{i}}{\partial x_i}&= \frac{\partial }{\partial x_i} \left( \overline{c u}_i-\overline{c} \ \overline{u}_i \right) , \end{aligned}$$
(7)

which is the SGS concentration flux. The overbar denotes the spatial filter of width equal to the computational cell size. The characteristic cell size is estimated as \(\overline{\Delta }= ( {\Delta }_x {\Delta }_y {\Delta }_z )^{1/3}\), where \(\Delta _i\) is the extension of the cell in the i-direction (notice that in this case the overbar does not represent the application of the filter to the computational grid).

The momentum equations are closed using the dynamic Lagrangian model proposed by Meneveau et al. (1996). First, the SGS kinematic momentum fluxes are expressed by means of the Smagorinsky (1963) model

$$\begin{aligned} \tau _{ij} - \frac{\delta _{ij}}{3}\tau _{kk}&= -2 \nu _{\textsc {sgs}} \overline{S_{ij}}, \end{aligned}$$
(8)
$$\begin{aligned} \nu _{\textsc {sgs}}&= c_s^2 \overline{\Delta }^2 |\overline{S_{ij}}|, \end{aligned}$$
(9)

where \(\nu _{\textsc {sgs}} \) is the SGS viscosity, \(c_s\) is the Smagorinsky constant, \( \overline{S_{ij}} = (1/2)[(\partial \overline{u}_i / \partial x_j) + (\partial \overline{u}_j / \partial x_i)]\) is the resolved strain rate tensor, and \(|\overline{S_{ij}}|=\sqrt{2\overline{S_{ij}} \ \overline{S_{ij}}}\) is its magnitude.

Subsequently, the constant \(c_s^2\) is determined by comparison between two scales of motion. The information on small scales is obtained by means of the grid filter of width \(\overline{\Delta }\), while for the larger scales an additional filter of width \(\widehat{\Delta } = 2 \overline{\Delta }\) is used. Following the dynamic Lagrangian procedure, the value of the Smagorinsky constant is computed, minimizing the modelling error along the streamlines (in order to stabilize the computation). In the end, the constant reads

$$\begin{aligned} c_s^2= \frac{ \mathcal {I}_{LM}({x}, t)}{\mathcal {I}_{MM}({x},t)}, \end{aligned}$$
(10)

where the additional quantities \(\mathcal {I}_{LM}\) and \(\mathcal {I}_{MM}\) are approximated by a sequence specified by recursion on index n (which in practice is the simulation timestep). The sequences specified by recursion are

$$\begin{aligned} \left\{ \begin{array}{l} \mathcal {I}_{LM}^{n+1}({x}) = \epsilon [L_{ik} M_{ik}]^{n+1}({x}) + (1- \epsilon ) \mathcal {I}_{LM}^{n} ({x}-\overline{{u}}^n \Delta t) \\ \mathcal {I}_{LM}^{0}({x}) = c^2_{s,0} [M_{ik} M_{ik}]^0({x}) \end{array} \right. \end{aligned}$$
(11)

and

$$\begin{aligned} \left\{ \begin{array}{l} \mathcal {I}_{MM}^{n+1}({x}) = \epsilon [M_{ik} M_{ik}]^{n+1}({x}) + (1- \epsilon ) \mathcal {I}_{MM}^{n} ({x}-\overline{{u}}^n \Delta t) \\ \mathcal {I}_{MM}^{0}({x}) = [M_{ik} M_{ik}]^0({x}) \end{array} \right. \end{aligned}$$
(12)

with the first iteration marked by the 0 index, \(c^2_{s,0}=0.0256\) the standard Smagorisky constant (squared), \(\Delta t\) the simulation timestep, and

$$\begin{aligned} \epsilon = \frac{\Delta t / \textsc {t}^n}{1 + \Delta t / \textsc {t}^n}, \quad \textsc {t}^n = 1.5 \; \overline{\Delta } (\mathcal {I}_{LM}^n \mathcal {I}_{MM}^n)^{-1/8}. \end{aligned}$$
(13)

Additionally, the tensors \(L_{ik}\) and \(M_{ik}\) are defined as

$$\begin{aligned} L_{ik} = \widehat{\overline{u}_i \overline{u}}_k - \widehat{\overline{u}}_i\widehat{\overline{u}}_k, \quad M_{ik} = 2 \overline{\Delta } ^2 \left( \widehat{|\overline{S}|\overline{S}_{ik}}-4|\widehat{\overline{S}}|\widehat{\overline{S}}_{ik} \right) . \end{aligned}$$
(14)

A numerical clipping is applied to \(\mathcal {I}_{LM}\), \(\mathcal {I}_{MM}\) in order to enforce positive values of the constant.

The concentration equation is closed adopting the Reynolds analogy for the SGS diffusivity. The SGS concentration flux reads

$$\begin{aligned} h_{i} = -\lambda _{\textsc {sgs}} \frac{\partial \overline{c}}{\partial x_i}, \end{aligned}$$
(15)

where \(\lambda _{\textsc {sgs}} = \nu _{\textsc {sgs}} / Sc_t \), with \(Sc_t\) the turbulent Schmidt number. Further details on the SGS model adopted here can be found in Cintolesi et al. (2015).

Notice that the dynamic Lagrangian model is an SGS model particularly suitable for reproducing anisotropic and localized turbulence (which characterizes the present case) because the turbulent content of the flow is determined cell by cell, taking information from the resolved scales of motion. Moreover, the model empirical content is reduced since the Smagorinsky constant is not set a priori.

3.2 Settings and Physical Parameters

The building Reynolds number is set to \({ Re} = 2.0 \times 10^{4}\), which is greater than the critical values, ensuring the Reynolds number independence for a square canyon flow (Michioka et al. 2011; Chew et al. 2018). Hence, the simulation of the smooth-roof case reproduces the same fluid dynamics as in real cases. In the obstacle-roof case the Reynolds number independence has not been yet assessed. The fluid is driven by an imposed pressure gradient, that is dynamically computed to ensure that the average velocity at the domain top face is equal to \(U_{\mathrm{Ref}}=0.2\) m s\(^{-1}\). The kinematic viscosity is \(\nu =1 \times 10^{-5}\) m\(^{2}\) s\(^{-1}\), and the turbulent Schmidt number (the ratio between turbulent viscosity and turbulent diffusivity of pollutant concentration) is \(Sc_t = 0.85\).

The boundary conditions are: at the fluid vertical faces, cyclic conditions for all the variables to reproduce the infinite extension in streamwise and spanwise directions; at the solid walls, no-slip conditions for velocity and zero-gradient for pressure and concentration; at the domain top face, zero-gradient for velocity and zero-value for pressure and concentration. The zero-gradient condition for the velocity at the top boundary is required to correctly reproduce the removal of the pollutant from the domain; in particular, this implies that \(u_z\) can be non-zero at the top boundary.

The simulations were initialized with a uniform velocity field of \(u_x=0.1\) m s\(^{-1}\) and zero values for other variables. To speed-up the convergence to a statistical steady state, a preliminary simulation using the Smagorinsky SGS model was performed. Once the field was developed, the simulation was re-run using the dynamic Lagrangian SGS model. After the steady state established, the simulation was run for a period of \(6 \tau \) where the statistics were collected.

3.3 Computational Mesh

An overview of the computational grid for the obstacle-roof case is given in Fig. 1a. The grids of both cases were generated using the utility snappyHexMesh of Open- FOAM (version 6.0) , by setting up a background mesh refined in selected areas with an iterative procedure. The background mesh was composed of \(144 \times 144 \times 80\) cells in the x,  y,  z direction, respectively; they were slightly stretched to have a denser distribution at the roof level (\(z/H=1\)). The refinement procedure was the following: first, the cells within the solid bodies were removed (buildings and obstacles, if present); second, the mesh near the solid walls was refined (cell dimensions divide by factor two); third, supplementary layers of stretched cells were added near the wall surfaces. Preliminary simulations were run using coarse refinement, which was subsequently increased to ensure direct resolution of the boundary layer in the entire domain and avoid the use of wall functions. To this end, the final meshes ensure that the first computational point near the walls is placed at \(y^+ = u_\tau y / \nu < 1\), where y here is the distance normal to the wall and \(u_\tau \) is the friction velocity. To check the adequacy of the background mesh used, an additional simulation for the smooth-roof case is performed increasing the resolution of the background mesh of \(20\%\) in all directions. Mean and root-mean-square (r.m.s.) profiles of velocity components (along selected vertical lines) obtained from the original and refined mashes collapse on the same profiles, showing no detectable differences. In the obstacle-roof case, the mesh was also refined in the area behind the obstacles to better capture the wake flow (see Fig. 1a). The final mesh for the smooth-roof case counts 1,592,562 cells, while for the obstacle-roof case the total number of cells is 4,087,700.

3.4 Numerical Set-Up

The simulation is carried out using the open-source software Open FOAM (version 6.0). A home-made numerical solver was created to implement the mathematical model described above. The pimpleFoam basic code was customized to also solve pollutant concentration equations, while the dynamic Lagrangian model was implemented ex novo. Such a turbulent model has been previously tested and validated (Cintolesi et al. 2019).

The resolution algorithm is based on the Pressure-Implicit with Splitting of Operators (PISO) proposed by Oliveira and Issa (2001). Additional details are reported in the official Open FOAM (version 6.0) documentation. The discretization in time uses an implicit Euler backward schemes while in space a second-order central difference scheme is used for all terms except the concentration advective term. For such a term, we use a Gauss-Gamma scheme proposed by Jasak et al. (1999) with \(\gamma = 0.2\). This is a bounded version of central-difference scheme for unstructured meshes based on the normalized variable diagram. Overall, the resolution accuracy is of second-order. A numerical clipping on concentration is imposed to prevent non-physical (negative) values that can appear in a negligible number of cells at the street level during the development of the flow. After the statistical steady-state regime is reached, the negative fluctuations practically disappear and the plots obtained without clipping do not show detectable discrepancies with respect to the simulation with clipping.

4 Validation

The smooth-roof case simulation is validated against existing laboratory experiments. The first- and second-order velocity statistics are compared with three datasets present in the literature: Brown et al. (2000), Li et al. (2007), and Chew et al. (2018). In Chew et al. (2018), a water-channel experiment is performed using seven canyons with \(W=0.12\) m, taking the measurements at canyons 3–6. The inflow was stabilized by a combination of honeycombs, plastic tubes, and wire mesh, while a layer of ceramic marbles accelerated the development of the incoming flow profile. Velocity measurements were made with acoustic Doppler velocimeters, with a sampling volume of 0.06 m\(^3\). Li et al. (2007) used a laboratory water flume containing seven canyons of width \(W=0.1\) m. Two-colour laser Doppler anemometry was used for velocity measurements at the central canyon, at a vertical plate centred in the spanwise direction. Brown et al. (2000) carried out a wind-tunnel experiment including six canyons at \(W=0.15\) m. Irwin spires and roughness elements at the inlet simulate a neutral atmospheric boundary layer. A pulsed-wire anemometer was used for velocity measurements; the flow reaches the equilibrium at the 3rd–4th canyons, where the measurements were taken. The concentration profile is validated against the LES of Michioka et al. (2011), and the measurements of Meroney et al. (1996) and Pavageau and Schatzmann (1999).

Statistics are computed in time (over a period of \(6 \tau \)) and in space along the spanwise direction, taking advantage of the spatial homogeneity. If \(\psi \) is a generic variable, \(\langle \psi \rangle \) indicates the time–spatial average, \(\psi ' = \psi - \langle \psi \rangle \) is the fluctuation, and \(\psi _{\mathrm{rms}} = \sqrt{\langle \psi '^2 \rangle }\) is the root-mean-square (r.m.s.) value.

4.1 First-Order Statistics

The velocities are plotted along three vertical lines within the canyon near the upwind building (\(x/H = 0.75\)), at the centreline (\(x/H=1.00\)), and near the downwind building (\(x/H=1.25\)).

Fig. 2
figure 2

Non-dimensional mean velocity components along three vertical lines: \(x/H = 0.75, 1.00, 1.25\). Black solid line, present LES. Symbols are measurements: blue \(+\), Chew et al. (2018); black \(\circ \), Li et al. (2007); red \(\square \), Brown et al. (2000)

Full size image

Figure 2a shows the non-dimensional mean streamwise velocity component at the selected vertical lines. The present LES results practically collapse to the reference data. In the canyon internal region (\(z/H<1.0\)), a weak primary vortex develops (see Fig. 8) with maximum velocities of one order of magnitude lower than the freestream velocity \(U_{\mathrm{Ref}}\). Two inflection points delimiting the shear layer atop the canyon are detectable; as expected, the distance between the points increases near the downwind building (Letzel et al. 2008). In the external region (\(z/H>1.0\)), a logarithmic boundary layer is established.

Figure 2b displays the non-dimensional mean vertical velocity. The experimental profiles exhibit different results at the centreline in the canyon internal region: Brown et al. (2000) measures almost zero velocities, while Li et al. (2007) and Chew et al. (2018) report positive values and the maximum at \(z/H \approx 0.5\). The former indicates that the centre of the primary vortex is close to the canyon centre (\(x/H=z/H=0.5\)); the latter reveals that the centre is shifted downwind and upward. The present LES slightly underestimates the velocity magnitude near the upwind building, being in agreement with the profile of Brown et al. (2000) at the canyon centreline, and is in good agreement with measurements near the downwind building.

Table 1 Summary of settings used in laboratory experiments of urban canyon with unity aspect ratio
Full size table

There is no clear consensus in the literature about the genesis of the vortex displacement: in general, numerical simulations do not display such a feature, whereas laboratory investigations exhibit conflicting results summarized in Table 1. Li et al. (2008) suggested that the displacement is due to the interference of the experimental apparatus boundaries: an insufficient spanwise extension can produce motions parallel to the canyon axis altering the two-dimensional flow, while a limited height can generate a pressure-suction effect that pushes the vortex centre upward. This hypothesis seems to be supported by previous studies reported in Table 1; the vortex displacement appears in small domains (height \(\le 5 H\) and width \(< 9 H\)) while for larger experimental domains it does not appear. To scrutinize the existence of such boundary effects, two additional test simulations are performed setting as wall the top boundary (face \(z/H=3\)) and the spanwise boundaries (faces \(y/H=0,2\)). The non-slip condition for u and the zero-gradient condition for p is imposed at the added walls. The Spalding wall function is used since the mesh is not fine enough to directly resolve the boundary layer on these boundaries (Spalding 1961). The other settings do not change. Neither simulation exhibits relevant displacement of the primary vortex centre. The vertical profiles show a minimal variation with respect to the ones reported in Fig. 2 (the comparison is not reported). Therefore, the top and lateral solid boundaries do not produce a relevant effect on the mean flow profile in the current geometry. Additional simulations are required to detect the cause of vortex displacement; however, the authors suggest that the displacement may be caused by the inflow condition or some alterations of the external flow (see the unexpected downward vertical velocity at the centreline for \(z/H>1\), Fig. 2b).

Fig. 3
figure 3

Non-dimensional mean concentration \(C^*\) at the upwind building facade (\(x/H=0.5\)) centreline, and downwind building facade (\(x/H=1.5\)). Simulation: black solid line, present LES; red dashed line, LES by Michioka et al. (2011). Measurements: blue \(\diamond \), Pavageau and Schatzmann (1999); green \(*\), Meroney et al. (1996)

Full size image

Figure 3 presents the non-dimensional mean concentration profile \(C^*\) defined by Eq. 2, at the centreline of the upwind building facade (\(x/H=0.5\)) and downwind building facade (\(x/H=1.5\)). In the external region, the reference profiles rapidly decrease to zero because the pollutant is realized in a single canyon, and clear fluid flows atop the canyon. The LES profiles assume non-zero values that decrease almost linearly with height as a result of the removal from the upwind canyons. In the internal region, present LES slightly underestimates the concentration at the upwind building facade, while it is in a satisfactory agreement with the references at the centreline and at the downwind facade. The discrepancy at the upwind facade is possibly attributed to the different pollutant source (linear for the references, planar for the present case), which can influence the removal by the primary canyon vortex near the upwind building (Pavageau and Schatzmann 1999). The planar pollutant source is chosen in this study because it is estimated to be a more realistic model of traffic pollution.

4.2 Second-Order Statistics

Fig. 4
figure 4

Non-dimensional mean TKE (a), resolved kinematic momentum flux \(\langle u'w' \rangle \) and velocity r.m.s. (b) along three vertical lines: \(x/H = 0.75, 1.00, 1.25\). Black solid line, present LES. Symbols are measurements: blue \(+\), Chew et al. (2018); black \(\circ \), Li et al. (2007); red \(\square \), Brown et al. (2000)

Full size image

Figure 4a shows the resolved turbulence kinetic energy (\(\mathrm{{TKE}}, e= 0.5 \langle u'_i u'_i \rangle \)) at the selected vertical lines along with the resolved kinematic momentum flux \(\langle u'w' \rangle \) at the canyon centre line. It is worth noting that the experiments of Chew et al. (2018) and Brown et al. (2000) do not reproduce the well-documented TKE peak at the roof level (\(z/H=1\)) generated by the turbulent shear layer from the upwind roof (Louka et al. 2000; Kastner-Klein et al. 2004; Letzel et al. 2008), whereas they exhibit high values in the canyon external region. These high values may be due to transport by the background flow rather than production by the interaction with buildings, which is limited in the turbulent shear layer near atop the canyon. To the best of the authors’ knowledge, there is not a clear explanation in the literature of the dynamics underlying the formation of this high TKE in the external flow. The TKE profiles are in excellent agreement with all the references in the internal region and reproduce the laboratory experiment of Li et al. (2007) in the external region. Surprisingly, the kinematic momentum flux of the present LES at the centreline follow the profile reported by Chew et al. (2018), which presents the expected roof-level minimum. This can suggest that the above-discussed higher value of TKE in the external region is possibly due to large spanwise fluctuations, unexpected in two-dimensional flow. The negative values of kinematic momentum flux reveal the turbulent motion feeds the mean reverse flow from the external to the internal region.

Figure 4b reports the spanwise and vertical velocity r.m.s. values, compared with the dataset of Li et al. (2007): there is a very good agreement of the simulation profiles with the reference at all the locations.

Overall, the simulation satisfactorily reproduces the reference profiles, whilst some discrepancy among the laboratory data are detected and discussed.

5 Results

The obstacle-roof case is now analyzed with respect to the smooth-roof case. Unless otherwise specified, the statistics for the former are made along the spanwise direction for a direct comparison with the latter.

5.1 The Obstacle-Roof Case Dynamics

The presence of the obstacles channels the flow atop the roof in the gap between two consecutive obstacles. To discuss this dynamic, the time-average velocity in two \(x-y\) planes is shown: the behind-obstacle plane, that cuts an obstacle at the mid-width, and the behind-gap plane, that passes through the gap centre.

Fig. 5
figure 5

Streamlines of the time-averaged velocity at the behind-obstacle plane, cutting the obstacle at the mid-width, and distribution of the non-dimensional time-averaged vertical velocity

Full size image

Figure 5 displays the time-averaged streamlines and the non-dimensional mean vertical velocity distribution at the behind-obstacle plane. It is worth highlighting the dynamics around the obstacle: at the upwind face, the flow deviates above and to the side of the obstacle, and a weak stagnation point appears at the bottom corner. At the downwind face, a positive vertical velocity drives the fluid upward from the canyon. Such velocity injects in the external region part of canyon air that escapes from the canyon. Within the canyon, three corners host recirculation regions that are also present in the smooth-roof case and are later discussed (see Fig. 8). The streamlines at the behind-gap plane are very similar to those of the smooth-roof case, hence, they are not shown. The gap centre is sufficiently far from the obstacles to not be affected by strong perturbation (see also Fig. 7).

Fig. 6
figure 6

Non-dimensional time-averaged velocity components at three locations: \(x/H=0.6,1.00,1.4\). Solid line: smooth-roof case; dashed lines: profiles at behind-obstacle plane; dotted lines: profiles at behind-gap plane

Full size image

Figure 6 presents the profiles of non-dimensional time-averaged velocity components at three locations: \(x/H=0.6,1.00,1.4\). The profiles at the behind-obstacle and behind-gap plane are compared; the profiles of the smooth-roof case are also reported. At \(x/H=0.6\), the streamwise velocity component behind the obstacle is largely decreased, while behind the gap there is almost no difference with the smooth-roof velocity. The vertical velocity behind-obstacle shows a peak in the vicinity of the obstacle top and slightly higher values in the upper part of the canyon (\(z/H>0.75\)) compare to the smooth-roof case. Additionally, the profile behind the gap exhibits slightly higher values. At \(x/H=1.0\), the influence of the obstacle is attenuated: all the curves practically collapse on the same profile, even if the streamwise velocity component of the obstacle-roof case is slightly weaker compared to the smooth-roof case in the external region. At \(x/H=1.4\), the streamwise velocity component does not present significant variation. The vertical velocities of the obstacle-roof case show small negative values at the roof level and the almost zero values at the street-level. The former is possibly due to the increase of the primary vortex height, while the latter is caused by the different horizontal extension of the recirculation vortex at the bottom-downwind corner (\(x/H=1.5, \ z/H=0\)). These features are reported in Sect. 5.2 (see Fig. 8).

Fig. 7
figure 7

Resolved TKE made non-dimensional by \({ {e}}_{\mathrm{Ref}}=0.5 U_{\mathrm{Ref}}^2\). Distribution at the horizontal plane \(z/H=1.05\)

Full size image

Figure 7 shows the TKE distribution over a horizontal plane cutting the obstacle at the mid-height \(z/H=1.05\). The higher values are triggered at the lateral faces of the obstacles and propagate streamwise in the wake flow. Behind the obstacles, there is a zone of very low TKE with a horizontal extension of about \(0.5<x/H<0.6\). This is a low-turbulence region where the upward flow drives the fluid from the internal region to the external one (see Fig. 5).

5.2 Main Flow and Turbulence Features

In this section, the mean flow behaviour is discussed irrespectively of the spanwise variation introduced by the obstacles. Thus, all the quantities are averaged in the y-direction.

Fig. 8
figure 8

Mean velocity streamlines in the internal region of the canyon: the smooth-roof case (left) and the obstacle-roof case (right). The average is in time and in the y-direction. Same system of reference as in Fig. 1b: horizontal axis is the x-direction, vertical axis is the z-direction

Full size image

Figure 8 reports the streamlines in the internal region. In both the configurations, the canyon is dominated by the clockwise primary vortex whose centre \(\gamma \) is close to the geometrical centre of the canyon: \(\gamma =(1.030,0.546)\) in smooth-roof case, \(\gamma =(1.036,0.528)\) in obstacle-roof case. Three stable, low-velocity, recirculation zones are detectable at the canyon corners, except in the top-downwind corner where the principal external flow impinges and creates a stagnation point. In the obstacle-roof configuration the primary vortex reaches the height of 0.04H, thanks to the reduction of the horizontal velocity by the obstacle. The height of the vortex is defined as the maximum height of the interface between the vortex and the ambient flow. The counterclockwise recirculation zone at the bottom-upwind corner (\(x/H=1.5, \ z/H=0\)) is slightly reduced, while the size of the zone at the top-upwind corner (\(x/H=1.5, \ z/H=1\)) decreases by half in the horizontal. The recirculation zone at the top-upwind corner is created by the sharp and highly energetic (high values of TKE) shear layer leaving the upwind roof that acts as a barrier for the up-willing flow. The wake flow triggered by the obstacles destabilizes such layers, reducing the barrier effect. The recirculation at the bottom-downwind corner slightly increases the horizontal extension (see also discussion of Fig. 6). The stagnation point at the top-downwind corner is lightly shifted upward of 0.01H.

Fig. 9
figure 9

Non-dimensional mean \(\mathrm{TKE}\) and resolved kinematic momentum flux \(\langle u'w' \rangle \) along three vertical lines: \(x/H = 0.75, 1.00, 1.25\). Solid lines, smooth-roof case; dashed lines, obstacle-roof case

Full size image

Figure 9 presents the non-dimensional vertical profiles of resolved TKE and the kinematic momentum flux for the two cases simulated. Overall, the profiles in the internal region practically overlap each other, both for the TKE and the kinematic momentum flux.

The TKE in the obstacle-roof case assumes higher values in the external region and does not exhibit the peak at the roof level, except near the upwind building where there is a localized climax due to the obstacle wake flow. The maximum value is slightly shifted upward with respect to the smooth-roof case. The obstacles destabilized the sharp shear layer that characterized the smooth-roof case, providing a comparable level of TKE, which is transported in the external region. It is worth noting that the TKE centreline profile is qualitatively similar to those reported by Chew et al. (2018) and Brown et al. (2000) in Fig. 4. This supports our hypothesis that experimental results have been altered by an additional turbulent content in the external flow, possibly triggered by the devices (e.g., Irwin spires and roughness elements) use to generate a turbulent inflow.

The kinematic momentum flux displays similar values for the two cases also in the external region. Differences are detectable at the obstacle-top level \(0.95<z/H<1.15\), where the obstacle-roof case has lower values at the centreline and near the downwind building. Near the upwind building, the shear triggered at the obstacle top \(z/H=1.1\) and at the gap \(z/H=1.0\) generates two local maxima. Similar to the TKE, in the obstacle-roof case the roof-level peak is inhibited.

Instantaneous turbulent structures are also described. For this purpose, it is customary to analyze the spanwise vorticity or the pressure fluctuations. Following the discussion by Dubief and Delcayre (2000), the Q-criterion investigation is here preferred. This is based on the study of the second invariant of the velocity gradient

$$\begin{aligned} Q = \frac{1}{2}\left( \Omega _{ij} \Omega _{ij} - S_{ij} S_{ij} \right) , \end{aligned}$$
(16)

where \(\Omega _{ij}= 0.5(\partial u_i/\partial x_j - \partial u_j/\partial x_i)\) is the rotation strain rate, Q is interpreted as the balance between the rotation and strain rate, and positive isosurfaces are entitled to be delimiters of turbulent eddies.

Fig. 10
figure 10

Turbulent structures visualized by the isosurface of scalar Q, together with the instantaneous velocity magnitude at a vertical slice in background: the smooth-roof case (left) and the obstacle-roof case (right). The isosurface corresponds to a suitable threshold and is three-dimensional, i.e., all structures in spanwise direction are shown

Full size image

Figure 10 visualizes the turbulent structures of the two cases by means of the Q-criterion. In the smooth-roof case, the structures are sparse, and limited in a horizontal zone of about \(0.9<z/H<1.5\) and near the downwind building facade. In the obstacle-roof case, the structures are more dense, and numerous smaller structures are generated in the proximity of the obstacles. They penetrate more in the external region, spreading in a large horizontal zone approximately delimited by \(0.9<z/H<2.0\).

5.3 Interface Turbulent Exchange

The interface between the internal canyon region and the external background flow (the horizontal plane \(z/H=1.0\)) is a critical area for pollution removal.

Fig. 11
figure 11

Profiles of non-dimensional mean vertical velocity, TKE, the \(\langle u w \rangle \) budget along the canyon–atmosphere interface \(z/H=1.0\) for the smooth-roof case (solid lines) and the obstacle-roof case (dashed lines)

Full size image

Figure 11 shows the non-dimensional mean profiles of vertical velocity, TKE, and the velocity flux budget

$$\begin{aligned} \langle u w \rangle = \langle u \rangle \langle w \rangle + \langle u'w' \rangle , \end{aligned}$$
(17)

which is the sum of the mean flux and the kinematic momentum flux. The vertical velocity of the smooth-roof case is almost zero in the upwind part (\(x/H<1.0\)) because of the large streamwise velocity that separates the internal and the external region. In the downwind part (\(x/H>1.0\)) is characterized by the strong positive peak due to the stagnation point at the downwind building facade. The obstacle-roof case exhibits positive values in the upwind part and negative values in the downwind part given by the rise of the primary vortex, and no peak because of the upward shifting of the stagnation point (see also Fig. 8). The TKE in the obstacle-roof case has lower values as a result of the vertical dispersion of the TKE, except near the upwind building where the obstacles are more effective in triggering turbulent motion. The mean vertical fluxes are greater in the obstacle-roof case, consistent with the vertical velocity profile. The kinematic momentum fluxes are similar in the two cases, but near the upwind building (\(x/H<0.75\)) the obstacle-roof case displays larger negative values associated with more intense events of ejection and sweeps. Overall, a more intense vertical exchange by the mean and turbulent flow is detectable in the obstacle-roof case.

Fig. 12
figure 12

Contour plot of instantaneous non-dimensional fluctuations of resolved velocity components at the horizontal plane \(z/H=1.05\). Planar velocity vectors are also reported

Full size image

Figure 12 shows the contourplot of velocity component fluctuations, along with the planar velocity vectors, at a horizontal plane \(z/H=1.05\). The areas of negative values (associated with ejection and sweep events) are more numerous and less extended in the obstacle-roof case, and they are more localized in the region \(x/H<0.75\) behind the obstacle. Hence, this can be seen as a region of higher turbulent exchange.

Fig. 13
figure 13

Joint probability density function of fluctuations of the resolved the velocity components at three points at the roof level \(z/H=1.0\), at spanwise centreline \(y/H=1.0\) (i.e., between obstacles in the obstacle-roof case), and at the streamwise locations \(x/H=0.6, 1.0, 1.4\). Time series data from virtual probes

Full size image

Figure 13 displays the joint probability density function of the streamwise and vertical velocity fluctuations (\(u'\) and \(w'\)) for the obstacle-roof case and the smooth-roof case. They are computed at the roof level \(z/H=1.0\) and at three horizontal locations \(x/H=0.6,1.0,1.4\), where virtual probes registered a time series of velocity data. Following the quadrant analysis, the first quadrant (\(u'>0, w'>0\)) records outward interaction, the second (\(u'<0, w'>0\)) ejection events, the third (\(u'<0, w'<0\)) inward interaction, the fourth (\(u'>0, w'<0\)) sweep events.

In general, ejections and sweeps are dominant with respect to outward/inward interaction in both cases. The former transports low-momentum fluid upward and the latter transports high-momentum fluid downward; hence, they significantly influence the momentum fluxes and the pollutant removal (Michioka et al. 2011).

At \(x/H=0.6\), the obstacle-roof case shows a larger probability of ejection events with respect to the smooth-roof case. This is mainly due to the turbulent wakes generated by the obstacles, which decrease streamwise velocity and boost vertical turbulent motions (see also Fig. 12 that shows the turbulent structures localized behind the obstacles). At \(x/H=1.0\), the extreme values of fluctuations are reduced when obstacles are present. This is possibly due to the upward shifting of the interface between the primary vortex and the ambient flow (see discussion in Sect. 5.2): in the obstacle-roof case, the point (\(x/H=1.0, z/H=1.0\)) is located within the vortex where turbulence is reduced, whereas in the smooth-roof case, it is located at the turbulent interface. At \(x/H=1.4\), the differences between the two cases are less pronounced, but the obstacle-roof case exhibits a slightly lower probability of turbulent events with respect to the smooth-roof case. Overall, the obstacle-roof case exhibits an increase of turbulent fluctuations at the roof level behind the obstacles (\(x/H=0.6\)), while it shows a decreasing at the canyon centreline (\(x/H=1.0\)) and near the upwind building (\(x/H=1.4\)). This is in agreement with the TKE profiles in Fig. 11 and the turbulence structure distribution displays in Fig. 12. Such behaviours can be ascribed to the rise of the turbulent interface between the primary canyon vortex and the ambient flow, which shifts the TKE peak upward as seen in Fig. 9.

5.4 Pollutant Removal from the Canyon

Table 2 The mean recirculation quotient \(\Phi \) defined in Eq. 18, the convective transfer coefficient \(\Omega \) of pollutant defined in Eq. 19, and its positive and negative parts
Full size table

The pollutant concentration is used to quantify the overall fluid recirculation in the canyon. Following Göthe et al. (1988), the mean recirculation quotient is defined as

$$\begin{aligned} \Phi = { \int _{\mathrm{canyon}} \langle C^* \rangle \mathrm{d}V }\Big /{ \int _{\mathrm{domain}} \langle C^* \rangle \mathrm{d}V }, \end{aligned}$$
(18)

where \(\langle C^* \rangle \) is the non-dimensional concentration averaged in time and space, integrated on the canyon volume (numerator) and on the entire domain (denominator); see also Fig. 1. The quotients for the two cases simulated are reported in Table 2: in the roof-obstacle case, the recirculation quotient \(\Phi \) is sensibly smaller than in the smooth-roof case. This indicates that obstacles induce a stronger dispersion of the pollutant from the canyon to the atmosphere. The estimation of the pollutant concentration within the canyon shows that the obstacle-roof case has a concentration \(34\%\) lower that the smooth-roof case.

Fig. 14
figure 14

Mean non-dimensional concentration profile (left) and concentration r.m.s. (right) along three selected vertical lines, as in Fig. 2. Solid line, smooth-roof case; dashed lines, obstacle-roof case

Full size image

Figure 14 shows the vertical profiles of mean concentration \(C^*\) and its r.m.s value. With respect to the smooth-roof case, the obstacle-roof case exhibits a lower pollutant concentration within the canyon and a more homogeneous distribution in the vertical direction as a consequence of the higher turbulent mixing. The r.m.s. profiles are very similar and display larger values at the street level and a local peak at the roof level.

To scrutinized the pollution transfer through the roof interface, the convective transfer coefficient of the pollutant (Barlow and Belcher 2002) is computed as

$$\begin{aligned} \Omega = \frac{1}{H C_{\mathrm{Ref}} U_{\mathrm{Ref}}} \int _{\mathrm{roof}} \langle c w \rangle \mathrm{d}x. \end{aligned}$$
(19)

The removal \(\Omega ^+\) and re-entrainment \(\Omega ^-\) components are also estimated using the Eq. 19 where the positive and negative part function, respectively, is applied to w. Table 2 displays the values of the coefficients. The convective transfer coefficient (per unit of time) is \(7\%\) higher in the obstacle-roof case, where the removal component is slightly higher while the re-entrainment component is slightly lower. The rate of re-entrainment can be measured by the ratio \(\Omega ^-/\Omega ^+\) (Cheng and Liu 2011), which shows that in the obstacle-roof case the re-entrainment is reduced about \(1.6\%\). This is in accordance with the elevation of the primary vortex boundaries and the separation points at the downwind building facade, which are detected in the obstacle-roof case (see discussion of Fig. 8).

Fig. 15
figure 15

Contourplot of vertical turbulent concentration flux \(\langle c'w' \rangle \) made non-dimensional

Full size image

Figure 15 shows the distribution of vertical turbulent mass fluxes \(\langle c'w' \rangle \) (made non-dimensional). In both cases, two regions of high values are detectable at the street level (where the pollutant is released) and at the canyon–atmosphere interface. As expected, the smooth-roof case exhibits lower values with respect to the obstacle-roof case, and the region of larger fluxes is located approximately atop the roof in a central zone \(0.70<x/H<1.30\). In the obstacle-roof case, the zone of high fluxes at roof level is more extended both in the vertical and horizontal directions: it encompasses the entire canyon width, denoting a more efficient vertical mixing at the interface.

6 Conclusions

A numerical study of idealized urban canyons with and without obstacles on the building roof is performed to investigate the effects on pollutant removal of natural turbulence increases at roof level. This simplified configuration reproduces the fundamental physical processes governing real urban canyons, while it allows for detailed analysis of the dominant removal mechanisms.

Large-eddy simulation is used, along with the dynamic Lagrangian subgrid model, which is particularly accurate to reproduce anisotropic and localized turbulence. The paradigmatic case of periodic arrays of infinitely long square canyons is investigated. Two simulations are performed considering a smooth-roof (the building roof is a flat surface) and an obstacle-roof, where a series of solid obstacles are placed on the roof upwind side. Pollution is modelled by a passive scalar released at the canyon bottom surface.

The simulations are successfully validated against several laboratory and numerical datasets. Some discrepancies between the reference data concerning the displacement of the primary vortex are highlighted and discussed. The main effect of the obstacles is to destroy the sharp shear layer arising at the roof level, thus increasing the mixing between the canyon interior and the surrounding flow. Such an effect is particularly evident near the upwind building, where the ascending flow in the canyon transported the largest part of the pollutant upward. The principal flow features generated by the roof obstacles are:

  • The turbulent interface between the canyon primary vortex and the ambient flow is shifted upward, and its maximum height increases of \(4\%\). Also, the internal recirculation regions become smaller.

  • In the region near the upwind building at the roof level, the turbulent mixing increases: a higher level of TKE is generated, the probability of ejection events increases, and the vertical turbulent fluxes intensify. It is worth noticing that this is a crucial area where the pollutant accumulates transported by the primary vortex.

  • The ambient flow exhibits higher turbulent content (i.e., higher TKE) above the canyon (\(z/H > rsim 1.3\) at centreline) with peaks in correspondence to the vortex-ambient turbulent interface.

In terms of pollution removal, there is a reduction of the mean pollutant concentration within the canyon of \(34\%\) in the roof-obstacle case with respect to the smooth-roof case. Also, a slight increase in the convective transfer coefficient and a moderate reduction of re-entrainment rate are detected.

Finally, the presence of the obstacles at the building roof considerably increases the removal of pollutants from urban canyons. In practical situations, modern buildings are more and more often flat and free of obstacles; deliberately placing obstacles to increase turbulence at the roof level can improve ventilation and air quality in new urban areas. Although in the real case many other wind directions are possible, the authors focus on the perpendicular one because it is more studied in the literature and is more suitable to reveal the pollutant removal mechanisms. Moreover, real cities can exhibit a complex and irregular roof configuration which can influence the pollutant dispersion in different ways. In this work, a simple configuration is investigated with the aim of demonstrating that rooftop infrastructures can greatly influence pollutant removal. The evidence discussed in this study may be expanded in future works and form the basis of further applications in a realistic environment. To this end, subsequent studies should analyze the effectiveness of the presence of obstacles in more complex settings, for example, by varying the wind direction or the arrangement of obstacles.