A one-dimensional model of turbulent flow through “urban” canopies (MLUCM v2.0): updates based on large-eddy simulation

In mesoscale climate models, urban canopy flow is typically parameterized in terms of the horizontally averaged (1-D) flow and scalar transport, and these parameterizations can be informed by computational fluid dynamics (CFD) simulations of the urban climate at the microscale. Reynolds averaged Navier–Stokes simulation (RANS) models have previously been employed to derive vertical profiles of turbulent length scale and drag coefficient for such parameterization. However, there is substantial evidence that RANS models fall short in accurately representing turbulent flow fields in the urban roughness sublayer. When compared with more accurate flow modeling such as large-eddy simulations (LES), we observed that vertical profiles of turbulent kinetic energy and associated turbulent length scales obtained from RANS models are substantially smaller specifically in the urban canopy. Accordingly, using LES results, we revisited the urban canopy parameterizations employed in the one-dimensional model of turbulent flow through urban areas and updated the parameterization of turbulent length scale and drag coefficient. Additionally, we included the parameterization of the dispersive stress, previously neglected in the 1-D column model. For this objective, the PArallelized Large-Eddy Simulation Model (PALM) is used and a series of simulations in an idealized urban configuration with aligned and staggered arrays are considered. The plan area density (λ_p) is varied from 0.0625 to 0.44 to span a wide range of urban density (from sparsely developed to compact midrise neighborhoods, respectively). In order to ensure the accuracy of the simulation results, we rigorously evaluated the PALM results by comparing the vertical profiles of turbulent kinetic energy and Reynolds stresses with wind tunnel measurements, as well as other available LES and direct numerical simulation (DNS) studies. After implementing the updated drag coefficients and turbulent length scales in the 1-D model of urban canopy flow, we evaluated the results by (a) testing the 1-D model against the original LES results and demonstrating the differences in predictions between new (derived from LES) and old (derived from RANS) versions of the 1-D model, and (b) testing the 1-D model against LES results for a test case with realistic geometries. Results suggest a more accurate prediction of vertical turbulent exchange in urban canopies, which can consequently lead to an improved prediction of urban heat and pollutant dispersion at the mesoscale.

1 Introduction

Mesoscale meteorology is of particular interest for urban climate analysis: many weather phenomena that directly impact human activities occur at this scale, and the effects of urban roughness, heat, pollutant, and moisture on the atmospheric boundary layer (characterized as urban boundary layer) have important mesoscale implications. Accordingly, mesoscale modeling is a powerful tool for the analysis of urban climate and further prediction and management of urban heat and pollution.

In mesoscale models, urban climate variables on timescales of hours to days depend on multiple spatial scales from the street scale to synoptic scales. Given contemporary computational resources, however, it is not feasible to explicitly resolves building shapes (O(1–100 m)) and at the same time span a domain large enough to assess mesoscale impacts on the urban boundary layer (UBL; O(10–100 km)). Therefore, mesoscale models must parameterize the subgrid-scale exchanges of momentum, pollutant, moisture, and heat across the urban canopy layer (UCL) and UBL interface (Fig. 1).

These “subgrid”-scale urban processes may be classified as hydrodynamic (flow) or thermal (e.g., radiation, convection, conduction). In the case of the former (focus of this study), the flow near the surface is being treated with approaches of varying complexity. The simplest and oldest is the bulk transfer approach, with the Monin–Obukhov similarity theory (Monin and Obukhov, 1954) to account for varying atmospheric stability. However, this approach accounts only for surface–atmosphere exchange and the effects on the overlying atmosphere. Canopies (e.g., forests, urban neighborhoods) result in a new atmospheric layer of importance: the roughness sublayer (Rotach, 1993) or its subset, the canopy layer (Oke, 1976). It is the flow in this layer that directly impacts the wind, air temperature, and pollutant levels to which urban dwellers are exposed.

In the past 2 decades, urban canopy models (UCMs) have been developed to approximate the flow and thermal exchanges within and above neighborhoods and to couple with mesoscale models. Single-layer UCMs (Masson, 2000; Kusaka et al., 2001; Kanda et al., 2005; Bueno et al., 2013) have only one layer within the canopy and focus on the overall exchange of heat, momentum, and moisture with the overlying atmospheric model. Moreover, they typically parameterize the exchange of momentum using MOST (Monin–Obukhov similarity theory) and use simple empirical relations to diagnose canopy wind speed. Multi-layer UCMs (Martilli et al., 2002) have several layers within the canopy (Fig. 1) and permit a more process-based treatment of canopy physics. However, they are computationally expensive as they employ prognostic equations for both momentum and turbulent kinetic energy (TKE) solved with “urban canopy parameterization” or UCP (Martilli et al., 2002; Dupont et al., 2004; Santiago and Martilli, 2010). For instance, Santiago and Martilli (2010) presents a one-dimensional (column) model of vertical exchange of momentum and turbulent kinetic energy based on the k−l turbulence closure scheme (1.5 order). This model employs horizontally averaged microscale computational fluid dynamics (CFD) simulations based on Reynolds-averaged Navier–Stokes (RANS) with the standard k−ε turbulence model to determine required input parameters to the column model (drag coefficients and turbulent length scales as a function of height), and it is designed to predict the hydrodynamic component of multi-layer UCMs. Similarly, Simón-Moral et al. (2014) employed CFD simulations of idealized urban configurations and updated the parameterization of drag and turbulent length scale based on the horizontal heterogeneity caused by the variation in streamwise and spanwise streets. Krayenhoff et al. (2015) further extended the column model to include the effects of tree foliage on mean wind and turbulent kinetic energy in urban canopies. Subsequently, Krayenhoff et al. (2020) added temperature, humidity, and buoyancy effects to the Krayenhoff et al. (2015) flow model (Santiago et al., 2014) and combined it with previously developed models on radiation (Krayenhoff et al., 2014) and thermal (Martilli et al., 2002) balance for a comprehensive representation of trees at the street level.

The combined multi-layer urban canopy model, called BEP-Tree (Krayenhoff et al., 2020), is the first multi-layer (column) model of urban flow and energy exchange at the neighborhood scale that includes the radiation and dynamic effects of trees in the street canyon. However, when compared with detailed spatially sampled measurements over a 2 km by 2 km area in the Sunset neighborhood in Vancouver, Canada (Crawford and Christen, 2014), results indicated that the model strongly overestimated daytime air temperature. Krayenhoff (2014) concluded that the underestimation of vertical exchange of heat is what results in a significantly higher canopy air temperature calculated. Additionally, the study reported that large differences persist with or without trees and for several days of simulation; therefore, the underestimation cannot be attributed to the parameterization of trees or anomalies in the observations. Recent work by Krayenhoff et al. (2020) demonstrates that larger turbulent length scales (based on the current large-eddy simulations, LES, analysis) markedly improve pedestrian-level air temperature predictions compared to measurements.

