Mean flow and turbulence in unsteady canopy layers

Abstract Non-stationarity is the rule in the atmospheric boundary layer (ABL). Under such conditions, the flow may experience departures from equilibrium with the underlying surface stress, misalignment of shear stresses and strain rates, and three-dimensionality in turbulence statistics. Existing ABL flow theories are primarily established for statistically stationary flow conditions and cannot predict such behaviours. Motivated by this knowledge gap, this study analyses the impact of time-varying pressure gradients on mean flow and turbulence over urban-like surfaces. A series of large-eddy simulations of pulsatile flow over cuboid arrays is performed, programmatically varying the oscillation amplitude $\alpha$ and forcing frequency $\omega$. The analysis focuses on both longtime-averaged and phase-dependent flow dynamics. Inspection of longtime-averaged velocity profiles reveals that the aerodynamic roughness length $z_0$ increases with $\alpha$ and $\omega$, whereas the displacement height $d$ appears to be insensitive to these parameters. In terms of oscillatory flow statistics, it is found that $\alpha$ primarily controls the oscillation amplitude of the streamwise velocity and Reynolds stresses, but has a negligible impact on their wall-normal structure. On the other hand, $\omega$ determines the size of the region affected by the unsteady forcing, which identifies the so-called Stokes layer thickness $\delta _s$. Within the Stokes layer, phase-averaged resolved Reynolds stress profiles feature substantial variations during the pulsatile cycle, and the turbulence is out of equilibrium with the mean flow. Two phenomenological models have been proposed that capture the influence of flow unsteadiness on $z_0$ and $\delta _s$, respectively.

