Articles | Volume 17, issue 7
Model description paper
05 Apr 2024
Model description paper |  | 05 Apr 2024

A one-dimensional urban flow model with an eddy-diffusivity mass-flux (EDMF) scheme and refined turbulent transport (MLUCM v3.0)

Jiachen Lu, Negin Nazarian, Melissa Anne Hart, E. Scott Krayenhoff, and Alberto Martilli

In recent years, urban canopy models (UCMs) have been used as fully coupled components of mesoscale atmospheric models as well as offline tools to estimate temperature and surface fluxes using atmospheric forcings. Examples include multi-layer urban canopy models (MLUCMs), where the vertical variability of turbulent fluxes is calculated by solving prognostic momentum and turbulent kinetic energy (TKE, k) using mixing length scale (l) and drag parameterizations. These parameterizations are based on the well-established 1.5-order kl turbulence closure theory and are often informed by microscale fluid dynamics simulations. However, this approach can include simplifications such as assuming the same diffusion coefficient for momentum, TKE, and scalars. In addition, the dispersive stresses arising from spatially averaged flow properties have been parameterized together with the turbulent fluxes despite being controlled by different mechanisms. Both of these assumptions impact the quantification of the turbulent exchange of flow properties and subsequent air temperature predictions in urban canopies. To assess these assumptions and improve corresponding parameterization, we analyzed 49 large-eddy simulations (LES) for idealized urban arrays, encompassing variable building height distributions and a comprehensive range of urban densities (λp[0.0625,0.64]) seen in global cities. We find that the efficiency of turbulent transport (numerically described via diffusion coefficients) is similar for scalars and momentum but is 3.5 times higher for TKE. Additionally, parameterizing the dispersive momentum flux using the kl closure was a source of error, while scaling with the pressure gradient and urban morphological parameters appears more appropriate. In response to these findings, we propose two changes to the previous version of MLUCM: (a) separate characterization for turbulent diffusion coefficient for momentum and TKE and (b) introduction of an explicit physics-based “mass-flux” term to represent the fraction of the dispersive momentum transport directly induced from buildings as an amendment to the existing “eddy-diffusivity” framework. The updated one-dimensional model, after being tuned for building height variability, is further compared against the original LES results and demonstrates improved performance in predicting vertical turbulent exchange in urban canopies.

1 Introduction

The urban canopy is a unique and complex land cover type in climate models (Oke et al.2017). As climate models become more capable of high-resolution simulations, the challenge resides in developing and validating urban canopy models (UCMs) to account for finer-scale interactions with complex urban forms. Applications such as running UCMs to accurately describe urban surface fluxes to atmospheric boundary layers (Schoetter et al.2020) and UCMs as offline models coupling with mesoscale forcing to simulate in-canopy processes (Redon et al.2020) can then benefit from such infrastructure developments.

The development of UCMs relies on the similarity between flow through the urban canopy layer (UCL) and other more commonly investigated types of flow. For example, Macdonald (2000) modified the exponential wind profile from flow over vegetative canopy (Finnigan2000) based on the aerodynamic roughness properties of the urban surface to predict the urban wind speed that is commonly used in single-layer urban canopy models (SLUCM (Kusaka et al.2001)). The 1.5-order turbulence closure scheme (Bougeault and Lacarrere1989) used for turbulence parameterization in the planetary boundary layer (PBL) was also extended to the UCL, facilitating the development of the multi-layer urban canopy model (MLUCM (Nazarian et al.2020)). The urban form and fabric, however, feature distinct behaviors that require modifications to length scales (Cheng and Porté-Agel2021) and form drag (Sützl et al.2021b) that are uniquely developed to relate urban morphological parameters (Lu et al.2023b). In addition, some physical processes commonly seen in urban environments, such as thermal instability, also result in unique flow characteristics in urban canopies that have been separately evaluated in urban canopy models (Santiago et al.2014; Simón-moral et al.2016; Nazarian and Kleissl2015).

The modeling of urban surfaces develops along with the capability of computational power and availability of realistic building morphology (Sirko et al.2021). As the community is pushing the resolution of regional climate models to the sub-kilometer scale (Sützl et al.2021a), there is an urgent need for a more detailed characterization of finer(subgrid)-scale (and higher-order) urban phenomena (such as building wake and canyon re-circulation (Oke et al.2017)). The most common way to enable UCMs to capture these factors is by relating the variation of urban geometric parameters to the corresponding flow statistics based on computational fluid dynamics (CFD) models, where the large-eddy simulation (LES) offers a feasible way to obtain the higher moment of flow statistics.

The development of UCMs includes designing modules to predict flow properties such as wind speed (Castro2017), turbulent kinetic energy (TKE) (Christen et al.2009), and scalars (Lim et al.2022) that interact with urban structures differently. For example, increasing the frontal density of roughness elements was found to amplify the turbulent momentum transport but hinder the turbulent scalar transport (Li and Bou-Zeid2019). However, efforts to identify and characterize the dissimilarity of the transport mechanisms of flow properties are rare due to the lack of experimental evidence to support quantitative characterizations. Therefore, momentum, TKE, and scalar diffusion efficiency in modeling practice are typically assumed to be the same (Martilli et al.2015; Nazarian et al.2020), introducing errors that may degrade model performance.

Multi-layer UCMs are fundamentally more versatile than single-layer UCMs in characterizing the urban effects (Hendricks et al.2020; Lipson et al.2023) and in-canopy processes, and the use of UCMs has evolved into different implementations (Kondo et al.2005; Yuan et al.2019; Santiago and Martilli2010b). The building effects parameterization (BEP (Martilli et al.2002)) model and developments to its flow parameterization through the multi-layer urban canopy models (MLUCMs (Santiago and Martilli2010a)) have been directed toward multiple applications such as simulation of building–tree interactions (BEP-tree (Krayenhoff et al.2020)) and indoor building energy exchanges with outdoor climate (BEP-BEM (Salamanca et al.2010)). Simulations with MLUCMs are particularly sensitive to details of the parameterization that influence the strength of turbulence intensity within the canopy that facilitates turbulent transport of flow properties. Thus, accurate quantification of TKE is essential; however, previous work leading to the development of MLUCM v2.0 has shown an underestimation of canopy TKE (Nazarian et al.2020). Moreover, canopy TKE must be accurately quantified for simplified neighborhoods before the addition of the flow effects arising from the added complexity of natural urban areas, such as sub-facet-scale structures (overhangs, awnings, irregular building shapes, etc.) or trees (Krayenhoff et al.2020).

Arising from spatially averaged flow properties (that is core to MLUCM development), dispersive fluxes illustrate the transport of variables by time-mean structures smaller than the averaging grid size, constituting another unique urban phenomenon (Poggi and Katul2008b). The strength of the dispersive flux is highly related to the horizontal (e.g., Lu et al.2023b) and vertical (e.g., Xie et al.2008) urban structures and exhibits substantial spatial variability (Harman et al.2016) that precludes generalized characterizations. Therefore, dispersive flux has only been pragmatically considered in Simón-moral et al. (2016) based on the K-theory framework as an increment to its turbulent counterpart. In practice, this approach resulted in an enhanced eddy diffusivity and implicitly assumed the correlation between dispersive flux and turbulence intensity and the vertical gradient of the flow property. However, the physical basis of such a correlation is still unclear and may not comprehensively improve the predictability of UCMs. Varying with the shape of roughness elements, the dispersive component can be about 50 % of the total transport over non-cubic blocks (Li and Bou-Zeid2019). Furthermore, realistic urban neighborhoods induce even higher dispersive contributions due to the more substantial spatial flow variability (Akinlabi et al.2022; Giometto et al.2016). Consequently, the lack of dispersive stress characterization restricts the study of canopy flows to relatively less complex and homogeneous flow, such as on gentle topography (e.g., Finnigan and Belcher2004) with slowly varying packing density (Ross2012).

Based on the aforementioned fundamental deficiencies embedded in the multi-layer model, a better performance of UCMs on more complex flows requires a finer and more physically based characterization of dispersive fluxes. In this research, based on a recently developed urban flow dataset (Lu et al.2023a), we aim to assess and improve the prediction of the TKE and momentum in the urban canopy through a more accurate parameterization of turbulent transport in the MLUCM (Nazarian et al.2020) in two ways:

  • separate characterization of transport efficiency of different flow properties.

  • introduce a physics-based modeling of the dispersive transport of momentum.

This paper first describes the flow dataset using the Parallelized Large-Eddy Simulation Model (PALM) in Sect. 2.1. Then different components of MLUCM v2.0 (Nazarian et al.2020) are revisited for the conventional eddy diffusion module in Sect. 3.1 and for the drag and dissipation module in Sect. 3.3. A novel physically based mass-flux parameterization is proposed based on the pressure gradient scaling and the geometric condition of the neighborhood in Sect. 3.2. In Sect. 4, the proposed changes for MLUCM are tested and compared across 49 urban arrays with different horizontal (staggered, aligned) and vertical (height variabilities) arrangements. Section 5 summarizes the findings of this study and provides perspectives for future developments of the multi-layer model.

