Large-Eddy Simulation Analyses of Heated Urban Canyon Facades

: Thermal convective ﬂows are common phenomena in real urban canyons and strongly affect the mechanisms of pollutant removal from the canyon. The present contribution aims at inves-tigating the complex interaction between inertial and thermal forces within the canyon, including the impacts on turbulent features and pollutant removal mechanisms. Large-eddy simulations reproduce inﬁnitely long square canyons having isothermal and differently heated facades. A scalar source on the street mimics the pollutant released by trafﬁc. The presence of heated facades triggers convective ﬂows which generate an interaction region around the canyon-ambient interface, characterised by highly energetic turbulent ﬂuxes and an increase of momentum and mass exchange. The presence of this region of high mixing facilitates the pollutant removal across the interface and decreases the urban canopy drag. The heating-up of upwind facade determines favourable convection that strengthens the primary internal vortex and decreases the pollutant concentration of the whole canyon by 49% compare to the isothermal case. The heating-up of the downwind facade produces adverse convection counteracting the wind-induced motion. Consequently, the primary vortex is less energetic and conﬁned in the upper-canyon area, while a region of almost zero velocity and high pollution concentration (40% more than the isothermal case) appears at the pedestrian level. Finally, numerical analyses allow a deﬁnition of a local Richardson number based on in-canyon quantities only and a new formulation is proposed to characterise the thermo-dynamics regimes.


