3.2.1. Longtime-averaged velocity profile

Profiles of the longtime-averaged streamwise velocity are shown in figure 5. Flow unsteadiness leads to a horizontal shift of profiles in the proposed semilogarithmic plot. This behaviour is distinct from the ‘two-log’ profile in flow over sand grain roughness (Fredsøe et al. Reference Fredsøe, Andersen and Sumer1999; Yang et al. Reference Yang, Tan, Lim and Zhang2006; Yuan & Madsen Reference Yuan and Madsen2015), and also in stark contrast to the one previously observed in current-dominated pulsatile flow over aerodynamically smooth surfaces, where the longtime-averaged field is essentially unaffected by flow unsteadiness (Tardu & Binder Reference Tardu and Binder1993; Scotti & Piomelli Reference Scotti and Piomelli2001; Tardu Reference Tardu2005; Manna et al. Reference Manna, Vacca and Verzicco2012; Weng et al. Reference Weng, Boij and Hanifi2016). As also apparent from figure 5, departures from the statistically stationary flow profile become more significant for larger values of $\alpha$ and $\omega$.

Figure 5. Wall-normal profiles of longtime-averaged streamwise velocity $\langle \bar {u}_1\rangle _l$ of high-amplitude cases (a) and low-amplitude cases (b). Line colour specifies the forcing frequency: navy blue, $\omega T_{h}=0.05{\rm \pi}$; dark green, $\omega T_{h}=0.125{\rm \pi}$; light green, $\omega T_{h}=0.25{\rm \pi}$; yellow-green, $\omega T_{h}=1.25{\rm \pi}$. Here $\langle \bar {u}_1\rangle$ from the SS case, represented by the red dashed line, is included for comparison. Black solid line indicates the slope of $1/\kappa$, with $\kappa =0.4$. The shaded area highlights the fitting region for the estimation of the roughness length scale $z_0$.

Variations in the aerodynamic roughness length $z_0$ and displacement height $d$ parameters with $\omega T_h$ are shown in figure 6 for the considered canopy. These parameters are evaluated via the Macdonald et al. (Reference Macdonald, Griffiths and Hall1998) approach, where $d$ is the barycentre height of the longtime-averaged pressure drag from the urban canopy, and $z_0$ is determined via curve fitting. More specifically,

(3.2)\begin{equation} d = \frac{\int^h_0 \langle \bar{D} \rangle_l (x_3)x_3 \,{\rm d}\kern0.7pt x_3 }{\int^h_0 \langle \bar{D} \rangle_l (x_3 )\,{\rm d}\kern0.7pt x_3 }, \end{equation}

where the wall-normal distribution of the instantaneous canopy pressure drag $D$ is obtained by taking an intrinsic volume average of the pressure gradient, i.e.

(3.3)\begin{equation} D(x_3,t) =\frac{1}{V_f}\int_{x_3-\delta z/2}^{x_3+\delta z/2}\int_0^{L_2} \int_0^{L_1}\frac{1}{\rho}\frac{\partial p}{\partial x} \,{\rm d}\kern0.7pt x_1 \,{\rm d}\kern0.7pt x_2 \,{\rm d}\kern0.7pt x_3. \end{equation}

Figure 6. Normalized aerodynamic roughness length $z_0$ (a) and displacement height $d$ (b). Different colours correspond to different oscillation amplitudes: blue, $\alpha =0.2$; red, $\alpha =0.4$. The black square symbol denotes the reference SS case.

Note that, in principle, one should also account for the SGS drag contribution in (3.3); in this work, we omit SGS contributions because they are negligible when compared with the total drag – a direct result of the relatively small aerodynamic roughness length that is prescribed in the wall-layer model (see § 2.1).

Here $z_0$ is solved by minimizing the root mean square error between the longtime-averaged velocity and the law of the wall with $\kappa =0.4$ in the $x_3 \in [2h,6h]$ interval, i.e.

(3.4)\begin{equation} E=\left\| \langle \bar{u}_1\rangle _{l}-\frac{u_{\tau}}{\kappa}\log\left(\frac{x_3-d}{z_0}\right)\right\|_{2} . \end{equation}

The fitting interval is highlighted in figure 5. The estimated $z_0$ was found to be poorly sensitive to variations in the fitting interval within the considered range of values. Cheng et al. (Reference Cheng, Hayden, Robins and Castro2007) argued that MacDonald's method is accurate when surfaces are characterized by a low packing density – a requirement that is indeed satisfied in the considered cases.

As apparent from figure 6, $d$ is poorly sensitive to variations in both $\alpha$ and $\omega$ (variations across cases are within the $\pm 3\,\%$ range). This behaviour can be explained by considering that for flows over sharp-edged obstacles, such as the ones considered herein, flow separation consistently takes place at the sharp edges, irrespective of variations in $\alpha$ and $\omega$. This results in a similar net pressure distribution on the faces of the cuboids across all the considered cases. This can be readily observed in figure 4. The $d$ parameter is here evaluated as the integral of the pressure gradient field over the surface area, so the above considerations provide a physical justification for the observed behaviour. Note that this finding might not be generalizable across all possible roughness morphologies. For instance, Yu et al. (Reference Yu, Rosman and Hench2022) showed that separation patterns in pulsatile flows over hemispheres feature a rather strong dependence on $\alpha$ and $\omega$, yielding corresponding strong variations in $d$.

Contrary to $d$, the $z_0$ parameter is strongly impacted by flow unsteadiness, and its value increases with $\alpha$ and $\omega$. Bhaganagar (Reference Bhaganagar2008) reported a similar upward shift of velocity profile in his simulations of pulsatile flow over transitionally rough surfaces at a low Reynolds number. She attributed the increase in $z_0$ to the resonance between the unsteady forcing and the vortices shed by roughness elements, which is induced when the forcing frequency approaches that of the vortex shedding. However, such an argument does not apply to the cases under investigation, since we observed no spurious peaks in the temporal streamwise velocity spectrum. Rather, the increase in $z_0$ stems from the quadratic relation between the phase-averaged canopy drag and velocity, as elaborated below.

3.2.2. Phase-averaged drag-velocity relation

