Global Intermittency and Collapsing Turbulence in the Stratified Planetary Boundary Layer

Abstract

Direct numerical simulation of the turbulent Ekman layer over a smooth wall is used to investigate bulk properties of a planetary boundary layer under stable stratification. Our simplified configuration depends on two non-dimensional parameters: a Richardson number characterizing the stratification and a Reynolds number characterizing the turbulence scale separation. This simplified configuration is sufficient to reproduce global intermittency, a turbulence collapse, and the decoupling of the surface from the outer region of the boundary layer. Global intermittency appears even in the absence of local perturbations at the surface; the only requirement is that large-scale structures several times wider than the boundary-layer height have enough space to develop. Analysis of the mean velocity, turbulence kinetic energy, and external intermittency is used to investigate the large-scale structures and corresponding differences between stably stratified Ekman flow and channel flow. Both configurations show a similar transition to the turbulence collapse, overshoot of turbulence kinetic energy, and spectral properties. Differences in the outer region resulting from the rotation of the system lead, however, to the generation of enstrophy in the non-turbulent patches of the Ekman flow. The coefficient of the stability correction function from Monin–Obukhov similarity theory is estimated as $\beta \approx 5.7$ in agreement with atmospheric observations, theoretical considerations, and results from stably stratified channel flows. Our results demonstrate the applicability of this set-up to atmospheric problems despite the intermediate Reynolds number achieved in our simulations.

1 Introduction

The characteristics of a planetary boundary layer (PBL) crucially depend on its density stratification. In the absence of humidity and advection, the stratification is caused by heating or cooling at the surface—mostly as a consequence of radiative processes. If the surface cools sufficiently, turbulence is often observed to cease partially or even entirely (Mahrt 1999; Wiel et al. 2012). This cessation of turbulence and the associated decoupling of the PBL from the surface impose challenges for mixing formulations in numerical weather prediction (NWP) and climate modelling. Enhanced mixing formulations need often to be used in the mixing parametrizations for boundary layers of NWP and climate models to prevent a decoupling of the atmosphere from the surface (Wiel et al. 2012). Being heuristically formulated and tuned for the performance of NWP models, these enhanced mixing formulations lack a physical basis and cause warm biases at the surface under very cold conditions (Tjernstrom et al. 2005). A better understanding of the underlying dynamics and physical processes, especially in the very stable limit, could hence contribute to alleviate and ultimately overcome these problems of mixing formulations under stable stratification (Mahrt 1999). Conceptual studies of the very stable boundary layer (McNider et al. 1995; Derbyshire 1999; Wiel and Moene 2002) have led to qualitative models of turbulence collapse and global intermittency. In this work, we address fundamental aspects of wall-bounded stably-stratified turbulence and their implications for the stably stratified PBL (the stable boundary layer—SBL).

Often, SBLs are classified into three regimes (Mahrt et al. 1998; Garg et al. 2000; Sun et al. 2012). First, in the weakly stable regime, temperature behaves almost as a passive scalar, and the PBL structure is indistinguishable from the neutral reference: the weakness of the temperature gradient limits the turbulent exchange of heat. Consequently, if stratification is strengthened slightly, the turbulent heat flux increases monotonically as a function of the stratification. Second, in the intermediately stable regime, the turbulent heat flux stagnates if stratification is strengthened; an increased temperature gradient is compensated by a decrease of vertical velocity fluctuations. Third, in the very stable regime, stratification drastically alters the turbulence structure of the SBL to the degree that the weakness of turbulent motion limits the vertical heat exchange: The turbulent heat flux decreases with strengthening stratification. If strong enough, stratification can—locally or globally—lead to the absence of turbulence.

Although the very stable regime is commonly observed in the atmosphere (Mahrt et al. 1998; Ha et al. 2007), there is still a lack of a general framework for the SBL incorporating that very stable regime (Wiel and Moene 2012; Mahrt 2014). Monin–Obukhov similarity theory (MOST, Obukhov 1971) lacks the ability to properly reproduce turbulent fluxes under weak-wind conditions (Ha et al. 2007), i.e. in the very stable regime. From atmospheric observations, it is unclear if stratification can become strong enough to suppress turbulent mixing entirely (Mauritsen and Svensson 2007), and it proves problematic to locally classify a very stable PBL as turbulent or non-turbulent. If turbulence is treated as an on–off process, a runaway cooling at the surface is often seen in PBL models and large-eddy simulations (LES) applied under very stable conditions (Jiménez and Cuxart 2005; Wiel et al. 2012; Huang et al. 2013).

In fact, it is a well-accepted hypothesis that the cessation of turbulence is not an on–off process but rather a complex transition beginning with the local absence of turbulence in an otherwise turbulent boundary layer. Such a local break-down of turbulence is referred to as global intermittency (Mahrt 1999). It has been shown from observations that in a sufficiently stable environment, globally intermittent turbulence can be triggered by a variety of external disturbances including orographic obstacles (Acevedo and Fitzjarrald 2003), solitary and internal gravity waves (Sun et al. 2004) and non-turbulent wind oscillations, such as nocturnal low-level jets (Sun et al. 2012). Whether global intermittency can occur in a SBL without these external triggering mechanisms remains, however, unclear.

A quantitative understanding of strongly stable and globally intermittent turbulence from observations has proved difficult. In particular, accurate flux measurements are hard to obtain with standard methods, and the above mentioned processes often interact: under atmospheric conditions, the entire range of scales, orographic complexity, interaction with the surface, and radiative processes are always present. Thus, it is hard to isolate the signal of a single process as is sometimes necessary for a basic physical understanding. We choose here Ekman flow over a smooth wall—a much simplified configuration. This choice enables a systematic and quantitative study of the intermediately and very stable regimes of turbulence in a simplified and well-defined set-up.

While the study of the weakly stratified limit of simplified SBL configurations is well accomplished by LES (Beare et al. 2006; Huang and Bou-Zeid 2013), the study of very stable cases remains a challenge (Saiki et al. 2000; Jiménez and Cuxart 2005). In particular, the treatment of quasi-laminar patches under very stable conditions is problematic within the conceptual framework of LES. Hence, we consider direct numerical simulation (DNS). Following the pioneering work by Coleman et al. (1990), neutrally stratified Ekman flow has been subject to a number of studies (Coleman 1999; Shingai and Kawamura 2004; Miyashita et al. 2006; Spalart et al. 2008, 2009; Marlatt et al. 2010). Since Coleman et al. (1992) and Shingai and Kawamura (2002), who studied Ekman flow under weak to moderate stratification in rather small domains, however, no studies have come to the authors’ attention.

Instead, it is common practice to use stratified channel flow as a surrogate for the stratified Ekman boundary layer, which is possible due to the analogy between the surface layer of channel flow and that of Ekman flow. In channel flow, oscillations on a period that is large when compared with the large-eddy turnover time, were observed; but no intermittency was found at moderate Reynolds numbers (Nieuwstadt 2005). More recently, global intermittency is simulated in channel flow, and it was proven that too-small domains lead to the formation of artificial flow regimes manifest for instance in the low-frequency oscillation of global statistics (García-Villalba and del Álamo 2011; Flores and Riley 2011). Whereas García-Villalba and del Álamo 2011 use a fixed-temperature boundary condition, Flores and Riley (2011) impose a constant buoyancy flux at the surface. Flores and Riley (2011) observe large-scale intermittency linked to the collapse of turbulence. In contrast to channel flows, Ekman flow has no symmetry in the spanwise direction, which is known to cause large-scale structures in the neutrally stratified limit (Shingai and Kawamura 2004). Whether these structures affect the collapse of turbulence remains unclear. Moreover, Jiménez et al. (2009) found in a non-rotating configuration that the outer flow of boundary layers and channel flows are intrinsically different. Hence, we expect Ekman flow to differ from channel flow as well—we address here the question of how much.

Utilizing our simplified set-up, we complement existing studies of stably stratified flows relevant to the atmospheric boundary layer. We use domains large enough to capture the spatio-temoral structure of a globally intermittent flow, and the stratification is increased to the extreme limit of laminarization of the flow. Our results show that this approach allows for novel insights into the problem of stably stratified boundary layers and enables direct comparison between DNS and observations of the SBL. We address the following research questions:

• (1) Does our simplified set-up reproduce all three regimes of turbulence in terms of a bulk Richardson number as the only control parameter? If so: which are the differentiating bulk properties among these regimes?

• (2) What is the character of the very stable regime? Can we study the turbulence collapse and global intermittency with our simplified set-up – even if the forcing usually made responsible as a trigger for intermittent events is absent?

