2.1 Numerical simulations set-up

To study mean–turbulence interactions, we employ the LES technique, which resolves the large scales of motion and parameterizes eddies smaller than the filter size. In the code, the mean pressure gradient is represented as an equivalent horizontal geostrophic wind vector with components $U_{g}$ and $V_{g}$ , as shown in figure 1. All of the results are normalized by the maximum magnitude of the geostrophic wind vector, defined as $G=8~\text{m}~\text{s}^{-1}$ , and the BL depth $z_{i}$ ( ${\approx}1000$  m in our simulations). The molecular terms are ignored due to the very high $Re$ regime assumption. This also requires the use of a rough wall model, and as such the viscous length scale and Stokes-layer thickness are smaller than the grid–filter size, and are not being resolved in our LES.

The LES code solves the filtered incompressible continuity and Navier–Stokes equations in dimensionless form (rescaled by $G$ and $z_{i}$ ), where a hydrostatic approximation is assumed to hold for the mean flow in the vertical direction (mean vertical pressure gradient and gravity cancel out):

The tilde denotes filtering at the grid/filter scale; $x_{i}$ is the position vector; $\tilde{u} _{i}=(\tilde{u} ,\tilde{v},\tilde{w})=(\tilde{u} _{1},\tilde{u} _{2},\tilde{u} _{3})$ is the resolved velocity vector; $\unicode[STIX]{x1D70F}_{ij}$ is the deviatoric part of the subgrid scale (SGS) stress tensor $\unicode[STIX]{x1D70E}_{ij}$ ; $p^{\ast }$ is a modified pressure into which the resolved and SGS TKE are lumped (see § 2.3); $f_{c}$ is the Coriolis parameter, the dimensional value of which is taken as $1.394\times 10^{-4}$  Hz (corresponding to $73^{\circ }$ latitude); and $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. For all cases except one, $\unicode[STIX]{x1D70F}_{ij}$ is parameterized using the Lagrangian scale-dependent dynamic (LASD) Smagorinsky model, which is detailed and validated for neutral flows in Bou-Zeid, Meneveau & Parlange (Reference Bou-Zeid, Meneveau and Parlange2005). The static Smagorinsky model (with an imposed coefficient that is damped near the surface as described by Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005) is used in one of the cases to assess the influence of the SGS closure on the TKE budget computations.

Figure 1.

Schematic view of the simulations: $t$ is time; $\unicode[STIX]{x1D70F}_{f}$ is the period of the variable pressure forcing; $\unicode[STIX]{x1D70C}$ is density; $f_{c}$ is the Coriolis parameter; $z_{0}=0.01$  m is the surface roughness length; $L_{z}=1000$  m is the domain height.

The code uses a pseudo-spectral scheme in the horizontal directions to compute derivatives and for test filtering in the dynamic SGS model. A second-order finite-difference scheme is used in the vertical direction, necessitating a staggered grid. Therefore, the vertical velocity $w$ in our simulations is computed and stored halfway in between the other variables ( $u$ , $v$ and $p$ ). This means that $w$ grid points begin at the wall and end at the top of the domain at $z=L_{z}$ . The other variables are computed on a grid that starts at $\text{d}z/2$ above the wall and continues up to $L_{z}+\text{d}z/2$ . This implies that some vertical interpolation will be necessary to compute some of the TKE budget terms later. The vertical velocity, $w$ , is set to 0 at the surface. The top boundary conditions for the velocity are no stress, $\unicode[STIX]{x1D70F}_{i3}=0$ , and an impermeable wall, $w=0$ . The boundary conditions for the pressure field at both the top and bottom surfaces are found from consideration of the vertical momentum balance equation across these boundaries, to be consistent with the zero-vertical-velocity condition, yielding $\unicode[STIX]{x2202}\tilde{p}^{\ast }/\unicode[STIX]{x2202}x_{3}+\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{j3}/\unicode[STIX]{x2202}x_{j}=0$ , i.e. zero net force. We impose the surface shear stresses through a local-similarity theory formulation: the filtered values of the horizontal velocities (at twice the grid scale) are used to compute the stresses through a local log-law to ensure that the average stresses over the wall are close to the classic mean log-law predictions (see Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005) for more details, performance and validation of the wall model used). Time advancement uses the second-order Adams–Bashforth scheme. The code has been validated in spatially variable flows (e.g. Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005); Li, Bou-Zeid & Anderson (Reference Li, Bou-Zeid and Anderson2016)) that impose similar challenges for the SGS model to the temporally variable ones simulated here; as such, no further validation is deemed necessary. The reader is referred to Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005) and Momen & Bou-Zeid (Reference Momen and Bou-Zeid2016) for more details about the code as used in the present simulations.

Table 1.

Details of the LES cases. Here, $\unicode[STIX]{x0394}t$ is the time-averaging period (2–10 min) and ‘duration’ denotes the dimensional simulation time. The variable forcing, non-dimensionalized by $G=8~\text{m}~\text{s}^{-1}$ , has the general form $U_{g}=\cos (2\unicode[STIX]{x03C0}t/\unicode[STIX]{x1D70F}_{f});V_{g}=\cos (2\unicode[STIX]{x03C0}t/\unicode[STIX]{x1D70F}_{f})$ . In unsteady cases, $u_{\ast }$ is calculated during the last four inertial oscillations, $\unicode[STIX]{x1D70F}_{EBL}=12.5$  h, of the simulations. The domain height $z_{i}=1000$  m in all of our simulations.

