Budgets of Second-Order Turbulence Moments over a Real Urban Canopy

This study analyses budgets of second-order turbulence moments over a real urban canopy using large-eddy simulation. The urban canopy is representative of the City of Boston, MA, United States and is characterized by a significant height variability relative to the mean building height. The budgets of double-averaged Reynolds-stress components, scalar fluxes, and scalar variances are examined with a focus on the importance of the dispersive terms above the mean building height. Results reveal the importance of the wake (dispersive) production term, in addition to the shear production term, in the turbulence kinetic energy (TKE), streamwise velocity variance and scalar variance budgets well above the mean building height. In this region, the turbulent and dispersive transport terms are smaller than the production and dissipation terms. Nonetheless, the dispersive transport terms in the TKE and scalar variance budgets can be as important as their turbulent counterparts. The subgrid-scale dissipation term is the main sink in the TKE, vertical velocity variance and scalar variance budgets. In the momentum and scalar flux budgets, the pressure-strain correlation term and the pressure gradient-scalar interaction term are the significant sink terms, respectively. Our analysis highlights the complexity associated with the budgets of second-order turbulence moments over real urban canopies and has important implications for developing urban parameterizations for weather and climate models.

2016; Auvinen et al. 2017), the quantification of pedestrian thermal comfort Nazarian and Lee 2021), and the estimation of building energy consumption (Zhao and Magoulès 2012;Javanroodi et al. 2022). During the past decades, substantial efforts have been devoted to characterizing turbulence and turbulent transport inside the so-called urban roughness sublayer (RSL) (Rotach 1999;Masson 2006;Fernando et al. 2010). The established paradigm is that the urban RSL spans from the street level to roughly two to five times the mean building height (Oke et al. 2017); in this layer, turbulence and its associated transport are strongly influenced by the individual urban roughness elements (e.g., buildings, trees) and are thus both vertically and horizontally inhomogeneous (Rotach 1993;Oikawa and Meng 1995;Kastner-Klein et al. 2001;Britter and Hanna 2003).
The grid resolution in numerical weather prediction (NWP) models is too coarse to explicitly resolve individual roughness elements (Skamarock et al. 2008). To account for the impact of the roughness elements on the resolved-scale exchange processes between the urban canopies and the atmosphere, it is hence common to homogenize the governing equations at the horizontal grid resolution of NWP models (Chen et al. 2011). The horizontal grid resolution of NWP models is typically a few kilometers (the neighborhood scale), over which some degree of statistical homogeneity in the canopy morphology and resulting flow statistics might be expected (Britter and Hanna 2003). Given that the urban canopy domain is not simply connected (or it is multiply connected), the coarse-graining operation has to be based on the volume-averaging theorem (Whitaker 1967), whose theoretical and implementation details for flow over rough surfaces are discussed in Nikora et al. (2007), Mignot et al. (2009), Xie and Fuka (2018) and Schmid et al. (2019). When time-and spatialaveraging (hereafter double-averaging) operations are performed in a multiply-connected domain, additional terms arise in the averaged equations besides the turbulent fluxes, namely the dispersive fluxes (Mahrt 1987). The turbulent fluxes are caused by temporal deviations from the temporally-averaged flow, while the dispersive fluxes arise from the spatial correlations of temporally-averaged flow quantities over the spatial averaging scale. The dispersive fluxes and related dispersive terms in the budget equations for turbulence moments remain poorly understood and their contributions to the flow dynamics are often overlooked.
Previous studies of turbulence moment budgets in flow over urban canopies have primarily focused on the turbulence kinetic energy (TKE) budget (Louka et al. 2000;Bou-Zeid et al. 2009;Christen et al. 2009;Giometto et al. 2016;Blackman et al. 2017;Tian et al. 2021;Blunn et al. 2022). Much less is known about the budgets of momentum and scalar fluxes, as well as scalar variances over urban canopies. This is in sharp contrast to the attention these budgets received in studies over vegetative canopies (Meyers and Baldocchi 1991;Dwyer et al. 1997;Katul et al. 2009Katul et al. , 2013Viana Parente Lopes et al. 2021;Watanabe et al. 2021). Moreover, motivated by the need to understand the physical system in its simplest form, the majority of previous studies have considered idealized urban canopy configurations such as arrays of aligned or staggered cuboids Yakhot et al. 2006;Blackman et al. 2017;Tian et al. 2021); these canopies are characterized by a few length scales and hence lend themselves to analytical treatment. As shown in recent work, the dynamics of turbulent transport in idealized conditions might profoundly differ from their real-world counterparts (Giometto et al. 2016;Inagaki et al. 2017;Auvinen et al. 2020;Akinlabi et al. 2022). To bridge this knowledge gap, we conduct a budget study over a real urban canopy. Specifically, we focus on quantifying the double-averaged budgets for second-order turbulence moments. We propose to use large-eddy simulations (LESs) because field studies of turbulence budgets have been restricted to one or few locations (Rotach 1993;Christen et al. 2009;Santiago and Martilli 2010). The LES technique has been applied to study the budgets of second-order turbulence moments over vegetation (Dwyer et al. 1997;Yue et al. 2008;Viana Parente Lopes et al. 2021) and urban canopies (Bou-Zeid et al. 2009;Giometto et al. 2016;Tian et al. 2021). However, compared to previous work, our contribution is novel because we examine for the first time the momentum and scalar fluxes, as well as scalar variance budgets over a real urban environment.
In what follows, a standard notation is used where x i = (x, y, z) are the Cartesian coordinates (i.e., x, y, z represent the streamwise, spanwise and vertical directions, respectively), and u, v, w are the streamwise, spanwise and vertical velocity components (resolved by LES), respectively; and s represents a passive scalar (e.g., the concentration of pollutants etc.). The Einstein summation convention for repeated indices is used. The overbar (.) and angular brackets . denote time and spatial (volume) averaging, respectively. Double-averaging (DA) refers to taking the average in time first and then in space. The prime and double prime denote temporal and spatial deviations, respectively. Namely, X = X − X is the temporal fluctuation of X (i.e., deviations from the temporally-averaged X ) and X = X − X is the spatial deviation of X from its spatial average X . This paper is organized as follows: Sect. 2 provides the theoretical framework and presents the double-averaged budgets of second-order moments in a multiply-connected domain; the large-eddy simulation model and the simulated case are presented in Sect. 3; Sect. 4 presents the analysis of second-order moment budgets and conclusions are drawn in Sect. 5.