2 Methodology

In this study, we used two methods to assess the turbulent transport behavior of the urban canopy flow, i.e., computational fluid dynamics (CFD) simulations and an offline multi-layer urban canopy model. In Sect. 2.1, with the help of LES simulations recording the third-order moments in turbulent flows, we evaluate differences in turbulent transport of momentum, TKE, and scalars. Following the differences, we further revisit the turbulent parameterization in Sects. 3.1 and 3.2 followed by a comprehensive evaluation of model performance against LES in Sect. 4.

2.1 Large-eddy simulation (LES) setup and dataset

Conventionally, two configurations of idealized building arrays, “aligned” and “staggered” (Coceal et al.2006), are employed to approximate horizontal arrangements of real neighborhoods. The aligned array simulates the situation in natural communities where streets align with the prevailing wind direction. The staggered array, by contrast, is mainly employed to calibrate UCM (Santiago and Martilli2010b; McNorton et al.2021) in that it has no significant corridors (i.e., streets) aligned with the wind, potentially resembling the frequent situation where wind and street directions do not align. Both choices simplify realistic urban neighborhoods, but staggered arrays provide a closer approximation to average conditions in real cities (Nazarian et al.2020). Accordingly, the staggered configuration was selected to obtain urban wind flow in determining parameterizations in MLUCM.

Simulations are conducted in the Parallelized Large-Eddy Simulation Model (PALM, version r4554) (Maronga et al.2020) with λp ranging from 0.0625 to 0.64, representing urban densities seen in global cities. The numerical setup follows (Lu et al.2023a; Nazarian et al.2020; Lu et al.2023b), which has been validated against direct numerical simulation (DNS) (Coceal et al.2007) and wind tunnel experiments (Brown et al.2001) but will be listed here for completeness.

The computational domain is discretized equidistantly using second-order central differences (Piacsek and Williams1970) with a horizontal and staggered Arakawa C-grid in the vertical direction. The momentum field is solved using the filtered Navier–Stokes equation with the Boussinesq approximation with the time integration following a minimal storage scheme (Williamson1980). The covariance terms from the filtering procedure were parameterized using a 1.5-order closure scheme after (Deardorff1980). The pressure perturbation was calculated in Poisson's equation and was solved by the FFTW scheme (Frigo and Johnson1998). Apart from solving momentum and pressure fields, a passive scalar was introduced with an initial surface value c0=0.06 and a constant negative gradient of -5×10-5 m−1 to form the initial profile. At the same time, the total scalar concentration was maintained by two Neumann boundary conditions with surface emission equal to 0.0001 and a sink at the domain top equal to 0.0001×(1-λp,0), where λp,0 indicates urban density at the surface.

The detailed design of building arrays can be found in Lu et al. (2023a) and is briefly explained here: while maintaining the same volume of roughness elements (i.e., same mean building height, Hmean=16 m), we consider three height distributions, i.e., Hstd[0m,2.8m,5.6m], which result in different maximum building height Hmax. Additionally, we consider three types of urban arrays: (a) continuous distribution of building height with a relatively low Hmax; (b) clustered configuration with two clusters (three buildings) of low and tall building blocks; and (c) high-rise configuration, featuring one tall and a cluster (three) of low building blocks. Each configuration is then run with seven urban packing densities, i.e., λp[0.0625,0.1111,0.16,0.25,0.35,0.4444,0.64], which yields a total of 49 simulation cases. Details of the configurations considered are shown in Table 1.

Table 1Dataset details for 49 urban arrays discussed in this study. Four building height arrangements (uniform, continuous, clustered, and high-rise) with two standard deviations of height configurations are considered. The maximum Hmax and minimum Hmin building heights and a three-dimensional view of an example (λp=0.25) are shown below each category.

Download Print Version

3 Assessment of the multi-layer urban canopy model (MLUCM)

In mesoscale climate models, flow through canopies is often simulated at scales much larger than the typical surface processes, such as turbulent eddies and obstacle wakes. A common approach to derive parameterizations is to apply a time averaging over an interval longer than the scale of the slowest eddies and a spatial averaging over a length larger than the typical spatial deviations in the flow (Raupach and Shaw1982; Finnigan1985). For momentum, the Reynolds decomposition is first applied to the three-dimensional instantaneous equations that decompose mean flow quantities (ϕ) from their fluctuating components (ϕ, time or ensemble averaging, ϕ=ϕ+ϕ). Then, flow properties are spatially averaged to match a grid cell of a mesoscale model (horizontal averaging, ϕ=ϕ+ϕ̃, where ϕ̃ is termed a dispersive component).

In urban canopies where the volume fraction occupied by obstacles is not negligible, the intrinsic average (Schmid et al.2019) that excludes the roughness fraction in the spatial averaging is more favorable given that additional terms compensating for the volume fraction occupied by roughness elements in the canopy are included Blunn et al. (2022). The friction Reynolds number (Re=uτHT/ν, where uτ∼0.21 m s−1 is the friction velocity, HT is the domain height, and ν is the kinematic viscosity) of the flow dataset is large enough (106) to neglect viscous stress on momentum and molecular flux on scalars. The resulting one-dimensional representation of the momentum and passive scalar equation assuming horizontal homogeneity, negligible Coriolis effect, and buoyancy force reads


where ρ is the air density, u and w are the streamwise and vertical velocity so that uw and ũw̃ form the spatially averaged turbulent and dispersive Reynolds stress, respectively. The first term on the RHS of both equations represents vertical transport events, while the second term is a term that compensates for volume change in the intrinsic averaging, where γ(z)=1-λp(z) is the fluid fraction at the vertical index z. The third term, fi and ci, accounts for external sources or sinks, such as pressure gradient driving the flow and surface emission of pollutants. The fourth term of Eq. (1) represents the pressure drop between windward and leeward building facets (Santiago and Martilli2010b), which accounts for momentum sink due to building-induced drag rising from spatial averaging over the momentum equation. The multi-layer model is designed to relate terms in Eq. (1) to the variation of urban surface, where the vertical diffusion of momentum (Sects. 3.1 and 3.2) and drag (Sect. 3.3) are parameterized.

3.1 Characterization of the eddy diffusivity for different flow statistics

The Reynolds averaging over the momentum equation yields nonlinear terms that must be parameterized to close the equation. The most common approach in parameterizing the Reynolds stress in UCMs is based on K theory (as also seen in models compared in Best and Grimmond2015, Lipson et al.2023) to represent the efficiency of turbulent transports by employing the eddy-diffusivity coefficient together with the vertical gradient of the flow properties (for momentum, Km=-uwuz). In a variety of ways specifying the eddy diffusivity, the 1.5-order closure (Martilli et al.2002) assuming the transport strength related to the turbulence intensity is commonly employed for urban canopy flow as follows:

(3) K m = C k l k k 1 / 2 ,

where Ck is a model constant, lk is a turbulent length scale, and k is TKE that has to be prognostically obtained during the simulation following Eq. (4).

(4) k ¯ t = - k w + k ̃ w ̃ z - k w + k ̃ w ̃ γ γ z - u i u j u i x j - u i u j ̃ u i ̃ x j - D k - ε ¯

Similar to the momentum equation (Eq. (1)), the first two terms on the RHS represent the transport of TKE, followed by two terms for shear and dispersive TKE production, respectively. Dk represents the source of TKE generated through the interaction with the buildings and the airflow, while the last term ε¯ is the viscous dissipation of TKE that can be parameterized as

(5) ε ¯ = C ε k ¯ 3 / 2 l ε ,

where Cε and lε are the model constant and the dissipation length scale, respectively, and have a relationship with turbulent length scale based on the turbulent viscosity Cμ (Santiago and Martilli2010b) that measures the relative strength of viscous dissipation over turbulent transport:

(6) C k l k = C μ l ε C ε .

To date, most models assume the same value for the diffusion coefficient in the momentum, TKE, and scalar equations, possibly due to the lack of data on TKE and scalar budget terms (note that the explicit evaluation of TKE transport and diffusivity, in particular, is only possible from numerical simulations that resolve third-order moments, e.g., LES and direct numerical simulations (DNS)). However, previous LES analyses of urban flow indicated a different transport efficiency of varying flow properties (Nazarian et al.2020). Accordingly, we distinguish between three diffusion coefficients as