Major drivers of non-stationarity in the ABL include time-varying horizontal pressure gradients, associated with non-turbulent motions ranging from submeso to synoptic scales, and time-dependent thermal forcings, induced by the diurnal cycle or by cloudinduced time variations of the incoming solar radiation (Mahrt & Bou-Zeid 2020).These conditions often result in departures from equilibrium turbulence, with important implications on time-and area-averaged exchange processes between the land surface and the atmosphere.The first kind of non-stationarity was examined in Mahrt (2007Mahrt ( , 2008)); Mahrt et al. (2013), which showed that time-variations of the driving pressure gradient can enhance momentum transport under strong stable atmospheric stratifications.The second kind of non-stationarity was instead analyzed in Hicks et al. (2018) making use of data from different field campaigns, and showed that the surface heat flux can change so rapidly during the morning and late afternoon transition that the relations for equilibrium turbulence no longer hold.
Numerical studies have also been recently conducted to study how exchange processes between the land surface and the atmosphere are modulated by non-stationarity in the ABL.In their study, Edwards Edwards et al. (2006) conducted a comparison of a prevailing single-column model based on equilibrium turbulence theories with the observations of an evening transition ABL, as well as results from large-eddy simulations (LES).Their findings emphasized the inadequacy of equilibrium turbulence theories in capturing the complex behavior of ABL flows during rapid changes in thermal surface forcing.This breakdown of equilibrium turbulence was particularly notable during the evening transition period, which is known for its rapid changes in thermal forcing.Momen & Bou-Zeid (2017) investigated the response of the Ekman boundary layer to oscillating pressure gradients, and found that quasi-equilibrium turbulence is maintained only when the oscillation period is much larger than the characteristic time scale of the turbulence.The majority of the efforts have focused on atmospheric boundary layer flow over modeled roughness, where the flow dynamics in the roughness sublayer-that layer of the atmosphere that extends from the ground surface up to about 2-5 times the mean height of roughness elements (Fernando 2010)-are bypassed, and surface drag is usually evaluated via an equilibrium wall-layer model (see, e.g., Momen & Bou-Zeid (2017)), irrespective of the equilibrium-theory limitations outlined above.It hence remains unclear how unsteadiness impacts flow statistics and the structure of atmospheric turbulence in the roughness sublayer.Roughness sublayer flow directly controls exchanges of mass, energy, and momentum between the land surface and the atmosphere, and understanding the dependence of flow statistics and structural changes in the turbulence topology on flow unsteadiness is hence important in order to advance our ability to understand and predictively model these flow processes.
This study contributes to addressing this knowledge gap by focusing on nonstationarity roughness sublayer flow induced by time-varying pressure gradients.Unsteady pressure gradients in the real-world ABL can be characterized by periodic and aperiodic variations in both magnitude and direction.In this study, we limit our attention to a pulsatile streamwise pressure-gradient forcing, consisting of a constant mean and a sinusoidal oscillating component.This approach has two major merits.First, the temporal evolution of flow dynamics and associated statistics, as well as structural changes in turbulence, can be easily characterized thanks to the time-periodic nature of the flow unsteadiness.Second, the time scale of the pulsatile forcing is well defined, and can hence be varied programmatically to encompass a range of representative flow regimes.
Pulsatile turbulent flows over aerodynamically smooth surfaces have been the subject of active research in the mechanical engineering community because of their relevance across a range of applications; these include industrial (e.g., a rotating or poppet valve) and biological (blood in arteries) flows.The corresponding laminar solution is an extension of Stokes's second problem (Stokes 1901), where the modulation of the flow field by unsteady pressure gradient is confined to a layer of finite thickness known as the "Stokes layer".The thickness of the Stokes layer δ s is a function of the pulsatile forcing frequency ω, i.e., δ s = 2l s , where l s = 2ν/ω is the so-called Stokes length scale and ν is the kinematic viscosity of the fluid.In the turbulent flow regime, it has been found that the characteristics of the pulsatile flow are not only dependent on the forcing frequency, but also on the amplitude of the oscillation.Substantial efforts have been devoted to investigating this problem, from both an experimental (Ramaprian & Tu 1980;Tu & Ramaprian 1983;Ramaprian & Tu 1983;Mao & Hanratty 1986;Brereton et al. 1990;Tardu & Binder 1993;Tardu 2005) and a computational perspective (Scotti & Piomelli 2001;Manna et al. 2012Manna et al. , 2015;;Weng et al. 2016).Scotti & Piomelli (2001) drew an analogy to the Stokes length and proposed a turbulent Stokes length scale, which can be expressed in inner units as where u τ is the friction velocity based on the surface friction averaged over pulsatile cycles and ν t is the so-called eddy viscosity.Parameterizing the eddy viscosity in terms of the Stokes turbulent length scale, i.e., ν t = κu τ l t , where κ is the von Kármán constant, and substituting into (1.1)one obtains When l + t is large (e.g., when ω → 0), the flow is in a quasi-steady state.Under such a condition, the flow at each pulsatile phase resembles a statistically stationary boundary layer flow, provided that the instantaneous friction velocity is used to normalize statistics.As ω increases and l + t becomes of the order of the open channel height L + 3 , the entire flow is affected by the pulsation, i.e., time lags occur between flow statistics at different elevations, and turbulence undergoes substantial structural changes from its equilibrium configuration.When l + t < L + 3 /2, the flow modulation induced by the pulsation is confined within the Stokes layer δ + s = 2l + t .Above the Stokes layer, one can observe a plug-flow region, with the turbulence being frozen to its equilibrium configuration and simply advected by the mean flow pulsation.A few years later, Bhaganagar (2008) conducted a series of direct numerical simulations (DNS) of low-Reynolds-number pulsatile flow over transitionally rough surfaces.She found that flow responses to pulsatile forcing are generally similar to those in smooth-wall cases, when the roughness size is of the same order of magnitude as the viscous sublayer thickness.The only exception is that, as the pulsation frequency approaches the frequency of vortex shedding from the roughness elements, the longtime averaged velocity profile deviates significantly from that of the steady flow case due to the resonance between the pulsation and the vortex shedding.In the context of a similar flow system, Patil & Fringer (2022) also reported a comparable observation in their DNS study.
In addition to the work of Bhaganagar (2008); Patil & Fringer (2022), pulsatile flow over small-scale roughness, e.g., sand grain roughness, has also been studied extensively in the oceanic context, i.e., combined current-wave boundary layers, which play a crucial role in controlling sediment transport and associated erosion in coastal environments (Grant & Madsen 1979;Kemp & Simons 1982;Myrhaug & Slaattelid 1989;Sleath 1991;Soulsby et al. 1993;Mathisen & Madsen 1996;Fredsøe et al. 1999;Yang et al. 2006;Yuan & Madsen 2015).The thickness of the wave boundary layer, which is the equivalent of the Stokes layer in the engineering community, is defined as where u τ,max = √ τ max , and τ max denotes the maximum of kinematic shear stress at the surface during the pulsatile cycle.Within the wave boundary layer, mean flow and turbulence are controlled by the nonlinear interaction between currents and waves.Above this region, the modulation of turbulence by waves vanishes.The wall-normal distribution of the averaged velocity over the pulsatile cycle deviates from the classic logarithmic profile, and is characterized by a "two-log" profile, i.e., the velocity exhibits a logarithmic profile with the actual roughness length within the wave boundary layer and a different one characterized by a larger roughness length further aloft (Grant & Madsen 1986;Fredsøe et al. 1999;Yang et al. 2006;Yuan & Madsen 2015).Such behavior was first predicted by a two-layer time-invariant eddy viscosity model by Grant & Madsen (1979), followed by many variants and improvements (Myrhaug & Slaattelid 1989;Sleath 1991;Yuan & Madsen 2015).
On the contrary, pulsatile flow at high Reynolds numbers over large roughness elements, such as buildings, has received far less attention.Yu et al. (2022) conducted a series of LES of combined wave-current flows over arrays of hemispheres, which can be seen as a surrogate of reefs near the coastal ocean.They focused only on low Keulegan-Carpenter numbers KC ∼ O(1 − 10), where KC is defined as the ratio between the wave excursion U w T and the diameter of the hemispheres, and U w and T are the wave orbital velocity and the wave period, respectively (Keulegan et al. 1958).However, conclusions from Yu et al. (2022) cannot be readily applied to pulsatile flow over urban-like roughness (i.e., large obstacles with sharp edges), mainly because different surface morphologies yield distinct air-canopy interaction regimes under pulsatile forcings (Carr & Plesniak 2017).
Motivated by this knowledge gap, this study proposes a detailed analysis on the dynamics of the mean flow and turbulence in high-Reynolds-number pulsatile flow over idealized urban canopies.The analysis is carried out based on a series of LES of pulsatile flow past an array of surface-mounted cuboids, where the frequency and amplitude of the pressure gradient are programmatically varied.The LES technique has been shown capable of capturing the major flow features of pulsatile flow over various surface conditions (Scotti & Piomelli 2001;Chang & Scotti 2004).The objective of this study is to answer fundamental questions pertaining to the impacts of the considered flow unsteadiness on the mean flow and turbulence in the urban boundary layer: (i) Does the presence of flow unsteadiness alter the mean flow profile in a longtimeaveraged sense?If so, how do such modifications reflect in the aerodynamic surface parameters?
(ii) To what extent does the unsteady pressure gradient impact the overall momentum transport and turbulence generation in a longtime-averaged sense within and above the canopy?(iii) How do the phase-averaged mean flow and turbulence behave in response to the periodically varying pressure gradient?How are such phase-dependent behaviors controlled by the oscillation amplitude and the forcing frequency?
This paper is organized as follows.Section 2 introduces the numerical algorithm and the setup of simulations, along with the flow decomposition and averaging procedure.Results are presented and discussed in §3.Concluding remarks are given in §4.

