The structure of turbulence in unsteady flow over urban canopies

Abstract

The topology of turbulent coherent structures is known to regulate the transport of energy, mass and momentum in the atmospheric boundary layer (ABL). While previous research has primarily focused on characterizing the structure of turbulence in stationary ABL flows, real-world scenarios frequently deviate from stationarity, giving rise to nuanced and poorly understood changes in the turbulence geometry and associated transport mechanisms. This study sheds light on this problem by examining topological changes in ABL turbulence induced by non-stationarity and their effects on momentum transport. Results from a large-eddy simulation of pulsatile open channel flow over an array of surface-mounted cuboids are examined. The analysis reveals that the flow pulsation triggers a phase-dependent shear rate, and the ejection-sweep pattern varies with the shear rate during the pulsatile cycle. From a turbulence structure perspective, it is attributed to the changes in the geometry of hairpin vortices. An increase (decrease) in the shear rate intensifies (relaxes) these structures, leading to an increase (decrease) in the frequency of ejections and an amplification (reduction) of their percentage contribution to the total momentum flux. Furthermore, the size of the hairpin packets undergoes variations, which depend on the geometry of the constituting hairpin vortices, yet the packet inclination preserves its orientation throughout the pulsatile cycle. These observations reinforce the important role non-stationarity holds in shaping the structure of ABL turbulence and the momentum transport mechanisms it governs.

1. Introduction

Coherent turbulence structures, also known as organized structures, control the exchange of energy, mass and momentum between the Earth's surface and the atmosphere, as well as within engineering systems. In wall-bounded flows, these structures have been shown to carry a substantial fraction of the mean shear stress (Lohou et al. 2000; Katul et al. 2006), kinetic energy (Carper & Porté-Agel 2004; Huang, Cassiani & Albertson 2009; Dong et al. 2020) and scalar fluxes (Li & Bou-Zeid 2011; Wang et al. 2014; Li & Bou-Zeid 2019). It hence comes as no surprise that substantial efforts have been devoted to their characterization across many fields. These structures are of practical relevance in applications relating to biosphere–atmosphere interaction (Raupach, Coppin & Legg 1986; Pan, Chamecki & Isard 2014), air quality control (Michioka, Takimoto & Sato 2014), urban climate (Christen, van Gorsel & Vogt 2007), oceanography (Yang & Shen 2009) and energy harvesting (Ali et al. 2017), to name but a few.

Previous studies on coherent structures in atmospheric boundary-layer (ABL) flows have mainly focused on the roughness sublayer (RSL) and the inertial sublayer (ISL) – the lower portions of the ABL. These layers host physical flow phenomena regulating land–atmosphere exchanges at scales relevant to weather models and human activities (Stull 1988; Oke et al. 2017). The RSL, which extends from the surface up to 2 to 5 times the average height of the roughness elements, is characterized by flow heterogeneity due to the presence of these elements (Fernando 2010). In the RSL, the geometry of turbulence structures is mainly determined by the underlying surface morphology. Through field measurements and wind tunnel data of ABL flow over vegetation canopies, Raupach, Finnigan & Brunet (1996) demonstrated that coherent structures near the top of a vegetation canopy are connected to inflection-point instabilities, akin to those found in mixing layers. Expanding on the framework of this mixing-layer analogy, Finnigan, Shaw & Patton (2009) employed conditional averaging techniques to show that the prevalent eddy structure in the RSL is a head-down hairpin vortex followed by a head-up one. This pattern is characterized by a local pressure peak and a strong scalar front located between the hairpin pair. More recently, Bailey & Stoll (2016) challenged this observation by proposing an alternative two-dimensional roller structure with streamwise spacing that scales with the characteristic length suggested by Raupach et al. (1996).

Extending the mixing-layer analogy to the urban RSL has proven challenging. In a numerical simulation study, Coceal et al. (2007) discovered the absence of Kelvin–Helmholtz waves, which are a characteristic of the mixing-layer analogy, near the top of the urban canopy. This finding, corroborated by observations from Huq et al. (2007), suggests that the mixing-layer analogy is not applicable to urban canopy flows. Instead, the RSL of urban canopy flows is influenced by two length scales; the first is dictated by the size of individual roughness elements such as buildings and trees, and the second by the imprint of large-scale motions above the RSL. The coexistence of these two length scales can be observed through two-point correlation maps (Castro, Cheng & Reynolds 2006; Reynolds & Castro 2008) and velocity spectra (Basley, Perret & Mathis 2019). However, when the urban canopy has a significant aspect ratio between the building height $h$ and width $w$, such as $h/w > 4$, the momentum transport in the RSL is dominated by mixing-layer-type eddies, as shown by Zhang et al. (2022).

The ISL, located above the RSL, is the geophysical equivalent of the celebrated law-of-the-wall region in high Reynolds number turbulent boundary-layer (TBL) flows. In the absence of thermal stratification effects, mean flow in the ISL displays a logarithmic profile, and the momentum flux remains approximately constant with height (Stull 1988; Marusic et al. 2013; Klewicki et al. 2014). Surface morphology has been shown to impact ISL turbulence under certain flow conditions, and this remains a topic of active research. Volino, Schultz & Flack (2007) highlighted the similarity of coherent structures in the log region of TBL flow over smooth and three-dimensional rough surfaces via a comparison of velocity spectra and two-point correlations of the fluctuating velocity and swirl. Findings therein support Townsend's similarity hypothesis (Townsend 1976), which states that turbulence dynamics beyond the RSL does not depend on surface morphological features, except via their role in setting the length and velocity scales for the outer flow region. The said structural similarity between TBL flows over different surfaces was later confirmed by Wu & Christensen (2007) and Coceal et al. (2007), where a highly irregular rough surface and an urban-like roughness were considered, respectively. However, Volino, Schultz & Flack (2011) later reported pronounced signatures of surface roughness on flow structures beyond the RSL in a TBL flow over two-dimensional bars. Similar observations were also made in a TBL flow over a surface characterized by cross-stream heterogeneity (Anderson et al. 2015a), thus questioning the validity of Townsend's similarity hypothesis. To reconcile these contrasting observations, Squire et al. (2017) argued that structural similarity in the ISL is contingent on the surface roughness features not producing flow patterns significantly larger than their own size. If the surface-induced flow patterns are larger than their own size, then they may control flow coherence in the ISL. For example, cross-stream heterogeneous rough surfaces can induce secondary circulations as large as the boundary-layer thickness, which profoundly modify momentum transport and flow coherence in the ISL (Barros & Christensen 2014; Anderson et al. 2015a).

Although coherent structures in cases with significant surface-induced flow patterns necessitate case-specific analyses, researchers have extensively worked towards characterizing the topology of turbulence in cases that exhibit ISL structural similarity. These analyses have inspired scaling laws (Meneveau & Marusic 2013; Yang, Marusic & Meneveau 2016; Hu, Dong & Vinuesa 2023) and the construction of statistical models (Perry & Chong 1982) for TBL turbulence. In this context, the hairpin vortex packet paradigm has emerged as the predominant geometrical model (Christensen & Adrian 2001; Tomkins & Adrian 2003; Adrian 2007). The origins of this model can be traced back to the pioneering work of Theodorsen (1952), who hypothesized that inclined hairpin or horseshoe-shaped vortices were the fundamental elements of TBL turbulence. This idea was later supported by flow visualizations from laboratory experiments (Bandyopadhyay 1980; Head & Bandyopadhyay 1981; Smith et al. 1991) and high-fidelity numerical simulations (Moin & Kim 1982, 1985; Kim & Moin 1986). In addition to providing evidence for the existence of hairpin vortices, Head & Bandyopadhyay (1981) also proposed that these vortices occur in groups, with their heads describing an envelope inclined at $15^\circ \unicode{x2013}20^\circ$ with respect to the wall. Adrian, Meinhart & Tomkins (2000) expanded on this idea, and introduced the hairpin vortex packet paradigm, which posits that hairpin vortices are closely aligned in a quasi-streamwise direction, forming hairpin vortex packets with a characteristic inclination angle of $15^\circ \unicode{x2013}20^\circ$. Nested between the legs of these hairpins are low-momentum regions, which extend approximately 2–3 times the boundary-layer thickness in the streamwise direction. These low-momentum regions are typically referred to as large-scale motions (Smits, McKeon & Marusic 2011). Flow visualization studies by Hommema & Adrian (2003) and Hutchins et al. (2012) further revealed that ABL structures in the ISL are also organized in a similar manner.

Of relevance for this work is that previous studies on coherent structures have predominantly focused on (quasi-)stationary flow conditions. However, stationarity is a rare occurrence in both ABL and engineering flow systems (Lozano-durán et al. 2020; Mahrt & Bou-Zeid 2020). As discussed in the recent review paper by Mahrt & Bou-Zeid (2020), there are two major drivers of non-stationarity in the ABL. The first involves temporal variations of surface heat flux, typically associated with evening transitions or the passage of individual clouds (Grimsdell & Angevine 2002). The second kind corresponds to time variations of the horizontal pressure gradient driving the flow, which can be induced by modes associated with propagating submeso-scale motions, mesoscale disturbances and synoptic fronts (Monti et al. 2002; Mahrt 2014; Cava et al. 2017). Previous studies have demonstrated that non-stationarity significantly affects flow statistics in the ABL, and can result in deviations from equilibrium turbulence Hicks et al. (2018) reported that, during morning and late afternoon transitions, the rapid change in surface heat flux disrupts the equilibrium turbulence relations. Additionally, several observational studies by Mahrt (2007, 2008) and Mahrt et al. (2013) demonstrated that time variations in the driving pressure gradient can enhance momentum transport under stable atmospheric stratifications. Non-stationarity is also expected to impact the geometry of turbulence in the ABL, but this problem has not received much attention thus far. This study contributes to addressing this knowledge gap by investigating the impact of non-stationarity of the second kind on the topology of coherent structures in ABL turbulence and how it affects the mechanisms controlling momentum transport. The study focuses on flow over urban-like roughness subjected to a time-varying pressure gradient. To represent flow unsteadiness, a pulsatile pressure gradient with a constant average and a sinusoidal oscillating component is selected as a prototype. In addition to its practical implications in areas such as wave-current boundary layers, internal-wave induced flows and arterial blood flows, this flow regime facilitates the analysis of coherent structures, owing to the periodic nature of flow statistics.