The phase-averaged canopy drag $\langle \bar {D}\rangle$ and the local phase- and intrinsic-averaged velocity $\langle \bar {u}_1 \rangle$ at $x_3/h\approx 0.8$ are shown in figure 7, and are representative of corresponding quantities at different heights within the UCL. Results from cases with three lower frequencies and that from the SS case cluster along a single curve, highlighting the presence of a frequency-independent one-to-one mapping between $\langle \bar {D}\rangle$ and $\langle \bar {u}_1 \rangle$. As apparent from figure 8(a), at these three forcing frequencies, the interaction between the wind and the canopy layer is in a state of quasiequilibrium, i.e. $\langle \bar {D}\rangle$ is in phase with $\langle \bar {u}_1 \rangle$. Moreover, the shape of the aforementioned curve generally resembles the well-known quadratic drag law, which is routinely used to parameterize the surface drag in reduced-order models for stationary flow over plant and urban canopies (Lettau Reference Lettau1969; Raupach Reference Raupach1992; Macdonald et al. Reference Macdonald, Griffiths and Hall1998; Coceal & Belcher Reference Coceal and Belcher2004; Katul et al. Reference Katul, Mahrt, Poggi and Sanz2004; Poggi, Katul & Albertson Reference Poggi, Katul and Albertson2004). This finding comes as no surprise, given the quasiequilibrium state of the $\langle \bar {D}\rangle -\langle \bar {u}_1 \rangle$ relation for the three lower forcing frequencies.

Figure 7. Phase-averaged canopy drag $\langle \bar {D} \rangle$ at $x_3/h\approx 0.8$ as a function of the local phase- and intrinsic-averaged velocity $\langle \bar {u}_1 \rangle$ (a) and the same $\langle \bar {D} \rangle -\langle \bar {u}_1 \rangle$ plot but with the time lag between $\langle \bar {D} \rangle$ and $\langle \bar {u}_1 \rangle$ removed (b). Different colours are used to denote different forcing frequencies: navy blue, $\omega T_{h}=0.05{\rm \pi}$; dark green, $\omega T_{h}=0.125{\rm \pi}$; light green, $\omega T_{h}=0.25{\rm \pi}$; yellow-green, $\omega T_{h}=1.25{\rm \pi}$. Different symbols are used to distinguish between oscillation amplitudes: triangle, $\alpha =0.2$; cross, $\alpha =0.4$. Red square represents the SS case. Red dashed line represents $\langle \bar {D} \rangle =C_d\lambda _f \langle \bar {u}_1 \rangle | \langle \bar {u}_1 \rangle |$ with $C_d$ obtained from the SS case.

Figure 8. Time evolution of $\langle \bar {D} \rangle$ (black solid lines) and $\langle \bar {u}_1 \rangle$ (red dashed lines) at $x_3/h\approx 0.8$ from the HL (a) and LVH (b) cases. The acronyms of LES runs are defined in table 1.

On the other hand, as shown in figure 8(b), results from the highest-frequency cases (LVH and HVH) exhibit an orbital pattern, which stems from the time lag between $\langle \bar {D} \rangle$ and $\langle \bar {u}_1 \rangle$. Artificially removing this time lag indeed yields a better clustering of data points from the highest-frequency cases along the quadratic drag law (see figure 7b).

These findings suggest the following parameterization for $\langle \bar {D} \rangle$:

(3.5)\begin{equation} \langle \bar{D} \rangle(x_3,t)= C_d(x_3) \lambda_f \langle \bar{u}_1 \rangle \left | \langle \bar{u}_1 \rangle \right |(x_3,t+\Delta t) , \end{equation}

where $C_d$ is a sectional drag coefficient that is constant in time and does not depend on $\alpha$ or $\omega$, and $\Delta t$ accounts for the time lag between $\langle \bar {D} \rangle$ and $\langle \bar {u}_1 \rangle$, which instead does indeed depend on $\alpha$ and $\omega$. Note that, throughout the considered cases, the wall-normal-averaged drag coefficient $\int ^{h}_0 C_d \,{\rm d}\kern0.7pt x_3/h\approx 0.9$ – a value that is similar to those previously reported ones for stationary flow over cube arrays (Coceal & Belcher Reference Coceal and Belcher2004) (note that the exact value depends on the formula used to define $C_d$).

Morison, Johnson & Schaaf (Reference Morison, Johnson and Schaaf1950) developed a semiempirical model relating the phase-averaged drag generated by obstacles in an oscillatory boundary layer to a given phase- and intrinsic-averaged velocity – a model that has been extensively used in the ocean engineering community to evaluate drag from surface-mounted obstacles (Lowe, Koseff & Monismith Reference Lowe, Koseff and Monismith2005; Yu, Rosman & Hench Reference Yu, Rosman and Hench2018; Yu et al. Reference Yu, Rosman and Hench2022). The Morison model assumes that the total force applied to the fluid by obstacles consists of a quadratic drag term and an inertial term; the latter accounts for the added mass effect and the Froude–Krylov force arising as a direct consequence of the unsteady pressure field. While it might seem plausible to utilize the Morison model to evaluate surface drag as a function of phase-averaged velocity at a given $x_3$ for the cases under consideration, unfortunately, this approach is not applicable. As shown in Patel & Witz (Reference Patel and Witz2013), the Morison model provides relatively accurate evaluations of obstacle drag when the phased-averaged acceleration at different $x_3$ are in phase. This is not the case in this study, where substantial phase lags between phase-averaged accelerations at different wall-normal locations substantially degrade the accuracy of such a model. This behaviour can be easily inferred from the results in § 3.3.

In the following, we will make use of (3.5) to derive an alternative phenomenological surface-drag model for the considered flow system.

3.2.3. Mapping roughness length variability to longtime-averaged flow statistics

Here $z_0$ and $d$ are input parameters of surface flux parameterizations that are routinely used in numerical weather prediction, climate projection and pollutant dispersion models (see, e.g. Skamarock et al. Reference Skamarock, Klemp, Dudhia, Gill, Barker, Duda, Huang, Wang and Powers2008; Benjamin et al. Reference Benjamin2016). These models are typically based on Reynolds-averaged Navier–Stokes closures, and feature time steps that can go from one hour up to several days. When departures from stationarity occur at a time scale that is much smaller than the time step of the model, model predictions are essentially longtime-averaged quantities, and the validity of surface flux parameterizations based on flow homogeneity and stationarity assumptions may break down. As a first step towards addressing this problem, this section proposes a phenomenological model relating $z_0$ to longtime-averaged pulsatile-flow statistics.