Introduction
The ongoing worldwide urbanisation and increasing population of the urban environments have raised awareness towards a series of environmental issues having adverse effects on public health [1,2]. Air pollution due to vehicular emissions is recognised as one of the major health threats [3], as a consequence of a large number of sources in a poorly ventilated environment and high level of exposure for the urban residents [4][5][6]. Understanding the flow circulation and ventilation within the urban canopy is paramount to predicting population exposure and developing solutions to improve residents' livability and comfort, thus requiring a morphology revision of the canopy environment [7] or the introduction of novel techniques to reduce pollutant concentrations [8].
Most of the recent investigations on the flow circulation and ventilation processes in urban canopies have delved into the characterisation of the inertial flow regime under isothermal conditions (i.e., driven by the mechanical force exerted by the free-stream wind velocity above the canopy) and less so on the thermal effects. Thermal stratification and convection are caused by the buoyancy force induced by differential surface heating by solar irradiation and contribute to the modification of the canopy flow and turbulent regimes, especially under weak inertial forcing [9]. The horizontal temperature gradient establishing between differently heated building facades forces a thermally-driven circulation within the canopy that can modify or superimpose the already-established inertially-driven vortex [9][10][11][12][13]. Contextually, turbulent energy and momentum transport is also modified [14][15][16][17], playing a significant role in the canyon pollutant dispersion and scalar mixing [17,18]. Additional buoyancy forces arise from the street and building surfaces [19] and rooftop heating [20], enhancing the convective forcing according to the solar path [14,21], and thus modifying the heat transfer [22] and thermally-related ventilation [23]. Boundary-layer unstable stratification further enhances the in-canyon ventilation due to the interaction of the canopy vortex with the large convective eddies above [24].
Different non-dimensional coefficients have been introduced to discern inertial from thermal flow regimes within an urban canyon. Among others, Dallman et al. [9] derived a buoyancy parameter as the ratio between buoyancy force generated by opposite building facades at different surface temperature and the free-stream velocity to determine if the in-canyon circulation is inertially or thermally driven. Magnusson et al. [11] refined this parameter by further characterising the mixed-phase condition, where both mechanical and thermal forces impact the observed circulation. Most commonly, the canopy-flow regimes are investigated in terms of a bulk Richardson number where g is the gravity acceleration, β is the thermal expansion coefficient, ∆T is the characteristic temperature difference, H is the mean building height, and U 0 free-stream velocity. Such a quantity compares the buoyancy force per unit mass gβ∆T between canyon surfaces (e.g., street, building facades, rooftops) with the wind-driven mechanical force per unit mass U 2 re f /H. The results of the buoyancy-mechanical force balance is oftentimes addressed in terms of convection regimes. For Ri 1, the in-canyon flow regime is inertially-driven and convection is suppressed (forced convection regime), while for Ri ∼ O(1) both thermal and inertial forces become significant (mixed convection regime). The third regime of natural convection is obtained at Ri 1 when buoyancy force dominates inertia (condition typically achieved under wind calm). Unlike the buoyancy parameter, the Richardson number does not have precise thresholds, but the convective regimes seem to depend on local atmospheric and morphology conditions. For instance, Li et al. [25] observed a direct transition between forced and natural convection at |Ri|= 20 in asymmetric urban canyons while Allegrini et al. [26] observed a mixed-phase at Ri ≥ 1.
In the context of non-isothermal urban canyon flows, the use of non-dimensional parameters provides a useful tool to discern the mechanical and thermal force balance. However, the characteristic temperature difference and free-stream velocity required by the Richardson number are not unequivocally determined in the literature, somehow limiting the comparability between different studies. As revealed by Zhao et al. [27], the characteristic temperature difference has been computed between the air at the rooftop and street level, windward and leeward building facades or a heated surface and ambient air, thus limiting the inter-study comparability. Moreover, in field experiments the canyon and above atmosphere are rarely isothermal; thus, the ambient temperature is a critical parameter too. This aspect is typically overtaken in the laboratory and numerical experiments where an isothermal ambient flow is perturbed by heating a target surface. The free-stream velocity is also not unequivocally determined: if in field studies the computation of this quantity undergoes the limitation of the instrumental arrangement (i.e., how above the canyon we have managed to install an anemometer), there is not a unique definition for laboratory and numerical studies.
A limited number of studies have used the convective-regime classification given by the Richardson number to investigate the buoyancy force and thermally-driven flows in urban canyons, with no definitive conclusions on their impact on the observed regimes. Given the flow inertia, the intensity and direction of the buoyancy force have determined an intensification/enlargement or weakening/shrink of the inertial vortex [28,29], or the formation of one/more counter-rotating vortex/vortexes both from field measurements [30] and numerical simulations [31]. However, in the same study by Louka et al. [31], a single recirculation vortex has been observed from field data, with the heating effect of the build-ing facade being confined in a thin layer close to the same building. This second outcome is in agreement with other studies (e.g., [14,[32][33][34]) which only observed an intensification of the upwards motion at the heated facade. Xie et al. [35] also found that an increasing aspect ratio leads to the formation of multiple vortexes within the canyon, in contradiction to Magnusson et al. [11] where the number of vortexes decreased to only one when introducing a buoyancy force.
The present study aims at providing further insight into the mixed-convection regime of a single infinitely-long urban canyon of unity aspect ratio. This study will delve into the thermo-fluid dynamics and turbulence of the in-canyon region and impacts on the pollutant removal, as obtained from a stationary background flow and varying buoyancy forces. Thus, we will not explore the full extent of the boundary-layer processes above it nor the variations in the background flow structure. Highly-resolved large-eddy simulations (LESs) are employed to reproduce the non-linear interaction among the in-canyon inertial velocities induced by the surrounding flow, the natural convection arising from heated building facades, and the mechanisms of pollution removal. Numerically, this required a careful choice of numerical schemes to realistically reproduce the coupling between inertial and thermal turbulence. Physically, the simultaneous presence of these elements generates complex thermo-fluid dynamics and turbulence features which are analysed in details. Also, the in-canyon flow regimes are investigated to further characterise the potential issues related to pollutant concentration removal and pedestrian exposure.
This paper is organised as follows: Section 2 describes the case study of this investigation, while Section 3 is devoted to the numerical simulation approach adopted. In Section 4 the simulation settings are presented leading to the simulation-approach validation in Section 5. Section 6 provides the major results of this investigation. Finally, Section 7 draws the conclusions.

Case Study Description
Numerical simulations reproduce the paradigmatic case of periodic arrays of infinitely long square canyons. The ambient wind is perpendicular to the canyon axis. A pollutant concentration is released from the bottom of the canyon with constant flux, while the facades of the buildings bordering the canyon heat up differently and generate natural convective motions (see Figure 1a). Two cases are analysed: (1) the favourable convection case, where the upwind facade is heated up and the other surfaces are at ambient temperature. In such a case, the convective velocity supports the velocity induced by ambient wind and the in-canyon circulation is strengthened. (2) The adverse convection case, where the downwind facade is heated up and the other surfaces are at ambient temperature. In this case, the convection flow opposes the ambient-induced velocity and the overall internal circulation is altered. The neutral case, where all the solid surfaces are at ambient temperature, is also reported as a reference. Figure 1 displays the case study geometry and computational domain dimensions. Two square buildings of height H = 1.0 m and width W = 1.0 m delimit a canyon with an aspect ratio H/W = 1. The computational domain consists of a rectangular region that is 2H wide (centred in the canyon), 3H height (starting from the bottom canyon level), and 2H deep along the canyon axis direction which is a direction of flow homogeneity. Since the attention is focused on the in-canyon dynamics, the height of the domain has been limited to 3H following the example of similar studies [36,37]. The origin of the system of reference is placed at the bottom left corner of the computational domain; the vertical, streamwise and spanwise direction are the z, x and y direction, respectively.
The ambient flow is perpendicular to the canyon axis, i.e., directed in the x-direction, and the average free-stream velocity at the top boundary is U 0 . The source of pollutant concentration c is at the street level, in a region 0.7H large and centred in the canyon. Such a region delimited the carriageway where motor vehicles release pollutant and do not extend the full width of the canyon to account for the presence of footpaths. The natural convection arises at the vertical building facade within the canyon. The buildings are at ambient temperature, while the temperature of the upwind facade (T 1 ) and downwind facade (T 2 ) are changed to produce two different horizontal thermal gradients. It should be noticed that in real case configurations, also the building roofs and the street can be heated up by solar radiation. Possibly, a mixed configuration where the rooftop, one building facade and the street are heated up, is also more common. The present study aims to analyse the effects of a horizontal thermal gradient, while additional setups can be investigated in future dedicated studies.

Non-Dimensional Parameters
The ambient wind fluid dynamics is characterised by the Reynolds number based on building height: where ν is the kinematic viscosity of air. The natural convection is described by the Grashof number, which is the ratio between the buoyancy force and the viscous force: The Richardson number can be interpreted also as the ratio of the buoyancy term to the flow shear term, and it is introduced to state the relative importance of the buoyancy force and the atmospheric wind: where U b = gβ∆TH is the characteristic buoyancy velocity. Additional parameters are the Prandtl number Pr = ν/α, the ratio between kinematic viscosity and thermal conductivity, and the Schmidt number Sc = ν/λ, the ratio between the kinematic viscosity and concentration diffusivity. Quantities are made non-dimensional by means of the buildings height H for length and the free stream velocity U 0 for velocities. The ambient flow time scale is defined as τ = H/U 0 . Temperature is normalised as T * = (T − T 0 )/∆T where T 0 is the ambient reference temperature and ∆T = |T 1 − T 2 | is the difference of temperature between the two building facade. Concentration is made non-dimensional as C * = c/C 0 and the reference concentration C 0 = q/U 0 is defined from the concentration flux at the street level q = λ ∂c ∂n , where λ is the molecular diffusivity and n is the street-normal direction (see also Section 4.3).

Simulation Approach
The thermo-fluid dynamics proprieties of air are summarised in Table 1. The modelling approach is detailed in the following sections.

Mathematical Model
The air is modelled as an incompressible fluid and the Boussinesq approximation is adopted to account for the buoyancy force. The governing equations read: where u i are the velocity components (i = x, y, z), p = (p/ρ 0 − g i x i ) is the deviation from the kinematic hydrostatic pressure, g i is the gravity acceleration vector, T is temperature, α is the molecular thermal conductivity, ρ 0 and T 0 are the reference density and temperature, respectively. The pollutant concentration is modelled as follows: with λ the molecular pollutant diffusivity. Being a passive scalar, c is just advected and diffused by the flow.

Large-Eddy Simulation Approach
In the LES approach, the scales of motion larger than the computational cell width are directly resolved, while the smallest scales are modelled (for example see Rodi et al. [38]). The numerical resolution of the governing equations over a discrete mesh acts as an implicit spatial filter (denoted with an overlying bar), which decomposes the velocity field in a filtered and a residual component: u i = u i + u i . Applying the filter to momentum Equation (6), an extra term appears and is modelled through the dynamic Lagrangian model proposed by Meneveau et al. [39]. Details on the mathematical model implemented can be found in Cintolesi et al. [40] and are not here reported for the sake of conciseness. Such a model is particularly suitable for reproducing localised and anisotropic turbulence because it can determine the sub-grid scale (SGS) contribution to the resolved flow at each point of the computational domain, extracting information from the instantaneous velocity field. Hence, this is convenient for simulating the present case, where the largest part of the turbulent motion is localised at the canyon-atmosphere interface. It was also shown that it can adjust the value of the SGS viscosity ν SGS near the solid boundaries, overcoming a well know shortcoming of the Smagorisnky model which uses damping function as the Van Driest function [41] to correctly reproduce the turbulent energy dissipation near walls [39]. Moreover, the empirical content and case-dependency of the model are reduced, since the value for the model parameter c 2 s is not specified as in the original Smagorinsky model but inherently computed.
A model for SGS contribution is also required for temperature and concentration. When the implicit spatial filter is applied to temperature Equation (7), a SGS thermal flux appears and is modelled by the flux-gradient hypothesis: where the last equation is the flux-gradient model and α SGS is the sub-grid scale thermal diffusivity. This is expressed using the Reynolds analogy as α SGS = ν SGS /Pr t , where Pr t = 0.7 is the turbulent Prandt number. In a very similar way, the filtering of the concentration Equation (8) leads to a SGS concentration flux and a SGS pollutant diffusivity that is computed as λ SGS = ν SGS /Sc t , with Sc t = 0.85 the turbulent Schmidt number. Additional details on the SGS model for active and passive scalar are reported in Cintolesi et al. [42,43].

Algorithm and Numerical Schemes
The numerical simulations are carried out using OpenFOAM [44], open-source software for computational fluid dynamics. The solver pimpleFoam has been customised to include the resolution of temperature and concentration Equations (7) and (8). The SGS model described in Section 3.2 has been originally implemented and tested in several previous studies [40,43,45].
The solution algorithm is the Pressure-Implicit with Splitting of Operators (PISO) as it is implemented in OpenFOAM (see the software documentation [44]). The buoyancy term in momentum Equation (6) is added as an explicit source term. The discretisation schemes are of the second-order accuracy: the implicit Euler backward scheme is used for time derivative, the central-difference scheme is used for spatial derivative. To stabilised the simulation, the convective term in temperature and concentration equations are discretised using the Gauss-Gamma scheme (γ = 0.2) proposed by Jasak et al. [46], which is a bounded central-difference scheme based on the normalised variable diagram. A numerical clipping is enforced on T and c to prevent nonphysical values given by numerical instabilities. The clipping acts on a limited number of computational cells during the initial development of the flow and once the statistically steady state is reached, it is practically not activated.

Case Settings
This section reports the numerical settings of the cases under consideration, along with the physical parameters adopted and compared with typical field-experiment values. Table 1 summarises the non-dimensional numbers that characterised the flow regime, together with the characteristic velocity and temperature difference, and the physical parameters adopted in this investigation. The Reynolds number based on the building high is set to be Re = 2 × 10 4 which is larger than the critical value ensuring the fluid-dynamic Reynolds independence for square urban canyons [47,48]. Thus, the velocity field (suitably scaled) reproduces the same configuration obtained at higher Re that are characteristic of real atmospheric conditions. Nevertheless, caution should be paid on this point since no studies have verified the validity of the criterion of independence of the motion from Reynolds number in the presence of natural convection. The temperature gradient imposed between the two vertical building facades is set to be ∆T = 2.5 K, leading to a Grashof number Gr = 8.17 × 10 8 , indicating that the convective flow along the vertical building surfaces is in a turbulent regime. The Richardson number is Ri = 2.05 suggesting that convective forces introduce non-negligible effects in the overall circulation.

Physical Parameters Set-Up in View of Field-Experiment Evidences
The physical parameters for this case study are derived because of evidence from several field experiments devoted to the investigation of convective forces in urban canyons. Table 2 summarises the typical values of Reynolds, Grashof and Richardson numbers in several urban canopies subjected to different free-stream velocities and temperature differences. Field studies encompass real cities and controlled environments (e.g., ad-hoc canopies build with containers or homogeneous concrete bodies positioned in unperturbed open spaces), with different mean building heights H and aspect ratios H/W. It is paramount to be aware that free-stream velocity and temperature differences are not equally determined in these studies due to the diverse instrumentation equipment and setup of the field experiments. The free-stream velocity U 0 is defined as the speed of the background flow unperturbed by the canopy. This measurement is not commonly available in an urban environment, where is typically replaced by the wind speed measured at the highest available instrumented level above the canyon, which is necessarily different in each selected study. Conversely, the identification of the temperature difference ∆T is a matter of various definition adopted by individual authors. The ∆T is estimated using the temperature of a heated surface (e.g., building facades or street) measured with thermocouples or thermo-cameras and a reference air temperature from thermocouples or probes. The major issue with the temperature-difference definition is the computation of the reference air temperature. Most studies [12][13][14]31,49] in Table 2 chose the air temperature measured above the canyon rooftop as the reference (at the same elevation used for measuring the free-stream velocity). The remaining studies [30,34,50] in Table 2 used the mean canyon temperature as the reference. These differences open up a discussion on the necessity to define a clear and unique Richardson number, suitable for real urban environments where the canopy and the surrounding atmosphere are rarely isothermal. The absence of consensus leads to ambiguities in the set-up of physical and non-dimensional parameters for studying urban canyon flows, limiting the research comparability. Compared to Table 2, the Richardson number used in the present study is within the typical range observed in field experiments, while Reynolds and Grashof numbers are smaller by orders of magnitudes, which is due to the reduced ambient wind velocity and building height compare to real case settings.

