A Note on Friction Velocity and Viscous Effect for Idealized Urban Canopy Flows

Friction velocity is one of the important scaling parameters in atmospheric boundary layer studies. However, several definitions of friction velocity exist in the literature: e.g. estimated from the total drag force or used only pressure drag, etc. In this study, a series of large-eddy simulation (LES) calculations were carried out to evaluate the impact of various definitions on the friction velocity for idealized urban geometries, i.e. staggered array of cubes with different packing densities and wind directions. We first compared the normalized velocity fields with the literature data for the case with a packing density of 25%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$25\%$$\end{document}. The results show that the LES data normalized by the friction velocity derived from the Reynolds stress using the extrinsic spatial average is more consistent with the direct numerical simulation data. Furthermore, when varying the wind direction, the distribution of Reynolds stress and pressure drag show significant change in streamwise and spanwise directions. We further found that for packing density of 44.44%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$44.44\%$$\end{document}, the frictional drag accounts for more than a quarter of the total drag, and even higher than the pressure drag in parallel wind direction. This leads to the deviation of friction velocity estimated from the pressure drag and that calculated from the total drag force up to 33%. Such characteristic of viscous effect challenges the assumption widely used in wind tunnel experiments and urban canopy parameterizations that the contribution of viscous force is negligible, especially for ultra-dense arrays.


Introduction
Friction velocity is widely used in urban-type boundary layer studies to normalize the velocity fields (Cheng and Castro 2002;Castro et al. 2006;Coceal et al. 2006Coceal et al. , 2007aKono et al. 2010;Leonardi and Castro 2010;Tomas et al. 2016;Blackman et al. 2017;Li and Bou-Zeid 2019;Tian et al. 2021;Xu et al. 2021;Akinlabi et al. 2022;Li and Katul 2022). Although the definitions of friction velocity in literature are different, the normalized numerical results are in satisfactory agreement with the experimental data. It must be noted that the experimental constraints often prevent measuring the viscous force on the solid wall leads to the friction velocity in wind-tunnel experiments is generally derived by measuring the pressure drag only (e.g. Cheng and Castro 2002;Castro et al. 2006;Cheng et al. 2007;Perret et al. 2019). Although the DNS work of Coceal et al. (2006) has confirmed that the viscous drag is negligibly small, the pressure drag shown in Leonardi and Castro (2010) is gradually reduced with the packing density (λ p ) increases, which implies that the viscous contribution may not be entirely insignificant for high packing density configurations (Leonardi et al. 2007).
Apart from using the pressure drag, existing studies have many other ways to define the friction velocity: e.g. derived from the maximum value of the Reynolds stress. However, this value depends largely on how spatial averaging is performed, namely the extrinsic spatial average on the comprehensive volume, or the intrinsic spatial average in which only the fluid volume is included. Xie and Fuka (2018) pointed out that the extrinsic spatial average approach yields a continuous vertical momentum flux within and above the urban canopy. However, Schmid et al. (2019) claimed that the continuity of the averaged quantities depends on the boundary condition at the top of the cube. Only when the quantity is zero at the fluid and solid interface, the intrinsic spatial average undergoes a step change. ( Blunn et al. (2022)) found that if the solid fraction is not accounted for in the intrinsic mixing-length closure, the error of the predicted mean velocity up to 17%. Note that it is rarely reported in the literature on which spatial average process should be used to calculate the friction velocity, and moreover (very) few studies devotes to the cases of oblique wind direction. Recently, several numerical studies (Tian et al. 2021;Akinlabi et al. 2022) performed on the staggered cube arrays, in which the friction velocity is derived from the total drag force, i.e. sum of the pressure drag and frictional drag. They showed that the normalized mean flow are slightly different from the experimental data of Castro et al. (2006). In addition, some authors (Coceal et al. 2006;Leonardi and Castro 2010;Claus et al. 2012) have pointed out that the wall friction velocity reduced by a factor √ 1 − d/H can yield the most satisfactory fit to the experimental velocity data.
As reported in the aforementioned literature, the definition of friction velocity is not unique around the urban boundary layer studies, which leaves an open question regarding the statistics normalization: to what extent does the friction velocity calculated from these definitions differ from that calculated from the total drag force? This note therefore focuses on this issue and attempts to provide a detailed comparison. We propose to set up and run the building-resolved large-eddy simulation model over a range of constant height staggered cube arrays with different packing densities and wind directions. The definitions of friction velocity are introduced in Sect. 2. The numerical methods and physical models are described in Sect. 3 and the results are presented in Sect. 4. Finally, conclusions are summarized in Sect. 5.