The longtime-averaged friction velocity can be written as

(3.6)\begin{equation} u^2_{\tau}=\int^{h}_0 \langle \bar{D} \rangle_l \,{\rm d}\kern0.7pt x_3=\int^{h}_0 C_d \lambda_f \langle \overline{\langle \bar{u}_1 \rangle \left | \langle \bar{u}_1 \rangle \right |}\rangle_l\,{\rm d}\kern0.7pt x_3 , \end{equation}

where $\langle \bar {D} \rangle _l$ is the longtime-averaged surface drag. Note that the wall-normal structure of $C_d$ is approximately constant in the UCL (not shown here), except in the vicinity of the surface, where local contributions to the overall drag are, however, minimal due to the small value of $\langle \bar {u}_1 \rangle$. Thus, it is reasonable to assume $C_d$ is constant along the wall-normal direction. Also, depending on $\alpha$ and $\omega$, the flow within the canopy might undergo a local reversal in the phase- and intrinsic-averaged sense, meaning that $\langle \bar {u}_1 \rangle < 0$ at selected $x_3$ locations. Assuming that there is no flow reversal within the UCL, i.e. $\langle \bar {u}_1 \rangle \ge 0$ in $x_3 \le h$, (3.6) can be written as

(3.7)\begin{equation} u_{\tau}^2= C_d \lambda_f \left( \int^{h}_{0}{\langle \bar{u}_1\rangle_l^2}\,{{\rm d}\kern0.7pt x_3} + \int^{h}_{0}{\langle\overline{\langle \bar{u}_1\rangle_o ^2}\rangle_l}\,{{\rm d}\kern0.7pt x_3} \right) , \end{equation}

where the second term on the right-hand side of (3.7) is identically zero for the SS case. Equation (3.7) essentially states that an unsteady canopy layer requires a lower longtime-averaged wind speed to generate the same drag of a steady canopy layer since quadratic drag contributions are generated by flow unsteadiness (the second term on the right-hand side of (3.7)). Note that $\sqrt {\int ^{h}_0 ({\cdot })^2\, {\rm d}\kern0.7pt x_3/h }$ is an averaging operation over the UCL based on the $L_2$ norm. Rearranging terms in (3.7) leads to

(3.8)\begin{equation} \langle \bar{u}_1\rangle_{l,{avg}} = \sqrt{\frac{u^2_{\tau}}{C_d \lambda_f h} -\frac{1}{h} \int^{h}_{0}{\langle\overline{\langle \bar{u}_1\rangle_o ^2}\rangle_l}\,{{\rm d}\kern0.7pt x_3} } , \end{equation}

and

(3.9)\begin{equation} \langle \bar{u}_1\rangle_{l,{avg}}^{{SS}} = \sqrt{\frac{u^2_{\tau}}{C_d \lambda_f h}} , \end{equation}

for the pulsatile cases and the SS case, respectively. Here $({\cdot })_{{avg}}$ denotes the canopy-averaged quantity, and $({\cdot })^{{SS}}$ represents a quantity pertaining to the SS case.

As discussed in § 3.2.1, flow unsteadiness yields a shift of the $\langle \bar {u}_1\rangle _{l}$ profile with negligible variations in the $d$ parameters when compared with the stationary flow with the same $u_\tau$. In terms of the law-of-the-wall, this behaviour can be described as a variation in $z_0$, i.e.

(3.10)\begin{equation} \langle \bar{u}_1\rangle_{l,{avg}}^{{SS}}-\langle \bar{u}_1\rangle_{l,{avg}}=\frac{u_\tau}{\kappa}\log \left( \frac{x_3-d}{z_0^{{SS}}}\right)-\frac{u_\tau}{\kappa}\log \left( \frac{x_3-d}{z_0}\right) , \end{equation}

where (3.10) is valid for any $x_3$ in the logarithmic region. The shift in the velocity profile is approximately constant for $x_3 \in [0,L_3]$, so one can write

(3.11)\begin{equation} \langle \bar{u}_1\rangle_l^{{SS}}-\langle \bar{u}_1\rangle_l \approx \langle \bar{u}_1\rangle_{l,{avg}}^{{SS}}-\langle \bar{u}_1\rangle_{l,{avg}} , \end{equation}

and substituting (3.8)–(3.10) into (3.11) finally yields

(3.12)\begin{equation} z_0=z_0^{{SS}} \exp \left[ \kappa \left(\sqrt{\frac{1}{C_d \lambda_f h}}-\sqrt{\frac{1}{C_d \lambda_f h}-\frac{\int^{h}_{0}{\langle\overline{\langle \bar{u}_1\rangle_o ^2}\rangle_l}\,{{\rm d}\kern0.7pt x_3}/h}{u^2_{\tau}}}\right)\right] . \end{equation}

Equation (3.12) is a diagnostic model relating variations in the $z_0$ parameter to the UCL phase- and intrinsic-averaged velocity variance – a longtime-averaged quantity. The $z_0$ estimates from (3.12) are compared with LES results in figure 9, using $C_d = 0.9$. It is apparent that the proposed model is able to accurately evaluate $z_0$ for most of the considered cases. For the LVH, HH and HVW runs, $z_0$ is overestimated by the model; these departures are attributed to the presence of flow reversal in the UCL, which contradicts the model assumptions.

Figure 9. Comparison between $z_0$ estimated via (3.12) and $z_0$ from LES. Symbols and colours correspond to those used in figure 7.

Equation (3.12) highlights that, in the absence of flow reversal, $z_0$ can be described as a monotonically increasing function of the $\langle \bar {u}_1\rangle _o$ variance in the UCL. As explained at the beginning of this section, this finding is of high relevance from a flow modelling perspective, because it relates a longtime-averaged flow statistic to the $z_0$ parameter. Note that $z_0^{{SS}}$ can be accurately evaluated using any existing parameterization for stationary ABL flow over aerodynamically rough surfaces, including the Lettau (Reference Lettau1969), Raupach (Reference Raupach1992) and Macdonald et al. (Reference Macdonald, Griffiths and Hall1998) models. Further, $C_d = 0.9$ (see discussion in § 3.2.2) and $\lambda _f$ and $h$ are morphological parameters that are in general a priori available. In weather forecasting and climate models, $\langle \bar {u}_1\rangle _o$ is an SGS quantity and would hence have to be parameterized as a function of longtime-averaged statistics that the model computes or from available in-situ or remote sensing measurements.