• (3) Does Ekman flow differ from other configurations used to study the PBL such as channel-flow surrogates? How do large scales in the flow interact with stratification?

2 Formulation

Following Coleman et al. (1992) we solve the governing flow equations in the Boussinesq limit,

$$\begin{aligned} \dfrac{\partial u_i}{\partial t}&= -u_j \dfrac{\partial u_i}{\partial x_j} + \nu \dfrac{\partial ^2u_i}{\partial x_{j}^2} - \dfrac{\partial \pi }{\partial x_i} + f \epsilon _{i3k}(u_k - G \delta _{k1}) + b \delta _{i3}, \end{aligned}$$

(1a)

$$\begin{aligned} \dfrac{\partial b}{\partial t}&= -u_j \dfrac{\partial b}{\partial x_j} + \kappa \dfrac{\partial ^2b}{\partial x_{j}^2},\end{aligned}$$

(1b)

$$\begin{aligned} \dfrac{\partial u_i}{\partial x_i}&= 0, \end{aligned}$$

(1c)

where $u_i$ are the velocity components, $b$ is the buoyancy, $\nu$ is the kinematic viscosity, $f$ is the Coriolis parameter and $\kappa$ is the diffusivity. The buoyancy $b$ is linearly related to the potential temperature $\theta$ as $b(x,y,z)=g\theta (x,y,z)/\Theta _0$ with the acceleration due to gravity $g$ and $\Theta$ is an appropriately chosen background temperature. In the non-turbulent flow above the turbulent boundary layer, the Coriolis force exerted on the geostrophic wind $\mathbf {G}=G\hat{e}_x$ balances the large-scale pressure gradient on an f-plane. We refer to the geostrophic wind direction as streamwise and to the direction opposite to the large-scale pressure gradient as spanwise. The boundary conditions are no-slip at the surface and free-slip at the top. The system is statistically homogeneous in the horizontal directions, and we denote the average of a quantity $\left( a\right)$ along horizontal planes by $\left\langle \left( a\right) \right\rangle$. If $b=0$, i.e. under neutral stratification, the system is also statistically steady and $\overline{\left( a\right) }$ denotes an average over time.

2.1 Scaling of the Neutrally Stratified System

In the neutral case, the boundary-layer dynamics are governed by the quantities $\{G,f,\nu ,\kappa \}$ once turbulence has fully developed and the flow fields have sufficiently de-correlated from its initial conditions. Following previous studies (Coleman et al. 1992; Spalart et al. 2009; Marlatt et al. 2010), we replace in the dimensional analysis the Coriolis parameter $f$ by the laminar Ekman-layer depth $D\equiv \sqrt{2\nu f^{-1}}$. We obtain the Reynolds and Prandtl numbers

$$\begin{aligned} Re=\frac{GD}{\nu }, \quad Pr=\frac{\nu }{\kappa }, \end{aligned}$$

(2a)

and note that $Re\propto \nu ^{-1/2}$. We fix here $Pr\equiv 1$, and hence the steady solution in terms of a statistical description of turbulence of this system is only a function of $Re$.

Once the flow becomes turbulent, the laminar length scale $D$ no longer describes the flow appropriately. Instead, we use the boundary-layer depth scale $\delta \equiv u_\star /f$, where $u_\star$ is the friction velocity defined below. The following parameters characterize the turbulent flow,

$$\begin{aligned} u_{\star }^2 = \left. \nu \frac{\partial \sqrt{\langle u \rangle ^2 + \langle v \rangle ^2}}{\partial z} \right| _{z=0}, \quad \delta = u_\star / f \quad \text {and} \quad Re_\tau = \frac{u_\star \delta }{\nu }. \end{aligned}$$

(2b)

In contrast to channel flows, $u_\star$ in the Ekman layer cannot be known a priori but only a posteriori, and $u_\star$ depends weakly on $Re$ (Spalart 1989). We follow common practice, and study the flow in terms of an inner and an outer layer. In the inner layer, i.e. the surface layer, we choose the wall unit $\nu /u_\star$ and the friction velocity $u_\star$ for normalization; normalized quantities are denoted by a superscript $+$. In the outer layer quantities are normalized by $u_\star$ and $\delta$, and correspondingly normalized quantities are denoted by a superscript $-$. In the outer layer we choose, because of its physical meaning, the inertial period $2\pi f^{-1}$ as the reference time scale.

2.2 Imposing Stratification: Initial and Boundary Conditions

For the stratified cases we use a fully-turbulent, statistically-steady, and neutrally-stratified Ekman flow as initial condition for the velocity fields. Here, we study the problem only for a fixed surface buoyancy; for a discussion of the impact of flux boundary conditions and more complex set-ups we refer to Flores and Riley (2011) and Wiel et al. (2012). We define the buoyancy difference $B_0$ between the surface and the far field (cf. Coleman et al. 1992); this new parameter combines into the Froude number $Fr=G^2/ (B_{0} D)$. $Fr$, however, incorporates the laminar length scale $D$, loosing its relevance in turbulent flow. Stratification can be expressed more appropriately in terms of a global bulk Richardson number,

$$\begin{aligned} Ri_B&\equiv \frac{B_0 \delta _{\mathrm{neutral}}}{G^2}. \end{aligned}$$

(3a)

Due to our choice of a Dirichlet boundary condition,

$$\begin{aligned} b(x,y,z=0,t>0) = 0 \quad \text {and}\quad b(x,y,z=z_\mathrm{top},t>0) = B_0, \end{aligned}$$

(3b)

and—using $\delta _\mathrm{neutral}$ as reference length scale—$Ri_B$ is an inviscid external parameter. We use $Ri_B$ in the following to classify our simulations according to their stratification.

The Obukhov length $L_O$ (Obukhov 1971) is an alternative measure to characterize stratification. In particular, $L_{O}^+$, the ratio of the Obukhov length $L_O$ with the wall unit $\nu /u_\star$,

$$\begin{aligned} L_{O}^+ = L_O \frac{u_\star }{\nu } \equiv \left( -\frac{\theta _0}{g} \frac{u_\star ^3}{\overline{\theta w}|_{z=z_0}}\right) ^+ = \left( \frac{u_\star ^3}{u_\star b_\star }\right) ^+ \quad \text {with}\quad u_\star b_\star = \nu \left. \frac{\partial \langle b\rangle }{\partial z}\right| _{z=0} \end{aligned}$$

(3c)

determines the character of turbulence (Flores and Riley 2011): $L_{O}^+$ measures the biggest possible scale separation in a stratified flow (Flores and Riley 2011), and is therefore an appropriate Reynolds number in a stratified environment. For a smooth wall, this parameter can also be interpreted in terms of the gradient Richardson number

$$\begin{aligned} Ri_G&\equiv \frac{\partial _z B|_{z=0}}{(\partial _z U)^2|_{z=0}} = \nu \frac{b_\star }{u_{\star }^3} = \left( L_{O}^{+}\right) ^{-1}, \end{aligned}$$

(3d)

where we have used the condition $Pr=1$, as defined above. Therefore, $L_{O}^+$ also contains information about the stability character in the near-wall region. Large $L_{O}^+$ implies a small $Ri_G$ and hence turbulence can develop in the lower part of the SBL. In our flow, both $L_O$ and $Ri_G$ are not external parameters to the problem, but they are time-dependent measures describing the evolution of the system. For that reason we use $Ri_B$ as control parameter.

Besides the strength of stratification, the profile of stratification imposed as the initial condition has a significant impact during the initial phase: the initial profile determines the duration of the initial transient. We focus here on the cases where almost the entire stratification concentrates initially within the viscous sub-layer of the flow, mimicking a sudden cooling of the surface. Following Coleman et al. (1992), we choose as the initial condition for the buoyancy

$$\begin{aligned} b(x,y,z,t_0) = 1-\mathrm{erf}\left[ \frac{z/D}{2 a}\right] = 1-\mathrm{erf}\left[ \frac{z}{2a}\left( \frac{\delta _\mathrm{neutral}}{\sqrt{2Re_\tau }} \right) ^{-1}\right] , \end{aligned}$$

(4)