Computational Domain
The computational domain is discretised with an unstructured and staggered mesh, which is depicted in Figure 1b. The mesh consists of 1,592,562 cells; it is refined near the solid boundaries to ensure the direct resolution of the wall boundary layer and slightly stretched near the canyon-atmosphere interface (z/H = 1) where high turbulent motions are expected [40]. Preliminary simulations showed that the first computational point near the walls is placed within y + = u τ y/ν < 1 (note that to maintain the traditional notation, here y is the wall-normal coordinate and u τ is the friction velocity); thus, the wall boundary layer is well resolved [51]. Also, the adequacy of spatial discretization was checked: an additional simulation has been performed over a finer grid having 20% more cells in all directions. The first-and second-order statistics of velocity field along selected lines, computed using the present and the refined mesh, do not show differences. The time steps are dynamically computed to ensure the maximum Courant-Friedrichs-Lewy number CFL max ≤ 0.5, which preserves the numerical stability and accuracy (e.g., Trivellato and Raciti Castelli [52]).

Initial and Boundary Conditions
The flow is driven by a pressure gradient dynamically computed in a way that the average streamwise velocity at the top boundary is maintained at a fixed value U 0 . At the top boundary, the zero-gradient conditions are prescribed for velocity and temperature, while a zero value is set for pressure and concentration. These conditions mimic a free interface that allows the removal of pollutant and temperature from the domain. At solid walls, a non-slip condition is enforced for velocity and zero-gradient is imposed to pressure. The pollutant is released at the canyon street level (z/H = 0) in the area 0.65 < x/H < 1.35, corresponding to the carriageway for motorcycles. In this area, it is imposed a constant flux of concentration q = λ ∂c ∂n , where n is the street-normal direction. At the remaining solid boundaries, a zero-gradient condition is imposed. Regarding temperature, different constant values are set at the various walls: the horizontal boundaries are at the atmospheric temperature T 0 (building roof and canyon street); the vertical boundaries composed by the upwind (x/H = 0.5) and downwind (x/H = 1.5) building facade have a fixed temperature T 1 and T 2 respectively (see Figure 1). In the favourable convection case, the upwind facade is heated up at T 1 = T 0 + ∆T while the downwind facade is at ambient temperature T 2 = T 0 . Conversely, in the adverse convection case, the upwind facade is at ambient temperature T 1 = T 0 and the downwind facade is at T 2 = T 0 + ∆T. Lastly, an isothermal configuration T 1 = T 2 = T 0 is imposed in the neutral case. The computational domain is made virtually infinite in the y direction and periodic in the x direction, by imposing a cyclic boundary condition for all variables.
The initial conditions are constant weak streamwise velocity u x = 0.5U 0 m/s and zero values for others quantities. To accelerate the development of the flow dynamics, a preliminary simulation has been run to reach the statistical steady state for velocity and pressure, without resolving the temperature and concentration fields. Subsequently, the temperature and concentration have been transported as a passive scalar by an instantaneous configuration of the velocity field. Finally, the resolution of the complete model is activated and the coupling set of equations are resolved till a statistical steady state configuration is reached, after which the statistics are collected for a period of 60τ non-dimensional time.
Statistics are computed by averaging in time (over a period of 6τ) and along the spanwise direction (y, which is a direction of homogeneity). The average in time and space of a generic variable ψ is denoted by angular bracket ψ ; the fluctuations are denoted ψ = ψ − ψ ; the root-mean-square (RMS) are indicated and computed as ψ rms = ψ 2 .

