## Abstract

High-resolution large-eddy simulations of the Antarctic very stable boundary layer reveal a mechanism for systematic and periodic intermittent bursting. A nonbursting state with a boundary layer height of just 3 m is alternated by a bursting state with a height of ≈5 m. The bursts result from unstable wave growth triggered by a shear-generated Kelvin–Helmholtz instability, as confirmed by linear stability analysis. The shear at the top of the boundary layer is built up by two processes. The upper, quasi-laminar layer accelerates due to the combined effect of the pressure force and rotation by the Coriolis force, while the lower layer decelerates by turbulent friction. During the burst, this shear is eroded and the initial cause of the instability is removed. Subsequently, the interfacial shear builds up again, causing the entire sequence to repeat itself with a time scale of ≈10 min. Despite the clear intermittent bursting, the overall change of the mean wind profile is remarkably small during the cycle. This enables such a fast erosion and recovery of the shear. This mechanism for cyclic bursting is remarkably similar to the mechanism hypothesized by Businger in 1973, with one key difference. Whereas Businger proposes that the flow acceleration in the upper layer results from downward turbulent transfer of high-momentum flow, the current results indicate no turbulent activity in the upper layer, hence requiring another source of momentum. Finally, it would be interesting to construct a climatology of shear-generated intermittency in relation to large-scale conditions to assess the generality of this Businger mechanism.

## 1. Introduction

This study presents a mechanism for shear-generated intermittent turbulence in the very stable boundary layer (VSBL) based on a high-resolution large-eddy simulation (LES) study, which is representative for conditions on the Antarctic plateau (van der Linden et al. 2019). Here, with intermittent turbulence, we refer to “global intermittency” as defined by Mahrt (1999), where periods of “quiescent” flow are interrupted by sudden bursts of turbulence. It is shown that shear is built up and eroded in a natural, cyclic manner at the top of the boundary layer. The high-shear flow then generates unstable waves that lead to turbulent bursting, which in turn erode the shear itself by which they are generated. Subsequently, a relatively “quiescent” period follows during which shear is built up again. In this study, we will show that this shear-generated intermittency on short time scales (≈10 min) is a systematic and periodic feature in our VSBL, and how it contributes to the steady-state VSBL over longer time scales (>1 h).

Shear-generated intermittent bursting is a frequently observed phenomenon within the weak-wind, stable boundary layer (SBL) (see, e.g., Nappo 1991; Mahrt 1999). In spite of its omnipresence, the reasons behind such intermittent flow have remained unclear as no general dominating mechanism has been identified (Mahrt 1999). Analyzing intermittent flow (from observations) is difficult because the turbulent intensity within the background flow is extremely weak. As such, the effects of local heterogeneities or case-specific disturbances are amplified, and the cause or origin of the burst is easily obscured. Multiple triggers of these events have been identified in literature, for example, density currents and solitary waves (Sun et al. 2004, 2012), spatially dependent (de)coupling depending on local topography (Acevedo and Fitzjarrald 2003), or the interplay between radiative surface cooling and pressure-gradient induced mixing (van de Wiel et al. 2002). Other frequently observed causes are unstable internal gravity waves resulting from the Kelvin–Helmholtz instability (see, e.g., Gossard et al. 1970; Finnigan et al. 1984; De Baas and Driedonks 1985; Coulter 1990; Nappo 1991; Blumen et al. 2001).

Recently, Petenko et al. (2019) showed that successive wave disturbances frequently occur over periods exceeding several hours during the polar winter of 2012 at Dome C, Antarctica. It therefore appears to be a systematic feature of the long-lived Antarctic SBL. Using high-resolution sodar echogram observations, they were able to observe the finescale structure of such wave events and estimate both their characteristic temporal and spatial scales. In particular, they show that shear-generated wave disturbances occur under stationary conditions in periodic wave trains lasting 4–6 min even at supercritical bulk Richardson numbers (Ri_{b} > 0.25).

Nearly half a century ago, Joost Businger proposed a mechanism by which such shear-generated bursts could occur in the VSBL even at supercritical Richardson numbers (see Businger 1973). He conjectured that, if the supercritical Richardson number is reached at a particular height, vertical transfer of momentum and heat is blocked. Locally, shear at this height builds up as the wind below is decelerated (vertical momentum flux divergence) and wind above is accelerated (convergence). The shear between the upper layer and lower layer increases until the flow becomes hydrodynamically unstable, and a burst of momentum and heat toward the surface can occur. The shear is rapidly reduced and the flow becomes quiescent again until the next burst. We will show that Businger’s mechanism is largely applicable except for one component. Whereas the lower layer is, indeed, decelerated by momentum divergence, momentum convergence is not the major cause of flow acceleration in the upper layer (as turbulent activity is very weak). Instead, the acceleration above is caused by the combination of the pressure force and wind turning due to the Coriolis force rather than by the assumed momentum convergence.

The favorable conditions in the Antarctic winter at the Dome C site may prove to be key in identifying the mechanism for shear-generated intermittent bursts and its periodic occurrence. During these Antarctic winter months (June to August), the SBL at Dome C may reach long periods (lasting for several days) of “steady state” during which the wind and temperature profiles do not change significantly over time (Vignon et al. 2017b; Baas et al. 2019). The potential of the arctic regions has already been recognized before (see, e.g., Dabberdt 1970; Grachev et al. 2005, 2008) as they may serve as “natural laboratories” for the study of the SBL—in particular—during the winter months when the daily cycle is absent.

By analysis of observations in combination with high-resolution LES, van der Linden et al. (2019) show that, in the Antarctic, a thermal steady state is possible when the turbulent cooling of the SBL as a whole is balanced by heating through large-scale subsidence (see also Vignon et al. 2018; Baas et al. 2019). The close correspondence between their LES results and the observations encourages the use of LESs for in-depth process studies. In contrast to the observations, within LESs, the boundary conditions and forcings can be fully controlled. Therefore, the LES approach is an attractive complementary tool to study the Antarctic SBL and the associated intermittency found by Petenko et al. (2019). However, they can only be considered complementary as simulations remain heavily idealized compared to complex reality with respect to, e.g., their forcing and surface boundary conditions (Bosveld et al. 2014). Such simulations have also been used before to study shear-generated instabilities in realistic settings, for example, based on CASES-99 (Zhou and Chow 2014) and the Beaufort Sea Arctic Stratus Experiments (Na et al. 2014).