Numerical procedure
A suite of LES is performed using an extensively validated in-house code (Albertson & Parlange 1999a,b;Bou-Zeid et al. 2005;Chamecki et al. 2009;Anderson et al. 2015;Fang & Porté-Agel 2015;Li et al. 2016;Giometto et al. 2016).The code solves the filtered continuity and momentum transport equations in a Cartesian reference system, which read ∂u i ∂x i = 0 , (2.1) where u 1 , u 2 , and u 3 are the filtered velocities along the streamwise (x 1 ), lateral (x 2 ), and wall-normal (x 3 ) direction, respectively.The advection term is written in the rotational form to ensure kinetic energy conservation in the discrete sense (Orszag & Pao 1975).ρ represents the constant fluid density, τ ij is the deviatoric component of the subgridscale (SGS) stress tensor, which is evaluated via the Lagrangian scale-dependent dynamic (LASD) Smagorinsky model (Bou-Zeid et al. 2005).The LASD model has been extensively validated in wall-modeled simulations of unsteady atmospheric boundary layer flow (Momen & Bou-Zeid 2017;Salesky et al. 2017) and in the simulation of flow over surface-resolved urban-like canopies (Anderson et al. 2015;Li et al. 2016;Giometto et al. 2016;Yang 2016).Note that viscous stresses are neglected in the current study; this assumption is valid as the typical Reynolds number of the ABL flows is Re ∼ O(10 9 ), and the flow is in the fully rough regime.p * = p + 1 3 ρτ ii + 1 2 ρu i u i is the modified pressure, which accounts for the trace of SGS stress and resolved turbulent kinetic energy.The flow is driven by a spatially uniform but temporally periodic pressure gradient, i.e., (2.3) where f m denotes the mean pressure gradient.α p is a constant controlling the amplitude of the forcing, and ω represents the forcing frequency.δ ij is the Kronecker delta tensor.Periodic boundary conditions apply in the wall-parallel directions, and free-slip boundary conditions are employed at the upper boundary.The lower surface is representative of an array of uniformly distributed cuboids, which serves as a surrogate of urban landscapes.This approach has been commonly adopted in studies of ABL, as this type of surface morphology is characterized by a limited number of length scales, making it amenable to analytical treatment (Reynolds & Castro 2008;Cheng & Porté-Agel 2015;Tomas et al. 2016;Basley et al. 2019;Omidvar et al. 2020).Such an approach is justified on the basis that one should first study a problem in its simplest setup before introducing additional complexities.Nonetheless, it is crucial to acknowledge that the introduction of randomness in roughness alters the flow characteristics and the generation of turbulence in rough-wall boundary layer flows.Xie et al. (2008) studied flows over random urban-like obstacles and found that turbulence features in the roughness sublayer are controlled by the randomness in the roughness.Giometto et al. (2016) conducted an LES study and highlighted that roughness randomness enhances the dispersive stress in the roughness sublayer.Chau & Bhaganagar (2012) carried out a series of DNS of flow over transitionally rough surfaces, demonstrating that different levels of roughness randomness lead to distinct turbulence structures in the near-wall region and subsequently affect turbulence intensities.
Spatial derivatives in the wall-parallel directions are computed via a pseudo-spectral collocation method based on truncated Fourier expansions (Orszag 1970), whereas a second-order staggered finite difference scheme is employed in the wall-normal direction.A second-order Adams-Bashforth scheme is adopted for time integration.Nonlinear advection terms are de-aliased via the 3/2 rule (Canuto et al. 2007;Margairaz et al. 2018).Roughness elements are explicitly resolved via a discrete-forcing immersed boundary method (IBM), which is also referred to as the direct forcing IBM in the engineering fluid mechanics community (Mohd-Yusof 1996;Mittal & Iaccarino 2005;Fang et al. 2011).The IBM was originally developed in Mohd-Yusof (1996) and first introduced to ABL studies by Chester et al. (2007).Since then, the IBM has been extensively validated in subsequent studies (e.g.Graham & Meneveau 2012;Cheng & Porté-Agel 2015;Giometto et al. 2016;Anderson 2016;Yang & Anderson 2018;Li & Bou-Zeid 2019).Specifically, an artificial force F i drives the velocity to zero within the cuboids, and an inviscid equilibrium logarithmic wall-layer model (Moeng 1984;Giometto et al. 2016) is applied over a narrow band centered at the fluid-solid interface to evaluate the wall stresses.As shown in Appendix A, for flow over cuboids in the fully rough regime, the use of an equilibrium wall model does not impact the flow field significantly.Xie & Castro (2006) reached a similar conclusion for a comparable flow system.The incompressibility condition is then enforced via a pressure-projection approach Kim & Moin (1985).
Figure 1 shows a schematic of the computational domain.The size of the domain is , and L 3 = 8h, where h is the height of the cuboids.The planar and frontal areas of the cube array are set to λ p = λ f = 0.1.An aerodynamic roughness length of z 0 = 10 −4 h is prescribed at the cube surfaces and the lower surface via the wall-layer model.With the chosen value of z 0 , the SGS pressure drag is a negligible contributor to the overall momentum balance (Yang & Meneveau 2016).The domain is discretized using a uniform Cartesian grid (N 1 , N 2 , N 3 ) = (576,192,128) where each cube is resolved by (n 1 , n 2 , n 3 ) = (8, 8, 16) grid points.As shown in Appendix B, this resolution yields flow statistics-up to second-order moments-that are poorly sensitive to grid resolution.

Averaging operations
Given the time-periodic nature of the flow system, phase averaging is the natural approach to evaluate flow statistics in pulsatile flows (Scotti & Piomelli 2001;Bhaganagar 2008;Weng et al. 2016;Önder & Yuan 2019).Phase averaging can be best understood as a surrogate of Reynolds ensemble averaging for time-periodic flows.The phase and intrinsic volume average (Schmid et al. 2019) (hereafter referred to as phase-average for brevity) of a quantity of interest θ can be defined as (2.4) where • denotes the phase averaging operation, V f is a thin fluid slab of thickness δ 3 in the x 3 direction, N p denotes the number of the pulsatile cycles over which the averaging operation is performed, and T = 2π/ω is the time period of the pulsatile forcing.A given instantaneous quantity θ can be decomposed as where (•) denotes a departure of the instantaneous value from the corresponding phaseaveraged quantity.A phase-averaged quantity can be further decomposed into a longtime average and an oscillatory component with zero mean, i.e., This work relies on the Scotti & Piomelli (2001) approach to analyze the flow system; in this approach, an oscillatory quantity θ is split into two components: one corresponding to the flow oscillation at the forcing frequency (fundamental mode), and one which includes contributions from all of the remaining harmonics, i.e., where A θ and φ θ are the oscillatory amplitude of the fundamental mode and the phase lag with respect to the pulsatile forcing, respectively.These components are evaluated via minimization of e θ 2 at each x 3 .