3.2.4. Longtime-averaged resolved Reynolds stress

This section shifts the attention to longtime-averaged resolved Reynolds stresses. For all of the considered cases, contributions from SGS stresses account for $<1\,\%$ of the total phase-averaged Reynolds stresses and are hence not discussed.

Throughout the boundary layer, $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ profiles are indistinguishable from the SS one (not shown), indicating a weak dependence of such a quantity on $\alpha$ and $\omega$. Above the UCL, the divergence of $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ balances the longtime-averaged driving pressure gradient $f_m$, i.e.

(3.13)\begin{equation} -\frac{\partial \langle\overline{u_1^\prime u_3^\prime}\rangle_l}{\partial x_3}+f_m=0 . \end{equation}

Since $f_m$ does not vary across the considered cases and $\langle \overline {u_1^\prime u_3^\prime }\rangle _l(L_3) = 0$, a collapse of $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ profiles in this region was to be expected from the mathematical structure of the governing equations. Within the UCL, the longtime-averaged momentum budget reads

(3.14)\begin{equation} -\frac{\partial \langle\overline{u_1^\prime u_3^\prime}\rangle_l}{\partial x_3}+f_m-\langle\bar{D}\rangle_l=0 . \end{equation}

Given that $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ profiles collapse in this region, $\langle \bar {D}\rangle _l$ is also expected to feature a weak dependence on $\alpha$ and $\omega$, which in turn explains the weak variations in the $d$ parameter that were observed in § 3.2.1.

Figure 10 depicts wall-normal profiles of the longtime-averaged resolved turbulent kinetic energy, which is defined as

(3.15)\begin{equation} \langle\bar{k}\rangle_l=\tfrac{1}{2}(\langle\overline{u_1^\prime u_1^\prime}\rangle_l+\langle\overline{u_2^\prime u_2^\prime}\rangle_l+\langle \overline{u_3^\prime u_3^\prime}\rangle_l) . \end{equation}

Such a quantity features a significant increase in the UCL when compared with its values farther away from the surface, which as discussed in Schmid et al. (Reference Schmid, Lawrence, Parlange and Giometto2019), is due to dispersive contributions caused by the canopy geometry. Flow unsteadiness has a relatively more important impact on such a quantity when compared with the $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ component, especially for flow in the UCL and for cases with a high oscillation amplitude. An increase in $\alpha$ generally results in more pronounced departures of $\langle \bar {k}\rangle _l$ profiles from the SS one. This behaviour can be best explained by considering the shear production terms in the budget equation for $\langle \bar {k}\rangle _l$, i.e.

(3.16)\begin{equation} \langle \bar{\mathcal{P}}\rangle_l={-}\left\langle\overline{\langle \overline{u_1^\prime u_3^\prime} \rangle \frac{\partial \langle \bar{u}_1 \rangle }{\partial x_3}}\right\rangle_l=\underbrace{-\langle \overline{ u_1^\prime u_3^\prime}\rangle_l \frac{\partial\langle \bar{u}_1\rangle_l }{\partial x_3}}_{\langle \bar{\mathcal{P}}\rangle_{l,1}}\underbrace{-\left\langle\overline{\langle\overline{ u_1^\prime u_3^\prime }\rangle_o \frac{\partial \langle\bar{u}_1\rangle_o }{\partial x_3}}\right\rangle_l} _{\langle \bar{\mathcal{P}}\rangle_{l,2}} , \end{equation}

where $\langle \bar {\mathcal {P}}\rangle _{l,1}$ represents the work done by $\langle \overline { u_1^\prime u_3^\prime }\rangle _l$ onto the longtime-averaged flow field, and ${\langle \bar {\mathcal {P}}\rangle _{l,2}}$ is the longtime average of the work done by the oscillatory shear stress $\langle \overline { u_1^\prime u_3^\prime }\rangle _o$ onto the oscillatory flow field. Figure 11(a) shows that $\langle \bar {\mathcal {P}}\rangle _{l,1}$ is poorly sensitive to variations in $\alpha$ and $\omega$. This behaviour stems from the constancy of $\langle \overline { u_1^\prime u_3^\prime }\rangle _l$ and ${\partial \langle \bar {u}_1\rangle _l }/{\partial x_3}$ across cases (the latter can be inferred from the systematic shift of $\langle \bar {u}_1\rangle _l$ profiles in figure 5). Conversely, $\langle \bar {\mathcal {P}}\rangle _{l,2}$ from high-amplitude cases are generally larger than those from low-amplitude ones, mainly due to the higher $\langle \overline { u_1^\prime u_3^\prime }\rangle _o$ and $\langle \bar {u}_1\rangle _o$ values. Discrepancies in $\langle \bar {\mathcal {P}}\rangle _{l,2}$ among high-amplitude cases are larger than those among low-amplitude ones, which ultimately yields the observed variability in $\langle \bar {k}\rangle _l$.

Figure 10. Longtime-averaged resolved turbulent kinetic energy $\langle \bar {k}\rangle _l=(\langle \overline {u_1^\prime u_1^\prime }\rangle _l+\langle \overline {u_2^\prime u_2^\prime }\rangle _l+\langle \overline {u_3^\prime u_3^\prime }\rangle _l)/2$ from high-amplitude cases (a) and low-amplitude cases (b). Line colours correspond to those used in figure 5.

Figure 11. Normalized shear production terms of $\langle \bar {k}\rangle _l$: $\langle \bar {\mathcal {P}}\rangle _{l,1}$ (a) and $\langle \bar {\mathcal {P}}\rangle _{l,2}$ (b). Symbols and line colours correspond to those used in figure 7.