In this study, we aim to investigate the factors contributing to the underestimation of vertical exchange of heat and momentum in the multi-layer column model. We hypothesize that the following factors may be responsible:

• RANS simulations as the basis for 1-D parameterization. Given the simplified assumption of the turbulent flow in the RANS models, it is likely that the turbulent length scales derived from the RANS-CFD model are a culprit.

• Contribution of dispersive stress. The dispersive stress has been neglected in the parameterization and formulation of multi-layer model, though it has been shown (Coceal et al., 2006) that it contributes to the total turbulent flux at urban canopy level, specifically for higher urban densities.

• Idealized versus realistic configurations. So far, the parameterizations are derived for the simplified “urban” arrays with uniform height, while mesoscale models aim to represent the impact of real urban neighborhoods.

• Thermal effects. It is possible that the simplistic representation of thermal effects on vertical turbulent heat transport further contributes to the underestimation of turbulent exchange.

Figure 2 Flow chart of the present study. Three-dimensional LES simulations are performed and vigorously tested in an idealized configuration of buildings. Then, using the spatially averaged profiles, urban canopy parameterizations of the multi-layer (column) model are revisited. The updated multi-layer (1-D) model is then evaluated against the UCPs in Santiago and Martilli (2010), the 3-D simulation results, as well as LES results in a realistic configuration by Giometto et al. (2017).

Considering that (a) the underestimation of vertical exchange of momentum is also seen in neutral cases and (b) the height variability in the Sunset neighborhood in Vancouver, which was used in Krayenhoff et al. (2020) for model evaluation, is relatively small, we focus on the first two factors in this analysis. Accordingly, for a more robust assessment of the urban canopy parameterization in the column model (focusing on the turbulent length scales, in particular), we employ a LES model for a more accurate representation of turbulent flow (Xie and Castro, 2009; Salim et al., 2011; Gousseau et al., 2011; Nazarian et al., 2018a, b) and aim to include the contribution of dispersive stress.

Figure 2 summarizes the structure of the present study. Sect. 2 describes the methodology to achieve the objectives in the 1-D multi-layer urban canopy model, and in Sect. 2.2.2, the LES model and setup are rigorously tested to ensure the fidelity of the results. Subsequently, drag coefficients and turbulent length scales are derived from LES of idealized arrays of cubes at varying densities for neutral conditions (Sect. 3.3 and 3.2). Finally, the column model with updated parameters is evaluated against horizontally averaged LES results in both idealized and realistic configurations (Sect. 3.4). Section 4 further summarizes the findings of this study and maps out future developments of the multi-layer model.

2 Numerical methods

2.1 One-dimensional column k−l model

The momentum equation in mesoscale models undergoes two averaging processes (Martilli and Santiago, 2007; Santiago and Martilli, 2010). First, the Reynolds decomposition is applied to the momentum equation such that the mean flow quantities are separated from fluctuating turbulent parameters (time or ensemble averaging, $u=\stackrel{\mathrm{‾}}{u}+u^{\prime }$). Second, quantities are spatially averaged over volumes that can be compared to a grid cell of a mesoscale model (horizontal-averaging, $\stackrel{\mathrm{‾}}{u}=\langle \stackrel{\mathrm{‾}}{u}\rangle +\stackrel{\mathrm{̃}}{u}$). Additionally, assuming (1) horizontal homogeneity (and hence, zero mean vertical velocity due to the assumed incompressibility), (2) negligible Coriolis effect, and (3) negligible buoyancy effects, the equation for the horizontal momentum is presented as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }\langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial t}}=-{\displaystyle \frac{\partial \mathit{\rho }\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle }{\partial z}}-{\displaystyle \frac{\partial \mathit{\rho }\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle }{\partial z}}-〈{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{P}}{\partial x_i}}〉+\mathit{\nu }\langle {\mathrm{\nabla }}^{\mathrm{2}}\stackrel{\mathrm{̃}}{u}\rangle ,\end{array}$$

where u and w are the streamwise and vertical velocity components, P is the pressure, and ρ is the air density (assumed to be constant here). In this equation and onward, 〈ψ〉 and $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ denote the spatial and time average of parameter ψ, respectively, and ψ^′ and $\stackrel{\mathrm{̃}}{\mathit{\psi }}$ are the departure of the instantaneous parameter ψ from the time or ensemble mean and the deviation of the mean quantity $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ from its spatial average, respectively (i.e., $\stackrel{\mathrm{̃}}{\mathit{\psi }}=\stackrel{\mathrm{‾}}{\mathit{\psi }}-\langle \stackrel{\mathrm{‾}}{\mathit{\psi }}\rangle$ and ${\mathit{\psi }}^{\prime }=\mathit{\psi }-\stackrel{\mathrm{‾}}{\mathit{\psi }}$). More information on the averaging techniques can be found in Martilli and Santiago (2007).

Accordingly, the first term on the right-hand side (RHS) of Eq. (1) is the spatial average of the time-averaged turbulent fluxes, while the second term is the dispersive stress (Raupach and Shaw, 1982; Martilli and Santiago, 2007), which accounts for the transport due to time-averaged structures smaller than the size of the averaging volume. Additionally, the third and fourth terms indicate the spatially averaged acceleration due to the pressure gradient, as well as the spatial average of dispersive viscous dissipation (viscous drag), respectively.

To parameterize the contribution of the spatially averaged turbulent momentum flux (first RHS term in Eq. 1), a K-theory approach is used:

$$\begin{array}{}\text{(2)}& \langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle =-K_{\text{m}}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial z}},\end{array}$$

where K_m is the diffusion coefficient for momentum using a k−l closure (Martilli et al., 2002) as

$$\begin{array}{}\text{(3)}& K_{\text{m}}=C_k{l}_k\langle \stackrel{\mathrm{‾}}{k}{\rangle }^{\mathrm{1}/\mathrm{2}},\end{array}$$

where C_k is a model constant for momentum, l_k is a length scale, and k is the TKE. C_kl_k is parameterized in the column model based on the CFD results (further detailed in Sect. 3.3).

Figure 3 Plan view of configurations used for LES analyses, representing “staggered” (a) and “aligned” (b) arrays of buildings. Note that computational domains consist of 5×3 and 6×3 (N_x by N_y) buildings for the aligned and staggered configurations, respectively, and only a subsection of this domain is shown here.