* See § 4 for details on the resonance in these flows.

It should be noted that usually (as long as $Ro$ is not too large) the Coriolis frequency characterizes the response time of the mean velocity in geophysical flows with rotation, and thus the natural inertial time scale of the EBL would be $\unicode[STIX]{x1D70F}_{EBL}=2\unicode[STIX]{x03C0}/f_{c}\approx 12.5$  h (Tennekes & Lumley Reference Tennekes and Lumley1972). The largest turbulence time scale can also be estimated via $\unicode[STIX]{x1D70F}_{t}=z_{i}/u_{\ast }\leqslant 2$  h, where $u_{\ast }$ is the friction velocity ( ${\approx}0.2~\text{m}~\text{s}^{-1}$ ). This is the turnover time of the largest eddies which will be approximately as large as the domain height $z_{i}$ , but the full spectrum of eddies would be characterized by a range of time scales reaching down to the Kolmogorov time scale (although in LES the smallest scales are not resolved). We will examine the Ekman flow and turbulent BL behaviour for different forcing time periods, relative to the mean and turbulence time scales (i.e. under various values of $\unicode[STIX]{x1D70F}_{f}/\unicode[STIX]{x1D70F}_{EBL}$ and $\unicode[STIX]{x1D70F}_{f}/\unicode[STIX]{x1D70F}_{t}$ ). The time scale that emerges from the non-dimensionalization of the equations ( $z_{i}/G\approx 125$  s) is not the most relevant to the physics of the present flow (although it is related to the eddy turnover time). Hence, we normalize the simulation time using the inertial time scale of the EBL, $\unicode[STIX]{x1D70F}_{EBL}$ , which is constant in all of the cases.

Table 1 describes the eight simulations that we perform, as well as their acronyms which will be used hereafter. The temporal resolution of all cases is 0.15 s (i.e. equivalent to a non-dimensional value of $1.2\times 10^{-3}$ when normalized by $z_{i}/G$ ) and is of the order of the turnover time of the smallest resolved eddies; it is primarily set to maintain numerical stability and satisfy the Courant–Friedrichs–Lewy condition (Madabhushi & Vanka Reference Madabhushi and Vanka1991). The term $\unicode[STIX]{x0394}t$ in table 1 stands for the temporal averaging period used in the analyses; some temporal averaging is useful to obtain converged statistics, but its time span must be kept much shorter than any of the relevant time scales of the problem in these unsteady simulations. The Ekman boundary layer thickness ( $u_{\ast }/f$ ) in most of our simulations is less than 1000 m, and thus we resolve this height in such cases. The largest Ekman height occurs in the resonance case, in which it can reach up to 4000 m. However, this height is defined as the height at which the velocity vector is aligned with the geostrophic wind, and its influence on the flow dynamics is not clear (in the real ABL, for example, an inversion could prevent mixing up to that height). To verify whether the domain should extend up to the maximum $u_{\ast }/f$ , we conducted another LES of the resonance case with a doubled domain height ( $L_{z}=2000$  m), while maintaining the same resolution by doubling the number of vertical grid points. The results indicated that the differences between the two simulations were not significant, and therefore our basic domain height of 1000 m is sufficient. Furthermore, since in our simulations the imposed pressure gradients are unsteady, the friction velocity will be time-dependent, and thus $u_{\ast }/f$ will also vary in time and with the exerted forcing (see table 1). For instance, this length scale varies between 200 and 700 m in the last two inertial oscillations of the fast and very fast oscillation cases. Therefore, the $L_{z}$ , which sets the upper limit on the size of turbulent eddies, is preferred over the $u_{\ast }/f$ for normalization of the height, and we set it to 1000 m in all of our simulations.

All of the simulations were performed on a $128^{3}$ grid cell numerical mesh. Thus, the vertical spatial resolution of the domain was 1000/128 m and the horizontal resolutions were $2\unicode[STIX]{x03C0}\times 1000/128$  m. To test for grid convergence, we also ran some of our simulations (RO and VFO) with a higher resolution ( $256^{3}$ ), and the results obtained were similar to the basic $128^{3}$ resolution (not shown); therefore, the educed dynamics and conclusions are not influenced by the numerical set-up.