Dimensional analysis and suite of simulations
Having fixed the domain size, the aerodynamic surface roughness length, and the spatial discretization, the remaining physical parameters governing the problem are: (i) the oscillation amplitude α, defined as the ratio between the oscillation amplitude of u 1 at the top of the domain and the corresponding mean value; (ii) the forcing frequency ω; (iii) the friction velocity based on the mean pressure gradient u τ = √ f m L 3 ; and (iv) the height of roughness elements h.The latter is a characteristic length scale of the flow in the urban canopy layer (UCL).System response is studied in time and along the x 3 coordinate direction, so (v) the wall-normal elevation x 3 and (vi) time t should also also be included in the parameter set.Note that the viscosity is not taken into account since the flow is in the fully rough regime.Choosing u τ and h as repeating parameters, a given normalized longtime (Y ) and phase-average ( Y ) quantity of interest can hence be written as respectively, where f and g are universal functions.(2.8) show that Y and Y only depend on two dimensionless parameters, namely α and ωT h .T h = h/u τ is the turnover time of the largest eddies in the UCL and can be best understood as a characteristic time scale of the flow in the UCL.ωT h ∼ T h /T is hence essentially the ratio between the turnover time of the largest eddies in the UCL (T h ) and the pulsation time period (T ).Also, note that given their identical mathematical formulation, the normalized ω can be seen as an equivalent Strouhal number which has been used as a non-dimensional parameter to characterize pulsatile flows in early works on pulsatile flow.Also, note that ωT h can be seen as an equivalent Strouhal number, which has been used as a nondimensional parameter to characterize pulsatile flows in earlier studies of this flow system, such as Tardu et al. (1994), given their identical mathematical formulation.Four different values of ωT h are considered in this study, namely wT h = {0.05π,0.125π, 0.25π, 1.25π}.For u τ ≈ 0.1 ms −1 and h ≈ 10 m-common ABL values for these quantities (Stull 1988)-the considered ωT h set encompasses time scales variability from a few seconds to several hours, which are representative of submeso-scale phenomena (Mahrt 2009;Hoover et al. 2015).
In this study, we aim to narrow our focus to the current-dominated regime, i.e., 0 < α < 1. Larger values of α would lead to a wave-dominated regime, which behaves differently.The specific values (α = {0.2,0.4}) have been chosen because they are sufficiently large to lead to an interesting flow response, yet sufficiently far from the wave-dominated regime; this enables us to focus on departures from stationarity induced by flow pulsation in the current-dominated regime.Due to the lack of a straightforward relation between the oscillatory pressure gradient (α p ) and α, α p has been tuned iteratively to achieve the desired α.As shown in table 1, despite the optimization process, there are still discernible discrepancies between the target α and its actual value.As also shown in previous works (Scotti & Piomelli 2001;Bhaganagar 2008), α is highly sensitive to variations in α p , which makes it challenging to obtain the desired α values.
The suite of simulations and corresponding acronyms used in this study are listed in table 1.A statistically stationary flow case (α = 0) is also carried out to highlight departures of pulsatile flow cases from the steady state condition.Simulations with pulsatile forcing are initialized with velocity fields from the stationary flow case; the ωT h = {0.25π,1.25π} and ωT h = {0.05π,0.125π} cases are then integrated in time over 200T L3 and 400T L3 , respectively, where T L3 = L 3 /u τ is the turnover time of the largest eddies in the domain.This approach yields converged phase-averaged flow statistics.The size of the time step δt is chosen to satisfy the Courant-Friedrichs-Lewy stability condition (uδt)/δ x 0.05, where u is the maximum velocity magnitude at any given spatial location and time during a run, and δ x is the grid stencil in the computational domain.Instantaneous three-dimensional snapshots of the velocity and pressure fields are collected every T /16 for the ωT h = {0.25π,1.25π} cases and every T /80 for the ωT h = {0.05π,0.125π} cases, after an initial 20T L3 transient period.

Instantaneous velocity field
To gain insights into the instantaneous flow field, figure 2 displays the streamwise fluctuating velocity within and above the UCL from the HM case at two different x 3 planes and phases.The chosen two phases t = 0 and t = π/2 correspond to the end of the deceleration and acceleration period, respectively.
It is apparent from figure 2(a, b) that meandering low-momentum (u < 0) streaks are flanked by adjacent high-momentum (u > 0) ones.Comparing these two, turbulence structures at t = π/2 appear smaller in size in both streamwise and spanwise directions.Additionally, within the UCL, apparent vortex shedding occurs on the lee side of the cubes at t = π/2, while it is less pronounced at t = 0.This suggests that flow unsteadiness substantially modifies the flow field during the pulsatile cycle and is expected to impact the flow statistics substantially.3.2.Longtime-averaged statistics