To calculate the spatially averaged TKE in Eq. (3), a prognostic equation is then solved where the same assumptions as Eq. (1) are made. The resulting equation is

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{k}\rangle }{\partial t}}=& -\langle \stackrel{\mathrm{‾}}{{u^{\prime }}_i{u^{\prime }}_j}\rangle {\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{u_i}\rangle }{\partial x_j}}-{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{k^{\prime }{w}^{\prime }}\rangle }{\partial z}}-{\displaystyle \frac{\partial \langle \stackrel{\mathrm{̃}}{k}\stackrel{\mathrm{̃}}{w}\rangle }{\partial z}}\\ & -〈\stackrel{\mathrm{̃}}{{u^{\prime }}_i{u^{\prime }}_j}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_i}}{\partial x_j}}〉-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{p^{\prime }{u^{\prime }}_i}\rangle }{\partial x_i}}-\langle \stackrel{\mathrm{‾}}{\mathit{\epsilon }}\rangle .\end{array}\end{array}$$

By (a) parameterizing the shear production $\left(-\langle \stackrel{\mathrm{‾}}{{u^{\prime }}_i{u^{\prime }}_j}\rangle \frac{\partial \langle \stackrel{\mathrm{‾}}{u_i}\rangle }{\partial x_j}\right)$ and turbulent transport terms $\left(-\frac{\partial \langle \stackrel{\mathrm{‾}}{k^{\prime }{w}^{\prime }}\rangle }{\partial z}\right)$ with K theory (Eq. 3) and (b) assuming K_m to be same for TKE and momentum (i.e., $K_{\text{m}}=-\frac{\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }\rangle }}{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle /\partial z}=-\frac{\langle \stackrel{\mathrm{‾}}{k^{\prime }{w}^{\prime }\rangle }}{\partial \langle \stackrel{\mathrm{‾}}{k}\rangle /\partial z}$), the TKE equation in the 1-D model is described as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{k}\rangle }{\partial t}}=& K_{\text{m}}\left[{\left({\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{v}\rangle }{\partial z}}\right)}^{\mathrm{2}}\right]+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\text{m}}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{k}\rangle }{\partial z}}\right)\\ & -{\displaystyle \frac{\partial \langle \stackrel{\mathrm{̃}}{k}\stackrel{\mathrm{̃}}{w}\rangle }{\partial z}}-〈\stackrel{\mathrm{̃}}{{u^{\prime }}_i{u^{\prime }}_j}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_i}}{\partial x_j}}〉+D_k-\langle \stackrel{\mathrm{‾}}{\mathit{\epsilon }}\rangle ,\end{array}\end{array}$$

where D_k is the source of $\langle \stackrel{\mathrm{‾}}{k}\rangle$ generated through the interaction with the buildings and the air flow and $\langle \stackrel{\mathrm{‾}}{\mathit{\epsilon }}\rangle$ is the viscous dissipation rate computed as

$$\begin{array}{}\text{(6)}& \langle \stackrel{\mathrm{‾}}{\mathit{\epsilon }}\rangle =C_{\mathit{\epsilon }}{\displaystyle \frac{\langle \stackrel{\mathrm{‾}}{k}{\rangle }^{\mathrm{3}/\mathrm{2}}}{l_{\mathit{\epsilon }}}}.\end{array}$$

C_ε and l_ε here are the model constant and the length scale of dissipation, respectively. In Santiago and Martilli (2010), l_ε in Eq. (6) is derived from the CFD-modeled turbulent kinetic energy and dissipation (Eq. 6), and using the RANS model constant for turbulent viscosity (C_μ), the turbulent length scale (l_k) is calculated as

$$\begin{array}{}\text{(7)}& C_k{l}_k=C_{\mathit{\mu }}{\displaystyle \frac{l_{\mathit{\epsilon }}}{C_{\mathit{\epsilon }}}}.\end{array}$$

Accordingly, to solve prognostic Eqs. (1) and (5), two main parameterizations should be provided. First, the turbulent length scales (and consequently the dissipation length scales) are parameterized based on the CFD results of $\langle \stackrel{\mathrm{‾}}{k}\rangle$, $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle$, and $\frac{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial z}$ at different heights in the UCL (Eqs. 2 and 3). Second, the drag term due to buildings is parameterized as follows. In the momentum equation (Eq. 1), the drag force is introduced as a sink of momentum, given that buildings are not explicitly resolved and the averaging air volume is not connected (i.e., containing porosities representing the volume of the buildings). Accordingly, the drag at height z is parameterized (Santiago and Martilli, 2010):

$$\begin{array}{}\text{(8)}& {\left.{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}〈{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{P}}{\partial x}}〉\right|}_z=D=S\left(z\right)C_{\text{d}}\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle |\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle |,\end{array}$$

where S(z) is sectional building area density (square meters of area facing the wind per cubic meter of outdoor air volume), $\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle$ is the spatially averaged mean wind speed, and C_d is the sectional drag coefficient for buildings.

Additionally, analogous to the momentum equation, the source term of TKE due to the conversion of mean kinetic energy into turbulent kinetic energy by the presence of buildings is parameterized as

$$\begin{array}{}\text{(9)}& -〈\stackrel{\mathrm{̃}}{{u^{\prime }}_i{u^{\prime }}_j}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_i}}{\partial x_j}}〉+D_k=S\left(z\right)C_{\text{d}}|\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle {|}^{\mathrm{3}}.\end{array}$$

Similarly, the parameterization of drag induced by tree foliage and the interaction with the buildings can be considered, detailed in Krayenhoff et al. (2015). Using the LES results, we revisit the parameterization of length scales and drag coefficient induced by buildings and discuss the consideration of dispersive stresses in Sect. 3.3 and 3.2.

2.2 Three-dimensional large-eddy simulation model

The LES results are used as a superior method to RANS models for evaluating turbulence characteristics and dispersion behavior in urban canopies (Xie and Castro, 2006; Salim et al., 2011; Nazarian and Kleissl, 2016). A PArallelized Large-Eddy Simulation Model (Raasch and Schröter, 2001; Letzel et al., 2008; Maronga et al., 2015) is employed here, which solves the following: filtered incompressible Boussinesq equations, the first law of thermodynamics, passive scalar equation, and the equation for subgrid-scale (SGS) turbulent kinetic energy. The subgrid-scale fluxes are parameterized using the 1.5-order Deardorff flux–gradient relationships (Deardorff, 1980), which use the SGS-TKE equation to calculate eddy viscosity. The Temperton algorithm for the fast Fourier transform (FFT) is also used to solve the Poisson equation for the perturbation pressure. A more detailed description of PALM can be found in Maronga et al. (2015).