The first two cases in table 1 are to elucidate the baseline steady-state time variations of the ( $x$ – $y$ ) planar-averaged TKE budget and to depict the resolved TKE profiles. We also use them to assess the quality of the SGS modelling and how it influences the closure of the TKE budget. Cases 3–6 are performed to study the dynamics of the mean and turbulence with different forcing variability time scales, from slow to very fast. The forcing magnitude changes in these cases (3–6); its direction, however, remains always aligned along $y$ (inducing a geostrophic wind along  $x$ ). In case 7, on the other hand, the pressure-gradient magnitude is fixed but its orientation rotates in time (i.e. the pressure is varying with time in both the $x$ - and $y$ -directions). The last case is with an unsteady mean pressure forcing but without a Coriolis force (i.e. $f_{c}=0$ ); we ran it to investigate the influence of (lack of) rotation on the mean and turbulence response. The long simulations are primarily designed to understand the underlying dynamics of the Ekman flows (ensure independence from initial conditions and reproducibility of the statistics over multiple cycles). Although the forcing conditions are only expected to persist in nature for shorter periods ( ${\sim}$  h), their dynamics can also be investigated from our cases. Such extended periodic forcings are not uncommon in engineering applications (e.g. the UNR case, which has a similar TKE budget dynamics to the other cases with Coriolis, as we show later).

2.2 Large-eddy simulation initialization

A 50 h $64^{3}$ node steady simulation is performed to initiate the spin-up. Then, the outputs are linearly interpolated into the $128^{3}$ node domain and used as initial values for cases 1 and 2. The rest of the simulations (except UNR) use the final output of case 1 as their initial condition. A similar approach is used for initializing the UNR case, starting from an independent 50 h $64^{3}$ node simulation.

2.3 Resolved TKE equation

The pressure field that is computed in our incompressible LES code is not a thermodynamic variable, but rather a dynamic pressure that is modified to include the resolved and SGS TKE components that act normal to a surface (like pressure):

This is done to improve the conservation of resolved energy and to allow the modelling of the deviatoric part of the SGS stress only. If we substitute this pressure expression into (2.2), and apply the Reynolds decomposition and averaging, the conventional forms of the mean kinetic energy $\overline{E}=\overline{u}_{i}\overline{u}_{i}/2$ and the resolved turbulent kinetic energy $\overline{e}=\overline{u_{i}^{\prime }u_{i}^{\prime }}/2$ arise:

where the overbar denotes the ensemble average and the prime denotes the turbulent fluctuation from the ensemble average. The hydrostatic approximation for the mean allows us to omit the gravity and $\unicode[STIX]{x2202}\overline{p}/\unicode[STIX]{x2202}z$ terms in (2.4) (hence only the $\unicode[STIX]{x1D6FF}_{i1}+\unicode[STIX]{x1D6FF}_{i2}$ in the pressure term). It should be noted that all $u_{i}$ , $p$ and $S_{ij}$ in the above equations are the resolved components, and we will omit the tildes for notational simplicity. Since the SGS kinetic energy $\unicode[STIX]{x1D70E}_{kk}$ would require a model that is distinct from our model for the deviatoric part $\unicode[STIX]{x1D70F}_{ij}$ , and since it makes a small contribution to the terms in (2.4) and (2.5) relative to its resolved counterpart away the wall (this is verified using a simple model for $\unicode[STIX]{x1D70E}_{kk}$ in appendix A), we neglect the terms that include this component rather than risk introducing modelling errors that could be of the same order as the magnitude of the terms involving $\unicode[STIX]{x1D70E}_{kk}$ . Furthermore, in our TKE budget calculations, we can safely assume no subsidence, $\overline{w}=0$ , and periodic statistical homogeneity in both horizontal directions; hence, equations (2.4) and (2.5) reduce to

We compute the TKE budget in terms of the resolved variables and define the following expressions for the terms in (2.6) and (2.7) for the rest of the paper as follows:

where $P$ , $PT$ , $RT$ , $ST$ , $T$ and $\unicode[STIX]{x1D716}$ are respectively the shear production, pressure transport (correlation), resolved turbulent transport, SGS turbulent transport, total TKE transport and SGS dissipation. Our LES code computes each of these terms by averaging in space and time. To obtain the averaged quantities on the right-hand side of (2.8), first, a spatial average over the horizontal directions is applied to acquire a single instantaneous profile for each term (to define the ensemble average). Then, all of these TKE budget profiles are calculated instantaneously (this involves products and vertical derivatives), and finally the profiles of the budget terms are averaged in time over the short $\unicode[STIX]{x0394}t$ intervals indicated in table 1. Time averaging before computing the terms in (2.8) was found to cause significant errors in this unsteady framework due to the time variation of the parameters that need to be Reynolds-averaged. Hence, only spatial averaging is used as a surrogate for Reynolds averaging, although the resulting unsteady statistics are then averaged over $\unicode[STIX]{x0394}t$ . A similar procedure is applied to produce velocity and other time series.

It should be noted that because of the staggered set-up of the numerical domain, we linearly interpolated all terms from the $w$ nodes into $u{-}v{-}p$ nodes when needed, and subsequently applied the multiplications and differentiations in (2.8). This could be a source of error in the results, but as the numerical grids become finer this error would be expected to decrease.

4.1 Turbulence intensity variations