Simulation Approach Validation
The accuracy of the simulation approach and the numerical solver has been assessed in Cintolesi et al. [40]. The neutral case, i.e., in absence of a temperature gradient and natural convection, has been reproduced and successfully validated against several experimental and numerical results from the literature. For the sake of completeness and to make this study as self-consistent as possible, the validation of the neutral case is here briefly reported and discussed. The simulation has been validated against laboratory experiments and one numerical simulation reported in the literature. Five datasets have been used for the comparison: Li et al. [36] and Chew et al. [48] perform a water-channel experiment considering seven canyons, collecting measurements at the central canyon. The former uses a two-colour laser Doppler anemometer, the latter uses an acoustic Doppler velocimeter with a sampling volume of 0.06 m 3 . Brown et al. [53] used a wind tunnel with six canyons and take the measurements with a pulsed-wire anemometer at the 3rd-4th canyons, where the flow reaches the equilibrium. Michioka et al. [47] perform a laboratory experiment in a wind tunnel with 24 canyons, collecting velocity measurements with a laser Doppler velocimeter. They also present an LES provided with the dynamic SGS model proposed by Germano [54] and Lilly [55]. Meroney et al. [56] and Pavageau and Schatzmann [57] provide concentration profiles obtained in laboratory experiments. Figure 2 shows first-and second-order statistics for the streamwise velocity and the mean concentration profile along the canyon vertical centerline (x/H = 1). To be consistent with the reference datasets, the velocity has been made non-dimensional through U Re f = 0.17 m/s, which is the mean streamwise velocity at the height of z/H = 2. The profiles of the mean streamwise velocity collapse to the same line, while the RMS profile slightly underestimates the resolved fluctuations predicted by the experiments but is in agreement with the simulation of Michioka et al. [47]. The peak of the resolved turbulent kinematic momentum fluxes u w at the canyon-atmosphere interface (z/H = 1) is lightly weaker, while the pollutant concentration within the canyon is satisfactorily reproduced. Given the validation performance, the simulation approach is considered adequate to reliably reproduce the case study under investigation.  [47]. Symbols are laboratory measurements: black •, Li et al. [36]; red , Brown et al. [53]; light green , experiment Michioka et al. [47]; blue +, Chew et al. [48]; dark green , Meroney et al. [56]; violet * , Pavageau and Schatzmann [57].