Pulsatile flows share some similarities with oscillatory flows, i.e. flow driven by a time-periodic pressure gradient with zero mean. Interestingly, in the context of oscillatory flows, several studies have been devoted to the characterization of coherent structures. For instance, Costamagna, Vittori & Blondeaux (2003) and Salon, Armenio & Crise (2007) carried out a numerical study on transitional and fully turbulent oscillatory flow over smooth surfaces, and observed that streaky structures form at the end of the acceleration phases, then distort, intertwine and eventually break into small vortices. Carstensen, Sumer & Fredsøe (2010) performed a series of laboratory experiments on transitional oscillatory flow, and identified two other major coherent structures, namely, cross-stream vortex tubes, which are the direct consequences of inflectional-point shear layer instability, and turbulent spots, which result from the destruction of near-wall streaky structures such as those in stationary flows. Carstensen, Sumer & Fredsøe (2012) observed turbulent spots in oscillatory flows over sand-grain roughness, suggesting that the presence of such flow structures is independent of surface types, and it was later highlighted by Mazzuoli & Vittori (2019) that the mechanism responsible for the turbulent spot generation is similar over both smooth and rough surfaces. Although the primary modes of variability in oscillatory flows are relatively well understood, the same cannot be said for pulsatile flows. A notable study by Zhang & Simons (2019) on wave-current boundary layers, a form of pulsatile flow, revealed phase variations in the spacing of streaks during the wave cycle. However, a detailed analysis of this particular behaviour is still lacking.

To investigate the structure of turbulence in current-dominated pulsatile flow over surfaces in fully rough aerodynamic flow regimes, we conduct a wall-modelled large-eddy simulation (LES) of flow over an array of surface-mounted cuboids. This study builds on the findings of a companion study that was recently accepted for publication in the Journal of Fluid Mechanics, focusing on the time evolution of flow statistics in pulsatile flow over a similar surface (Li & Giometto 2023). By contrasting findings against a corresponding stationary flow simulation, this study addresses these specific questions: (i) Does flow unsteadiness alter the topology of coherent structures in a time-averaged sense? (ii) How does the geometry of coherent structures evolve throughout the pulsation cycle? (iii) What is the effect of such modifications on the mechanisms governing momentum transfer in the ABL? Answering these questions will achieve a twofold research objective: first, contributing to a better understanding of coherent patterns in pulsatile flow over complex geometries, and second, shedding light on how these patterns regulate momentum transfer.

This paper is organized as follows. Section 2 outlines the numerical procedure and the simulation set-up. First- and second-order statistics are presented and discussed in § 3.1. Section 3.2 focuses on a quadrant analysis, whereas §§ 3.3 and 3.4 interpret the flow field in terms of two-point correlations and instantaneous flow behaviour. Further insight into the time evolution of the turbulence topology is proposed in § 3.5 via conditional averaging. Concluding remarks are given in § 4.

2. Methodology

2.1. Numerical procedure

Simulations are carried out via an in-house LES algorithm (Albertson & Parlange 1999a,b; Giometto et al. 2016). The LES algorithm solves the spatially filtered momentum and mass conservation equations, namely

(2.1)\begin{gather} \frac{\partial u_i}{\partial t} + u_j \left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right) ={-} \frac{1}{ \rho}\frac{\partial P}{ \partial x_i} - \frac{\partial \tau_{ij}}{\partial x_j} - \frac{1}{\rho} \frac{\partial P_\infty}{ \partial x_1}\delta_{i1}+F_i , \end{gather}

(2.2)\begin{gather}\frac{\partial u_i}{\partial x_i} = 0, \end{gather}

where $(u_1,u_2,u_3 )$ represent the filtered velocities along the streamwise $x_1$, cross-stream $x_2$ and wall-normal $x_3$ directions, respectively. The rotational form of the convective term is used to ensure kinetic energy conservation in the discrete sense in the inviscid limit (Orszag & Pao 1975). Also, $\tau _{ij}$ is the deviatoric part of the kinematic subgrid-scale (SGS) stress tensor, parameterized via the Lagrangian scale-dependent dynamic Smagorinsky model (Bou-Zeid, Meneveau & Parlange 2005). The flow is assumed to be in the fully rough aerodynamic regime, and viscous stresses are not considered. Further, $P=p+\rho \frac {1}{3}\tau _{ii}+\rho \frac {1}{2} u_i u_i$ is a modified pressure, which accounts for the trace of SGS stress and resolved turbulent kinetic energy, and $\rho$ is a constant fluid density. The flow is driven by a spatially uniform, pulsatile pressure gradient in the $x_1$ direction, namely ${\partial P_\infty }/{ \partial x_1} = -\rho f_m[ 1+\alpha _p \sin (\omega t) ]$, where the $f_m$ parameter controls the magnitude of the temporally averaged pressure gradient, $\alpha _p$ controls the forcing amplitude and $\omega$ the forcing frequency; $\delta _{ij}$ in (2.1) denotes the Kronecker delta tensor.

Periodic boundary conditions apply in the wall-parallel directions, and a free-slip boundary condition is imposed at the top of the computational domain. The lower surface consists of an array of uniformly distributed cuboids, which are explicitly resolved via a discrete forcing immersed boundary method (IBM) (Mittal & Iaccarino 2005). The IBM approach makes use of an artificial force $F_i$ to impose the no-slip boundary condition at the solid–fluid interfaces. Additionally, it utilizes an algebraic equilibrium wall-layer model to evaluate surface stresses (Piomelli 2008; Bose & Park 2018). The algorithm has been extensively validated for the simulation of flow in urban environments (see, e.g. Tseng, Meneveau & Parlange 2006; Chester, Meneveau & Parlange 2007; Giometto et al. 2016).

Spatial derivatives in the wall-parallel directions are computed via a pseudo-spectral collocation method based on truncated Fourier expansions (Orszag 1970), whereas a second-order staggered finite differences scheme is employed in the wall-normal direction. Since dealiasing errors are known to be detrimental for pseudo-spectral discretization (Margairaz et al. 2018), nonlinear convective terms are de-aliased exactly via the $3/2$ rule (Canuto et al. 2007). The time integration is performed via a second-order Adams–Bashforth scheme, and the incompressibility condition is enforced via a fraction step method (Kim & Moin 1985).

2.2. Simulation set-up

Two LESs of flow over an array of surface-mounted cubes are carried out. The two simulations only differ by the pressure forcing term: one is characterized by a pressure gradient that is constant in space and time (CP hereafter), and the other by a pressure gradient that is constant in space and pulsatile in time (PP).

The computational domain for both simulations is sketched in figure 1. The size of the box is $[0,L_1]\times [0,L_2]\times [0,H]$ with $L_1=72h$, $L_2=24h$ and $H=8h$, where $h$ denotes the height of the cubes. Cubes are organized in an in-line arrangement with planar and frontal area fractions set to $\lambda _p = \lambda _f = 0.\bar {1}$. The relatively high packing density along with the chosen scale separation $H/h = 8$ support the existence of a well-developed ISL and healthy coherent structures in the considered flow system (Castro 2007; Coceal et al. 2007; Zhang et al. 2022). In terms of horizontal extent, $L_1/H$ and $L_2/H$ are larger than those from previous works focusing on coherent structures above aerodynamically rough surfaces (Coceal et al. 2007; Xie, Coceal & Castro 2008; Leonardi & Castro 2010; Anderson, Li & Bou-Zeid 2015b) and are sufficient to accommodate large-scale motions (Balakumar & Adrian 2007). An aerodynamic roughness length $z_0=10^{-4} h$ is prescribed at the cube surfaces and the ground via the algebraic wall-layer model, resulting in negligible SGS drag contributions to the total surface drag (Yang & Meneveau 2016). The computational domain is discretized using a uniform Cartesian grid of $N_1 \times N_2 \times N_3 = 576 \times 192 \times 128$, so each cube is resolved via $8 \times 8 \times 16$ grid points. Such a grid resolution yields flow statistics that are poorly sensitive to grid resolution in both statistically stationary and pulsatile flows at the considered oscillation frequency (Tseng et al. 2006; Li & Giometto 2023).

Figure 1. Side and planar views of the computational domain (a,b, respectively). The red dashed line denotes the repeating unit.

