Evidence for Raupach et al.'s mixing-layer analogy in deep homogeneous urban-canopy flows

Abstract

The mixing-layer analogy is due to Raupach, Finnigan & Brunet (Boundary-Layer Meteorol., vol. 25, 1996, pp. 351–382). In the analogy, the flow in the roughness sublayer of a homogeneous deep vegetation canopy boundary layer is analogous to a plane mixing layer rather than a surface layer. Evidence for the analogy includes the inflected velocity profile, which resembles the velocity profile in a plane mixing layer, and, most notably, the following estimate as a result of the Kelvin–Helmholtz instability: $\varLambda _x=8.1L_s$, where $\varLambda _x$ is the spacing of the large-scale eddies, and $L_s$ is the shear length. The mixing-layer analogy has been very successful in vegetation canopy flow research, but has received only limited support in urban-canopy flow research. This work revisits Raupach et al.'s mixing-layer analogy, and we present the evidence for the mixing-layer analogy in urban-canopy flows: the exponential velocity profile in the canopy layer, i.e. $(U-U_c)/(U_h-U_c)=\exp (z/L_m)$, and $L_m\sim [(U_h/U_c-1)(U_h/U_c+3)]^{-1}$. Here, $z$ is the vertical coordinate, $L_m$ is the attenuation length and is a measure of the largest eddy in the canopy layer, $U_h$ is the wind speed at the canopy crest and $U_c$ is the velocity in the inactive layer. We conduct direct numerical simulations of various deep homogeneous urban-canopy flows and test the above two scalings. We also discuss why Raupach et al.'s analogy has not seen as many successes in urban-canopy flows as in vegetation canopy flows.

1. Introduction

Canopy flows are common in the lower part of the atmosphere, and they have been the topic of many reviews (Finnigan 2000; Barlow & Coceal 2009; Fernando 2010; Belcher, Harman & Finnigan 2012; Nepf 2012; Brunet 2020). Figure 1(a) shows the layered structure of a vegetation canopy flow. It consists of the roughness sublayer, the logarithmic layer and the wake layer. The definition of the roughness sublayer differs from one community to another (Finnigan 2000; Fernando 2010; Flack & Schultz 2014; Chung et al. 2021; Flack & Chung 2022). Here, the roughness sublayer contains the canopy layer. The behaviour of the mean flow plays an important role in defining the layers in a rough-wall boundary layer. The mean velocity profile in the logarithmic layer follows the logarithmic law of the wall (Raupach & Thom 1981), and the mean flow in the wake layer abides by Townsend's outer layer similarity (Townsend 1976; Jiménez 2004). Lastly, Raupach et al.'s mixing-layer analogy gives an estimate of the flow in the canopy sublayer (Raupach, Finnigan & Brunet 1989, 1996).

Figure 1. (a) A sketch of a vegetation canopy flow. The flow consists of the wake layer, the logarithmic layer and the roughness sublayer. Here, $\delta$ is the height of the boundary layer above the canopy, $h$ is the height of the canopy and $U_h$ is the velocity (wind speed) at the top of the canopy. The inflected velocity profile in the roughness sublayer gives rise to K–H-type waves above the canopy. (b) A sketch of a plane mixing layer. The inflected velocity profile gives rise to K–H instability. Here, ${\rm \Delta} U$ is the velocity difference between the two streams. A vertical ($z$) coordinate is defined such that the origin is at the top of the canopy, which corresponds to the centre of the mixing layer. The flow is symmetric with respect to $z=0$.

1.1. Mixing-layer analogy and vegetation canopy flows

First, we explain what the mixing-layer analogy is. Raupach et al. (1989) noticed that the inflected velocity profile in the roughness sublayer is similar to that in a plane mixing layer and hypothesized that the roughness sublayer is analogous to a plane mixing layer. This hypothesis was elaborated in Raupach et al. (1996), according to which the inflected velocity profile in the roughness sublayer gives rise to two-dimensional Kelvin–Helmholtz (K–H) waves.

An extensive review of the mixing-layer literature and the K–H instability falls outside the scope of this work, but a summary of the basics helps the discussion that follows. Consider a plane mixing layer, as sketched in figure 1(b). The inflected velocity gives rise to a shear length scale

(1.1)\begin{equation} L_{s,ML}=\frac{{\rm \Delta} U/2}{{{{\rm d}}U}/{{{\rm d}}z}|_{z=0}}, \end{equation}

where ${\rm \Delta} U$ is the velocity difference between the two streams. The factor 2 is included here to simplify the derivations in the following sections. Linear stability analysis gives an estimate of K–H waves’ spacing, i.e. $\varLambda _x=7L_s$ to $10L_s$.

Now, we consider the canopy flow. The velocity profile in the roughness sublayer is also inflected and should give rise to a similar shear length scale

(1.2)\begin{equation} L_s=\left.\frac{U}{{{{\rm d}}U}/{{{\rm d}}z}}\right|_{z=0}, \end{equation}

and the spacing of the transverse vortices in the roughness sublayer should be $\varLambda _x=7L_s$ to $10L_s$. Raupach et al. (1996) measured vertical velocity correlations in a number of canopy flows and concluded that

(1.3)\begin{equation} \varLambda_x=8.1L_s.\end{equation}

Equation (1.3) was found to work well in Novak et al. (2000), Dupont & Brunet (2008), Dupont & Patton (2012), Huang, Cassiani & Albertson (2009), Bailey & Stoll (2013) and Gavrilov et al. (2013) (to name a few) for homogeneous dense vegetation canopies, and is the cornerstone of the mixing-layer analogy. Raupach et al.'s mixing-layer analogy is also evidenced by similar flow structures in other boundary-layer flows and mixing layers (Finnigan 2000; Finnigan, Shaw & Patton 2009; Huang et al. 2009).

Deviations from (1.3) are found in stably stratified canopy flows (Brunet & Irvine 2000), in non-homogeneous canopies (Thomas & Foken 2007), in crop canopies (Py, de Langre & Moulia 2006; Dupont et al. 2010; Tschisgale et al. 2021) and in sparse canopies (Novak et al. 2000; Huang et al. 2009), where the flow is akin to a surface layer rather than a mixing layer. However, these deviations are not unexpected since the mixing-layer analogy is not meant for these flows.

1.2. Mixing-layer analogy and urban-canopy flows

While the vegetation canopy flow community has actively scrutinized the mixing-layer analogy, the urban-canopy flow community has not given nearly as much attention to Raupach et al. (1996). This is partly because the mixing-layer analogy is meant for deep homogeneous canopies, but deep urban canopies are not common and are found only in metropolitan areas (Barlow & Coceal 2009). Figure 2 is a sketch of a deep homogeneous urban-canopy flow. Surface roughnesses in an urban canopy are buildings, which are modelled as rectangular prisms. Comparing figures 1(a) and 2, we see that a homogeneous deep vegetation canopy flow and a homogeneous deep urban-canopy flows share a lot in common. They both have a wake layer and a logarithmic layer above the roughness sublayer. Moreover, they both have an inflected mean velocity within the roughness sublayer. It would seem that Raupach et al.'s mixing layer analogy would work well for the flow in urban-canopy flows.

Figure 2. A sketch of an urban-canopy flow. The flow consists of the wake layer, the logarithmic layer and the roughness sublayer. The roughness sublayer contains the canopy occupied region.

However, Raupach et al.'s mixing-layer analogy has not seen much success in the urban canopy flows. In the following, we explain why. Consider the urban canopy in figure 2. Coceal et al. (2007) argued that the inflection in the velocity profile at the top of the canopy is a result of the no-slip condition at the top surface of each individual rectangular prism rather than the roughness canopy as a whole. Coceal et al. (2007) further argued that, in an urban canopy, mixing layers are only found locally, and there are no large-scale K–H waves in urban-canopy flows. Coceal et al.'s (2007) argument is supported by the available empirical evidence. Huq et al. (2007) searched their experimental data but could not find K–H rollers despite extensive visualization attempts. Figure 3 shows the instantaneous vortical structures in a few urban canopies flows obtained from direct numerical simulations. The flow configuration is a half-channel, and the canopy consists of rectangular prisms with the aspect ratio $h/w_r=4$, where $h$ and $w_r$ are the height and width of the prisms, respectively. The ground coverage density is $\lambda _p=0.06$, 0.11, 0.25 and 0.39 in figure 3(a–d), covering a large range of $\lambda _p$ values. Like the previous authors (Kanda, Moriwaki & Kasamatsu 2004), we cannot find large-scale K–H waves. This becomes more clear if we compare figures 4 and 5, where we show the premultiplied spectra in a canopy flow and a plane channel flow. The two spectra look very much alike, and there is no energy accumulation at any specific wavelength, particularly at the wavelength that corresponds to K–H instability. It is worth mentioning that we choose to show the results in figures 3 and 4 here rather than in the result section because of the following. First, the results in figures 3 and 4 explain the difficulty of Raupach et al.'s mixing-layer analogy in urban-canopy flows and motivate the present work. Second, similar results like the ones in figures 3 and 4 are already extensively available in the present literature and should not need much explanation.

