3.1. Overview of flow statistics

Li & Giometto (Reference Li and Giometto2023) 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. (Reference Pokrajac, Campbell, Nikora, Manes and McEwan2007), $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. Reference Coceal, Dobre, Thomas and Belcher2007).

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 (Reference Scotti and Piomelli2001), 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 (Reference Li and Giometto2023). 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 (Reference Wallace, Eckelmann and Brodkey1972), and has thereafter been routinely employed to characterize the structure of turbulence across a range of flow systems (Wallace Reference Wallace2016). 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. Reference Adrian, Meinhart and Tomkins2000); $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. (Reference Coceal, Dobre, Thomas and Belcher2007) and Finnigan (Reference Finnigan2000) 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 (Reference Kanda, Moriwaki and Kasamatsu2004) and Castro et al. (Reference Castro, Cheng and Reynolds2006) 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 Reference Dennis and Nickels2011; Guala et al. Reference Guala, Tomkins, Christensen and Adrian2012). 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 (Reference Zhang and Simons2019) 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 Reference Marusic and Monty2019).

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 (Reference Krogstad and Antonia1994). 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. Reference Takimoto, Inagaki, Kanda, Sato and Michioka2013) and $\bar {R}_{11}=0.5$ (Volino et al. Reference Volino, Schultz and Flack2007; Guala et al. Reference Guala, Tomkins, Christensen and Adrian2012). 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. Reference Chauhan, Hutchins, Monty and Marusic2013; Salesky & Anderson Reference Salesky and Anderson2020). 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. (Reference Coceal, Dobre, Thomas and Belcher2007), 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. (Reference Hutchins, Chauhan, Marusic, Monty and Klewicki2012) and Chauhan et al. (Reference Chauhan, Hutchins, Monty and Marusic2013). 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 Reference Adrian2007) can be used to provide an interpretation for these findings. Note that alternative paradigms, such as that proposed by Del Alamo et al. (Reference Del Alamo, Jimenez, Zandonade and Moser2006), 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. Reference Adrian, Meinhart and Tomkins2000; Christensen & Adrian Reference Christensen and Adrian2001; Dennis & Nickels Reference Dennis and Nickels2011) to ABL field measurements (Hommema & Adrian Reference Hommema and Adrian2003; Morris et al. Reference Morris, Stolpa, Slaboch and Klewicki2007) and numerical simulations (Lee, Sung & Krogstad Reference Lee, Sung and Krogstad2011; Eitel-Amor et al. Reference Eitel-Amor, Örlü, Schlatter and Flores2015). 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. Reference Zhou, Adrian, Balachandar and Kendall1999; Ganapathisubramani et al. Reference Ganapathisubramani, Hutchins, Hambleton, Longmire and Marusic2005; Volino et al. Reference Volino, Schultz and Flack2007; Guala et al. Reference Guala, Tomkins, Christensen and Adrian2012; Hutchins et al. Reference Hutchins, Chauhan, Marusic, Monty and Klewicki2012). 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 (Reference Stanislas, Perret and Foucaut2008) and Elsinga et al. (Reference Elsinga, Poelma, Schröder, Geisler, Scarano and Westerweel2012). 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 (Reference Tomkins and Adrian2003) 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. Reference Zhou, Adrian, Balachandar and Kendall1999).

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. Reference Hutchins, Chauhan, Marusic, Monty and Klewicki2012), 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 Reference Blackwelder1977). 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. (Reference Zhou, Adrian, Balachandar and Kendall1999) 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 (Reference Dennis and Nickels2011) 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 (Reference Watanabe2004) used a scalar microfront associated with a sweep event. Shortly after, Finnigan et al. (Reference Finnigan, Shaw and Patton2009) 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. (Reference Coceal, Dobre, Thomas and Belcher2007), 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 (Reference Dennis and Nickels2011) 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. Reference Adrian, Meinhart and Tomkins2000; Ganapathisubramani, Longmire & Marusic Reference Ganapathisubramani, Longmire and Marusic2003). 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. (Reference Zhou, Adrian, Balachandar and Kendall1999) and Adrian (Reference Adrian2007), 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. Reference Zhou, Adrian, Balachandar and Kendall1999). In this case, the conditionally averaged swirling field corresponds to a triggering event at $x_3^e/h=2$. Zhou et al. (Reference Zhou, Adrian, Balachandar and Kendall1999) 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 Reference Hunt, Wray and Moin1988) and the $\lambda _2$ criterion (Jeong & Hussain Reference Jeong and Hussain1995), are expected to result in qualitatively similar vortex structures (Chakraborty, Balachandar & Adrian Reference Chakraborty, Balachandar and Adrian2005).

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. Reference Zhou, Adrian, Balachandar and Kendall1999), persist throughout the pulsatile cycle. According to Adrian (Reference Adrian2007), 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.