Friction Velocity Definitions
Friction velocity u * (also known as u τ ), as a characteristic velocity, is widely used in atmospheric boundary layer studies. It represents the effect of wind stress on the underlying surface, and varies with the nature of roughness (Kaimal and Finnigan 1994). In this section, we introduce four definitions of friction velocity that were previously used in literature.

By Total Drag Force
The friction velocity defined here is derived from the time-averaged total drag force (see Kono et al. (2010); Tian et al. (2021) and Akinlabi et al. (2022)), which is represented by: where ρ is the reference density and τ = F total /A total , A total is total plan area and F total is the total drag force, i.e. sum of pressure drag on the cube (D pc ), frictional drag on the cube (D f c ) and on the ground (D f g ).

By Pressure Drag
Except using the Karman momentum integral equation (Kármán 1921), in some wind-tunnel experiments (e.g. Castro et al. (2006); Cheng et al. (2007); Herpin et al. (2018); Perret et al. (2019)), the authors assume that the contribution of viscous force to total drag is negligible with the pressure drag (or called form drag) being the major contributor. Hence, the friction velocity is estimated by measuring the pressure drag only, which gives:

By Zero-Plane Displacement Height
Some authors (e.g. Coceal et al. 2006;Leonardi and Castro 2010;Claus et al. 2012) pointed out that the friction velocity reduced from u τ by a factor √ 1 − d/H yields the most satisfactory fit to the experimental velocity data of Cheng and Castro (2002), i.e.: where H is the domain height and d demotes the zero-plane displacement height. Following (Jackson 1981), d is the effective height at which pressure drag acts (Coceal et al. 2006;Perret et al. 2019), where h is the height of cubes, as indicated in Fig. 1b. Leonardi and Castro (2010) pointed out that, the difference between u * d and u τ is less than 4% in their DNS cases.

By Reynolds Stress
Some studies (e.g. Böhm et al. 2013;Cheng and Porté-Agel 2015;Blackman et al. 2017) assume a constant momentum flux layer, and the friction velocity u * is estimated using the maximum value of the Reynolds stress near the canopy top, which is defined as: where overbar and angular brackets denote the time and spatial averaging, respectively. In fact, the aforementioned definition uses only the longitudinal component of the Reynolds stress, which is primarily introduced in flat plate boundary-layer study of Prandtl (1949) where the Reynolds stress is parallel to the mean flow direction. However, this definition may not true for oblique atmospheric flows (Weber 1999). In Landau and Lifshitz (1987), they discussed the near-wall turbulence, and suggested that the horizontal Reynolds stress should be used to define the friction velocity, i.e.: It is worth noting that there has been continuous debate in the literature on intrinsic or extrinsic spatial average should be chosen for the urban-type modeling (Kono et al. 2010;Böhm et al. 2013;Castro 2017;Xie and Fuka 2018;Schmid et al. 2019;Sützl et al. 2021). The discussion mainly focus on the spatial averaging only in fluid region or over the whole domain including the blocks part. In general, the spatial averaging of a flow property ϕ over a volume V can be expressed as: When V corresponds to the entire domain volume including both fluid (V f ) and block (V b ) regions, this is extrinsic spatial average and defined as: When one considers only the fluid volume (V = V f ), then the spatial average is intrinsic spatial average, and defined as: There is no difference between these two average processes above the canopy because V b = 0 there. The main difference occurs inside the canopy, namely: where = V f /V is the fluid volume fraction. In what follows, these two spatial average processes are employed on Reynolds stress, which yields u * I (intrinsic spatial average) and u * E (extrinsic spatial average), respectively.