Turbulence typically has a smaller time scale than the mean flow; thus, it is expected that the turbulence intensity (TI) should adjust itself to the variations in the mean very rapidly. That is, turbulence is usually near a quasi-equilibrium state with the mean flow, and TI should remain constant for a given flow. However, if the forcing of the mean varies at time scales commensurate with those of turbulence, the adjustment might not be sufficiently fast. This could result in variations in the turbulence statistics, even the ones normalized by mean variables such as TI. This departure from quasi-equilibrium is illustrated in figure 4 for the varying pressure-gradient simulations: TI fluctuates strongly compared with the steady-state simulation. It should be noted that since the simulation time of the steady case SF is $0.8\unicode[STIX]{x1D70F}_{EBL}$ , its averaged value is extrapolated up to $t=8.0\unicode[STIX]{x1D70F}_{EBL}$ in the figure for comparison with the longer unsteady runs. Since the forcing is not changing in this case, the variations in TI are expected to remain similar to the simulated 10 h.

Figure 4.

Turbulence intensity for five simulations, (a) in the middle of the EBL and (b) near the wall. The peaks of the FO case in (b) rise up to 0.8.

The TI in the slow-oscillation case displays a periodic behaviour in which the magnitude of the oscillations is approximately three quarters of the mean value at $z/z_{i}=0.59$ (figure 4 a), despite the fact that the forcing time scale in this case is two times the inertial time scale. If one compares the TI minima and maxima for SO with the mean horizontal velocity $M$ in figure 5(a), opposite trends can be observed (as one peaks the other troughs). The very-fast-oscillation case also shows variations in the TI times series (up to 20 % of the mean) at $z/z_{i}=0.59$ . However, the strongest time variability occurs in the fast-oscillation FO case (this will be explained in the next subsection).

Figure 5.

The $M$ (a,c,e,g) and TKE budget (b,d,f,h) time series for case 3 SO (a,b), 4 RO (c,d), 5 FO (e,f) and 6 VFO (g,h).

Figure 4 also indicates that the TI variations near the surface ( $z/z_{i}=0.055$ ) are lower for the SO, RO and SF cases than near the middle of the domain ( $z/z_{i}=0.59$ ). This is due to the quicker response (smaller effective $\unicode[STIX]{x1D70F}_{t}$ ) of the smaller eddies near the surface to the imposed variations and their faster equilibration. On the other hand, when the imposed forcing time scale is of the order of or smaller than the turbulence time scale (FO and VFO), the eddies might not find adequate time to respond, and therefore the variations in TI do not decrease as height decreases. Nonetheless, the results suggest that while the FO case depicts strong variability as expected, the VFO case seems not to depart strongly from equilibrium conditions. These high variations in the surface layer emerge when $M$ gets close to zero while the TKE, which decays slowly and can also be transported from above, remains high in the EBL.

These simulations imply that the TI is rarely constant in unsteady flows; it changes according to the mean pressure variations even when those variations occur at diurnal time scales ( ${\sim}24~\text{h}=1~\text{day}$ ) two times larger than the inertial period of the EBL in the paper ( ${\sim}12$  h). As the forcing time scale increases even further, a quasi-steady TI would eventually be expected. These variations could either imply non-equilibrium, or they could be due to the changes in the $Ro$ and the effective $Re$ of the LES (note that the LES has an effective $Re$ imposed by the SGS model and/or the grid resolution, despite the fact that we ignore the molecular terms). We will revisit this point in § 4.4.

4.2 Turbulent kinetic energy budget variations

The unsteady forcing alters the mean and turbulent velocity time series differently, as can be concluded from the TI plots in the previous subsection. We seek to understand how these different responses of the mean and turbulence are modulated by their dynamics and their interaction. To that end, the time series of the dimensionless mean wind velocity at three heights and the dimensionless TKE budget terms at two heights are plotted. Figure 5(a,b) exhibits case SO. One can observe that all of the TKE budget terms, at a given height, peak and trough at the same time. Furthermore, the peaks in the shear production correspond in time to the peaks in $M$ , albeit with a small lag relative to the mean wind (see § 6), and a similar behaviour is observed in the troughs. This indicates that the turbulence and mean dynamics are well correlated, although they are not perfectly synchronized.

The next simulation, presented in figure 5(c,d) is the resonant-forcing case ( $\unicode[STIX]{x1D70F}_{forcing}\approx \unicode[STIX]{x1D70F}_{EBL}$ ), which illustrates the increased TKE levels as the mean resonates and increases beyond the geostrophic wind speeds ( $M/G>1$ , figure 5 c). As detailed in Momen & Bou-Zeid (Reference Momen and Bou-Zeid2016), this can be conceptually explained by considering the damped-oscillator (mass–spring–damper) model for the response of the mean flow in the EBL to a variable forcing. In that paper, we showed that the mean dynamics can be simply represented by two coupled second-order ordinary differential equations (ODEs) for the means of $u$ and $v$ , like a damped oscillator (underdamped for the EBL). As such, the mean dynamics are characterized by a single natural time scale $\unicode[STIX]{x1D70F}_{EBL}$ , and they damp or magnify the magnitude of the forcing (transmissibility ratio) depending on the ratio $\unicode[STIX]{x1D70F}_{f}/\unicode[STIX]{x1D70F}_{EBL}$ .