Longtime-averaged velocity profile
Profiles of the longtime-averaged streamwise velocity are shown in figure 3. Flow unsteadiness leads to a horizontal shift of profiles in the proposed semi-logarithmic plot.This behavior is distinct from the "two-log" profile in flow over sand grain roughness (Fredsøe et al. 1999;Yang et al. 2006;Yuan & Madsen 2015), 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 1993;Tardu 2005;Scotti & Piomelli 2001;Manna et al. 2012;Weng et al. 2016).As also apparent from figure 3, departures from the statistically stationary flow profile become more significant for larger values of α and ω.
Variations in the aerodynamic roughness length z 0 and displacement height d parameters with ωT h are shown in figure 4 for the considered canopy.These parameters are evaluated via the Macdonald et al. (1998) approach, where d is the barycenter height of the longtime-averaged pressure drag from the urban canopy, and z 0 is determined via curve fitting.More specifically, 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., Note that, in principle, one should also account for the SGS drag contribution in (3.2); in this work, we omit SGS contributions because they are negligible when compared to the total drag-a direct result of the relatively small aerodynamic roughness length that is prescribed in the wall-layer model (see §2.1).z 0 is solved by minimizing the mean square error (MSE) between the longtime-averaged velocity and the law of the wall with κ = 0.4 in the x 3 ∈ [2h, 6h] interval, i.e., The fitting interval is highlighted in figure 3. 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. (2007) 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 4, d is poorly sensitive to variations in both α and ω (variations across cases are within the ±3% range).This behavior can be explained by considering that for flows over sharp-edged obstacles, such as the ones considered herein, flow separation patterns are poorly sensitive to variations in α and ω, resulting in a rather constant total volume of wake regions and momentum deficits across cases, and constant longtime-averaged pressure on the surfaces of the cuboids.The d parameter is here evaluated as integral of the pressure gradient field over the surface area, so the above considerations provide a physical justification for the observed behavior.Note that this finding might not be generalizable across all possible roughness morphologies.For instance, Yu et al. (2022) showed that separation patterns in pulsatile flows over hemispheres feature a rather strong dependence on α and ω, yielding corresponding strong variations in d.
Contrary to d, the z 0 parameter is strongly impacted by flow unsteadiness, and its value increases with α and ω.Bhaganagar (2008) 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.

Phase-averaged drag-velocity relation
The phase-averaged canopy drag D and the local phase-averaged velocity u 1 at x 3 /h ≈ 0.8 are shown in figure 5, 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 frequencyindependent one-to-one mapping between D and u 1 .As apparent from figure 6(a), at these three forcing frequencies, the interaction between the wind and the canopy layer is in a state of quasi-equilibrium, i.e., D is in phase with u 1 .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 1969;Raupach 1992;Katul et al. 2004;Poggi et al. 2004;Macdonald et al. 1998;Coceal & Belcher 2004).This finding comes as no surprise, given the quasi-equilibrium state of the D − u 1 relation for the three lower forcing frequencies.
On the other hand, as shown in figure 6(b), results from the highest-frequency cases (LVH and HVH) exhibit an orbital pattern, which stems from the time lag between D and u 1 .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 5b).
These findings suggest the following parameterization for D : where C d is a sectional drag coefficient that is constant in time and does not depend on α nor ω, and ∆t accounts for the time lag between D and u 1 , which instead does indeed depend on α and ω.Note that, throughout the considered cases, the wall-normalaveraged drag coefficient h 0 C d dx 3 /h ≈ 0.9-a value that is similar to those previously reported ones for stationary flow over cube arrays (Coceal & Belcher 2004) (note that the exact value depends on the formula used to define C d ).  1. Morison et al. (1950) developed a semi-empirical model relating the phase-averaged drag generated by obstacles in an oscillatory boundary layer to a given phase-averaged velocity-a model that has been extensively used in the ocean engineering community to evaluate drag from surface-mounted obstacles (Lowe et al. 2005;Yu et al. 2018Yu et al. , 2022)).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.One might argue that the Morison model could also be used to evaluate surface drag as a function or the phase-averaged velocity at a given x 3 for the cases under consideration, but unfortunately, this is not the case.As shown in Patel & Witz (2013), 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 important phase lags between phase-averaged accelerations at different wallnormal locations substantially degrade the accuracy of such a model.This behavior can be easily inferred from figure 16.
In the following, we will make use of (3.4) to derive an alternative phenomenological surface-drag model for the considered flow system.