The turbulent fluxes, vertical gradient, and eddy diffusivity of momentum, TKE, and scalars are shown in Fig. 1 for staggered-uniform configurations for cubic buildings. The roughness flow over other building shapes such as cuboids (Blunn et al.2022; Li and Bou-Zeid2019), cylinders (Agbaglah and Mavriplis2019), and spheres (Ye et al.2016) exhibits different transport behavior. However, most buildings are a combination of cuboids (Oke et al.2017), and non-urban roughness flow is not suitable to inform urban parameterization. Therefore, the transport behavior over non-cuboid roughness is outside of the current scope. Note that the profiles over other designs with building height variability (not shown here) share similar characteristics. Therefore, we chose the uniform-height layout to better reveal the difference among flow properties. Within the canopy, the turbulent fluxes of momentum and TKE are predominantly negative, reflecting downgradient transports driven by their differing production and destruction mechanisms. The turbulent fluxes of passive scalars are positive due to the surface emission being the only source and exhibiting a minor variation across different densities. The eddy diffusivity of scalars (Kc) and momentum (Km) maintains a similar shape as a result of the relatively simpler mechanism in their production, destruction, and transport, which justifies the simplification (KcKm) that is consistent with Li and Bou-Zeid (2019). However, the eddy diffusivity of TKE (Kk) demonstrates a greater complexity with some negative values at the ground and canopy interface that marks the counter-gradient transport of TKE indicating two vital local shear production regions: (a) the canopy-top shear zone is relatively well-studied and is induced by clearings of roughness elements (Giometto et al.2016) and the (b) near-ground shear zone that is unique in medium-dense layouts (λp≤0.25) where flow reversal appears and counteracts with the pressure gradient, inducing a significant shear source of TKE (Lu et al.2023a). This lack of consideration of the near-ground flow features partially explains the underestimation of local TKE seen in Nazarian et al. (2020) and can be mitigated by imposing a height-dependent model constant (e.g., Glazunov et al.2021), but is beyond the scope of this study that aims to refine the turbulent transport characterization within the whole canopy. Apart from the near-ground and canopy-top inversion, Kk shares a similar shape with Km but presents a significantly larger magnitude, which indicates that models using unified K underestimate the transport of TKE and require closer characterization. The distinct behavior of Km in the densest case (purple lines in Fig. 1c) comes from the positive gradient of streamwise velocity in the middle of the canopy and positive Reynolds stress in the lower canopy that is hard to characterize. Such a high density is not rare, especially in some European and South American cities. However, realistic neighborhoods with extremely high density typically exhibit varied horizontal arrangements, leading to different flow characteristics from idealized building blocks (Lu et al.2023b).

Figure 1Vertical profiles of turbulent flux (first column), vertical gradient (second column), and eddy diffusivity (third column) of momentum, TKE, and scalar flux from Eq. (7) evaluated over uniform-staggered configurations. The yellow-shaded range indicates the mean building height. The x axis is limited to show most profiles for better presentability. Values are off the axis limit at the canopy top due to the strong shear zone.


Following the parameterization of length scales in Nazarian et al. (2020), we average the profile of eddy diffusivities shown in Fig. 1 up to Hmean to reveal the integral behavior of transport behaviors across configurations. Figure 2 presents the scatter of the vertically averaged eddy diffusivity for momentum Km,T and TKE Kk,T without considering dispersive flux, and these include the dispersive contribution (Km,T+D and Kk,T+D). Characterization of dispersive fluxes within the K-theory framework enhances the exchange rate of flow properties, which agrees with previous studies such as Akinlabi et al. (2022) where the dispersive flux is mostly downward.

The magnitude of eddy diffusivity for TKE is ∼3.5 times higher than that for momentum (note the different vertical ranges for momentum and TKE in Fig. 2). This considerable difference demonstrates a significant contrast to the parameterization strategy in Nazarian et al. (2020) that assumed Km,T+D=Kk,T+D, which results in underestimations of TKE within the canopy. In addition, unlike Km,T+D, which generally decreases with urban density before the layout becomes very dense, Kk,T+D exhibits a less decisive trend to density. The magnitude difference and variations over densities between Km,T+D, and Kk,T+D urge for a distinction in their parameterizations.

Figure 2Canopy-averaged eddy diffusivity for momentum and TKE with (Km,T+D and Kk,T+D) and without (Km,T and Kk,T) additional consideration of the dispersive component. Black lines indicate second-order polynomial regressions over the means of seven configurations under the same densities where the regression equation is shown at the bottom of each subplot.


The eddy diffusivity was not directly parameterized under the 1.5-order closure framework developed in Martilli et al. (2002) but as a turbulent length scale correlated with the dissipation length scale in Eq. (6) (i.e., Cklk=CμlεCε). Here we follow the same approach as in Nazarian et al. (2020) to parameterize the turbulent viscosity (Cμ) and dissipation length scale (lεCε) to obtain the turbulent length scale (Cklk). Considering that the parameterization for dissipation length scale is the same across the scenarios considered in Table 2, the difference in the transport efficiency demonstrated in Sects. 3.1 and 3.2 can be expressed as the difference in the turbulent viscosity (Cμ), as shown in Fig. 3.

Figure 3 shows the turbulent viscosity calculated in 49 idealized urban arrays binned over densities for better presentation, where the original parameterization with seven uniform-height cases (Nazarian et al.2020) is shown in green lines. Generally, Cμ increases with density over sparse and medium layouts, which indicates the overweighting trend of viscous dissipation of TKE over the turbulent transport that is consistent with the transition of the flow regime (Oke et al.2017). Meanwhile, denser layouts amplify the impact of building height variability, leading to a higher standard deviation among binned scenarios, which is consistent with Lu et al. (2023a). Turbulent viscosity for momentum lumping the dispersive flux Cμ,Km,T+D shows a good agreement with that from Nazarian et al. (2020) except for the peaking at higher density. Due to the separate explicit characterization of the building-induced dispersive momentum transport, Cμ,Km,T+S is slightly lower than Cμ,Km,T+D but shows a similar shape. The turbulent viscosity of TKE (Cμ,Kk,T+D) is about 3.5 times higher than that of momentum (note the difference on the y-axis limit) and shows exponential growth over density.

Figure 3Canopy-averaged turbulent viscosity for momentum and TKE. Panel (a) shows the turbulent viscosity of momentum where the red and black dots indicate evaluations adding the dispersive and turbulent fluxes (Cμ,TKm) and sampling out the fraction of building-induced dispersive flux (Cμ,SKm) and lumping the residue, respectively. Panel (b) shows the turbulent viscosity of TKE, lumping the dispersive component. Data points are binned by their ground density λp(z=0), where the vertical line attached to each square indicates the standard deviation among the data points binned. Green lines indicate the parameterization from Nazarian et al. (2020).


3.2 Characterization of the dispersive momentum transport

The momentum transport in urban canopy flow is also modulated by its spatial variability (Eq. 1) induced by the rigid volume of roughness elements that is hard to characterize due to its extreme heterogeneity (Poggi and Katul2008b). Nevertheless, flow heterogeneity is not unique to the urban boundary layer but is prevalent in the real-world atmosphere layer due to the surface heating and has been well characterized in PBL parameterizations (e.g., Han and Bretherton2019). In the context of convective boundary layer parameterization, the eddy diffusivity (ED) approach based on K theory (as discussed in Sect. 3.1) has been successful in representing neutral and surface boundary layers that promote more homogeneous and isotropic turbulence (Fig. 4a). On the other hand, the flow heterogeneity induced by convection from thermal forcing and moisture phase change has been successfully captured by the mass-flux (MF) approach (Siebesma et al.2007). In the effort to model both atmospheric conditions coexisting in PBL, an eddy-diffusivity mass-flux (EDMF) scheme has been developed by partitioning flow fluctuations into contributions from local turbulent mixing (ED) and non-local coherent rising and descending air parcels (MF, Lopez‐Gomez et al.2020; Lu2019).

This section explores this alternative approach to characterize the dispersive momentum flux, other than the pragmatic approach developed in Nazarian et al. (2020); Simón-moral et al. (2016), by considering dispersive transports as mass-flux components under the framework of the EDMF scheme. The dispersive transports of momentum (Fig. 4b) in the urban boundary layer are induced by coherent motions near the building facets and present a strong local distribution that resembles the patchy coherent structures induced by thermal forcing in PBL. Therefore, the similarity in flow heterogeneity between UCL and PBL implies a phenomenological analogy, which favors the introduction of the EDMF approach to multi-layer UCMs. Following the formulation in Siebesma et al. (2007), the dispersive flux of momentum is firstly decomposed into individual plumes (ai) enclosing coherent motions with strength ui-uwi-w and then sorted into upward (positive, with a fraction of au) and downward (negative, with a fraction of ad) contributions. The dispersive flux of momentum then reads

(8) u ̃ w ̃ = i a i u i - u w i - w = a u u ̃ w ̃ u + a d u ̃ w ̃ d a u + a d = 1 .

Note that, for simplicity, the entrainment process (Abdella and Petersen2000) that represents scalar transport between subdomains (between dispersive structures and the turbulent environment) is not considered here. Using the mass-flux approach to characterize dispersive momentum flux requires the fraction (au+ad=1) and strength (ũw̃u, ũw̃d) of those coherent flow motions to be parameterized based on the variation of the urban arrangement. An example of the spatial distribution of turbulent and dispersive momentum flux is shown in Fig. 4 for the S-highrise-SD56 configuration of density λp=0.25 sampled at z=15 m.