Figure 3. Visualization of instantaneous vortical structures in urban-canopy flows. Shown here are $Q$ iso-surfaces (Hunt, Wray & Moin 1988) coloured by their distances from the ground. The roughnesses are rectangular prisms, whose dimensions are ${w_r} \times {w_r}\times 4{w_r}$. Here, the vertical coordinate is such that the origin is at the ground. The surface coverage densities are (a) $\lambda _p=0.06$, (b) $\lambda _p=0.11$, (c) $\lambda _p=0.25$ and (d) $\lambda _p=0.39$. The flows correspond to R3L06/11/25/39A (see § 4 for details of the direct numerical simulations).

Figure 4. Premultiplied energy spectra of the vertical velocity, i.e. $k_xk_yE_{ww}$ at a distance $z^{+}=15$ from the canopy crest. The results in panels (a–d) correspond to the flows in figure 3(a–d).

Figure 5. Premultiplied vertical velocity energy spectra, i.e. $k_xk_yE_{ww}$ in plane channel flow of $Re_{\tau }=550$ (where the friction Reynolds number $Re_\tau$ is defined based on the friction velocity and the half channel height) at $z^{+}=60$ above the wall. Figure reproduced using data reported in del Alamo et al. (2004).

Although (1.3) does not hold for homogeneous urban-canopy flows, we cannot hasten to the conclusion that Raupach et al.'s mixing-layer analogy fails. Here, we explain why. We begin by reviewing the historical background of Raupach et al. (1989) and Raupach et al. (1996). Before the establishment of the mixing-layer analogy, the conventional view was that the roughness sublayer of a vegetation canopy is a superposition of standard surface layer and small-scale turbulence generated in the wakes of plant elements (Raupach & Thom 1981). In the 1970s, it becomes evident that this view is false. For one, the velocity profile is inflected. For another, the velocity in the roughness sublayer is controlled by a single length scale rather than many scales as in a surface layer. The mixing-layer analogy was proposed in the 1980s to account for these facts. Raupach et al. (1996) put an emphasis on (1.3), but the foundations of the mixing layer analogy are the inflected velocity and the single length scale in the canopy layer, both of which show their presence in homogeneous deep urban-canopy flows. Following this logic, the mixing layer analogy should apply to urban-canopy flows. The challenge, however, is to find mathematical formulations one can put to the test.

1.3. This work

The objective of the present work is to find mathematical formulations that one can test to confirm/refute the mixing-layer analogy and then put these formulations to the test in direct numerical simulations (DNSs) of urban-canopy flows. DNS resolves all scales in a turbulent flow and is the most accurate computational fluid dynamics tool. A limitation of DNS is its high cost (Yang & Griffin 2021), and because of its high cost, the use of DNS has been limited to low and moderate Reynolds number flows (Moin & Mahesh 1998). Canopy flows in the lower part of the atmosphere are at high Reynolds numbers and are usually studied using wall-modelled large-eddy simulation (WMLES) (Giometto et al. 2016; Yang 2016b; Zhu et al. 2017; Zhu & Anderson 2018; Fan et al. 2021; Yang & Ge 2021; Zhang et al. 2021a,b). However, a WMLES does not resolve the flow close to the top surfaces of roughness elements and the thin shear layers downstream of the roughness elements, which are critical to the mixing layer analogy.

The discussion so far has focused on topics that are closely registered near the mixing layer analogy, and we have not covered the previous research on homogeneous deep urban canopy. Urban-canopy flow is, in general, a very well-researched topic (Barlow & Coceal 2009). Idealized urban canopies have been studied experimentally (Cheng & Castro 2002; Huq et al. 2007; Inagaki & Kanda 2008; Hagishima et al. 2009), computationally (Kanda et al. 2004; Coceal et al. 2006; Leonardi & Castro 2010; Lee, Sung & Krogstad 2011; Anderson, Li & Bou-Zeid 2015b; Cheng & Porté-Agel 2015; Li & Bou-Zeid 2019) and theoretically (Yang et al. 2016; Chung et al. 2021). However, the existing studies are, by and large, limited to low aspect ratio roughness elements, and deep canopy flow remains an insufficiently explored territory (Sadique et al. 2017; MacDonald et al. 2018). A more in-depth review of the urban-canopy literature falls outside the scope of this work, but in anticipation of § 2, reviewing the layered structures of the urban-canopy layer would be beneficial. Figure 2 is a sketch of the layered structure of the canopy layer. The roughness elements are rectangular prisms of the same height. The flow in the active layer directly exchanges momentum and energy with the flow above the canopy. The flow in the inactive layer, on the other hand, does not, at least not directly. The flow in figure 2 is pressure driven. The pressure force pushes the flow through the roughness in the inactive layer as in a porous medium. A ground layer emerges between the inactive layer and the ground connecting the finite fluid velocity in the inactive layer to 0 velocity at the ground level. In a zero-pressure-gradient (ZPG) deep canopy flow, the velocity in the inactive layer is 0 because of a lack of a driving force, and the inactive layer and the ground layer become one, as sketched in figure 6(b). Figure 6(c) is a sketch of a shallow urban canopy, where the flow in the canopy is affected by the boundary layer above the canopy as well as the ground. Evidently, the flow in a shallow urban canopy is more akin to a surface layer, where the mixing-layer analogy would not apply. This work focuses on the pressure-driven deep canopy flows, i.e. flow sketched in figure 6(a).

Figure 6. Sketches of urban-canopy flows: (a) pressure-driven deep urban-canopy flow, (b) ZPG deep urban-canopy flow, (c) shallow urban canopy.

The rest of the paper is organized as follows. In § 2, we show that the exponential profile is closely related to Raupach et al.'s mixing-layer analogy. In § 3, we derive mathematical formulations that can be tested to confirm/refute Raupach et al.'s mixing layer analogy. These mathematical formulations are put to the test in § 5 after we detail our DNS set-up in § 4. Further discussion of the results is included in § 6, followed by conclusions in § 7.

2. The exponential profile

In § 3, we will show that the velocity following an exponential scaling in the urban canopy is direct evidence for Raupach et al.'s mixing-layer analogy. In anticipation of § 3, we review Inoue (1963) and Cionco's (1965) original formulation of the exponential profile, explain why it has seen limited success (MacDonald 2000; Yang 2016a; Yang et al. 2016) and propose a new formulation.

2.1. Inoue and Cionco's original formulation

We begin by deriving Inoue and Cionco's original formulation. Consider a fully developed atmospheric boundary layer over a deep urban canopy. The momentum equation reads

(2.1)\begin{equation} -\frac{{{{\rm d}} \langle \overline{u'w'} \rangle }}{{{{\rm d}}z}} - {C_d}{U^{2}}=0, \end{equation}

in the canopy (Inoue 1963; Cionco 1965), where $\overline {u'w'}$ is the Reynolds shear stress, and $C_d U^{2}$ is the drag force exerted by the urban roughness. Here and throughout the paper, $u$, $v$ and $w$ are the velocity in the streamwise ($x$), spanwise ($z$) and wall-normal ($z$) directions, respectively, $\bar {{\cdot }}$ denotes time averaging, $\langle {\cdot } \rangle$ denotes comprehensive/superficial averaging in the $x$ and $y$ directions (Mignot, Barthelemy & Hurther 2009; Xie & Fuka 2018; Schmid et al. 2019). The averaging includes the solid volume within which the fluid velocity is assumed to be 0. Here, $C_d$ is the drag coefficient (and depends on the canopy layout), and $U=\langle \bar {u} \rangle$ is the double-averaged velocity in the $x$ direction. The mean flow advection, the viscous diffusion and the pressure force are neglected since the flow is fully developed, and is not subjected to any mean pressure gradient at a high Reynolds number. Invoking the eddy viscosity model and mixing length model, we have