Results
The three introduced configurations are herein numerically reproduced and discussed. The neutral case where the building facades are at the same ambient temperature and therefore no natural convection is generated. The favourable case where the upwind building facade is heated up and the convective motion develops in the same direction of the inertially-driven motion of the principal vortex within the canyon. The adverse case where the downwind building facade is heated up, triggering a vertical convective velocity that opposes the motion of the main canyon vortex generated by the overlying ambient wind. See Section 4.3 for a complete description of the numerical settings.

Overall Canyon Dynamics
A first insight into the general circulation within the canyon is given through the spatial distribution of the principal variables (averaged and made non-dimensional) at a vertical plane y/H = 1. A quantitative analysis is reported in the subsequent Section 6.2.  (Figure 3c) shows a quite different configuration: the convective velocity at the downwind facade counteracts the downward flow from ambient wind by decreasing its intensity and deflecting its trajectory towards the canyon centre. The result is that the primary vortex is less energetic and confined in the upper-canyon (z/H > 0.2). The two recirculation regions located in the bottom corners of the canyon increase their extension and are displaced upward, while an area characterised by an almost zero velocity appears at the canyon bottom (z/H < 0.2), generating a zero-velocity region (see also discussion of Figure 5). The interface is a straight line starting at the upwind building top corner and intercepting the downwind building facade at z/H = 0.94. Figure 3d-f reports the pollutant concentration distribution and contour lines. All the cases have a maximum concentration near the street, where the source of the pollutant is located. The neutral and favourable cases (Figure 3d,e) exhibit an accumulation of concentration in the bottom-left corner, in correspondence of the recirculation region, and near the upwind facade as a consequence of the transport by the primary vortex. Instead, the adverse case (Figure 3f) accumulates the larger part of the concentration in the zerovelocity region, while decreasing with an almost vertical gradient. This is due to the thermal stratification which characterises the zero-velocity region, see Figure 4. Table 3 depicts the distribution of the pollutant concentration within the entire domain, the canyon (z/H < 1.0), and the zero-velocity region (z/H < 0.2), using the concentration density Ψ and the main recirculation quotient Φ. The concentration density for each region is computed as where V region is the volume of the region under consideration; the main recirculation quotient introduced by Göthe et al. [58] is defined as On average, the concentration density in the domain is higher in the neutral case which is less effective in dispersing the pollutant across the top boundary. Within the canyon, the favourable case shows lower values thanks to the enhanced internal recirculation, while in the zero-velocity region the adverse case exhibits a very high level of concentration. This zero-velocity region approximately coincides with the pedestrian level defined by Buccolieri et al. [59] which corresponds to about 1 to 2 m above the street, and where a high concentration of pollutants particularly affects human health. The main recirculation quotient shows that the neutral and adverse cases are equivalent in the removal of pollutant concentrations from the canyon, with the favourable case being more effective. Also, it points out that in the adverse case the pollutant accumulates near the street, whilst the upper part has a level of concentration similar to the favourable case. Table 3. Non-dimensional concentration density Ψ and the main recirculation quotient Φ within three selected regions: the entire domain (0 < z/H < 3), the canyon region (0 < z/H < 1.0), and the zero-velocity region (0 < z/H < 0.2). Percentage are computed with respect to the entire domain concentration.  The interaction region has a larger vertical extension, spanning within 1.2 < z/H < 0.8 and enlarging downstream. The negative peak is localised at the top corner of the downwind building, where the vertical convective flow perturbs the ambient wind. The descent-flow region is slightly detached from the downwind facade, reduced in the vertical direction (0.2 < z/H < 0.7), and directed along the anti-diagonal of the square canyon. Figure 4a,b displays the temperature distributions for the cases with convection. The favourable case (Figure 4a) shows that temperature is transported from the upwind facade to the interaction and descent-flow regions. The adverse case (Figure 4b) is characterised by a vertical stratification of the canyon. A strong thermal stratification is detectable in the zero-velocity region z/H < 0.2 near the street (see also discussion of Figure 3c), while in the upper-canyon (z/H > 0.2) the high temperature is localised near the downwind facade and weakly transported toward the upwind building by the primary vortex. When the warm air reaches the upwind cold facade, an additional weak vertical convection is generated, leading to an expansion of the left-bottom recirculation region along the facade (see Figure 3c). This contributes to limit the horizontal extent of the main vortex, displacing its centre downwind. Figure 4c-f presents the vertical and horizontal turbulent thermal fluxes. They can be interpreted through the eddy-viscosity hypothesis (e.g., see Pope [60]):