For the PP case, the forcing frequency is set to $\omega T_h={\rm \pi} /8$, where $T_h = h/u_\tau$ is the averaged turnover time of characteristic eddies of the urban canopy layer (UCL) and ${u}_{\tau }=\sqrt {f_m H}$ the friction velocity. This frequency selection is based on both practical and theoretical considerations. Realistic ranges for the friction velocity and UCL height are $0.1 \leq {u}_{\tau } \leq 0.5 \ {\rm m}\ {\rm s}^{-1}$ and $3 \leq h \leq 30 \ {\rm {m}}$ (Stull 1988). Using such values, the chosen frequency corresponds to a time period $24 \leq T \leq 4800 \ {\rm {s}}$, where $T=2{\rm \pi} /\omega = 16 T_h$. This range of time scales pertains to sub-mesoscale motions (Mahrt 2009; Hoover et al. 2015), which, as outlined in § 1, are a major driver of atmospheric pressure gradient variability. From a theoretical perspective, this frequency is expected to yield substantial modifications of coherent structures within the ISL. The chosen frequency results in a Stokes layer thickness $\delta _s=5h$, encompassing both the RSL and the ISL. Within the Stokes layer, turbulence generation and momentum transport undergo significant modifications during the pulsation cycle, as demonstrated in Li & Giometto (2023). Moreover, the oscillation period $T$ is comparable to the average lifespan of eddies in the ISL of the considered flow system, as elaborated below. Coceal et al. (2007) showed that, in flow over rough surfaces, the characteristic length scale of ISL eddies ($\ell$) is bounded below by $h$, thus yielding $\min {(\ell )} \sim h$. Based on Townsend's attached-eddy hypothesis, $\ell \sim x_3$, which results in $\max {(\ell )} \sim H$. The time scale associated with ISL eddies is $T_\ell = \ell /u_\tau$, so that $\min {(T_\ell )} \sim h/u_\tau = T_h$ and $\max {(T_\ell )} \sim H/u_\tau = T_H$. The modest scale separation characterizing our set-up ($H = 8h$) yields a modest separation of time scales in the ISL, and considering $T \approx T_H$, one can conclude that the time scale of ISL eddies is comparable to $T$. With $T_\ell \approx T$, flow pulsation will considerably modify the structure of ISL turbulence and drive the flow out of equilibrium conditions. This is because changes in the imposed pressure gradient occur at a rate that enables turbulent eddies to respond. This behaviour can be contrasted with two limiting cases: with $T_\ell \gg T$, turbulence is unable to respond to the rapid changes in the environment and is advected in a ‘frozen’ state, i.e. it does not undergo topological changes. With $T_\ell \ll T$, ISL eddies do not ‘live’ long enough to sense changes in the environment, and maintain a quasi-steady state throughout the pulsatile cycle. In terms of forcing amplitude, such a quantity is set to $\alpha _p=12$ for the PP case; this amplitude is large enough to induce significant changes in the coherent structures with the varying pressure gradient while avoiding mean flow reversals.

Both simulations are initialized with velocity fields from a stationary flow case and integrated over $400 T_{H}$, corresponding to $200$ pulsatile cycles for the PP case. Here, $T_{H}= H/{u}_{\tau }$ refers to the turnover time of the largest eddies in the domain. The time step ($\delta t$) is set to ensure $\max {(CFL)}=u_{max} \delta t/\delta \approx 0.05$, where CFL denotes the Courant–Friedrichs–Lewy stability condition, $u_{max}$ is the maximum velocity magnitude at any point in the domain during the simulation and $\delta$ is the smallest grid stencil in the three coordinate directions. The initial $20 T_{H}$ are discarded for both the CP and PP cases (transient period for the PP case), which correspond to approximately 10 oscillation periods, after which instantaneous snapshots of velocities and pressure are saved to disk every $0.025T_{H}$ ($1/80$ of the pulsatile cycle).

2.3. Notation and terminology

For the PP case, $\overline {({\cdot })}$ denotes an ensemble averaging operation, performed over the phase dimension and over repeating surface units (see figure 1), i.e.

(2.3)\begin{align} \bar{\theta} (x_1,x_2,x_3,t)&= \frac{1}{N_p n_1 n_2}\sum^{N_p}_{n=1} \sum^{n_1}_{i=1}\sum^{n_2}_{j=1} \theta(x_1+il_1,x_2+jl_2,x_3,t+nT) ,\nonumber\\ &\qquad 0\leq x_1 \leq l_1,\quad 0\leq x_2 \leq l_2,\quad 0\leq t \leq T , \end{align}

where $\theta$ is a given scalar field and $n_1$ and $n_2$ are the number of repeating units in the streamwise and cross-stream directions, respectively. Using the usual Reynolds decomposition, one can write

(2.4)\begin{equation} \theta (x_1,x_2,x_3,t)= \bar{\theta} (x_1,x_2,x_3,t)+\theta^\prime(x_1,x_2,x_3,t) , \end{equation}

where $({\cdot })^\prime$ denotes a fluctuation from the ensemble average. For the CP case, $\overline {({\cdot })}$ denotes a quantity averaged over time and repeating units. An ensemble-averaged quantity can be further decomposed into an intrinsic spatial average and a deviation from the intrinsic average (Schmid et al. 2019), i.e.

(2.5)\begin{equation} \bar{\theta} (x_1,x_2,x_3,t)= \langle \bar{\theta} \rangle (x_3,t)+ \bar{\theta} ^{\prime \prime}(x_1,x_2,x_3,t) . \end{equation}

Note that, for each $x_3$, the intrinsic averaging operation is taken over a thin horizontal ‘slab’ $V_f$ of fluid, characterized by a thickness $\delta _3$ in the wall-normal ($x_3$) direction, namely

(2.6)\begin{equation} \langle \bar{\theta} \rangle (x_3,t)= \frac{1}{V_f}\int_{x_3-\delta_3/2}^{x_3+\delta_3/2}\int_0^{l_2} \int_0^{l_1}\bar{\theta} (x_1,x_2,x_3^{\prime},t) \,{\rm d}\kern0.7pt x_1 \,{\rm d}\kern0.7pt x_2 \,{\rm d}\kern0.7pt x_3^{\prime} . \end{equation}

Further, any phase-averaged quantity from the PP case consists of a long-time-averaged component and an oscillatory component with a zero mean, which will be hereafter denoted via the subscripts $l$ and $o$, respectively, i.e.

(2.7)\begin{equation} \bar{\theta} (x_1,x_2,x_3,t)= \bar{\theta}_{l} (x_1,x_2,x_3)+ \bar{\theta}_{o} (x_1,x_2,x_3,t), \end{equation}

and

(2.8)\begin{equation} \langle \bar{\theta} \rangle (x_3,t)= \langle \bar{\theta} \rangle_{l} (x_3) + \langle \bar{\theta} \rangle_{o} (x_3,t) . \end{equation}

As for the CP case, the long-time and ensemble averages are used interchangeably due to the lack of an oscillatory component. In the following, the long-time-averaged quantities from the PP case are contrasted against their counterparts from the CP case to highlight the impact of flow unsteadiness on flow characteristics in a long-time-average sense. Oscillatory and phase-averaged quantities are analysed to shed light on the phase-dependent features of the PP case.

3. Results

3.1. Overview of flow statistics

Li & Giometto (2023) have proposed a detailed analysis of pulsatile flow over an array of surface-mounted cuboids, discussing the impact of varying forcing amplitude and frequency on selected flow statistics. Here, we repropose and expand upon some of the findings for the chosen oscillation frequency and amplitude that are relevant to this work.

Figure 2(a) presents the wall-normal distributions of the long-time-averaged resolved Reynolds shear stress $\langle \overline {u^\prime _1 u^\prime _3} \rangle _l$ and dispersive shear stress $\langle \bar {u}^{\prime \prime }_1 \bar {u}^{\prime \prime }_3 \rangle _l$. Note that SGS components contribute less than $1\,\%$ to the total Reynolds stresses and are hence not discussed. From the figure, it is apparent that flow unsteadiness does not noticeably affect the $\langle \overline {u^\prime _1 u^\prime _3} \rangle _l$ profile, with local variations from the statistically stationary scenario being within a $3\,\%$ margin. On the contrary, flow pulsation within the UCL leads to pronounced increases in $\langle \bar {u}^{\prime \prime }_1 \bar {u}^{\prime \prime }_3 \rangle _l$, with local surges reaching up to a fivefold increase. However, despite this increase, the dispersive flux remains a modest contributor to the total momentum flux in the UCL. Figure 2(b) displays the long-time-averaged resolved turbulent kinetic energy $k_l = \langle \overline {u^\prime _i u^\prime _i} \rangle _l/2$ and wake kinetic energy $k_{w,l} = \langle \bar {u}^{\prime \prime }_i \bar {u}^{\prime \prime }_i \rangle _l/2$. Both $k_l$ and $k_{w,l}$ from the PP case feature modest (${<}5\,\%$) local departures from their CP counterparts, highlighting a weak dependence of both long-time-averaged turbulent and wake kinetic energy on flow unsteadiness. Also, the RSL thicknesses $(x_3^R)$ for the CP and PP cases are depicted in figure 2. Following the approach by Pokrajac et al. (2007), $x_3^R$ is estimated by thresholding the spatial standard deviation of the long-time-averaged streamwise velocity normalized by its intrinsic average, namely

(3.1)\begin{equation} \sigma=\frac{\sqrt{\langle(\bar{u}_{1,l}-\langle\bar{u}_1\rangle_l)^2\rangle}}{\langle\bar{u}_1\rangle_l} , \end{equation}