(2.2)\begin{equation} -\langle \overline{u'w'} \rangle=\nu_t\frac{{{\rm d}}U}{{{\rm d}}z}= \left(l_m^{2}\left|\frac{{{\rm d}}U}{{{\rm d}}z}\right|\right)\frac{{{\rm d}}U}{{{\rm d}}z}, \end{equation}

where $\nu _t$ is the turbulent viscosity, and $l_m$ is the mixing length and is approximately a constant in the urban canopy (Bai, Meneveau & Katz 2012; Forooghi, Yang & Abkar 2020). Plugging (2.2) into (2.1), we have

(2.3)\begin{equation} \frac{{{\rm d}} }{{{{\rm d}}z}}\left[{\left( {l_m\frac{{{{\rm d}}U}}{{{{\rm d}}z}}} \right)^{2}}\right] = {C_d}{U^{2}}, \end{equation}

and one can easily verify that the following exponential profile is a solution to (2.3)

(2.4)\begin{equation} U = U_h \exp\left(a\frac{z}{h}\right)= U_h \exp\left(\frac{z}{L_m}\right). \end{equation}

Here, we put the origin of the $z$ coordinate at the top of the urban canopy, and the canopy occupies $-h< z<0$, $U_h$ is the velocity at top of the canopy (following the conventional notation), $L_m=h/a$ is the attenuation length scale and $a$ is the attenuation factor.

In the above derivation, Inoue (1963) and Cionco (1965) invoked the following assumptions in addition to a constant $C_d$ and a constant $l_m$. First, the viscous diffusion is negligible. Second, the mixing length $l_m$ is a constant. Third, the flow is not subjected to any mean pressure gradient. The exponential profile will not be valid if any of the above assumptions is not true.

2.2. A few caveats

There are a number of caveats when applying (2.4), which are often ignored (Castro 2017). First, (2.4) is meant for ZPG flows and is not valid for pressure-driven flows. To make this clear, let us consider a deep urban canopy, for which $h\gg L_m$. Let us consider a $z$ location that is sufficiently far away from the canopy crest and also sufficiently far away from the bottom wall, i.e. a $z$ location in the inactive layer. This should be possible if the canopy is deep. There, the vertical momentum flux is small, and according to (2.4), we should have $U=U_h\exp (z/L_m)\approx 0$. However, this is true only if it is a ZPG urban-canopy flow. When the mean pressure gradient is non-zero, the pressure force pushes the fluid through the urban canopy. The moving fluid gives rise to a non-zero drag force, which then balances the driving pressure force, leading to the following Reynolds-averaged momentum equation:

(2.5)\begin{equation} - \frac{{{{\rm d}}P}}{{{\rm d}\kern0.7pt x}} = \rho{C_d}{{U_c}^{2}}, \quad -h\ll z\ll 0, \end{equation}

where $P$ is the double-averaged pressure, and $U_c$ is the non-zero double-averaged velocity in the inactive layer. In light of the discussion above, the exponential profile should really be

(2.6)\begin{equation} \frac{U-U_c}{U_h-U_c}=\exp\left(\frac{z}{L_m}\right), \end{equation}

for pressure-driven urban-canopy flows. It is worth noting that (2.6) is a physical ansatz rather than a solution of the governing equation.

The second caveat is: (2.4) is valid only sufficiently far away from the ground, i.e. for $z$ such that $L_m\ll (z+h)$ (note that $z=0$ at the canopy top and $z=-h$ at the ground). This requirement arises because (2.4) is valid if and only if the mixing length $l_m$ is a constant. Since the ground constrains the mixing length (Townsend 1976; Marusic & Monty 2019; Yang & Meneveau 2019; Zhang et al. 2021b), (2.4) fails close to the ground. In a pressure-driven deep canopy boundary layer, a ground layer emerges close to the wall, and the fluid velocity decreases from $U_c$ to $0$ on the ground (Huang et al. 2009).

The third caveat is: (2.4) is valid only in regions where the viscous term is small. In vegetation canopy flows, this is barely an issue. The viscous term is small except very close to the ground, where (2.4) is not valid anyway. This, however, is not true for rectangular roughness arrays with smooth top surface. The no-slip top surfaces of these rectangular roughness elements give rise to a sharp velocity gradient at the top of the canopy, and viscous diffusion becomes the dominant term in the momentum equation especially when the roughness packing density is sufficiently high (Xu et al. 2021). For these rectangular roughness arrays, the origin of the exponential profile cannot be at the top of the canopy, and a displacement $d_c$ must be incorporated

(2.7)\begin{equation} \frac{U}{U_h}=\exp\left(\frac{z+d_c}{L_m}\right). \end{equation}

2.3. A new formulation

Taking into consideration these caveats in § 2.2, we propose a new formulation for homogeneous deep canopy flows,

(2.8)\begin{equation} \frac{U-U_c}{U_h-U_c}=\exp\left(\frac{z+d_c}{L_m}\right), \quad {\rm for}\ L_m\ll(z+h),\ z+d_c<0, \end{equation}

where, again, the coordinate is such that $z=0$ at the top of the canopy, $U_c$ accounts for a non-zero $-{{\rm d}}P/{{\rm d}\kern0.7pt x}$ and is such that $-\rho C_dU_c^{2}-\rho \, {{\rm d}}P/{{\rm d}\kern0.7pt x}=0$, $U_h$ is the double-averaged velocity at the top of the canopy and is a normalization velocity scale, $d_c$ accounts for the sharp velocity gradient at the top of canopy when the canopy is a rectangular roughness array of the same height and $L_m$ is the attenuation length scale. For a ZPG atmospheric boundary layer over vegetation canopies, $U_c=0$, $d_c\approx 0$, and (2.8) reduces to Inoue and Cionco's original formulation in (2.4). For dense and deep urban canopies, $d_c\approx 0$, $U_c\approx 0$, $z+d_c\approx z$ in most of the exponential layer, and the new formulation (2.8) reduces to Inoue and Cionco's original formulation as well. The new formulation (2.8) contains the following parameters: $U_c$, $d_c$ and $L_m$, and is applicable sufficiently far away from the ground and for $z<-d_c$.

3. Mixing-layer analogy in deep urban-canopy flow

In this section, we derive mathematical relations that can be tested to confirm or refute Raupach et al.'s mixing-layer analogy in the context of deep urban-canopy flows.

In a plane mixing layer, where the velocities of the two streams are $U_1$ and $U_2$, the velocity follows the universal scaling (Pope 2000)

(3.1)\begin{equation} \frac{{U - {U_2}}}{{{U_1-U_2}}} = f\left( {\frac{{z}+l}{{{{L_s}}}}} \right), \end{equation}

where $U_1>U_2$, ${L_s}$ is the shear length scale, $l$ is a displacement length and $f$ is a generic function (Rogers & Moser 1994). Invoking the mixing-layer analogy, the velocity in the urban canopy should follow a similar scaling. In other words,

(3.2)\begin{equation} \frac{{U - {U_c}}}{{{U_h} - {U_c}}} = f \left( {\frac{{z+l}}{{{L_s}}}} \right), \end{equation}

where $z$ is the vertical coordinate. Comparing the exponential profile in (2.8) and the ansatz in (3.2), we notice that the exponential profile is direct evidence for the mixing-layer analogy: the universal function $f$ in (3.2) is the exponential function $\exp ({\cdot })$ in (2.8), the displacement $l$ is $d$ and $L_s$ is $L_m$.

The most challenging part of applying the exponential profile is $L_m$. In the following, we estimate $L_m$. To do that, we write the force balance in a canopy flow

(3.3)\begin{equation} \int_{{-}h}^{0} \rho{C_d}U{{(z)}^{2}}\, {{\rm d}}z = {-} ({h+\delta}) \frac{{{{\rm d}}p}}{{{\rm d}\kern0.7pt x}}, \end{equation}

where the drag balances the driving pressure force, and we have neglected the skin friction. It is worth noting that (3.3) is a direct consequence of the mean momentum equation and can be obtained by integrating the mean momentum equation from the ground to the top boundary. Plugging (2.6) into (3.3), we have

(3.4)\begin{equation} {\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right)^{2}} \frac{{{L_m}}}{2}\left[ {1 - \exp \left( -{2\frac{{ h}}{{{L_m}}}} \right)} \right] + 2\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right){L_m}\left[ {1 - \exp \left(-{\frac{{h}}{{{L_m}}}} \right)} \right] = \delta. \end{equation}

In a deep canopy, $h \gg L_m$, $\exp (-h/L_m)\approx 0$ and (3.4) becomes