Unlike the first two cases in the figure for which both $\unicode[STIX]{x1D70F}_{EBL}$ and $\unicode[STIX]{x1D70F}_{f}$ are much larger than $\unicode[STIX]{x1D70F}_{t}$ , the last two cases (figures 5 e,f and 5 g,h) demonstrate decaying turbulence with a fast (4 times the natural EBL frequency, FO) and a very fast (37.5 times the natural EBL frequency, VFO) variable pressure forcing. In the FO case (as well as RO), we observe that the peaks and troughs of the mean and TKE production time series remain well correlated, with some lag. In the very-fast-forcing scenario (faster than $\unicode[STIX]{x1D70F}_{t}$ ), the production, transport and dissipation terms decay with the mean wind velocity, but there appear to be fast variations in the mean that are not reflected in the TKE budget terms.

As illustrated in figure 5, when $\unicode[STIX]{x1D70F}_{f}\geqslant \unicode[STIX]{x1D70F}_{EBL}$ , the mean adjusts to the forcing and varies with a single time scale $\unicode[STIX]{x1D70F}_{f}$ . Since that implies that $\unicode[STIX]{x1D70F}_{f}\gg \unicode[STIX]{x1D70F}_{t}$ , the turbulence has time to adjust and follow the mean. However, when $\unicode[STIX]{x1D70F}_{f}<\unicode[STIX]{x1D70F}_{EBL}$ , the dynamics of the mean continues to exhibit variability at the forcing scale $\unicode[STIX]{x1D70F}_{f}$ , but now it also varies at its natural time scale $\unicode[STIX]{x1D70F}_{EBL}=12.5$  h. Since a damped oscillator damps input frequencies much larger than its natural frequencies (the higher the frequency the stronger the damping), in the last two cases the mean decreases in time because it is not able to fully absorb the energy provided by the forcing. Moreover, in these cases, $\unicode[STIX]{x1D70F}_{f}\sim \unicode[STIX]{x1D70F}_{t}$ , and therefore turbulence might not adjust rapidly enough. Turbulence, unlike the mean, has a broader spectrum and therefore a broader range of characteristic time scales (corresponding to the range of eddy sizes). It thus acts as a broadband filter that filters the very fast oscillations occurring at a time scale $\unicode[STIX]{x1D70F}_{f}=0.2\unicode[STIX]{x1D70F}_{t}$ in the VFO case, but responds to the slower decrease in the mean. In the FO case, the mean variability is at a time scale ${\geqslant}\unicode[STIX]{x1D70F}_{f}=1.5\unicode[STIX]{x1D70F}_{t}$ , and turbulence responds (and does not filter) to this variability as well as to the slower decrease in the mean. However, it is unlikely to be in equilibrium given that $\unicode[STIX]{x1D70F}_{f}\sim \unicode[STIX]{x1D70F}_{t}$ (we will investigate equilibrium and its criteria further in § 4.4).

Figure 6.

Space–time graphs of the dimensionless shear production (multiplied by $z_{i}/G^{3}$ ) for cases (a) 3 SO, (b) 4 RO, (c) 5 FO and (d) 6-VFO.

To further examine the turbulence dynamics in the cases of figure 5, spacetime graphs of the shear-production terms are plotted in figure 6. The production oscillations are correlated with the mean velocity fluctuations, and the magnitude of this term decreases in height (except perhaps some parts of the FO case which will be discussed in § 4.4). Another important feature to observe in these graphs is the phase lag of the sinusoidal behaviour of the TKE budget at different heights. One could notice approximately 10–30 min phase lag, for example, between the time series extrema at $z/z_{i}=0.055$ and 0.234 in figure 5 (the higher elevations lag the lower ones), which is related to the faster response of the smaller eddies near the surface compared with the larger eddies further aloft, as well as to the time required for the shear-produced turbulence at the surface to be transported to higher levels. For instance, a lag of approximately 30 min between the two shown heights can be observed in the production and dissipation terms of figure 5(d) at $t=3.2\unicode[STIX]{x1D70F}_{EBL}$ and the minima of figure 5(b) at $t=3.6\unicode[STIX]{x1D70F}_{EBL}$ . These lags are more apparent in figure 6, where they are illustrated by the vertical tilt of the patterns of strong production, which are particularly clear in the FO case (figure 6 c). We will quantitatively calculate this lag between two heights in § 6.

The role of the variability in the forcing direction, at a constant magnitude, is also explored. This is shown in figure 7, which depicts the velocity (figure 7 a) and TKE budget (figure 7 b) time series of the 2D fast oscillation. As in the 1D FO case, various TKE terms are directly responding to the mean velocity oscillations, and the mean decreases with time due to the filtering of the forcing by the damped-oscillator dynamics. Therefore, time variabilities in the magnitude or the direction of the forcing seem to result in comparable responses and mean–turbulence interactions.

Figure 7.

The $M$ (a) and TKE budget (b) time series for case 7 FO_2D. The TKE budget terms are as defined in the legend of figure 5(h). The error bars on the production term are the residuals of the budget, reflecting the uncertainty in the budget terms.

4.3 The role of the Coriolis force