where $b$ is the buoyancy normalized with the surface boundary value $B_0$, $Re_\tau$ is defined in Eq. 2b and $a = 0.15$ is the non-dimensional thickness. Also, $\delta _\mathrm{neutral}$ refers to the value of the neutrally stratified case used for the initialization of the velocity fields. This particular choice of the initial condition through the concentration of the entire buoyancy gradient into the surface layer bears the dilemma that we can have very strong stratification close to the surface even if $Ri_B$, which we expect to control the long-time evolution of the system, is sub-critical. In that case, the initial transient can contain phases wherein the turbulence is shut off through a very efficient cut of production in the buffer layer; this initial phase is most appropriately characterized by $Ri_G$. Thereafter molecular mixing slowly diminishes the buoyancy gradient until the flow becomes unstable again. Although the focus of this work is the long-time evolution of the system, the ratios of these transition time scales to the integral time scale of turbulence, $2\pi f^{-1}$, are relevant in the context of the atmosphere. For instance, they determine whether a fully-developed boundary layer has time to reach it quasi-steady state over the course of a night (or other externally set time scales). In the Appendix we show that our key findings are, at least in the range considered here, independent of the choice of initial condition.

2.3 The Numerical Method

The governing equations are integrated using a high-order finite-difference algorithm on a structured, collocated grid. The time advancement is carried out with a low-storage fourth-order Runge–Kutta scheme (Williamson 1980). Derivatives are computed using a compact scheme described in Lele (1992) such that the accuracy in the interior of the domain is sixth order with overall fourth-order accuracy. We solve the pressure-Poisson equation using a Fourier decomposition along the horizontal directions, which results in a set of second-order differential equations along the vertical coordinate (Mellado and Ansorge 2012). Boundary conditions are no-slip at the wall and free-slip with a Rayleigh damping layer of 10 points and an e-folding time of one inertial period at the top. The grid is stretched moving upwards, and the top boundary is placed at $z\approx 3\delta _{\mathrm{neutral}}$. Adequacy of both horizontal and vertical resolution has been assured for the neutrally stratified simulations using approximately twice the vertical resolution and twice the horizontal resolution for small test simulations at each $Re$ (not discussed here).

2.4 Set-up of Numerical Simulations

A series of simulations is carried out as shown in Table 1 to identify the parameter range of $Ri_B$ in which global intermittency occurs. From Coleman et al. (1992), a configuration of parameters within the weakly stable regime is known, and we use this as a starting point in the weakly stable regime for our case $\mathtt W015S$. With the series of cases W031–S620 the stratification is increased until turbulence collapses.

The impact of both the vertically and horizontally finite domain sizes has been investigated. We find that large domains (at least on the order of $G/f$) are necessary for size-independent fluxes in the outer layer. This is caused by large-scale structures (discussed in Sects. 3 and 5) contributing to the turbulent flux. Hence, the domain was chosen much larger than in previous studies, as interactions between such large-scale coherent structures and the stratification are expected to affect the flow.

Table 1 Set-up of numerical simulations and their parameters

3 The Neutrally Stratified Ekman Layer

This section introduces the neutrally stratified cases. These simulations are used as initial condition and reference for the stratified Ekman layer discussed in Sects. 4 and 5. We show that, despite the moderate Reynolds number $Re$, we can distinguish between the inner and outer layer, the latter being characterized by external intermittency. In the overlap region where the velocity profile is approximately logarithmic, we find both a large-scale structure originating from the outer layer and small-scale hairpin vortices stemming from the buffer layer. We also ascertain the degree of $Re$-independence to better differentiate the effect of stratification by considering three Reynolds numbers $Re=\{500,\,750,\,1000\}$.

3.1 Conventional Statistics

Since the seminal work of Coleman et al. (1990), the wall friction velocity $u_\star$ as well as the turning angle $\alpha$ of the surface stress with respect to the geostrophic velocity are commonly used to compare simulations of neutrally stratified Ekman flows and provide a first estimate of their dependency on $Re$. Figure 1 compares the values of $u_\star$ and $\alpha$ with the semi-empirical theory of Spalart (1989) as well as previous work, and shows that they are within the uncertainty range of the data for our three cases.

Fig. 1

Dependency on $Re$ of the wall friction velocity $u_\star$ (a) and the turning angle $\alpha$ of the surface wind with respect to $\mathbf {G}$ (b) for the neutrally stratified configuration. We compare here also with available data from other groups (Miyashita et al. refers to Miyashita et al. 2006; Coleman et al. and CFS90 refers to Coleman et al. 1990; Coleman 1999; Spalart et al. 2008). The dashed lines show the higher order theory derived in Spalart (1989). The best fit for the available data is obtained with $A=5.9$, $B=0.1$ and $C_5=-30$ following the nomenclature of Coleman et al. (1990). Statistical uncertainty is estimated from the variance in the time series of the respective quantities

Further support for the relatively weak dependency on $Re$ is obtained by the values of the vertically integrated turbulent kinetic energy (TKE, $e$) and the viscous dissipation rate ($\varepsilon$) shown in Table 2 for three Reynolds numbers, where $e$ and $\varepsilon$ are defined as follows,

$$\begin{aligned} e&\equiv \overline{\left\langle \frac{u'_iu'_i}{2}\right\rangle }, \end{aligned}$$

(5a)