where the threshold is taken as 1 %. An alternative method to evaluate $x_3^R$ involves using phase-averaged statistics instead of long-time-averaged ones in (3.1). Although not shown, such a method yields similar predictions (with a discrepancy of less than $5\,\%$). Both $\langle \bar {u}^{\prime \prime }_1 \bar {u}^{\prime \prime }_3 \rangle _l$ and $k_{w,l}$ reduce to less than $1\,\%$ of their peak value above $x_3^R$. From figure 2, one can readily observe that flow unsteadiness yields a modest increase in the extent of the RSL, with an estimated $x_3^R$ not exceeding $1.5h$ in both cases. Hereafter, we will hence assume $x_3^R=1.5h$. As discussed in § 1, RSL and ISL feature distinct coherent structures. Specifically, the structures in the RSL are expected to show strong imprints of roughness elements, whereas those in the ISL should, in principle, be independent of surface morphology (Coceal et al. 2007).

Figure 2. (a) Long-time-averaged shear stresses from the PP (black) and CP (red) cases. Resolved Reynolds shear stress $\langle \overline {u^\prime _1 u^\prime _3} \rangle _l$, solid lines; dispersive shear stress $\langle \bar {u}^{\prime \prime }_1 \bar {u}^{\prime \prime }_3 \rangle _l$. (b) Long-time-averaged turbulent and wake kinetic energy from the PP (black) and CP (red) cases. Resolved turbulent kinetic energy $k_l = \langle \overline {u^\prime _i u^\prime _i} \rangle _l/2$, solid lines; wake kinetic energy $k_{w,l}=\langle \bar {u}^{\prime \prime }_i \bar {u}^{\prime \prime }_i \rangle _l/2$, dashed lines. Dashed-dotted horizontal lines denote the upper bound of the RSL $(x_3^R)$.

The response of selected first- and second-order flow statistics to flow unsteadiness is depicted in figure 3. In figure 3(a), an oscillating wave is evident in the oscillatory shear rate $\partial \langle \bar {u}_1 \rangle _o /\partial x_3$. This wave, generated at the canopy top due to flow unsteadiness, exhibits a phase lag of ${\rm \pi} /2$ relative to the pulsatile pressure forcing. Such a wave propagates in the positive vertical direction while being attenuated and diffused by turbulent mixing. It is noteworthy that the propagation speed of the oscillating shear rate is to a good degree constant, as suggested by the constant tilting angle along the $x_3$ direction of the ${\partial \langle \bar {u}_1 \rangle _o}/{\partial x_3}$ contours. As apparent from figure 3(b,c), the space–time diagrams of the oscillatory resolved Reynolds shear stress $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and oscillatory resolved turbulent kinetic energy $k_o=\langle \overline {u^\prime _i u^\prime _i}\rangle _o/2$ are also characterized by decaying waves travelling away from the RSL at constant rates. The speeds of these waves are similar to that of the corresponding oscillating shear rate, which can be again inferred by the identical tilting angles in the contours. There is clearly a causal relation for this behaviour: above the UCL, the major contributors of shear production terms in the budget equations of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and $k_o$ are

(3.2) \begin{equation} \langle\overline{\mathcal{P}}\rangle_{13,o}= -2\langle \overline{u_3^\prime u_3^\prime}\rangle_l \frac{\partial \langle \overline{u}_1 \rangle_o}{\partial x_3}-2\langle \overline{u_3^\prime u_3^\prime}\rangle_o \frac{\partial \langle \overline{u}_1 \rangle_l }{\partial x_3} ,\end{equation}

and

(3.3) \begin{equation} \langle\overline{\mathcal{P}}\rangle_{k,o} = -\langle \overline{u_1^\prime u_3^\prime}\rangle_l \frac{\partial \langle \overline{u}_1 \rangle_o}{\partial x_3}-\langle \overline{u_1^\prime u_3^\prime}\rangle_o \frac{\partial \langle\overline{u}_1\rangle_l}{\partial x_3}, \end{equation}

respectively. As the oscillating shear rate travels upwards away from the UCL, it interacts with the local turbulence by modulating $\langle \bar {\mathcal {P}}\rangle _{13,o}$ and $\langle \bar {\mathcal {P}}\rangle _{k,o}$, ultimately yielding the observed oscillations in resolved Reynolds stresses. On the other hand, no pulsatile-forcing-related terms appear in the budget equations of resolved Reynolds stresses. This indicates that the oscillating shear rate induced by the pulsatile forcing modifies the turbulence production above the UCL, rather than the pressure forcing itself. A similar point about pulsatile flows was made in Scotti & Piomelli (2001), where it was stated that ‘[$\ldots$]in the former [pulsatile flow] it is the shear generated at the wall that affects the flow’. It is worth noting that such a study was, however, based on pulsatile flow over smooth surfaces and at a relatively low Reynolds number.

Figure 3. Space–time diagrams of (a) oscillatory shear rate ${\partial \langle \bar {u}_1 \rangle _o}/{\partial x_3}$, (b) oscillatory resolved Reynolds shear stress $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and (c) oscillatory resolved turbulent kinetic energy $k_o=\langle \overline {u^\prime _i u^\prime _i}\rangle _o/2$ from the PP case. Results are normalized by $u_{\tau }$ and $h$. Horizontal dashed lines highlight the top of the UCL.

In addition, a visual comparison of the contours of ${\partial \langle \overline{u}_1\rangle_o }/{\partial x_3}$ and $-\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ highlights the presence of a phase lag between such quantities throughout the flow field. Further examination of this phase lag can be found in Li & Giometto (2023). During the pulsatile cycle, the turbulence is hence not in equilibrium with the mean flow. This is the case despite the fact that neither the pulsatile forcing nor the induced oscillating shear wave significantly alters the long-time-averaged turbulence intensity, as evidenced in figure 2. To gain further insight into this behaviour, the next section examines the structure of turbulence under this non-equilibrium condition.

3.2. Quadrant analysis

The discussions will first focus on the impact of flow pulsation on the $u_1^\prime u_3^\prime$ quadrants, with a focus on the ISL. This statistical analysis enables the quantification of contributions from different coherent motions to turbulent momentum transport. The quadrant analysis technique was first introduced by Wallace, Eckelmann & Brodkey (1972), and has thereafter been routinely employed to characterize the structure of turbulence across a range of flow systems (Wallace 2016). The approach maps velocity fluctuations to one of four types of coherent motions (quadrants) in the $u_1^\prime -u_3^\prime$ phase space, namely

(3.4)\begin{equation} \begin{cases} Q1: & u_1^\prime>0, u_3^\prime>0 , \\ Q2: & u_1^\prime<0, u_3^\prime>0 , \\ Q3: & u_1^\prime<0, u_3^\prime<0 , \\ Q4: & u_1^\prime>0, u_3^\prime<0 . \end{cases} \end{equation}

The quantities $Q$2 and $Q$4 are typically referred to as ejections and sweeps, respectively. They are the main contributors to the Reynolds shear stress, and constitute the majority of the events in boundary layer flows. Ejections are associated with the lift-up of low-momentum fluid by vortex induction between the legs of hairpin structures, whereas sweeps correspond to the down-draft of the high-momentum fluid (Adrian et al. 2000); $Q$1 and $Q$3 denote outward and inward interactions, and play less important roles in transporting momentum when compared with $Q$2 and $Q$4. Coceal et al. (2007) and Finnigan (2000) showed that the RSL of stationary flows is dominated by ejections in terms of the number of events, but the overall Reynolds stress contribution from sweep events exceeds that of ejections. This trend reverses in the ISL. This behaviour is indeed apparent from figure 4, where ejection and sweep profiles are shown for the CP case (red lines).

Figure 4. (a) Relative contribution to $\overline {u_1^\prime u_3^\prime }$ by events in each quadrant summed over the wall-parallel planes and the whole sampling time period and (b) relative number of events in each quadrant from the PP case (black) and CP (red) as a function of $x_3$. Cross: outward interaction; triangles: ejection; diamonds: inward interaction; circles: sweep.

We first examine the overall frequency of events in each quadrant and the contribution of each quadrant to the resolved Reynolds shear stress. For the considered cases, the contribution to $\overline {u_1^\prime u_3^\prime }$ and the number of the events of each quadrant are summed over different wall-parallel planes and over the whole sampling time period (i.e. these are long-time-averaged quantities). Results from this operation are also shown in figure 4. What emerges from this analysis is that flow pulsation does not significantly alter the relative contribution and frequency of each quadrant. Some discrepancies between CP and PP profiles can be observed immediately above the UCL, but do not sum to more than 4 % at any given height.

A more interesting picture of the flow field emerges if we consider the phase-dependent behaviour of ejections and sweeps. Hereafter, the ratio between the numbers of ejections and sweeps is denoted by $\gamma _\#$, and the ratio of their contribution to $\overline {u_1^\prime u_3^\prime }$ by $\gamma _c$. As outlined in the previous section, turbulent fluctuations are defined as deviations from the local ensemble average. Consequently, both the frequency of occurrences and the contribution to $\overline {u_1^\prime u_3^\prime }$ from each quadrant are influenced by two main factors: the relative position to the cube within the repeating unit and the phase in the PP case. This dual dependency extends to $\gamma _\#$ and $\gamma _c$ as well. Conversely, in the CP case, $\gamma _\#$ and $\gamma _c$ are only functions of the spatial location relative to the cube. Figure 5(a,c) presents $\gamma _\#$ up to $x_3/h=2$ at a selected streamwise/wall-normal plane for the PP and CP cases, respectively. The chosen plane cuts through the centre of a cube in the repeating unit, as shown in 5(b). In the cavity, the ejection-sweep pattern from the PP case is found to be qualitatively similar to its CP counterpart throughout the pulsatile cycle (compare panels (a,c) in figure 5). Specifically, a preponderance of sweeps characterizes a narrow region on the leeward side of the cube (the streamwise extent of this region is ${\lessapprox }0.3h$), whereas ejections dominate in the remainder of the cavity. As also apparent from figure 5(a), the streamwise extent of the sweep-dominated region features a modest increase (decrease) during the acceleration (deceleration) time period. During the acceleration phase, the shown above canopy region $(h< x_3<2h)$ transitions from an ejection-dominated flow regime to a sweep-dominated one, and vice versa as the flow decelerates. This transition initiates just above the cavity, characterized by a higher occurrence of sweeps during the acceleration phase and a predominance of ejections in the deceleration period. This continues until both phenomena are distributed throughout the RSL. While not discussed in this work, it is worth noting that the trend observed for $\gamma _c$ is precisely the inverse.