(3.5)\begin{equation} \frac{1}{2}{\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right)^{2}}{{L_m}} + 2\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right){L_m} = \delta. \end{equation}

Rearranging (3.5), we arrive at an estimate of the shear length scale

(3.6)\begin{equation} \frac{{{L_m}}}{{\delta}} \approx \frac{2}{{\left( {{U_h}/{U_c} - 1} \right)\left( {{U_h}/{U_c} + 3}\right)}}, \end{equation}

which should be true irrespective of the Reynolds number and other parameters. Plugging in (2.8) into (3.3) gives a longer expression but leads to the same conclusion after simplification as we will show in § 6.2.

Equations (3.2) and (3.6) are mathematical formulations that can be tested to confirm or refute Raupach et al.'s mixing-layer analogy in the context of deep canopy flows.

4. Direct numerical simulation

We conduct DNSs to test (3.2) and (3.6). In this section, we present the details of the DNS set-up. The pseudo-spectral code LESGO is employed for our DNSs. The code solves the incompressible Navier–Stokes equations in a periodic half-channel. It uses a pseudo-spectral method in the streamwise and spanwise directions and a second-order finite difference method in the wall-normal direction. A second-order Adams–Bashforth method is employed for time advancement. The flow is pressure driven. Canopy roughness is resolved using an immersed boundary method (Chester, Meneveau & Parlange 2007). Many have used this code for boundary-layer flow simulations, including flows over complex terrain (Anderson & Chamecki 2014; Anderson et al. 2015a), vegetative canopies (Chester et al. 2007; Bai et al. 2012) as well as urban canopies (Cheng & Porté-Agel 2015; Giometto et al. 2016; Yang et al. 2019).

We consider deep canopy flows. The roughnesses are rectangular prisms. Two roughness arrangements are considered, namely, aligned and staggered, as sketched in figure 7(a,b). Here, $s$ is the spacing of the prisms, and ${w_r}$ is the width of the prisms. In most cases, $h/{w_r}=4$ and $\delta /{w_r}=4$, 6. The spacing $s/{w_r}$ is 1.6, 2, 3, or 4. The surface coverage density, $\lambda _p={w_r}^{2}/s^{2}$, is 0.39, 0.25, 0.11 or 0.06. Figure 7(c) shows a side view of the flow domain.

Figure 7. (a,b) Top views of a repeating tile for the aligned/staggered configuration, (c) side view of the canopy and (d) perspective view of the canopy (here, $h/{w_r}=4$, $\lambda _p=0.25$).

Similar configurations are considered by Sharma & Garcia-Mayoral (2020) and MacDonald et al. (2018) in their numerical studies of deep roughness canopies. Figure 7(d) is a perspective view of a specific case, i.e. $h/{w_r}=4$, $\delta /{w_r}=4$ and $\lambda _p=0.25$.

In addition to the canopy geometry, the flow is controlled by the Reynolds number. The Reynolds number of a deep canopy boundary layer can be defined in a few ways (MacDonald et al. 2018; Xu et al. 2021), and the nominal Reynolds number and the effective Reynolds number are two. The nominal Reynolds number reads

(4.1)\begin{equation} Re_{\tau ,N} = \frac{{{u_{\tau ,N}}(\delta + h)}}{\nu}, \end{equation}

where $u_{\tau,N}=\sqrt {(\delta +h)(-{{\rm d}}P/{{\rm d}\kern0.7pt x})/\rho }$ is the nominal friction velocity, and $\nu$ is the kinematic viscosity. The effective Reynolds number reads

(4.2)\begin{equation} Re_{\tau ,E} = \frac{{{u_{\tau ,E}}\delta }}{\nu }, \end{equation}

where $u_{\tau,E}=\sqrt {\delta (-{{\rm d}}P/{{\rm d}\kern0.7pt x})/\rho }$ is the effective friction velocity. Since the canopy occupies a substantial portion of the computational domain, the nominal Reynolds number will be larger than the effective Reynolds number. Considering that the boundary layer interacts with only the top part of the canopy, the effective Reynolds number is a more appropriate measure of the ratio between the largest and the smallest scales in the flow – although the computational cost scales with the nominal Reynolds number. Three Reynolds numbers are considered, namely $Re_{\tau,E}=280$, 367 and 560. These effective Reynolds numbers correspond to $Re_{\tau,N}=790$, and 1580.

Table 1 shows the details of our DNSs, where we tabulate the effective and the nominal Reynolds numbers, the domain sizes, the canopy and the boundary-layer heights, the grid numbers in the three Cartesian directions and the number of rectangular prisms in the streamwise and the spanwise directions. The nomenclature of the cases is as follows: $R[Re_{\tau,E}/100]L[100\lambda _p][{A/S}]$, where $A$ is for ‘aligned’ and $S$ is for ‘staggered’. Here and throughout the paper, the superscript ‘$+$’ denotes normalization by $u_{\tau,E}$ and $\nu /{u_{\tau,E}}$.

Table 1. DNS details. Here, $L_x$, $L_y$ and $L_z$ are the domain sizes in the streamwise, spanwise and vertical directions, respectively; ${w_r}$ is the width of the prisms; $h$ is the canopy height; $\delta$ is the height of the boundary layer above the canopy; $N_x$, $N_y$ and $N_z$ are the numbers of the grid points in the streamwise, spanwise and vertical directions, respectively; $n_x$ and $n_y$ are the numbers of prisms in the $x$ and $y$ directions. The nomenclature of the cases is as follows: $R[Re_{\tau,E}/100]L[100\lambda _p][{A/S}]$, where A is for ‘aligned’ and S is for ‘staggered’.

Next, we explain our choices of domain size and grid resolution. The boundary layer interacts with only the top part of the canopy, and therefore the height of the boundary layer above the canopy is a more appropriate measure of the domain size than the half-channel height. In our DNSs, the spanwise extents of computational domains are such that $L_y\ge 4\delta$. The streamwise domain sizes are slightly different for flows at different Reynolds numbers. We use a long streamwise domain for the R3 and R5 DNSs, i.e. $L_x=8\delta$ and $9\delta$, respectively. A shorter streamwise domain is used for the R2 DNSs, i.e. $L_x=4\delta$. In the following, we compare our domains with those in the existing literature. Lozano-Durán & Jiménez (2014) studied the effects of domain size on flow statistics and concluded that the streamwise and the spanwise domain sizes must be $L_x>2{\rm \pi} \delta$ and $L_y>3\delta$. The authors argued that a smaller domain might affect the first- and second-order statistics. Lozano-Durán & Jiménez (2014) considered plane channel flow. The requirement for canopy flow, or rough-wall channel flow, is different and is less stringent. Leonardi & Castro (2010) used a domain of size $L_x\times L_y\times L_z=8{w_r}\times 6{w_r}\times 6{w_r}$. Coceal et al. (2006) studied the effects of domain size on the flow over cube arrays, and concluded that $L_x\times L_y\times L_z=4{w_r}\times 4{w_r}\times 4{w_r}$ is sufficient if one is interested in only low-order flow statistics. These two domains translate to $L_x\times L_y\times L_z=1.6\delta \times 1.2\delta \times \delta$ and $L_x\times L_y\times L_z=1.33\delta \times 1.33\delta \times \delta$, both of which are notably smaller than the minimum channel in Lozano-Durán & Jiménez (2014). In the more recent work, Chung et al. (2015) and MacDonald et al. (2017) concluded that, if the purpose is to study the roughness’ aerodynamic properties (which is our purpose), one may well employ a minimal-span channel whose size is $L_y \ge \max (100\nu /u_{\tau },k/0.4,\lambda _{r,y})$, where $k$ is the (effective) roughness height, $\lambda _{r,y}$ is the spanwise length scale of roughness elements. In light of the discussion above and the recent DNS studies (Sharma & Garcia-Mayoral 2020), we may conclude that our domains are larger than needed. The purpose of employing a large domain for our higher Reynolds number DNSs is to resolve the K–H waves – if they exist. We have also conducted a DNS where we double the streamwise domain size for R2L11A. The case is referred to as R2L11A(L). In figure 8(a), we compare the mean flows in R2L11A and R2L11A(L), and there is barely any difference. Here and in later sections, all statistics are streamwise and spanwise averaged and temporally averaged for 20 large-eddy turnover times. We can confirm a flow's statistical convergence by examining the momentum budget, which we will do in § 6.4.

