Numerical Simulations of Boundary-Layer Airflow Over Pitched-Roof Buildings

Arrays of buildings with pitched roofs are common in urban and suburban areas of European cities. Large-eddy simulations are performed to predict the boundary-layer flows over flat and pitched-roof cuboids to gain a greater understanding of the impact of pitched roofs on urban boundary layers. The simulation methodology is validated for an array of flat roof cuboids. Further simulations show that changes in the type of grid conformity have a negligible effect on the mean flow field and turbulent stresses, while having a visible, but small, effect on the dispersive stresses for a given packing density. Comparisons are made for flat and 45∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$45^\circ $$\end{document} pitched roof cuboid arrays at packing densities of 16.7% and 33.3%. The interactions between pitched-roof buildings and their effect on the urban boundary layer are considerably different to those of flat-roof buildings. The pitched roofs at a packing density of 33.3% leads to significant changes in the mean flow field, the Reynolds stresses, and the aerodynamic drag. Further work investigates the effects of changes in turbulence level and atmospheric thermal stratification in the approaching flow. Importantly, in comparison to a flat-roof array, the pitched-roof one at a packing density of 33.3% evidently increase the friction velocity and greatly reduces the effects of stable stratification conditions and changes in inflow turbulence level.


Introduction
Past research in urban atmospheric problems has shown that coupling may exist between atmospheric airflow at the city scale (O(10 km)), neighbourhood scale (O(1 km)), and street scale (O(0.1 km)) ( Barlow et al. 2017;Fernando 2010). This spatial coupling effect is par-ticularly strong around urban forms such as groups of tall buildings (Han et al. 2017;Fuka et al. 2018;Hertwig et al. 2019), or around steep topographic changes. To understand the importance of coupling effects we require a mapping between urban morphology and flowfield characteristics, such as mean velocity, turbulent statistics, momentum and heat fluxes. Acquiring such data for the varied range of geometries encountered in urban areas challenges both numerical and experimental simulations in many ways.