2.2.1 Setup of LES simulations

A series of neutral simulations is considered for idealized urban-like configurations with aligned (Fig. 3a) and staggered (Fig. 3b) arrays of identical cubes. The plan area density (${\mathit{\lambda }}_{\text{p}}=A_{\text{p}}/A_{\text{T}}$) is varied from 0.0625 to 0.44 for both configurations to span a wide range of urban densities (from sparsely developed to compact midrise neighborhoods, respectively), where A_p and A_T are the plan area and total area of roughness elements, respectively. Similar to Santiago and Martilli (2010), the obstacles are cubes, such that λ_p=λ_f. Total height in the simulation domain is 7.4H, where H is the building height (16 m), and the wind direction is in the x direction and perpendicular to the array (Fig. 3). The simulations are performed over arrays of 5×3 and 6×3 (N_x by N_y) buildings for the aligned and staggered configurations, respectively. In all simulations, the canyon height is resolved by 32 grids and the same uniform grid resolution is used in x and y directions (0.0312H). In the vertical direction (z), a uniform grid resolution is used up to 4H and grid spacing is gradually increased thereafter. Periodic boundary conditions are employed in horizontal directions (x and y axes) to simulate an infinite array. The specifications of the geometric configuration, domain height, and grid resolution are motivated by detailed sensitivity analyses in Yaghoobian et al. (2014) and Nazarian et al. (2018a) to ensure that the large eddies influencing the canopy flow are resolved. Additionally, Sect. 2.2.2 further discusses the validity of simulation setups for the parameterization of canopy flow.