Figure 8. The mean velocity profile in R2L11A compared with the results (a) in a larger domain: R2L11A(L), (b) with two coarsened resolutions: R2L11A(M), R2L11A(C) and (c) with two refined resolutions: R2L11A(F1), R2L11A(F2).

In addition, the vertical domain must be sufficiently high so that the top boundary does not affect the flow near the canopy crest. This usually requires the height of the domain above the canopy, denoted as $\delta$, to be sufficiently larger than the roughness’ characteristic length scale. For deep canopies, the canopy height, denoted as $h$, is not the roughness’ characteristic length scale because the lower part of the canopy does not directly interact with the flow above the canopy. A good measure of the deep canopy's characteristic length scale is the attenuation length scale $L_m$. For the cases studied in this work, $L_m/\delta \lesssim 0.2$, and therefore the top boundary should not affect the flow near the canopy crest. Further evidence will be presented in § 6.3.

Last, we explain our choice of the grid. We use a uniform grid spacing in the streamwise and the spanwise direction and slightly stretch the grid in the wall-normal directions away from the top of the canopy and the ground, where the mean velocity gradient is large. The grid resolution is such that ${\rm \Delta} x^{+}={\rm \Delta} y^{+}<4.5$ and ${\rm \Delta} z^{+}\approx 2.5$. The wall unit is defined using the kinematic viscosity and the nominal friction velocity $u_\tau =\sqrt {D/\rho }$, where $D$ is the drag force per unit area. In a turbulent rough-wall boundary layer, most drag is form drag, and therefore the defined wall unit is an underestimate (and therefore overly conservative estimate) of the viscous scale at the top of the roughness. A more physical measure of the viscous length scale is the Kolmogorov length scale. Figure 9 shows the grid resolution measured by the Kolmogorov length scale. We see that ${\rm \Delta} x/\eta ={\rm \Delta} y/\eta$ is around 1.5 at most places and nowhere exceeds 2.5, and ${\rm \Delta} z/\eta$ is around 1 at most places and nowhere exceeds 2. This grid resolution is quite high, and we see no apparent Gibbs oscillations in our DNSs. The non-physical Gibbs oscillation has a less significant impact on the simulations with DNS or wall-resolved large-eddy simulation (LES) rather than WMLES, since the velocity in DNS and wall-resolved LES transitions more gradually to the no-slip condition near the wall (Li, Bou-Zeid & Anderson 2016).

Figure 9. Ratio of the grid spacing to the Kolmogorov scale: (a–c) ${\rm \Delta} x/\eta$, (d–f) ${\rm \Delta} z/\eta$, and (a,d) R2L06/11/25/39A/S, (b,e) R3L06/11/25/39A/S, (c,f) R5L11/25A/S. The solid/dashed lines represent the aligned/staggered configurations respectively.

We also conduct a grid convergence study. Specifically, we coarsen R2L11A's grid by a factor of 1.25 and 1.67 in the horizontal directions. The two cases are R2L11A(M) and R2L11A(C), where M stands for medium and C stands for coarse. Furthermore, we refine R2L11A's grid by a factor 1.5 in the streamwise and the spanwise direction, and stretch the wall-normal grid such that the wall-normal grid resolution at the roughness crest is refined by a factor of 1.25 and 3. These two cases are R2L11A(F1) and R2L11A(F2). We compare the mean flow results in these five calculations in figure 8(b,c). R2L11A(C) yields a slightly lower $U$ near the top boundary, but the profiles in R2L11A and R2L11A(M) collapse. The R2L11A(F1), R2L11A(F2) and the R2L11A results collapse. These results suggest that the standard grid resolution, i.e. the grid resolution used for R2L11A, is sufficiently fine. It is, however, worth noting that this conclusion applies only to low-order statistics, high-order statistics may well require finer grid resolution (Yang et al. 2021). We now compare our grid resolution with those in the literature. The typical grid resolution requirement for flat-plate boundary-layer flow is ${\rm \Delta} x^{+} \approx 12$, ${\rm \Delta} y^{+} \approx 6$, and ${\rm \Delta} z_{max}^{+} \approx 8$ (Kim, Moin & Moser 1987; Lee et al. 2011; Lozano-Durán & Jiménez 2014; Xu et al. 2021). The grid resolution requirement for rough-wall boundary layers is similar. The resolution is ${\rm \Delta} x^{+}={\rm \Delta} y^{+}=19$ in Leonardi & Castro (2010), ${\rm \Delta} ^{+}=$7.8 to 15.6 in Coceal et al. (2006) and ${\rm \Delta} z_{max}^{+}=8$ to 13 in MacDonald et al. (2018). Compared with these previous studies, our grid resolutions are excessive. The purpose of employing such a fine grid resolution is to fully resolve the thin shear layer at the top of the canopy.

5. Evidence for the mixing-layer analogy in deep canopy flows

The basis of Raupach et al.'s mixing-layer analogy is the inflected velocity profile. Figure 10 shows the time- and plane-averaged streamwise velocity in our DNSs. We see that the velocity profiles are inflected at the canopy crest. Although this is trivial, the results in figure 10 is the first evidence for the mixing-layer analogy. In addition to the fact that velocity profiles are inflected, we make the following observations. Firstly, an inactive layer emerges below the active layer (see figure 6 for the definition of the inactive layer), where the mean velocity is approximately constant and equals $U_c$. The fact that an inactive layer emerges shows that the canopies are deep (Poggi et al. 2004; Nepf 2012; Nikora, Nikora & O'Donoghue 2013). Secondly, as more roughness leads to more resistance in the canopy layer, $U_c$ decreases as the surface coverage density increases for a given roughness arrangement. Also, as the staggered arrangement leads to more blockage than the aligned arrangement, the staggered arrangement results in a smaller $U_c$ than the aligned arrangement for a given surface coverage density. Thirdly, the reader may notice that the velocity at the top of the half-channel, $U(\delta )/U_h$, is larger in the staggered cases than in the aligned cases. This counter-intuitive result is due to smaller $U_h^{+}$ in the staggered cases: although $U(\delta )^{+}$ is also smaller in the staggered cases, an even smaller $U_h^{+}$ results in a larger $U(\delta )/U_h$ in the staggered cases.

Figure 10. Time- and plane-averaged streamwise velocity profiles. (a) R2 cases, (b) R3 cases, (c) R5 cases. The solid and dashed lines are the aligned and staggered configuration results, respectively. The ranges of the $x$ and $y$ axes are kept the same among the three figures.

In addition to the inflected velocity profile, (2.8), if confirmed, is evidence for the mixing-layer analogy. In the following, we compare the velocity profile in the canopy layer with (2.8). Figure 11 shows the double-averaged velocity in the canopy layer. The velocity is scaled according to $(U-U_c)/(U_h-U_c)$, where $U_c$ is the velocity in the inactive layer and is readily available in the DNSs. We see that $\log ((U-U_c)/(U_h-U_c))\sim z/{w_r}$ in the active layer, which is encouraging. In addition, we observe the following. Firstly, the staggered arrangement results in more attenuated velocity profiles than the aligned arrangement for a given surface coverage density. Secondly, the velocity overshoots near the wall. This overshoot is quite peculiar and does not seem to have been reported by any other authors before. However, since this paper does not concern the ground layer, we will leave this overshoot for future investigation.

Figure 11. Mean velocity profiles in the canopy layer. (a) R2 case, (b) R3 cases, (c) R5 cases. Again, the solid and dashed lines represent the aligned and staggered cases.

Next, we fit for $L_m$ and $d_c$ in the active layer. Figure 12 shows the scaled velocity as a function of $(z+d_c)/L_m$. We see that the mean velocity follows the exponential profile closely from approximately $z+d_c=-2L_m$ to $z+d_c=-0.4L_m$. The extent of the active layer varies from one case to another as the sizes of the viscous layer and the ground layer vary from one case to another. The viscous layer occupies more space at higher packing densities, and the ground layer is thicker at lower packing densities. It is also worth noting that the Reynolds number does not play an important role here, and (3.2) and (3.6) are valid irrespective of the Reynolds number. It is worth noting that whether DNS results are Reynolds number independent is not a concern to us. Our concern is whether (3.2) and (3.6) are valid at all (turbulent) Reynolds numbers. The results in figure 12 is the second evidence for the mixing-layer analogy.

Figure 12. Mean velocity profiles in the canopy layer. The velocity is scaled according to (2.8). The black solid line represents (2.8). ‘Exp’ in the legend stands for ‘exponential profile’.