Further insight into the problem can be gained by looking at the normal components of the longtime-averaged resolved Reynolds stress tensor, which are shown in figure 12. In this case, it is apparent that increases in the oscillation frequency lead to a decrease in $\langle \overline {u_1^\prime u_1^\prime }\rangle _l$ and an increase in $\langle \overline {u_2^\prime u_2^\prime }\rangle _l$ and $\langle \overline {u_3^\prime u_3^\prime }\rangle _l$ within the UCL – a behaviour that is especially apparent for the high-amplitude cases (figure 12d–f). These trends can be best understood by examining the pressure-strain terms from the budget equations of the longtime-averaged resolved Reynolds stresses, i.e.

(3.17)\begin{equation} {\mathsf{R}}_{ij} = \left\langle\overline{\frac{p^\prime}{\rho} \left(\frac{\partial u^\prime_j}{\partial x_i}+\frac{\partial u^\prime_i}{\partial x_j}\right)}\right\rangle_l. \end{equation}

Figure 12. Longtime-averaged resolved normal Reynolds stresses $\langle \overline {u_1^\prime u_1^\prime }\rangle _l$, $\langle \overline {u_2^\prime u_2^\prime }\rangle _l$ and $\langle \overline {u_3^\prime u_3^\prime }\rangle _l$ from low-amplitude cases (a–c) and high-amplitude cases (d–f). Line colours correspond to those used in figure 5.

Here $R_{ii}$ are responsible for redistributing kinetic energy among the longtime-averaged normal Reynolds stresses (Pope Reference Pope2000), and are shown in figure 13. With the exception of the very near surface region ($x_3 \lessapprox 0.2$) where no clear trend can be observed, increases in $\omega$ and $\alpha$ yield a decrease in ${\mathsf{R}}_{11}$ and an increase in ${\mathsf{R}}_{22}$ and ${\mathsf{R}}_{33}$, which justify the observed isotropization of turbulence in the UCL.

Figure 13. Pressure redistribution terms ${\mathsf{R}}_{ii}$: (a–c) low-amplitude cases; (d–f) high-amplitude cases. Line colours correspond to those used in figure 5.

3.3.1. Oscillation amplitude impacts on the oscillatory fields

Flow statistics from the LL and HL simulations are here examined to study the impact of the oscillation amplitude $(\alpha )$ on the oscillatory velocity and resolved Reynolds stresses. Only the LL and HL runs are discussed, as they were found to be representative of the observed behaviours for the other low and high-amplitude cases, respectively.

Figure 14 contrasts the profiles of oscillatory velocity $(\langle \bar {u}_1\rangle _o)$ from the LL and HL runs at four equispaced phases during the pulsatile cycle. The $\langle \bar {u}_1\rangle _o$ at the top of the domain is controlled by the $\alpha$ parameter, so it is natural to use $u_{\tau }\alpha$ as a normalization constant to study the problem. Manna et al. (Reference Manna, Vacca and Verzicco2012) investigated pulsatile open channel flow over an aerodynamically smooth surface, and showed that using such a normalization constant is indeed convenient as it leads to a collapse of $\langle \bar {u}_1\rangle _o$ profiles across cases with different amplitudes, even in the presence of strong flow reversal. This indicates that, at a given forcing frequency, the amplitude of the oscillatory velocity within the domain is proportional to that at the top of the domain. In this work, we show that such a scaling works well also in the presence of aerodynamically rough surfaces, as evidenced by the excellent collapse of $\langle \bar {u}_1\rangle _o/(u_{\tau }\alpha )$ profiles in figure 14. This is not trivial, especially when considering the different scaling of surface drag between aerodynamically smooth and rough walls in the presence of flow unsteadiness.

Figure 14. Profiles of the oscillatory velocity $\langle \bar {u}_1\rangle _o$ from the LL (blue) and HL (red) cases at different phases of the pulsatile cycle. Profiles are $T/4$ apart.

The oscillatory resolved turbulent kinetic energy can be defined as

(3.18)\begin{equation} \langle\bar{k}\rangle_o=\tfrac{1}{2}\big(\langle\overline{u_1^\prime u_1^\prime}\rangle_o+\langle\overline{u_2^\prime u_2^\prime}\rangle_o+\langle\overline{u_3^\prime u_3^\prime}\rangle_o\big) . \end{equation}

As shown in figures 15 and 16, both the oscillatory resolved Reynolds shear stress $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and $\langle \bar {k}\rangle _o$ scale with $u_{\tau }^2 \alpha$. Although not shown, the three oscillatory normal Reynolds stresses also obey such a scaling, suggesting that any change in $\langle \bar {k}\rangle _o$ is proportionally distributed to the three normal Reynolds stresses during the pulsatile cycle. As discussed in the following paragraphs, the scaling of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and $\langle \bar {k}\rangle _o$ is a direct consequence of the mild nonlinearity in the production of $\langle \bar {k}\rangle _o$ and $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$. Subtracting the budget equation of $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ from that of $\langle \overline {u_1^\prime u_3^\prime } \rangle$, one obtains

(3.19)\begin{align} \frac{\partial \langle\overline{u_1^\prime u_3^\prime}\rangle_o}{\partial t}= & \underbrace{-2\langle\overline{u_3^\prime u_3^\prime}\rangle_o \frac{\partial \langle\bar{u}_1 \rangle_l}{\partial x_3}}_{\langle \bar{\mathcal{P}}_{13}\rangle_{o,l1}} \underbrace{-2\langle\overline{u_3^\prime u_3^\prime}\rangle_l \frac{\partial \langle\bar{u}_1 \rangle_o}{\partial x_3}}_{\langle \bar{\mathcal{P}}_{13}\rangle_{o,l2}}\nonumber\\ & \underbrace{-2\langle\overline{u_3^\prime u_3^\prime}\rangle_o \frac{\partial \langle\bar{u}_1 \rangle_o }{\partial x_3}+2\left\langle\overline{\langle\overline{u_3^\prime u_3^\prime}\rangle_o \frac{\partial \langle\bar{u}_1 \rangle_o }{\partial x_3}}\right\rangle_l}_{\langle \bar{\mathcal{P}}_{13}\rangle_{o,nl}}+\cdots, \end{align}