Figure 5. (a) Ratio between the numbers of ejections to sweeps ($\gamma _\#$) from the PP case on a streamwise/wall-normal plane. (b) Location of the selected streamwise/wall-normal plane (red dashed line) within a repeating unit. (c) Value of $\gamma _\#$ from the CP case on the same plane. Black dashed lines denote $x_3/h=1.5$, which is the upper limit of the RSL.

Shifting the attention to the ejection-sweep pattern in the ISL, which is indeed the main focus of this study, figure 6 shows the intrinsic average of $\gamma _c$ and $\gamma _\#$ in the $x_3/h = \{2, 3, 4 \}$ planes. These quantities are hereafter denoted as $\langle \gamma _c \rangle$ and $\langle \gamma _\# \rangle$, respectively. The use of $\langle \gamma _c \rangle$ and $\langle \gamma _\# \rangle$ instead of $\gamma _c$ and $\gamma _\#$ to characterize the ejection-sweep pattern in the ISL can be justified by the fact that the spatial variations in $\gamma _\#$ and $\gamma _c$ in the wall-parallel directions vanish rapidly above the RSL, as apparent from figure 5. This is in line with the observations of Kanda, Moriwaki & Kasamatsu (2004) and Castro et al. (2006) that the spatial variations in $\gamma _\#$ and $\gamma _c$ are concentrated in the RSL for stationary flow over an urban canopy. Further, as shown in figure 6, the ejection-sweep pattern varies substantially during the pulsatile cycle. For instance, at a relative height of $x_3/h=2$, even though the contribution from ejections to $\overline { u_1^\prime u_3^\prime }$ dominates in a long-time-average sense ($\langle \gamma _c \rangle _l> 1$), sweep contributions prevail for $\omega t \in [0,{\rm \pi} /2]$. Interestingly, at a given wall-normal location, this ejection-sweep pattern appears to be directly controlled by the intrinsic- and phase-averaged shear rate ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$. This is particularly evident when $\langle \gamma _c \rangle$ and $\langle \gamma _\# \rangle$ are plotted against ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ (refer to figure 7). As ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ increases at a given $x_3$, the corresponding $\langle \gamma _c \rangle$ increases whereas $\langle \gamma _\# \rangle$ decreases, highlighting the presence of fewer but stronger ejections events. Maxima and minima of $\langle \gamma _c \rangle$ and $\langle \gamma _\# \rangle$ approximately coincide with the maxima of ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$. This observation is consistent across the considered planes. As discussed in the next sections, such behaviour can be attributed to time variations in the geometry of ISL structures.

Figure 6. (a–c) Intrinsic-averaged ratio of contributions to $\overline { u_1^\prime u_3^\prime }$ from ejections and sweeps ($\langle \gamma _c \rangle$); (d–f) intrinsic-averaged ratio of ejections to sweeps ($\langle \gamma _\# \rangle$); (g–i) intrinsic- and phase-averaged shear rate ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ from the PP case at three wall-normal locations within the ISL (a,d,g) $x_3/h=2$, (b,e,h) $x_3/h=3$ and (c,f,i) $x_3/h=4$ as a function of phase. Black dashed lines denote long-time-averaged values, whereas solid red lines represent corresponding quantities from the CP case.

Figure 7. (a) Values of $\langle \gamma _c \rangle$ and (b) $\langle \gamma _\# \rangle$ vs ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ at $x_3/h=2$ (blue), $x_3/h=3$ (green) and $x_3/h=4$ (magenta).

3.3. Spatial and temporal flow coherence

To gain a better understanding of the extent and organization of coherent structures in the ISL, this section analyses two-point velocity autocorrelation maps. These flow statistics provide information on the correlation of the flow field in space, making it an effective tool for describing spatial flow coherence (Dennis & Nickels 2011; Guala et al. 2012). For the PP case, the phase-dependent two-point correlation coefficient tensor $\bar {R}_{ij}$ can be defined as

(3.5)\begin{equation} \bar{R}_{ij}(\varDelta_1,\varDelta_2,x_3,x_3^*,t)=\frac{\langle \overline{ u_i^\prime(x_1,x_2,x_3^*,t) u_j^\prime(x_1+\varDelta_1,x_2+\varDelta_2,x_3,t)} \rangle}{\sqrt{\langle \overline{u_i^\prime u_i^\prime}\rangle (x_3^*,t)\langle \overline{ u_j^\prime u_j^\prime}\rangle (x_3,t)}} , \end{equation}

where $\varDelta _i$ is the separation in the wall-parallel directions, $x_3^*$ represents a reference wall-normal location and $t$ denotes the phase. In the CP case, the flow is statistically stationary, and therefore $\bar {R}_{ij}$ is not a function of $t$, i.e. $\bar {R}_{ij} = \bar {R}_{ij,l}$.

Figure 8 compares $\bar {R}_{11,l}$ for the PP and CP cases over the $x_3^*/h = \{1.5, 2, 3, 4 \}$ planes. In both cases, $\bar {R}_{11,l}$ features an alternating sign in the cross-stream direction, signalling the presence of low- and high-momentum streaks flanking each other in the cross-stream direction. The cross-stream extent of long-time-averaged streaks can be identified as the first zero crossing of the $\bar {R}_{11,l}$ contour in the $\varDelta _2$ direction. Based on this definition, figure 8 shows that flow unsteadiness has a modest impact on such a quantity. This finding agrees with observations from Zhang & Simons (2019) for pulsatile flow over smooth surfaces. Further, although not shown, the streamwise and cross-stream extent of streaks increases linearly in $x_3$, suggesting that Townsend's attached-eddy hypothesis is valid in a long-time-average sense (Marusic & Monty 2019).

Figure 8. Long-time-averaged two-point correlation coefficient tensor $\bar {R}_{11,l}$ at (a) $x_3^*/h=1.5$, (b) $x_3^*/h=2$, (c) $x_3^*/h=3$ and (d) $x_3^*/h=4$. Black lines correspond to the PP case, and red lines to the CP one. Here, $\bar {R}_{11,l}=0.6$ and $\bar {R}_{11,l}=0.3$ are denoted by solid lines, and dashed lines represent $\bar {R}_{11,l}=0$.

Turning the attention to the phase-averaged flow field, figure 9 shows the time variation of the cross-stream streaks extent, which is identified as the first zero crossing of the $\bar {R}_{11}=0$ field in the cross-stream direction. The linear $x_3$-scaling of the streak width breaks down in a phase-averaged sense. Such a quantity indeed varies substantially during the pulsatile cycle, diminishing in magnitude as ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ increases throughout the boundary layer. Interestingly, when ${\partial \langle \bar {u}_1 \rangle }/{\partial x_3}$ reaches its maximum at $\omega t \approx {\rm \pi}$ and $x_3/h \approx 1.5$, the cross-stream extent of streaks approaches zero, suggesting that streaks may not be a persistent feature of pulsatile boundary-layer flows.

Figure 9. Time evolution of (a) the cross-stream streak width normalized by $h$ and (b) $\partial \langle \bar {u}_1 \rangle /\partial x_3$. The cross-stream width is identified as the first zero crossing of the $\bar {R}_{11}=0$ field.

To further quantify topological changes induced by flow pulsation, we hereafter examine variations in the streamwise and wall-normal extent of coherent structures. Such quantities will be identified via the $\bar {R}_{11}=0.3$ contour, in line with the approach used by Krogstad & Antonia (1994). Note that the choice of the $\bar {R}_{11}$ threshold for such a task is somewhat subjective, and several different values have been used in previous studies to achieve this same objective, including $\bar {R}_{11}=0.4$ (Takimoto et al. 2013) and $\bar {R}_{11}=0.5$ (Volino et al. 2007; Guala et al. 2012). In this study, the exact threshold is inconsequential as it does not impact the conclusions. Figure 10 presents $\bar {R}_{11,l}$ contours in the streamwise/wall-normal plane for $x_3^*/h = \{ 1.5,2,3,4 \}$. The jagged lines at $x_3/h \approx 1$ (the top of the UCL) bear the signature of roughness elements. The dashed lines passing through $x_3^*$ identify the locus of the maxima in $\bar {R}_{11,l}$ at each streamwise location. The inclination angle of such lines can be used as a surrogate for the long-time-averaged tilting angle of the coherent structure (Chauhan et al. 2013; Salesky & Anderson 2020). It is clearly observed that, at each reference wall-normal location, the tilting angle of long-time-averaged structures is similar between the PP case and CP case. The tilting angle in both cases decreases monotonically and slowly from $15^\circ$ at $x_3^*/h=1.5$ to $10^\circ$ at $x_3^*/h=4$ – a behaviour that is in excellent agreement with results from Coceal et al. (2007), even though a different urban canopy layout was used therein. Further, the identified tilting angle is also similar to the one inferred from real-world ABL observations in Hutchins et al. (2012) and Chauhan et al. (2013). On the other hand, long-time-averaged coherent structures in the PP case are relatively smaller than in the CP case in both the streamwise and wall-normal coordinate directions. Discrepancies become more apparent with increasing $x_3^*$. Specifically, the difference in the streamwise extent of the long-time-averaged structure from the two cases increases from $2\,\%$ at $x_3^*/h=1.5$ to $15\,\%$ at $x_3^*/h=4$. Corresponding variations in the wall-normal extent are $2\,\%$ and $4\,\%$.

Figure 10. Value of $\bar {R}_{11,l}$ in the streamwise/wall-normal plane of the PP (black) and CP (red) cases. Results correspond to four reference wall-normal locations: (a) $x_{3}^*/h=1.5$, (b) $x_{3}^*/h=2$, (c) $x_{3}^*/h=3$ and (d) $x_{3}^*/ h=4$. Contour levels (solid lines) range from $0.2$ to $0.5$ with increments of $0.1$. Dashed lines denote the locus of the maximum correlation at each streamwise location. The slopes of the dashed lines represent the tilting angles of the structures.

More insight into the mechanisms underpinning the observed behaviour can be gained by examining the time evolution of such structures for the PP case in figure 11. When taken together with figure 9(b), it becomes clear that both the streamwise and the wall-normal extents of the coherent structures tend to reduce with increasing local $\partial \langle \bar {u}_1 \rangle /\partial x_3$. Compared with the streamwise extent, the wall-normal extent of the coherent structure is more sensitive to changes in $\partial \langle \bar {u}_1 \rangle /\partial x_3$. For example, at $x_3^*/h=4$, we observe an overall $15\,\%$ variation in the wall-normal extent of the coherent structure during a pulsation cycle, whereas the corresponding variation in streamwise extent is $8\,\%$. Further, the flow field at the considered heights appears to be more correlated with the flow in the UCL for small $\partial \langle \bar {u}_1 \rangle /\partial x_3$, thus highlighting a stronger coupling between flow regions in the wall-normal direction. Interestingly, the tilting angle of the coherent structure remains constant during the pulsatile cycle, as shown in figure 12.

Figure 11. Time evolution of $\bar {R}_{11}=0.3$ in the streamwise/wall-normal plane. Line colours denote the contours corresponding to different $x_3^*$ planes: $x_3^*/h=1.5$ (black), $x_3^*/h=2$ (blue), $x_3^*/h=3$ (green) and $x_3^*/h=4$ (magenta). Dots highlight the location of the reference plane.

Figure 12. The locus of the maximum $\bar {R}_{11}$ at four phases: $\omega t=0$ (solid lines), $\omega t={\rm \pi} /2$ (dashed lines), $\omega t={\rm \pi}$ (dashed dotted lines) and $\omega t=3{\rm \pi} /2$ (dotted lines). Line colours denote different reference elevations: $x_3^*/h=1.5$ (black), $x_3^*/h=2$ (blue), $x_3^*/h=3$ (green) and $x_3^*/h=4$ (magenta).

Next, we will show that the hairpin vortex packet paradigm (Adrian 2007) can be used to provide an interpretation for these findings. Note that alternative paradigms, such as that proposed by Del Alamo et al. (2006), may offer different interpretations of the results, but are not discussed in this work. The validity of such a paradigm is supported by a vast body of evidence from laboratory experiments of canonical TBL (Adrian et al. 2000; Christensen & Adrian 2001; Dennis & Nickels 2011) to ABL field measurements (Hommema & Adrian 2003; Morris et al. 2007) and numerical simulations (Lee, Sung & Krogstad 2011; Eitel-Amor et al. 2015). This formulation assumes that the dominant ISL structures are hairpin vortex packets, consisting of a sequence of hairpin vortices organized in a quasi-streamwise direction with a characteristic inclination angle relative to the wall. These structures encapsulate the low-momentum regions, also known as ‘streaks’. The structural information obtained from the two-point correlation has been considered to reflect the averaged morphology of the hairpin vortex packets (Zhou et al. 1999; Ganapathisubramani et al. 2005; Volino et al. 2007; Guala et al. 2012; Hutchins et al. 2012). Specifically, in this study, the observed changes in $\bar {R}_{11,l}$ between the CP and PP cases and of $\bar {R}_{11}$ contours during the pulsatile cycle reflect corresponding changes in the geometry of vortex packets in a long-time- and phase-averaged sense. That is, as $\partial \langle \bar {u}_1 \rangle /\partial x_3$ increases, the phase-averaged size of vortex packets is expected to shrink, and, in the long-time-averaged sense, the vortex packets are smaller than their counterparts in the CP case. However, upon inspection of $\bar {R}_{11}$ in figure 11, it is unclear whether the observed change in packet size is attributable to variations in the composing hairpin vortices or the tendency for packets to break into smaller ones under high $\partial \langle \bar {u}_1 \rangle /\partial x_3$ and merge into larger ones under low $\partial \langle \bar {u}_1 \rangle /\partial x_3$. To answer this question, we will next examine the instantaneous turbulence structures and extract characteristic hairpin vortices through conditional averaging. Also, the constant tilting angle of the structure evidenced in figure 12 during the pulsatile cycle indicates that, no matter how vortex packets break and reorganize and how individual hairpin vortices deform in response to the time-varying shear rate, the hairpin vortices within the same packet remain aligned with a constant tilting angle.

3.4. Instantaneous flow structure

Figure 13(a,b) shows the instantaneous fluctuating streamwise velocity $u_1^\prime$ at $x_3/h=1.5$ from the PP case. The chosen phases, $\omega t={\rm \pi} /2$ and $\omega t={\rm \pi}$, correspond to the local minimum and maximum of $\partial \langle \bar {u}_1 \rangle /\partial x_3$, respectively (see figure 6g). Streak patterns can be observed during both phases. As shown in figure 13(a), at low $\partial \langle \bar {u}_1 \rangle /\partial x_3$ values, instantaneous $u_1^\prime$ structures intertwine with neighbouring ones, and form large streaks with a cross-stream extent of approximately $5h$. Conversely, when $\partial \langle \bar {u}_1 \rangle /\partial x_3$ is large, the streaks are shrunk into smaller structures, which have a cross-stream extent of approximately $h$. This behaviour is consistent with the observations we made based on figure 9.

Figure 13. (a,b) Instantaneous fluctuating streamwise velocity $u_1^\prime$ normalized by ${u}_{\tau }$ at $x_3=2h$; (c,d) wall-normal swirl strength $\lambda _{s,3}$ of the PP case at $x_3=2h$; (a,c) $\omega t={\rm \pi} /2$; (b,d), $\omega t={\rm \pi}$. Shaded regions in (c,d) highlight the low-momentum ($u_1^\prime <0$) regions. The instantaneous flow fields correspond to the same pulsatile cycle. Green solid lines highlight the background location of the cuboids.

Further insight into the instantaneous flow field can be gained by considering the low-pass filtered wall-normal swirl strength $\lambda _{s,3}$, shown in figure 13(c,d). The definition of the signed planar swirl strength $\lambda _{s,i}$ is based on the studies of Stanislas, Perret & Foucaut (2008) and Elsinga et al. (2012). The magnitude of $\lambda _{s,i}$ is the absolute value of the imaginary part of the eigenvalue of the reduced velocity gradient tensor $\boldsymbol{\mathsf{J}}_{jk}$, which is

(3.6)\begin{equation} \boldsymbol{\mathsf{J}}_{jk}=\begin{pmatrix} {\partial u_{j}}/{\partial x_{j}} & {\partial u_{j}}/{\partial x_{k}}\\ {\partial u_{k}}/{\partial x_{j}} & {\partial u_{k}}/{\partial x_{k}} \end{pmatrix},\quad i \neq j \neq k , \end{equation}

with no summation over repeated indices. The sign of $\lambda _{s,i}$ is determined by the vorticity component $\omega _i$. Positive and negative $\lambda _{s,i}$ highlight regions with counterclockwise and clockwise swirling motions, respectively. To eliminate the noise from the small-scale vortices, we have adopted the Tomkins & Adrian (2003) idea and low-pass filtered the $\lambda _{s,i}$ field (a compact top-hat filter) with support $h$ to better identify instantaneous hairpin features. As apparent from this figure, low-momentum regions are bordered by pairs of oppositely signed $\lambda _{s,3}$ regions at both the considered phases; these counter-rotating rolls are a signature of hairpin legs. Based on these signatures, it is also apparent that hairpin vortices tend to align in the streamwise direction. Comparing panels (c,d) in figure 13, it is clear that, as $\partial \langle \bar {u}_1 \rangle /\partial x_3$ increases, the swirling strength of the hairpin's legs is intensified, which in turn increases the momentum deficits in the low-momentum regions between the hairpin legs. This behaviour leads to a narrowing of low-momentum regions to satisfy continuity constraints. Also, it is apparent that a larger number of hairpin structures populates the flow field at a higher $\partial \langle \bar {u}_1 \rangle /\partial x_3$, which can be attributed to hairpin vortices spawning offsprings in both the upstream and downstream directions as they intensify (Zhou et al. 1999).

Figure 14 displays a $u_1^\prime$ contour for the PP case at a streamwise/wall-normal plane. Black dashed lines feature a tilting angle $\theta =12^\circ$. It is evident that the interfaces of the low- and high-momentum regions, which are representative instantaneous manifestations of hairpin packets (Hutchins et al. 2012), feature a constant tilting angle during the pulsatile cycle. This behaviour is in agreement with findings from the earlier $\bar {R}_{11}$ analysis, which identified the typical tilting angle of coherent structures as lying between $10^\circ$ and $15^\circ$, depending on the reference wall-normal location. We close this section by noting that, while the instantaneous flow field provides solid qualitative insight into the structure of turbulence for the considered flow field, a more statistically representative picture can be gained by conditionally averaging the flow field on selected instantaneous events. This will be the focus of the next section.

Figure 14. Instantaneous fluctuating streamwise velocity $u_1^\prime$ in a streamwise/wall-normal plane during a pulsatile cycle. Black dashed lines denote the $12^\circ$ structural tilting angle of the coherent structure. Green solid lines represent the canopy layer top.

3.5. Temporal variability of the composite hairpin vortex

This section aims at providing more quantitative insights into the temporal variability of the individual hairpin structures, and elucidating how variations in their geometry influence the ejection-sweep pattern (§ 3.2) and the spatio-temporal coherence of the flow field (§ 3.3). To study the phase-dependent structural characteristics of the hairpin vortex, we utilize the conditional averaging technique (Blackwelder 1977). This technique involves selecting a flow event at a specific spatial location to condition the averaging process in time and/or space. The conditionally averaged flow field is then analysed using standard flow visualization techniques to identify the key features of the eddies involved. By applying this technique to the hairpin vortex, we can gain valuable insights into its structural attributes and how they vary over time.

In the past few decades, various events have been employed as triggers for the conditional averaging operation. For example, in the context of channel flow over aerodynamically smooth surfaces, Zhou et al. (1999) relied on an ejection event as the trigger, which generally coincides with the passage of a hairpin head through that point. More recently, Dennis & Nickels (2011) considered both positive cross-stream and streamwise swirls as triggers, which are indicative of hairpin heads and legs, respectively. In flow over homogeneous vegetation canopies, Watanabe (2004) used a scalar microfront associated with a sweep event. Shortly after, Finnigan et al. (2009) noted that this choice might introduce a bias towards sweep events in the resulting structure and instead used transient peaks in the static pressure, which are associated with both ejection and sweep events.

Here, we adopt the approach first suggested by Coceal et al. (2007), where the local minimum streamwise velocity over a given plane was used as the trigger. It can be shown that this approach yields similar results as the one proposed in Dennis & Nickels (2011) and that it is suitable for the identification of hairpin vortices in the ISL. The conditional averaging procedure used in this study is based on the following operations:

1. (i) Firstly, at a chosen $x_3^e$, we identify the set of locations $(x_1^e,x_2^e)$ where the instantaneous streamwise velocity is $75\,\%$ below its phase-averaged value. This is our ‘triggering event’. Such an operation is repeated for each available velocity snapshot.

2. (ii) Next, for each identified event, the fluctuating velocity field at the selected $x_3^e$ plane is shifted by $(-x_1^e,-x_2^e)$. After this operation, all identified events are located at $(x_1^\prime,x_2^\prime ) = (0,0)$, where $(x_1^\prime,x_2^\prime )$ is the new (translated) coordinate system.

3. (iii) Lastly, the shifted instantaneous velocity fields are averaged over the identified events and snapshots, for each phase.

The end result is a phase-dependent, conditionally averaged velocity field that can be used for further analysis. Figure 15 shows a wall-parallel slice at $x_3/h=2$ of the conditionally averaged fluctuating velocity field in the same plane as the triggering event. Counter-rotating vortices associated with a low-momentum region in between appear to be persistent features of the ISL throughout the pulsation cycle. Vortex cores move downstream and towards each other as $\partial \langle \bar {u}_1 \rangle /\partial x_3$ increases, and the vortices intensify. This behaviour occurs in the normalized time interval $\omega t \in [{\rm \pi} /2, {\rm \pi}]$. Instead, when $\partial \langle \bar {u}_1 \rangle /\partial x_3$ decreases, the cores move upstream and further apart. Such behaviour provides statistical evidence of the behaviour depicted in figure 13(c,d) for the instantaneous flow field. Note that the composite counter-rotating vortex pair in the conditionally averaged flow field is, in fact, an ensemble average of vortex pairs in the instantaneous flow field. Thus, the spacing between the composite vortex pair cores ($d_{\omega }$) represents a suitable metric to quantify the phase-averaged widths of vortex packets in the considered flow system. Figure 16 presents $d_{\omega }$ evaluated with the triggering event at $x_3^e/h=\{1.5, 2, 3, 4 \}$. The trend in $d_{\omega }$ is similar to that observed in figure 9(a) for the first zero crossing of $\bar {R}_{11}$, which is an indicator of the streak width. The explanation for this behaviour is that low-momentum regions are generated between the legs of the hairpins, justifying the observed linear scaling of the streak width with the cross-stream spacing of hairpin legs.

Figure 15. Vector plot of the conditionally averaged fluctuating velocity (PP case) over the $x_3/h=2$ wall-parallel plane. The flow has been conditioned on a local minimum streamwise velocity event in the same plane. Colour contours represent the wall-normal swirling strength $\lambda _{s,3}$. Green dots identify the cores of the counter-rotating vortices.

Figure 16. Spacing between the composite vortex pair cores $d_{\omega }$, corresponding to local minimum streamwise velocity events at $x_3^e/h=1.5$ (black line), $x_3^e/h=2$ (blue line), $x_3^e/h=3$ (green line) and $x_3^e/h=4$ (magenta line).

Figures 17 and 18 depict a conditionally averaged fluctuating velocity field, which is obtained with a triggering event at $x_3^e/h=2$, in the $x^\prime _2=0$ plane and the $x_1^\prime =-h$ plane, respectively. Note that the $x^\prime _2=0$ plane corresponds to the centre plane, and the $x_1^\prime = -h$ cross-section is located $h$ upstream of the triggering event. From figure 17, a region of positive $\lambda _{s,2}$ can be identified immediately above and downstream the location of the triggering event, i.e. $(x_1^\prime,x_2^\prime,x_3^e)=(0,0,2h)$. This $\lambda _{s,2}>0$ region can be interpreted as the head of the composite hairpin vortex (Adrian et al. 2000; Ganapathisubramani, Longmire & Marusic 2003). As $\partial \langle \bar {u}_1 \rangle /\partial x_3$ increases, the vortex structure is deflected downstream and $\lambda _{s,2}$ increases, leading to enhanced upstream ejection events. This behaviour is also apparent from figure 6, where the overall contribution from ejection events to $\langle \overline {u_1^\prime u_3^\prime }\rangle$ increases, while the number of ejection events reduces, highlighting enhanced individual ejection events. The deflection of the hairpin head in the downstream direction is caused by two competing factors. The first is the increase in $\langle \overline {u_1^\prime u_3^\prime }\rangle$, which leads to the downstream deflection. The second factor is the enhancement of the sweep events, which induce an upstream deflection. The first factor outweighs the second, thus yielding the observed variations in the hairpin topology.

Figure 17. Time evolution of the conditionally averaged fluctuating velocity field of the PP case in the streamwise/wall-normal plane $x_2^*/h=0$ given a local minimum streamwise velocity event at $x_3^e/h=2$. Colour contours represent the cross-stream swirling strength $\lambda _{s,2}$. Red and blue lines mark the $\lambda _{s,2}=0.1$ and $\lambda _{s,2}=-0.1$ contours, respectively.

Figure 18. Time evolution of the conditionally averaged fluctuating velocity field in figure 17 in a cross-stream/wall-normal plane $x^\prime _1=-h$. Colour contours represent the streamwise swirling strength $\lambda _{s,1}$. Red and blue lines mark $\lambda _{s,1}=0.1$ and $\lambda _{s,1}=-0.1$, respectively. Green dots identify the cores of the counter-rotating vortices.

Figure 18 shows the response of hairpin legs to changing $\partial \langle \bar {u}_1 \rangle /\partial x_3$ in a cross-stream plane at $x_1^\prime = -h$. A pair of counter-rotating streamwise rollers are readily observed, which, as explained before, identify the legs of the composite hairpin vortex. It also further corroborates our analysis, highlighting that the spacing between the legs reduces from $\approx 5h$ at $\omega t={\rm \pi} /2$ to $\approx 2h$ at $\omega t={\rm \pi}$. This also provides a justification for findings in §§ 3.3 and 3.4. Further, the swirling of the hairpin legs, which is quantified with $\lambda _{s,1}$ and $\lambda _{s,3}$ in the wall-normal/cross-stream and wall-parallel planes, respectively, intensifies with increasing $\partial \langle \bar {u}_1 \rangle /\partial x_3$. Interestingly, when $\partial \langle \bar {u}_1 \rangle /\partial x_3$ approaches its peak value at $\omega t={\rm \pi}$, a modest albeit visible secondary streamwise roller pair is induced by the hairpin legs at $x_2^\prime =\pm 3$. This suggests that the hairpin vortex not only generates new offsprings upstream and downstream, as documented in Zhou et al. (1999) and Adrian (2007), but also in the cross-stream direction when it intensifies. The intensification of hairpin legs creates counter-rotating quasi-streamwise roller pairs between the hairpin vortices adjacent to the cross-stream direction. These roller pairs are lifted up due to the effect of the induced velocity of one roller on the other according to the Biot–Savart law, and the downstream ends of the rollers then connect, forming new hairpin structures.

A more comprehensive picture is provided by isocontours of the conditionally averaged swirling magnitude $\lambda _s = 0.1$ shown in figure 19; $\lambda _s$ is the imaginary part of the complex eigenvalue of the velocity gradient tensor (Zhou et al. 1999). In this case, the conditionally averaged swirling field corresponds to a triggering event at $x_3^e/h=2$. Zhou et al. (1999) pointed out that different thresholds of the iso-surface result in vortex structures of similar shapes but different sizes. The value $\lambda _s=0.1$, in this case, strikes the best compromise between descriptive capabilities and surface smoothness. Note that other vortex identification criteria, such as the Q criterion (Hunt, Wray & Moin 1988) and the $\lambda _2$ criterion (Jeong & Hussain 1995), are expected to result in qualitatively similar vortex structures (Chakraborty, Balachandar & Adrian 2005).

Figure 19. Time evolution of the conditionally averaged swirling field $\lambda _s$ of the PP case given a local minimum streamwise velocity event at $x_3^e=2h$. The shown iso-surfaces are for $\lambda _s=0.1$.

The extents of the conditional eddy in figure 19 vary substantially from roughly $10h\times 8h \times 5h$ at relatively low $\partial \langle \bar {u}_1 \rangle /\partial x_3$ ($\omega t={\rm \pi} /2$), to $6h\times 6h \times 3h$ at high $\partial \langle \bar {u}_1 \rangle /\partial x_3$ ($\omega t={\rm \pi}$). During the period of decreasing $\partial \langle \bar {u}_1 \rangle /\partial x_3$, i.e. $0<\omega t<3/4{\rm \pi}$ and ${\rm \pi} <\omega t<2{\rm \pi}$, the conditional eddy resembles the classic hairpin structure in the stationary case, where two hairpin legs and the hairpin head connecting the hairpin legs can be vividly observed. The sizes of the hairpin legs increase with decreasing $\partial \langle \bar {u}_1 \rangle /\partial x_3$, and so does their spacing, which is in line with our prior observations based on figure 18. One possible physical interpretation for the change in the size of hairpin legs is that the reduction in swirling strength of the hairpin head resulting from a decrease in $\partial \langle \bar {u}_1 \rangle /\partial x_3$ weakens the ejection between the hairpin legs, as shown in figure 17. As a result, the swirling strength of the legs decreases, causing an increase in their size due to the conservation of angular momentum. Conversely, during the period of increasing $\partial \langle \bar {u}_1 \rangle /\partial x_3$ ($3/4{\rm \pi} <\omega t<{\rm \pi}$), the hairpin structure is less pronounced. The conditional eddy features a strengthened hairpin head, and the intensified counter-rotating hairpin legs move closer to each other and ultimately merge into a single region of non-zero swirling strength, as apparent from figure 19. Moreover, downstream of the conditional eddy, a pair of streamwise protrusions, known as ‘tongues’ (Zhou et al. 1999), persist throughout the pulsatile cycle. According to Adrian (2007), these protrusions reflect the early stage of the generation process of the downstream hairpin vortex. These protrusions would eventually grow into a quasi-streamwise vortex pair and later develop a child hairpin vortex downstream of the original one.

In summary, the proposed analysis reveals that the time-varying shear rate resulting from the pulsatile forcing affects the topology and swirling intensity of hairpin vortices. As the shear rate increases (decreases), hairpin vortices tend to shrink (grow) with a corresponding enhancement (relaxation) of the swirling strength. These variations in hairpin geometry are responsible for the observed time-varying ejection-sweep pattern (figure 6). Ejection events primarily occur between the hairpin legs, which become more widely spaced as the vortices grow and less spaced as they shrink. Therefore, a decrease in hairpin vortex size due to an increasing shear rate reduces the number of ejection events, while an increase in vortex size due to the decreasing shear rate leads to an increased number of ejections. Moreover, the intensification (relaxation) of hairpin vortices at high (low) shear rates results in enhanced (attenuated) ejection events between the hairpin legs, as evidenced by figures 17 and 18. This enhancement and attenuation of ejection events is also corroborated by results from figure 6, which indicated that high (low) shear rates decrease (increase) the number of ejection events but increase (decrease) their contribution to $\overline {u_1^\prime u_3^\prime }$. From a flow coherence perspective, this physical process also explains the observed time evolution of $\bar {R}_{11}$ (see figures 9 and 11), which is a statistical signature of hairpin packets. Changes in the size of individual hairpin vortices in response to the shear rate directly influence the dimensions of hairpin packets, as the latter are composed of multiple individual hairpin structures.

4. Conclusions

In this study, the structure of turbulence in pulsatile flow over an array of surface-mounted cuboids was characterized and contrasted with that in stationary flow regimes. The objective was to elucidate the effects of non-stationarity on turbulence topology and its implications for momentum transfer.

Flow unsteadiness was observed to not significantly alter the long-time-average profiles of turbulent kinetic energy and resolved Reynolds shear stress, but it marginally increased the height of the RSL. In the context of quadrant analysis, it was found that flow unsteadiness does not noticeably alter the overall distribution within each quadrant. However, the ejection-sweep pattern exhibited an apparent variation during the pulsation cycle. Flow acceleration yielded a large number of ejection events within the RSL, whereas flow deceleration favoured sweeps. In the ISL, it was shown that the ejection-sweep pattern is mainly controlled by the intrinsic- and phase-averaged shear rate $\partial \langle \bar {u}_1 \rangle /\partial x_3$ rather than by the driving pressure gradient. Specifically, the relative contribution from ejections increases, but their frequency of occurrence decreases with increasing $\partial \langle \bar {u}_1 \rangle /\partial x_3$. The aforementioned time variation in the ejection-sweep pattern was later found to stem from topological variations in the structure of ISL turbulence, as deduced from inspection of the two-point streamwise velocity correlation function and the conditionally averaged flow field.

Specifically, the geometry of hairpin vortex packets, which are the dominant coherent structures in the ISL, was examined through the analysis of two-point velocity correlation to explore its long-time-averaged and phase-dependent characteristics. Flow unsteadiness was found to yield relatively shorter vortex packets in a long-time-average sense (up to 15 % shorter). From a phase-averaged perspective, the three-dimensional extent of hairpin packets was found to vary during the pulsation cycle and to be primarily controlled by $\partial \langle \bar {u}_1 \rangle /\partial x_3$, while their tilting angle remained constant throughout. A visual examination of instantaneous structures also confirmed such behaviour: the size of low-momentum regions and spacing of the hairpin legs encapsulating them were found to change with $\partial \langle \bar {u}_1 \rangle /\partial x_3$, while the hairpin vortices remained aligned at a constant angle during the pulsation cycle.

Further insight into phase variations of instantaneous hairpin structures was later gained using conditional averaging operations, which provided compelling quantitative evidence for the behaviours previously observed. Specifically, the conditional-averaged flow field revealed that the size and swirling intensity of the composite hairpin vortex vary considerably with $\partial \langle \bar {u}_1 \rangle /\partial x_3$. When $\partial \langle \bar {u}_1 \rangle /\partial x_3$ increases to its peak value, the swirling strength of the hairpin head is intensified, yielding strengthened ejections upstream of the hairpin head and a downstream deflection of the hairpin head. As the hairpin head intensifies, there is a corresponding increase in the intensity of the hairpin legs, coupled with a reduction in the spacing between them. This development accounts for the noted decrease in the extent of the ejection-dominated region. In other words, individual ejections become stronger and are generated at a reduced frequency as the shear rate increases, which provides a kinematic interpretation and justification for the observed time variability of the quadrant distribution. Such a process, needless to say, is reversed when the shear rate decreases.

Findings from this study emphasize the significant influence that departures from statistically stationary flow conditions can have on the structure of ABL turbulence and associated processes. Such departures are typical in realistic ABL flows and have garnered growing attention in recent times (Mahrt & Bou-Zeid 2020). While the study focuses on a particular type of non-stationarity, its results underscore the importance of accounting for this flow phenomenon in both geophysical and engineering applications. The modification of turbulence structures due to flow unsteadiness has a substantial effect on exchanges between the land and the atmosphere, as well as on the aerodynamic drag experienced by vehicles. This underlines the necessity for concerted efforts to fully characterize these modifications. From a modelling perspective, empirical insights obtained from this study hold promise for guiding the evolution of more advanced wall-layer model formulations (Piomelli 2008). These models are routinely used in weather and climate forecasting, as well as in aerospace and mechanical engineering applications, facilitating the assessment of area-aggregate exchanges between solid surfaces and the adjacent fluid environment. A recurrent shortcoming of operational wall-layer models lies in their reliance on assumptions of statistical stationarity, overlooking flow unsteadiness and state-dependent turbulence topology information (Monin & Obukhov 1954; Piomelli 2008; Skamarock et al. 2008). This represents an important area for improvement. Past investigations have proposed pathways to integrate turbulence topology information into wall-layer model predictions, leveraging parameters like the vortex packet inclination angle and size (Marusic, Kunkel & Porté-Agel 2001; Marusic, Mathis & Hutchins 2010). These approaches open a fruitful avenue for assimilating the insights derived from this study into wall-layer model infrastructures.