Figure 13 compares the data with (3.6). The range of $(U_h/U_c-1)(U_h/U_c+3)$ is more than a decade. A higher surface coverage density leads to a larger $U_h/U_c$, which, in turn, leads to a larger $(U_h/U_c-1)(U_h/U_c+3)$. For a given surface coverage density, the staggered arrangement leads to a larger $U_h/U_c$ and a larger $(U_h/U_c-1)(U_h/U_c+3)$ than the aligned arrangement. The data follow (3.6) closely for $(U_h/U_c-1)(U_h/U_c+3)>6$, and the result is the third piece of evidence for the mixing-layer analogy. Disagreements are found for $(U_h/U_c-1)(U_h/U_c+3)<5$, specifically, for $R2L06A/S$ and $R3L06A/S$. In these L06 cases, the active layer intrudes to the wall, and the deep canopy assumption breaks.

Figure 13. The attenuation length scale as a function of $(U_h/U_c-1)(U_h/U_c+3)$. We use different colours for different Reynolds numbers, i.e. blue for R2, red for R3 and yellow for R5. The open and solid symbols represent aligned and staggered roughness arrangements, with squares for L06A/S, circles for L11A/S, left triangle for L25A/S and right triangle for L39A/S.

6. Discussion

Having presented evidence for the mixing-layer analogy, we now scrutinize some of the claims in the previous sections. In § 6.1, we will explain why (1.3) fails. In § 6.2, we discuss a few extensions of the derivation in § 3. We will show that the canopies in our DNSs are deep in § 6.3 and that our data are statistically converged in § 6.4 . Our results are compared with those in the vegetation canopy flow literature in § 6.5. Last, we will report the roughness properties in § 6.6.

6.1. Revisiting the methodology in Raupach et al. (1996)

We follow Raupach et al. (1996) and compare our data with (1.3). The results will add to the discussion in § 1 and serve as a motivation for the present work.

Raupach et al. (1996) argued that the vertical velocity is less susceptible to inactive eddies and therefore is more indicative of the active eddies in the roughness sublayer than the streamwise component. We follow the previous authors and examine the vertical velocity. Figures 14 and 15 show the contours of the instantaneous vertical velocity fluctuations at $z^{+}=5$ and $z^{+}=15$, respectively, in the R3 cases. The results in the R2 and R5 cases are similar and are not shown here for brevity. Both $z^{+}=5$ and 15 are within the roughness layer, but the wakes behind the individual roughness elements are less discernible at $z^{+}=15$ than at $z^{+}=5$. To prevent any interplay between the wake flow and the K–H type waves (if they exist), we focus on the data at $z^{+}=15$. There, we see no clear evidence for K–H-type instability, i.e. no K–H rollers.

Figure 14. Contours of the instantaneous vertical velocity fluctuation at $z^{+}=5$. The shear length scales are indicated above each panel. (a) R3L06A, (b) R3L11A, (c) R3L25A, (d) R3L39A.

Figure 15. Same as figure 14 but for $z^{+}=15$.

Next, we follow Raupach et al. (1996) and examine the velocity correlation. The two-point correlation reads

(6.1)\begin{equation} {R_{\phi_1\phi_2}(z,{z_{ref}},{r_x})} = \frac{{\langle {\phi_1' (x,y,{z_{ref}},t)\phi_2'(x + {r_x},y,z,t)} \rangle }}{{\sqrt {\langle {\phi_1'{{(x,y,{z_{ref}},t)}^{2}}} \rangle \langle {\phi_2'{{(x,y,z,t)}^{2}}} \rangle }}}, \end{equation}

where $z_{ref}$ is the reference vertical location, $r_x$ is the streamwise separation distance, $\phi '$ is the velocity fluctuation in time. Here, $z_{ref}^{+}=15$. Figures 16 and 17 show the two-point correlation of the streamwise and the vertical velocities, i.e. $R_{uu}$ and $R_{ww}$, in the R3 cases. The statistical object $R_{uu}$ is included here for comparison purposes. We see from figure 16 that the contour lines are inclined towards the positive $x$ direction (Adrian, Meinhart & Tomkins 2000; Heisel et al. 2018; Zhang et al. 2021b). Raupach et al. (1996) measured the spacing between two adjacent eddies in the roughness sublayer according to

(6.2)\begin{equation} {\varLambda _x} = 2{\rm \pi} \alpha {L_w}. \end{equation}

Here, $\alpha =1.8$ is a constant, and $L_w$ is the correlation length scale computed as follows:

(6.3)\begin{equation} L_{w}(z_{ref},z) = \int_0^{\infty} {R_{ww}({r_x},{z_{ref}},z)} \,{{\rm d}}{r_x}. \end{equation}

Because of the finite computational domain, $r_x$ is limited to $L_x/2$, and we employ the following equation for $L_w$:

(6.4)\begin{equation} L_{w}(z_{ref}^{+},z) = \frac{1}{2}\int_{{-}L_x/2}^{L_x/2} {R_{ww}({r_x},{z_{ref}^{+}},z)} \, {{\rm d}} {r_x}, \end{equation}

again, $z_{ref}^{+}=15$. Figure 18 shows the correlation length scale $\varLambda _x$ as a function of the shear length scale $L_s$. We see that the data do not support (1.3) or the presence of K–H rollers. We also examined the data reported in Coceal et al. (2007), Kanda et al. (2004), Yang et al. (2019) and Xu et al. (2021), but could not find K–H rollers. Huq et al. (2007) carried out extensive visualization attempts in experiments, and they could not find K–H rollers either. The above is consistent with the existing literature: there are no K–H rollers in urban-canopy flows. Although we cannot exclude the existence of K–H rollers in deep canopy flows, the chance that K–H rollers can be found in deep urban-canopy flows is slim. Following this line of thinking, it would be even harder to find K–H rollers in shallow urban canopies and above closely packed buildings of the same height. The mixing-layer analogy is not meant for the former, and the surface layer analogy is more appropriate for the latter.

Figure 16. Contours of the two-point correlation of the streamwise velocity ($R_{uu}$) in streamwise/vertical plane. The reference location is at $z_{ref}^{+}=15$, i.e. slightly above the roughness canopy. (a) R3L06A, (b) R3L11A, (c) R3L25A, (d) R3L39A.

Figure 17. Same as figure 16 but for the two-point correlation of the vertical velocity ($R_{ww}$).

Figure 18. The correlation length scale $\varLambda _x$ as a function of the shear length scale $L_s$. The solid line represents $\varLambda _x=8.1L_s$, while the horizontal dashed line represents $\varLambda _x/{w_r}=0.7$. The symbols are results obtained from R3 cases (square: R3L06A, circle: R3L11A, left triangle: R3L25A, right triangle: R3L39A).

We remark on the results in this subsection. Firstly, the failure of (1.3) in deep urban-canopy flows and the success of the mixing-layer analogy in deep vegetation canopy flows motivated this work. Secondly, the failure of (1.3) does not invalidate the mixing-layer analogy. Raupach et al. proposed the mixing-layer analogy as an alternative to the surface layer analogy. The analogy has value because it gives us more insights into the flow in the roughness sublayer than the surface layer analogy. While (1.3) has been the cornerstone of the mixing-layer analogy, the analogy is more than (1.3) (Brunet 2020). Thirdly, the failure of (1.3) does not invalidate the arguments that led to it. Raupach et al. argued that the inflected velocity profile would give rise to K–H instability, which then results in (1.3). The velocity profile is an inflected one in urban-canopy flows, and K–H instability surely plays a role in urban-canopy flows (Brunet 2020). Equation (1.3) fails because K–H instability is not the dominant mechanism in deep urban-canopy flows. Fourthly, the failure of (1.3) does not invalidate the shear length scale $L_s$ or the correlation length scale $\varLambda _x$ as important length scales in the roughness sublayer. The correlation length scale is a measure of the largest eddy in the flow and is undoubtedly an important length scale. Figure 19 compares the shear length scale $L_s$ and the attenuation length scale $L_m$. We see that the two are well correlated.

Figure 19. The shear length scale $L_s$ as a function of the attenuation length scale $L_m$. We use different colours for different Reynolds numbers: blue for R2, red for R3 and yellow for R5. The open and solid symbols represent the aligned and staggered configurations. Different symbols represent different surface coverage densities. Here, squares: L06A/S, circles: L11A/S, left triangles: L25A/S, right triangles: L39A/S. The solid and dashed lines represent $L_s=0.7L_m\pm 0.05\delta$.

6.2. Further discussion on the formulation in § 3

We discuss two extensions of the derivation in § 3.