Volume Averaging
The volume averaging operation is carried out over horizontal slab of thickness z. Two types of volume averaging need to be distinguished. The first is intrinsic averaging (Nikora et al. 2007), where the averaging volume includes the ambient air only. The second is extrinsic (or superficial) averaging (Schmid et al. 2019), where averaging is performed over the entire horizontal slab (i.e., including the volume occupied by solid elements such as buildings and trees). Intrinsic averaging is widely used in the literature to characterize flow fields over vegetation canopies (Wilson and Shaw 1977;Raupach and Shaw 1982), gravel beds (Nikora et al. 2007), rigid canopies (Raupach et al. 1991;Coceal et al. 2006;Xie et al. 2008) and real urban canopies (Giometto et al. 2016;Akinlabi et al. 2022). The intrinsic averaging operation is the natural approach for this study for two reasons. The first is that the resulting statistics are more representative of typical values inside the fluid. The second is that intrinsically-averaged dispersive fluxes are zero for a constant velocity field (due to the zero spatial deviation of the constant velocity field from its mean) but the superficially-averaged dispersive fluxes are not necessarily zero.
For a volume (V ) centered at location x i that composes of fluids (with volume V f ) and solid elements (with volume V s ), we define the intrinsic average of a temporally averaged variable F as: (1)

Double-Averaged Budgets of Second-Order Turbulence Moments
Double-averaged budgets are obtained by first averaging the flow field in time and then in space using the intrinsic volume-averaging operations on temporally averaged fields. The time-averaged budget equations have been extensively studied and can be found in classical textbooks (Stull 1988;Garratt 1992). Therefore, in the following, we only briefly discuss the volume-averaging rules for flow in urban canopies (Schmid et al. 2019). We define a timeand intrinsically-averaged quantity ϕ , where the intrinsic and time-averaging operations commute. In the following, the quantity ϕ is termed the double-averaged ϕ. However, based on the volume averaging theorem (Whitaker 1967), the intrinsically-averaged spatial gradient of ϕ is not equal to the spatial gradient of ϕ , but rather: where α p = 1 − λ p and λ p is the plan area fraction defined as the fraction of space occupied by the solid elements in a given averaging volume. The function α p is needed to account for the change of the fluid volume V f with height, which is important for real urban canopies (Giometto et al. 2016). For generality, α p is written as a function of x i in Eq. 2 (and Eq. 5 below) following Schmid et al. (2019). The surface integral represents the effect of the solid-fluid interface and is zero when ϕ is any of the velocity components due to the no-slip boundary conditions. A f s (x) is the solid-fluid interface contained in the averaging volume V (x) while n i is the unit normal vector of A f s pointing from the fluid phase into the solid phrase. For more details about the surface integral in Eq. 2, readers are referred to Mignot et al. (2009) and Schmid et al. (2019). With these rules, the budget equations for double-averaged second-order moments can be obtained using the following procedure. First, we average each term in time and space. Second, we switch the order of spatial averaging and differentiation following Eq. 2. Third, we expand the spatial averaging of the product of ϕ i and ϕ j following: Note that ϕ = 0 due to the averaging rules. The spatial averaging of the product of ϕ i and the gradient of ϕ j reads: In this case, ∂ϕ j ∂ x i does not disappear, as can be seen from Eq. 2, namely: When applying the double averaging procedure to analyzing LES outputs, we make further simplification by assuming horizontal homogeneity at scales beyond the spatial averaging scale (i.e., ∂ ϕ ∂ x = ∂ ϕ ∂ y = 0), stationarity (i.e., ∂ ϕ ∂t = 0), and no large-scale subsidence, (i.e., w = 0). Furthermore, because of the assumption of horizontal homogeneity at scales larger than the spatial averaging scale, α p becomes also only a function of z and thus only ∂α p ∂z = dα p dz is non-zero. In Appendix 1, we follow the above-mentioned procedure to derive the budget equations for double-averaged second-order moments.

Reynolds Stress
The budget equation for double-averaged Reynolds stress tensor reads: where P s ik is the shear production term, P w ik is the wake (dispersive) production term, P m ik is the rate of work of the temporally averaged velocity fluctuations against the shear production (given that α p varies with height, P m ik = 0), which is called the form induced production term hereafter, T t ik is the turbulent transport term, T d ik is the dispersive transport term, T p ik is the pressure transport term, S P ik is the pressure-strain correlation or pressure redistribution term, D ik is the SGS transport term and ε ik is the SGS dissipation term. The SGS thirdorder velocity correlations and SGS pressure-strain correlation are negligible above the mean building height; they will be incorporated into the budget residual term. In general, Eq. 6 shows that each component of the Reynold stress is produced by shear production (P s ik ), wake production (P w ik ) and form induced production (P m ik ), transported by turbulent transport (T t ik ), dispersive transport (T d ik ), pressure transport (T p ik ) and SGS transport (D ik ), and dissipated by ε ik . It is assumed that ε ik takes place at small scales where the local isotropy of the Kolmogorov hypothesis prevails. As a result, the Reynold stress dissipation is modeled as an isotropic tensor and its deviatoric part is incorporated into the pressure-strain correlation term (making this term a sink for this case). The isotropic part of ε ik is computed using the SGS model. This procedure is standard in second-order turbulence closure modelling .

TKE
The budget equation for TKE is simply the trace of Eq. 6 multiplied by 1/2. It reads: where the definition for P s T K E and D T K E are similar to the definitions given in Eq. 6 and ε T K E is the TKE dissipation. The pressure-strain correlation is not in the TKE budget because it only acts to redistribute energy between the components. In general, Eq. 7 shows that TKE is produced by shear production P s T K E , wake production P w T K E and P m T K E , redistributed by turbulent transport T t T K E , dispersive transport T d T K E , pressure transport T p T K E and SGS transport D T K E and finally dissipated by the work of SGS stresses onto the resolved field ε T K E (Christen et al. 2009;Giometto et al. 2016).