Mapping roughness length variability to longtime-averaged flow statistics
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. 2008;Benjamin et al. 2016).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 where D 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 u 1 .Thus, it is reasonable to assume C d is constant along the wall-normal direction.Also, depending on α and ω, the flow within the canopy might undergo a local reversal in the phase-averaged sense, meaning that u 1 < 0 at selected x 3 locations.Assuming that there is no flow reversal within the UCL, i.e., u 1 0 in z h, (3.5) can be written as where the second term on the right-hand side of (3.6) is identically zero for the SS case.
(3.6) 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.6)).Note that h 0 (•) 2 dx 3 /h is an averaging operation over the UCL based on the L 2 norm.Rearranging terms in (3.6) leads to For the pulsatile cases and the SS case, respectively.Here (•) avg denotes the canopyaveraged quantity, and (•) SS represents a quantity pertaining to the SS case.
As discussed in §3.2.1, flow unsteadiness yields a shift of the u 1 profile with negligible variations in the d parameters when compared to the stationary flow with the same u τ .In terms of the law-of-the-wall, this behavior can be described as a variation in z 0 , i.e., where (3.9) is valid for any x 3 in the logarithmic region.The shift in the velocity profile is approximately constant for x 3 ∈ [0, L 3 ], so one can write and substituting (3.7)-(3.9)into (3.10)finally yields ) is a diagnostic model relating variations in the z 0 parameter to the UCL phaseaveraged velocity variance-a longtime-averaged quantity.z 0 estimates from (3.11) are compared against LES results in figure 7, 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 HHVW 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.(3.11) highlights that, in the absence of flow reversal, z 0 can be described as a monotonically increasing function of the u 1 variance in the UCL.As explained at the beginning of this section, this finding is important from a flow modeling perspective, because it relates a longtime-averaged flow statistic to the z 0 parameter.Note that z SS 0 can be accurately evaluated using any existing parameterization for stationary ABL flow over aerodynamically rough surfaces, including the Lettau (1969), Raupach (1992), and Macdonald et al. (1998) models.Further, C d = 0.9 (see discussion in §3.2.2) and λ f and h are morphological parameters that are in general a-priori available.In weather forecasting and climate models, u 1 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.

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, u 1 u 3 profiles are indistinguishable from the SS one (not shown), indicating a weak dependence of such a quantity on α and ω.Above the UCL, the divergence of u 1 u 3 balances the longtime-averaged driving pressure gradient f m , i.e., Since f m does not vary across the considered cases and u 1 u 3 (L 3 ) = 0, a collapse of u 1 u 3 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 Given that u 1 u 3 profiles collapse in this region, D is also expected to feature a weak dependence on α and ω, which in turn explains the weak variations in the d parameter that were observed in §3.2.1.
Figure 8 depicts wall-normal profiles of the longtime-averaged resolved turbulent kinetic energy, which is defined as Such a quantity features a relatively rapid increase in the UCL, which as discussed in Schmid et al. (Schmid et al. 2019), is due to dispersive contributions caused by the canopy geometry.wh Also note that the rapid decrease in the near-surface region signals the presence of variations over scales of variability smaller than the vertical grid stencil (Coceal et al. 2007).Flow unsteadiness has a relatively more important impact on such a quantity, especially for flow in the UCL and for cases with a high oscillation amplitude.An increase in α generally results in more pronounced departures of k profiles from the SS one.This behavior can be best explained by considering the shear production terms in the budget equation for k, i.e., .15)where P 1 represents the work done by u 1 u 3 onto the longtime-averaged flow field, and P 2 is the longtime average of the work done by the oscillatory shear stress u 1 u 3 onto the Figure 10.Longtime-averaged resolved normal Reynolds stresses u 1 u 1 , u 2 u 2 , and u 3 u 3 from low-amplitude cases (a, b, c) and high-amplitude cases (d, e, f ).Line colors correspond to those used in figure 3.
oscillatory flow field.Figure 9(a) shows that P 1 is poorly sensitive to variations in α and ω.This behavior stems from the constancy of u 1 u 3 and ∂u 1 /∂x 3 across cases (the latter can be inferred from the systematic shift of u 1 profiles in figure 3).Conversely, P 2 from high-amplitude cases are generally larger than those from low-amplitude ones, mainly due to the higher u 1 u 3 and u 1 values.Discrepancies in P 2 among high-amplitude cases are larger than those among low-amplitude ones, which ultimately yields the observed variability in k.
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 10.In this case, it is apparent that increases in the oscillation frequency lead to a decrease in u 1 u 1 and an increase in u 2 u 2 and u 3 u 3 within the UCL-a behavior that is especially apparent for the high-amplitude cases (figure 10(e,d,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., R ii are responsible for redistributing kinetic energy among the longtime-averaged normal Reynolds stresses (Pope 2000), and are shown in figure 11.With the exception of the very near surface region (x 3 0.2) where no clear trend can be observed, increases in ω and α yield a decrease in R 11 and increase in R 22 and R 33 , which justify the observed isotropization of turbulence in the UCL.

Oscillatory fields
This section shifts the focus to the time evolution of velocity and resolved Reynolds stresses during the pulsatile cycle.

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 (α) 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 behaviors for the other low and high-amplitude cases, respectively.
Figure 12 contrasts the profiles of oscillatory velocity ( u 1 ) from the LL and HL runs Profiles of the oscillatory resolved turbulent kinetic energy k = ( u 1 u 1 + u 2 u 2 + u 3 u 3 )/2 from the LL (blue) and HL (red) cases at different phases of the pulsatile cycle.Profiles are T /4 apart.
at four equispaced phases during the pulsatile cycle.0 < t < T /2 and T /2 < t < T are flow acceleration and deceleration periods, respectively.u 1 at the top of the domain is controlled by the α parameter, so it is natural to use u τ α as a normalization constant to study the problem.Manna et al. (2012) 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 u 1 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 u 1 /(u τ α) profiles in figure 12.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.The oscillatory resolved turbulent kinetic energy can be defined as As shown in figure 13 and 14, both the oscillatory resolved Reynolds shear stress u 1 u 3 and k scale with u 2 τ α.Although not shown, the three oscillatory normal Reynolds stresses also obey such a scaling, suggesting that any change in k is proportionally distributed to the three normal Reynolds stresses during the pulsatile cycle.As discussed in the following paragraphs, the scaling of u 1 u 3 and k is a direct consequence of the mild non-linearity in the production of k and u 1 u 3 .Subtracting the budget equation of u 1 u 3 from that of u 1 u 3 , one obtains, where P 13,l1 and P 13,l2 are the linear production terms of u 1 u 3 , while P 13,nl is the nonlinear production.Similarly, the budget equations of k can be written as where P k,l1 and P k,l2 are the linear production terms of k, while P k,nl is the nonlinear production.As shown in figure 15, the nonlinear production term is substantially smaller than the sum of the corresponding linear productions for both u 1 u 3 and k.Given that u 1 u 3 , u 3 u 3 , and ∂u 1 /∂x 3 from the LL and HL cases are similar, when u 1 ∼ u τ α and u 3 u 3 ∼ u 2 τ α, the total production of u 1 u 3 is whereas that of k is This in turn leads to the observed u 1 u 3 ∼ u 2 τ α and k ∼ u 2 τ α.This scaling is expected to fail under two conditions.First, when α is sufficiently large, the contribution of P k,nl and P k,nl can no longer be neglected.Second, when the variations in u 1 u 3 or u 3 u 3 among cases with different α are so large that the linear production terms do not scale with u 2 τ α any more, which is what occurs within the UCL for the LVH and HVH cases (see figure 10), (3.20) and (3.21) do not hold.
The above analysis has shown that the α parameter primarily controls the amplitude of the oscillatory flow quantities, but has little impact on the wall-normal profiles of those quantities, which are instead controlled by ω, as will be shown in the following section.This behavior is also expected to hold in the smooth-wall setup (Manna et al. 2012), although the mechanism responsible for generating drag over aerodynamically-smooth surfaces is quite distinct.

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.
Figure 16 and 17 present the time evolution of u and ∂ u/∂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 ∂ u/∂x 3 is typically confined within the UCL.This behavior 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 ≈ 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 π/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 ∂ u 1 /∂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 (δ s ) is smaller than the domain height L 3 .Above the Stokes layer, the slope of u 1 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 δ s > L 3 .Figure 18 and 19 depict the time evolution of u 1 u 3 and u 1 u 1 , respectively.Although the contours of u 2 u 2 and u 3 u 3 are not shown, these quantites vary in a similar fashion as u 1 u 1 during the puslatile 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 behavior to ∂ u 1 /∂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.
A visual comparison of the tilting angles in the contours of oscillatory resolved Reynolds stresses and ∂ u 1 /∂x 3 reveals that the oscillating waves in these quantities feature similar propagation speeds.This behavior closely resembles the one observed in smooth-wall cases (Scotti & Piomelli 2001;Manna et al. 2015).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.18) and (3.19), thus inducing the observed oscillations in the resolved Reynolds stresses.To further quantify the propagation speeds of the oscillating waves in ∂ u 1 /∂x 3 and oscillatory resolved Reynolds stresses, figure 20 present the phase lag of those quantities with respect to the pulsatile forcing.For a quantity θ, the propagation speed c θ is defined based on the slope of the phase lag, i.e., Table 2 summarizes the wave propagation speeds for ∂ u 1 /∂x 3 , − u 1 u 3 , u 1 u 1 , u 2 u 2 , and u 3 u 3 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 ω.
Three other observations can be made from figure 20.First, throughout the considered cases, there appears a marked phase lag of roughly π/6 between ∂ u 1 /∂x 3 and (negative) u 1 u 3 , indicating a deviation from the Boussinesq eddy viscosity assumption.Weng et al. (2016) reported a similar finding, and they attributed such behavior to non-equilibrium effects arising when the time period of the pulsatile forcing is short compared to 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 u 1 u 1 .Third, the lifetimes of oscillating waves in the resolved Reynolds stresses and shear rate, which is 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 ω.They are considerably shorter than those in smooth-wall cases, which are typically larger than one oscillation period (Scotti & Piomelli 2001;Manna et al. 2015;Weng et al. 2016).

Scaling of the Stokes layer thickness
δ s is a quantity of interest across many applications, 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 where u τ,max is the friction velocity based on the maximum phase-averaged wall stress during the pulsatile cycle (Grant & Madsen 1979).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 ω, whereas τ max is determined by both α and ω.Rather, the scaling of δ s in the current cases is a trivial extension of the model first introduced by Scotti & Piomelli (2001), as discussed next.
Let us recall from §1 that the Stokes layer thickness for turbulent pulsatile flow over aerodynamically-smooth surface (Scotti & Piomelli 2001) is defined as ) . (3.24) 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 ν 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 behavior can be easily captured by augmenting δ 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.d is a plausible choice of the offset since it captures the limiting behavior of the flow system as the canopy packing density varies.For instance, in the limit of λ p → 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 λ p → 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 Note that, in the limit of ω → 0, the Stokes layer no longer exists, rendering (3.25) invalid.At this limit, as previously stated in Scotti & Piomelli (2001), T is much larger than the turbulence relaxation time.As a result, the turbulence maintains a quasi-equilibrium 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 21 compares the predictions of (3.25) against LES results.Note that only lowamplitude cases are shown, since, as mentioned earlier, δ s only depends on ω.The upper limit of the Stokes layer is identified as the location where σ 2 k is 1% of its maximum, where 2 dt/T of the LL (navy blue), LM (dark green), LH (light green), and LVH (yellow green) cases.Circle denotes δs estimated via (3.25).Triangle represents the location where σ 2 k is reduced to 1% of its maximum value.