Figure 4An example of horizontal and vertical field sampled at 93.8 % of the mean building height (z=15 m and y=110 m, respectively), for the S-highrise-SD56 configuration of density λp=0.25. The left panels (a, c) show the turbulent flux of momentum ((a) uw/uτ2=-0.6052) that forms the “background” transport events, with the title showing the spatial average. The right panels (b, d) show the dispersive flux that forms “hot spot” transport events ((b) adũw̃d/uτ2=-0.2122, auũw̃u/uτ2=0.0253), with the contribution from upward and downward transport indicated in the title. Regions enclosed by red and blue lines indicate values sampled with a 50 % CDF. Vectors in the bottom left (c) panel show the mean velocity direction, whereas those in the bottom right show the direction of the dispersive velocity (ũ=u-u).


Unlike the turbulent momentum flux, which shows a uniform downward transport with an indecisive correlation with the local roughness elements, a portion of the dispersive counterpart is directly connected to the local roughness, whereas the rest remains relatively homogeneous. Therefore, a proper sampling filter should be applied to enclose these regions before gathering flow properties for parameterization. Again, we adopt a common sampling approach in the development of the EDMF (Couvreux et al.2010; Sušelj et al.2012) based on the joint cumulative probability density function (CDF) of the dispersive velocities (ũ,w̃).

(9) CDF ( ( u ̃ , w ̃ ) ) = u ̃ min u ̃ max w ̃ min w ̃ max PDF ( u ̃ , w ̃ ) d u ̃ d w ̃ = 1

The CDF can be adjusted from [0,1] that corresponds to different integration ranges of dispersive streamwise and vertical velocity (Eq. 9). By lowering the CDF, the sampled regions are adjusted smaller only to enclose regions with outstanding strength. We tested a range of cut-off values of CDF[0.3,0.4,0.5,0.6,0.7] (test results are shown in the Appendix and in the Supplement), and a cut-off value here is determined as 0.5 for all cases considered. Figure 4b shows the sampled upward and downward transports enclosed by red and blue lines, respectively.

The sampling filter aims to capture the building-induced component of the dispersive momentum flux that K theory cannot characterize. Here, we test the correlation between dispersive streamwise velocity (ũ) and vertical velocity (w̃) to show how the urban structures modify the statistical distribution of dispersive fluxes and the corresponding characterization fulfilled by the sampling. Figure 5 shows the 2D CDF distribution for S-highrise-SD56 configuration of density λp=0.25. The distribution exhibits a more eccentric behavior along the vertical range, with

  1. One tail growing stronger and thinner at the third quadrant (ũ<0,w̃<0) representing the blockage effects of the frontal area of the roughness (Fig. 4d), which is well-captured by the sampling strategy (red squares in Fig. 5). The height distribution of building blocks can explain the thinner trend and reflect how the dispersive flux responds to the urban vertical structure.

  2. A relatively symmetric elongated region at the second (ũ<0,w̃>0) and fourth quadrant (ũ>0,w̃<0) representing the leeward flow cavity (Fig. 4b, characterized by the sampled flux at most heights) and flow divergence at the lateral facets (represented by the sampled flux at 2.5 m), respectively. The elongation of both sides reflects a more irregular distribution of downward dispersive flux with heights, which is consistent with the analysis in Lu et al. (2023a).

Figure 5Contour plot showing the two-dimensional histograms of the CDF of the dispersive velocities at 2.5, 5, 7.5, 10, 15, and 22 m. The integral of the PDF (from outside to the mean value) within the area is delineated by isolines on the logarithmic scale. The black/red squares depict the averaged dispersive flux over the sampled regions to represent the strength of the gradient and counter-gradient transport parameterized in MLUCM v3.0.


After subtracting dispersive structures from the sampled regions, the rest of the relatively homogeneous regions may safely resort back to the lumping approach and update Eq. (8) with Eq. (10):

(10) u ̃ w ̃ = a u u ̃ w ̃ u co = 0.5 + a d u ̃ w ̃ d co = 0.5 + u ̃ w ̃ residue .

We speculate the reason that dispersive flux has long been ill represented in UCMs is attributed to the inappropriate scaling with the turbulent flux (e.g., Poggi and Katul2008b), which usually yields a small fraction without recognizing these two processes are controlled by very different physical processes. Alternatively, the dispersive perturbations are caused by the pressure fields associated with the topography or the changed resistance to flow through the canopy (Finnigan et al.2015; Finnigan and Shaw2008), which implies a possibility of scaling its strength to the external pressure gradient. This was only made possible recently due to the increased computational power and awareness of intrinsic urban heterogeneity (Boeing2019; Lu et al.2023b). Therefore, we hypothesize the strength of sampled flow motions can be parameterized by pressure gradient and urban morphological parameters as Eq. (11):

(11) a u , d u ̃ w ̃ u , d = f z , p x , λ p , H max , mean , std .

Figure 6Vertical profiles of sampled downward and upward dispersive momentum flux (gray-scale lines) for seven building height configurations covering λp[0.0625,0.4444]. The range of mean building height is marked in yellow shades with the transparency following the λp. Diminished red lines indicate the strength of the external pressure gradient constant with height. Dashed red lines show the parameterization from Eq. (12).


To gain a direct impression of how the sampled dispersive motions scale with the pressure gradient, we present the profile of the vertical gradient of the sum of two sampled terms from Eq. (10) in Fig. 6. The variation of the dispersive momentum transport exhibits a strong dependence on the flow regime but decouples with the turbulence level (not shown here; see Figs. 9–11 from Lu et al. (2023b) for spatial pattern of TKE as well as turbulent and dispersive momentum flux). The sampled dispersive motion is of a similar magnitude to the only driver in our experiments – pressure gradient collectively leaving such a term unattended is equivalent to only applying half of the pressure gradient to the flow and inevitably degrades model predictability. Considering that realistic urban geometries will exhibit increased complexity and dispersive flux (Giometto et al.2016), explicit parameterization is necessary. The little variation in the profiles for each density in Fig. 6 indicates building height distribution does not induce significant variance over the sampled dispersive flux among cases within the Hmean. Instead, urban density dictates the variation of the sampled flux, whose profile is different between sparse (the first row in Fig. 6) and denser (the second row in Fig. 6) layouts: for the three sparse layouts, a relatively linear increasing trend crossing zero with heights is developed. Denser layouts lead to strong overlapping of building wakes, which helps to maintain a constant strength profile of similar magnitude to the pressure gradient, which indicates a trend of a skimming flow regime (Oke et al.2017). Above the Hmean, sampled fluxes quickly converge to zero except for the high-rise configurations over denser layouts (λp≤0.25) and at the Hmax of each configuration, demonstrating a minor magnitude but complex pattern across configurations. Therefore, for simplicity, we only provide a parameterization for the sampled fluxes below Hmean that scales with the external pressure gradient dpdx as follows:


The above parameterization retains the physical significance of upward and downward dispersive motion demonstrated in Fig. 6. By definition, the scaling coefficients to the external pressure gradient are designed to be lowered over denser layouts, consistent with previous findings where dispersive momentum flux is more significant over sparse arrangements (e.g., Lu et al.2023b).

3.3 Drag and dissipation parameterization

The presence of roughness elements on the urban surface implicitly modifies turbulent length scales (Li et al.2020) that regulate the transport mechanisms of flow properties. Explicitly, the buildings block and divert the flow, inducing pressure deficits between their windward and leeward facets that collectively serve as a sink for momentum (Eq. 1) and a source for TKE (Eq. 4). In the following, we revisit the modeling strategy for drag and dissipation in Nazarian et al. (2020) and add the necessary adjustments to reflect the building height variability. By assuming the strength of building effects is proportional to the area facing the wind per cubic meter of outdoor air volume S(z) and the spatially averaged mean wind speed U, and a drag coefficient Cd encodes the variation of building arrangements, the drag parameterization reads


Figure 7Drag parameterization based on the equivalent drag coefficient Cdeq. Panel (a) shows the error bar plot where black squares mark the binned value of Cdeq among seven building height configurations, and the vertical line attached to each square indicates the standard deviation for each bin. The green dots indicate Cdeq evaluated over uniform-staggered cases from Nazarian et al. (2020). Panel (b) shows the scatter plot of the ratio between Cdeq and the drag coefficient evaluated above the mean building height where the dashed line indicates the ratio of 1. The green-shaded area indicates the corresponding parameterization in the plot.


Accordingly, characterizing drag impacts converts to evaluating the drag coefficient Cd that comes in different forms in the literature (e.g., Santiago et al.2013). The equivalent drag coefficient Cdeq evaluated following Santiago and Martilli (2010b) requires fewer inputs but maintains a similar performance, as demonstrated by its capacity in MLUCM v2.0 (Nazarian et al.2020). It will be retained in the present study. Using Cdeq for the evaluation assumes a constant profile of drag that ensures the integral of the drag over the whole urban canopy is equivalent to that evaluated from LES simulations:

(14) C deq = - 1 ρ H 0 H Δ p ( z ) d z 1 H 0 H u ( z ) | u ( z ) | d z ,

where Δp(z)=pf(z)-pb(z) is the pressure deficit between the windward and leeward facets of each building block. Due to the low variance of Cdeq over the impact of height variability except for denser layouts λp≥0.35, we binned all seven cases over each density in Table 2 to show their mean and standard deviation in Fig. 7 together with that from Nazarian et al. (2020) for uniform-staggered arrays. The excellent agreement with Cdeq in the present study indicates the parameterization in Nazarian et al. (2020) presents great resilience to the height variability, especially over the sparse layouts. Therefore, updating the constants in the original parameterization requires only a minor change. However, height variability not only extends the vertical range of the urban canopy but also complicates the vertical structure of the drag coefficient and requires extra characterization. Figure 7b shows the ratio between the drag coefficient above the mean building height Hmean and Cdeq that covers a range that is generally larger than unity and shows a decreasing trend with densities, which can be captured by considering the morphological condition of the urban surface from Lu et al. (2023a).

The last term parameterized is the dissipation length scale from Eq. (5). As indicated in the left panel of Fig. 8, the impact of height variability lε/Cε below the mean building height only demonstrates a very minor variation for each density range. The variation is more evident above Hmean, with a clear shape following the sectional density, and was discussed in Lu et al. (2023a). Despite the fact that applying a sophisticated dissipation and drag parameterization might lead to better model performance (Glazunov et al.2021; Castro2017; Sützl et al.2021a), it is beyond the scope of the present study, which is focused on refining the transport characterization of UCM. Therefore, we decided to keep the original parameterization from MLUCM v2.0 (Nazarian et al.2020). For completeness, we briefly describe the parameterization strategy: by segmenting the urban canopy flow into the mixing layer, transition region, and wall-bounded layer (Coceal and Belcher2004), the parameterization reads


where α1=5.5 is a revised value based on 49 simulations in the present study, and d=Hmeanλp0.15 is the displacement height parameterized following Krayenhoff et al. (2015). Equation (15a) depicts a mixing layer type of flow with a constant dissipation efficiency and Eq. (15b) represents the transition to a fully wall-bounded flow in Eq. (15c) where the dissipation efficiency decreases with height based on d2 and α2 parameterized as follows:


Figure 8The left panels ((a)(c),  (e)(g)) show the vertical profiles of dissipation length scales lε/Cε up to 2.5Hmean where the yellow-shaded range indicates mean building height. The right panel (d) shows the binned vertical averaged dissipation length scale lε/Cε in the same manner as Fig. 7 overlaid with the parameterization and with that from Nazarian et al. (2020).


4 MLUCM performance with updated parameterizations

In this section, we evaluate the predictability of two modifications to the original MLUCM developed in Nazarian et al. (2020), where the combination of changes proposed in Sects. 3.1 and 3.2 and the original model form four scenarios, as shown in Table 2.

Table 2Four scenarios tested in the present study with the combination of the EDMF approach (Sect. 3.2) and differentiation in the TKE diffusion (Sect. 3.1 ). Here, 1D-Original aims to set a baseline by employing a similar parameterization from Nazarian et al. (2020) with an adaptation for the height variability from Sect. 3.3. Therefore, the 1D-Original does not differentiate turbulent TKE flux from turbulent momentum flux and lumps the dispersive momentum flux into the calculation of Km (To Km). 1D-Kk considers the differentiation of TKE transport from momentum as discussed in Sect. 3.1. By contrast, 1D-EDMF only considers the “EDMF” parameterization for dispersive momentum flux while 1D-KkEDMF considers both modifications. The bold font in the table indicates modified components from the original MLUCM.

Download Print Version | Download XLSX

Figure 9Vertical profiles of velocity (u/uτ), turbulent kinetic energy (k/uτ2), and turbulent momentum flux (uw/uτ2) obtained with the four scenarios in Table 2 and LES results for sparse (λp=0.1111), medium (λp=0.25), and dense (λp=0.4444) layouts.


Figure 9 shows the vertical profiles of horizontally averaged velocity, TKE, and turbulent momentum flux calculated with the one-dimensional models with different transport characterization and compared with the corresponding LES results for the high-rise configuration with Hstd=5.6 m. The profile comparison for other configurations can be found in the Supplement. We select sparse (λp=0.1111), medium (λp=0.25), and dense (λp=0.4444) layouts to reflect different flow regimes for comparison.

The prediction of streamwise velocity does not present a great variation, where 1D-Kk-EDMF exhibits the highest value. This is explained by the newly introduced mass-flux term as an explicit source of momentum and differentiation of Kk from Km enhancing the momentum transport from the atmosphere to the canopy.

For TKE estimation, 1D-Kk and 1DKk-EDMF with differentiation of Kk from Km (Table 2) demonstrate a significant improvement on the TKE profile for medium and dense layouts below Hmean. The introduction of EDMF is phenomenologically equivalent to imposing a stronger pressure gradient, which leads to higher velocity and explains the better agreement with LES for the 1D-EDMF and 1D-Kk-EDMF scenarios. As a result, 1D-Kk-EDMF yields an even higher TKE due to the introduction of a mass-flux term to parameterize dispersive momentum flux as an explicit source leading to higher TKE production. The correction of TKE not only leads to better modeling of the velocity profile but also enhances the prediction of urban fluxes such as the overestimation of daytime air temperature by correcting the vertical exchange of heat (Krayenhoff2014).

The variation in the performance of scenarios is not pronounced across different height variability, which further consolidates the binning strategy for parameterization in the present study. To evaluate the comprehensive performance of scenarios over different height configurations, Fig. 10 presents the root mean square error (RMSE) between vertical profiles (ranging from z=[0-3Hmean]) of LES and four MLUCM scenarios averaged over seven configurations (one uniform and six variable height, as shown in Table 1) over seven densities.

Figure 10Difference in the root mean square error (RMSE) between the original MLUCM v2.0 (original) and three improved scenarios (positive values indicate better performance from scenarios). The RMSE is evaluated between the one-dimensional model from Table 2 against LES across seven densities in combination with seven height configurations from Table 1. Each column shows the averaged RMSE difference across seven height configurations.


The new parameterizations with refined transport characterization represent an overall improvement over the previous multi-layer model (Nazarian et al.2020). All three scenarios with modifications lead to better predictability for velocity except for the two sparse layouts. The predictability for TKE is substantially better from the 1D-Kk and 1D-KkEDMF scenarios, except for two sparse layouts that can be further improved with a better characterization of the TKE budget terms such as dissipation and wake production (Blunn et al.2022). The turbulent momentum flux is diagnostically obtained in MLUCM based on the kl closure (uw=-Cklkk1/2uz) that is highly related to streamwise velocity and TKE. However, it exhibits a consistently contrasting performance trend compared to TKE, which leads to less improvement in denser layouts. This pattern can be attributed to both improvements (Table 2) enhancing the vertical exchange of flow properties within the canopy, which results in a smaller vertical gradient of streamwise velocity. Having the most remarkable performance, 1D-KkEDMF with both modifications from Sects. 3.1 and 3.2 is desirable for MLUCM v3.0.

5 Summary and conclusions

This study refined the characterization of the transport of flow properties for the multi-layer urban canopy model via a separate parameterization of TKE diffusion and introduced a “mass-flux” term for the dispersive transport of momentum. The updated multi-layer model demonstrates an improved performance by correcting the underestimation of turbulent exchange and velocity below the mean building height Hmean.

By analyzing 49 LES simulations over staggered urban arrays of uniform and variable building height and comprehensive density coverage, we found the turbulent exchange rate of momentum is similar to that for scalars but 3.5 times lower than TKE. However, the previous model (MLUCM v2.0, Nazarian et al.2020) assumed a unified transport efficiency for momentum and TKE, which curbs the transport of the flow above Hmean into the canopy and causes an underestimation of in-canopy TKE and wind speed. The accurate characterization of TKE transport into the lower canopy becomes critical as it also controls the turbulent exchange of heat and moisture between buildings and vegetation.

Regarding applications, the enhanced model, which addresses the underprediction of in-canopy TKE, has the potential to alleviate the overestimation of daytime air temperature and diurnal temperature range (Krayenhoff et al.2020). Furthermore, the improved prediction of urban fluxes reinforces the advantages of MLUCM over simpler bulk schemes, emphasizing their ability to realistically simulate vertical diffusion processes in mesoscale atmospheric models (Hendricks et al.2020). The proposed refinements are not exclusively constrained to urban applications but are also applicable to broader turbulence modeling. For example, the explicit parameterization of dispersive momentum flux could benefit turbulence modeling over other types of roughness elements such as coral reefs (Davis et al.2021) and vegetative canopies (Finnigan2000). In parallel, emphasizing the distinction between TKE transport efficiency and its momentum counterpart underscores the importance of gathering higher flow moments for modeling general atmospheric flow (Ghonima et al.2017).