The flow is driven by a pressure gradient of magnitude $\mathit{\tau }=\mathit{\rho }{u}_{\mathit{\tau }}^{\mathrm{2}}/H_T$, where u_τ is the total wall friction velocity and H_T is the total domain height (7.4H). The corresponding u_τ is ≈0.21 m s⁻¹, which results in ${\mathit{Re}}_H=UH/\mathit{\nu }\approx {\mathrm{10}}^{\mathrm{6}}$. We note that calculating u_τ using the surface kinematic momentum fluxes in the horizontal directions (i.e., $u_{\mathit{\tau }}=({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}$) yields the same value. For the top boundary condition for momentum, a zero-gradient (free-slip) boundary condition that enforces a parallel flow is used.

2.2.2 Model evaluation: validity of the simulation setups

The PALM model is widely used and has been validated against various experimental measurements (Maronga et al., 2015; Park et al., 2012; Yaghoobian et al., 2014; Nazarian et al., 2018b). However, since the parameterization of the multi-layer model requires a high accuracy of results for the turbulent flow characteristics, we extended our analysis to evaluate the validity of the simulation setups for this study and further compare the results with the RANS model (Sect. 2.2.3). First, we compared the profiles of turbulent kinetic energy and Reynolds stresses with wind tunnel experiment data as well as other available LES studies. Our velocity and Reynolds stress showed good agreement when compared with the LES results of Kanda et al. (2004) (not shown), and the quadrant analysis showed good agreement of the flow structures and coherent structures when LES was compared with the direct numerical simulations of Coceal et al. (2007) in Nazarian et al. (2018b). Lastly, we compared the TKE profiles obtained with the LES results with the wind tunnel experiment of Brown et al. (2001) for a 3-D building array with aligned configurations and observed good agreement in the shape of the profiles and TKE above the canyon, while underestimation of TKE within the building levels is seen (Fig. 4). Such underestimation of TKE compared to measurements in the canopy was also reported in Giometto et al. (2016) for a realistic urban configuration. Additionally, since the exact value of friction velocity was not available in the experimental dataset, the velocity at 3H is used for this comparison, which may further contribute to the discrepancy. A direct comparison between LES and RANS demonstrates that RANS underestimates TKE even further compared to the wind tunnel results (Sect. 2.2.3).

Figure 4 Comparison of the TKE profile at the center of the canyon with experimental results of Brown et al. (2001) for a 3-D building array with aligned configuration (11×7 obstacles). The aspect ratio of the wind tunnel experiment and numerical simulations are set to 1 ($H/W=\mathrm{1}$), resulting in the skimming flow regime (Oke, 2002). The domain height in the numerical simulations was set to 8H to be compatible with the experimental setup as well as the numerical results of Santiago et al. (2007). Vertical profiles along the center line of the last three street canyons (indicated by M, O, and Q here) are compared with the ensemble-averaged canyon in the LES simulations. More information regarding the experiment configuration and comparison with numerical results can be found in Brown et al. (2001) and Santiago et al. (2007).

Second, in order to ensure the accuracy of our LES analysis, the choice of simulation setups is rigorously evaluated here, and a series of sensitivity analyses are performed to compare the profiles (time- and ensemble-averaged) of mean flow, TKE, and velocity covariances based on the (1) geometrical configuration (size and height of the domain), (2) grid resolution, and (3) run time parameters (spin-up time, sampling frequency, and time-averaging interval).

We find the domain height to be critical for both staggered and aligned arrays. The domain size of 4H previously used in the RANS simulations of Santiago and Martilli (2010) is insufficient for the present LES analysis, as it modifies the vertical profile of Reynolds stresses and accordingly TKE above the building height. Therefore, following a series of sensitivity analyses, 7.4H is used to ensure that the top boundary condition (i.e., the lack of solution for the entire boundary layer) minimally affects simulation results in the roughness sublayer. Similarly, the choice of domain size (number of buildings in an array) is critical. In the streamwise direction, a sufficient number of buildings should be included in the computational domain to resolve the large eddies influencing the canopy flow (Inagaki et al., 2012; Coceal et al., 2006). Similarly, Yaghoobian et al. (2014) compared 3×3, 5×3, and 5×5 arrays of aligned buildings and found 5×3 to be the best compromise between accuracy and computational cost. In this analysis, we extended the domain to an array of 6×4 aligned cubes and found insignificant differences in the vertical profile of turbulent parameters. Therefore, 5×3 and 6×3 arrays of cubes are selected for aligned and staggered configurations, respectively. Additionally, the grid resolution of 32 grid cells per H is used following the grid sensitivity analysis done by Yaghoobian et al. (2014) and Nazarian et al. (2018b) that showed lower grid resolution (such as 0.05H or 20 grid cells per canopy height H) to be insufficient for resolving the wall flow.

Regarding the run time calculations, three main parameters are evaluated. First the volume-averaged results are monitored throughout the runs, and the spin-up time (i.e., the initial time interval that is discarded in the subsequent analysis) is chosen to be 3 h, corresponding to 125 eddy turnover time ($T=H/u_{\mathit{\tau }}$). This initial time interval is necessary to reach quasi-steady behavior in TKE, velocity variances, and friction velocity at the surfaces. Second, the choice of sampling frequency is evaluated by comparing the vertical profiles of TKE and Reynolds stresses for 10, 20, and 50 time step sampling frequencies (time step size is 2 s). TKE results are influenced when a low frequency (50 time steps) is used. However, there is no significant change in the TKE profile between 10 and 20 time steps, though the computational cost is affected. Hence, results are saved every 20 time steps. The last and most important factor in the run time parameters is the time-averaging intervals. Coceal et al. (2006) and Nazarian et al. (2018b) have shown that in order to filter the formation of roll-like circulations over the urban-like configurations, the results should be averaged over a time period of 200–400 eddy turnovers. Therefore, in this analysis, the results are averaged over 6 h, which corresponds to 250T.

2.2.3 Model evaluation: comparison with RANS

Santiago and Martilli (2010) used vertical profiles of flow properties calculated from Reynolds-averaged Navier–Stokes (RANS) simulations of idealized arrays of buildings (setups similar to Sect. 2.2.1) to parameterize the 1-D column model (Sect. 2.1). The RANS model used for the urban canopy parameterization showed good agreement when evaluated against direct numerical simulation (DNS) and wind tunnel results for flow over aligned cube arrays (Santiago et al., 2008; Simón-Moral et al., 2014) and wind tunnel results for canopies of “vegetation” (Santiago et al., 2013b; Krayenhoff et al., 2015). When compared to the large-eddy simulation results (Fig. 5), the streamwise velocity as well as Reynolds stress at the building height calculated in RANS shows agreement with the LES results. For the vertical profile of Reynolds stress ($\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle$), the aligned configuration results in a better agreement within the canyon, while the above-canopy results are mainly dominated by the domain height (which is set at 4H in RANS, significantly lower than 7.4H in LES). However, when the vertical variation in normalized turbulent kinetic energy ($\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$) is compared to wind tunnel experiments (not shown) and LES (Fig. 5), RANS substantially underestimates turbulent kinetic energy in the urban canopy layer. Similar behavior is previously reported when the distribution of TKE obtained by LES and RANS k−ε models are compared with the measurements in a realistic urban configurations (Antoniou et al., 2017) and wind tunnel experiments (Xie and Castro, 2006). Additionally, Krayenhoff et al. (2020) suggested that the 1-D model of Santiago and Martilli (2010) contributed to underestimation of the venting in UCL, and the discrepancy has been traced to turbulent length scales derived from the RANS simulations. These new findings, and the recent advancements in the high-performance computing, motivate a revisit of these parameterizations with a more accurate flow model such as LES that has been shown to be superior in representing the turbulent flow statistics (Xie and Castro, 2009; Salim et al., 2011; Gousseau et al., 2011).

Figure 5 Vertical profiles of spatially averaged velocity $\left(\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}\right)$, TKE $\left(\langle k\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$, and turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$ obtained from LES simulations (PALM model described in Sect. 2.2) compared with the RANS simulation (Santiago and Martilli, 2010). “A” and “S” correspond to aligned and staggered arrays of buildings, and different urban densities (λ_p) are considered.

Figure 6 Vertical profiles of normalized velocity $\left(\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}\right)$, turbulent kinetic energy $\left(\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$, and turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$ obtained from LES results, spatially and time- averaged for different λ_p and configurations. $\langle \stackrel{\mathrm{‾}}{k}\rangle$ from the LES results is calculated based on the turbulent variances at the resolved scale as well as modeled subgrid-scale TKE. “A” in this graph indicates “aligned” configuration, while “S” stands for “staggered”.

3 Results and discussions

3.1 Large-eddy simulations: vertical profile of mean flow and dispersive stress

Figure 6 displays the vertical profiles of flow parameters ($\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}$, $\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$, and $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$) spatially averaged within the roughness sublayer for two urban configurations (aligned and staggered) and varying urban density (λ_p). It is evident that the flow profiles are significantly influenced by the urban configuration. Overall, average wind speed, and consequently the turbulent momentum flux and TKE, are significantly larger for aligned arrays of cubes with streamwise flow aligned with urban street canyon. However, it is worth noting that in real cities, the aligned configuration with a 0^∘ wind angle may not be the most representative of the flow field. Real cities experience a range of wind directions relative to the street grid, and many cities do not have a grid but rather streets of many orientations. Our simulations (similar to many urban CFD simulations) represent buildings with two street directions oriented perpendicular to each other, with streamwise flow oriented perpendicular to one set of building faces. The aligned version of this setup represents a special case relative to real cities: those scenarios where wind direction is aligned with one of the two street directions. The staggered version of this setup, conversely, presents no major corridors (i.e., streets) aligned with the wind that do not include building drag. As such, we believe that the staggered configuration better represents the impacts of real cities on urban canopy flow under a variety of wind directions. Any choice here is a simplification of reality, and the choice of a regular staggered array provides a closer approximation to average conditions in real cities in our estimation.

Another investigation made here is regarding the significance of dispersive fluxes in urban canopy parameterizations. In the formulation of the multi-layer urban canopy model, the dispersive transport processes are neglected so far (Santiago and Martilli, 2010; Krayenhoff et al., 2015), while in fact they are non-negligible in many real urban configurations (Giometto et al., 2016). The variability of the spatially averaged dispersive stress obtained from LES for varied urban configuration and packing density and the contribution of $\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle$ to the total turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle +\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle \right)$ is represented in Fig. 7. It is observed that the dispersive stress can in fact be in the same order of the total momentum flux. Hence, given the importance of the dispersive term in the momentum budget, the subsequent analysis seeks to represent the effects of dispersive motions in the column model by driving the parameterization from the 3-D results of $\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle$ together with $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle$ (Sect. 3.3). Note that positive values of dispersive fluxes within the canopy, of similar magnitude to the turbulent stress, implies that the flux is countergradient, indicating downward transport of slow air.

3.2 Drag parameterization