Conclusions
This paper has examinined the impact of flow pulsation on longtime-and phaseaveraged flow statistics in an open-channel flow over urban-like roughness.A series of LES of pulsatile flow past an array of cuboid elements has been carried out, programmatically varying the oscillation amplitude (α) and frequency (ω).The forcing frequencies have been chosen as a multiple of the characteristic frequency of turbulence in the UCL and encompass a range of values representative of submesoscale motions (Mahrt & Bou-Zeid 2020).The main findings and contributions of this study are outlined below.
(i) Flow pulsation leads to an increase of the z 0 parameter educed from longtimeaveraged u 1 profiles, with larger α and ω values yielding a larger z 0 .On the contrary, d was found to be insensitive to variations in α and ω.The increase of z 0 was shown to be a direct consequence of the quadratic relation between the phase-averaged canopy drag D and the phase-averaged velocity u 1 , and this relation was leveraged to construct a phenomenological model for z 0 .The proposed model takes surface information and the variance of the phase-averaged velocity in the UCL as input parameters and captures the impact of flow unsteadiness on the z 0 parameter in the absence of flow reversal.
(ii) The wall-normal distributions of the longtime-averaged shear stress and canopy drag are unaltered by the flow unsteadiness.In contrast, the same cannot be said for longtime-averaged resolved normal Reynolds stresses, especially in the UCL.In particular, k profiles were found to be relatively more sensitive to variations in α via the longtime-averaged shear production of k.The highest frequency cases were also characterized by a relatively more isotropic turbulence field in the UCL, owing to a more efficient kinetic energy redistribution by the pressure-strain terms.
(iii) The oscillation amplitudes of phase-averaged streamwise velocity and resolved Reynolds stresses scale with α.This behavior is due to the fact that the nonlinear production terms of u 1 u 3 and k are of relatively modest magnitude when compared to the linear ones.Increasing the pulsation amplitude might lead to more substantial contributions from nonlinear production terms and break down this scaling.
(iv) For each case, profiles of oscillatory shear rate and resolved Reynolds stresses are characterized by oscillating waves which are advected away from the UCL at a constant speed while also being dissipated.ω is found to determine both the speeds of the oscillating waves and the extent of these waves, which identifies the Stokes layer thickness (δ s ).More specifically, δ s was found to increase with decreasing ω, whereas the wave speed increased with ω.The scaling of δ s has also been discussed, and findings have been used to propose a model for δ s .
All in all, flow pulsation is found to have a significant impact on both longtimeaveraged and phase-averaged flow statistics, with nuanced dependencies on oscillation amplitude and frequency.The observed enhancement of the longtime-averaged surface drag, the isotropization of turbulence in the UCL, and the presence of a Stokes layer, amongst others, are expected to have important implications on the exchange of mass, energy, and momentum between the land surface and the atmosphere, as well as affect our ability to model these processes in weather forecasting and climate models.These models typically rely on surface flux parameterizations and theories that are based on flow stationarity assumptions and are not able to capture these behaviors correctly (see, e.g., Stensrud 2007).The proposed phenomenological models for z 0 and δ s , as well as the identified scaling of phase-averaged flow statistics, contribute to advancing our understanding of flow unsteadiness in the ABL and offer a pathway for the development of improved surface flux parameterizations.Given the massive parameter space of unsteady ABL flow processes, it is also essential to acknowledge that several questions remain unanswered and deserve further investigation.For example, what is the impact of different types of periodic and aperiodic flow unsteadiness on turbulence statistics and topology?How are these variations in the structure of turbulence impacting landatmosphere exchange rates of momentum, energy, and mass?Can prevailing surface flux parameterizations be modified to account for these impacts?Addressing these questions will be the subject of future studies.
Declaration of Interests.The authors report no conflict of interest.length scale.The grid resolution of DNS400 is (n 1 , n 2 , n 3 ) = (64, 64, 128) per cube.Such a grid resolution ensures that the ratio between the grid size ∆ = 3 √ ∆ 1 ∆ 2 ∆ 3 and the Kolmogorov scale η does not exceed 2, which has been proven sufficient for DNS of flow over fully rough surfaces (Zhang et al. 2022).
In both simulations, the contribution from the tangential stresses at the cube facets and lower surface to the total surface drag remains below 1%, confirming that the flow is in the fully rough regime.Figure 22 and 23 display the phase-averaged velocity ( u 1 ) and turbulent kinetic energy ( k /u 2 τ ), respectively.Profiles from the LES400 case are in good agreement with corresponding DNS quantities, with the maximum error in the LES400 profiles relative to those from the DNS400 case being approximately 1% and 6% for ( u 1 ) and k , respectively.The minor mismatches in k can be partly explained by the fact that the SGS contribution to k is zero for the LES400 case.It is suggested that, although the equilibrium assumption does not hold in a strict sense, the use of an equilibrium wall-layer model does not result in a noticeable impact on model results for the considered unsteady flow cases.