$$\begin{aligned} \varepsilon&\equiv \nu \overline{\left\langle \frac{\partial u'_i}{\partial x_j}\left( \frac{\partial u'_i}{\partial x_j} +\frac{\partial u'_j}{\partial x_i}\right) \right\rangle }. \end{aligned}$$

(5b)

The reader is reminded that $\langle \cdot \rangle$ and $\overline{\left( \cdot \right) }$ denote the average over a horizontal plane and the time respectively (note that $\varepsilon$ is the dissipation rate based on the viscous stress tensor (Pope 2000, Eq. 5.160), and not the pseudo- dissipation). We find that the viscous dissipation rate scales independent of $Re$ when normalized with the geostrophic forcing $G$ whereas the energy does so approximately when normalized with $u_{\star }^2$. This further illustrates only slight changes in the organization of the flow while the viscosity varies by a factor of four in this study (cf. Fig. 2).

Fig. 2

Hodograph (a) and streamwise wind-speed profile in semi-logarithmic space (b) for the neutrally stratified cases N500 (red), N750 (blue) and N1000 (orange). In panel a the laminar hodograph is shown as a black dashed line. The levels $z^+=15$ and $z^-=0.12$ are marked by dots in the hodographs to illustrate the increase of scale separation from $Re=500$ to $Re=1{,}000$. In panel b Viscous law of the wall ($u^+=y^+$) and the logarithmic law $u^+=\kappa ^{-1}\,\hbox {ln} ( y^+) + A$ with $\kappa =0.41$ and $A=5.0$ are shown as dashed black lines. As in the hodograph, the level $z^-=0.12$ is marked by a dot in the corresponding wind-speed profile

Table 2 Global statistics as a function of the Reynolds numbers for the neutrally stratified configuration

Velocity hodographs are shown in Fig. 2a. The dependency on $Re$ is also relatively small, in particular the change from $Re=750$ to $Re=1{,}000$ is much smaller than the change from $Re=500$ to $Re=750$.

Attempts have been undertaken to obtain the logarithmic law for Ekman flow and associated constants describing the mean wind-speed profile within the surface layer at Reynolds numbers as high as $Re=2{,}828$ (Spalart et al. 2008, 2009). We find an approximately constant slope of the velocity $U^+$ around $z^+=30$ for the three Reynolds numbers considered here (Fig. 2b); this agrees with Coleman (1999) and Miyashita et al. (2006). In accordance with Tennekes (1973), the height of departure from the logarithmic law for the three cases coincides with the height at which the velocity begins to turn significantly (see markers at $z^-=0.12$ in Fig. 2b). This logarithmic variation supports the analogy with channel flows, which have been studied in great detail.

Beyond the mean profiles, similarity with channel flow in the inner layer is also found in the TKE budget terms as confirmed by Fig. 3a. Irrespective of $Re$, the production of TKE peaks in the buffer layer, around $z^+=12$, and so does the removal by turbulent transport. At the wall, all energy is provided by diffusive downward transport of energy that dissipates locally. The upward transport of TKE away from the production region is caused by turbulent convection. Above $z^+=30$, the TKE budget is dominated by a balance between production and dissipation.

Fig. 3

a Budgets of the evolution equation for $e$ in the inner layer of our simulations N1000 (solid) and N500, N750 (opaque). Circles show data from channel flow [Pope 2000, Fig. 7.34, p. 314]). b Relative contribution in the outer layer. The terms are normalized with their sum of squares at each height

Contributions to the TKE budget in the outer layer, normalized such that at any level the sum of their squares equals one, are shown in Fig. 3b. The change from the production-dominated to the transport- dominated regime, both balanced by viscous dissipation, occurs at $z^-\approx 0.5$. This is about 20 % lower than in a non-rotating boundary layer (Pope 2000, Fig. 7.34), and might hint at the fact that the outer scale $\delta$ has to be adapted by an order-unity constant when quantitatively comparing Ekman flow to other non-rotating flows. Otherwise, the normalized profiles qualitatively agree with those in a non-rotating boundary layer.

3.2 External Intermittency and Conditional Statistics

Despite the qualitative and quantitative agreement with channel flows in many statistics, Ekman flow is not bound by an upper solid wall, i.e. it is an external flow. This causes the co-existence of strongly vortical patches adjacent to much less vortical ones in the outer layer, a property termed external intermittency. External intermittency is widely studied for non-rotating boundary layers since the seminal work by Corrsin and Kistler (1955). They introduced the intermittency function defined as

$$\begin{aligned} \gamma (z) = \overline{\langle H(\omega ^2 - \omega ^2_{0})\rangle } , \end{aligned}$$

(6)

where $\omega$ is the vorticity magnitude, $H$ is the Heaviside function, and $\overline{\langle \cdot \rangle }$ denotes averaging along planes and in time. $\gamma (z)$ is the fraction of the domain at a given height $z$ exceeding a threshold of enstrophy $\omega _{0}^2$, and $\gamma (z)$ is known to be a useful measure of the turbulent area fraction in the outer part of external flows (Kovasznay et al. 1970).

The evaluation of $\gamma (z)$ involves defining $\omega ^2_{0}$ for the turbulent/non-turbulent interface, which is to some degree arbitrary. In many free flows, external intermittency is characterized by very sharp gradients in the enstrophy $\omega ^2$ between the turbulent and non-turbulent patches, and this feature deems $\gamma (z)$ quite insensitive to the choice of $\omega _0$, at least in a certain range (Kovasznay et al. 1970; Bisset et al. 2002; Mellado et al. 2009). Here, we choose $\omega _0=\omega _{rms}(\delta )$, the root-mean-square value (r.m.s.) of vorticity $\sqrt{\overline{\langle \omega '^2 \rangle }}$ at $z=\delta$ as a reference vorticity for the distinction between turbulent and non-turbulent regions for several reasons. First, this level—according to classical definitions of the boundary-layer height such as $\delta _{95}$ (Table 2)—is well outside the part that is generally considered turbulent. Second, we find $\overline{\left\langle u_i u_i\right\rangle } \propto z^{-4}$ for $0.75\lesssim z^-\lesssim 2$ (not shown), which is a signature of potential flow above a turbulent boundary layer (Phillips 1955). Third, the resulting profile $\gamma (z)$ (Fig. 4) is similar to that found in non-rotating boundary layers (Kovasznay et al. 1970).

Fig. 4

a Intermittency factor versus height for the Reynolds number $Re=500$ (red) and $Re=1{,}000$ (orange) where the threshold $\omega _0$ expressed in terms of the vorticity r.m.s. at $z=\delta$ is varied by a factor of 4. b Turbulent kinetic energy conditioned to turbulent and non-turbulent patches of the flow for $\omega _{0}= \omega _\mathrm{rms}(\delta )$

For Ekman flow, it is readily seen from the enstrophy equation that vortex tilting of planetary vorticity $f\hat{e}_z$ constitutes a source $f \partial u/\partial z$ of streamwise vorticity $\omega _x$ (and hence a source $2f\omega _x \partial u/\partial z$ of streamwise enstrophy $\omega _{x}^2$). This source is absent in a non-rotating reference frame. Hence, small gradients in the mean velocity profile $\partial \langle u \rangle / \partial z$ will, even in a purely non-turbulent flow, generate mean vorticity that grows until balanced by dissipation. In comparison to the non-rotating flows mentioned above, this enstrophy generation leads to a larger variation in $\gamma (z)$ if the threshold is doubled and halved with respect to $\omega _0$ (Fig. 4a). The $Re$-independency of $\langle u(z)\rangle$ in the outer layer suggests that this mechanism is independent of the Reynolds number since the term $\langle \omega _x\rangle f \partial \langle u\rangle /\partial z$ scales inviscidly. This $Re$-independency is corroborated by the small sensitivity of the profile to $Re$ seen in Fig. 4. We conclude that this vortex tilting, irrespective of $Re$, is a fundamental mechanism in Ekman flow that renders the outer layer different from non-rotating external flows.

The definition of a discriminator between turbulent and non-turbulent regions of the flow allows for the use of conditional statistics. It enables a separation of turbulent from non-turbulent contributions to bulk quantities of the flow that is useful for a process-oriented study of the flow based on fundamental principles. Such a separation is important, in particular in the outer layer where the variation of mean properties between turbulent and non-turbulent patches can contribute significantly to the variances (Pope 2000, Eq. 5.306); in terms of the fluctuation kinetic energy $e$, this decomposition reads as

$$\begin{aligned} e = \gamma e_T + (1-\gamma ) e_N + \gamma (1-\gamma )\left( \langle u_T\rangle -\langle u_N\rangle \right) ^2, \end{aligned}$$

(7)

where the subscripts $T$ and $N$ stand for the average of the respective quantity over the turbulent and non-turbulent regions—as determined by a point-wise application of the intermittency conditioning–, respectively. In our flow we find $\gamma (1-\gamma )\approx 0$ below $z^-\approx 0.5$ and above $z^-\approx 1.5$ (Fig. 4a): Below $z^-\approx 0.5$, the turbulent fraction is very close to one, whereas above $z^-\approx 1.5$, the total variance is largely explained by the non-turbulent fluctuations (Fig. 4b). In between those extremes, $\gamma (1-\gamma ) >0$. Hence, both the turbulent as well as the non-turbulent variances, and thus also their cross term $\gamma (1- \gamma )(\langle u_T\rangle - \langle u_N\rangle )^2$ contribute to the total variance.

3.3 Flow Visualizations

The vertical structure suggested by $\gamma (z)$ is consistent with a visual inspection of the flow enstrophy fields (Fig. 5). The strongly vortical regions adjacent to the surface (Fig. 5a–b) are typical of wall-bounded flows and indicate the level of the buffer layer. In the buffer layer, vorticity is mainly associated with so-called surface streaks (Fig. 5c). Ejections of turbulent fluid are seen as excursions of the white colours into higher levels of the boundary layer (Fig. 5a, b).

Fig. 5

a, b Streamwise–vertical cross-sections showing the magnitude of the gradient of concentration of a passive scalar originating from the surface. Grey scale varies from $10^3 \varLambda _{RO}^{-1}$ (white) to $10^{-1}\varLambda _{RO}^{-1}$ (black) where $\varLambda _{RO}=G/f$. Block grey-shading in a indicates the region shown b. c–g Horizontal cross-sections. Wind magnitude in the buffer layer ($z^+\!\approx \! 15$ c), the upper part of the logarithmic layer ($z^+\!\approx \! 100$, $z^-\!\approx \! 0.11$ d), and the outer layer ($z^-\approx 0.75$, e). c Only the subset marked by white shading in d and e. f, g Plot the scalar gradient at $z^+\approx 40$ in the whole domain f and for the grey shaded square from f illustrating the hairpin vortices g. Shown in this figure is case N1000L, which has been carried out in the computational domain rotated by approximately $-\alpha$: The horizontal planes in frames c–g are rotated such that the geostrophic velocity vector is pointing from left to right as shown in the sketch

The change of organization in the turbulent flow when moving upwards is illustrated by horizontal cross-sections of wind magnitude (Fig. 5c–e): In the buffer layer (Fig. 5c), the flow is dominated by surface streaks aligned with the mean wind at that level, which is anti-parallel to the force exerted on the fluid by surface shear stress $\tau _\mathrm{wall}$. In the fully turbulent part of the outer layer (Fig. 5d), the turbulence is modulated at large scale that is rotated by about $20^\circ$–$30^\circ$ clockwise with respect to the geostrophic wind. At this height, the small-scale structure appears as noise. At higher levels (Fig. 5e), the boundary layer is externally intermittent, because turbulence at those levels is mainly provided by strong ejections from lower levels happening only sporadically. Such generated turbulent structures in the outer layer of the flow are long-lived because of their relatively large extent and the weakness of turbulent dissipation at these large scales.

Horizontal planes in the quasi-logarithmic layer of the flow (at $z^+\approx 40$, Fig. 5f–g) show—consistent with the intermittency function $\gamma (z^+=40)\approx 1$—that the field is homogeneously turbulent, and that the dominating small-scale structures are hairpin vortices typical of the logarithmic layer (Adrian 2007). Their intensity is modulated at a large scale that is rotated about $20^\circ$–$30^\circ$ clockwise with respect to the geostrophic wind vector, similarly to the large-scale organization observed before in Fig. 5d. This large-scale organization is typical of wall-bounded flows (Marusic et al. 2010; Adrian 2007), and there remains considerable controversy about the role of these large-scale structures in the inner layer (Jiménez 2012). In our case, they have a clear organization that can be attributed to some large-scale instability inherent to the flow (Barnard 2001). We consider the existence of such large-scale structures a fundamental property of turbulent Ekman flow and expect that they are crucial when the flow is exposed to stable stratification, as discussed in Sect. 5.

4 Turbulence Regimes and Stability

When the flow is exposed to a suddenly cooled surface (Sect. 2.2), there is an initial interval of time in which the buoyancy mainly mixes by viscous diffusion, and a thin strongly stratified layer develops at the surface. This, when compared to the inertial period, short initial transient in our simulations ends at $t^-\approx 0.2$, and this transient is followed by a much slower recovery (Fig. 6a, b). Despite our fixed-temperature boundary condition, the recovery is similar to the evolution of a nighttime boundary layer with residual turbulence in the outer layer (Coleman et al. 1992). We focus on this slow recovery and, in this section, classify our simulations as weakly, intermediately, or very stable following the classification introduced in Sect. 1. Note that, throughout this slow evolution, turbulence in the entire column does not adapt immediately to a change in surface friction, and both $u_\star$ and $\delta$ evolve on a time scale $\delta /u_\star =1/f$; for this reason, we use $u_{\star ,\mathrm{neutral}}$ and $\delta _{\mathrm{neutral}}$ instead as reference scales to normalize the results and compare with the neutral case.

Fig. 6

Panel a shows the temporal Temporal evolution of vertically integrated turbulent kinetic energy $E(t)$ (solid) and averaged wall friction velocity $u_\star (t)$ (dashed) normalized with the corresponding neutral reference $E_\mathrm{neutral}$ and $u_{\star ,\mathrm{neutral}}$, respectively. b Same as a, but for the streamwise enstrophy $\langle \varOmega _{x}^2\rangle (t)$. Panel c plots Vertically integrated terms of the TKE budget equation at $t^-=tf/2\pi =0.5$ for cases W015S, W030S, I150, S310S and S620 normalized with the shear production rate of the neutral reference

4.1 Classification

The time evolution of the vertically integrated TKE (Fig. 6a) suggests the following classification of our simulations:

• (1) Weakly stable: integrated TKE changes slightly (10–20 %, cases W015S, W031S) with respect to the neutral configuration

• (2) Intermediately stable: integrated TKE significantly (50 %) decreases and subsequently recovers (case I150, S310S)

• (3) Very stable: integrated TKE is diminished nearly entirely, and subsequently recovers (case S620).

This classification is supported by the time series of integrated enstrophy as well as the vertically-integrated budget of TKE at $t^-=0.5$ (Fig. 6b, c): in particular, the buoyancy flux is a maximum among all cases in case I150 supporting its identification as intermediately stable (Mahrt et al. 1998). Both the shear production (blue bars) and the buoyancy flux (black bars) change drastically when the most stable case S620 is considered. In this case the terms in the TKE budget as well as the TKE itself reduce to ${}\approx 5~\%$ of the neutral reference. This reduction of order one in both the turbulence production and buoyancy flux with respect to the neutral reference illustrates that the buoyancy flux is limited by the absence or weakness of turbulent motion and not by the pure strength of buoyancy destruction $\int \langle bw\rangle \mathrm{d} z$.

After this initial decrease or even breakdown of turbulence, such as that simulated in case S620, the turbulence intensity recovers because the buoyancy difference across the boundary layer is fixed (Sect. 2.2). Hence, the stratification close to the surface decreases to compensate for a deepening of the stratified layer. Concomitantly, the fluctuation kinetic energy reaches and exceeds that of neutral stratification; such a recovery happens in all our cases (Fig. 6a). The magnitude of this recovery is even larger if the energy is normalized with the instantaneous friction velocity $u_\star (t)$ as $u_\star (t)$ decreases when the flow is exposed to stratification (Fig. 6a). The reasons for this recovery are complex, and for an explanation we refer the reader to Sect. 5. None of the cases reaches an equilibrium on time scales of one or two inertial periods. In the Appendix we study the associated time scales, and we conclude that, even under weak stratification, no equilibrium of the entire turbulent boundary layer is reached over the course of a night.

A possible origin of this increase in fluctuation kinetic energy with respect to the non-stratified configuration is non-vortical or, when compared to turbulent eddies, weakly vortical large-scale modes—possibly gravity waves—in the flow. This is confirmed by the fact that, unlike the fluctuation kinetic energy, the integrated streamwise vorticity r.m.s. remains below the neutral reference (Fig. 6b). We investigate this and the organization of the flow further in Sect. 5 by means of the intermittency function, conditional statistics, and spectra. The remainder of this section is devoted to a more detailed description of the turbulence regimes that we have identified.

4.2 The Weakly Stable Regime

The boundary layer forming in the weakly stable regime ($Ri_B \lesssim 0.05$; W015S, W031S) is very similar to that found in the neutral reference (Fig. 8a, d; Sect. 3). Turbulent mixing efficiently weakens the stratification and a quasi-neutral weakly-stable boundary layer forms, which after $tf/2\pi \approx 1$ enters into a relatively small amplitude oscillation. As expected and found elsewhere (Monin 1970; Ha et al. 2007; Sun et al. 2012), the weakly stable boundary layer is well described when considered as a perturbation of the neutrally stratified one. TKE is altered most strongly in the outer layer (not shown) as also found by Coleman et al. (1992) and García-Villalba and del Álamo (2011), and the hodograph is barely distinguishable from that of the corresponding neutral case N500 (Fig. 7c, red line).

Fig. 7

Vertical profiles of $e$ (a) and normalized buoyancy frequency $N/f$ (b) blue case I150, orange case S620. d The in-plane Reynolds stress (solid: $t^-=1$; dashed: $t^-=2.0$; dash-dotted: $t^-=3.0$). Thin solid lines show the initial condition for the respective case. Panel c shows Hodographs after one inertial period ($t^-=1$) and those from the neutral cases repeated from Fig. 2a

Fig. 8

Horizontal cross-sections of enstrophy in a sub-domain spanning $\approx (10 \times 10)\delta _\mathrm{neutral}^2$ at $tf/2\pi \approx 1$. a, d $Ri_B=0.015$, case W015S; b, e $Ri_B=0.15$, case I150; c, f $Ri_B=0.62$, case S620. a–c The inner layer $z^+\approx 15$, d–f the outer layer around $z^-=0.56$

4.3 The Intermediately Stable Regime

In the intermediately stable regime at $Ri_B=0.15$, TKE is reduced by $\approx 50~\%$ throughout the initial transient, and the integrated buoyancy flux at $t^-=1$ is the maximum of all simulations carried out within this study (Fig. 6c). The increase in the integrated buoyancy flux is one order of magnitude smaller than the reduction in TKE (Fig. 6a) and shear production (Fig. 6c) with respect to the neutral reference. This illustrates that the main impact of buoyancy on the flow is not the direct destruction of TKE but a decrease in the shear-induced production, in particular of the stress $\langle u'w'\rangle$ (Jacobitz et al. 1997, p. 243). Profiles of shear production (not shown) indeed confirm this explanation, in agreement with Jacobitz et al. ’s study of stably stratified shear flow.

In contrast to the simulations attributed to the weakly stable regime, the hodograph in the intermediately stable regime ($Ri_B=0.15$, case I150) departs significantly from the neutral reference case. It lies in between the hodographs from the neutrally stratified case and a laminar one (Fig. 7c).

After an initial decay, the fluctuation kinetic energy recovers slowly on a time scale of a few inertial periods (blue curve in Fig. 7a). If expressed in terms of $f^{-1}$, the time scale of this slow oscillation matches the time scale for recovery observed by Nieuwstadt (2005) in a stably stratified channel flow. Turbulence intensity recovers across the entire boundary layer, and concomitantly the depth of the stratified layer increases (sequence of blue lines in Fig. 7b). This increase in depth of the stratified layer is compensated by weakening stratification in the surface layer ($z^-\lesssim 0.1$). Eventually, during this recovery, the fluctuation kinetic energy increases beyond the neutral reference both above $z^-\approx 0.5$ and in the production region (Fig. 6a) as also observed by Nieuwstadt (2005).

In agreement with recent work on channel flow (Flores and Riley 2011) and in contrast to the findings of Nieuwstadt (2005), the simulated boundary layer is globally intermittent. A local break-down of turbulence is evident from Figs. 8 and 9 showing quasi- laminar patches in a turbulent environment. These quasi- laminar patches extend through the entire vertical fluid column in an otherwise turbulent flow, and have a distinct alignment clockwise to the geostrophic flow. Hence, we identify this state with that of global intermittency in the sense of Mahrt (1999). We note that our $Ri_B$ defined in terms of $\delta _{\mathrm{neutral}}$ as an external control parameter is smaller than the Richardson number defined in terms of the depth $\delta (t)$ of the SBL ($\delta (t)\approx 0.5\delta _{\mathrm{neutral}}$ at $tf/2\pi =1$; see Fig. 7b). Hence, the occurrence of global intermittency in this particular case at $Ri_B=0.15$ agrees with the observation that global intermittency often occurs if $\delta B_0/G^2\approx 1$.

Fig. 9

Flow visualizations at $tf/2\pi \approx 2$. a Enstrophy iso-surface of $\omega ^2\approx \left( 20 \omega _\mathrm{rms}(\delta )\right) ^2$ coloured by horizontal wind speed in the range $0.4<\sqrt{u^2+v^2}/G<1.15$. b Wall shear stress (colour coded from black (low) to white (high). a, b Gives a perspective looking from above. c Iso-surface of non-dimensional buoyancy $b/B_0=0.8$ ($b(z=0)=0,\quad b(z=z_{top})=B_0$) coloured by streamwise wind speed. Low wind speeds (blue) tend to occur close to the surface whereas high wind speeds occur further up in the flow. d Schematic of flow organization as illustrated in a–c of this figure: streaks aligned with the surface shear stress are represented by the red dashed line. The orientation of the large-scale structures in the outer layer is shown by the blue line. (The red lines in a and b indicated the location of the vertical cross-section presented in Fig. 11; the white line indicates the alignment of the large-scale structure discussed in the main text)

4.4 The Very Stable Regime

Under very strong stability (case S620), the turbulence dies out nearly completely since the production region is eliminated. The hodograph (Fig. 7a) is close to that of the corresponding laminar Ekman flow. In fact, the eddy diffusivity estimated from the laminar fit to the velocity profiles (not shown) is $1.01\nu$. This re-laminarization in the inner layer is seen in Fig. 8: the turbulence with relatively high enstrophy magnitudes in panel (a) is replaced in panel (c) by large-scale roll-like structures aligned parallel to the wall-shear stress, i.e. is rotated by $45^\circ$ counter-clockwise with respect to the geostrophic wind vector. This initial re-laminarization is followed by a recovery of turbulence as seen in the time series of TKE and enstrophy (Fig. 6a, b). The recovery of turbulence is similar to that observed in the intermediately stable case. This recovery, however, takes longer, and while the enstrophy levels off around $60~\%$ of the neutral value, the normalized TKE grows beyond one.

At the beginning of the recovery of TKE, the maximum of TKE associated with the peak shear production in the buffer layer is eroded (Fig. 7a), that is, the production region of turbulent stress is eliminated. Around $z^-=0.25$, the turbulence intensity is reduced even more than at the peak of production; this illustrates the absence of vertical turbulent exchange across the buffer layer and a decoupling of the flow inside this surface layer from upper layers of the flow. Above the decoupled surface layer ($z^-\gtrsim 0.5$), turbulence is affected less strongly by stratification and decays slowly from its fully-turbulent initial state between $t^-\approx 1$ and $t^-\approx 2$. Such slowly-decaying residual turbulence is common to nighttime boundary layers cooled from below (Stull 1988), which illustrates the appropriateness and the relevance of these simulations for the study of such cases.

The decoupling of the outer layer from the surface layer is an important consequence of very strong stratification, and there has been debate on whether the decoupling produced by boundary-layer schemes in NWP models (Derbyshire 1999; Acevedo et al. 2012) and LES (Saiki et al. 2000; Jiménez and Cuxart 2005) is an artifact of the turbulence subgrid model. From our data, we conclude that, at bulk Richardson numbers of order one ($Ri_B=0.62$ for this particular case), a decoupling is possible—at least for an intermediate Reynolds number. This is in accordance with Wiel and Moene (2012). Our estimate, in contrast to estimates from NWP or LES, is not subject to uncertainties in subgrid schemes, but, similar to stratified shear flow (Jacobitz et al. 1997), the particular value of a critical Richardson number for decoupling might depend on $Re$.

5 Flow Organization

By comparison between the neutral case N500 and the intermediately stable case I150 at $Ri_B=0.15$, we show in this section that the recovery of turbulence in the intermediate cases is accompanied by a large-scale organization of the flow. This organization efficiently couples the outer and inner layers of the flow and hence is the main mechanism governing the spatio-temporal structure of global intermittency in the intermediately stable regime. We present a measure to quantify this global intermittency. For that purpose, we extend the classical concept of external intermittency, that is, the alternation of turbulent and non-turbulent patches of fluid in the outer layer: If there are laminar patches of fluid extending from the outer layer down to the surface layer, we identify the flow as globally intermittent.

5.1 Spatial Variability

We have seen in Fig. 5 how the spatial organization of the outer layer leaves a footprint in the quasi-logarithmic layer even in the neutral case. A conspicuous manifestation of this organization is evident in Fig. 9 for the intermediately stable case I150: strongly convoluted parts of an enstrophy iso-surface penetrate deep into the outer layer of the flow where the wind velocity reaches values comparable to the geostrophic velocity. These bulges coexist with much smoother patches of low velocities marking regions very close to the surface (Fig. 9c), where the enstrophy is almost entirely caused by the background shear. Signatures of streaks are also present in these smoother patches; streaks, however, seem to be pushed downwards, which prevents ejections into the boundary layer and locally limits buoyancy and momentum exchange.

The angle along which the laminar and turbulent structures are evidently oriented is estimated to be around $23^\circ$ clockwise with respect to the geostrophic wind (white line in Figs. 9, 11). This is the same orientation that is observed for the large-scale outer-layer structures under neutral stratification (Sect. 3; Fig. 5). The wavelength $\lambda$ of these smooth patches has been extracted from visualizations and was measured as $2.5 < \lambda /\delta < 4.2$ with $\delta$ the boundary-layer depth scale as introduced in Eq. 2b. In fact, also the spectrum of the vertical velocity $w$ averaged over the layer $0.62 < z^- < 0.65$ (Fig. 10) shows three isolated local maxima of the energy density corresponding to wavelengths of $\lambda /\delta \approx \{2.68, 4.14, 5.26\}$.

Fig. 10

Turbulent energy spectrum of the neutral case (N500; contours) and the intermediate case (I150; shading) averaged between $0.62 < z^- < 0.65$ around $t^-=1$. Wavelengths $\lambda ^-$ are normalized with $\delta =u_\star /f$ for the neutral case, and wavelengths $\lambda ^{LM}$ are normalized with $L_O=u_{\star }^2 b_{\star }^{-1}$ for the stable case I150 (shading). Following García-Villalba and del Álamo (2011), we use $\lambda ^{LM}:=\lambda /L_O$ where $L_O=1.07\delta$ at $t^-=1$

Fig. 11

Vertical cross-section along the red line shown in Fig. 9a, b at $tf/2\pi \approx 2$: a logarithm of enstrophy colour-coded linearly from black to white in the range $-8 < \zeta < 13$; b streamwise velocity $U/G$ colour-coded from black to white in the range $0.12 < u/G < 1.11$; c buoyancy colour-coded from black to white in the range $0.96 > b/B_0 > 0.04$; d vertical velocity colour-coded from red to blue in the range $-0.43 < w/G < 0.43$. The vertical bars indicate the length scales $\delta _\mathrm{neutral}, \delta _{95}$ and $0.1\delta _\mathrm{neutral}$

The organization of the globally intermittent flow is qualitatively summarized in Fig. 9d: the surface streaks (red, dashed) are parallel to the wall shear stress; their separation distance is much smaller than the boundary-layer depth (Fig. 8). Even adjacent to the wall, the flow is not fully turbulent everywhere but the turbulence is patchy. There are laminar patches extending through the entire vertical column of the flow. The fully turbulent and quasi-laminar patches in the inner and outer layer are organized along lines oriented as indicated by the solid blue line in Fig. 9d. The angle with the mean flow is ${} \approx {23}^\circ$ clockwise, which is similar to the neutral case (Sect. 3). That is, under stable stratification, the orientation of the large-scale structures in the outer layer sets the spatial organization of the flow even in the inner layer by determining the patterns of global intermittency. The wavelength of these large-scale structures, which are already present under neutral conditions, is of the order of several boundary-layer depths (2–5$\delta$, a fraction $1/10\,\hbox {to}\,1/4$ of the Rossby radius); these structures need to be resolved for a realistic representation of this complex flow.

5.2 External and Global Intermittency

A quantitative measure of external intermittency in the flow is provided by the intermittency factor $\gamma$ introduced in Sect. 3.2. Compared to the neutral case (Fig. 12, grey line), the intermediate case (blue line) has a much lower intermittency factor in the outer layer, which indicates the absence of turbulence in this region. This absence occurs between $z^-\approx 0.5$ and $z^-\approx 1.0$, that is, precisely at those levels where the TKE increases beyond the neutral reference (Fig. 7a) while the integrated streamwise enstrophy r.m.s. (Fig. 6b) does not. This apparent contradiction is resolved if the flow field is inspected visually: in the outer layer of the neutral case, the turbulent field is characterized by a variety of small vortices with rather small vertical velocities (Figs. 5, 8d). In the intermediate case, after one inertial period, these structures are replaced by the large-scale structures discussed in Sect. 5.1 above. We show streamwise-vertical intersects approximately perpendicular to the large-scale structure in the flow (red and white lines in Fig. 9) in Fig. 11. This large-scale structure carries TKE by means of the comparatively large values of vertical velocity: there are large regions with vertical velocities as large as ${} \approx {\pm }\,\,0.1G$ (Fig. 11d). Regions of large vertical velocities are connected by corresponding positive and negative streamwise (Fig. 11b) and spanwise (not shown) anomalies in between them, which indicates a roll-like structure approximately aligned with the large-scale organization discussed above. These large-scale structures are responsible for the large increase of TKE (about 25 %, cf. Fig. 6a) beyond the neutral reference and thus resolve the above noted apparent contradiction.

The intermittency factor, i.e. the area fraction occupied by turbulent patches, decreases even close to the surface ($0.1<z^-<0.2$; see inset of Fig. 12a) with respect to the neutral reference. Given the visually large portion of quasi-laminar patches (Fig. 8), the change of only ${\approx }{5}~\%$ in $\gamma$ appears, however, quite small. This problem with the recognition of quasi-laminar patches in the inner layer is due to the high background enstrophy. This enstrophy stems from the mean velocity gradient contributing a large portion of the mean enstrophy in the inner layer of the flow. This limits the change in $\gamma$. In fact, this method, introduced to measure global intermittency in the outer layer, only detects laminar patches of fluid originating from the outer layer and penetrating deep into the inner layer of the flow. In contrast, the method cannot detect non-turbulent regions where turbulence has collapsed locally in the near-wall region since those are characterized by high mean enstrophy. Nonetheless, we consider the intermittency factor useful in this context as it properly represents the state of the flow: a complete depletion of turbulence in the upper part of the outer layer as well as globally intermittent turbulence in the lower part of the outer layer and in the inner surface layer.

Fig. 12

a Intermittency factor profiles for the intermediately stable case I150 (blue $t^-=\{1,2,3\}$) and the very stable case S620 (orange $t^-=\{1,1.75\}$). Inset shows the same focused on the lower part of the boundary layer. b TKE conditioned to turbulent and non-turbulent patches as in Fig. 4 (grey) and for $t^-=3$ of case I150 (blue)

The TKE conditioned to turbulent and non-turbulent patches of the flow at $Ri_B=0.15$ (case I150) is plotted in Fig. 12b. Beyond $z^-\approx 0.7$ the contribution of the non-turbulent region in the stable case increases with respect to the neutral reference, and this increase partly explains the overshoot in the integrated TKE. We also note that, according to our partitioning, the turbulent contribution to the fluctuation kinetic energy $e$ does not only increase above $z^-\approx 0.7$ but also in the lower parts of the outer layer, around $z^-\approx 0.2$. This is likely due to the problems of the intermittency function with the recognition of global intermittency close to the surface as described in the previous paragraph: contributions of the cross term in Eq. 7 are aliased onto the turbulent partition. Further work on conditioning methods to study global intermittency is necessary.

6 Discussion

6.1 Relation to Monin–Obukhov Similarity Theory

Agreement of our numerical experiments with Monin–Obukhov Similarity theory (MOST, Obukhov 1971) supports their relevance for atmospheric conditions, despite the difference in $Re$. A common formulation for the gradient of velocity in the stratified surface layer is

$$\begin{aligned} \frac{\kappa z^+}{u_\star }\frac{\partial U}{\partial z^+}&= \varPhi _M. \end{aligned}$$

(8a)

Under the assumption of a stability correction function of the form $\varPhi _M(\zeta )\equiv 1+\beta \zeta$, with $\zeta \equiv z/L_O$, we obtain

$$\begin{aligned} \int \limits _{z_0}^{z} \varPhi _M \frac{\partial {z^{'+}}}{z^{'+}} = \int \limits _{z_0}^z \mathrm{d} \ln (z'^+) + \varPhi _M-1 - \beta \frac{z_0}{L_O}. \end{aligned}$$

(8b)

The term $\beta z_0/L_O$ is the stability correction at the height $z= z_0$, the lower end of the surface layer and in the following it is assumed $\varPhi _M(z_0/L_O)=0$. Hence, using the common assumption $z_{0}^+=1$,

$$\begin{aligned} \varPhi _M - 1&= \kappa \left( U^+(z) - \kappa ^{-1}\ln {z^+}\right) = \beta \zeta . \end{aligned}$$

(8c)

Eq. 8a implies a logarithmic profile $U^+=\kappa ^{-1}\ln {z^+}$ if $\varPhi _M=1$. Instead, we use the profile of the corresponding neutral reference case, and estimate the stability correction as

$$\begin{aligned} \varPhi _M-1= \kappa \left( U^+(z) - U_{\mathrm{neutral}}^+(z)\right) \end{aligned}$$

(8d)

with the von-Kármán constant $\kappa =0.4$. Under weak stratification, we find good agreement between MOST and the measured stability correction from our data (red circles in Fig. 13). If the stratification is increased to $\zeta >0.05$, scatter increases and errors are on the order of $10~\%$ for $Ri_B=0.15$ and of the order of $20~\%$ for $Ri_B=0.3$: the limit of applicability for MOST is reached. In the very stable case where we have seen that turbulence has ceased nearly entirely, the measured and predicted stability corrections differ by about 100 % indicating that MOST is inapplicable under very strong stratification. Indeed, we have shown above that, in the latter case, the profile is very well approximated by a laminar one. A least-square fit for the coefficient $\beta$ in $\varPhi _M$ gives $\varPhi _M-1=5.68\zeta - 0.04$. We explain the small offset $0.04$ of $\varPhi _M-1$, which has a theoretical value of $\varPhi _M(0)\approx 0$, by the neglect of the lower boundary condition $\beta z_0/L_O$. The fit explains $98.5~\%$ of the variance considering data from cases with $Ri_B \le 0.31$. This agrees well with data obtained from atmospheric measurements ($\beta =5.3$; Högström 1996) and channel-flow DNS ($\beta =4.5$; Wiel et al. 2008).

Fig. 13

Stability correction function $\varPhi _M -1 := \kappa \left( U^+(z) - U_{\mathrm{neutral}}^+(z)\right)$ evaluated from three stably stratified cases with respect to the neutral reference profiles. The Obukhov length is calculated from the fluxes at the surface. The dashed black line shows the linear fit $(\varPhi _M-1)=5.7\zeta$

6.2 Global Intermittency

In the light of our results, the occurrence of global intermittency might be explained as follows: once turbulence cannot be fully sustained, the turbulence decays and the relative size of the non-turbulent region is determined by a bulk property of the system, like a bulk Richardson number in the case studied here. We thus proved that global intermittency can arise from a global constraint on the flow, for instance the exceedance of the maximum sustainable heat flux (Wiel et al. 2012; Wiel and Moene 2012). Local perturbations, such as surface heterogeneities—or large-scale dynamics as in our case—, simply determine the spatio-temporal distribution of global intermittency (Sun et al. 2012, 2004; Acevedo and Fitzjarrald 2003), but they are not needed as a trigger.

While previous work mainly considers local or external mechanisms as a trigger for global intermittency in the SBL, we show that this process is intrinsic to a stratified atmosphere. This has important implications for turbulence models applied under very stable stratification: global intermittency effects could be incorporated into turbulence closures for Reynolds-averaged NWP models once its dependency on $Ri$ and $Re$ is quantitatively understood—without depending strongly on the local details of the flow such as surface heterogeneities.

Owing to the absence of heterogeneities in our numerical set-up, the spatio-temporal pattern of global intermittency close to the surface is caused by a large-scale structure in the outer layer of the flow. Indeed, Chung and Matheou (2012) and Brethouwer et al. (2012) observe a very similar phenomenon in homogeneously stratified sheared turbulence respectively rotating Couette flow, in both of which no local or coherent external perturbations are present. This occurrence of global intermittency and its linkage to large-scale forcing from the outer layer is consistent with Mauritsen and Svensson (2007) who, based on observational data, suggest internal gravity waves as a cause of non-zero turbulent fluxes under very stable stratification.

The agreement with observations in terms of $Ri_B$ and the occurrence of global intermittency suggests that this mechanism is relevant at the atmospheric scale. If our simplified set-up is considered in the phase space spanned by its non-dimensional parameters $Ri_B$ and $Re$, there are two ways in which a transition from turbulent to laminar can happen: either through stronger stratification, that is changing only the parameter $Ri_B$, or via a decrease in the mean shear, affecting both $Ri_B$ and $Re$. Once in the fully turbulent regime, it is unlikely that the fundamental character of this transition changes. Such a change would in fact require turbulence in the SBL to be caused by a different instability than in our simulations. The quantitative dependency of this mechanism on $Re$ has, however, to be elucidated in future work. We expect certain properties of this transitions such as a critical Richardson number to depend on $Re$ similarly to $u_\star$ and $\alpha$ (Sect. 3) and as observed in stably stratified shear and channel flow (Jacobitz et al. 1997; Flores and Riley 2011). Nonetheless, we have shown that this simplified set-up is suited to study dynamics of the stable and very stable stratified PBL. This analogy over a cascade of complexity—ranging from the canonical flow problem of stable channel flow via rotating Cuoette and stable Ekman flow to an atmospheric boundary layer—encourages further investigation of the fundamental aspects of stably stratified turbulence in rotating reference frames.

7 Conclusions

We have defined a framework to investigate the SBL using direct numerical simulation of turbulent Ekman flow over a smooth surface with a fixed temperature. This set-up is described in terms of a Reynolds and a Richardson number solely.

In the neutral limit, the flow depends only weakly on the Reynolds number as we demonstrate by a comprehensive flow description including hodographs, the integral value of the turbulence dissipation rate as well as vertical profiles of velocity, TKE, and external intermittency. We thereby show that the analogy between the surface layer in Ekman flow and channel-flow also applies to the budget of TKE, and not only to the logarithmic law of the mean velocity. Through a conditioning of the TKE to the turbulent and non-turbulent patches, we demonstrate that the alternation between those contributes significantly to the velocity variance in the outer layer of the flow.

Considering our stably stratified cases, we conclude the following:

• (1) Stably stratified Ekman flow is suited to study aspects of the SBL. The regimes of stratified turbulence are reproduced varying a single parameter: the bulk Richardson number. Characteristics of these turbulence regimes, such as hodographs and TKE profiles, compare well with atmospheric observations and theoretical considerations. We estimate the Monin–Obukhov stability correction for stable stratification as $\varPhi _M-1 \approx 5.7z/L_O$ in agreement with data from channel-flow and atmospheric observations.

• (2) The analogy of the surface layer with channel flow holds beyond qualitative aspects. The recovery and overshoot of integrated TKE as well as the laminar patches observed in our simulations are similar to their equivalents in channel flow.

• (3) A large portion of fluctuation kinetic energy is carried by velocity fluctuations associated with regions of potential flow, i.e. non-turbulent patches. These fluctuations associated with potential-flow regions occur at a large scale, and are possibly coherent motions associated with a time scale of several large-eddy turnover times. If such coherency occurs in the atmosphere, it has implications both for the flux measurements in the field and modelling of such flows with LES: the averaging period of flux measurements, respectively the grid size of an LES, would have to be chosen to allow for either resolution or complete filtering of these structures.

• (4) Global intermittency can occur without external perturbations of the flow: in our cases global intermittency is simulated despite the absence of finite-size triggers from synoptic conditions, low-level jets, or surface heterogeneities. Global intermittency is hence intrinsic to a stable atmospheric boundary layer beyond a certain stability, and we suggests it cannot be treated as an on–off process in time, but should rather be seen as happening in time and space.

Appendix: Initial Conditions and Transition Time Scales of the System

In Sect. 2.2 boundary conditions are discussed in detail, and transition time scales are identified as fundamental properties of the flow. The question how these time scales depend on the initial condition, i.e. are they determined by the strength of stratification only or also by the way in which it is imposed, we address through an additional set of simulations labeled with the prefix IC (Table 3). In these simulations the stratification is strengthened gradually beginning with the state of case W030S at $t^-\approx 2$. At this time, the case W030 has reached the state of a weakly SBL with strong stratification at the bottom and a capping inversion around $0.5\delta$. For cases IC150S and IC038S the state of this case at $t^-\approx 2$ is used as initial condition, and stratification is strengthened (which corresponds to multiplying the buoyancy profile by a constant larger than 1) such that $Ri_B=0.15$ respectively $Ri_B=0.038$. Similarly, the cases IC050S and IC075S were spun off IC038S and IC050S respectively (Table 3; Fig. 14a).

Table 3 Set-up of numerical simulations and their parameters for the study of initial conditions (see Table 1 for definition of the main cases; the cases S620 and I150 are repeated from that table)

Fig. 14

a Sketch illustrating the way in which initial conditions where studied; coloring matches panels b and c of this figure. b Vertically integrated TKE for the series of cases (IC038S, IC050S, IC075S, IC150S) described in a and Table 3. c Same as b but for vertically integrated spanwise vorticity RMS $\varOmega _{x}^2$

We find that our simulations vary on three time scales (Figs. 6, 14b):

• First, an initial adaption of the system to the strong perturbation from the surface that is seen in a reduction of both turbulence intensity and enstrophy magnitude in all cases. The strength of this initial response is stronger for stronger stratification and scales with the Brunt–Väisälä frequency at the surface, that is, $\sqrt{Ri_B}$. From our simulations we find for this initial transient

  $$\begin{aligned} \frac{1}{E_{\mathrm{neutral}}\sqrt{Ri_B}}\frac{\partial E}{\partial t^-} \approx -4, \end{aligned}$$

  (9)

  where the derivative is approximated as $(E(t')-E(0))/t'$ for $t'^-=0.2$ and $E(t')=\int e(t',z) dz$. This transient is only observed for the cases in which the stratification was concentrated into the viscous sublayer, and to some extent in the case IC150S, where the stratification is increased by a factor of five.

• Second,the fluctuation kinetic energy recovers on the time scale of several inertial periods starting between $t^-=0.5$ and $t^-=1$. During this recovery, the enstrophy fluctuation (Fig. 14b) grows until it reaches an equilibrium at a level close to, but smaller than the neutral reference. The time scale of this enstrophy recovery, discussed in Sect. 4, is about two inertial periods independent of stratification strength. The fluctuation kinetic energy initially grows at a larger rate than the enstrophy and is out of equilibrium until the end of our simulations; it enters a slow oscillation on the time scale of several inertial periods. In contrast to the enstrophy the fluctuation kinetic energy exceeds the neutral reference in all cases simulated here. This indicates significant contributions from weakly or non-vortical motions. The physical mechanism behind this is the excitation of waves in a stratified fluid which efficiently extract energy from the mean flow. In our cases, the significant contribution of wave energy of the flow is independent of the way in which the initial condition is imposed. We find, however, larger growth rates of TKE for stronger stratification.

• Third, all cases vary at the inertial time scale $2\pi /f$. Inertial oscillations are evident in time series of horizontal velocities (not shown). In the global statistics, a distinguished signal on the inertial time scale is only visible for the strongest stratified case. This is related to enhanced shear during certain phases in the inertial oscillation of the mean velocity.

We conclude that transients on the order of inertial periods are present in the SBL, even under weak stratification. To the extent covered by our simulations, this finding does not depend on the way in which the initial conditions are imposed. Hence, the SBL does not equilibrate over the course of a night. This also means that, irrespective of the initial condition, our simulations show the qualitative behaviour discussed in Sects. 4 and 5.