It is known that the sectional drag coefficient depends on the packing density and the configuration of the array with a strong dependency on height, such that $C_{\text{d}}=C_{\text{d}}(z,{\mathit{\lambda }}_{\text{p}})$ (Macdonald, 2000; Santiago et al., 2008; Santiago and Martilli, 2010). However, as indicated by Santiago and Martilli (2010), height-dependent parameterization of drag coefficients is challenging due to the high variability of C_d close to the ground due to small $\langle \stackrel{\mathrm{‾}}{u}\rangle$ as well as the lack of experimental information on the vertical profiles of this property inside the urban canopy. Accordingly, Santiago and Martilli (2010) proposed the following calculation of equivalent drag coefficient that is constant with height in the urban canyon, considering that “when integrated in the whole urban canopy, the drag force must be equal to that computed by the CFD simulations”.

$$\begin{array}{}\text{(10)}& C_{\text{deq}}={\displaystyle \frac{\frac{-\mathrm{1}}{\mathit{\rho }H}{\int }_{\mathrm{0}}^H\mathrm{\Delta }\langle \stackrel{\mathrm{‾}}{p\left(z\right)}\rangle \text{d}z}{\frac{\mathrm{1}}{H}{\int }_{\mathrm{0}}^H\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle |\langle \stackrel{\mathrm{‾}}{u\left(z\right)}\rangle |\text{d}z}}\end{array}$$

Following this method, the drag coefficient parameterization using the LES results is shown in Fig. 8. C_deq is computed by means of the ratio between the horizontally averaged mean pressure deficit around an obstacle and the square of the horizontally averaged mean velocity around the obstacle (Eq. 10). C_deq depends on the configuration (aligned or staggered) and packing density (λ_p) of the array shown here. Comparing the LES and RANS results, the trends in C_deq with λ_p are in good agreement, but as previously demonstrated by Simón-Moral et al. (2014), RANS tends to overestimate the value of C_deq.

Figure 7 (a) Vertical (spatially and temporally averaged) profiles of normalized dispersive stress $\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$ and (b) the contribution of dispersive stress to the total turbulent momentum flux $\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle /\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle +\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle \right)$. “A” in this graph indicates “aligned” configuration, while “S” stands for “staggered”.

Figure 8 Variation of C_deq with urban packing density λ_p for the staggered configuration and compared with the RANS results of Santiago and Martilli (2010). The fitted line indicates the parameterization proposed for the 1-D model. Additional data points (in green) are added for parameterization corresponding to λ_p=0 (zero building-induced drag) and λ_p>0.44 (where the drag coefficient reaches a constant value for high urban packing densities).

3.3 Length scale parameterization

In this section, the length scales obtained from the spatially averaged LES results and the k−l turbulence closure theory for the urban canopy parameterization are discussed. Combining Eqs. (2) and (3), the turbulent length scale C_kl_k is traditionally calculated only considering the Reynolds stress $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle$ (Eq. 11a). Here, following the discussions in Sect. 3.1, we recalculate turbulent length scale using total momentum fluxes that include turbulent dispersive flux $\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle$, shown as C_kl_kM in Eq. (11b). Note that in the column model, C_kl_k (l_ε∕C_ε) is parameterized instead of l_k (l_ε) in order to avoid the uncertainties regarding the values of C_k (C_ε) proposed in the literature.

$$\begin{array}{}\text{(11a)}& {\displaystyle }& {\displaystyle }\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle =-C_k{l}_k\langle \stackrel{\mathrm{‾}}{k}{\rangle }^{\mathrm{1}/\mathrm{2}}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial z}}\text{(11b)}& {\displaystyle }& {\displaystyle }\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle +\langle \stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\rangle =-C_k{l}_{kM}\langle \stackrel{\mathrm{‾}}{k}{\rangle }^{\mathrm{1}/\mathrm{2}}{\displaystyle \frac{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle }{\partial z}}\end{array}$$

Figure 9 shows the vertical profile of turbulent length scale for varied urban densities (panel a) and the canopy-averaged length scale obtained from RANS and LES with or without considering the dispersive term (panel b). Two observations are made. First, we observe that the turbulence length scale is larger for the LES results within the building canopy, specifically when dispersive stress is included. This is further explained by the significant difference in the TKE profile between LES and RANS shown in Fig. 5. Second, the length scale calculated using LES does not vary monotonically with urban packing density (λ_p) but rather follows the behavior of roughness length (Grimmond and Oke, 1999). This can be explained due to the varying flow regimes from the isolated (λ_p=0.0625) to wake interference (λ_p=0.25) and skimming (λ_p=0.44) flow. As noted by Grimmond and Oke (1999), “as the density increases so does the roughness of the system, but a point comes where adding new elements merely serves to reduce the effective drag of those already present due to mutual sheltering. This reduces the effective height of the canopy for momentum exchange.” Accordingly, we observe that the drag coefficient (Fig. 8) plateaus with increasing density. Similarly, the non-monotonic behavior in LES-derived turbulent length scales (that resolve the turbulent flux of momentum and energy across larger scales of motions as opposed to RANS) can be attributed to different flow regimes. The LES results suggest that the largest scales of turbulence (i.e., the most turbulent organization) are produced in the wake interference regime. As the turbulent length scale is a measure of the efficiency of vertical transport, the higher l_k values indicate higher vertical transport of momentum (including turbulent and dispersive) for the same TKE and vertical flow gradient. For intermediate λ_p, mainly in the wake interference regime, the presence of the buildings favor the formation of organized motions, likely at the scale of the buildings, that enhance the vertical transport. For higher λ_p values, in the skimming flow, these movements are suppressed, and for isolated buildings they are less strong. Given the time-averaged representation of the turbulent field, RANS is not able to reproduce these effects, resulting in the monotonic decrease in derived turbulent length scale with urban density.

Figure 9 (a) Vertical profiles of turbulent length scale calculated using LES simulations with the dispersive stress included for aligned (A) and staggered (S) urban configurations. (b) Variation in normalized turbulent length scale averaged within the building canyon ($z/H<\mathrm{1}$) with plan area density (λ_p) for staggered configurations. Data points represent (1) the RANS simulation (Santiago and Martilli, 2010), (2) the LES simulation without dispersive stress included, and (3) the LES simulation with dispersive stress included.

To assess the dissipation length scale l_ε∕C_ε, Eq. (6) is used assuming that the dissipation is only happening at the subgrid scale, and therefore is correlated with the SGS-TKE (Maronga et al., 2015). Figure 10a demonstrates the vertical profile of l_ε∕C_ε for staggered urban configurations with variable packing density. The vertical profiles of l_ε∕C_ε exhibit similar characteristics compared to the RANS results of Santiago and Martilli (2010): inside the canopy, the length scale is mostly constant with height (specifically for λ_p≥0.25) as it is controlled by the shear layer (H−d, where d is the zero-plane displacement height), above the canopy the dissipation length scale increases with height, and the lower values of l_ε∕C_ε close to the ground (particularly for lower λ_p) and at building height correspond to the locations of maximum dissipation in the urban canopy. This is likely due to the fact that dissipation depends only on small-scale motions and therefore is less affected by larger structures induced by the presence of the buildings.