The final case investigated (figure 8) illustrates the influence of the Coriolis force on the dynamics of the EBL. This simulation has the same variable forcing as the RO case, but the mean flow responds very differently and no resonance is observed. This distinct behaviour is due to the fact that the 12.5 h forcing in RO matches the mean EBL time scale ( $\unicode[STIX]{x1D70F}_{EBL}=2\unicode[STIX]{x03C0}/f_{c}\approx 12.5$  h), while it does not represent any special frequency for the UNR case where the Coriolis effect is not present. Hence, such a forcing time scale in the RO case leads to the resonance phenomena and to an increase in the magnitude of the velocity and TKE budgets, whereas in the UNR case the velocity decreases exponentially. This exponential decay can be exactly predicted by the solutions of the two coupled ODEs which form a reduced model of the Reynolds-averaged Navier–Stokes (RANS) equations, as derived in Momen & Bou-Zeid (Reference Momen and Bou-Zeid2016) by setting $f_{c}=0$ (which results in a mass–damper system with no spring). Despite the significant differences in the mean response, the TKE budget terms peak and fall with a small lag relative to $M$ as in the EBL simulations, suggesting that turbulence is only influenced by the mean variability it experiences, regardless of the driving forcing variability that modulates the mean. This is consistent with the observation that in the MKE equation (2.6), the forcing variability appears through $\overline{u_{i}}\unicode[STIX]{x2202}\overline{p}/\unicode[STIX]{x2202}x_{i}$ , while in the TKE equation (2.7) the forcing variability does not appear directly. Thus, the forcing pressure gradient interacts with the mean, which then interacts with turbulence. However, turbulence is not directly modulated by the variability in the pressure forcing. Therefore, these cases indicate that in general the turbulence and the mean interactions are similar in flows with and without rotation (as long as $Ro$ is not too small, since rotation can damp turbulence at small $Ro$ ), but the dynamics of the mean is strongly affected by the Coriolis force. In other words, the Coriolis force only directly influences the mean dynamics (2.4), while it does not explicitly appear in the TKE budget equation (2.5) and has only an indirect influence on the TKE. Hence, the main difference between the UNR case (analogous to a mass–damper system) and other cases in this paper (analogous to mass–spring–damper systems) is in their mean response to the imposed forcing, while their TKE budget dynamics are expected to be similar according to (2.5), which is confirmed by our LES results.

Furthermore, the simulations of this case are in agreement with the LES results of other papers, which have investigated oscillatory flows driven by harmonic pressure gradients without rotation. For instance, the production and dissipation terms in figure 8(b) show one complete sinusoidal cycle at each imposed pressure-gradient sinusoidal oscillation (without the appearance of other oscillatory time scales such as the $\unicode[STIX]{x1D70F}_{EBL}$ ). This is similar to the results of Salon et al. (Reference Salon, Armenio and Crise2007), although in our simulation the TKE budgets, as well as $M$ , are decreasing from the initial conditions to reach the statistically steady-state oscillations (where the crests of $M$ do not vary between cycles).

Figure 8.

The $M$ (a) and TKE budget (b) time series for case 8 UNR. The solid lines in (b) show the budgets for $z/z_{i}=0.234$ and the lines with dots are for $z/z_{i}=0.055$ . The error bars on the production term are the residuals of the budget, reflecting the uncertainty in the budget terms.

4.4 Turbulence quasi-equilibrium and normalized TKE budget

We define the ‘quasi-equilibrium with the mean’ as the conditions where the turbulence statistics of a flow at any point in time are equivalent to those that would be generated if the current mean values at that same point in time were steady. This definition can be directly assessed from the present simulation results. For instance, figure 9(a) displays the mean velocity and the production and dissipation time series for both steady-state and slow-oscillation cases. As figure 9(a) indicates, when the mean velocities of both simulations are the same, the production and dissipation terms have the same values (the oscillating and straight horizontal lines intersect), implying that the turbulence is in quasi-equilibrium with the mean. In such conditions, if the turbulence statistics are normalized properly with the mean, they will collapse independent of the imposed variations. Examples of this non-dimensionalization with the mean velocity for three cases are shown in figure 9(c). We will use this normalization to develop a turbulence closure model based on the mean variables in § 6. Nevertheless, in the FO case, the production term departs from its steady-state equivalent values (figure 9 b), and hence the turbulence is out of equilibrium with the mean, and normalization with the mean yields a highly variable times series (figure 9 d).

Given that turbulence varies with the unsteady mean, and some mean-normalized quantities like TI vary in time, we will also normalize the TKE budget terms by a turbulence statistic of order 2 (TKE production) to further examine the departure of turbulence from quasi-equilibrium. Production–dissipation balance is often used to assert turbulence equilibrium (Pope Reference Pope2000); we will also examine the transport terms here.

When the forcing time period is larger than the turbulence time scale $\unicode[STIX]{x1D70F}_{f}\gg \unicode[STIX]{x1D70F}_{t}$ , if quasi-equilibrium holds, it is expected that the production-normalized TKE budget terms would only be functions of height, independent of the forcing function that is imposed on them. The SO and RO cases shown in figure 9(e,f) indeed seem to be close to a quasi-equilibrium, since the variations of the time series of the TKE budget terms are moderate (dissipation/production ratio varies by less than 20 % of its mean signal). This implies that although the TI of these cases varies in time, as shown in figure 4, the TKE budget dynamics equilibrates fast and maintains quasi-steady budget profiles. In other words, the relative shapes of the budget profiles in these cases remain the same, while their magnitude might change in time; the profiles will collapse if properly normalized (e.g. by shear production as we did here, or by the mean velocity as in figure 9 c). This is expected, since the turbulence time scale is smaller than the imposed forcing time period, such that turbulence should have enough time to adapt itself to the exerted slow variations.