Appendix B. Resolution sensitivity analysis
To identify grid resolution requirements for this study, a sensitivity analysis has been conducted for the stationary flow case, i.e., α = 0.The domain size for this analysis is (36h, 12h, 4h), and we have studied the convergence of u 1 , u 1 u 1 , and u 1 u 3 profiles as the grid stencil is progressively reduced.Three grid resolutions have been considered, namely (4, 4, 8), (8, 8, 16), and (12, 12, 24) on a percube basis.Note that the reduced domain size may have an impact on the evaluated flow statistics, but this serves the purpose of this analysis, since we are here only interested in quantifying relative variations of selected profiles across grid resolutions.Other numerical and physical parameters of the grid-sensitivity analysis simulations are set equal to the ones used in the main simulations.
As apparent from figure 24, profiles from the (8, 8, 16) case are essentially matching corresponding ones from the (12, 12, 24) case, indicating that the chosen grid resolution for the pulsatile channel flow analysis is sufficient to yield grid-independent flow statistics (up to second order).

Figure 1 .
Figure 1.Side (a) and planar view (b) of the computational domain.

Figure 2 .
Figure 2. Instantaneous streamwise fluctuating velocity field at streamwise/cross-stream plane x3 = 2h (a, b) and x3 = 0.75h (c, d ) from the HM case.Panels a and c correspond to t = 0, whereas panels b and d to t = π/2.

Figure 3 .
Figure3.Wall-normal profiles of longtime-averaged streamwise velocity u1 of high-amplitude cases (a) and low-amplitude cases (b).Line color specifies the forcing frequency: navy blue, ωT h = 0.05π; dark green, ωT h = 0.125π; light green, ωT h = 0.25π; yellow-green, ωT h = 1.25π.u1 from the SS case, represented by the red dashed line, is included for comparison.Black solid line indicates the slope of 1/κ, with κ = 0.4.The shaded area highlights the fitting region for the estimation of the roughness length scale z0.

Figure 5 .
Figure 5. Phase-averaged canopy drag D at x3/h ≈ 0.8 as a function of the local phase-averaged velocity u1 (a) and the same D − u1 plot but with the time lag between D and u1 removed (b).Different colors are used to denote different forcing frequencies: navy blue, ωT h = 0.05π; dark green, ωT h = 0.125π; light green, ωT h = 0.25π; yellow-green, ωT h = 1.25π.Different symbols are used to distinguish between oscillation amplitudes: triangle, α = 0.2; cross, α = 0.4.Red square represents the SS case.Red dashed line represents D = C d λ f u1 | u1 | with C d obtained from the SS case.

Figure 6 .
Figure 6.Time evolution of D (black solid lines) and u1 (red dashed lines) at x3/h ≈ 0.8 from the HL (a) and LVH (b) cases.The acronyms of LES runs are defined in table 1.

Figure 7 .Figure 8 .
Figure 7.Comparison between z0 estimated via (3.11) and z0 from LES. Symbols and colors correspond to those used in figure 5.

Figure 9 .
Figure 9. Normalized shear production terms of k: P1 (a) and P2 (b).Symbols and line colors correspond to those used in figure 5.

Figure 13 .
Figure 13.Profiles of the oscillatory resolved Reynolds shear stress u 1 u 3 from the LL (blue) and HL (red) cases at different phases of the pulsatile cycle.Profiles are T /4 apart.

Figure 15 .
Figure 15.Normalized production terms for u 1 u 3 (a) and for k (b) at x3/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.

Figure 16 .
Figure 16.Space-time diagrams of u1/uτ from the LL (a), LM (b), LH (c), and LVH (d ) cases.Horizontal dashed lines identify the top of the UCL.

Table 1 .
List of LES runs