Entire Domain Canyon Region Zero-Velocity Region
where κ i is a positive constant, implying that negative (positive) turbulent flux increases (decreases) the temperature diffusion along the i−direction.

First-and Second-Order Statistics over Lines
The effects of favourable and adverse convection are quantitatively analysed comparing the profiles over selected vertical and horizontal lines. Figure 5 shows the mean non-dimensional vertical and horizontal velocity components along three vertical lines near the upwind building (x/H = 0.75), at the centerline (x/H = 1.00), and near the downwind building (x/H = 1.25). The streamwise velocity ( Figure 5 upper panel) reveals that both convective cases produce a higher streamwise velocity above the canyon (z/H > 1) compared to the neutral case. The latter shows an almost zero velocity at the rooftop level, and the development of a logarithm-like profile for the ambient wind. Such a behaviour can be ascribed to the strong horizontal shear at the roof level, which acts as a barrier between the canyon and ambient dynamics, mimicking the presence of a solid boundary. This effect is more evident at x/H = 0.75 while it is damped at x/H = 1.25 where the turbulence destabilises the sharp horizontal shear layer. The vertical convection within the canyon destroys the horizontal shear, increasing the turbulent transport of momentum and the mixing at the interface (see Figure 3); hence, the ambient wind velocity is more homogeneously distributed and reaches a constant value U 0 above the canyon. The favourable case shows a larger streamwise velocity above the canyon concerning the adverse case because it is more effective in vertical mixing, which is also amplified by the periodicity along the x-direction. Overall, turbulence reduces the drag of the urban canopy and allows more energetic wind velocities. Within the canyon (z/H < 1), the favourable case strengthens the circulation established by the principal vortex, increasing the flow intensity and generating more pronounced horizontal velocity peaks near the street and rooftop levels. On the contrary, the adverse case exhibits a very weak positive velocity in the zero-velocity region, whilst the upper-canyon (z/H > 0.2) is dominated by the (weakened and spatially constrained) principal vortex: negative velocity are detectable in 0.2 < z/H < 0.8, positive velocity for z/H > 0.8. The vertical velocity profiles ( Figure 5 bottom panel) highlight that in the favourable case the primary vortex velocity is reinforced if compared to the neutral case, and that a weak ascending motion is triggered also above the rooftop (see plot over lines x/H = 0.75, 1.25). The centerline (X/H = 1.00) reports a positive peak in the bottom half of the canyon which is due to the reduction of the bottom-right reticulation region (see Figure 3e). The adverse case displays almost zero velocity in the zero-velocity region and weak values in the upper-canyon coherently with the three-vortex circulation described in Section 6.1. Overall, the profiles of the vertical and horizontal velocity components of the neutral case show intermediate values between the favourable (higher values) and adverse (lower values) cases. This is mainly due to the alteration of the principal vortex, which is enhanced by favourable convection and hampered by adverse convection. x/H=0.75 x/H=1.00 x/H=1.25  Figure 6 shows the non-dimensional mean vertical velocities along three selected horizontal lines within the canyon. The adverse convection reverses the profiles. Near the downwind facade, the positive velocity is localised in a thin layer near the building surface. The descending flow from the principal vortex has a similar magnitude as in the neutral case at z/H = 0.75, while becomes weaker at the mid-height z/H = 0.5 and disappears at the lower level z/H = 0.25. Close to the upwind facade, a peak of weak negative velocity can be generated by the temperature gradient given by the thermal stratification of the canyon.  Figure 7 displays the temperature and concentration profiles along three vertical lines. Temperature profiles (top panel) inform that the favourable case is practically not stratified, having an almost constant temperature within and above the canyon. Exceptions are two narrow regions at the street and interface level. The peak of temperature at z/H = 1 is due to the transport of temperature by the main canyon vortex, which drives air from the heated facade to the interface region. Instead, the adverse case is strongly stratified in the zero-velocity region and exhibits a moderate vertical temperature gradient in the upper canyon. Outside the canyon, the stratification is almost absent. The overall temperature of the canyon is higher than in the favourable case because the internal air mixing enhances the heat exchange with the hot boundaries. Concentration profiles (bottom panel) shows that convection within the canyon largely improves the pollution removal: both convective cases decrease the concentration levels of the 50% at the canyon mid-high (z/H = 0.5) compare to the neutral case. This is ascribed to the large turbulent mixing triggered by convective motion at the canyon-ambient interface, which improves considerably the exchanges between the canyon and its surroundings. The profiles of the favourable and adverse cases are comparable, but some discrepancies are detectable: the former has a lower level of concentration near the street (lines x/H = 1.00, 1.25) and shows a weak peak at the interface as a consequence of the pollutant transport by the primary vortex. The latter has larger values of concentration in the zero-velocity region, which decreases smoothly with height. Near the upwind facade (x/H = 0.75) the concentration near the street of the favourable case is similar to the adverse case because of the corner recirculation region that entraps part of the pollutant.  in the zero-velocity region u u are almost zero and smoothly increase until the maximum located at the interface, which is slightly shifted downwards concerning the favourable case. This is caused by the deflection of the interface described in Figure 3. The adverse case displays larger fluctuations in the external region (z/H > 1) and a slightly high peak value at the interface. This can be ascribed to the stronger perturbation of the wind flow produced by the adverse convection.