We also revealed the distinct horizontal distribution of dispersive momentum flux, which is induced by flow heterogeneity responding to windward (quadrant 3, upward), lateral (quadrant 4, downward), and leeward (quadrant 2, downward) flow patterns. In response to this complication, we applied a sampling filter to the cumulative density function of the dispersive transport to segment the building-induced contribution (MF) that scales well with the pressure gradient and the rest of the relatively undisturbed component that can be securely lumped into the turbulent counterpart in the 1.5-order turbulent closure (ED). The EDMF framework, having successfully modeled the planetary boundary layer with flow heterogeneity induced by thermal effects, was demonstrated to be also favorable for the representation of flow heterogeneity caused by rigid building volumes. The new framework enables an explicit consideration of the impact of flow heterogeneity in MLUCM, improving our ability to calibrate the model for simulations of complex canopy flow.

With adaptations to the parameterizations of drag coefficient and dissipation length scale to account for building height variability, the model configuration that includes both the updated eddy diffusivity of TKE (Kk) and the EDMF scheme yields an improved performance relative to MLUCM v2.0. However, analysis of the flow profiles and length scales indicates future investigation is needed to address the following aspects:

  1. The turbulent length scales show a great variance over heights for layouts with non-uniform building height. A height-dependent parameterization of turbulent length scale and drag profile is necessary to correctly simulate features of urban flow such as the logarithmic or exponential behavior of the wind speed profile (Castro2017).

  2. The dispersive transport of TKE and scalars can be considered similarly following the dispersive momentum flux parameterization. In the original implementation of EDMF in PBL parameterizations (Siebesma et al.2007), the sampled air parcel also exchanges flow properties with the turbulent environment and other air parcels, which could be considered in future developments.

  3. Staggered building arrays were found to be less representative of flow over realistic dense urban neighborhoods (Lu et al.2023b). The wind direction (Santiago et al.2013) and thermal effects (Simón-moral et al.2016) also affect dynamic urban flow parameterizations. Further simulation for model calibration should account for these factors for a more realistic flow characterization.

  4. Similar to the scarcity of evaluations of the eddy diffusivity of flow properties other than momentum, the variation of the TKE budget terms is still poorly considered in urban flow models. Accordingly, a more comprehensive analysis that addresses the relationships between the urban morphology and the TKE budget may provide further improvements.

  5. The applicability of the 1.5-order kl turbulence closure in MLUCM is restricted to isotropic turbulence flow, where the momentum exchange rate is directly correlated with turbulence strength. However, this scaling becomes questionable even in neutral conditions and diminishes under thermal stratification, where turbulence tends to be highly anisotropic (Stiperski and Calaf2018). In response to these limitations, employing a higher-order closure that diagnostically assesses momentum transfer (Wilson1988; Ayotte et al.1999) or considering alternative scaling approaches for poorly performing regions (Sun et al.2020) may contribute to the development of more effective urban canopy models in the future.

Appendix A: Sampling of the dispersive momentum flux

The sampling procedure to extract the components of building-induced dispersive fluxes was conducted as follows: (a) separating the dispersive fluxes into gradient (downward) and counter-gradient (upward) components, (b) sorting these two fluxes into a descending order in terms of the magnitude, (c) lowering the cumulative density function of the two fluxes from the largest allows a smaller subdomain to be sampled until subdomains only enclose transport events connected to buildings, and (d) spatially averaging the sampled dispersive fluxes over subdomains for both gradient and counter-gradient dispersive momentum flux. Figure A1 shows a range of cut-off values (co[0.3,0.4,0.5,0.6,0.7]) for CDF over an example for λp=0.25. The cut-off value here is determined as 0.5 for all cases considered (an example is provided in Fig. 4b, where blue (red) lines enclose the sampled downward (upward) transports).

The three-dimensional distribution of the sampled dispersive structures is shown in Fig. A2 S-highrise-SD56 for λp[0.0625,0.4444] where the dispersive momentum flux concentrates in the immediate vicinity downstream of the cubes. Structures near the ground contributed significantly to negative stresses (green structures) where the recirculating flows dominate. In response to this distribution (note the dispersive velocity vectors in Fig. 4d), the parameterization for the vertical variation of sampled dispersive fluxes is stronger over the lower part of the canopy over sparse layouts, which is in general agreement with Poggi and Katul (2008a). In the upper levels of the canopy, the counter-gradient stresses (red structures) are more pronounced at the windward facets due to flow reflection.

Figure A1Sampling regions of dispersive fluxes for five (co[0.3,0.4,0.5,0.6,0.7]) cut-off CDF values as discussed in Sect. 3.2. Higher cut-off values for CDF lead to broader sampling regions and lower mean strength in the mass-flux term.


Figure A2Three-dimensional distribution of the dispersive coherent structures for cases with S-highrise-SD56 covering λp[0.0625,0.4444]. Green and red regions represent gradients and counter-gradient dispersive transport sampled with a 50 % CDF.


Code availability

The source code and the supporting data of three scenarios for the 1D Multi-layer Urban Canopy Model are publicly available at (Lu et al.2023c) under GPL 3.0 license.

Data availability

The dataset is available in the Supplement.


The supplement related to this article is available online at:

Author contributions

JL, NN, MAH, ESK, and AM collectively developed and planned the study. JL ran the large-eddy simulations designed by NN and modified the model with the help of NN and AM. JL carried out the result analyses and wrote the paper with significant input and critical feedback from NN, MAH, ESK, and AM.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Simulations were undertaken with the assistance of resources and services from the Australian Government's National Collaborative Research Infrastructure Strategy (NCRIS), with access to computational resources provided by the National Computational Infrastructure (NCI), which is supported by the Australian Government through the National Computational Merit Allocation Scheme. We also thank the anonymous reviewers for their constructive and valuable comments.

Financial support

This research has been supported by the Australian Research Council Centre of Excellence for Climate Extremes (Grant CE170100023).

Review statement

This paper was edited by Yongze Song and reviewed by two anonymous referees.


Abdella, K. and Petersen, A.: Third-Order Moment Closure Through A Mass-Flux Approach, Bound.-Lay. Meteorol., 95, 303–318,, 2000. a

Agbaglah, G. and Mavriplis, C.: Three-dimensional wakes behind cylinders of square and circular cross-section: early and long-time dynamics, J. Fluid Mech., 870, 419–432,, 2019. a

Akinlabi, E., Maronga, B., Giometto, M. G., and Li, D.: Dispersive Fluxes Within and Over a Real Urban Canopy: A Large-Eddy Simulation Study, Bound.-Lay. Meteorol., 185, 93–128,, 2022. a, b

Ayotte, K. W., Finnigan, J. J., and Raupach, M. R.: A Second-Order Closure for Neutrally Stratified Vegetative Canopy Flows, Bound.-Lay. Meteorol., 90, 189–216,, 1999. a

Best, M. J. and Grimmond, C. S. B.: Key Conclusions of the First International Urban Land Surface Model Comparison Project, B. Am. Meteorol. Soc., 96, 805–819,, 2015. a

Blunn, L. P., Coceal, O., Nazarian, N., Barlow, J. F., Plant, R. S., Bohnenstengel, S. I., and Lean, H. W.: Turbulence Characteristics Across a Range of Idealized Urban Canopy Geometries, Bound.-Lay. Meteorol., 182, 275–307,, 2022. a, b, c

Boeing, G.: Urban spatial order: street network orientation, configuration, and entropy, Appl. Network Sci., 4, 67,, 2019. a

Bougeault, P. and Lacarrere, P.: Parameterization of Orography-Induced Turbulence in a Mesobeta–Scale Model, Mon. Weather Rev., 117, 1872–1890,<1872:POOITI>2.0.CO;2, 1989. a

Brown, M., Lawson, R., DeCroix, D., and Lee, R.: Comparison of Centerline Velocity Measurements Obtained Around 2D and 3D Building Arrays in a Wind Tunnel, article, 1 July 2001, New Mexico, University of North Texas Libraries, UNT Digital Library, (last access: 2 April 2024), 2001. a

Castro, I. P.: Are Urban-Canopy Velocity Profiles Exponential?, Bound.-Lay. Meteorol., 164, 337–351,, 2017. a, b, c

Cheng, W.-C. and Porté-Agel, F.: A Simple Mixing-Length Model for Urban Canopy Flows, Bound.-Lay. Meteorol., 181, 1–9,, 2021. a

Christen, A., Rotach, M. W., and Vogt, R.: The Budget of Turbulent Kinetic Energy in the Urban Roughness Sublayer, Bound.-Lay. Meteorol., 131, 193–222,, 2009. a