Buildings with Pitched Roofs
Most of the published experimental and numerical studies of idealized buildings (e.g. Cheng and Castro 2002;Stoesser et al. 2003;Coceal et al. 2006;Xie and Castro 2006) have focused on arrays of cuboid blocks with flat roofs, placed in an aligned or staggered arrangement with uniform or non-uniform heights. Such studies are most relevant to urban flows over heavily populated urban areas where most of the building blocks are likely to be cuboid in shape. However, in the surrounding urban and suburban areas of cities (e.g. European cities) which experience regular heavy precipitation, whether in the form of rain or snow, most of the residential houses have pitched roofs. For example, in 2008 it was recorded that 93% of the dwellings in England had pitched roofs (Department for Communities and Local Government 2008). It is therefore of great interest to assess the effect of having pitched, rather than flat roofs, on the local flowfield and large-scale boundary-layer flow.
There are a small number of published papers (Barlow et al. 2004;Yassin 2011;Ferrari et al. 2016;Nosek et al. 2016Nosek et al. , 2017Llaguno-Munitxa et al. 2017;Badas et al. 2017;Woodward et al. 2021) in which different roof shapes and their effect on the flow field and dispersion has been investigated. Most of these studies (Barlow et al. 2004;Yassin 2011) have investigated two-dimensional street canyons. A few studies (Holmes 1994;Tominaga et al. 2015;Ozmen et al. 2016;Fouad et al. 2018;Woodward et al. 2021) have also examined a single isolated building with a pitched roof and these have provided an understanding of the mean surface pressure around the building (Ginger and Letchford 1995;Oliveira and Younis 2000;Tominaga et al. 2015) which has supported the development of building regulation codes. Tominaga et al. (2015) investigated the flow field and surface pressure on an isolated building for one aspect ratio and three roof pitch angles using particle image velocimetry and computational fluid dynamics (CFD). They suggested that the flow pattern around a building with a pitched roof changes critically at a roof angle of around 20 • . Holmes (1994) found that a 30 • pitched roof on a single tropical house had a considerable effect on the mean roof pressures. The mean pressure coefficient on the upwind half of the roof was positive, while on the leeward side the flow did not re-attach on the roof, resulting in a near uniform surface pressure. On the other hand, Reardon and Holmes (1981) found that roof pitches up to 10 • in winds normal to the ridge yielded a small separation bubble at the leading edge with high shear-layer curvatures, which were associated with high suction pressure, and rapid pressure recovery downwind to the flow reattachment position. They suggested these low pitch roofs to be called "aerodynamically flat". These papers suggest that modest changes in the angle of pitched roofs may have strong effects on the flowfield around the building and building surface pressures. Fouad et al. (2018) investigated isolated pitched-roof buildings to explore the potential for obtaining useful design data using CFD methods by comparing their results to the Eurocode and ASCE10 standards. They concluded "the application of CFD techniques show great potential to offer very good wind design data for structures with shapes not listed in the existing codes." While building regulations are based on considering a single isolated building taking into account urban roughness, it is probable that the flow within the canopy of an array of pitched roof buildings will be considerably different to that around an isolated building, with some consequent effects on the surface pressure distribution. The limited studies on the effect of interference between buildings have focused on two buildings (e.g. Bailey and Vincent 1943;Holmes 2007). Bailey and Vincent (1943) suggested that having a maximum of two buildings upstream of the instrumented model was sufficient to be representative of built up areas. Holmes (2007) concluded that the building spacing was a key parameter. It is now well established that two building rows upstream of the measurement location is insufficient for a simulation of fully-developed internal boundary layer (IBL) over a very large array (e.g. Hanna et al. 2002;Sessa et al. 2018). For arrays of building blocks placed in a fully developed turbulent boundary layer in the scenarios typically studied (e.g. low-rise buildings and small roughness length upstream of the array), it requires a fetch length of about ten average building heights (e.g. Hanna et al. 2002;. For arrays of buildings placed in a smooth laminar boundary layer, the required fetch length is much greater. It is to be noted that while generating near equilibrium and fully developed urban boundary layer flows is necessary for reducing the uncertainty of the research data, such conditions are not generally established in real urban atmospheric airflows. Two parameters frequently used to characterise the morphology of urban areas are the frontal and plan solidities (e.g. Placidi and Ganapathisubramani (2015)). The latter is more often termed the 'packing density' λ p . The effect of packing density in arrays of simplified buildings with flat roofs has been intensively studied (e.g. Cheng et al. 2007;Ganapathisubramani 2015, 2018), and it has been shown that their total drag coefficient based on the freestream velocity is a function of the packing density. The drag coefficient function has minimums at λ p = 0 or 100%, and a maximum in the lower half of the packing density range. Cheng et al. (2007) showed that the drag coefficient of a staggered array of cubes was greater than for an aligned one at both λ p = 6.25% and 25%. The difference was greater at the smaller packing density. These results show that any study into the effect of pitched roofs should consider the effect of packing density.

Research Needs
The existing body of literature shows a clear need to investigate the effect of pitched roofs in building arrays at different frontal solidities, packing densities and aspect ratios with fetches sufficient to develop realistic urban boundary layers; to demonstrate improved methods for predicting the airflow and dispersion at neighbourhood scales; and to obtain spatially averaged momentum and scalar flux data which may be used to develop parametrizations for mesoscale models.
In addition to the above there is the question of thermal stratification effects in urban areas. These attracted little attention until one decade ago (Sessa et al. 2020). The studies published since then have focused on arrays of idealised cuboid-shape buildings (e.g. Boppana et al. 2014;Kanda and Yamao 2016;Tomas et al. 2016;Marucci et al. 2018;Marucci and Carpentieri 2020;Sessa et al. 2020), or real urban areas in which the building geometries have been simplified to have flat roofs (e.g. Xie et al. 2013). They have concluded that even weakly stable stratification conditions (i.e. the bulk Richardson number ≤ 1) have a considerable impact on the urban boundary-layer turbulence and significantly impacted pollutant dispersion.
To our knowledge there has been no research into the nature of the flow over arrays of pitched-roof buildings in thermally stratified conditions. The sparse literature which exists though shows that a pitched roof generates a very different flows within and above the canopy top in neutral stratification, with more three-dimensional flow structures than flows over flatroof buildings. The question arises as to whether the changes in thermal stratification lead to similar changes in the flow over pitched-roof buildings as those for flat roofs.
Rather than using conformal meshes to accurately capture the geometric details of buildings (such as pitched roofs), two widely used urban CFD codes PALM4U (Maronga et al. 2020;Krč et al. 2020) and uDALES (Heus et al. 2010;Grylls et al. 2020;Suter et al. 2021) use non-conformal Cartesian meshes. This simple treatment is attractive as it leads to significant improvements in computational efficiency, but raises the question of what the penalty in accuracy is compared to using conformal (body fitted) meshes? This question is also addressed by the study.

Outline of the Current Work
This study aims to address three knowledge gaps: (1) how the urban boundary layer is affected by having an array of buildings with pitched roofs, (2) how the effect of pitched roofs varies with the packing density of the array, and (3) how the effect of having pitched rather than flat roofs changes as the conditions vary from neutral to stable stratification.
Section 2 presents details of the numerical methods, Sect. 3 describes the cases and simulation settings. Section 4.1 details the verification and validation of the CFD model. Section 4.2 examines the effect of using conformal and non-conformal meshes, Sect. 5 presents a comparison between idealised block arrays with flat and pitched roofs, which includes spatially-averaged velocities, Reynolds stresses and surface pressures. Section 5.2 reports the effect of packing density in neutral stratification. Section 6 reports the atmospheric thermal stratification effects on airflow over pitched-roof buildings. Section 7 summarises the conclusions from the research undertaken.

Governing Equations
The flow which develops over arrays of blocks is innately unsteady and so is best resolved by adopting a large-eddy-simulation (LES) approach (e.g. Kanda et al. 2004;Xie and Castro 2006;Castro 2017). The LES models developed in this study were implemented was implemented within the open-source package OpenFOAM version 2.1.1. Eqs. 1 and 2 are respectively the filtered continuity and Navier-Stokes equations: where x i and t are respectively spatial coordinate and time, u i and p are respectively the resolved or filtered velocity and pressure, ρ is the density, ν is the kinematic viscosity, τ i j is the Subgrid-scale (SGS) stress. δ i1 is the Kronecker-delta and ∂ P/∂ x i is the pressure gradient or body force term which drives the flow when periodic boundary conditions are used in LES. The term f δ i3 is the body force in the vertical direction due to applied thermal buoyancy and is calculated by using the Boussinesq approximation based on the estimated temperature from Eq. 3. The mixed time scale SGS model (Inagaki et al. 2005) was used to avoid using the near-wall damping functions as used in the Smagorinsky SGS model. Nevertheless, Xie and Castro (2006) suggested that the flow is very much building block-scale dependent and that the results are not sensitive to the precise nature of the SGS model, subject to a requirement that the grid sufficiently resolves the inertial range of the turbulent spectra. The filtered transport equation of temperature is: where T is the resolved-scale temperature, D is the molecular diffusivity of temperature, D r is the SGS diffusivity given by ν r /Pr r , where ν r is the SGS kinematic viscosity, Pr r is the subgrid Prandtl number which is set to 0.9 (Xie et al. 2013).

Settings of Atmospheric Thermal Stratification
Atmospheric thermal stratification effects are set by using the bulk Richardson number. This offers a convenient approach for quantifying the strength of stratification due to the vertical temperature gradient in numerical simulations. The bulk Richardson number is defined as where g is acceleration due to gravity, δ is the domain height, T 0 is the mean ground temperature, T ∞ is the mean freestream temperature and U ∞ is the mean freestream wind speed.

Boundary Conditions
The LES model must be used with boundary conditions which are appropriate to the problem. For urban flow simulations the inflow and outflow faces typically have either periodic boundary conditions (PBCs) (e.g. Coceal et al. 2006) or synthetic turbulence inflow conditions (STI) (e.g. . The LES approach with PBCs is based on the assumption that the simulated domain is a repeat unit of a much greater region. This approach can substantially reduce the computational cost. However, when modelling a "real" urban geometry, the PBC approach could lead to an inaccurate representation of the flow in some scenarios, such as a change in roughness (Xie and Castro 2009;Tomas et al. 2016Tomas et al. , 2017Sessa et al. 2018Sessa et al. , 2020. Using the PBC approach in non-neutral conditions requires additional precautions, such as adding additional forcing terms into the governing equations that should not influence the turbulence statistics (Boppana et al. 2014;Grylls et al. 2020).
The STI method has been well documented Sessa et al. 2018) and simulations using it have compared well to data from experiments (Marucci et al. 2018;Marucci and Carpentieri 2020). The STI approach requires vertical profiles of mean velocity, Reynolds stresses, and integral length scales as inputs. To apply thermal stratification at the inlet STI requires vertical profiles of mean temperature, and the variance and integral length of temperature fluctuations as additional inputs.
We use the PBC and STI approaches for two different sets of simulation cases. The choice of boundary conditions was driven by the physics to be simulated, and so STI boundary conditions were used for examining the effect of thermal stratification. To enable the results to be compared with past work on block arrays with flat roofs, the same STI inputs were used for the neutral and stable stratification as in Sessa et al. (2018Sessa et al. ( , 2020) (specifically Fig. 4 in Sessa et al. 2020). This meant that the prescribed vertical profiles of mean velocity, Reynolds stresses, mean temperature and temperature fluctuations for simulating Ri = 0, 0.2, 0.5and1 were identical to those used in Sessa et al. (2018Sessa et al. ( , 2020, which matched the experimental data in Marucci et al. (2018). Moreover, the same prescribed inflow data were used for Ri = 0.2, 0.5 and 1, because as stated in Sessa et al. (2020), using the same settings allows the effect of thermal stratification to be isolated.

Spatial Averaging
The spatial averaging method was adopted for processing the outputs from all simulations. Using spatially averaging over a whole (x − y) plane (e.g. Fig. 3), a resolved instantaneous flow quantity φ in LES can be further decomposed into space-time average φ (Eq. 6), spatial variation of the time average Φ , and resolved turbulence fluctuation φ which is the deviation of the resolved instantaneous quantity φ from the time average φ, where the time-averaged quantity φ is equivalent to Φ for simplicity, such asū i being equivalent to U i . In particular, the comprehensive spatial average (Xie and Fuka 2018) was selected. This method includes the solid regions in the total volume, but ascribes zero values to the output quantities in them (i.e. the volume occupied by the building blocks shown as the grey area in Fig. 1d), while the solid regions within the domain are included in the total volume. Equation 6 defines the comprehensive spatial average: where denotes spatial average, φ denotes the quantity to be spatially-averaged, S f denotes the fluid area at height z. For the periodic boundary condition (PBC) case, S c denotes the total area at height z for the entire horizontal plane, i.e. x ∈ (0, 12H ), y ∈ (−6H , 6H ) (see Fig. 3a). For the synthetic turbulence inflow (STI) cases, S c denotes the area at height z for a horizontal plane across the entire span (y ∈ (−6H , 6H )), and between two x coordinates a and b separated by 2H (see Fig. 3b). The comprehensive spatial average has the advantage that it produces a smooth change of total momentum flux across the canopy interface . In comparison, the intrinsic spatial average, which only includes the fluid regions in the averaging, leads to a discontinuity in total momentum flux at the canopy height. The fact that the comprehensive spatial average gives a smooth change in momentum flux is potentially advantageous in developing a parametrization of momentum flux, or other quantities, for a meso-scale model.

Array Geometries
The idealized-morphology studied here takes the form of a simplified urban-like geometry. Arrays of simple blocks, i.e. cuboids with either flat or pitched roofs, are studied. The baseline cuboid has dimensions 2H × H × H , where H = 1 m. The Reynolds number based on the block width H and the velocity at height H in the upstream boundary layer is about 7400, or about 830 based on the friction velocity u τ (Coceal et al. 2006;Cheng and Castro 2002), estimated using an extrapolation from the linear turbulent shear stress profile above the canopy (e.g. Figure 6). Given the different block heights, and dimensionless domain heights in the work, we choose to use u τ to non-dimensionlise quantities, and to focus on the physics within and immediately above the canopy. Based on their research on arrays of cubes at a smaller Reynolds number, Xie and Castro (2006) state that "Reynolds number dependency (if it does exists) is very weak for such flows, except no doubt very close to solid walls. Turbulence generated by urban-like obstacles (with sharp edges) is large scale dominated and hardly at all dependent on the much-smaller-scale viscous dominated processes on the body surfaces". Flows over arrays of pitched-roof obstacles with sharp edges are expected very weakly dependent on the Reynolds number too.
Two different packing densities (i.e. λ p = 33.3% and 16.7%) are simulated. The baseline flat roof cases of height H were validated against Castro et al. (2017). Further cases of height 1.5H with a flat roof or a 45 • pitched roof, were then compared with the baseline case. Figure 1 shows the geometries of the three different blocks used and diagrams illustrating the packing densities.

Meshes and Test Cases
The meshes used in this study were created using SnappyHexMesh which is far more flexible and efficient than standard structured-mesh generators. It was primarily used here to assess its suitability for application to much larger domains in follow-on studies. Nevertheless, this study utilised the capability of SnappyHexmesh to create uniform grids, grids with multiple levels of resolution and non-conformal and conformal meshes.
For the baseline flat-roof cases uniform Cartesian grids (identified as (U), Table 1) were first created with a resolution of H /16. This was in accordance with the minimum level of resolution suggested in Xie and Castro (2006) and Coceal et al. (2006), which was adopted for all cases, except for those with specifically defined grids. Cartesian meshes with three levels of grid resolution (identified as (3R), Table 1) were created, as used in Xie and Castro (2006), for testing mesh sensitivity. These had a resolution H /16 up to z = 2.5H , a resolution H /8 between z = 2.5 H and z = 8.5 H , and finally a resolution H /4 from 8.5H up to the top of the domain at 12H . The time step was chosen following the guidance that the maximum Courant number was less than unit. For the PBC cases (Table 1), the initialisation period was about 600T p , and the averaging period was at least 700T p , where T p was based on the block height and the freestream velocity. For the STI cases (Table 1), the initialisation period was 50T p , and the averaging period was 200T p . It took from 240 to 300 wall-clock hours on 200 cores to complete one simulation case.
All the Cartesian meshes created for the flat roof cases were conformal. However, if a simple Cartesian mesh is created for a pitched roof, the mesh is non-conformal as shown in Fig. 2. To address the question of what impact this has, three meshing strategies were used for x coordinate is perpendicular to the long side of the cuboid and the long street, and z coordinate is in the vertical direction, d packing density λ p = 33.3% (as in Fig. 1), e packing density λ p = 16.7% where the width of the long street marked with "2" is tripled compared to "1" in (d) the pitched-roof cases: uniform conformal-grids (UC), Uniform non-conformal grids (UNC) and a 3-level resolution-conformal grid (3RC).
The majority of the cases were run at neutral stability conditions with PBC at the inflow and outflow (Table 1). For these cases the domain had dimensions of 12H × 12H × 12H , following Castro et al. (2017) as shown in Fig. 3a. For the cases with thermally stratified inflows, the STI boundary condition was applied (Table 1). The domain was sized to match that used in Sessa et al. (2018Sessa et al. ( , 2020, which was 29H × 12H × 12H shown in Fig. 3b. In both cases periodic boundary conditions were applied to the spanwise boundaries, and a stress free condition was set for the top boundary. Table 1 summarises the neutral stability cases simulated with the baseline packing density of 33.3%, the neutral cases run to examine the effect of changing the packing density from

Validation of Flow Simulation Around a Cuboid Array
To validate the methodology simulations were first made of cases 1 H(U) and 1 H(3R) with PBC at λ p = 33.3%. Figure 4 shows vertical profiles of the mean streamwise velocity component U , streamwise normal stress u u , vertical normal stress w w , and Reynolds shear stress −u w at the centre of the long street marked "1" in Fig. 1d. Overall, the differences in these turbulent statistics between the cases 1 H(U), 1 H(3R) and the LES data in Castro et al. (2017) are negligible within the canopy.
The mean streamwise velocity profiles sampled at typical locations (not shown) also agree well with those in Castro et al. (2017). The two cases 1 H(U) and 1 H(3R) yield accurate mean flow fields within the street canyon. The turbulent stress profiles of the case 1 H(U) are in agreement with those in Castro et al. (2017) over the entire vertical extent. Within the street canyon, the turbulent stresses for case 1 H(3R) are in reasonable agreement with those in Castro et al. (2017). Nevertheless, the discrepancy in the Reynolds stresses above z = 2.5H between the case 1 H(3R) and Castro et al. (2017) is visible, with its approximate maximum 5%. This difference is due to a greater proportion of the Reynolds stresses being modelled (unresolved) by the SGS model above z = 2.5H due to the coarser resolution for the case 1 H(3R). The modelled portion of the mixed time scale SGS model to the total Reynolds stress can be estimated (e.g. Xie et al. 2004). This is not of a primary concern of this paper. Despite the difference, the mean streamwise velocity is accurately predicted for the case 1 H(3R), again suggesting that the mixed time scale SGS model performs well for the coarser mesh.
In summary, the results suggest that the two conformal (body fitted) grids 1 H(U) and 1 H(3R) are able to produce satisfactory time-averaged velocities and turbulent statistics, e.g. Reynolds stresses. This provides the baseline of the mesh structure, from which the assessment of mesh sensitivity for the pitched-roof cases is carried out in Sect. 4.2.

Assessment of Conformal and Non-conformal Meshes
This section examines the effect of having conformal or non-conformal grids on the flow predictions obtained for a block array with pitched roofs. This is done by comparing spatially averaged profiles for the 1.5H45 • (UC), 1.5H45 • (UNC) and 1.5H45 • (3RC) cases at a packing density of 33.3% with PBC. Figure 5 shows that all three cases produced consistent vertical profiles of mean streamwise velocity and root-mean-square velocity fluctuation components. Overall agreement between the three cases is satisfactory, especially for the mean velocity profile in Fig. 5a. All cases predict similar peak root-mean-square velocity fluctuations with differences of less than 5% within and above the canopy. At the top of the domain, the increases in axial and particularly spanwise turbulence intensities are due to the so-called 'splatting' of the eddies on the top boundary, where the vertical velocity is constrained to be zero. This is not important, because as discussed in  the domain height is great enough to avoid any top boundary condition effect on the regions of interest, i.e. within and immediately above the canopy. It can be seen that at z/h = 2 where the interface between the H /16 and H /8 grid regions occurs, the axial fluctuating velocity root-mean-square peak in Fig. 5b and a Reynolds shear stress peak in Fig. 6a are visible. These only have a very local effect and do not evidently impact the flow regions within and immediately above the canopy. Figure 6 shows spatially averaged vertical Reynolds shear stress u w and dispersive shear stress −U W . There are negligible differences between the three cases within the canopy. The peak dispersive stress occurs just below the canopy top in all cases. The 1.5H45 • (UNC) case with non-conformal grid and 1.5H45 • (3RC) case with coarse grid above the canopy produce a small but visible discrepancy (less than 2% of u 2 τ ) from height z/h = 2 to z/h = 5, compared to the case 1.5H45 • (UC).
The discrepancy in the dispersive shear stress profile between the conformal and nonconformal grids is visible, with an approximate maximum of 3% of u 2 τ at z/h = 3. This small discrepancy is more likely owing to the uncertainty of the calculation of the dispersive shear stress than grid conformity. Span-width size structures above the canopy require significant averaging time for obtaining a 'converged' dispersive shear stress. These data confirm the satisfactory accuracy of using the non-conformal grids for the global quantities and those far away from the obstacles. Nevertheless, the accuracy of the near-wall quantities, e.g. surface pressure, is still a question for the use of a non-conformal grid, such as adopted by the PALM4U (Maronga et al. 2020) and uDALES (Heus et al. 2010;Suter et al. 2021) codes.
We speculate that the inaccurate representation of the local flow details in the vicinity of the pitched roof in the non-body fitted grid case 1.5H45 • (UNC) could affect the prediction of local skin friction and other flow parameters. Nevertheless, Figs. 5 and 6 suggest that the overall discrepancy between the spatially averaged mean velocity, Reynolds stresses and dispersive stresses due to using a non-conformal rather than a conformal grid is small. Figure 6a shows a linear decrease in Reynolds shear stress from the canopy top, reducing to zero at the domain top. The three profiles of the non-dimensional total shear stress (not shown) including the drag due to the blocks linearly increase within the canopy to almost 1 at the ground, which is expected because a constant body force ∂ P/∂ x is set to drive the flow. This confirms the satisfactory accuracy of the simulations (see Xie and Fuka 2018). The mesh resolution and numerical settings for the flat roof cases validated in Sect. 4.1, were kept the same for the pitched roof cases. All these ensure reliable results for the pitched roof cases.

Comparison of the Turbulence Statistics and Aerodynamics of Flat and 4• Pitched Roofs
This section compares flow fields, turbulence statistics, surface pressure coefficients and drag coefficients between arrays of cuboids with flat and 45 • pitched roofs. It aims to provide an understanding of how having a pitched roof affects these quantities and to thereby quantify the importance of accurately accounting for the impact of pitched roofs on the boundary-layer flow. Two packing densities λ p = 16.7 and 33.3% of uniform cuboid arrays with and without pitched roofs were simulated. Spatially-averaged mean velocities and Reynolds stress profiles were examined to understand the combined effects of a pitched roof and packing density.

Effect of Pitched Roof on Urban Canopy Flows
It is known that the upper boundary conditions with a large ratio of the boundary-layer thickness to the building height, such as δ/h ≈ 10, yields a negligible effect on the flow and turbulence within and above the canopy (e.g. Coceal et al. 2006;. It is worth identifying any differences of flow and turbulence over pitched-roof and flat-roof buildings in a thick boundary layer. Figure 7a shows a comparison of the spatially averaged vertical profiles of axial mean velocity between the pitched roof case 1.5H45 • (UC), and the flat roof cases 1.5 H(U) and 1 H(U), in which the vertical coordinate is normalized by the canopy height h, and the profiles of root-mean-square velocity fluctuation components are normalised by the friction velocity u τ . Within the canopy there are only small differences between the three cases, as they are of the same packing density λ p = 33.3%. The two flat roof cases show a much greater inflection of mean streamwise velocity at the canopy height (Fig. 7b), which is associated with an overall "smoother" canopy top resulting a weaker shear immediately above the canopy, compared to the pitched roof case. The pitched roof generates far more drag on the flow. This can be seen from the larger momentum deficit in Fig. 7a, and from the values of normalised friction velocity in Table 2, which show a significant increase for the pitched roof compared to the flat roof at λ p = 33.3%. Note a square of the normalised friction velocity (u 2 τ /U 2 3h ) is denoted "drag coefficient" (Boppana et al. 2014). Here the mean streamwise velocity U 3h at z = 3h is chosen as the reference velocity. Because the dimensionless domain height L z /h is different for different cases (see Fig. 15), the mean velocity at the domain top is not chosen as the reference velocity. Table 2 shows that the normalised friction velocity of the pitched-roof cuboid array is approximately 30% greater than the 1.5H flat roof case, and approximately 21% greater than the the 1H flat roof case. This strongly suggests that an improved urban canopy model by using more advanced parametrizations should take account of the pitched-roof effects.
We now examine the flowfields in detail to understand the reasons for the increase in friction velocity. Figure 7d shows far higher vertical velocity fluctuations at the canopy height for the pitched roof in the same freestream wind; it is to be noted again that the friction velocity u τ for the pitched roof case at λ p = 33.3% is significantly greater than the flat roof cases. This could be considered similar to a random height array , which produces greater drag and turbulence mixing as well (i.e. than uniform height blocks with flat roof). The axial, spanwise and vertical velocity fluctuation root-mean-square data are shown in Fig. 7b, c and d respectively. Within the canopy, the dimensionless root-mean-square velocity fluctuations are significantly different. The behaviour of the vertical profiles for the two flat roof cases are similar above the canopy, whereas the deeper canopy evidently suppresses the turbulent  Figure 8a, b for the flat roof cases show skimming flow regimes (Oke 1988), where the bulk of the flow does not enter the street canyon, and a stable circulation flow is formed in the canyon. A stagnation at z/h = 0.94 (Figs. 11c and 13a) confirms this. Fig. 8c for the pitched roof exhibits a flow in the wake interference regime (Oke 1988), where the downward flows is deflected down the windward roof of the next block downstream, at the stagnation of z = 0.85h (Figs. 11e and 13b). Dispersive stress within the canopy of the pitched-roof case 1.5H 45 • (UC) (Fig. 14) is much greater than that of the flat-roof case 1.5H (U ) in the same incoming flow. Again, this suggests that the former generates a less stable circulation in the canyon, and is in the wake interference regime, whereas the latter is in the skimming region. Figure 8c shows that the shear layer at the canopy top impinging on the windward pitched roof, which explains the greater vertical velocity fluctuations at the canopy height in the same freestream wind (Fig. 7d). Much greater vertical mean velocity magnitude (not shown) was also observed at the canopy height for the pitched roof case.
The observations above are reinforced by Fig. 9 which shows very different normalised instantaneous axial velocity contours for the 1.5 H flat and pitched roof cases. The flat roof creates an obvious interface at the canopy top, where a thin shear layer forms; whereas the pitched roof generates a far thicker shear layer from the apex of the roof, and a less visible interface. Vorticity contours (not shown) for both cases show the same phenomenon. Figure 10 shows how the pitched roof changes the spanwise flow around the cuboid. The development of the secondary flows evident in Fig. 10a, b has been well documented for uniform height flat-roof cuboids with various aspect ratios (Willingham et al. 2014;Vanderwel and Ganapathisubramani 2015;Tomas et al. 2017). However, the pitched roof effectively destroys any secondary flow (Fig. 10c). This is due to the large mean upwards flow at the apex of the roof, which is largely uniform across the span. The pitched roof also greatly enhances the axial rotation seen at the corner of the apex. Figure 11 shows pronounced differences in the windward and leeward surface pressure distributions normalised by ρu 2 τ on the three blocks. The windward side plots show similar pressure distributions for both the flat-roof cases, 1 H(U) and 1.5 H(U) (Fig. 11a, c respectively). The pitched-roof case, however, shows a considerably different surface pressure distribution. The vertical face from z = 0 to 2h/3 shows a similar low pressure region to the two flat roof cases, but with a lower pressure which extends much deeper into the canopy and further across the span.
On the leeward side similar distributions are again observed for the flat-roof cases (Fig. 11b,  d). But the leeward side of the pitched roof case (Fig. 11e) displays a much more uniform distribution across the span, due to the more uniform mean flow across the span shown in Fig. 10c. Figure 12 shows contours of pressure difference between the windward and leeward sides for the three cases 1 H(U), 1.5 H(U), and 1.5H45 • (UC). This figure shows a similar surface pressure difference distribution to that in Tsutsumi et al. (1992) with the highest volume ratio (packing density). The data collected also follows the same pattern in that, the surface pressure difference increases from low to high in going from the centre span to edge. The pressure difference peak occurs at y = 0 and almost the canopy top for the flat-roof block, whereas the pressure difference peaks on the pitched roof occurs at y/H = 0.15 and 1.85, and z/h = 0.9. The pressure difference on the pitched roof shows a high uniformity across the span, which probably results from a dominant flow along the roof surface towards the apex of the roof. Figure 13a shows vertical profiles of spanwise averaged surface pressure on windward and leeward sides. The profiles for two flat-roof cases are similar in shape; whereas the pitched roof creates significantly different shaped profiles above z = 2h/3. A noticeable inflection is present in the pitched roof profiles at the transition between side and roof, which occurs at slightly different vertical locations on windward and leeward sides, this is likely due to not having data at the transition point. The maximum surface pressure at z/h = 0.9 on the windward side, denotes the stagnation. Stagnation is not visible on the windward sides of the flat-roof cuboids, which occurs too close to the canopy top to be resolved in the current mesh. Figure 13b shows normalised spanwise-averaged pressure difference. A considerable difference is evident between the flat-roof case 1.5 H(U) and the pitched-roof case 1.5H45 • (UC) above z = 2h/3. The form drag due to the pressure difference between the windward and leeward faces is dominant. Xie and Castro (2006) show that the integrated pressure difference between the windward and leeward of an array of cubes is approximately 90% of the total body force imposed, which is equal to an integral of ρu 2 τ over the ground surface. The pressure difference for the pitched-roof case exhibits an inflection at z = 2h/3 as in Fig. 13a, accompanied by an abrupt increase in the pressure difference above z = 2h/3, which peaks at 0.9h and decreases towards the canopy top. The 1.5 H(U) case exhibits a lower pressure difference across the depth of the canopy, compared to the shorter cuboid case 1 H(U). A sharp increase occurs above z = 0.8h for case 1.5 H(U), and one occurs above z = 0.7h for case 1 H(U), which are consistent with the pressure profiles on the windward sides for the respective flatroof cases in Fig. 13b. These changes in shape are likely associated with the heights of the recirculation centres within the long street shown in Fig. 10a, b. Again, these suggest that a one-dimensional urban canopy model that resolves the aerodynamic drag distribution within the canopy should take account of the pitched-roof effects. Figure 14a shows dimensionless spatially-averaged Reynolds shear stress profiles. Within the canopy the pitched roof generates far less Reynolds shear stress. The profiles for the two flat roof cases are similar, with the lower cuboid case generating slightly more Reynolds shear stress within the canopy. When focusing on the flow immediately above and within the canopy, it is more appropriate to normalise quantities by the friction velocity. Table 2 shows more details of u τ /U 3h , where U 3h is the mean streamwise velocity at z = 3h. When focusing on the flow far above the canopy, it is more appropriate to normalise the flow quantities by the free stream velocity U ∞ , or the velocity high above the canopy e.g. z = 3h. We noticed that the peak Reynolds shear stress normalised by U 3h is about 100% higher for the pitched roof case (not shown), suggesting that the pitched roof vastly increases the drag, mixing and re-entrainment of the flow into the canopy. Within the canopy, Reynolds shear stress normalised by U 3h shows the same behaviour to that normalized by the friction velocity u τ . Considering only the packing density λ p = 33.3% against a flow region map based on flat-roof data, the flows fall into the the skimming region. However, the entrainment into the street canyon is less-energetic due to the pitched roof because the mean flow in the shear layer has a non-negligible upwards component and some portion of flow is directed upwards instead of into the street canyon. Taking account of this, the flow over the pitched-roof array with λ p = 33.3% falls into the wake interference regime.
The pitched roof creates a more uniform mean flow field across the span of the cuboid. The local peak Reynolds shear stress differs by less than 7% across the span of the cuboid at the long street location, i.e. station "1" in Fig. 1, with a maximum occurring at the center of the span. The flat roof cuboid generates the maximum Reynolds shear stress at 0.1H from the spanwise edge of the cuboid. Figure 14b shows spatially-averaged vertical dispersive stress profiles, which are the mean flow contributions to vertical momentum flux. The pitched roof creates significantly larger dispersive stresses at the canopy height compared to the flat roof -it is to be noted the friction velocity u τ is much greater for the pitched roof case. At z = 2 h and 1.5 h, the flat roof cuboids generate larger dispersive stresses, whereas the pitched roof generate negligible dispersive stress at the same height. We speculate that the pitched roof generates more turbulent fluctuations and mixing at the canopy height resulting in less significant mean flow variations. Table 2 shows an evident increase of normalised friction velocity u τ /U 3h for the pitched roof compared to the flat roof at packing density λ p = 33.3%, suggesting a significant increase of drag coefficient u 2 τ /U 2 3h . Nevertheless, Table 2 shows no evident difference in u τ /U 3h at packing density λ p = 16.7%, suggesting that the impact of the pitched roof is highly dependent on the packing density. Figure 15 shows spatially-averaged mean streamwise velocity and Reynolds shear stress profiles for the three block geometries at packing densities of 33.3% and 16.7%. To emphasise the variation of the canopy drag for different blocks in a given wind, we deliberately use U 3h to scale the quantities. With the exception of the 1.5H45 • (3RC) 16.7% case, the differences between the mean velocity profiles within and immediately above the canopy are hard to discern. Oke (1988) suggests that a threshold line between the wake interference and skimming regimes is h/W ≈ 0.65, where h/W is the canyon aspect ratio, and W is the canyon width. These data are based on a low ratio of building height to width, e.g. h/W = 1 (see Fig. 1, where W = H for λ p = 33.3%). Oke (1988) also suggests being prudent when using such a threshold, given the uncertainties. As discussed in Sect. 5.1, the flow over the flat roofed arrays at λ p = 33.3% and h/W = 1.5 is in the skimming regime, whereas at λ p = 16.7% (a) (b) Fig. 15 Effects of packing density and roof shape on a spatial-and time-averaged axial velocity profile, and b spatial-and time-averaged Reynolds shear stress and h/W = 0.725 it is in the wake interference regime. Again, these are supported by the flow pattern, surface pressure and skin friction data.

Packing Density Effects Accounting for Roof Shape
The peak values of Reynolds shear stress for the higher packing density cases are all lower than the peaks for the lower packing density.
It is crucial to note that the pitched roof increases the total drag coefficient substantially at the higher packing density λ p = 33.3%, whereas at the lower packing density λ p = 16.7% the effect of the roof on the drag coefficient is small. This confirms the observation in Sect. 5.1 that the pitched roof alters the flow regime at λ p = 33.3%, but it does not at λ p = 16.7%. This suggests that the flow regime map against the canyon aspect ratio needs to take account of roof shape, and requires a more systematic study of the issue, such as covering a certain range of λ p and h/W , but as a starting point, a specific case is of interest. Given the significant effect of the pitched roof, the width W r between one apex of pitched roof to the next apex is used to calculate the canyon aspect ratio h/W r . For the pitched roof at λ p = 33.3%, the aspect ratio h/W r = 0.725 falls in the wake interference regime. This is consistent with the h/W = 0.725 falling in the wake interference regime for flat-roof arrays at λ p = 16.7%. Figure 16 shows spatially-averaged mean streamwise velocity and Reynolds shear stress within and immediately above the canopy. The spatially-averaged mean velocity profiles are similar in shape, but the arrays with the lower packing density produce a smoother transition in the mean flow at the canopy top. As with the values above the canopy shown in Fig. 15, there is a clear difference in the values of Reynolds shear stress within the canopy for the two packing densities. It is noticeable though that the pitched roof reduces the Reynolds shear stress in most of the regions within the canopy, for both packing densities even though it substantially enhances the the Reynolds shear stress at the canopy height when λ p = 33.3%. We speculate that the 45 • slope on the windward side of the block convects most of the (a) (b) Fig. 16 Effects of packing density and and roof shape within the canopy on the a spatial-and time-averaged axial velocity profile and the b spatial-and time-averaged Reynolds shear stress turbulence upwards and downstream, whereas the straight vertical windward side of the flatroof blocks convects a large fraction of the turbulence downwards into a large circulation within the long street canyon as shown in Fig. 8.

Effects of Atmospheric Thermal Stratification
The study of flow over pitched roofs in neutral conditions as shown previously is a crucial stepping stone, albeit that the occurrence of neutral atmospheric stability conditions is very rare. In this section the effects of stable stratification are examined for pitched-roof cuboids. Again, it is to be noted that the same turbulent and temperature inflow statistics were used for the neutral and stable stratification as in the flat roof cases reported in Sessa et al. (2020). The imposed turbulent kinetic energy at the inlet in the neutrally stratified case is approximately twice that in the stably stratified cases. Data from the STI domain is shown as spatially averaged between a(17H ) − b(19H ) (Eq. 6). Figure 17 shows laterally averaged vertical mean velocity and Reynolds stress profiles at x = 18.5H for four stratified conditions. The mean velocity profile is negligibly affected by all the stratification conditions or the inflow turbulence level, whereas for flat-roof cuboids in Sessa et al. (2020) (i.e. Fig. 7) the effect was more evident. Figure 17 shows an evident effect of stratification and inflow turbulence level on the turbulent statistics. The reduction of streamwise normal stress (Fig. 17b) in going from Ri = 0.2 to Ri = 0.5 for the pitched roof case is approximately 7%, in contrast to a reduction of 25% for the flat-roof case (Sessa et al. 2020). For Ri = 0.2 and Ri = 1 the differences are 16% and 50% for the pitched and flat roof cases, respectively. These suggest the effect of a pitched roof opposes that typically associated with the stable stratification. The evident differences in the streamwise and vertical stresses between the neutral and the case Ri = 0.2 (in Fig. 17) are mainly due to the large difference of the inflow turbulence level. Nevertheless, these differences for the pitched roof cuboids are far smaller than for the flat-roof cases shown in Fig. 8 in Sessa et al. (2020). This suggests that the flow created around pitched roof cuboids leads to an evident reduction in the effect of inflow turbulence level, when compared to the flows over flat-roof cuboids for the same inflow conditions.
In neutral conditions the flat-roof generates around 10% stronger streamwise fluctuations than the pitched roof, whereas it generates around 10% less vertical fluctuations than the pitched roof (see Fig. 17, and Sessa et al. 2020). This suggests that the pitched roof generates more three-dimensional turbulent eddies above the canopy, as opposed to the thin shear layers generated by the flat roof at the canopy height. Thermal stratification normally inhibits vertical momentum transport while pitched roofs promote it by pushing flow upwards or downwards, so they have the opposite effect. The more energetic vertical turbulent fluctuations at the top of the pitched-roof cuboids reduce the effects of the stratification and the inflow turbulence level. Figure 18 shows contours of vertical stress w w across a vertical plane located at y = −1.5H (see Fig. 3) for Ri = 0, 0.2, 0.5 and 1. It is evident that the stable stratification reduces the strength of fluctuations across the entire fetch of the domain. Figure 18 also shows that the growth of the IBL is suppressed in the stratified flow. However, such a stratification effect is evidently less effective in suppressing the IBL when compared to its effect on the flow over flat-roof cuboids (Sessa et al. 2020), as we discussed in the preceding paragraphs. Sessa et al. (2018) developed an approach to defining the IBL interface by identifying the abrupt change in gradient of the vertical Reynolds stress profiles. Following the methodology of Sessa et al. (2018), we estimated the IBL depth for the pitched-roof array. Figure 19 shows spatially-averaged vertical normal stress profiles over an area of 2H (streamwise) × 12H (lateral) centred at x lee = 4H , 6H , …, and 22H , marked with 1, 2,…, and 10, respectively. The spatial average is defined in Eq. 6, where a = x lee − H , and b = x lee + H . The colour lines with symbols denote the IBL interfaces. It is evident that the thermal stratification suppresses the IBL growth. The IBL depth in Fig. 19 is evidently greater than that over a flat-roof array with the same input settings, including stratification condition Sessa et al. (2020). This confirms our early concluding remarks that compared to the flat roof, the pitched roof evidently enhances the mixing, significantly increases the aerodynamic drag, and reduces the thermal stratification effect.

Conclusions
Large-eddy simulations have been carried out to simulate the flow of urban boundary layers over idealized arrays of cuboids with and without pitched roofs, under neutral and stable stratification conditions. The reliability and accuracy of the results has been ensured by conducting rigorous evaluation tests, including examining mesh sensitivity and body conformity.
The results for a baseline flat roof cuboid array were first validated against previous data. The same mesh and numerical settings were then preserved as far as possible to ensure the reliability of flow simulations made for the 45 • pitched roof cases.
The non-conformal Cartesian meshes used in some CFD codes, such as PALM4U and uDALES, provide high computational efficiency, but their users might have doubts with Fig. 19 Spatially-averaged vertical normal stress profiles over an area of 2H (streamwise) × 12H (lateral) centred at x lee = 4H , 6H ,…, and 22H , marked with 1, 2,…,and 10, respectively. The coloured lines with symbols denote the IBL interfaces. a, e Ri = 0, b, f Ri = 0.2, c, g Ri = 0.5 and d, h Ri = 1 regard to the accuracy, such as in the near wall regions. This study has examined their accuracy of the global quantities, such as the spatially averaged mean velocities, second order turbulent statistics and dispersive stresses. We conclude that these meshes produce accurate spatially-averaged quantities compared to those from conformal (body-fitted) meshes.
Compared to that produced by an array of flat roof cuboids, introducing an array with 45 • pitched roofs leads to significant changes in the mean flow field, the Reynolds stresses, and drag, and increases the turbulent momentum flux at the canopy height by approximately 50%. The magnitude of these changes strongly suggests that it is important to account for the effect of roof shape in urban arrays.
The results have shown that the pitched roof case exhibits a different flow regime to that of the flat roof cases at packing density λ p = 33.3%. The pitched roof, due to decreased packing density in the top third of the canopy, is in the wake-interference regime, as opposed to the skimming flow regime experienced by the flat roof cuboids. At λ p = 16.7% both the flat-roof and pitched roof cases are in the wake-interference regime. The normalised spatially averaged profiles of Reynolds stress over the flat-and pitched-roof buildings show a significantly smaller difference when compared to the results at λ p = 33.3%. This suggests that the pitched roof only plays a crucial role within a certain range of packing densities.
Although pitched roof and random height arrays are more common than flat-roofed buildings of uniform height in urban areas, current parametrizations are typically based on arrays of uniform-height cuboids with flat roofs. The need for taking the effect of pitched roofs into account is given more importance by the fact that the results showed a much stronger effect for a packing density of around λ p = 33.3% which is more common in suburban regions than λ p = 16.7%. The results suggest that omitting these effects in a high-resolution meso-scale models may introduce non-negligible errors into boundary-layer predictions over typical suburban and urban regions, where there are a high proportion of pitched roof or nonuniform height buildings. Frontal solidity and packing density provide adequate descriptions of arrays for developing parametrizations for the effects of flat roofed buildings in meso-scale models, adding adding a further parameter to differentiate between arrays with different roof geometries would improve the performance of these models.
Simulations show that having pitched roofs means that the effect of increasing stable stratification conditions (0 ≤ Ri ≤ 1) is very much reduced in comparison to the same incoming flow for flat roof cuboids. The relative difference in the vertical fluctuations between two weak stratification conditions (Ri = 0.2, and 1) is around 16% less for the pitched roof cases, while a reduction of 50% is found for flat roof ones (Sessa et al. 2020). We conclude that pitched roofs may greatly enhance the turbulent mixing in stable stratification conditions, and so improve the urban ventilation in the local environment. Consequently, pollution models which ignore the effects of pitched roofs may be overly pessimistic.