where $\langle \bar {\mathcal {P}}_{13}\rangle _{o,l1}$ and $\langle \bar {\mathcal {P}}_{13}\rangle _{o,l2}$ are the linear production terms of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$, while $\langle \bar {\mathcal {P}}_{13}\rangle _{o,nl}$ is the nonlinear production. Similarly, the budget equations of $\langle \bar {k}\rangle _o$ can be written as

(3.20)\begin{equation} \frac{\partial \langle\bar{k}\rangle_o}{\partial t}= \underbrace{-\langle\overline{u_1^\prime u_3^\prime}\rangle_o \frac{\partial \langle\bar{u}_1 \rangle_l}{\partial x_3}}_{\langle \bar{\mathcal{P}}_{k}\rangle_{o,l1}} \underbrace{-\langle\overline{u_1^\prime u_3^\prime}\rangle_l \frac{\partial \langle\bar{u}_1 \rangle_o}{\partial x_3}}_{\langle \bar{\mathcal{P}}_{k}\rangle_{o,l2}} \underbrace{-\langle\overline{u_1^\prime u_3^\prime}\rangle_o \frac{\partial \langle\bar{u}_1 \rangle_o }{\partial x_3}+\left\langle\overline{\langle\overline{u_1^\prime u_3^\prime}\rangle_o \frac{\partial \langle \bar{u}_1 \rangle_o }{\partial x_3}}\right\rangle_l}_{\langle \bar{\mathcal{P}}_{k}\rangle_{o,nl}}+\cdots , \end{equation}

where $\langle \bar {\mathcal {P}}_{k}\rangle _{o,l1}$ and $\langle \bar {\mathcal {P}}_{k}\rangle _{o,l2}$ are the linear production terms of $\langle \bar {k}\rangle _o$, while $\langle \bar {\mathcal {P}}_{k}\rangle _{o,nl}$ is the nonlinear production. As shown in figure 17, the nonlinear production term is substantially smaller than the sum of the corresponding linear productions for both $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and $\langle \bar {k}\rangle _o$. Given that $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$, $\langle \overline {u_3^\prime u_3^\prime }\rangle _l$ and ${\partial \langle \bar {u}_1\rangle _l}/{\partial x_3}$ from the LL and HL cases are similar, when $\langle \bar {u}_1\rangle _o \sim u_{\tau }\alpha$ and $\langle \overline {u_3^\prime u_3^\prime }\rangle _o\sim u_{\tau }^2 \alpha$, the total production of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ is

(3.21)\begin{equation} \langle \bar{\mathcal{P}}_{13}\rangle_{o}\approx \langle \bar{\mathcal{P}}_{13}\rangle_{o,l1}+\langle \bar{\mathcal{P}}_{13}\rangle_{o,l2} \sim \frac{u_{\tau}^3 \alpha}{h} , \end{equation}

whereas that of $\langle \bar {k}\rangle _o$ is

(3.22)\begin{equation} \langle \bar{\mathcal{P}}_{k}\rangle_{o}\approx \langle \bar{\mathcal{P}}_{k}\rangle_{o,l1}+\langle \bar{\mathcal{P}}_{k}\rangle_{o,l2} \sim \frac{u_{\tau}^3 \alpha}{h} . \end{equation}

This in turn leads to the observed $\langle \overline {u_1^\prime u_3^\prime }\rangle _o \sim u_{\tau }^2 \alpha$ and $\langle \bar {k}\rangle _o \sim u_{\tau }^2 \alpha$ scalings.

Figure 15. Profiles of the oscillatory resolved Reynolds shear stress $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ from the LL (blue) and HL (red) cases at different phases of the pulsatile cycle. Profiles are $T/4$ apart.

Figure 16. Profiles of the oscillatory resolved turbulent kinetic energy $\langle \bar {k}\rangle _o=(\langle \overline {u_1^\prime u_1^\prime }\rangle _o+\langle \overline {u_2^\prime u_2^\prime }\rangle _o+\langle \overline {u_3^\prime u_3^\prime }\rangle _o)/2$ from the LL (blue) and HL (red) cases at different phases of the pulsatile cycle. Profiles are $T/4$ apart.

Figure 17. Normalized production terms for $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ (a) and for $\langle \bar {k}\rangle _o$ (b) at $x_3/h=1.5$ from the LL (blue) and HL (red) cases. Solid lines denote the linear production terms, and dash lines represent the nonlinear production terms.

Equations (3.21) and (3.22) are expected to fail under two conditions. Firstly, when $\alpha$ is sufficiently large and the contribution of $\langle \bar {\mathcal {P}}_{13}\rangle _{o,nl}$ and $\langle \bar {\mathcal {P}}_{k}\rangle _{o,nl}$ can no longer be neglected. Secondly, when differences in $\langle \overline {u_1^\prime u_3^\prime }\rangle _l$ or $\langle \overline {u_3^\prime u_3^\prime }\rangle _l$ among cases with different $\alpha$ become large enough that the linear production terms cease to scale with $u_{\tau }^2 \alpha$, which occurs within the UCL for the LVH and HVH cases (see figure 12).

The above analysis has shown that the $\alpha$ parameter primarily controls the amplitude of selected oscillatory flow quantities, but has little impact on their wall-normal structure. In the next section, it will be shown that the wall-normal structure of quantities is instead controlled by $\omega$. This behaviour is also expected to hold in the smooth-wall set-up (Manna et al. Reference Manna, Vacca and Verzicco2012), although the mechanism responsible for generating drag over aerodynamically smooth surfaces is quite distinct.

3.3.2. Forcing frequency impacts on the oscillatory fields

This section discusses how the oscillatory velocity and resolved Reynolds stresses respond to variations in the forcing frequency. Only low-amplitude cases will be considered since conclusions can be generalized across the considered runs.