Additional Thoughts on the Use of a Local Richardson Number
Based on numerical simulations and laboratory experiments reported in the literature, it is noted that a clear definition of Richardson number based only on the quantities characteristic of the canyon region is somehow missing. Although non unequivocally defined, the classical formulation of the Richardson number requires the use of a freestream velocity and a reference temperature, both computed or measured outside the canopy. On a more theoretical level, we can argue about the effectiveness of these nondimensional parameters in representing the in-canyon flow regime. Driven respectively by the free-stream velocity and characteristic temperature differences, the inertially-and thermally-driven circulations involve the air volume confined within the canyon. Thus, the definition of a local non-dimensional parameter seems more appropriate although rather complex. In this context, the use of building-building or surface-canyon air temperature differences seem suitable choices: the former addresses the local source of the buoyancy force but it can become of critical computation when the street is also heated; the latter directly infers the buoyancy effect of the canyon air temperature, a fair approximation if the canopy air temperature is isothermal away from the surfaces (e.g., Louka et al. [31]). As for the reference velocity, a characteristic canyon velocity would address the typical flow intensity of the investigated in-canyon circulation. The research of a local parameter can be considered in view of the confined canyon region where the buoyancy force is produced and the effects of the wind force are declined as flow inertia. Examples of local Richardson numbers have been given by Allegrini et al. [26] and Fernando et al. [61]. Both computed a surface-canyon air temperature difference based on an in-canyon air temperature and introduced a characteristic canyon velocity. Allegrini et al. [26] computed the in-canyon quantities in the closest proximity of the heated surface, defining a Richardson number representative of the direct (near-surface) consequences of the buoyancy force exerted by the heated building facade on the canopy flow. Fernando et al. [61] used an inertiallydriven canyon velocity determined as the free-stream velocity scaled on the canyon aspect ratio discerned from the continuity equation (see Magnusson et al. [11]). This definition directly addresses the in-canyon effects of the two driving forces.
In this context, the authors propose an alternative formulation for a local Richardson number as follows. The temperature difference ∆T = (T 1 − T 2 ) involves the characteristic (main) temperature of the upwind building facade T 1 and the downwind building facade T 2 . Notice that this quantity can be negative. Such a ∆T is representative of a configuration where the thermal gradient within the canyon is mealy horizontal, typically observed when the sun heats up a single facade while the opposite one remains in the shadow. Different dynamics arise in presence of vertical thermal gradient, that appears, for example, when the street or the ambient wind has a different temperature with respect to the canyon buildings. In such cases, a vertical thermal gradient tends to impose a stable or unstable temperature stratification within the canyon, with convection from vertical facades being triggered when the hot/cool air is transported near them. Hence, the thermo-dynamics regimes produced by horizontal and vertical temperature gradients are qualitatively different and should be characterised by two different non-dimensional numbers. Since we are interested in vertical interaction between convection and induced flows, the w component along the mid-high horizontal line (thus avoiding boundary effects from the canyon top and bottom) in Figure 6 is scrutinised as characteristic inertial velocity. Specifically, we compute the integral vertical velocity in the left-and right-half canyon as which are equal thanks to the continuity equation. It is worth noticing that this definition is general and not depending on the particular Re here considered, by virtue of the principle of Reynolds independence ensured by the case setting (see Section 4.1). As an alternative approximated formulation, w Re f can be estimated locally by measuring the maximum and the minimum of the vertical velocity at the mid-high and approximating the integral of the left-and right-half canyon with triangles. The average between the two triangles is a first approximation of the reference velocity: in the present case, this estimation gives w Re f ≈ 0.055U 0 . Using these local formulations of the temperature difference and characteristic canyon velocity, the canyon Richardson number becomes The Ri c has a positive sign if the thermal gradient is directed in the streamwise direction, otherwise, it is negative. Following the setup of this study, a positive Ri c indicates a regime of favourable convection, while a negative value indicates adverse convection.
The critical values to switch from an inertia-dominated regime to a convection-dominated regime are to be determined separately for the two cases, which is reasonable since the two canyon dynamics are qualitatively different. In the present case, Ri c = ±5 × 10 2 for favourable (positive sign) and adverse (negative sign) convection. It is worth noticing that the definition of Ri here proposed can be further developed to take into account different geometries (e.g., various aspect ratios) and the different heating of the canyon surfaces (e.g., including vertical thermal gradient and surface roughness effects) to better describe the complex configurations that can rise in realistic urban setups.