Ideally, a full mechanistic analysis of intermittency directly from observations would be preferred. Unfortunately, measuring such bursts is complicated by the harsh, cold conditions that make accurate measurements of turbulent fluxes by standard sonic anemometers challenging (Vignon et al. 2017a). Also, the sodar echogram data, although measured at relatively high spatiotemporal resolution, are not easily transformed into quantitative fluxes (Petenko et al. 2019).

Somewhat surprisingly, the aforementioned LES results of van der Linden et al. (2019) indeed show the presence of periodic turbulent bursts in the VSBL similar to those reported by Petenko et al. (2019). Although the SBL is found to be in steady state with respect to its bulk quantities on an hourly basis, closer inspection reveals that the SBL is found to be periodically modulated by episodes of enhanced turbulence originating at the top of the boundary layer on time scales of ≈10 min. These events subsequently spread both upward and downward resulting in a temporarily larger boundary layer height and surface fluxes, respectively.

Here, we further investigate these top-down bursting events in the VSBL case of van der Linden et al. (2019) and show that they are the result of wave breaking after initial growth of a shear-generated instability. Using an extended simulation (viz., with a larger domain) at a high resolution (Δ = 0.08 m), the dominant wavelength is identified. The instability of this wave is confirmed by applying a linear stability analysis (LSA) on the background flow. Finally, we will identify the full intermittency cycle: the mechanism of wave growth, bursting, and erosion of the shear layer as well as the restoring mechanism to restore local shear again.

## 2. The steady Antarctic boundary layer?

In this section, we further investigate the LES case of the VSBL of van der Linden et al. (2019) to show the presence of intermittent turbulence. A short overview of their VSBL simulation can be found in appendix A. A comprehensive description of the observations and model formulation may be found in van der Linden et al. (2019).

Figure 1 shows the vertical profiles of the wind speed, potential temperature, kinematic temperature flux and the contributions to rate of change of potential temperature averaged over the final hour and the horizontal plane of the simulation. Van der Linden et al. (2019) show that on average (viz., averaged over simulation periods ≥1 h), a thermal steady state exists in which cooling of the boundary layer by vertical divergence of the kinematic temperature flux is balanced by subsidence heating of the air. The heating rate of subsidence has a maximum at approximately 3.75 m, and decreases to zero toward the surface and top of the domain where the imposed subsidence velocity and temperature gradient are zero, respectively. Apparently, the profile of the temperature flux “adapts” itself to the profile of subsidence heating, as the latter is a slower process (viz., the average temperature gradient changes over time scales longer than the typical time scale of turbulent mixing).

Although the LES case is found to reach a thermal steady state averaged over periods ≥1 h, closer inspection indicates a thermal steady state does not exist over averaging periods of approximately 10 min or shorter. Figure 2a presents the horizontally averaged kinematic temperature flux as function of time and height during the final simulation hour of the original LES case. The temperature flux exhibits clear periodic behavior in which events of enhanced temperature flux are superimposed on a relatively “quiescent” base state (i.e., a shallow SBL of depth *z* ≈ 2.5–2.9 m). These main bursting events appear to consistently start at the top of the boundary layer, and subsequently extent both upward and downward. After approximately 200 s, the enhancement of the temperature flux has largely disappeared, although some enhanced values are still observed near the surface <2 m. The time between the onset of these successive events is approximately 600 s. Similar patterns are also present in, for example, the horizontally averaged momentum fluxes and temperature variance. Conceptually, a “short” time scale of ≈200 s (or “fast” process) can be defined in which the bursts affects the mean flow, and a “long” time scale of ≈400 s (or “slow” process) in which the conditions favorable for the subsequent burst are created. As the magnitude of the bursts is relatively small [e.g., *O*(10) W m^{−2} in the heat flux], changes in the first-order statistics such as wind speed and temperature remain modest as well during an event; the standard deviations over the entire simulation hour are *σ*(*U*) < 0.04 m s^{−1} and *σ*(*θ*) = 0.29 K (not shown).

Figure 2b shows the temporal relation of the kinematic temperature flux at the surface (green) and at the top of the SBL (i.e., at a height of 2.72 m; purple). During a burst, the 2.72-m flux rapidly exceeds the magnitude of the surface flux. On the contrary, the variation in the surface temperature flux is <6%, which indicates that bursts barely reach the surface. For convenience, we will define two states according to these two fluxes. The bursting intervals are defined as the periods in which the magnitude of temperature flux at 2.72 m exceeds the value of the surface flux. These intervals are indicated by the dash–dotted lines in Fig. 2b.

The observed behavior of the temperature flux is consistent with the formation and breaking of traveling waves at the interface of the turbulent boundary layer and the air aloft, which is confirmed by vertical cross sections from the simulation (see section 3a). During the initial stages of the bursting event, (linear) wave perturbations form and grow in time until nonlinear effects become dominant and cause wave breaking. Subsequently, turbulent kinetic energy is generated at this interface, which causes the boundary layer to grow in height (see Fig. 2a). Relatively warm and fast air is entrained into the boundary layer resulting in a net transport of both energy and momentum toward the surface (cf. Fig. 2b). This resembles the “upside-down” boundary layer as observed during the CASES-99 experiment (Mahrt and Vickers 2002).

In addition to these main events (at *z* ≈ 2.5–2.9 m), a secondary event appears to be initiated above the turbulent SBL in response to the first events (see *z* ≈ 4.5 m, *t* ≈ 2800 s). This secondary event is weak compared to the main events and appears not to penetrate deep into the base state. Its peak values are about 20% of those of the main events. Such secondary events appear to occur sporadically in the simulation. It is unclear if they result from a separate instability or from residual turbulence of the main bursting events. The turbulence on average is weak or even absent at higher levels in the flow (*z* > 4 m), and the flow can be regarded as “quasi laminar” and decoupled from the surface layer (cf. with Banta et al. 2007). Therefore, residual turbulence ejected by the main events may take relatively long to dissipate. Due to its sporadic occurrence and weak impact, these secondary events are discarded in the main analysis.