Figures 18 and 19 present the time evolution of $\langle \bar {u}_1\rangle _o$ and $\partial \langle \bar {u}_1\rangle _o/\partial x_3$. Three distinct frequency regimes can be identified. The first regime corresponds to the highest amongst the considered forcing frequencies, i.e. the LVH case. For this flow regime, the oscillation in $\partial \langle \bar {u}_1\rangle _o/\partial x_3$ is typically confined within the UCL. This behaviour can be best explained when considering that the time period of the oscillation is comparable to the eddy turnover time of turbulence in the UCL, i.e. $T \approx T_h$, which is the characteristic time scale for ‘information transport’ within the UCL. At the three lower forcing frequencies, i.e. the LL, LM and LH cases, on the contrary, the interaction between the roughness elements and the unsteady flow induces an oscillation in the shear rate, which has a phase lag of roughly ${\rm \pi} /2$ with respect to the pulsatile forcing at the top of the UCL. This oscillating shear rate then propagates in the positive wall-normal direction while being progressively attenuated. The propagation speed of the oscillating shear rate appears to be constant for a given forcing frequency, which can be readily inferred by the constant tilting angle in the $\partial \langle \bar {u}_1\rangle _o/\partial x_3$ contours. The flow region affected by the oscillating shear rate defines the so-called ‘Stokes layer’. For cases with two moderate frequencies, i.e. the LM and LH cases, the Stokes layer thickness $(\delta _s)$ is smaller than the domain height $L_3$. Above the Stokes layer, the slope of $\langle \bar {u}_1\rangle _o$ is nominally zero over the pulsatile cycle, and the flow in such a region resembles a plug flow. On the contrary, in the LL case, the entire domain is affected by the oscillating shear rate, thus $\delta _s>L_3$.

Figure 18. Space–time diagrams of $\langle \bar {u}_1\rangle _o/u_\tau$ from the LL (a), LM (b), LH (c) and LVH (d) cases. Horizontal dashed lines identify the top of the UCL.

Figure 19. Space–time diagrams of $(h/u_\tau ) {\partial \langle \bar {u}_1\rangle _o}/{\partial x_3}$ from the LL (a), LM (b), LH (c) and LVH (d) cases.

Figures 20 and 21 depict the time evolution of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ and $\langle \overline {u_1^\prime u_1^\prime }\rangle _o$, respectively. Although the contours of $\langle \overline {u_2^\prime u_2^\prime }\rangle _o$ and $\langle \overline {u_3^\prime u_3^\prime }\rangle _o$ are not shown, these quantities vary in a similar fashion as $\langle \overline {u_1^\prime u_1^\prime }\rangle _o$ during the pulsatile cycle. These space–time diagrams confirm that the considered frequencies encompass three distinct flow regimes. For the LVH case, time variations of the oscillatory resolved Reynolds stresses are essentially zero above the UCL. In cases with three lower frequencies, oscillatory resolved Reynolds stresses exhibit a similar behaviour to $\partial \langle \bar {u}_1\rangle _o/\partial x_3$. Specifically, there appear oscillating waves propagating away from the UCL at a constant speed and meanwhile getting weakened. In the LM and LH cases, such oscillating waves are fully dissipated at the upper limit of the Stokes layer, above which the turbulence is ‘frozen’ and passively advected.

Figure 20. Space–time diagrams of $\langle \overline {u_1^\prime u_3^\prime }\rangle _o / u_\tau ^2$ from the LL (a), LM (b), LH (c) and LVH (d) cases.

Figure 21. Space–time diagrams of $\langle \overline {u_1^\prime u_1^\prime }\rangle _o/u_\tau ^2$ from the LL (a), LM (b), LH (c) and LVH (d) cases.

A visual comparison of the tilting angles in the contours of oscillatory resolved Reynolds stresses and $\partial \langle \bar {u}_1\rangle _o/\partial x_3$ reveals that the oscillating waves in these quantities feature similar propagation speeds. This behaviour closely resembles the one observed in smooth-wall cases (Scotti & Piomelli Reference Scotti and Piomelli2001; Manna et al. Reference Manna, Vacca and Verzicco2015). The physical interpretation is that when the oscillating shear rate is diffused upwards by the background turbulent flow, it interacts with the local turbulence via the mechanisms described in (3.19) and (3.20), thus inducing the observed oscillations in the resolved Reynolds stresses. To further quantify the propagation speeds of the oscillating waves in ${\partial \langle \bar {u}_1\rangle _o}/{ \partial x_3}$ and oscillatory resolved Reynolds stresses, figure 22 presents the phase lag of those quantities with respect to the pulsatile forcing. For a quantity $\theta$, the propagation speed $c_\theta$ is defined based on the slope of the phase lag, i.e.

(3.23)\begin{equation} c_\theta={-} \omega\frac{\partial x_3}{\partial \phi_\theta} . \end{equation}

Figure 22. Phase lag of ${\partial \langle \bar {u}_1\rangle _o}/{ \partial x_3}$ (red), $-\langle \overline {u_1^\prime u_3^\prime }\rangle _o$ (black), $\langle \overline {u_1^\prime u_1^\prime }\rangle _o$ (magenta), $\langle \overline {u_2^\prime u_2^\prime }\rangle _o$ (green) and $\langle \overline {u_3^\prime u_3^\prime }\rangle _o$ (blue) with respect to the pulsatile forcing from the LL (a), LM (b) and LH (c) cases.

Table 2 summarizes the wave propagation speeds for $\partial \langle \bar {u}_1\rangle _o /\partial x_3$, $-\langle \overline { u_1^\prime u_3^\prime }\rangle _o$, $\langle \overline { u_1^\prime u_1^\prime }\rangle _o$, $\langle \overline { u_2^\prime u_2^\prime }\rangle _o$ and $\langle \overline { u_3^\prime u_3^\prime }\rangle _o$ of each case. This again confirms that the oscillating waves propagate at a similar speed for the considered quantities. It is also noteworthy to point out that the speed of the propagating wave increases with $\omega$.

Table 2. Propagation speeds of oscillating waves in $\partial \langle \bar {u}_1\rangle _o /\partial x_3$ and oscillatory resolved Reynolds stresses.