Methodology
The present LES approach is performed using the filtered Navier-Stokes equations for incompressible flows, as defined by the following equations: where u i (for i = 1 (u), 2 (v), 3 (w)) are the filtered velocity components, p is the filtered pressure. ρ is the reference density and ν is the kinematic viscosity and τ i j = u i u j − u i u j is the subgrid scale (SGS) stress. In this work, the dynamic Smagorinsky model developed by Germano et al. (1991) and Lilly (1992) is chosen in order to allow the energy backscattering from subgrid scales. Note that the SGS contribution has been added to the resolved Reynolds stress to obtain the friction velocity, which accounts for up to about 2% of the total stress. The mesh sets to be uniform ( = h/32) in all directions from ground up to z = 1.5h. For z > 1.5h, the mesh size is constant in the horizontal directions but gradually stretched in the vertical direction. The average vertical mesh size is h/10 in the region 1.5 < z/h < 5 and h/2 in the region 5 < z/h < 8 with a maximum value of 2h/3 at the top of the domain. The number of cells in the streamwise, spanwise and vertical directions is N x = 512, N y = 384 and N z = 91, respectively. The PhD thesis of Tian (2018) have demonstrated that at least 1/128h resolution is able to have three mesh layers in the range z + < 5 in a smaller domain simulation with dimensions Length = W idth = Height = 4 h. However, considering the prohibitive computational cost, it is not possible to do this on the current size of domain (16h × 12h × 8h). Coceal et al. (2007b) have pointed out that the turbulence is mainly  Coceal et al. (2007b). Circles: LDV data from Castro et al. (2006). Squares: LDV data from Herpin et al. (2018). Triangles: PIV data from Blackman et al. (2017). Stars: PIV data from Blackman et al. (2017). Note that the different friction velocity definition of the x-axis label is only applicable to LES data dominated by large-scale eddies shed from the sharp edges of the cubes, rather than the thin boundary layers on the cube facets. Besides, good agreements between the simulation with 1/32h resolution and 1/64h resolution are observed in their study. In this work, the minimum non-dimensional wall distance z + = zu * ν+ν sgs at location P1 is 1.66, and at location P0 is 2.97 for the case of λ p = 25% with perpendicular wind direction, where u * is the local friction velocity determined in the viscous sublayer by u 2 * = (ν + ν sgs ) ∂ u ∂z . Note that the viscous stress is calculated by: Although there are only two mesh layers in the range z + < 5, we assume that the current grid resolution of 1/32h is fine enough to capture the velocity gradient inside the viscous sublayer. Furthermore, all simulations are run for an initial duration of 300T, where T = h/u τ is the turnover time for the typical eddies shed by the cubes. Statistics are then averaged over a further 500T, to ensure convergence, with a time-step of 0.002T. Coceal et al. (2006) and Claus et al. (2012) have pointed out that an averaging time of about 400T at least is necessary because of the possible existence of longitudinal rolls above the canopy, which move around over long time scales and give rise to non-zero residual dispersive stresses if the averaging time is too short. Each run in this study is performed using 160 processors on 'Taiyi' supercomputing cluster at the SUSTech and 600 h computing time.

Normalization by the Friction Velocity
Figure 2shows the vertical profiles of mean and standard deviation of streamwise velocity component for the case of λ p = 25% at location P1, normalized by the different friction velocities, respectively. From these comparisons, it can be clearly seen that the LES data normalized by the friction velocity derived from the Reynolds stress using the extrinsic spatial average process (u * E ) is more consistent with the DNS data than other definitions (see Fig. 2a5). Such good agreement is presumably because the friction velocity in DNS is obtained from the maximum Reynolds stress at the top of cubes, which yields a similar value to u * E . In addition, the strong velocity gradient above the canopy in the experimental data of Castro et al. (2006) is also well-matched by the LES data using u * E (see Fig. 2b5). However, Fig. 2a1, b1 show that the experimental mean velocity of Castro et al. (2006) are  Fig. 3 The contribution of pressure drag and frictional drag for the three investigated configurations, packing density of λ p = 6.25%, 25% and 44.44%, at wind direction θ = 0 • (a), θ = 45 • (b) and θ = 90 • (c), respectively underestimated by the LES data using the friction velocity u τ . It is noteworthy that comparing to the differences shown in the experimental data, the discrepancy in these normalized LES profiles is less significant. As reported in Tian et al. (2021), the growing misalignment of the mean velocity between the LES results and PIV measurements above the canopy may be attributed to the difference of relative boundary-layer height, while the deviation of σ u between LES data and LDV measurements within the canopy may be due to the lack of spatial resolution in numerical simulations.