## 3. Wave analysis

In this section, an in-depth analysis of the wave phenomenon is made. First, wave characteristics are diagnosed from the simulations. The dominant wavelength is extracted by spectral analysis of the vertical cross sections of the velocity field (see section 3a). Second, a linear stability analysis is applied in section 3b to show that the background flow is indeed unstable in time with respect to this wave perturbation, and that the wave growth enables turbulent bursting.

### a. Spectral analysis

To identify wave properties, such as the wavelength or amplitude, vertical cross sections of the simulation are analyzed. Before applying the Fourier transform to find the dominant wavelength from the horizontal velocity fields, first, a simple visual inspection is made. These suggest a wavelength of approximately 16–19 m in the original VSBL simulation of van der Linden et al. (2019) (not shown). Unfortunately, the full horizontal extent of the domain in their simulations amounts to only *L*_{x} = 19.2 m (with an isotropic grid spacing of Δ = 0.08 m). Therefore, the wavenumber bin resolution (i.e., its detectable change) Δ*k* is equal to 0.3272 m^{−1}, and accurate determination of the expected wavelength using spectral analysis is unfeasible.

To alleviate this problem, the original simulation is extended in both the horizontal directions according to the following procedure. First, five copies of the original simulation field at *t* = 23 h are pasted together in the *x* direction. Second, this “new” field is duplicated and joined in the *y* direction. Gaussian noise (*μ*_{G} = 0; *σ*_{G} = 0.02*σ*_{i}) is added as a random perturbation, where *σ*_{i} is the height-dependent standard deviation of the variable considered. The perturbation is added to ensure that turbulent fields will not be identical while keeping the averaged state unchanged. This is done for all three velocity components and the temperature. The simulation is restarted on the bigger domain with new domain sizes *L*_{x} = 96 m, *L*_{y} = 38.4 m and *L*_{z} = 19.2 m, and is allowed to freely evolve for 2 simulation hours. Only the second simulation hour is used for the analysis as the first simulation hour may be influenced by initial correlation between the individual field copies. As multiple wave cycles have passed during the first hour, it is assumed that these “memory effects” of the artificial initialization have disappeared after the first hour (cf. Fig. 2). The grid spacing is kept at 0.08 m. The extended simulation results in a wavenumber bin resolution of Δ*k*_{x} = 0.065 m^{−1} in the *x* direction (along the isobar) after application of the Fourier transform.

Figure 3 shows the perturbation of the *x* component of the velocity *u*′ at different times during a full cycle. Here, the velocity perturbation is defined as the difference between the local, instantaneous velocity and the horizontally averaged value. It is observed that at the top of the boundary layer with height of approximately 2.7 m, a wave pattern of alternating positive and negative velocity perturbations forms (cf. Fig. 3b). Subsequently, the wave amplitude grows in time and eventually breaks triggering more (vertical) turbulent mixing, which leads to an increase of the turbulent boundary layer (cf. Figs. 3c,d). During the later stages of the event, the wave patterns have disappeared and the boundary layer has grown to approximately 5.5 m with overall increased values of the velocity perturbation indicating an increase in turbulent activity (see Figs. 3e,f). Finally, the turbulent activity at the top of the boundary layer dissipates, and the boundary layer returns to its preburst state (Fig. 3g). Similar evolutions are observed for the perturbations of the cross-isobaric velocity component *υ*′, the vertical velocity component *w*′, and the potential temperature *θ*′ (not shown).