Figure 9.

The ‘quasi-equilibrium of the turbulence with the mean’ is examined for case 3 SO (a) and case 5 FO (b). Normalized TKE budget by the mean for cases where the turbulence is in quasi-equilibrium with the mean flow (c) and for case FO lacking quasi-equilibrium (d, the maximum rises up to 0.3). Relative TKE budget time series for cases 3 SO (e), 4 RO (f), 5 FO (g) and 6 VFO (h).

More surprisingly perhaps, when the forcing time period is much smaller than the time scale of the turbulence, $\unicode[STIX]{x1D70F}_{f}\ll \unicode[STIX]{x1D70F}_{t}$ , we still observe the same quasi-equilibrium behaviour. As figure 9(h) shows, the VFO case also exhibits the same dynamics in which the normalized TKE budget profile preserves its shape and does not change much in time. This behaviour is due to the fact that the imposed variable forcing is very fast, such that the turbulence cannot respond accordingly; its effect is thus to integrate or filter the very fast variations of the mean signal over its time scale $\unicode[STIX]{x1D70F}_{t}$ , as discussed above. Only the slow decay of the mean of this case modifies the TKE budget, and turbulence adapts to this gradual change. As we discussed before, the turbulence contains a spectrum of eddies with time scales ${\leqslant}\unicode[STIX]{x1D70F}_{t}$ . Therefore, there exist some smaller turbulent eddies that interact with this very fast oscillatory variation in the mean flow. That is the reason for the relatively higher variation of the normalized TKE budget of the VFO case in figure 9(h) compared with the SO and RO cases in figure 9(e,f), resulting in a weaker quasi-equilibrium state with the mean for the VFO in comparison with other slower cases shown in figure 9(c). Despite this fact, larger more energetic eddies cannot respond to this very fast oscillation, and this results in a filtered TKE production that remains near a quasi-equilibrium state. In other words, the most energetic eddies of the turbulence cannot respond to the very fast oscillations of the imposed forcing and the mean at time scales smaller than $\unicode[STIX]{x1D70F}_{t}$ , and hence it filters out the imposed fast fluctuations of the mean flow and only responds to the filtered signal, which is a slow decay here.

Nevertheless, if the forcing and turbulence time scales are comparable, $\unicode[STIX]{x1D70F}_{f}\approx \unicode[STIX]{x1D70F}_{t}$ , the TKE budget will no longer be in quasi-equilibrium. For instance, figure 9(g) shows that in the FO case, the forcing strongly alters the TKE budget profiles. Such behaviour is more pronounced when the mean velocity decreases, and consequently the shear production is reduced significantly, while other TKE terms do not decrease as fast as production. In the FO simulation, the production term becomes nearly zero even at the surface during some periods, as illustrated in figure 6(c). However, there is always a minimum level of shear production in the VFO (figure 6 d) and other cases (figure 6 a,b). Similar behaviour could also be observed when such oscillation occurs in both directions (as shown for the case FO_2D in appendix B, figure 17). It should be noted that in the FO case, the magnitude of the forcing changes, while in the FO_2D case, only the direction of the forcing varies. Yet, in both cases, the TKE budget profile changes and a quasi-equilibrium turbulence state is no longer observed. Moreover, this out of equilibrium phenomenon, and the subsequent failure of similarity laws like the log-law shown in the next section, will not only occur at this specific forcing frequency. It will occur over a range of forcing frequencies that have comparable time periods to the turbulence time scale; for example, when $\unicode[STIX]{x1D70F}_{f}=6.25$  h, an analogous EBL response will emerge (see appendix B for an illustration).

5.1 Weak to moderate departures from the log-law

The velocity profiles for the unsteady cases that display a moderate departure from the steady-state log-law are first examined at three distinct times in a single cycle, selected close to the end of the simulation: first at the trough of the cycle, then at a quarter period, and finally at a half-period corresponding to the peak. For cases with a moderate departure from the log-law, we plot the three distinct instantaneous profiles for the various cases in figure 11(a–d). The plots show that the normalized wind profile is not constant and varies according to the phase of the imposed unsteady forcing, indicating some moderate departures from the log-law similarity in cases VFO and UNR. The steady-state log-law is denoted by the solid black lines, while the dashed lines show the log-curves fitted to each instantaneous profile and are used to compute the parameters in (5.1); these parameters are provided in table 2 for all of the unsteady simulations. We will consider here that the log-law is not applicable in cases with more than 25 % change in the slope parameter or a tripling in the curvature parameter. These are ad hoc thresholds, but we will use them for a convenient distinction of moderate versus strong departure from the log-law profile.