Three different zones are then defined consistent with Santiago and Martilli (2010) and Krayenhoff et al. (2015) to parameterize l_ε∕C_ε: (a) inside the canopy ($z/H<\mathrm{1}$), l_ε∕C_ε is assumed to be constant with height; (b) well above the canopy ($z/H>\mathrm{1.5}$), where the behavior of l_ε∕C_ε is linear and the slope varies with λ_p; and (c) in the zone of transition ($\mathrm{1}\le z/H\le \mathrm{1.5}$) between the two previous zones:

Figure 10 (a) Vertical profiles of dissipation length scale (l_ε∕C_ε) for staggered (S) urban configuration and varying urban packing density (λ_p). l_ε∕C_ε is calculated based on dissipation at the subgrid scale (Maronga et al., 2015) and Eq. (6). (b) Vertical profiles of model constant for turbulent viscosity (C_μ) calculated based on Eq. (7).

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle }{l}_{\mathit{\epsilon }}/C_{\mathit{\epsilon }}={\mathit{\alpha }}_{\mathrm{1}}(H-d)z/H<\mathrm{1},\text{(12b)}& {\displaystyle }& {\displaystyle }{l}_{\mathit{\epsilon }}/C_{\mathit{\epsilon }}={\mathit{\alpha }}_{\mathrm{1}}(z-d)\mathrm{1}\le z/H\le \mathrm{1.5},\text{(12c)}& {\displaystyle }& {\displaystyle }{l}_{\mathit{\epsilon }}/C_{\mathit{\epsilon }}={\mathit{\alpha }}_{\mathrm{2}}(z-d_{\mathrm{2}})z/H>\mathrm{1.5},\end{array}$$

where α₁=4 is the revised value computed for all λ_p cases in LES, the displacement height (d) parameterization is taken from Krayenhoff et al. (2015) and Simón-Moral et al. (2014) as

$$\begin{array}{}\text{(13)}& {\displaystyle }d\left({\mathit{\lambda }}_{\text{p}}\right)=H{\mathit{\lambda }}_{\text{p}}^{\mathrm{0.15}},\end{array}$$

and finally α₂ and d₂ are parameterized as

$$\begin{array}{}\text{(14a)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{2}}\left({\mathit{\lambda }}_{\text{p}}\right)& {\displaystyle }=min\left(\mathrm{5},max\left(\mathrm{2},\mathrm{1.3}{\mathit{\lambda }}_{\text{p}}^{-\mathrm{0.45}}\right)\right),\text{(14b)}& {\displaystyle }{d}_{\mathrm{2}}\left({\mathit{\lambda }}_{\text{p}}\right)& {\displaystyle }=\mathrm{1.5}H\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{2}}}}\right)+d\left({\mathit{\lambda }}_{\text{p}}\right){\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{2}}}}.\end{array}$$

Note that α₁ (in-canopy) does not vary significantly with urban density, while α₂ (slope of l_ε∕C_ε above $Z/H=\mathrm{1.5}$) is a function of λ_p. Additionally, the parameterization of the dispersive length scale below building height (α₁) is slightly underestimated compared to the LES results to account for the localized maximum of dissipation close to the ground specifically for low urban densities (Fig. 11).

Comparing the turbulent (C_kL_kM) and dissipation (l_ε∕C_ε) length scales (Eq. 7), however, we find that the assumption of constant C_μ in the canopy does not hold in the LES results. Figure 10b demonstrates the vertical profile of C_μ calculated based on Eq. (7), and we observe the variability in in-canopy C_μ with λ_p. Accordingly, in addition to the dissipation-length-scale parameterization provided for the 1-D model, the C_μ value in the canopy is also parameterized based on λ_p (Fig. 11b), while above the canopy, C_μ=0.05 for all urban densities:

$$\begin{array}{}\text{(15)}& {\displaystyle }{C}_{\mathit{\mu }}=\left\{\begin{array}{ll}max(\mathrm{0.05},-\mathrm{1.6}{\mathit{\lambda }}_{\text{p}}^{\mathrm{2}}+\mathrm{0.75}{\mathit{\lambda }}_{\text{p}}+\mathrm{0.022})& z/H\le \mathrm{1}\\ \mathrm{0.05}& z/H>\mathrm{1}\end{array}\right..\end{array}$$

Figure 11 Variation in normalized dissipation length scale (l_ε∕C_εH) and model constant for turbulent viscosity (C_μ) averaged within the canopy ($z/H<=\mathrm{1}$) with plan area density (λ_p). The results are obtained using the staggered urban configuration (Fig. 3) and averaged in the canopy volume ($\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{Q}}$ indicates volume average of quantity Q).

3.4 Assessment of one-dimensional column model with LES and RANS results

The drag coefficient and length scales (Sect. 3.2 and 3.3) parameterizations derived from LES results are used to update the multi-layer (1-D) urban canopy model and evaluated here against (1) the RANS-derived multi-layer (1-D) model, (2) 3-D LES results with idealized configuration (present study), and (3) LES results with realistic urban configurations (Giometto et al., 2017).

Figure 12 shows the vertical profiles of horizontally averaged velocity, turbulent kinetic energy, and turbulent momentum flux calculated with the 1-D (multi-layer) model and compared with the LES results. The 1-D model results are calculated using previous parameterizations with RANS (Santiago and Martilli, 2010) as well as the updated LES parameterizations for λ_p=0.0625, 0.25, and 0.44. The vertical profile of $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$ obtained with the updated (LES) 1-D model shows improvement for all studied λ_p cases. For horizontally averaged $\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}$ and $\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$, however, the performance of the model is dependent on the λ value. Overall, the prediction of the horizontally averaged velocity is improved compared to 1D-RANS, particularly within the canopy ($z/H<\mathrm{1}$). However, a significant underestimation of wind speed is seen at the higher urban density. TKE profiles, on the other hand, are overestimated for λ_p=0.0625 while significantly improved for other cases. Additionally, despite the improvements with the new parameterization, the TKE close to the ground is still substantially underestimated for high-λ_p cases, indicating that there is underestimation of vertical turbulent transport deep in the canopy. This could be traced back to the parameterization of the TKE transport in the multi-layer model that assumes the same diffusion coefficient (K_m) for momentum and TKE equations, which does not hold in the LES results (not shown).