Vertical Scalar Flux
The budget equation for the vertical scalar flux reads: where the definitions for P s SF , P w SF , P m SF , T t SF and T d SF are similar to the definitions in Eq. 6. P S is the pressure gradient-scalar interaction (a de-correlation term) and D SF include SGS terms; ψ SF represents the surface integral terms arising from the volume averaging theorem. We point out that in deriving Eq. 6, the pressure term in the time-averaged Reynold stress tensor u i u k budget equation is split into two terms (see Appendix 1 for details). On the contrary, for scalar flux budget equation, this pressure term cannot be split. The SGS components of P S have been shown to be non-negligible in the scalar flux budget (Khanna 1998;Heinze et al. 2015). Due to difficulties in evaluating the SGS components of P S and the surface integral terms over complex urban geometry via the LES model used in our study (to be introduced later), this work incorporates them into the budget residual along with the SGS component of the turbulent transport term.

Scalar Variance
The budget equation for the scalar variance reads: where the definitions for P s SV , P w SV , P m SV , T t SV and T d SV for the scalar variance are similar to those in Eq. 6. D SV is the SGS transport for scalar, ε SV is the dissipation and ψ SV is the surface integral terms that capture the effect of the solid-fluid interface. ψ SV is incorporated into the budget residual due to the difficulties associated with evaluating it in the LES model used in our study.