To determine the dominant wavelength, the stages similar to Fig. 3b are selected from the final three events (out of a total of four) and analyzed. The first bursting event in the second simulation hour is discarded since it may be influenced by a secondary event (cf. Fig. 2a). Using a similar approach as Newsom and Banta (2003), the normalized power spectra at each height *z* are computed by taking the one-dimensional Fourier transform in the *x* direction of each cross section. Individual spectra are added and normalized by its maximum value. The Fourier components of the perturbation of a variable are indicated by the hat symbol. For example, $\theta ^k$ refers to the Fourier component at wavenumber *k*_{x} ≡ 2*π*/*λ*_{x} (with *λ*_{x} the wavelength in the *x* direction; *k*th mode) of the perturbation in the potential temperature *θ*′.

Figure 4 presents the normalized power spectra of both the vertical velocity component $w^k$ and potential temperature $\theta ^k$ as a function of both height *z* and wavenumber *k*_{x}. For clarity, only wavenumbers up to *k*_{x} = 1 m^{−1} are shown (out of a maximum of *k*_{x} = 39.22 m^{−1}) as the power at higher wavenumbers is negligible. The spectra of $w^k$ and $\theta ^k$ have their maxima at *k*_{x} = 0.3274 m^{−1} and *z* = 4.36 m, and at *k*_{x} = 0.3929 m^{−1} and *z* = 3.08 m, respectively. The location of the maximum of $w^k$ corresponds to a wavelength of *λ*_{x} = 19.2 m, while the location of the maximum of $\theta ^k$ corresponds to *λ*_{x} = 16.00 m. Hence, this analysis confirms the aforementioned visual inspection.

Both power spectra have an approximately equal horizontal extent. This indicates that the wave phenomenon is composed of multiple wavelengths in a narrow range. The vertical extent of the power spectrum of $w^k$ appears to be larger than that of $\theta ^k$. No explanation is found for this difference in height of the distribution. The vertical profiles of the power spectra at the dominant wavenumber are shown in Figs. 5a and 5b.

Visual inspection of the simulation field at *t* = 5100 s shows that the propagation direction of the primary wave events is *ϕ* ≈ 0° with respect to the isobars; that is, aligned with the *x* axis (not shown). For convenience, it is therefore taken as 0°. Unfortunately, no value of the complex phase speed can be calculated due to the limited frequency at which cross sections and simulation fields were saved, namely, 30 and 300 s, respectively.

### b. Linear stability analysis

Linear stability analysis provides information about the hydrodynamic stability of small perturbations (indicated by the prime) subject to a given background flow. Arbitrarily shaped perturbations of small amplitude are typically present in “quiescent,” nonturbulent background flows, and can be seen as a superposition of sinusoidal waves (Fourier decomposition). By LSA, one investigates if these wave components (modes) decay or grow in time (i.e., have a negative or positive growth rate). If all modes contained in the Fourier decomposition decay, the flow is said to be stable. However, if a number of modes grow (exponentially in time), it is assumed that the fastest growing mode of these will rapidly dominate over the others and continue to grow until secondary instability mechanisms cause that wave to break and overturn. An extensive overview on LSA can be found in Drazin and Reid (2004) and Kundu et al. (2012). Although LSA is traditionally used to investigate the stability of strictly laminar flows and predict their transition to turbulence (Kundu et al. 2012), the LSA approach has been stretched in its assumptions by applying it to flows that are not completely laminar, but are “smooth” with respect to their very weak turbulent activity. In those cases, LSA is used to analyze whether the mean flow (in a Reynolds-averaged sense) supports unstable wave modes that will lead to turbulence of more significant magnitude. Indeed, LSA has been employed with success to predict shear-generated instabilities in the “smooth,” but weakly turbulent SBL using the observed mean states (see, e.g., Finnigan et al. 1984; De Baas and Driedonks 1985; Newsom and Banta 2003).

#### 1) Method

Here, we briefly explain the implementation of the LSA. A detailed description can be found in appendix B. First, it is assumed that, at a given time, the wave perturbations propagate along one direction in the horizontal plane. This reduces the 3D problem to a 2D approximation. Note that, this assumption excludes the Coriolis force from the analysis. This simplification is motivated by the magnitude of the perturbation Coriolis term after linearization, which is negligible compared to the other terms. Second, we assume the flow to be inviscid. The velocity vector is then rotated over angle *ϕ*, which corresponds to an alignment of the flow with the propagation direction (here, *ϕ* ≈ 0°; section 3a). The mean 2D background states of wind speed and temperature are given by **U** = {*U*(*z*), 0} and Θ(*z*), respectively. Traveling-wave solutions are assumed for the wave disturbances. For example, for the vertical velocity component

where *k* is the wavenumber, $w^k\u2061(z)$ is the complex amplitude (profile) of the *k*th mode, *c*_{k} = *c*_{k,R} + *ic*_{k,I} is the complex phase speed. Additionally, *σ*_{k} = −*ikc*_{k} is introduced for convenience. For a mode to be unstable, the real part of *σ*_{k} has to be >0 s^{−1}. Our LSA model investigates the stability of a single mode solving for the unknown *σ*_{k}, and the corresponding profiles of the vertical velocity and temperature perturbations for a given *k* and *ϕ* of that mode. Apart from these, boundary conditions for the vertical velocity component have to be specified. Here, we require that the vertical velocity component is zero at the bottom ($w^k=0$; no penetration) and that the solution remains bounded for infinite height (viz., $w^k$ tends to a constant value). The latter boundary condition is approximated by $dw^k/dz=\u2212kw^k$ at the top of the computational domain. This condition automatically ensures that the solution has an exponential decrease to zero at infinite height, while recognizing that the actual boundary condition is imposed at finite height. Note that, both *d*^{2}*U*/*dz*^{2} and *d*Θ/*dz* tend to zero here (cf. De Baas and Driedonks 1985; Newsom and Banta 2003). The equation for the temperature perturbation can be eliminated by further substitution and would yield the classical Taylor–Goldstein equation (see, e.g., Newsom and Banta 2003), which is a second-order equation in $w^k$ requiring two boundary conditions. Here, this elimination is not done for convenience.

The system of equations is discretized in the vertical direction using *N*_{z} levels (the same as in the simulation), and transformed into a generalized eigenvalue problem with eigenvalue *σ*_{k} and eigenvector $[w^k,\theta ^k]T$. Solving the generalized eigenvalue problems gives 2*N*_{z} pairs of eigenvalues and eigenvectors of which *N*_{z} are independent. For each pair, its complex conjugate is also a valid solution with opposite growth rate (Kundu et al. 2012). The most unstable eigenvalue-eigenvector pair [largest Re(*σ*_{k})] is selected as it is expected to dominate the flow evolution.

#### 2) Result

We investigate the stability of waves with wavenumber and propagation direction set equal to *k* = 0.3929 m^{−1} and *ϕ* = 0°, respectively (see section 3a). The background profiles of wind speed and temperature are obtained by averaging the simulated profiles between *t* = 5400 and 5700 s from the extended domain simulation. This interval is approximately halfway between two successive bursts (based on the 2.72-m temperature flux) and is representative of the base state. The background wind speed profile is then projected onto the plane of propagation, which corresponds to setting *U*(*z*) = *u*(*z*) in our case. Using these parameters, this investigated mode is found to have the fastest growing eigenvalue *σ*_{k} = (0.0195 − 0.7899*i*) s^{−1}. This corresponds to complex phase speed components of *c*_{k,R} = 2.01 m s^{−1} and *c*_{k,I} = 0.05 m s^{−1} (see Fig. B1), where the subscripts *R* and *I* represent the real and imaginary parts, respectively. The wave speed *c*_{k,R} equals the speed of the background flow at *z* ≈ 2.92 m, so that the midplane of the wave does not move in a coordinate system moving with that flow speed. The *e*-folding time scale for exponential growth is ≈51 s (i.e., $k\u22121ck,I\u22121$). Although this time scale cannot be accurately determined from the simulations due to the limited output frequency of the cross sections, it appears to be reasonable compared to the time scale of the bursting event (cf. Fig. 2b). A strict comparison is not possible as the linear growth regime is violated relatively soon due to fast growth of the wave.

Figure 5 shows the normalized wave mode profiles of the vertical velocity component, temperature, vertical wave momentum flux and vertical wave temperature flux at the dominant wavenumber *k* = 0.3929 m^{−1} as inferred from the simulation (blue) and calculated by the LSA (red). Here, the vertical “fluxes,” resulting from the temporal growth of the wave amplitude, are calculated as the real part of the product of the variable considered and the complex conjugate of the vertical velocity component $w^k*$ (cf. Newsom and Banta 2003). This product is the generalization of the dot product for complex numbers. For temperature, this product represents that part of the temperature perturbation that is in phase with the perturbation of the vertical velocity component. Note that, for nongrowing (linear) waves [Re(*σ*_{k}) = 0 s^{−1}] this product is zero (viz., $\theta ^k$ lags 90° with respect to $w^k$), and, as such, no scalar or momentum is transported. However, for growing waves this product is nonzero. Physically, the vertical velocity does not change sign at the moment the densest (lightest) fluid is displaced through the midplane in a wave of which the amplitude grows in time. The presence of an in-phase component (nonzero product) follows from the LSA model equations [see Eq. (B7b) in appendix B]:

which shows that $Re\u2061(\theta ^kw^k*)\u22600$ if and only if Im(*c*_{k}) ≠ 0.

The calculated LSA profiles resemble those estimated from the simulation for all four variables to a high degree. Minor differences are mainly found near the surface, which are likely caused by some irregular motion (weak turbulence) of minor amplitude. The LSA-calculated profile for $||w^k||$ smoothly tends toward zero near the top of the domain, whereas the profile estimated from the simulations does not. Because turbulent activity is virtually absent in the upper half of the domain (cf. Fig. 1), this might indicate some wave activity there (possibly caused by minor reflections). As such, the domain is not large enough to fully exclude boundary effects, although these effects are believed to be minor.

A local minimum of $||w^k||$, and the maxima of $||\theta ^k||$ and the wave fluxes are present at *z* = 2.88 m coinciding with the inflection point of the velocity profile *U*(*z*). This height is a critical level of the flow: the real part of the phase speed *c*_{k,R} is equal to the local horizontal velocity at this height. The large, narrow peaks of the wave momentum and temperature fluxes indicate that large parts of $u^k$ and $\theta ^k$ are in phase with $w^k$ at this height, whereas they are out of phase near the surface and above the SBL (see Figs. 5c,d).

The profiles correspond to those found by De Baas and Driedonks (1985) (vertical velocity and temperature) and Newsom and Banta (2003) in a nondimensional form. The shape and structure of these profiles are consistent with a Kelvin–Helmholtz-type instability (Newsom and Banta 2003). This confirms that the wave formation and wave breaking (cf. Fig. 3) are indeed the result of a shear instability at the top of the SBL.

## 4. Mechanism behind the full cycle

In spite of the close correspondence between the LSA and the simulation results, the previous section merely confirms that the background flow is unstable for perturbations at the dominant wavenumber. It does, however, not reveal how the boundary layer responds during the burst and relaxes back to its base state. In this section, this process is analyzed by conditional averaging over the bursting and the nonbursting periods. First, the effect of the intermittent burst on the mean flow is shown. Second, the evolution of the boundary layer after a burst is presented and, in particular, it is revealed why the process of shear-generated intermittent bursts is periodic.

### a. Flow evolution during the burst

Figure 6 shows the flux profiles of momentum *F*(*u*_{i}) and temperature *F*(*θ*), and the contributions to the tendencies of the isobaric velocity component *u* and temperature *θ*. These values are conditionally averaged on the bursting states taken from the final simulation hour of the original VSBL simulation (see van der Linden et al. 2019). These contributions for *u* are the divergence of the total isobaric momentum flux and the *x* component of the Coriolis force, whereas the contributions for *θ* are the divergence of the kinematic temperature flux and heating by subsidence. The *x* component of the Coriolis force is given by *f*_{C}*υ*. This term does not “add” momentum (or energy) to the flow as the Coriolis force is always perpendicular to the wind vector. However, it can rotate the wind vector thereby transferring momentum (and energy) from the *y* direction to the *x* direction (and vice versa) in the case of a force imbalance. At the same time, the imposed pressure gradient force steadily adds momentum to the cross-isobaric direction (*y* direction). The bursting (nonbursting) state is defined as those time intervals in which the absolute value of the 2.72-m temperature flux is larger (smaller) than the absolute value of the surface temperature flux (cf. Fig. 2b). The total fraction of the time the SBL resides in the bursting (nonbursting) state is 26% (74%).

Both the isobaric momentum flux (*x* direction) and temperature flux exhibit large negative peaks centered around 2.72 m (cf. Fig. 5) exceeding the surface values. The averaged vertical extent of the peaks is approximately equal to 4 m and is dependent on the time during the burst: after the initial wave breaks, momentum and heat are progressively mixed in the vertical direction. As a result of the burst, turbulent kinetic energy is generated and the boundary layer height increases up to ≈5.5 m. Additionally, the base state becomes temporarily “coupled” to the layer above.

The *x* component of the Coriolis force *f*_{C}*υ* and the heating by subsidence both have a positive contribution to the tendencies of *u* and *θ*, respectively, and tend to zero for *z* > 6 m (see Figs. 6c,d). The contributions as a result of the flux divergences show a more complicated pattern: they are mainly positive in the lower layer and negative higher up. In total, a net acceleration and warming of the SBL occur below ≈2.9 m, whereas in the region 3–5 m (relatively) strong deceleration and cooling occur. The vertical transport by bursting, hence, reduces both the difference in the velocity magnitude and the temperature between the upper and lower layers. As the relative decrease of the shear is larger than the decrease in thermal gradient, the cause of instability is counteracted (see Fig. 8).

In addition to the isobaric velocity component, also changes in the cross-isobaric component occur. Although the profile of the total rate of change of the cross-isobaric velocity component exhibits a more complicated structure, its values are typically <50% of the total rate of change of *u* and have a relatively small contribution to the change of the total shear squared *S*^{2} (not shown).

### b. Flow evolution after the burst

In the nonbursting state, the flux and total rate-of-change profiles are markedly different than in the bursting state (see Fig. 7). The profiles of the momentum and temperature fluxes indicate that the main turbulent layer is now approximately 3 m in depth.

It is found that the profiles of the *x* component of the Coriolis force *f*_{C}*υ* and the heating by subsidence do not significantly differ in the nonbursting state as compared to the bursting state. The figures, however, do differ with respect to the turbulent flux contributions (see Figs. 7c,d). In absence of momentum and heat transport from above, the lower layer (*z* < 3 m) decelerates due to the surface friction (momentum flux divergence) and cools by the surface temperature flux (temperature flux divergence). This lower layer corresponds to the active turbulent layer in the nonbursting state, whereas the layer above can be regarded as “quasi laminar.” At the same time, this quasi-laminar layer experiences a net acceleration and warming by the Coriolis force and subsidence heating. The overall result is that the contrast between the lower and the upper layer increases with respect to the wind speed and temperature (i.e., an increase of the local shear and temperature gradient around *z* ≈ 3 m). The momentum transferred from the cross-isobaric to the isobaric direction by the Coriolis force is steadily replenished by the pressure force in the *y* direction (not shown).

Weighted by their respective fractions of occurrence, the deceleration and cooling, and the acceleration and warming balance in time. As such, both a steady state in the amount of momentum and a thermal steady state result, when averaged over, for example, times >1 h (see van der Linden et al. 2019, their Fig. 8). The intermittent bursts of the SBL, therefore, contribute to this thermal steady state in the presence of heating by subsidence. Periodically, they “entrain” relatively warm air heated by subsidence into the turbulent layer. Mirocha et al. (2005) already provided compelling evidence that warm air entrained into the boundary layer by subsiding motions balances a significant part of the turbulent heat flux near the surface in the observed Arctic clear-sky SBL. Similarly, a LES case based on this Arctic SBL (Mirocha and Kosović 2010) shows that the inclusion of subsidence resulted in a nearly thermal steady state. However, they did not report any (periodic) bursts within the SBL.

The impact of the bursting and the nonbursting phases on the mean quantities are summarized in Fig. 8. This figure shows the profiles of the total shear squared *S*^{2}, Brunt–Väisälä frequency *N*^{2}, and the gradient Richardson number Ri_{g} = *N*^{2}*S*^{−2} representative of different times during one cycle (just before a burst and after the burst) of the original VSBL simulation. The temporal variation in *S*^{2} and *N*^{2} result in clear changes of Ri_{g} over the shear layer during a cycle. Finally, a conceptual picture of the mechanism and its main actors is given in Fig. 9.

## 5. Discussion

### a. Comparison with suggested mechanisms

The current results suggest a systematic mechanism by which cyclic intermittent bursts are triggered by a Kelvin–Helmholtz instability at the interface of a shallow SBL and the quasi-laminar layer above. Similar mechanisms (or parts thereof) have been reported in literature. Yet, a comprehensive, observationally based explanation by which multiple intermittent bursts may occur successively or even periodically within an uninterrupted time span has not been given (Mahrt 2014). Indeed, systematic observations of such successive bursts may be difficult due to both observational limitations and nonstationarity of the SBL itself in the midlatitudes.

The mechanism identified in this study resembles the mechanism reported by Newsom and Banta (2003). They show that, just prior to the burst, the shear dominates the reduction of the Richardson number causing the flow to become locally unstable. In particular, the buildup of shear over a relatively small vertical extent triggers a Kelvin–Helmholtz instability. Furthermore, they find a net increase of Ri during the wave event as both shear and temperature gradient are mixed, and a small decrease of Ri after the wave event for which no cause is identified. This observation appears to correspond with our simulations, although a direct comparison is difficult due to observational limitations (e.g., determining gradients from discrete levels) and the number of events (one in their case).

Also, similarities and dissimilarities between the mechanism of van de Wiel et al. (2002) and the current mechanism are present. Van de Wiel et al. (2002) acknowledge the potential role of the ageostrophic pressure gradient (i.e., the effective pressure gradient in the direction of the mean wind) as a main external parameter governing intermittency in their bulk model. The main difference, however, is that their bulk model cannot capture the instability and the dynamics at the interface of the SBL and the quasi-laminar layer above, but considers a suppression of the turbulent activity of the SBL as a whole (governed by the bulk Richardson number). The present results, on the contrary, provide compelling evidence for a two-layer structure with separate dynamics: whereas no turbulence is present above the interface and the flow accelerates there, the SBL itself decelerates as a result of the surface friction in the nonbursting state. As such, it appears that the mechanism of van de Wiel et al. (2002) is less realistic. At the same time, our simulation imposes a fixed surface temperature via the boundary condition, whereas, in van de Wiel et al. (2002) the thermal balance of the surface is an active (dynamic) part of the system, which may allow for additional surface feedbacks not considered here.

Finally, the mechanism found in this study is remarkably close to the conjecture of Businger (1973). Here, we cite parts of his conjecture:

The point is that if

R_{fcr}is reached sometime … it will be reached first where the maximum value occurs at some height above, but relatively close, to the surface. As soon as this happens the turbulence will be dampened and a laminar layer will tend to form. This layer is an effective barrier for all the fluxes.… Under the laminar layer the transfer of momentum will continue down to the surface until the available momentum is depleted orR_{f}has become larger than critical. The result is that the wind diminishes and a period of calm sets in.… In the meantime, above the laminar layer momentum is still transferred downward whereas little heat is transferred. Consequently, the momentum increases in the upper part of the laminar layer because it cannot pass through this layer. A strong wind shear builds up and since there is no similar effect for the heat flux, Ri must decrease, eventually reaching a value below Ri_{cr}. This means that the laminar layer will gradually be eaten away by turbulence from above. Eventually the turbulence reaches the ground associated with a burst of momentum and heat. After this, the entire sequence of events may repeat itself.

However, the key difference is the actor that increases the momentum above the boundary layer during the nonbursting times. Whereas Businger suggested that momentum is transferred downward from higher up in the flow by stress convergence, the current results indicate momentum is increased by acceleration as a result of the pressure gradient and subsequent rotation by the Coriolis force, which are a rather constant factor in time. Apart from this difference, his conjecture is correct with regard to the origin of the burst, the deceleration in the bulk of the SBL and the possibility of periodicity. Finally, we acknowledge the fact that our studied long-lived ABL may differ from the midlatitude diurnal ABL, where additionally the nocturnal momentum budget might by influenced by, e.g., decaying convection and inertial oscillations.

### b. A systematic climatology of bursts?

The present study would largely benefit from a systematic climatology of bursts. Such climatology may clarify under which conditions successive or even periodic bursting events can occur. In our simulations, the external forcings (e.g., the geostrophic wind speed and subsidence profile) are kept constant, and the surface is homogeneous. As a result, the simulation reaches a steady state in which the bursts only marginally affect the background flow allowing a fast recovery and subsequent burst. A strict steady state is not expected to occur in the outdoor environment in which synoptic disturbances occur, but may be approached for several days in the polar regions.

Petenko et al. (2019) show that periods lasting several hours in which the SBL is perturbed by successive wave events, are frequent at the Dome C station (see their Fig. 9). They note that the large-scale weather conditions were stationary during these periods. However, the lack of accurate turbulent flux measurements (e.g., using eddy-covariance techniques) and the limited amount of measurement levels on the meteorological tower prevent the determination of the interactions between the mean flow and the wave events.

Another open question relates to the climatology of the event in relation to external forcings. In contrast to the present study, Petenko et al. (2019) seem to suggest that intermittency is more likely to occur within SBLs of depth 20–70 m than in very shallow SBLs of depth *z* ≈ 5 m. This implies a larger geostrophic forcing (i.e., near-surface large-scale pressure gradient), and/or a weaker subsidence warming as to allow a larger turbulent activity and a deeper SBL. At the same time, however, the time scales of successive event in their study corresponds to the time scale identified in the present study: 8–15 min in theirs as compared to 10 min in the current. Therefore, a one-to-one comparison of intermittency climatology with respect to forcings between observations and modeling is essential in order to generalize the present conclusions.

## 6. Conclusions

In this study, a mechanism for periodic shear-generated intermittent bursts is identified using high-resolution LES. This mechanism closely resembles the mechanism proposed by Businger (1973) differing only in the cause of acceleration above the SBL.

Van der Linden et al. (2019) simulate the VSBL based on observations of the Antarctic winter of 2015 from the Dome C station in related work. They show that the temperature flux divergence and heating by subsidence balance over time scales >1 h such that a steady-state SBL with depth ≈5.5 m is reached. Here, we find that the SBL is not in steady state over time scales <10 min, but is rather modulated by turbulent bursts, which enable the steady state over longer time scales.

Using an extended simulation domain, it is found that periodically wave perturbations form at the interface of a shallow SBL (2.5–2.9 m) and a quasi-laminar layer above (i.e., flow with negligible turbulent activity). The dominant wave is found to grow in time until it breaks resulting in increased turbulent activity and a temporary growth of the active turbulent layer. Spectral analysis shows that the wavelength of this dominant wave is 16–19 m.

A linear stability analysis confirms that small-amplitude waves of this wavelength are indeed unstable with respect to the mean wind and temperature profiles. Furthermore, the predicted perturbation profiles of the velocity components, temperature and fluxes correspond with those obtained by the spectral analysis. The shape of these perturbation profiles are indicative of the Kelvin–Helmholtz instability, which has been found to occur before in stable conditions (see, e.g., De Baas and Driedonks 1985; Newsom and Banta 2003).

The instability is created by an increase of the local shear at the interface that dominates over the increase in temperature gradient resulting in a decrease of Ri_{g} to a value <0.25, which is a prerequisite for instability to occur (see, e.g., Kundu et al. 2012). The interfacial shear is increased as a result of deceleration of the flow in the SBL by turbulent friction, and acceleration above by the combined action of the pressure forcing and the rotation by the Coriolis force. During the burst, these two layers become temporarily coupled and the momentum is exchanged; that is, the lower part accelerates and the higher part decelerates. The instability is mixed away by its own result and Ri_{g} becomes >0.25 at the interface. It is found, however, that the mean wind is only altered slightly by the burst and returns to its preburst state. As such, the flow is found to reside around its critical state, and a cyclic process of instability formation and bursting ensues. This is a (modified) Businger mechanism. Businger (1973) correctly proposed such intermittency could be periodic by the process described above with one exception. He stated that the momentum above the SBL is increased due to downward turbulent transfer. However, such transfer is not possible as a result of negligible turbulent activity above the SBL. It is important to note that, apart from turbulence, the potential impact of subsidence on the momentum budget is not considered in the present study, which can result in an additional supply of momentum.

The temperature dynamics follow a similar pattern. Prior to the burst, the SBL is cooled by the turbulent flux toward the surface and the quasi-laminar layer is heated by the subsidence heating. During the burst, the cooler air is mixed upward and the warmer air is mixed downward. It is this periodic mixing that explains the thermal steady state over time scales >1 h reported by van der Linden et al. (2019).

Although intermittent bursts are commonly observed in the SBL at both the mid- and high latitudes, the exact conditions leading to such bursts, and, in particular, successive (periodic) bursts remain elusive. At the same time, while the steady forcing conditions of the simulation allow periodic bursts to occur and the mechanism to be revealed, these conditions are just one realization of the SBL based on observations from the Antarctic winter and a sensitivity study in which these conditions, such as, the geostrophic wind speed, are systemically varied is recommended. A detailed climatology of shear-generated bursts in relation to the conditions in which they are found (e.g., mean wind or surface characteristics) would therefore be beneficial, and help to predict the time scales and vertical extent of such bursts among other things. Furthermore, realistic high-resolution simulations based on such climatological cases can clarify the contribution of bursts to vertical transfer of momentum and scalars.

## Acknowledgments

This work was sponsored by NWO Exact and Natural Sciences for the use of supercomputer facilities. We gratefully acknowledge funding by the European Research Council through an ERC-Consolidator grant (648666).

### APPENDIX A

#### Description of the LES Case

In the current study, the LES case for the VSBL of van der Linden et al. (2019) is used. Here, we briefly summarize the setup of their VSBL simulation. A detailed description of the observations, setup and results can be found in van der Linden et al. (2019). Furthermore, the used, open-source code MicroHH (http://microhh.org) is described in van Heerwaarden et al. (2017).

The subfilter-scale flux tensors are modeled by a Smagorinsky–Lilly-type eddy-viscosity model (Lilly 1962; Smagorinsky 1963) in which stratification effects are included (Lilly 1962; Mason 1989). Furthermore, the wall correction of Mason and Thomson (1992) is used for the length scale of the eddy viscosity. Surface fluxes are calculated using Monin–Obukhov similarity theory with the similarity functions of Högström (1988). Velocity boundary conditions for the horizontal components are no-slip at the surface and stress-free at the top, and no-slip at both surface and top for the vertical velocity. For temperature, Dirichlet conditions are used. Heating of the air by subsidence is calculated as the product of a constant linear subsidence profile (zero at the surface) and the domain averaged temperature gradient. Subsidence of momentum is not included in the current work for simplicity. However, it would be interesting to assess its potential impact on the intermittency mechanism in future work.

Simulations are initialized with constant temperature *θ*_{0} and constant velocity (*G*, 0, 0) in the *x*, *y*, and *z* directions, respectively. At the start of the simulation, the surface is cooled by 25 K after which cooling is stopped, and the simulation is continued to reach steady state. An overview of the parameters used in the VSBL case is given in Table A1.

### APPENDIX B

#### Derivation of the LSA

We consider the conservation equation of mass, the inviscid Navier–Stokes equation and the conservation equation of energy (written in temperature form) under the Boussinesq approximation in 2D:

where *u*_{i} are the velocity components in the *x* and *z* direction, *θ* is the potential temperature, *θ*_{0} is the reference temperature, *g* is the acceleration due to gravity, and *p* is the modified pressure.

We assume that our variables can be decomposed into their mean background states and a contribution due to perturbations indicated by a capital letter and a prime, respectively,

These expressions are inserted into Eq. (B1) and subsequently the mean state balance is subtracted. Additionally, products of perturbed quantities are assumed to be negligibly small and therefore removed. This results in a new set of *linearized* equations for the perturbed variables:

By taking the derivatives of Eqs. (B3b) and (B3c) with respect to *x* and *z*, respectively, adding them and applying Eq. (B3a), a Poisson equation for the pressure is obtained:

Subsequently, by taking the Laplacian (∇^{2}) of Eq. (B3c) and the *z* derivative of Eq. (B4), the pressure is eliminated. This results in a reduced set of equations for the vertical velocity perturbation and the temperature:

Next, traveling-wave solutions (complex Fourier components) are taken as ansatz, for example, for the vertical velocity (perturbation):

where *k* is the real wavenumber, $w^k\u2061(z)$ is the complex amplitude (profile) of the *k*th mode, *c*_{k} = *c*_{k,R} + *ic*_{k,I} is the phase speed, and *σ*_{k} = −*ikc*_{k} is the growth rate. A positive value of Re(*σ*_{k}) (or *c*_{k,I}) results in a growing wave mode in time indicating instability. Substitution of this ansatz and cancellation of the exponentials leads to (for each wave mode separately)

This set of equations is to be numerically solved for the unknown growth rate *σ*_{k}, and the corresponding profiles of the vertical velocity and temperature perturbations. To do this, a finite-difference approximation is used in which the amplitude profiles are discretized in *N*_{z} vertical levels [i.e., $w^k\u2061(z)$ is discretized as the vector $w^k$ of finite length *N*_{z}]. This transforms Eq. (B5) into a generalized eigenvalue problem with eigenvalue *σ*_{k} and eigenvector $[w^k,\theta ^k]T$:

in which the $A$, $B$_{11}, $B$_{12}, $B$_{21}, and $B$_{22} are block matrices of size *N*_{z} × *N*_{z} with *N*_{z} being the amount of vertical levels. These block matrices are given by

where $D$^{2} is the matrix for the finite-difference second derivatives, $I$ is the identity matrix, and **U**, **U**_{zz}, and **T**_{zz} are column vectors (size *N*_{z} × 1) of the discretized background velocity magnitude, second derivative of the velocity magnitude and the second derivative of the temperature, respectively. Note that, $B$_{12}, $B$_{21}, and $B$_{22} are diagonal matrices. The second derivatives are calculated by using a second-order central difference scheme.

For the configuration in this study, the boundary conditions are $w^k=0$ at *z* = 0 and $dw^k/dz=\u2212kw^k$ at *z* = *L*_{z} (top of the computational domain). The latter is an approximation for $w^k\u21920$ as *z* →∞ (see, e.g., Newsom and Banta 2003), which can be derived from the second-order differential equation (Taylor–Goldstein equation) resulting from elimination of the temperature perturbation and using the fact that both *d*^{2}*U*/*dz*^{2} and *d*Θ/*dz* tend to zero above the SBL. These boundary conditions for $w^k$ and its first derivative are imposed through modification of $D$^{2}. The first boundary condition does not require a change of $D$^{2}. The second is implemented by alteration of the trace element and subtrace element: $D2\u2061(Nz,Nz)=\u2061(\u22122\u22122\Delta zk)/\Delta z2$ and $D2\u2061(Nz\u22121,Nz)=2/\Delta z2$.

Figure B1 shows the growth rate as a function of the wavenumber for the case considered in the main text. The wavenumber of the unstable waves in the extended domain simulation as identified by the spectral analysis (indicated by the blue dotted line; see section 3a) is close to the global maximum given by the LSA. The difference in predicted growth rate is less than 4%. This minor discrepancy may be due to viscosity effects (not considered in the LSA), nonlinear growth or numerical approximations.

## REFERENCES

*Bound.-Layer Meteor.*

*Quart. J. Roy. Meteor. Soc.*

*J. Atmos. Sci.*

*Dyn. Atmos. Oceans*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Appl. Meteor.*

*Bound.-Layer Meteor.*

*Hydrodynamic Stability*

*J. Atmos. Sci.*

*J. Geophys. Res.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Fluid Mechanics*

*Tellus*

*Bound.-Layer Meteor.*

*Annu. Rev. Fluid Mech.*

*Bound.-Layer Meteor.*

*J. Atmos. Sci.*

*J. Fluid Mech.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Geophys. Res.*

*Bound.-Layer Meteor.*

*J. Atmos. Sci.*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*Bound.-Layer Meteor.*

*J. Atmos. Sci.*

*Bound.-Layer Meteor.*

*J. Atmos. Sci.*

*Geosci. Model Dev.*

*Bound.-Layer Meteor.*

*Quart. J. Roy. Meteor. Soc.*

*J. Adv. Model. Earth Syst.*

*J. Atmos. Sci.*

## Footnotes

Denotes content that is immediately available upon publication as open access.

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).

Workshop on Micrometeorology, D. A. Haugen, Ed., Amer. Meteor. Soc.