Firstly, we consider the limit $U_c\to 0$. This limit is attained if the buildings are very closely packed or if the boundary-layer flow is not driven by a pressure gradient. At this limit, the velocity in the canopy layer is $U/{U_h} = \exp (z/{L_m})$, i.e. (2.4), and the force balance reads

(6.5)\begin{equation} \int_{ - h}^{0} {\rho {C_d}U{{(z)}^{2}}\, {{\rm d}}z} = \rho {u_\tau }^{2}, \end{equation}

with $u_\tau =\sqrt {\tau _w/\rho }$ by definition. Plugging (2.4) into (6.5) leads to

(6.6)\begin{equation} {C_d}{U_h}^{2}\frac{{{L_m}}}{2}[ {1 - \exp ( - 2h/{L_m})} ] = {u_\tau }^{2}. \end{equation}

In a deep canopy, we have $h \gg L_m$, $\exp (-h/L_m)\approx 0$, and therefore (6.6) becomes

(6.7)\begin{equation} {L_m} \approx \frac{2}{{{C_d}{{( {{U_h}/{u_\tau }} )}^{2}}}}, \end{equation}

which suggests that $L_m$ is inversely proportional to ${( {{U_h}/{u_\tau }} )^{2}}$.

Secondly, we discuss the case if we had plugged (2.8) instead of (2.6) into (3.3). Plugging (2.8) into (3.3), the integration becomes

(6.8)\begin{align} &{\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right)^{2}} \frac{{{L_m}}}{2}\exp \left({\frac{{2{d_c}}}{{{L_m}}}} \right) \left[ {1 - \exp \left( { - \frac{{2h}}{{{L_m}}}} \right)} \right] \nonumber\\ &\quad + 2\left( {\frac{{{U_h}}}{{{U_c}}} - 1} \right){L_m}\exp \left( {\frac{{{d_c}}}{{{L_m}}}} \right)\left[ {1 - \exp \left( { - \frac{h}{{{L_m}}}} \right)} \right] = \delta. \end{align}

Here, $d_c$ measures the size of the viscous layer at the top of the canopy and is much smaller than $L_m$. It follows that $\exp ({d_c}/{L_m})\approx 1$ and (6.8) reduces to (3.4).

6.3. Canopy depth

A canopy is considered deep if the bottom wall does not affect the flow near and above the canopy crest.

Consider R2L11A(h1), R2L11A(h2), R2L1A (the last case is the standard case where $h/{w_r}=4$) and R2L11A(h6), where we vary the canopy depth from $h={w_r}$ to $h=6{w_r}$. Figure 20(a,b) shows the comprehensive/superficial average (Mignot et al. 2009; Xie & Fuka 2018; Schmid et al. 2019) of the velocity within and above the canopy. The R2L11A and the R2L11A(h6) results collapse in figure 20(a). All four profiles collapse near the canopy crest in figure 20(b). In addition to the velocity profiles, we may examine the roughness properties. Figure 21 shows the roughness function ${\rm \Delta} U^{+}$ and the location of the virtual wall $d/w$, as a function of the canopy depth. The general expectation is as follows (MacDonald et al. 2018). When the canopy is shallow, increasing the canopy depth increases the drag, which in turn increases the roughness function. Meanwhile, increasing the canopy depth pushes the virtual ground away from the canopy crest, and the roughness transitions from k-type to d-type (Jiménez 2004). Once the roughness is d-type, further increasing the canopy depth has no effect on the roughness function nor the location of the virtual ground as the boundary layer only interacts with the top part of the roughness. We see that ${\rm \Delta} U^{+}$ and $d$ increases as $h$ increases and stay a constant as $h$ increases beyond $h=4w$. Hence, we can conclude that the canopy in R2L11A is a deep one.

Figure 20. The time- and plane-averaged streamwise velocity profile (a) above the canopy and (b) within the canopy for $h/{w_r}=1, 2, 4$ and 6.

Figure 21. (a) The roughness function and (b) the location of the virtual wall in R2L11A and R2L11A(h1/h2/h6). We will explain how we measure ${\rm \Delta} U^{+}$ and $d$ in § 6.6.

The canopy in R2L11A is a deep one. This readily guarantees that the canopies are deep in the L11S cases (for R2, R3 and R5), the L25 cases (for R2, R3 and R5, for S and A) and the L39 cases (for R2, R3 and R5, for S and A) because of the following. A canopy is deep if $h\gtrsim L_m$; $L_m$ is a decreasing function of the surface coverage density, is larger for aligned roughness than for staggered roughness and is insensitive to the Reynolds number (at least for the Reynolds numbers investigated here). Hence, the fact that the canopy in R2L11A is deep ensures the canopies with higher surface coverage densities, higher Reynolds numbers and staggered roughness arrangement are all deep canopies. Nonetheless, it is worth noting that the deep canopy in R2L11A does not guarantee the canopies in the L06 cases are deep, which may well explain the discrepancies between the L06 data and our theory in figure 13.

6.4. Statistical convergence

The Reynolds-averaged $x$-momentum equation reads

(6.9)\begin{equation} \frac{{\rm d}}{{{\rm d}}z}\left(\nu\frac{{\rm d} \langle {U} \rangle}{{{\rm d}}z}-\langle \overline{u'v'} \rangle-\langle u''v'' \rangle\right)-\frac{1}{\rho} \frac{{\rm d} \langle \bar{p} \rangle}{{\rm d}\kern0.7pt x}-f=0, \end{equation}

where $f$ is the drag force (which is 0 outside the roughness occupied region), ${\rm d}/{{\rm d}\kern0.7pt x}$ and $d/dz$ are total derivative in the $x$ and $z$ directions, respectively (streamwise- and spanwise-averaged velocity and stresses are only functions of $z$), $\nu$ is the kinematic viscosity, $p$ is the pressure, $\bar {{\cdot }}$ denotes time average, $\langle {\cdot } \rangle$ denotes spatial average, $\phi '=\phi - \bar {\phi }$ is the turbulent fluctuation, $\phi ''=\bar {\phi }-\langle \bar {\phi } \rangle$ is the spatial variation and $\phi$ is a generic variable (note that $\overline {\phi ''_i\phi ''_j}\equiv \phi ''_i\phi ''_j$ for any $\phi _i$ and $\phi _j$). The terms on the left-hand side are the viscous diffusion term, the turbulent transport term, the dispersive stress term and the pressure-gradient term. Integrating (6.9) in the $z$ direction leads to

(6.10)\begin{equation} \nu\frac{{{\rm d}} \langle {U} \rangle}{{\rm d} y}-\langle \overline{u'v'} \rangle-\langle u''v'' \rangle={\rm Const.}+\frac{1}{\rho}\frac{{{\rm d}} \langle \bar{p} \rangle}{{\rm d}\kern0.7pt x} z, \end{equation}

outside the canopy. The total stress is a linear function of $z$ if the flow is statistically converged. Equation (6.10) is often used to check the statistical convergence of a numerical simulation: a simulation is statistically converged if the sum of the viscous, turbulent and dispersive fluxes is a linear function of $z$ (Oliver et al. 2014). Figure 22 shows the terms in (6.10) in the R3 cases. The results in the R2 and R5 cases are similar and are not shown here for brevity. The total stress is a linear function of $z$, and therefore the flow is statistically converged.

Figure 22. The terms in the Reynolds-averaged $x$-momentum equation. From the top to the bottom: $\lambda _p=0.06$, $0.11$, $0.25$ and $0.39$. Left/right: aligned/staggered. The solid line represents $1-z/\delta$.

In addition to the linear total stress, we observe the following. Firstly, the turbulent stress is the largest term outside the canopy. The dispersive stress decays to 0 slightly above the canopy, suggesting a thin roughness sublayer above the canopy. The viscous stress peaks at the canopy crest as a result of the inflection in the mean velocity profile, but is very small everywhere else. The above are consistent with previous DNS studies on cuboidal roughness (Leonardi & Castro 2010; Xu et al. 2021). In the canopy layer, the turbulent stress quickly decays away from the canopy crest, and the viscous stress is nearly negligible. Secondly, the dispersive stresses are notable in the L06 and the L11 cases but are much smaller in the L25 and the L39 cases. A non-zero dispersive stress indicates a non-homogeneous mean flow field. It certainly makes sense that the mean flow is inhomogeneous within the canopy. The fact that the dispersive stress is small within the canopy in the L25 and the L39 cases suggests that the flow is suppressed by the presence of the canopy roughness.

6.5. Further analysis