Large-Eddy Simulation Model
In this study, the PALM LES Model in revision 4901 (Maronga et al. 2015(Maronga et al. , 2020 is used. The PALM solver numerically integrates the filtered, non-hydrostatic Navier-Stokes equations in the Boussinesq-approximations form and the filtered transport equation for passive scalar concentration. Filtered transport equations for two thermodynamic variables, such as potential temperature and total water specific humidity, can be solved but are not used in this study. The filtering of these equations is carried out implicitly using the volume-based approach (Schumann 1975) and employs the 1.5-order SGS closure model of Deardorff (1980). Using a predictor-corrector method and iterative multigrid scheme (Hackbusch 1985), the mass conservation of the flow is enforced by solving a Poisson equation for pressure perturbation. The 5th order Wicker-Skamarock and the 2nd order central difference schemes are employed to discretize the advection and diffusion schemes. Temporal discretization is done with the 3rd order Runge-Kutta scheme. The computational domain is spatially discretized using the finite difference approach on Arakawa staggered C-grid (Arakawa and Lamb 1977). The PALM model explicitly resolves the solid obstacles using the masking method (Briscolini and Santangelo 1989) and hence does not need a parameterization to account for the effect of the solid obstacle on the flow dynamics.
The PALM model has been widely used to study flows over both idealized (Letzel et al. 2008;Park et al. 2012;Gronemeier and Sühring 2019;Nazarian et al. 2020;Blunn et al. 2022) and real urban canopies (Kanda et al. 2013;Park et al. 2015;Gronemeier et al. 2017) and it has been extensively validated (Fröhlich and Matzarakis 2020;Gronemeier et al. 2021;Resler et al. 2021). Heinze et al. (2015) used the PALM model to study second-order moment budgets in cloud topped boundary layers and found that the PALM model results agree with the results of other LES models except for the TKE dissipation rate. The disagreement in the TKE dissipation rate was attributed to truncation errors, which can be relatively large and lead to artificial dispersion (especially at high wavenumbers) when using low order schemes (Ghosal 1996;Giacomini and Giometto 2021). Uncertainties arising from these errors are captured in the residual of our budget analysis and discussed in the result section.

Case Description and Model Set-Up
In this study, we focus on an area of about 2.6 × 2.1 km 2 around Fenway-Kenmore square in the City of Boston, Massachusetts, USA (see Fig. 1a). This geographical region is the same as the one considered in Akinlabi et al. (2022). This region is located in the northern part of Boston. The chosen domain contains a dense arrangement of building blocks, an irregular distribution of narrow street canyons, a park in the northwest region, and the Charles River in the north. The northeastern part is a business district with many high-rise buildings of height above 100 m (e.g., the Prudential centre which is 227 m high), while the southwestern part is the home to several hospitals (Boston children hospital, Beth Israel Medical centre, Brigham and Women's hospital) and universities (Harvard school of public health, Emmanuel college, Simmons university and Massachusetts College of Pharmacy and Health Sciences) with moderately tall buildings of about 60-80 m. Figure 1 shows the map, the vertical profile of the plan area fraction, and the building height distribution, and its probability density function within the study area. The mean building height H is 18 m, and the standard deviation σ H is 16 m. The plan area fraction varies strongly with height and is 0.29 at the Fig. 1 a 3-D map of the area around Fenway-Kenmore square in the City of Boston, USA, with the Charles River in the north, the Brigham and Women's hospital (about 75 m) in the southwest, and the Prudential centre (maximum building height of 227 m) in the northeast. Imagery ©Google, b the vertical profile of the plan area fraction, c building heights in the study area and d the probability distribution function (PDF) of building heights ground level. The distribution of building heights in our study area differs from previous studies like Auvinen et al. (2020), whose building height distribution is relatively symmetric with σ H /H = 0.4 − 0.6. It is also different from the study by Giometto et al. (2016Giometto et al. ( , 2017 with a trimodal building height distribution and σ H /H = 0.4. In our study, the distribution is very skewed with σ H /H = 0.89. Vegetation is not included in our simulation, which is justified by its small plan area fraction (Giometto et al. 2016).
The domain is discretized in space using 864 × 720 × 360 grid points in streamwise, spanwise and vertical directions, respectively. A horizontal grid spacing of 3 m is used, which has been shown in our previous work (Akinlabi et al. 2022) to be adequate in resolving the buildings and street canyons in the domain. In the vertical direction, 3 m grid spacing is used up to 300 m. Above this height, we apply a grid stretching with a factor of 1.005 until a maximum value of 11 m is reached. This gives a domain height of 1.9 km. The boundarylayer height is δ/H = 70, which satisfies the δ/H 50 requirement (Jimenez 2004).
The flow is driven by a constant geostrophic wind U g = 3.5 m s −1 (an intermediate value to represent a weakly sheared flow) in the west-to-east direction and neutral stratification is assumed throughout the study. The no-slip wall boundary condition is imposed on all surfaces (including the roofs, ground, and building walls) whereas a free-slip condition is applied at the top of the domain. We apply an algebraic wall-layer model between the surface (including the roofs, ground, and building walls) and the first computational grid level. To account for the effects of low vegetation, structural details, and temporary structures, a SGS aerodynamic roughness length z 0,SGS = 0.01 m is used, which follows the recommendation of Basu and Lacser (2017) that z 0,SGS ≤ 0.02 × min( z). The value min( z) for this study is 1.5 m because the first computational grid node is positioned at 0.5 z. Cyclic boundary conditions are imposed in the lateral directions to simulate an infinite repetition of the study area. This setup is convenient as it does not require specification of an inflow boundary condition. The boundary condition for the passive scalar equation is a surface flux 0.05 kg m −2 s −1 , which is imposed on all surfaces (including the roofs, ground and building walls). The simulation runs for a spin-up period of 560T where T = H /u * to reach a steady state. Here T is interpreted as the eddy-turnover time for the largest eddies in the urban canopy (Coceal et al. 2006). The friction velocity u * = 0.3 m s −1 is computed from the total kinematic surface drag per unit floor area τ * , i.e., u * = √ τ * /ρ, where ρ is the air density (1 kg m −3 ) and τ * is the sum of the form and skin-friction drag (Kanda et al. 2013). The simulation is then pursued for another 240T to evaluate temporally averaged statistics, which has been verified to be long enough for the statistics to reach convergence (Akinlabi et al. 2022).

Double-Averaged Flow Statistics
We start our analysis by examining the double-averaged flow statistics. The streamwise velocity, vertical velocity, momentum flux, velocity variances, and the total pressure drag are normalized by the friction velocity u * . The scalar concentration is normalized with s * = w s 0 /u * where w s 0 = 0.05 kg m −2 s −1 is the surface scalar flux. The turbulent scalar flux is normalized with u * s * . The vertical height is normalized with the mean building height (H = 18 m). Figure 2 shows the normalized profiles of double-averaged streamwise and vertical velocities and their variances, TKE, and the logarithm of the scalar concentration and its variance. The streamwise velocity profile exhibits no inflection point, consistent with profiles presented by previous real urban canopy studies with large σ H values (Park et al. 2015;Inagaki et al. 2017;Akinlabi et al. 2022). The reason for this, as discussed by Makedonas et al. (2021) and Akinlabi et al. (2022), is the large spread of velocities below the height H max , indicating significant flow penetration caused by the large σ H . As a result, cities designed with large σ H could have higher mixing rates, which can positively impact urban air quality and natural ventilation (Makedonas et al. 2021). The profile of the streamwise velocity follows the conventional logarithmic form well above the urban canopy. However, closer to the buildings, the streamwise velocity profile deviates from the logarithmic form as it responds directly to the urban canopy (see   Kanda et al. (2013), while the displacement height is overestimated (i.e., z d /H = 2.22 using the equation in Kanda et al. 2013). The reason for the difference in displacement height estimates is beyond the scope of this work. Following our earlier work (Akinlabi et al. 2022), we identify z/H = 30 height as the RSL thickness (similar to the fitting range used above). This height corresponds to the 90th percentile of the dispersive flux profile. Dispersive fluxes for the urban canopy under consideration were examined in detail in Akinlabi et al. (2022) and will not be discussed here.
The double-averaged vertical velocity vanishes as expected from the use of periodic lateral boundary condition. The normalized TKE has its maximum value of 3 around z/H = 5 and decreases with increasing height with major contribution from the streamwise velocity variance (see Fig. 2b). The streamwise velocity variance peaks at z/H = 10 while the vertical velocity variance peaks at z/H = 3. The profile of the logarithm of scalar concentration is almost uniform with height, even though almost all the source of the scalar concentration is below z/H = 2 based on the boundary condition for the passive scalar. This uniformity indicates an intense mixing of passive scalar from urban surfaces where it is released to the atmosphere. This vigorous mixing may be caused by the significant flow penetration discussed above. Here, we show the profile of the logarithm of the scalar and its variance to highlight their variations better. The logarithm of normalized scalar variance has a maximum value of 6 at z/H = 1.
The turbulent momentum and scalar fluxes as well as the pressure drag are presented in Fig. 3. Only the resolved parts of the turbulent fluxes are presented since the subgrid-scale fluxes are less than 6% of the sum of resolved and subgrid-scale fluxes above z/H = 1. The turbulent momentum flux peaks at about z/H = 10 with magnitude twice as large as its value at z/H = 40. A similar behavior is observed for the turbulent scalar flux, which peaks at z/H = 4 (see Fig. 3b). The pressure drag, which is the major sink of momentum in the urban canopy, decreases with height from its surface value Fig. 3a).

Budgets for Second-Order Turbulence Moments
The budgets for second-order moments in the urban RSL are now discussed. Each term in the TKE, velocity variances and Reynold stress budgets are normalized by H /u 3 * . Those for scalar flux and scalar variance are normalized by H /(s * u 2 * ) and H /(u * s 2 * ), respectively. According to Akinlabi et al. (2022), the urban RSL for real urban canopies can extend much higher than the traditional definition (i.e., z/H = 2 − 5) which is primarily based on studies over idealized urban canopies. Using the height that corresponds to 90th percentile of the dispersive flux profile as the beginning of the inertial sublayer, they argued that the RSL extends to z/H = 30. Following Akinlabi et al. (2022), we will focus on the interval z < 30H . Three layers are considered. The first layer is the traditionally defined urban RSL (2H ≤ z ≤ 5H ) represented by the grey area in Figs. 4,5,6,7,8 and 9. We note that this layer roughly covers the region where the lowest atmospheric grid (about 30-100 m) in NWP and climate models with single-layer urban parameterizations often occurs. Above this interval, two additional layers are considered: 5H < z ≤ 12H (the second layer) and 12H < z ≤ 30H (the third layer). The second layer spans from the top of the traditionally defined urban RSL to H max , while the third layer are from H max to the top of the urban RSL. Averages of each budget term within each layer are presented in Tables 1, 2, 3, 4, 5 and 6.
The relatively small residual R in the computed budget of second-order moments when z is above H provides confidence in our numerical results. The residual terms contain all other SGS components of the budget terms such as the SGS third-order velocity correlations and the SGS pressure redistribution (see Heinze et al. 2015 for example). The residual below H (see Figs. 4,5,6,7,8,9) is primarily due to the spatial interpolation of variables in the near wall regions required to compute some of the budget terms; this leads to numerical truncation errors and degrades the quality of the computed budget. Hence, only the budgets at z/H > 1 will be analyzed.

TKE
Vertical profiles of terms in the TKE budget are shown in Fig. 4 while Table 1 shows the percentage contribution of each term to the total source (+) or sink (−) in the considered layer. The shear production P s T K E peaks at z/H = 1 where the strongest wind shear occurs and decreases with height. This agrees with previous studies of boundary layer flows over uniform strip or tree-like canopies (see Yue et al. 2008 andBöhm et al. 2013). Although P s T K E decreases with height for z/H > 1, its contribution to the total source increases with height because other production terms become even smaller with height (see Table 1). The wake production P w T K E is the production rate of TKE in the wakes of buildings (i.e., converting wake kinetic energy to TKE) through the interaction between the local turbulent stress and time-averaged strains. P w T K E also peaks at H and decreases to approximately zero above z/H = 15. Below H max , P w T K E ≈ 0.5 P s T K E , in agreement with previous studies of flow over real urban canopies (Giometto et al. 2016). This implies that P w T K E is nonnegligible over the urban canopy. Similar results have been presented in studies of flow over other regular canopies (Raupach et al. 1991). The form-induced production term P m T K E is negligible in our study (see Table 1). This result disagrees with Giometto et al. (2016), where P m T K E is found to be non-zero in the vicinity of the inflection layer, accounting for 16% of P s T K E . This disagreement is likely caused by the difference in σ H . The difference between Fig. 5 The streamwise velocity variance budget terms normalized by H /u 3 * . The grey region corresponds to 2H ≤ z ≤ 5H while the dashed horizontal line is the maximum building height H max the two studies seems to suggest that the importance of P m T K E decreases or even becomes negligible with large σ H . Note that other factors such as the plan area fraction λ p , the frontal area fraction λ f = A f /A total (A f is the product of the building width and height) might also be responsible for this disagreement. More detailed investigations of how P m T K E (as well as similar terms in the budgets of other second-order moments) respond to changes in the aforementioned parameters is beyond the scope of this analysis and is left for future work.
The transport terms are responsible for redistributing TKE vertically from regions of high production to others. They serve as local sources/sinks of TKE (Roth and Oke 1993). Within 2H ≤ z < 12H , the turbulent transport term T t TKE is negative and contributes to 10% of the total sink of TKE (see Table 1). It changes sign at z/H = 2 and z/H = 13, contributing 5% of the total source above z/H = 15. Our result agrees with studies of flow over urban canopies (Christen et al. 2009;Giometto et al. 2016) and field studies of flow over vegetation canopies (Leclerc et al. 1990;Shen and Leclerc 1997). T d T K E , T p T K E and D T K E are almost zero in the studied height ranges (see Table 1). The result of D T K E agrees with Yue et al. (2008).
The TKE dissipation rate ε T K E is a significant sink of TKE. We compute ε T K E as τ i j S i j = τ i j S i j − τ i j S i j where S i j is the filtered shear rate tensor while τ i j is the SGS stress tensor. Experimental studies compute ε T K E based on the energy spectra (e.g., Christen et al. 2009), but this approach is known to overestimate ε T K E Akinlabi et al. 2019). Here we found that ε T K E has a maximum value of 1.4 H /u 3 * at z/H = 1 (though it might increase even more within z/H < 1) and decreases with height until it balances TKE production at z/H > 20, after which P s T K E ≈ ε T K E . Local contributions of ε T K E to the total sink rate of TKE range between 86% at 5H ≤ z < 12H to 93% at z/H > 12.
Based on these results, we conclude that for the real urban canopy studied here, the shear production P s T K E , wake production P w T K E and dissipation of TKE ε T K E are the major players in the TKE budget. They need to be parameterized in large-scale meteorological models due to their significant contributions to the total source or sink of TKE in real urban canopy flows. The contributions of turbulent transport T t T K E to local TKE sources/sink are less than 15% with significant height variability.

Velocity Variances
In this section, we further examine the budgets of streamwise and vertical velocity variances. Tables 2 and 3 show the percentage contribution of each term to the total source (+) or sink (−) in the considered layers and the profiles are shown in Figs. 5 and 6 respectively. P s 11 and P w 11 are the key source terms in the budget of streamwise velocity variance, with peak values of 2.5H /u 3 * and 0.8H /u 3 * , respectively, at z/H = 1. The profiles of P s 11 and P w 11 are similar in shape to P s T K E and P w T K E , respectively. Since P s 33 = 0, the production of TKE due to Fig. 7 The momentum flux budget terms normalized by H /u 3 * . The grey region corresponds to 2H ≤ z ≤ 5H while the dashed horizontal line is the maximum building height H max shear occurs through the horizontal velocity components (i.e., P s T K E ≈ 0.5(P s 11 + P s 22 )). This explains why the profile of P s 11 is similar to that of P s T K E . Here it should be pointed out that these results might be altered by thermal stratification which is absent in our study. P w 33 is non-zero but contributes less than 10% to the total source rate of vertical velocity variance (see Table 3), explaining why the profile of P w 11 is similar to that of P w T K E . The form-induced production terms P m 11 and P m 33 are negligible. The anisotropy introduced by shear and wake productions is counteracted by the pressurestrain correlation terms S P 11 and S P 33 . S P 11 and S P 33 only act to redistribute TKE between the components returning turbulence to the isotropic state-a process known as "isotropization of turbulence" (Pope 2000;Hanjalić and Launder 2009). S P 11 is negative while S P 33 is positive throughout the considered height intervals. This implies that the vertical velocity variance grows at the expense of the streamwise velocity variance. The dissipation rates (ε 11 and ε 33 ) in the velocity variance budgets are determined based on the assumption of local isotropy at small scales i.e., ε 11 ≈ ε 33 ≈ 2 3 ε T K E . The percentage contribution of ε 11 to the total sink in the streamwise velocity variance budget is about 50% of the percentage contribution of S P 11 , with the sum of S P 11 and ε 11 nearly balancing the production terms. In the vertical velocity variance budget, S P 33 nearly balances ε 33 .
All the transport terms (i.e., T t 11 , T t 33 , T d 11 , T d 33 , T p 11 , T p 33 , D 11 and D 33 ) are much less critical in the velocity variance budgets. T t 11 and T t 33 make about 5-10% contribution to the total sink while the other transport terms are even smaller compared to other terms in the velocity  Tables 2, 3). In summary, the production, pressure-strain correlation and dissipation terms play significant roles in the velocity variance budgets.

Momentum Flux
Before we present the momentum flux budget, it is important to make a remark related to its interpretation. Unlike the velocity and scalar variances that are non-negative, the momentum flux can have either sign. As a result, any term in the momentum flux budget is treated a source term if it has the same sign as the momentum flux itself and a sink term if it has the opposite sign. To avoid any confusion, we multiply both side of the budget equation for double averaged Reynolds stress tensor with a negative sign so that a negative and positive term is a sink and source term, respectively.
Vertical profiles of terms in the budget of momentum flux are shown in Fig. 7 while the layer-wise percentage contribution of each term to the total layer source or sink are presented in Table 4. Approximate equilibrium exists between P s 13 and the pressure redistribution term S P 13 above z/H = 15, which agrees with Raupach et al. (1986). For all the considered layers, T t 13 has a sink contribution of around 4-11%, i.e., it is larger than other transport terms. However, the same term becomes a source below 2H . The significance of T t 13 over rough surfaces is not a new finding and has been reported by Maitani (1979) and Raupach (1981). Here we simply note that the momentum flux budget over real urban canopies has not been analyzed thus far. The closest comparison is the momentum flux budget for plant canopies based on measurements from Meyers and Baldocchi (1991). Our findings agree Fig. 9 The scalar variance budget terms normalized by H /(u * s 2 * ). The grey region corresponds to 2H ≤ z ≤ 5H while the dashed horizontal line is the maximum building height H max Table 1 Percentage contribution of P s T K E , D T K E and ε T K E to the total source and sink for the considered layers (+) and (−) denote a source and sink of TKE, respectively with Meyers and Baldocchi (1991) and Raupach et al. (1986) regarding the dominant role of the shear production term P s 13 and the pressure-strain correlation term S P 13 above the canopy.
(+) and (−) denote a source and sink of streamwise velocity variances, respectively (+) and (−) denote a source and sink of scalar fluxes, respectively Table 6 Percentage contribution of P s SV , P w SV , P m SV , T t SV , T d SV , D SV and ε SV to the total source and sink for the considered layers (+) and (−) denote a source and sink of scalar variances, respectively

Scalar Flux
The results for the scalar flux budget are shown in Fig. 8 and Table 5. Compared to the TKE, velocity variance, and momentum flux budgets, the scalar flux budget and the scalar variance budget to be discussed in the following section still have relatively large residuals at the lower heights since the SGS components of P S and the surface integral terms are incorporated into the budget residual. The residuals gradually decrease with height and become zero around 2H and 4H in the scalar flux budget and the scalar variance budget, respectively. Hence, the results below 2H and 4H for the scalar flux budget and the scalar variance budget, respectively, should be interpreted with caution. For the scalar flux budget, the terms P s SF , P w SF and P S have their extrema near the surface, where large gradients of s occur. P s SF and P S terms dominate the budget at z/H > 2 while P w SF is also important at z/H < 2, with the caveat that the residual remains large for z/H < 2. P s SF is positive (since s is a decreasing function of height) above z/H = 1 and is nearly balanced by the pressure gradient-scalar interaction P S, which acts to destroy scalar flux. T t SF is negative within the range 2H ≤ z ≤ 15, similar to the turbulent transport term for momentum flux. Above 12H , T t SF may be neglected since its contribution to the budget is only 5%. All other terms are rather small (less than 5% contribution) (see Table 5).
Similar to the momentum flux budget, the scalar flux budget over real urban canopies has not been analyzed so far. Hence, a direct comparison of our findings with previous results is not possible. Profiles of P s SF and T t SF agree with those in Coppin et al. (1986), in which scalars were emitted within a plant canopy in a wind-tunnel. Unfortunately, all other terms were not computed in their study. Our finding regarding the dominance of P s SF and P S also agrees with the simplified analysis of the scalar flux budget in the inertial sublayer (Garratt 1992).

Scalar Variance
In the budget of scalar variance, P s SV decreases with height from its peak near the surface. P w SV also gradually decreases from its peak value near the surface and becomes negligible at z/H = 6 (see Fig. 9). The form-induced production term P m SV is generally small. Hence, the production term P s SV is the major source term above z/H = 6, which is balanced by the scalar dissipation ε SV (see Table 6). ε SV is estimated similarly as the TKE dissipation ∂s ∂ x j where τ s, j is the SGS scalar flux, computed by the SGS model. P s SV and ε SV are dominant terms in the scalar variance budget, in agreement with findings from Coppin et al. (1986). Other terms are minor except the transport terms T t SV (below 4H ) and T d SV and D SV (below 2H ). However, we stress again that the residuals below 4H are significant and thus the results below 4H should be interpreted with caution.

Relative Importance of the Dispersive Terms to the Reynolds Terms
As discussed in the introduction, due to difficulties in their measurement and simulation, dispersive terms such as wake production and dispersive transport have received less attention than their Reynolds counterparts. Results in Sect. 4.2 indicate that these dispersive budget terms may be important, especially within the first layer (2 ≤ z/H ≤ 5).
In this section, we contrast the dispersive and Reynolds budget terms by examining the ratio of their absolute values for the entire urban canopy considered in our model domain. These are labelled as "Ref" in Figs. 11 and 12. The ratios have been averaged over the considered layers. We do not present this ratio for the streamwise and vertical velocity variances since P s 33 is zero and this ratio for the streamwise velocity variance is similar to that of TKE. The symbol ϑ v in Figs. 11 and 12 indicates that the profile ϑ is averaged within the given height range. Figure 11 shows the relative importance of wake production terms. For TKE, the ratio of wake production term to shear production term decreases from 0.5 at 2 ≤ z/H ≤ 5 to approximately zero at z/H = 15. For scalar variance, the ratio of wake production term to shear production term also decreases with height, from a value of 0.4 at 2 ≤ z/H ≤ 5 and approximately zero at z/H > 5. The ratios of wake production term to shear production term for momentum and scalar fluxes exhibit similar profiles: they increase from 2 < z/H ≤ 5 to their peak values at 5 < z/H ≤ 12 and then decrease with height. Peak values for momentum and scalar fluxes are however relatively small (about 0.15).
The ratio of dispersive to turbulent transport terms is presented in Fig. 12 for the considered budget equations. Even though the magnitude of transport terms is small relative to the production terms in general (see Sect. 4.2), the dispersive transport terms can be significant relative to their turbulent counterparts. The relative importance of dispersive transport of TKE increases from about 0.1 at 2 < z/H ≤ 5 to over 1 at 12 < z/H ≤ 15 and then decreases with height. The ratios T d 13 / T t 13 v and T d SF / T t SF v are less than 0.2 throughout the studied height ranges. For the scalar variance, the ratio of dispersive transport term to the turbulent transport decreases monotonically with height from the peak value of 1.7 at 2 < z/H ≤ 5.
In summary, the dispersive terms are more critical in the TKE and scalar variance budgets than in the flux budgets. Their ratios to the corresponding Reynolds terms can be about 0.5 to 1 in the TKE and scalar variance budgets. The next step is to determine the sensitivity of the relative importance of dispersive terms to different urban geometric parameters. To do this, we partition our model domain (see Fig. 1a) into four subdomains in the y-direction. The area for each subdomain is about 2.6 × 0.5 km 2 . For this part of the analysis, the intrinsic spatial averaging is carried out over each subdomain and hence the condition ϕ = 0 is satisfied. The building height distribution, the PDF of building heights and the plan area fraction in each subdomain is presented in Fig. 10. The ratios of the standard deviations of building height to the mean building heights (σ H /H ) are greater than 1 for subdomains 3 and 4 (subdomain 3 = 1.31, subdomain 4 = 1.04) but less than 1 for the subdomains 1 and 2 (subdomain 1 = 0.72, subdomain 2 = 0.65). Figure 11 shows the ratios of wake productions to shear productions for TKE, momentum flux, scalar flux, and scalar variance in the subdomains. The importance of P w T K E decreases with height in all subdomains. The ratio is larger than 1 for subdomain 4 for the layer 2 ≤ z/H ≤ 5 (see Fig. 11a). For the momentum/scalar fluxes, the ratios remain less than 0.15 for all subdomains. Still, subdomain 4 has the most significant values (see Fig. 11b, c). For the scalar variance, the importance of wake production decreases with height for all subdomains. Only in subdomain 4 is the ratio greater than 0.5 for the layer 2 ≤ z/H ≤ 5 (see Fig. 11d). All in all, these results suggest that the wake production, especially for TKE and scalar variance, can become significant in the vicinity of tall buildings, as in subdomain 4. The enhanced importance of wake production in subdomain 4 suggests that the wake production may depend on H max (or the ratio of H max and H ).
The ratios of dispersive transport to turbulent transport terms show a much wide range of variabilities and do not exhibit any generalizable behaviors across the 4 subdomains. For TKE, the ratio has the most significant value of 0.5 in subdomain 4 at the layer 2 ≤ z/H ≤ 5 (see Fig. 12a). For momentum flux, the ratio has the largest value of 6 in subdomain 4 in 12 < z/H ≤ 30. For scalar flux, the maximum value of the ratio is about 8 and again occurs in subdomain 4 at the layer 2 ≤ z/H ≤ 5. However, for scalar variance, the most considerable value of the ratio occurs in subdomain 2 at the layer 2 ≤ z/H ≤ 5 (see Fig. 12d). There seems to be no single parameter that controls the relative importance of the dispersive transport, at least over the real urban canopies studied here.

Conclusion
This study analyses budgets of double-averaged second-order turbulence moments over a real urban canopy using large-eddy simulation. We focus on the budgets above the mean building height, where residual terms are generally negligible. The TKE budget shows that shear production is the primary source of TKE, whereas dissipation is the primary sink. Interestingly, wake production is also an important contribution to the TKE budget.
The pressure-strain correlation terms play an essential role in the velocity variance budgets. These terms redistribute energy between velocity components, thereby driving turbulence to the isotropic state. Over the considered urban canopy, pressure-strain correlation terms are responsible for the growth of the vertical-velocity variance at the expense of the streamwise velocity variance, as commonly observed in shear flows.
Along with the shear production term, the pressure-strain correlation term plays a vital role in the budget of momentum flux, where turbulent and pressure transport terms appear to be of secondary importance. The budget of scalar flux is dominated by the shear production and pressure gradient-scalar interaction terms, while the turbulent transport appears to be of secondary importance. However, along with the shear production and the scalar dissipation terms, the wake production, and turbulent and dispersive transport terms are important for the budget of scalar variances in the 2 ≤ z/H ≤ 5 interval.
In addition to the above analysis, we also examine the relative importance of the dispersive terms to the corresponding Reynolds terms in our model domain and in a range of subdomains. To achieve this, our model domain is partitioned into four subdomains in the y-direction. For each case, the ratio of wake production to shear production and the ratio of dispersive transport to turbulent transport averaged over different z/H intervals, are examined for TKE, momentum flux, scalar flux and scalar variance budgets. The importance of wake production of TKE and scalar variances decreases with height, and this importance appears to depend on the maximum building height (or the ratio of maximum building height to the mean building height), although more investigations are needed to confirm this. Wake production is less significant for momentum and scalar flux budget equations. The dispersive transport terms can be significant relative to their turbulent counterpart, but we could not identify any trend of how these terms vary as a function of the morphological parameters over the considered urban canopies.
Results from this work have implications for both single-layer and multi-layer urban canopy parameterizations, which have been developed to represent the flow and transport within and above neighborhoods in NWP and global climate models. Both single-layer (Masson 2000; Kusaka et al. 2001) and multi-layer (Martilli et al. 2002;Schoetter et al 2020) urban canopy parameterizations often assume horizontal homogeneity for canopies and neglect dispersive fluxes and dispersive transport. Our work indicates that dispersive fluxes (and dispersive transport) over real urban canopies can be important even above the mean building height. For single-layer urban canopy parameterizations coupled to an atmospheric model, this finding raises the question of whether dispersive fluxes should be parameterized, in addition to turbulent fluxes, despite that the lowest atmospheric grid is often above the mean building height. For multi-layer urban canopy parameterizations, our study supports and complements recent work that emphasizes the importance of dispersive stress relative to turbulent stress and the role of wake production in the TKE budget over idealized urban canopies (Nazarian et al. 2020). Our results further highlight that multi-layer urban canopy parameterization should properly consider the dissimilarity between momentum and scalar transport over real urban canopies. Findings from this work are limited to neutrally stratified ambient conditions; further investigations are needed to examine the impact of thermal stratification (stable or unstable) on the considered flow statistics.

Appendix 1: Derivation of the Double-Averaged Second-Order Moment Budget Equations
In this appendix, we show how the budgets of double-averaged second-order moments are obtained. For simplicity, we only derive the budget equation for the double-averaged Reynold stress tensor (Eq. 5) , scalar flux (Eq. 8) and scalar variance (Eq. 9). Budget equations for other second-order moments can be obtained in a similar fashion. Note that all budget equations are derived for neutral conditions with the Boussinesq approximation.
The LES resolved Reynold stress tensor u i u k budget equation is given as: where is the SGS stress tensor and ν t represents the SGS eddy viscosity. The first term on the left-hand side of Eq. 10 represents the local change of u i u k while the second is the advection of u i u k . On the right-hand side, the first two terms are the production terms resulting from the interaction of the mean flow and turbulence while the third term can be interpreted as the transport of u i u k by turbulent fluctuations (i.e., the turbulent transport term). The fourth term represents the interaction of the fluctuating pressure and velocity fields while the last term is the SGS term. After some algebraic manipulation on the last term, we have: Now the pressure term is split into the pressure transport term (the fourth term on the right-hand side) and the pressure-strain correlation term (the fifth term on the right-hand side). The SGS term also includes four terms: the SGS diffusion terms (the sixth and seventh term on the right-hand side) and the SGS dissipation terms (the eighth and nineth term on the right-hand side). To facilitate derivations, we write the advection term on the left-hand side of the above equation in its flux form by invoking the Boussinesq approximation: Applying the intrinsic spatial averaging to the above equation and following the rules in Eqs. 2 and 4, we have: Note that the surface integral does not show up in the above equation due to the no-slip boundary conditions. We can further expand the advection term as: Substituting into Eq. 13 gives: When we apply Eq. 15 to diagnosing our LES outputs, further simplifications can be made. First, we assume horizontal homogeneity at scales beyond the spatial averaging scale, consistent with the doubly periodic boundary conditions used in our LES. This assumption does not imply that the flow is horizontally homogeneous at the grid cell scale. In fact, the dispersive terms would not exist if horizontal homogeneity at the grid cell scale was assumed. Rather, this assumption means that the volume-averaged flow fields (represented by . ) have no horizontal gradients. Therefore, ∂ . ∂ x = ∂ . ∂ y = 0. Second, we analyze the budgets at stationary conditions and hence ∂ . ∂t = 0. Third, due to the use of doubly periodic boundary condition and continuity, the mean vertical velocity is zero (i.e., no large-scale subsidence, w = 0). Furthermore, because of the assumption of horizontal homogeneity at scales beyond the spatial averaging scale, α p becomes also only a function of z and thus only ∂α p ∂z = dα p dz is non-zero. With these assumptions, we have: To facilitate our analysis, we group and name the terms as follows: where P s ik is the shear production term, P w ik is the wake (dispersive) production term, P m ik is the work of the temporally averaged velocity fluctuations against the shear production (given that α p varies with height, P m ik = 0, see Eq. 5), T t ik is the turbulent transport term, T d ik is dispersive transport term, T p ik is the pressure transport term, S P ik is the pressure-strain correlation, D ik is the SGS transport term and ε ik is the dissipation term. From the budget equation for the Reynold stress tensor above, we can obtain the momentum flux w u and velocity variances budget equations. Similar steps can be followed to obtain the scalar flux budget equation as follows. The LES resolved u i s budget equation is given as: where τ SG S s, j = −k S ∂s ∂ x j is the SGS scalar flux and k S represents the scalar diffusivity. All terms in Eq. 18 have similar definitions to those in Eq. 10. On the left-hand side of Eq. 18, we have the local change and advection terms. On the right-hand side, there are two production terms, a transport term, a pressure gradient-scalar interaction term and two SGS terms. After some algebraic manipulation of the last two terms and writing the advection term in its flux form, we have: Applying the intrinsic spatial averaging to the above equation and following the rules in Eqs. 2 and 4, we have: We can further expand the advection term as: Substituting into Eq. 20 gives: When we make further simplification to Eq. 22 by assuming horizontal homogeneity at scales beyond the spatial averaging scale, stationary conditions and no large-scale subsidence, we have: Focusing on the vertical scalar flux (namely, i = 3) and noticing that u j u i s + w u j s + w τ SG S s, j n j d A = 0 because u j = u j = 0 at the fluid-solid interface, the above equation becomes: where the definitions for P s SF , P w SF , P m SF , T t SF and T d SF are similar to those in Eq. 17. P S is the pressure gradient-scalar interaction and D SF include the SGS terms; ψ SF is the surface integral term that arises from the averaging theorem.
The scalar variance budget equation can be obtained using similar steps. The LES resolved s 2 budget equation is given as: All terms in Eq. 25 have similar definitions to those in Eq. 10. The left-hand side of Eq. 25 has the local change and advection terms while the right-hand side has the production term, transport term and SGS term. After some algebraic manipulation of the last term and writing the advection term in its flux:form, we have Applying the intrinsic spatial averaging to the above equation and following the rules in Eqs. 2 and 4, we have: We can further expand the advection term as: When we make further simplification to Eq. 29 by assuming horizontal homogeneity at scales beyond the spatial averaging scale, stationary conditions and no large-scale subsidence, we have: Noticing that (u j s 2 +s 2 u j )n j d A = 0 because u j = u j = 0 at the fluid-solid interface, the above equation becomes: where the definitions for P s SV , P w SV , P m SV , T t SV , T d SV and D SV are similar to those in Eq. 17; ε SV is the scalar dissipation term; ψ SV is the surface integral term that arises from the averaging theorem.