Three other observations can be made from figure 22. First, throughout the considered cases, there appears a marked phase lag of roughly ${\rm \pi} /6$ between $\partial \langle \bar {u}_1\rangle _o /\partial x_3$ and (negative) $\langle \overline {u_1^\prime u_3^\prime }\rangle _o$, indicating a deviation from the Boussinesq eddy viscosity assumption. Weng et al. (Reference Weng, Boij and Hanifi2016) reported a similar finding, and they attributed such behaviour to non-equilibrium effects arising when the time period of the pulsatile forcing is short compared with the local turbulence relaxation time so that the turbulence is not able to relax to its equilibrium state during the pulsation cycle. Second, the lack of phase lag among the oscillatory normal resolved Reynolds stresses implies that the oscillatory pressure-redistribution terms respond immediately to the change in $\langle \overline { u_1^\prime u_1^\prime }\rangle _o$. Third, the lifetimes of oscillating waves in the resolved Reynolds stresses and shear rate, which are inferred by the difference between the phase lags at the top of the UCL and at the upper limit of the Stokes layer, are no more than half of the oscillation time period, although they decrease with $\omega$. They are considerably shorter than those in smooth-wall cases, which are typically larger than one oscillation period (Scotti & Piomelli Reference Scotti and Piomelli2001; Manna et al. Reference Manna, Vacca and Verzicco2015; Weng et al. Reference Weng, Boij and Hanifi2016).

3.3.3. Scaling of the Stokes layer thickness

Across many applications, $\delta _s$ is a quantity of interest, since it defines the region where the turbulence and the mean flow are out of equilibrium. In such a region, established turbulence theories may fail to capture flow dynamics that are of relevance for, e.g. surface drag and scalar dispersion.

For the wave-current boundary flow, where the surface is typically transitionally rough, the wave boundary layer thickness – an equivalent of Stokes layer thickness – scales as

(3.24)\begin{equation} \delta_w \sim \frac{\kappa u_{\tau,max}}{\omega} , \end{equation}

where $u_{\tau,max}$ is the friction velocity based on the maximum phase- and intrinsic-averaged wall stress during the pulsatile cycle (Grant & Madsen Reference Grant and Madsen1979). Such a scaling argument is not valid for the current cases, even though the considered surface is also rough. As shown in § 3.3.1, normalized oscillatory velocity and resolved Reynolds stresses profiles collapse between cases with the same frequencies, implying that the Stokes layer thickness is only dependent on $\omega$, whereas $\tau _{max}$ is determined by both $\alpha$ and $\omega$. Rather, the scaling of $\delta _s$ in the current cases is a trivial extension of the model first introduced by Scotti & Piomelli (Reference Scotti and Piomelli2001), as discussed next.

Let us recall from § 1 that the Stokes layer thickness for turbulent pulsatile flow over an aerodynamically smooth surface (Scotti & Piomelli Reference Scotti and Piomelli2001) is defined as

(3.25)\begin{equation} \delta_s = 2\frac{\kappa u_\tau}{\omega}\left(1+\sqrt{1+ \frac{2\nu \omega}{\kappa^2 u_\tau^2}}\right) . \end{equation}

Here we apply two modifications to this model in order to make it applicable to the current rough-wall cases. First, given that the viscous stress is omitted, the molecular viscosity $\nu$ can be neglected. Also, in the current cases, the oscillating shear rate is generated within the UCL rather than at the bottom surface (as in the smooth-wall cases), and the extent of the oscillating shear rate propagation defines the thickness of the Stokes layer. This behaviour can be easily captured by augmenting $\delta _s$ by the displacement height ($d$). Specifically, we draw an analogy to smooth-wall cases by taking $d$ as the offset, since it is the virtual origin of the longtime-averaged velocity profile. The displacement height, $d$, is a plausible choice of the offset since it captures the limiting behaviour of the flow system as the canopy packing density varies. For instance, in the limit of $\lambda _p \rightarrow 0$ (very sparse canopies), $d = 0$, i.e. the oscillating shear rate grows from the bottom of the domain. On the contrary, in the limit of $\lambda _p \rightarrow 1$ (very dense canopies), $d = h$, i.e. the oscillating shear layer starts at the top of the UCL. Based on these considerations, a phenomenological model for the Stokes layer thickness is

(3.26)\begin{equation} \delta_s=\frac{4\kappa u_\tau}{\omega}+d . \end{equation}

Note that, in the limit of $\omega \rightarrow 0$, the Stokes layer no longer exists, rendering (3.26) invalid. At this limit, as previously stated in Scotti & Piomelli (Reference Scotti and Piomelli2001), $T$ is much larger than the turbulence relaxation time. As a result, the turbulence maintains a quasiequilibrium state, and the flow statistics are indistinguishable from those of the corresponding equilibrium canopy layer flows, if scaled with the instantaneous inner/outer units.

Figure 23 compares the predictions of (3.26) against LES results. Note that only low-amplitude cases are shown, since, as mentioned earlier, $\delta _s$ only depends on $\omega$. The upper limit of the Stokes layer is identified as the location where $\sigma ^2_{ \langle \bar {k}\rangle _o}$ is $1\,\%$ of its maximum, where

(3.27)\begin{equation} \sigma^2_{ \langle\bar{k}\rangle_o}= \frac{1}{T}\int^{T}_0 \left(\frac{1}{2}\big(\langle\overline{ u_1^\prime u_1^\prime }\rangle_o+\langle\overline{ u_2^\prime u_2^\prime }\rangle_o+\langle\overline{ u_3^\prime u_3^\prime }\rangle_o\big)\right)^2 \,{\rm d}t \end{equation}

is the time variance of $\langle \bar {k}\rangle _o$. From figure 23, it is apparent that the estimated $\delta _s$ compare very well with LES results. The estimation of $\delta _s$ for the LL case is not shown in figure 23 because it exceeds the height of the computational domain. Equation (3.26) can hence be used in future studies to identify the Stokes layer thickness for pulsatile flows over aerodynamically rough surfaces.

Figure 23. Time variance of $\langle \bar {k}\rangle _o$: $\sigma ^2_{ \langle \bar {k}\rangle _o}=\int ^{T}_0 \big(\tfrac {1}{2} \big(\langle \overline { u_1^\prime u_1^\prime }\rangle _o+\langle \overline { u_2^\prime u_2^\prime }\rangle _o+\langle \overline { u_3^\prime u_3^\prime }\rangle _o\big)\big)^2 \,{\rm d}t/T$ of the LL (navy blue), LM (dark green), LH (light green) and LVH (yellow green) cases. Circle denotes $\delta _s$ estimated via (3.26). Triangle represents the location where $\sigma ^2_{ \langle \bar {k}\rangle _o}$ is reduced to $1\,\%$ of its maximum value.