Figure 12 Comparison of the vertical profiles of velocity $\left(\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}\right)$, turbulent kinetic energy $\left(\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$, and turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$ obtained with the multi-layer (1-D) model with RANS and LES parameterization with the LES results for various λ_p.

Figure 13 further demonstrates the root mean square error (RMSE) of horizontally averaged velocity $\left(\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}\right)$, turbulent kinetic energy $\left(\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$, and turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$ compared with the LES results. RMSE is calculated for $z=\mathrm{0}-\mathrm{3}H$ for all λ_p cases studied here. It can be seen that the new parameterizations with LES (depicted in dark blue) represent an overall improvement compared to the previous multi-layer model derived from the RANS results. The RMSEs for $\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$ and $\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}$ are substantially lower in the updated multi-layer model with LES-derived parameters and formulations. For high-λ_p cases, however, the new parameterization underperforms in predicting the wind flow.

Figure 13 Root mean square error calculated for vertical profiles of velocity ($\frac{u}{u_{\mathit{\tau }}}$), TKE ($\stackrel{\mathrm{‾}}{k}/u_{\mathit{\tau }}^{\mathrm{2}}$), and Reynolds stress ($u^{\prime }{w}^{\prime })/u_{\mathit{\tau }}^{\mathrm{2}}$. The RMSE values are calculated using 1-D models with LES and RANS parameterization up to 3H.

Lastly, the 1-D multi-layer model is compared with the LES results of Giometto et al. (2017), which are conducted for a realistic urban neighborhood (Vancouver Sunset) in BC, Canada. The neighborhood characteristic in the modeled urban canopy subset (indicated as S1 in Giometto et al., 2017) is λ_p=0.34, and average building height is 6.6 m. The studied case in Fig. 14 represents a configuration without trees (given the fact that tree parameterization was not the focus of the current study). We observe that the updated parameterizations in the 1-D multi-layer model result in a substantial improvement compared to Santiago and Martilli (2010), specifically for the vertical profile of turbulent kinetic energy as well as wind speed above the building height. However, underestimation of wind speed and Reynolds stress in the street canopy is observed, which is likely attributed to the building configuration and wind direction considered in the realistic LES simulations. In Giometto et al. (2017), the urban configuration resembles evenly spaced aligned buildings with wind direction aligned with one of the primary street directions. This results in relatively linear profiles of wind speed and Reynolds stress in the canopy, which as discussed before, only represent one realization of urban canopy flow. Nonetheless, this demonstrates the need for assessing urban canopy parameterizations with various urban configurations and wind directions in the future. Additionally, underestimation of TKE deep in the canopy is seen again, further indicating that the current parameterization of the turbulent transport in the urban canopy is not adequate to determine $\langle \stackrel{\mathrm{‾}}{k}\rangle$ at the ground level, particularly in higher-λ_p cases.

Figure 14 Comparison of vertical profile of velocity $\left(\langle \stackrel{\mathrm{‾}}{u}\rangle /u_{\mathit{\tau }}\right)$, turbulent kinetic energy $\left(\langle \stackrel{\mathrm{‾}}{k}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$, and turbulent momentum flux $\left(\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}\rangle /u_{\mathit{\tau }}^{\mathrm{2}}\right)$ obtained with multi-layer (1-D) model with the (3-D) LES results of Giometto et al. (2017) for realistic urban configurations (G2017) as well as LES simulations discussed here for idealized configurations.

4 Summary and conclusions

The present study focused on updating the urban canopy parameterizations of drag coefficient and turbulent length scales using large-eddy simulations (LES) results, which is shown to be a superior numerical model for resolving the turbulent flow field compared to Reynolds-averaged Navier Stokes (RANS) previously used in multi-layer UCMs (Santiago and Martilli, 2010).

The detailed analyses of the spatially averaged turbulent field in urban configurations revealed the following: (1) LES results exhibit a significantly higher transport of TKE into the lower canopy compared to RANS; (2) dispersive fluxes are not negligible in the urban canopy, particularly in higher urban packing densities; and (3) the ratio between turbulent and dispersive length scale (commonly described by the model constant C_μ in multi-layer models) is not constant with λ_p at the canopy level. These findings motivated the revision of the UCPs to include dispersive fluxes and further parameterize turbulent length scale (C_kl_k) in addition to dissipation length scale (l_ε∕C_ε) through the parameterization of model constant C_μ.

We demonstrated that using LES results as the basis for parameterization, as well as the inclusion of dispersive stress, improves the performance of the multi-layer model, such that spatially averaged profiles of flow, and consequently the turbulent exchange in the urban canopy in realistic neighborhoods, can be predicted more accurately. Additionally, when the updated parameterizations were used in the BEP-Tree model (Krayenhoff, 2014), we observed improved performance compared to measurements taken across the diurnal cycle at three sites located in Vancouver (BC) and London (ON) in Canada and Salt Lake City (UT) in USA (Krayenhoff et al., 2020).

However, spatially averaged turbulent kinetic energy, $\langle \stackrel{\mathrm{‾}}{k}\rangle$, is still underestimated close to the ground for high λ_p values due to the underestimation of turbulent transport deep in the canopy. Preliminary analyses of turbulent transport in this study (not shown) reveal that the K-theory assumption that the diffusion coefficient K_m is the same for TKE and momentum equations (i.e., $K_{\text{m}}=-\frac{\langle \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }\rangle }}{\partial \langle \stackrel{\mathrm{‾}}{u}\rangle /\partial z}=-\frac{\langle \stackrel{\mathrm{‾}}{k^{\prime }{w}^{\prime }\rangle }}{\partial \langle \stackrel{\mathrm{‾}}{k}\rangle /\partial z}$) does not hold in the LES results. Accordingly, future work should revisit the multi-layer model formulations to assess (1) the parameterization of turbulent transport term in the 1-D TKE equation (Eq. 5) and (2) the distinction between the diffusion length scale of momentum and TKE.

Further analysis is also needed to fully evaluate the effects of idealized configurations in parameterizations and assess the impact of variable building heights and wind directions on turbulent length scales and drag parameterization. Santiago et al. (2013a) showed that a height-dependent drag coefficient is needed to capture the lateral effects within the canopy for oblique wind directions. To further account for the street and wind directions in realistic configurations, future work is needed to develop a methodology that derives dominant street directions over each grid cell and computes the drag coefficient as a function of height and the angle between street and wind direction above the canopy. Lastly, the current study focused on the momentum exchange without considering the role of thermal forcing on turbulent length scales. An updated parameterization of thermal effects (investigated by Krayenhoff et al., 2020) can also be evaluated using LES results.