Conclusions
The present contribution analysed and quantified the effects that heated walls in urban canyons have on the flow structures and pollutant removal, which is relevant for practical applications. The large-eddy simulation (LES) approach is adopted to reproduce periodic arrays of infinitely long canyons with unity aspect ratio, where the ambient wind is at Re = 2 × 10 4 and vertical building facades have been set to different temperatures. A flux of pollutant is released at the street level. The dynamic Lagrangian turbulent model by Meneveau et al. [39] has been adopted to reproduce the localised and anisotropic turbulence. Two convective cases have been analysed: the favourable convection case where the upwind facade is heated-up, and the adverse convection case where the downwind facade is heated-up. The Richardson number is set to Ri = 2.05, being in the range of values measured in real urban canyons. A neutral case with isothermal buildings is also presented as a reference and for validating the numerical approach.
Both the convective cases show an increase in turbulence levels, quantified by the turbulent kinematic momentum and thermal fluxes. The significant increase of turbulence in the interaction region (around the canyon-ambient interface) enhances the vertical mixing and destroys the rooftop horizontal shear layer, increasing the exchange of momentum and mass between the canyon and the surrounding ambient. As a major consequence, the velocity intensity above the canyon exhibits an almost linear profile. Moreover, the pollutant removal from the canyon is largely increased with respect to the neutral case.
In the favourable case, convection strengthens the internal primary clockwise vortex that strongly advects temperature and pollutants along canyon boundaries. The in-canyon pollutant concentration is reduced by 49% compare to the neutral case. In the adverse case, convection counteracts the descending flow of the ambient wind at the downwind facade, strongly reducing the internal circulation. The primary vortex is less energetic and mainly confined in the upper-canyon (z/H > 0.2). A region of almost zero velocity generates above the street (z/H < 0.2), having strong thermal stratification and high pollutant concentration values: 40% more than in neutral case and 165% more than in the favourable case. In real settings, this is a critical region since corresponds to the pedestrian level. Also, a thermal stratification arises within the canyon, which further inhibits the internal recirculation. To characterised the different regimes reached by the two convective cases, a definition of a local Richardson number based on in-canyon quantities only (which is lacking in the literature) is proposed.