Horizontal Contribution for Oblique Wind Directions
As shown in Table 1, when the flow direction is normal to the obstacle surface, the Reynolds stress in streamwise direction (P rest ) for three array types λ p = 6.25%, 25% and 44.44% account for 99.99%. This suggests that it is feasible to use only the streamwise component of Reynolds stress (−u w ) to calculate the friction velocity in case of perpendicular wind direction (θ = 0 • ). However, with the deflection of wind direction, it is found that the percentage of Reynolds stress in streamwise direction (P rest ) is gradually decreasing, while the contribution in spanwise direction (P resp ) is continuously growing. Hence, if Reynolds stress is selected to calculate the friction velocity in oblique wind direction, the horizontal Reynolds stress, i.e. [(−u w ) 2 + (−v w ) 2 ] 1 4 , is more appropriate. A similar trend is found in pressure drag on the cube (D pc ), that is, the percentage of pressure drag on the horizontal plane depends largely on the prevailing wind direction. Since the pressure drag is determined by the integral of local pressure times the surface area of the cube, it is clear from the data in Table 1 that the local pressure on the windward and leeward sides of the cube is closely related to the upstream airflow. It might therefore anticipate that the friction velocity derived from Eqs. (3) and (4) using only the pressure drag in streamwise direction would yield an underestimated value when the wind direction is oblique.

Non-negligible Viscous Effect
Figure 3a-c show the contribution of drag forces for the packing density of 6.25%, 25% and 44.44%, at wind directions θ = 0 • , θ = 45 • and θ = 90 • , respectively. It can be clearly seen that, for less dense arrays (λ p = 6.25% and λ p = 25%), the frictional drag has quantitatively different behavior. The frictional drag on ground (D f g ) is more significant than the frictional drag on cube (D f c ) for packing density of 6.25%, by contrast, D f g is only a few percent of the total drag (less than 4%) in the case of λ p = 25%. Note that, in most cases studied here,  Fig. 3c). In particular, the total frictional drag increased approximately threefold from λ p = 25% to λ p = 44.44%. This significant growth of the frictional drag may be attributed to the higher blockage in λ p = 44.44%, which increases the shelter volume in the flow and reduces the pressure drag on the cube ( (Raupach 1992)). It is worth noting that, scientific community widely accepted the view that the contribution of viscous force to total drag is negligible. However, as shown here, this hypothesis is only valid for less dense arrays, namely cases of packing density of 6.25% (isolated-wake flow) and 25% (wake-interference flow). In terms of the ultra-dense arrangements (skimming flow), the effect of viscous force, especially the frictional drag on the cube, is particularly significant.

Conclusion and Discussion
Friction velocity as one of the important scaling parameters is not uniquely defined in urban boundary-layer studies. In this work, we conducted LES calculations over a range of staggered cube arrays with three packing density, λ p = 6.25%, 25% and 44.44%, in different wind directions, respectively. A significant viscous effect is found in packing density of λ p = 44.44%, where the viscous force accounts for more than 25% of the total drag. This result challenges the widely used assumption in wind tunnel experiments and urban canopy parameterizations that the viscous contributions can be ignored. It is worthy noting that the viscous effect problems have not received much attention previously. A plausible reason is that, the experimental measurements and DNS data for less dense staggered configuration (λ p ≤ 25%) support the view that viscous contribution is not significant. Furthermore, for experimentalists, it is non-trivial to measure the viscous force on the solid wall, so that there is few work to directly measure the frictional resistance on the roughness element and ground. Therefore, this note puts forward the standpoint that viscous force is a non-negligible part in much dense arrangement, and there is a risk of misalignment in numerical validation when compared with the experimental data normalized by the friction velocity derived from pressure drag only.
In terms of friction velocity, it can be found in Table 1 that the friction velocity derived from total drag forces (u τ ) yields the largest value among all definitions, while the friction velocity estimated from the Reynolds stress using the intrinsic spatial average (u * I ) is larger than that using the extrinsic spatial average (u * E ), except for the cases of parallel wind direction (θ = 90 • ). Note that the wind direction θ = 90 • corresponds to the configuration that Coceal et al. (2006) termed 'aligned' array. Cheng and Castro (2002) pointed out that the viscous drag may not be entirely negligible for the aligned array because they found differences in friction velocities. In this work, a significant decrease observed between u τ and u τ p in the case of λ p = 44.44%, in which the maximum deviation exceeds 33%. Furthermore, it is interesting that for the case of λ p = 25% at perpendicular wind direction, the friction velocity u τ p , u * d , u * I and u * E are almost indistinguishable lying in the range of 0.37−0.38 (m/s). Such small difference (in friction velocity) makes the previous research did not pay much attention to the definition of friction velocity when validating with the literature data, and consequently ignoring the effect of viscous force.
Finally, it is worth emphasizing that the contribution of viscous force is not entirely insignificant in ultra-dense arrangements. Switching the definition of friction velocity from total drag force to pressure drag only is likely to yield unexpected deviation, which should be noted in future wind tunnel experiments and numerical modelings.