The mixing-layer analogy has received much attention in vegetation canopy flow research since Raupach et al.'s (1996) pioneering work. While it is worthwhile to repeat the analysis in Raupach et al. (1996) and re-examine (1.3), it is unreasonable to repeat all previous analyses and re-examine all previous conclusions in this work, and we have to limit the discussion to a few historically important quantities. In addition to the correlation length scale $\varLambda _x$, which we examine in § 6.1, and the time- and plane-averaged velocity, which we examine in § 5, $u'w'$ is a historically important flow quantity (as far as the mixing-layer analogy is concerned). The streamwise velocity fluctuation $u'$ is negatively correlated with the vertical velocity fluctuation $w'$ in a surface layer flow. The fact that $u'$ and $w'$ are negatively correlated in the roughness sublayer was often cited as the evidence for the surface layer analogy, since a first-order turbulent diffusivity can successfully predict the mean velocity profile in a similar way to the surface layer, see Brunet (2020) and the references therein. The surface layer analogy was refuted in the 1980s and replaced with the mixing-layer analogy in the 1990s, and since then more attention has been given to how $u'w'$ behaves differently in a surface layer and a roughness sublayer.

The following two quantities have been given considerable attention: S4/S2 and T4/T2, where S4 is the contribution of Q4 events (positive $u'$ negative $w'$) to the Reynolds stress $\langle u'w' \rangle$, S2 is the contribution of Q2 events (negative $u'$ and positive $w'$) to the Reynolds stress $\langle u'w' \rangle$, T4 and T2 are the time fractions of Q4 and Q2 events, respectively. Novak et al. (2000), Finnigan et al. (2009), Gavrilov et al. (2011), Marjoribanks et al. (2017), among others, considered vegetation canopies and reported that S4/S2 is larger than 1 within the canopy and close to 1 above the canopy, and that T4/T2 is smaller than 1 within the canopy and greater than 1 above the canopy. Figure 23 shows T4/T2 and S4/S2 in the R3 cases. Results in the R2 and R5 cases are similar and are not reported here for brevity. We see very similar results in urban-canopy flows as in vegetation canopy flows (Castro, Cheng & Reynolds 2006; Kanda 2006; Coceal et al. 2007; Inagaki et al. 2012; Anderson et al. 2015b; Li & Bou-Zeid 2019), where the Q4 events are found to be infrequent but strong.

Figure 23. Ratio of Q4 and Q2 events: (a,b) time fraction $T_4/T_2$ and (c,d) stress contribution $S_4/S_2$ for R3 cases with (a,c) aligned and (b,d) staggered configurations.

6.6. Roughness properties

Last, we report the roughness properties and examine the velocity profiles for the presence of a logarithmic region. Table 2 lists the roughness function and the location of the virtual wall for all DNSs. We fit the following logarithmic scaling for these roughness properties:

(6.11)\begin{equation} {U^{+} } = {\kappa^{ - 1}}\ln [{z^{+} } + {d^{+} }] + B - {\rm \Delta} {U^{+} }, \end{equation}

where $\kappa =0.4$ is the Kármán constant, $d$ is the zero-plane displacement, $B=5$ is the smooth-wall constant and ${\rm \Delta} U^{+}$ is the roughness function. We can also write the log law as follows:

(6.12)\begin{equation} {U^{+} } = \frac{1}{\kappa }\ln \frac{{z + d}}{{{z_0}}}, \end{equation}

where ${z_0} = ( {\nu /{u_\tau }} )\exp [ {\kappa ( { - B + {\rm \Delta} {U^{+} }} )} ]$. The exact procedure is as follows. We first determine the zero-plane displacement by examining the following diagnostic function (Lozano-Durán & Jiménez 2014):

(6.13)\begin{equation} {\varXi^{+} } = \kappa {(z^{+} + d^{+}) }\frac{{{\rm d}{U^{+} }}}{{{\rm d}{z^{+} }}}, \end{equation}

which should be equal to 1 in the logarithmic region at $z^{+} \approx 60$. Again, the $z$ coordinate is such that $z=0$ at the canopy top, and $d$ measures the distance from the canopy top to the virtual ground. We then fit the logarithmic law of the wall for the roughness function ${\rm \Delta} U^{+}$. Here, the location of the virtual origin is not computed according to Jackson (1981) as the location of the effective drag force, which is known to underestimate the displacement height for deep canopy flows (Huq & Rahman 2018; Xu et al. 2021). Figure 24 shows the scaled double-averaged velocity profiles. The results are very much as expected. A logarithmic layer emerges as we scale the wall-normal coordinate according to $(z+d)^{+}$ and the mean velocity using the friction velocity $u_\tau$. Except for the L39 case, the surface coverage densities investigated here are low, and the roughness function increases as the surface coverage density increases (Jiménez 2004). Furthermore, the staggered arrangement leads to a larger roughness function than the aligned arrangement.

Table 2. The zero-plane displacement and the roughness function of all DNS cases.

Figure 24. Double-averaged velocity profiles. (a) R2 cases, (b) R3 cases, (c) R5 cases. The solid black line is the logarithmic law of the wall, i.e. $U^{+}=\ln (y^{+})/\kappa +B$ with $\kappa =0.4$ and $B=5$.

The attenuation length scale $L_m$ plays a central role in the mixing-layer analogy. Here, it is our interest to see if the zero-plane displacement height and the roughness function are related to the attenuation length scale $L_m$. Figure 25 shows the zero-plane displacement height, measured from the canopy crest, the roughness function and the equivalent roughness height $z_0$ as a function of the attenuation length scale. The three quantities, i.e. $d/{w_r}$, ${\rm \Delta} U^{+}$ and $z_0$, are generally increasing functions of $L_m/\delta$ with most data lying within $d/w_r=5L_m/\delta \pm 0.3$, ${\rm \Delta} U^{+}=6.5L_m/\delta +10\pm 2$ and $z_0/\delta =0.25L_m/\delta \pm 0.015$.

Figure 25. (a) The zero-plane displacement ($d/w_r$), (b) the roughness function (${\rm \Delta} U^{+}$), (c) the roughness length scale ($z_0$) as a function of the attenuation length scale $L_m$. The solid and dashed lines are $d/w_r=5L_m/\delta \pm 0.3$ in (a), ${\rm \Delta} U^{+}=6.5L_m/\delta +10\pm 2$ in (b) and $z_0/\delta =0.25L_m/\delta \pm 0.015$ in (c).

7. Conclusions

Raupach et al. (1989, 1996) argued that the roughness sublayer is more analogous to a plane mixing layer rather than a surface layer. This argument led to the mixing-layer analogy. This work shows that the following are evidence for Raupach et al.'s analogy: first, a velocity profile that is inflected, second, a velocity profile that follows an exponential scaling in the active layer, i.e.

(7.1)\begin{equation} \frac{U-U_c}{U_h-U_c}=\exp\left(\frac{z+d_c}{L_m}\right), \end{equation}

and third, an attenuation length scale that follows

(7.2)\begin{equation} \frac{L_m}{\delta}=\frac{2}{({U_h/U_c-1})({U_h/U_c+3})}. \end{equation}

Here, $U$ is the double-averaged velocity, $U_c$ is the fluid velocity in the inactive layer and is non-zero when the flow is pressure driven, $U_h$ is the fluid velocity at the canopy crest, $z$ is the vertical coordinate and is such that $z=0$ at the canopy crest, $d_c$ is a displacement height and accounts for the viscous layer at the canopy crest and $L_m$ is the attenuation length scale.

In order to test the mixing-layer analogy, we conduct DNSs. We consider deep urban canopies with surface coverage densities between $\lambda =0.06$ and $0.39$, aligned and staggered roughness arrangements and effective Reynolds numbers $Re_{\tau, E}$ between 180 and 550. The data support the mixing-layer analogy.

In the light of the supporting evidence, a word of caution must be mentioned. The mixing-layer analogy has its applicable range and cannot be taken too literally (Brunet 2020). Even limited to deep homogeneous vegetation canopy, it is rather clear that the roughness sublayer and the plane mixing layer are two different flows. For one, the flows inside and above the canopy are asymmetric, whereas the plane mixing layer is centro-symmetric. For another, the roughness sublayer is bounded on one side and unbounded on the other. Historically, Raupach et al. proposed the mixing-layer analogy as an alternative to the surface layer analogy. The analogy has value because it is more reasonable than the surface layer analogy and because it gives a familiar baseline. The same is true in the context of deep urban-canopy flows.