Coceal, O. and Belcher, S. E.: A canopy model of mean winds through urban areas, Q. J. Roy. Meteor. Soc., 130, 1349–1372,, 2004. a

Coceal, O., Thomas, T. G., Castro, I. P., and Belcher, S. E.: Mean Flow and Turbulence Statistics Over Groups of Urban-like Cubical Obstacles, Bound.-Lay. Meteorol., 121, 491–519,, 2006. a

Coceal, O., Dobre, A., Thomas, T. G., and Belcher, S. E.: Structure of turbulent flow over regular arrays of cubical roughness, J. Fluid Mech., 589, 375–409,, 2007. a

Couvreux, F., Hourdin, F., and Rio, C.: Resolved Versus Parametrized Boundary-Layer Plumes. Part I: A Parametrization-Oriented Conditional Sampling in Large-Eddy Simulations, Bound.-Lay. Meteorol., 134, 441–458,, 2010. a

Davis, K. A., Pawlak, G., and Monismith, S. G.: Turbulence and Coral Reefs, Annu. Rev. Mar. Sci., 13, 343–373,, 2021. a

Deardorff, J. W.: Stratocumulus-capped mixed layers derived from a three-dimensional model, Bound.-Lay. Meteorol., 18, 495–527,, 1980. a

Finnigan, J.: Turbulence in plant canopies, Annu. Rev. Fluid Mech., 32, 519–571, 2000. a, b

Finnigan, J., Harman, I., Ross, A., and Belcher, S.: First-order turbulence closure for modelling complex canopy flows: First-Order Turbulence Closure for Canopy Flows, Q. J. Roy. Meteor. Soc., 141, 2907–2916,, 2015. a

Finnigan, J. J.: Turbulent Transport in Flexible Plant Canopies, Springer Netherlands, Dordrecht, 443–480, ISBN 978-94-009-5305-5,, 1985. a

Finnigan, J. J. and Belcher, S. E.: Flow over a hill covered with a plant canopy, Q. J. Roy. Meteor. Soc., 130, 1–29,, 2004. a

Finnigan, J. J. and Shaw, R. H.: Double-averaging methodology and its application to turbulent flow in and above vegetation canopies, Acta Geophys., 56, 534–561,, 2008. a

Frigo, M. and Johnson, S.: FFTW: an adaptive software architecture for the FFT, in: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. No.98CH36181),IEEE, Seattle, WA, USA, ISBN 978-0-7803-4428-0, vol. 3, 1381–1384,, 1998. a

Ghonima, M. S., Yang, H., Kim, C. K., Heus, T., and Kleissl, J.: Evaluation of WRF SCM Simulations of Stratocumulus‐Topped Marine and Coastal Boundary Layers and Improvements to Turbulence and Entrainment Parameterizations, J. Adva. Model. Earth Sy., 9, 2635–2653,, 2017. a

Giometto, M. G., Christen, A., Meneveau, C., Fang, J., Krafczyk, M., and Parlange, M. B.: Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface, Bound.-Lay. Meteorol., 160, 425–452,, 2016. a, b, c

Glazunov, A. V., Debolskiy, A. V., and Mortikov, E. V.: Turbulent Length Scale for Multilayer RANS Model of Urban Canopy and Its Evaluation Based on Large-Eddy Simulations, Supercomputing Frontiers and Innovations, 8,, 2021. a, b

Han, J. and Bretherton, C. S.: TKE-Based Moist Eddy-Diffusivity Mass-Flux (EDMF) Parameterization for Vertical Turbulent Mixing, Weather Forecast., 34, 869–886,, 2019. a

Harman, I. N., Böhm, M., Finnigan, J. J., and Hughes, D.: Spatial Variability of the Flow and Turbulence Within a Model Canopy, Bound.-Lay. Meteorol., 160, 375–396,, 2016. a

Hendricks, E. A., Knievel, J. C., and Wang, Y.: Addition of Multilayer Urban Canopy Models to a Nonlocal Planetary Boundary Layer Parameterization and Evaluation Using Ideal and Real Cases, J. Appl. Meteorol. Clim., 59, 1369–1392,, 2020. a, b

Kondo, H., Genchi, Y., Kikegawa, Y., Ohashi, Y., Yoshikado, H., and Komiyama, H.: Development of a Multi-Layer Urban Canopy Model for the Analysis of Energy Consumption in a Big City: Structure of the Urban Canopy Model and its Basic Performance, Bound.-Lay. Meteorol., 116, 395–421,, 2005. a

Krayenhoff, E. S.: A multi-layer urban canopy model for neighbourhoods with trees., Ph.D. thesis, University of British Columbia,, 2014. a

Krayenhoff, E. S., Santiago, J. L., Martilli, A., Christen, A., and Oke, T. R.: Parametrization of Drag and Turbulence for Urban Neighbourhoods with Trees, Bound.-Lay. Meteorol., 156, 157–189,, 2015. a

Krayenhoff, E. S., Jiang, T., Christen, A., Martilli, A., Oke, T. R., Bailey, B. N., Nazarian, N., Voogt, J. A., Giometto, M. G., Stastny, A., and Crawford, B. R.: Urban Climate A multi-layer urban canopy meteorological model with trees (BEP- Tree): Street tree impacts on pedestrian-level climate, Urban Clim., 32, 100590,, 2020. a, b, c

Kusaka, H., Kondo, H., and Kikegawa, Y.: A simple single-layer urban canopy model for atmospheric models: comparison with multi-layer and slab models, Bound.-Lay. Meteorol. 101, 329–358,, 2001. a

Li, Q. and Bou-Zeid, E.: Contrasts between momentum and scalar transport over very rough surfaces, J. Fluid Mech., 880, 32–58,, 2019. a, b, c, d

Li, Q., Bou‐Zeid, E., Grimmond, S., Zilitinkevich, S., and Katul, G.: Revisiting the relation between momentum and scalar roughness lengths of urban surfaces, Q. J. Roy. Meteor. Soc., 146, 3144–3164,, 2020. a

Lim, H. D., Hertwig, D., Grylls, T., Gough, H., Reeuwijk, M. v., Grimmond, S., and Vanderwel, C.: Pollutant dispersion by tall buildings: laboratory experiments and Large-Eddy Simulation, Exp. Fluids, 63, 92,, 2022. a

Lipson, M. J., Grimmond, S., Best, M., Abramowitz, G., Coutts, A., Tapper, N., Baik, J., Beyers, M., Blunn, L., Boussetta, S., Bou‐Zeid, E., De Kauwe, M. G., De Munck, C., Demuzere, M., Fatichi, S., Fortuniak, K., Han, B., Hendry, M. A., Kikegawa, Y., Kondo, H., Lee, D., Lee, S., Lemonsu, A., Machado, T., Manoli, G., Martilli, A., Masson, V., McNorton, J., Meili, N., Meyer, D., Nice, K. A., Oleson, K. W., Park, S., Roth, M., Schoetter, R., Simón‐Moral, A., Steeneveld, G., Sun, T., Takane, Y., Thatcher, M., Tsiringakis, A., Varentsov, M., Wang, C., Wang, Z., and Pitman, A. J.: Evaluation of 30 urban land surface models in the Urban‐PLUMBER project: Phase 1 results, Q. J. Roy. Meteor. Soc., 150, 126–169,, 2023. a, b

Lopez‐Gomez, I., Cohen, Y., He, J., Jaruga, A., and Schneider, T.: A Generalized Mixing Length Closure for Eddy‐Diffusivity Mass‐Flux Schemes of Turbulence and Convection, J. Adv. Model. Earth Sy., 12, e2020MS002161,, 2020. a

Lu, J.: Coherent Structures Sampling and a Stochastic Eddy-Diffusivity/Mass-Flux Parameterization for the Transition from Stratocumulus to Cumulus, Master's thesis, UC San Diego, (last access: 2 April 2024), 2019. a

Lu, J., Nazarian, N., Anne Hart, M., Krayenhoff, E. S., and Martilli, A.: Representing the effects of building height variability on urban canopy flow, Q. J. Roy. Meteor. Soc., 150, 46–67,, 2023a. a, b, c, d, e, f, g, h

Lu, J., Nazarian, N., Hart, M. A., Krayenhoff, E. S., and Martilli, A.: Novel Geometric Parameters for Assessing Flow Over Realistic Versus Idealized Urban Arrays, J. Adv. Model. Earth Sy., 15, e2022MS003287,, 2023b. a, b, c, d, e, f, g, h

Lu, J., Nazarian, N., Anne Hart, M., Krayenhoff, E. S., and Martilli, A: Code for “A one-dimensional urban flow model with an Eddy-diffusivity Mass-flux (EDMF) scheme and refined turbulent transport (MLUCM v3.0)” Submitted to GMD, Zenodo [code],, 2023c. a