For RO and SO simulations, the instantaneous ( $\unicode[STIX]{x1D6FE}_{d}$ , $\unicode[STIX]{x1D705}_{d}$ ) values varied in the ranges (1.0–2.0, 0.9–1.1) at different times (figures not shown, see table 2), implying that, based on our criteria, the log-law is ‘applicable’ in these cases in the sense that it can be a useful prognostic tool despite some inaccuracy. The VFO case, on the other hand, shows moderate departures from the log-law, as shown in figure 11(b). Table 2 shows that the curvature parameter in this case is approximately 2–3 times larger than the steady-state case, and the von Kármán constant increases or decreases based on the phase. The $\unicode[STIX]{x1D705}_{d}$ value is higher by approximately 25 % when $t=7.41\unicode[STIX]{x1D70F}_{EBL}$ and lower by 11 % when $t=7.64\unicode[STIX]{x1D70F}_{EBL}$ in the VFO (relative to a value of 1 which would mean that the slope is the same as the steady-state case). However, the trend of increasing wind speed with height persists, despite significant departure from the log-law at some instants. The unsteady pressure forcing with no Coriolis case also shows comparably moderate departures from the log-law profiles. Figure 11(d) and table 2 indicate a 10 % decrease and up to 43 % increase in the slope parameter, with approximately a doubling of the curvature parameter.

Figure 12.

Wind profiles for strong departures from the log-law: case 7 FO_2D (a,b) and case 5 FO (c–f).

Figure 13.

Space–time graphs of wind-profile evolutions for (a) case 5 FO and (b) case 6 VFO.

Table 2.

Log-law departures in seven cases; the cases with an overbar over the name are for averaged profiles. The strong departures from the log-law, for which it is not applicable based on our criteria, are shown in bold.

5.2 Strong departures from the log-law

The fast-oscillation simulations show strong log-law departures. For instance, the FO_2D case exhibits clear log-law failure, as shown in figure 12(b). The $\unicode[STIX]{x1D705}_{d}$ values in table 2 for this case are very far from unity, reaching $-86.4$ at $t=4.66\unicode[STIX]{x1D70F}_{EBL}$ (inverted velocity profile), with a sixfold increase in the curvature parameter. Furthermore, very significant departures are observed in the FO case. Figure 12(c–f) assesses the log-law validity in this case, here for multiple cycles and multiple times in each cycle. The wind profile is drastically different, and the log-law fails to describe it at most times. Figure 13 exhibits the space–time graph for part of the FO simulation (figure 13 a) compared with the VFO case (figure 13 b). Inverted velocity profiles can be observed at $t\approx 6.32\unicode[STIX]{x1D70F}_{EBL}$ and $t\approx 6.68\unicode[STIX]{x1D70F}_{EBL}$ in figure 13(a), while velocity profiles are not inverted in figure 13(b). Moreover, during some periods (e.g. $t\approx 6.28\unicode[STIX]{x1D70F}_{EBL}$ and $6.64\unicode[STIX]{x1D70F}_{EBL}$ in figure 13 a) the magnitude of the velocity is nearly zero throughout the EBL in the FO case. Table 2 shows that the curvature parameter increases significantly in this case. It becomes approximately seven times larger at $t=6.98\unicode[STIX]{x1D70F}_{EBL}$ and $t=7.30\unicode[STIX]{x1D70F}_{EBL}$ . The slope parameter is also modified substantially: at $t=6.66\unicode[STIX]{x1D70F}_{EBL}$ and $t=7.30\unicode[STIX]{x1D70F}_{EBL}$ , it even becomes negative, indicating that the wind velocity is decreasing with height (inverted velocity profiles). This phenomenon usually occurs when $M$ is close to a local minimum and the pressure-gradient force changes the wind direction.

5.3 Log-law validity for the averaged profiles

We also examined the cycle-averaged wind profiles to see whether the log-law is valid if the velocities are averaged in time. Figure 11(e) presents the VFO case averaged over a quarter of its forcing cycle (5 min), and figure 11(f) shows the profiles averaged over a quarter forcing cycle of the FO case (approximately 47 min). The last four rows of table 2 provide the curvature and slope parameters of these profiles averaged over a quarter and a full forcing cycle. The full-cycle-averaged values display less eccentricity from the steady-state reference than the instantaneous profiles, but they still exhibit large departures for the FO case, and consequently the log-law is still not applicable for the FO case (in practical applications, the forcing will not be cyclic and hence time averaging to recover a log-law is even less feasible). The full-cycle-averaged profiles of the VFO case match the log-law reasonably well.

To summarize, these strong departures in the instantaneous and averaged wind profiles indicate that the classic log-law fails to accurately describe the wind profile in these fast-oscillating unsteady cases. In particular, when the imposed forcing time scale is of the order of the turbulence time scale, where the interactions between the mean and turbulence are strongly out of equilibrium, these departures are significant. On the other hand, for the slow and resonant forcing, as well as for the very fast forcing, the log-law could offer a reasonable model of the wind profile (see supplementary material available at https://doi.org/10.1017/jfm.2017.76 for further information about wind velocity animations and the corresponding evolution of the TKE production and dissipation profiles under unsteady forcing).