Macdonald, R. W.: Modelling the mean velocity profile in the urban, Bound.-Lay. Meteorol. 97, 25–45,, 2000. a

Maronga, B., Banzhaf, S., Burmeister, C., Esch, T., Forkel, R., Fröhlich, D., Fuka, V., Gehrke, K. F., Geletič, J., Giersch, S., Gronemeier, T., Groß, G., Heldens, W., Hellsten, A., Hoffmann, F., Inagaki, A., Kadasch, E., Kanani-Sühring, F., Ketelsen, K., Khan, B. A., Knigge, C., Knoop, H., Krč, P., Kurppa, M., Maamari, H., Matzarakis, A., Mauder, M., Pallasch, M., Pavlik, D., Pfafferott, J., Resler, J., Rissmann, S., Russo, E., Salim, M., Schrempf, M., Schwenkel, J., Seckmeyer, G., Schubert, S., Sühring, M., von Tils, R., Vollmer, L., Ward, S., Witha, B., Wurps, H., Zeidler, J., and Raasch, S.: Overview of the PALM model system 6.0, Geosci. Model Dev., 13, 1335–1372,, 2020. a

Martilli, A., Clappier, A., and Rotach, M. W.: An Urban Surface Exchange Parameterisation for Mesoscale Models, Bound.-Lay. Meteorol., 104, 261–304,, 2002. a, b, c

Martilli, A., Santiago, J. L., and Salamanca, F.: On the representation of urban heterogeneities in mesoscale models, Environ. Fluid Mech., 15, 305–328,, 2015. a

McNorton, J. R., Arduini, G., Bousserez, N., Agustí‐Panareda, A., Balsamo, G., Boussetta, S., Choulga, M., Hadade, I., and Hogan, R. J.: An Urban Scheme for the ECMWF Integrated Forecasting System: Single‐Column and Global Offline Application, J. Adv. Model. Earth Sy., 13, e2020MS002375,, 2021. a

Nazarian, N. and Kleissl, J.: Urban Climate CFD simulation of an idealized urban environment: Thermal effects of geometrical characteristics and surface materials, Urban Clim., 12, 141–159,, 2015. a

Nazarian, N., Krayenhoff, E. S., and Martilli, A.: A one-dimensional model of turbulent flow through “urban” canopies (MLUCM v2.0): updates based on large-eddy simulation, Geosci. Model Dev., 13, 937–953,, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa

Oke, T. R., Mills, G., Christen, A., and Voogt, J. A.: Urban Climates, Cambridge University Press,, 2017. a, b, c, d, e

Piacsek, S. A. and Williams, G. P.: Conservation properties of convection difference schemes, J. Comput. Phys., 6, 392–405,, 1970. a

Poggi, D. and Katul, G. G.: The effect of canopy roughness density on the constitutive components of the dispersive stresses, Exp. Fluids, 45, 111–121,, 2008a. a

Poggi, D. and Katul, G. G.: Micro- and macro-dispersive fluxes in canopy flows, Acta Geophys., 56, 778–799,, 2008b. a, b, c

Raupach, M. R. and Shaw, R. H.: Averaging procedures for flow within vegetation canopies, Bound.-Lay. Meteorol., 22, 79–90,, 1982. a

Redon, E., Lemonsu, A., and Masson, V.: An urban trees parameterization for modeling microclimatic variables and thermal comfort conditions at street level with the Town Energy Balance model (TEB-SURFEX v8.0), Geosci. Model Dev., 13, 385–399,, 2020. a

Ross, A. N.: Boundary-layer flow within and above a forest canopy of variable density, Q. J. Roy. Meteor. Soc., 138, 1259–1272,, 2012. a

Salamanca, F., Krpo, A., and Martilli, A.: A new building energy model coupled with an urban canopy parameterization for urban climate simulations – part I. formulation, verification, and sensitivity analysis of the model, Theor. Appl. Climatol., 99, 331–344,, 2010. a

Santiago, J. L. and Martilli, A.: A Dynamic Urban Canopy Parameterization for Mesoscale Models Based on Computational Fluid Dynamics Reynolds-Averaged Navier-Stokes Microscale Simulations, Bound.-Lay. Meteorol., 137, 417–439,, 2010a. a

Santiago, J. L. and Martilli, A.: A Dynamic Urban Canopy Parameterization for Mesoscale Models Based on Computational Fluid Dynamics Reynolds-Averaged Navier–Stokes Microscale Simulations, Bound.-Lay. Meteorol., 137, 417–439,, 2010b. a, b, c, d, e

Santiago, J. L., Coceal, O., and Martilli, A.: How to Parametrize Urban-Canopy Drag to Reproduce Wind-Direction Effects Within the Canopy, Bound.-Lay. Meteorol., 149, 43–63,, 2013. a, b

Santiago, J. L., Krayenhoff, E. S., and Martilli, A.: Urban Climate Flow simulations for simplified urban configurations with microscale distributions of surface thermal forcing, Urban Clim., 9, 115–133,, 2014. a

Schmid, M. F., Lawrence, G. A., Parlange, M. B., and Giometto, M. G.: Volume Averaging for Urban Canopies, Bound.-Lay. Meteorol., 173, 349–372,, 2019. a

Schoetter, R., Kwok, Y. T., de Munck, C., Lau, K. K. L., Wong, W. K., and Masson, V.: Multi-layer coupling between SURFEX-TEB-v9.0 and Meso-NH-v5.3 for modelling the urban climate of high-rise cities, Geosci. Model Dev., 13, 5609–5643,, 2020. a

Siebesma, A. P., Soares, P. M. M., and Teixeira, J.: A Combined Eddy-Diffusivity Mass-Flux Approach for the Convective Boundary Layer, J. Atmos. Sci., 64, 1230–1248,, 2007. a, b, c

Simón-moral, A., Santiago, J. L., and Martilli, A.: Effects of Unstable Thermal Stratification on Vertical Fluxes of Heat and Momentum in Urban Areas, Bound.-Lay. Meteorol., 163, 103–121,, 2016. a, b, c, d

Sirko, W., Kashubin, S., Ritter, M., Annkah, A., Bouchareb, Y. S. E., Dauphin, Y., Keysers, D., Neumann, M., Cisse, M., and Quinn, J.: Continental-Scale Building Detection from High Resolution Satellite Imagery, arXiv [preprint], 2107.12283, 2021. a

Stiperski, I. and Calaf, M.: Dependence of near‐surface similarity scaling on the anisotropy of atmospheric turbulence, Q. J. Roy. Meteor. Soc., 144, 641–657,, 2018. a

Sun, J., Takle, E. S., and Acevedo, O. C.: Understanding Physical Processes Represented by the Monin–Obukhov Bulk Formula for Momentum Transfer, Bound.-Lay. Meteorol., 177, 69–95,, 2020. a

Sušelj, K., Teixeira, J., and Matheou, G.: Eddy Diffusivity/Mass Flux and Shallow Cumulus Boundary Layer: An Updraft PDF Multiple Mass Flux Scheme, J. Atmos. Sci., 69, 1513–1533,, 2012. a

Sützl, B. S., Rooney, G. G., Finnenkoetter, A., Bohnenstengel, S. I., Grimmond, S., and van Reeuwijk, M.: Distributed urban drag parametrization for sub‐kilometre scale numerical weather prediction, Q. J. Roy. Meteor. Soc., 147, 3940–3956,, 2021a. a, b

Sützl, B. S., Rooney, G. G., and van Reeuwijk, M.: Drag Distribution in Idealized Heterogeneous Urban Environments, Bound.-Lay. Meteorol., 178, 225–248,, 2021b. a

Williamson, J.: Low-storage Runge-Kutta schemes, J. Comput. Phys., 35, 48–56,, 1980. a

Wilson, J. D.: A second-order closure model for flow through vegetation, Bound.-Lay. Meteorol., 42, 371–392,, 1988. a

Xie, Z.-T., Coceal, O., and Castro, I. P.: Large-Eddy Simulation of Flows over Random Urban-like Obstacles, Bound.-Lay. Meteorol., 129, 1–23,, 2008. a

Ye, Q., Schrijer, F. F., and Scarano, F.: Geometry effect of isolated roughness on boundary layer transition investigated by tomographic PIV, International J. Heat Fluid Flow, 61, 31–44,, 2016.  a

Yuan, C., Shan, R., Zhang, Y., Li, X.-X., Yin, T., Hang, J., and Norford, L.: Multilayer urban canopy modelling and mapping for traffic pollutant dispersion at high density urban areas, Sci. Total Environ., 647, 255–267,, 2019. a

Short summary
This study enhances urban canopy models by refining key assumptions. Simulations for various urban scenarios indicate discrepancies in turbulent transport efficiency for flow properties. We propose two modifications that involve characterizing diffusion coefficients for momentum and turbulent kinetic energy separately and introducing a physics-based mass-flux term. These adjustments enhance the model's performance, offering more reliable temperature and surface flux estimates.