2.1 Applying finite-time stabilizing buoyancy flux

As discussed in the introduction, previous DNS of the Ekman layer have used either constant-flux stability or constant-temperature stability. The present study incorporates a hybrid boundary condition, i.e. a constant-flux stability for a finite time of $fT\approx 6$ and a constant-temperature stability afterwards. In other words, the applied buoyancy flux $B_{0}$ is constant over the time interval $T$ and the surface cools to a temperature $\unicode[STIX]{x1D703}_{0}(T)$ during this interval after which the surface temperature is held constant at $\unicode[STIX]{x1D703}_{0}(T)$ . The choice of $ft\approx 6$ as the interval of constant buoyancy flux is prompted by the finding of Gohari & Sarkar (Reference Gohari and Sarkar2017) that the initial adjustment of the Ekman layer to the imposed buoyancy flux occurs over approximately an inertial period. At high latitudes with large $f$ , night-time cooling can occur over as long as $ft\approx 6$ .

Different non-dimensional parameters can be used to characterize the stability imposed on the flow. The Obukhov length, $L$ , used to quantify the strength of the buoyancy flux ( $B_{0}$ is negative for a cooling flux) is defined by

where $\unicode[STIX]{x1D705}=0.41$ is the von Kármán constant. Turbulence at a height above $L$ is strongly affected by buoyancy. Using inner-layer scaling, the normalized Obukhov length becomes

while outer-layer scaling yields

where $\unicode[STIX]{x1D6FF}_{N}=u_{\ast _{N}}/f$ is the scale for the neutral Ekman boundary-layer height. The non-dimensional Obukhov length scale can also be interpreted in terms of the gradient Richardson number at the surface

where $N$ is the buoyancy frequency and $S$ is the mean vertical shear ( $Pr=1$ has been assumed). Therefore, $L^{+}$ contains stratification information. Note that, since $u_{\ast }$ is a time-dependent function in the stratified DNS cases, so are $L$ , $L^{-}$ and $L^{+}$ .

Bulk Richardson number, a stability measure well suited for constant-temperature stability, is defined as

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{0}=\unicode[STIX]{x1D703}_{\infty }-\unicode[STIX]{x1D703}_{0}$ is the difference between the temperature above the boundary layer and the surface temperature. In the present study, the surface temperature ( $\unicode[STIX]{x1D703}_{0}$ ) is a time-dependent variable during the time interval, $T$ , over which finite-time constant-flux stability is imposed. The value of $Ri_{b}$ is also initially dependent on time but it becomes constant at the end of the time interval, $T$ .

2.2 Criterion for turbulence collapse

Turbulence in the neutral boundary layer collapses initially if the imposed stability is sufficiently strong. There are several potential candidates for a turbulence-collapse criterion: $Ri_{b}$ , $L^{-}$ , and $L^{+}$ . Previous studies (Ansorge & Mellado Reference Ansorge and Mellado2014; Shah & Bou-Zeid Reference Shah and Bou-Zeid2014) have shown that the critical value of $Ri_{b}$ above which neutral Ekman flows undergo initial collapse of turbulence is a function of Reynolds number. Similarly, a $Re$ -dependent behaviour is observed for the critical value of $L^{-}$ (Flores & Riley Reference Flores and Riley2011; Gohari & Sarkar Reference Gohari and Sarkar2017). These findings are in agreement with previous studies of stratified channel flow at different Reynolds numbers, in which neither $(h/L)_{\mathit{cri}}$ ( $h$ is the half-width of channel) nor $Ri_{b}$ was found to be the control parameter of the initial turbulence collapse (Armenio & Sarkar Reference Armenio and Sarkar2002; Nieuwstadt Reference Nieuwstadt2005; Flores & Riley Reference Flores and Riley2011). Cases with the normalized Obukhov length, $L^{+}$ , as 650, 705, 683 and 714 exhibited turbulence collapse in the previous studies of Armenio & Sarkar (Reference Armenio and Sarkar2002), Nieuwstadt (Reference Nieuwstadt2005), Flores & Riley (Reference Flores and Riley2011) and Deusebio et al. (Reference Deusebio, Brethouwer, Schlatter and Lindborg2014), respectively. As explained by Flores & Riley (Reference Flores and Riley2011), $L^{+}$ is a measure of the dynamic range between the buoyancy constraint of $L$ on the large-eddy length and $\unicode[STIX]{x1D708}/u_{\ast }$ that represents the small-eddy length scale of the pre-existing near-wall turbulence. If $L^{+}$ becomes sufficiently small during the flow evolution (less than 100 according to Flores & Riley (Reference Flores and Riley2011) in a stratified channel), turbulence cannot be maintained in wall-bounded flows since buoyancy influences the entire dynamic range of turbulence scales. It is worth noting that the collapse criterion of $L^{+}=100$ was inferred by (Flores & Riley Reference Flores and Riley2011) from a stratified channel with an initial $L^{+}=683$ . The field measurements of CASES-99 campaign, a comprehensive collection of flow statistics taken for 28 nights with clear sky conditions (Poulos et al. Reference Poulos, Blumen, Fritts, Lundquist, Sun, Burns, Nappo, Banta, Newsom, Cuxart and Terradellas2002) provide further guidance. Banta et al. (Reference Banta, Mahrt, Vickers, Sun, Balsley, Pichugina and Williams2007) analysed the typical structure of five very stable nights from the CASES-99 campaign and discussed implications of turbulent exchange between the surface layer and the upper layers of the lower atmosphere. While the studied nights showed a typical $L^{+}\approx O(10^{4})$ , the lowest reported value of $L^{+}$ with turbulence was 667.

As noted in the introduction, we are interested in studying the stable regime where turbulence initially collapses. According to the above literature review, application of a cooling flux with $L_{\mathit{cri}}^{+}\approx 700$ to a neutral boundary layer is hypothesized to be sufficient for turbulence collapse and we conduct series A of the DNS with initial $L^{+}\approx 700$ . The degree of Reynolds independence of the $L_{\mathit{cri}}^{+}\approx 700$ condition for turbulence collapse is ascertained considering three Reynolds numbers, $Re_{\ast }=\{692,1112,1722\}$ . At each Reynolds number, a separate simulation is performed wherein a neutral Ekman boundary layer is allowed to evolve and reach a quasi-steady-state condition for approximately one inertial oscillation. At the beginning of the stratification phase, a constant stabilizing buoyancy flux, $B_{0}$ , corresponding to initial $L^{+}\approx 700$ is applied at the surface for a finite time of $fT\approx 6$ , and the surface temperature is held constant afterwards. We investigate the flow evolution subsequent to the imposed stratification in these three cases.

Table 1.

DNS parameters of simulated cases in series A. The case label in column 1 starts with N (neutral case) or $S$ (stable case) and has a number corresponding to the target value of friction Reynolds number, $Re_{\ast }$ . In the stable cases, a surface buoyancy flux, chosen to obtain a target $L^{+}\approx 700$ , is applied for a finite time of $fT\approx 6$ and then the surface temperature is held constant. Note that $\unicode[STIX]{x1D6FF}_{N}$ is a measured parameter calculated from averages over the last $ft=2$ of the neutrally stratified simulation at each Reynolds number. The numbers in bold signify the controlling parameters in the simulations.

2.3 Numerical model

The numerical model used to solve (2.1) and (2.2) is a modified version of that developed by Jalali, Rapaka & Sarkar (Reference Jalali, Rapaka and Sarkar2014) and Gohari & Sarkar (Reference Gohari and Sarkar2017). Major features of this solver are as follows: (i) a low storage, third-order Runge–Kutta method is used as the time marching algorithm; (ii) derivatives in the spanwise direction are computed using Fourier transforms, while the derivatives in the other two directions are computed using the central second-order finite difference method; and (iii) dealiasing of nonlinear terms is performed using the $3/2$ rule.

Figure 1.

Time-averaged statistics for neutrally stratified cases. (a) Mean velocity profile compared with the numerical studies of Shah & Bou-Zeid (Reference Shah and Bou-Zeid2014) and Shingai & Kawamura (Reference Shingai and Kawamura2004), and the logarithmic law: $M^{+}=(1/\unicode[STIX]{x1D705})\ln (z^{+})+B$ , where $M=\sqrt{\langle u\rangle ^{2}+\langle v\rangle ^{2}}$ is the mean horizontal velocity magnitude and $B=5.3$ is a constant. (b) The von Kármán constant, $\unicode[STIX]{x1D705}$ , computed from the neutrally stratified cases (see (3.1)). The typical range of von Kármán constant reported in the literature (0.36–0.41) is shown by the two horizontal dashed lines. The time averaging is performed over the last $ft\approx \unicode[STIX]{x03C0}$ .

Table 1 presents a summary of the simulated cases of series A. The stratified-flow simulations have the same initial value of $L^{+}\approx 700$ while $Re_{\ast }$ is varied among cases. An important requirement in DNS is the proper resolution of viscous effects near the surface which can extend up to $z^{+}=10$ . The viscous sublayer has 12 points in the neutral simulation with a non-dimensional vertical grid spacing, $\unicode[STIX]{x0394}z_{\mathit{min}}^{+}\leqslant 1$ . The grid is stretched moving upwards, but the grid spacing does not exceed $\unicode[STIX]{x0394}z^{+}=18$ at the top of the boundary layer in all cases. The streamwise resolution is $\unicode[STIX]{x0394}x^{+}\approx 8.7$ and the spanwise resolution is $\unicode[STIX]{x0394}y^{+}\approx 4.9$ similar to the previous DNS studies of wall-bounded turbulent flows (Flores & Riley Reference Flores and Riley2011; Ansorge & Mellado Reference Ansorge and Mellado2014; Gohari & Sarkar Reference Gohari and Sarkar2017). Boundary conditions correspond to no-slip and impermeability at the surface, and stress free and adiabatic at the top. To minimize spurious reflection of internal gravity waves in each case, a Rayleigh damping layer of 10 points is used at the top. In all the simulated cases, the domain size is similar to the one used in the recent studies of Shah & Bou-Zeid (Reference Shah and Bou-Zeid2014) and Gohari & Sarkar (Reference Gohari and Sarkar2017), which is observed to be sufficient to capture all the dynamically active turbulence scales.

4.1 Turbulence collapse and transitional Ekman layer

The stabilizing buoyancy flux leads to an initial period of collapsed turbulence, emergence of a low-level jet (LLJ) and the formation of a thermal inversion in all cases. The time evolution of normalized friction velocity in figure 4(a,c,e) shows that all the cases have an initial decrease of surface drag that is followed by an increase. The initial decay rate of $u_{\ast }/U_{\infty }$ is larger at higher Reynolds numbers with a typical drop of approximately $20\,\%$ in magnitude. Although the normalized friction velocity in the low- $Re$ case $S700$ follows the general trend of an increase after $ft\approx 4$ , its late-time behaviour is quite intermittent and the entire lower flank of the LLJ exhibits oscillations in the shear and the direction. In case $S1100$ , the LLJ oscillation near the surface is temporally less intermittent and there is an almost constant $u_{\ast }/U_{\infty }$ at a late time. The difference in the late-time behaviour of these cases is elaborated in § 4.2. In case $S1700$ with the higher $Re_{\ast }$ , the normalized friction velocity increases during the LLJ emergence; however, in contrast to the other cases, $u_{\ast }/U_{\infty }$ exhibits a trend of consistent decrease at late time with an oscillatory modulation at the inertial period of $ft=2\unicode[STIX]{x03C0}$ .

Figure 4.

The time evolution of: (a,c,e) normalized friction velocity ( $u_{\ast }/U_{\infty }$ ); and (b,d,f) integrated turbulent kinetic energy ( $E/(\unicode[STIX]{x1D6FF}_{N}u_{\ast _{N}}^{2})$ ), for stably stratified cases at $Re_{\ast }=692$ (a,b), $Re_{\ast }=1112$ (c,d) and $Re_{\ast }=1722$ (e,f). The finite-time buoyancy flux is imposed on a neutral Ekman layer at $ft=0$ for one inertial oscillation ( $ft\approx 6$ ). The letters correspond to the panels of figure 7.

There is a striking influence of Reynolds number on the time evolution of the integrated TKE (figure 4 b,d,f) defined by

In the higher- $Re$ case $S1700$ shown in figure 4(f), the TKE does not recover at long time while it does in the lower- $Re$ simulations. In all cases, the normalized $E$ initially decreases to nearly zero during a ‘decay period’, indicating a turbulence collapse within the interval of constant buoyancy flux. This confirms that the choice of $L_{\mathit{cri}}^{+}\approx 700$ results in collapsed turbulence during the initial transient, independent of Reynolds number.

The TKE recovery in case $S700$ (figure 4 b) is not monotonic and there is noticeable temporal intermittency in the TKE. The period of large-amplitude oscillations of the integrated TKE corresponds to the inertial period, $ft=2\unicode[STIX]{x03C0}$ . The oscillations are less pronounced in case $S1100$ with an intermediate $Re_{\ast }$ .

As noted above, case $S1700$ exhibits a decay period with collapsed turbulence, however, unlike case $S1100$ , collapse is not followed by a recovery period. This difference between the two cases is due to the stronger bulk stability of case $S1700$ that is an indirect effect of the change in $Re_{\ast }$ among cases. All cases have a similar initial $L^{+}$ but the larger value of $Re_{\ast }$ in case $S1700$ requires a larger surface buoyancy flux to achieve the same initial $L^{+}$ as will be discussed in § 4.3 and, therefore, the value of $Ri_{b}(t)$ (see figure 10 and its discussion) in case $S1700$ is larger at the end of the cooling period relative to the other cases. As a consequence, $Ri_{g}$ in the LLJ region becomes sufficiently large in case $S1700$ and prevents the recovery of turbulence, as described below.

Figure 5(a,b) shows the time–height plot of normalized TKE for case $S1100$ and $S1700$ , respectively. In case $S1100$ , there is virtually no near-surface turbulence during $2.5<ft<3.5$ . However, the region between $0.05\leqslant z^{-}\leqslant 0.1$ has residual TKE during this period. After $ft\approx 3.5$ , near-surface turbulence builds up. The vertical transport of fluctuating energy from the $0.05\leqslant z^{-}$ region into the near-surface layers ( $z^{-}\leqslant 0.02$ ) triggers the surface turbulence and leads to eventual recovery as was found by Gohari & Sarkar (Reference Gohari and Sarkar2017) to occur during the turbulence ‘rebirth’ episodes in the Ekman layer subject to constant-flux stability. For case $S1100$ , this vertical exchange between outer and inner layers occurs inside a LLJ region which is dynamically unstable. The entire region between $z=0$ and the $z$ -location of the square symbols on the dashed profiles of figure 5(d) is subcritical with $Ri_{g}<0.25$ . Therefore, the LLJ in case $S1100$ is receptive to fluctuations transported into that region and rebirth of turbulence is possible. Note that the $z$ -location of each square symbol marking the upper boundary of the subcritical region is obtained from the intersection (the first one traversing upward from $z=0$ ) of the $Ri_{g}(z^{-})$ profile in figure 5(c) of the corresponding case with the vertical line at $Ri_{g}=0.25$ .

Figure 5.

Comparison of early time evolution between case $S1700$ and case $S1100$ : (a,b) time–height plot of normalized turbulent kinetic energy ( $e/u_{\ast _{N}}^{2}$ ) in case $S1100$ and $S1700$ , respectively; (c) gradient Richardson number; (d) mean velocity magnitude. The upper boundary of the near-surface layer with $Ri_{g}\leqslant 0.25$ is marked with circles in case $S1700$ (solid lines) and with squares in case $S1100$ (dashed lines).

In contrast with case $S1100$ , the sheared lower and upper flanks of the LLJ in case $S1700$ have a supercritical value of $Ri_{g}>0.25$ that prevents instability even if there is a vertical transport of finite-amplitude fluctuations into the LLJ. Consequently, the near-surface turbulence cannot recover in this case. Although the LLJ in the velocity profile (black solid line of figure 5 d) has strong shear, the stratification is also strong so as to keep $Ri_{g}$ above 0.25. Note that the circles on the solid lines of figure 5(d) that mark the upper boundary of the subcritical near-surface region are lower than the nose of the LLJ and their offset from the nose progressively increases with increasing $ft$ . Thus, the near-surface layer becomes progressively decoupled from the LLJ mean shear aloft and laminar flow persists.

Figure 6.

Quantification of the initial collapse of turbulence. (a) The early time evolution of integrated TKE. Each dashed line is a linear fit to the initial decrease of TKE and the adjacent number denotes the slopes of the linear fit. (b) The slopes are plotted as a function of $Re_{\ast }$ along with the line that is the best fit to the data.

It is of interest to quantify the dependence of the time interval over which TKE collapses on the parameters $L^{+}$ and $Re_{\ast }$ in the problem. Figure 6(a) shows that $E$ decreases more rapidly as $Re_{\ast }$ increases. Furthermore, the initial decrease of $E$ is linear. The slope $\text{d}E/\text{d}t$ of the linear phase can be used to define a characteristic turbulence-collapse time scale, $t_{c}$ , as

since the initial value of $E$ (the height-integrated TKE) is proportional to $u_{\ast _{N}}^{2}\unicode[STIX]{x1D6FF}_{N}$ . The stable buoyancy flux introduces a length scale, the Obukhov length $L$ and the initial turbulent fluctuations introduce a velocity scale $u_{\ast _{N}}$ so that a natural time scale for the buoyancy effect is $t_{b}=L/u_{\ast _{N}}$ (Flores & Riley Reference Flores and Riley2011). Now,

where the last of the above equations follows from the definitions of $L^{+}$ by (2.6) and $Re_{\ast }$ by (2.4c ). We hypothesize that $1/ft_{c}$ is proportional to the inverse normalized buoyancy time scale $1/ft_{b}$ so that

The DNS data presented in figure 6(a) are obtained to compute the left-hand side of (4.4) and, as shown by figure 6(b), the magnitude of the decay rate is found to increase linearly with $Re_{\ast }$ in agreement with the prediction of (4.4) for constant $L^{+}$ . The coefficient of proportionality is found to be $1/4.06$ by linear regression of the DNS data of figure 6(b). Thus, a model for $t_{c}$ defined by (4.2) is

The corresponding values of $t_{c,model}=4.06L/u_{\ast _{N}}$ are denoted in figure 6(a) with circles. It is evident that $4.06L/u_{\ast _{N}}$ is a reasonable measure of the time interval required for turbulence to collapse in the simulation. It is worth noting that the DNS data do not show a sharp minimum of $E$ so that a precise estimate of the collapse time by the criterion of minimum $E$ is not possible.

Visualizations of normalized streamwise velocity fluctuation, $u^{\prime }/u_{\ast _{N}}$ , near the surface at $z^{+}=15$ (a–c) and in the outer layer at $z^{-}=0.2$ (d–f), are presented in figure 7 in order to illustrate spatio-temporal intermittency of stratified Ekman flows. These snapshots are taken during the initial decay time period and the corresponding times are indicated using dots in figure 4(b,d,f). Figure 7 indicates that the collapse process does not proceed uniformly in space near the surface, instead, turbulence persists in localized streaks contained within a region where fluctuations have subsided. Thus, there is spatio-temporal intermittency in the near-surface flow wherein turbulence coexists with regions of nearly laminar flow during the temporal decrease of TKE. In contrast, the intermittency is observed to be less strong in the outer layer where turbulence decays in an approximately homogeneous manner , i.e. the turbulent fluctuations in the outer layer are reduced with respect to the neutral case but are not preferentially distributed as streaks. It should be noted that the near-surface streaks are continually evolving and their growth or decay depends on the competition between destabilizing vertical mean shear and stabilizing stratification. For case $S700$ and case $S1100$ , these structures eventually develop into a spatially extended turbulent state that expands to fill the whole flow domain, however, these structures peter out in case $S1700$ and TKE becomes negligible as previously discussed.

Figure 7.

Spatial intermittency of near-surface fluctuations in the stratified Ekman layer. Streamwise velocity fluctuations, $u^{\prime }/u_{\ast _{N}}$ , on a horizontal plane: (a–c) correspond to a near-surface location, $z^{+}=15$ , and (d–f) to an outer-layer location, $z^{-}=0.2$ . (a,d) Case $S700$ at $ft=2.22$ ; (b,e) case $S1100$ at $ft=1.87$ ; and (c,f) case $S1700$ at $ft=1.21$ . The corresponding times are indicated with dots and letters in figure 4. Spatial intermittency refers to near-ground turbulence taking the form of isolated, organized turbulent patches within a background of negligible fluctuations.

Figure 8.

The time evolution of: (a) normalized Obukhov scale ( $L^{+}$ ); (b) buoyancy Reynolds number ( $Re_{b}$ ). Case $S1700$ has almost zero turbulent dissipation and $Re_{b}$ after $ft\approx 2$ .

Figure 9.

Profiles of statistics at long time against normalized height: (a) normalized velocity; (b) normalized potential temperature ( $u_{\ast _{N}}(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{\infty })/q_{0}$ ); (c) gradient Richardson number ( $Ri_{g}$ ); (d) normalized variance of temperature fluctuations ( $u_{\ast _{N}}\unicode[STIX]{x1D703}_{\mathit{rms}}/q_{0}$ ); (e) vertical buoyancy flux ( $\langle w^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle /q_{0}$ ); (f) normalized TKE ( $e/u_{\ast _{N}}$ ). The statistics are averaged during the last $ft\approx 2\unicode[STIX]{x03C0}$ .

Figure 8(a) illustrates the time series of normalized Obukhov scale. The value of $L^{+}$ , initially ${\approx}700$ for each case, reaches a minimum value ranging between 200 and 400, similar to the minimum value that is previously reported by Flores & Riley (Reference Flores and Riley2011), Deusebio et al. (Reference Deusebio, Brethouwer, Schlatter and Lindborg2014, Reference Deusebio, Caulfield and Taylor2015). At the point of minimum $L^{+}$ , the near-surface fluctuations exhibit spatial intermittency. $L^{+}$ eventually increases to values larger than the initial value. Case $S1700$ , although having larger $L^{+}$ at later times, is non-turbulent. This indicates that, although $L^{+}$ can provide guidance as to whether there is an initial turbulence collapse, it does not help predict whether turbulence remains collapsed at long time. The buoyancy Reynolds number, $Re_{b}$ , that has been used to distinguish between strong and weak turbulence in stratified flows (Smyth & Moum Reference Smyth and Moum2000; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Deusebio et al. Reference Deusebio, Caulfield and Taylor2015) is defined as follows,

where $\unicode[STIX]{x1D700}$ is the viscous dissipation. An overall value of $N^{2}=g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{0}/\unicode[STIX]{x1D6FF}_{N}$ is assumed to obtain the second equation in (4.6) and the surface value of $\unicode[STIX]{x1D700}$ is used in the computation of $Re_{b}$ . Figure 8(b) depicts the time evolution of $Re_{b}$ . By correlating $Re_{b}(t)$ to the fluctuations (TKE and spatial near-surface structures) it is possible to identify critical values of $Re_{b}$ that distinguish different flow states in stratified Ekman layers: fully turbulent when $Re_{b}\geqslant 180$ , intermittently turbulent when $90\leqslant Re_{b}\leqslant 180$ and non-turbulent when $Re_{b}\leqslant 90$ . These critical values are similar to the ranges previously reported by Deusebio et al. (Reference Deusebio, Caulfield and Taylor2015) for the boundary between intermittent and fully developed turbulence in stratified plane Couette flow ( $30\leqslant Re_{b}\leqslant 650$ ). It is worth noting that, although these values of $Re_{b}$ are helpful in classifying the state of fluctuations in Ekman flows given the value of dissipation rate, $Re_{b}$ is a derived turbulence-based quantity. We will explore the suitability of a parameter based on mean flow quantities as a criterion for turbulence recovery in § 4.3.

4.2 Stratified Ekman layer at late time

In the previous section, it was shown that the Ekman flow response to a constant finite-time stabilizing flux with initial $L^{+}\approx 700$ includes a decay period with collapsed turbulence and significant changes in the boundary-layer structure. During the decay period, the near-surface region has turbulent streaks in an otherwise quiescent flow. Turbulence recovers in the lower- $Re_{\ast }$ cases $S700$ and $S1100$ but the flow remains non-turbulent at the higher $Re_{\ast }=1700$ .

The flow state of the Ekman layer at a late time is quantified . Several statistics are time averaged over the last $ft\approx \unicode[STIX]{x03C0}$ of the simulations ( $27\leqslant ft\leqslant 30$ ) and the profiles are shown in figure 9. The velocity profiles in case $S1700$ (figure 9 a) exhibit the strongest LLJ with an overshoot in the streamwise velocity. The temporal evolution of velocity magnitude shows an LLJ between $0.05<z^{-}<0.15$ that oscillates in time (not shown), quite different from the quasi-steady behaviour of the other two stable cases and the neutral reference case. Although all cases develop stable stratification in the near-surface layer, the mean temperature profile (figure 9 b) exhibits quantitative and qualitative differences in the high- $Re$ case. Case $S1700$ develops the largest stable temperature gradient in the near-surface layer and, furthermore, the temperature varies linearly with height before smoothly merging to the constant-temperature region above. The other cases also have a near-surface layer that has stratification but there is turbulent mixing so that the maximum temperature gradient is not within the near-surface layer but resides above in the capping layer. Despite its stronger LLJ and shear, the strongly stratified inversion in the near-surface layer of case $S1700$ leads to a stable region with $Ri_{g}>0.25$ that extends from below the LLJ nose (figure 9 c) and is able to decouple the near-surface layer from aloft. In the other two cases, the region with $Ri_{g}>0.25$ starts higher up and, therefore, the LLJ has a subcritical $Ri_{g}$ allowing for shear instability. The root-mean-square of the temperature fluctuations (figure 9 d) peaks in the thermocline at the location where the temperature gradient also reaches its maximum.

The vertical distribution of TKE is shown in figure 9(f). The TKE profile of case $S1100$ has a notably two-layer structure with two local maxima, one below and the other above the LLJ nose level, and a local minimum of TKE at the LLJ nose height. This result is similar to the previous finding of Cuxart & Jiménez (Reference Cuxart and Jiménez2007). It is worth noting that the numerical study of Cuxart & Jiménez (Reference Cuxart and Jiménez2007) was initialized with the velocity and temperature profiles of an observed LLJ in contrast to the present study where the LLJ inherently forms as a response of the Ekman flow to the imposed stabilizing buoyancy flux. Finally, the vertical TKE distribution of case $S1700$ shows a non-turbulent near-surface layer with a small amount of fluctuating energy located in the outer region with oscillating LLJ shear ( $0.05<z^{-}<0.15$ ).

We investigated the robustness of the case $S1700$ result of a relaminarized boundary layer by conducting a new simulation where random perturbations were added to relaminarized state and examined if the boundary layer remains non-turbulent during its subsequent development. Initial velocity perturbations with fluctuation energy, $E/E_{\mathit{neutral}}=0.93$ , and concentrated in a region close to the surface, i.e. $z^{-}<0.1$ were added to the flow at $ft_{0}\approx 27$ . These fluctuations were generated using the model spectrum,

where $k$ is the spanwise wavenumber, $\unicode[STIX]{x1D700}$ is the dissipation rate, $\unicode[STIX]{x1D702}$ is the Kolmogorov scale and $L=0.5\unicode[STIX]{x1D6FF}_{N}$ . The details of this model spectrum and the non-dimensional functions, $f_{L}$ and $f_{\unicode[STIX]{x1D702}}$ , can be found in § 6.5.3 of Pope (Reference Pope2000). The magnitude of $L$ , $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D700}$ are taken from the neutral Ekman layer. The time evolution of TKE showed that the initially seeded perturbations were continuously damped and were eventually extinguished by $f(t-t_{0})\approx 0.5$ .

4.3 Criterion for recovery

As discussed earlier in § 4.1, although the initial $L^{+}\approx 700$ provides a stabilizing flux that is sufficiently strong to cause transient collapse of turbulence, it does not determine the final state of the boundary layer. Turbulence recovers in the lower- $Re_{\ast }$ cases $S700$ and $S1100$ but the flow remains non-turbulent in the higher $Re_{\ast }=1700$ case, although the late-time values of $L^{+}$ are similar among the cases. From the results of §§ 4.1 and 4.2, it appears that turbulence recovery occurs in the lower- $Re$ cases since the $Ri_{g}(z)$ profiles have subcritical value in the near-surface layer that allows shear instability to develop. In this section, we investigate how $Ri_{b}$ as an alternate choice of overall stratification can provide guidance as to whether turbulence recovers after its initial collapse, by providing a theoretical analysis for the evolution of $Ri_{b}$ that shows an increased $Ri_{b}$ when $Re_{\ast }$ increases at a fixed $L^{+}$ , a comparison of such analysis with DNS results, and providing a relation between $Ri_{b}$ and average $Ri_{g}$ .

Consider the equations governing horizontally averaged quantities denoted by $\langle .\rangle$ . Continuity and the impermeability condition leads to the result that $\langle u_{3}\rangle (z)$ is identically zero. Averaging of (2.2) over horizontal planes leads to

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{\infty }-\unicode[STIX]{x1D703}$ and superscript $^{\prime }$ denotes fluctuation with respect to the average. After using the following boundary conditions:

the height-integrated form of (4.8) provides the energy budget equation

Note that $q_{0}=-\unicode[STIX]{x1D706}\unicode[STIX]{x2202}_{x_{3}}\langle \unicode[STIX]{x1D703}_{0}\rangle =\unicode[STIX]{x1D706}\unicode[STIX]{x2202}_{x_{3}}\langle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{0}\rangle$ is a constant during the application of the constant buoyancy flux in series A. Therefore, integration of (4.10) in time yields

The left-hand side of (4.11) can be approximated as $\langle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{0}\rangle \unicode[STIX]{x1D6FF}_{N}/C_{1}$ , where $C_{1}$ is a coefficient ( $C_{1}=2$ , if $\langle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\rangle$ varies linearly with height similar to case $S1700$ ). Therefore, the equation given in (4.11) becomes

After substituting (4.12) for the surface temperature into the definition of the bulk Richardson number, (2.9), and using (2.6), we obtain the following expression that describes the initial time evolution of $Ri_{b}(t)$ :

The non-dimensional parameters in the problem are the initial Reynolds number ( $Re_{\ast }$ or $Re_{D}=\sqrt{2U_{\infty }/\unicode[STIX]{x1D708}f}$ ), the constant buoyancy flux ( $L^{+}$ ) and the time interval ( $fT$ ) over which the buoyancy flux is applied. In order to relate $Ri_{b}(t=T)$ to these non-dimensional parameters, (4.13) is rewritten as

where

The second equality of (4.14) is obtained by using the correlation $u_{\ast _{N}}/U_{\infty }\approx 1/\ln Re_{D}$ (Coleman et al. Reference Coleman, Ferziger and Spalart1990).

Figure 10.

The time evolution of bulk Richardson number, $Ri_{b}$ , in series A. The value of $Ri_{b}(ft=6)$ increases with $Re_{\ast }$ . Note that in series B to be discussed later in § 5.2, each case has a constant-temperature boundary condition imposed according to the value of $Ri_{b}(fT\approx 6)$ of the corresponding series A case.

Table 2.

Summary of the parameters in (4.14). The value of $C_{2}$ coefficient is calculated using the DNS data, i.e. figure 10 and table 1, showing a general value $C_{2}\approx 0.3{-}0.4$ .

The time-evolution relation given in (4.14) renders three scenarios under which the bulk stability of the Ekman layer measured by $Ri_{b}(T)$ becomes large when one of the three independent parameters $\{L^{+},fT,Re_{D}\text{ or }Re_{\ast }\}$ is varied, keeping the other two constant. First, if a stronger stabilizing buoyancy flux (reduced $L^{+}$ ) is imposed at the surface. Second, if the stabilizing buoyancy flux is imposed for a longer time interval, i.e. larger value of $fT$ . Third, if a stabilizing buoyancy flux is imposed at a larger Reynolds number, i.e. $Re_{D}$ or equivalently $Re_{\ast }$ is increased (the last term in 4.14 monotonically increases with $Re_{D}$ ). It is the third scenario that applies to DNS series A where $Re_{\ast }$ is increased keeping $L^{+}$ and $fT$ constant. From a physical viewpoint, an increased Reynolds number reduces the viscous scale of the Ekman layer ( $O(\unicode[STIX]{x1D708}/u_{\ast })$ ) and therefore, imposing a stabilizing buoyancy flux based on a constant $L^{+}$ (ratio of the Obukhov scale to the viscous scale) yields a reduced buoyancy length scale ( $L$ ) and a stronger stratification during the flow evolution. In agreement with this analysis, figure 10 shows that the value of $Ri_{b}(T)$ at the end of the constant-flux interval of $fT\approx 6$ increases with increasing value of initial $Re_{\ast }$ in series A.

The validity of equation (4.14) for predicting $Ri_{b}(fT)$ is examined against DNS results. Table 2 summarizes the independent parameters (initial $L^{+}$ and $Re_{D}$ ) for the conducted DNS simulations in series A, as well as the $C_{2}$ calculated after inserting $Ri_{b}(ft=6)$ and the values of the independent parameters in equation (4.14). The values of $C_{2}\approx 0.3{-}0.4$ vary in a narrow range, showing that equation (4.14) with a suitable choice of $C_{2}$ , i.e. 0.35, can be used to estimate $Ri_{b}(fT)$ prior to the simulation.

An increase in $Ri_{b}$ implies a larger value of average $Ri_{g}$ as shown below. A vertically averaged value of $Ri_{g}$ is introduced as

where,

Therefore, (4.16) becomes

Hence, a vertically averaged gradient Richardson number is proportional to $Ri_{b}$ , despite the boundary values of $Ri_{g}=(\unicode[STIX]{x1D705}L^{+})^{-1}$ at the surface and $Ri_{g}=0$ at the top of the layer being constant among the cases of series A. Therefore, an increased value of bulk Richardson number, $Ri_{b}$ , increases the average $Ri_{g}$ in the Ekman layer (see (4.14) and (4.18)), indicating that there is a greater propensity for supercritical $Ri_{g}$ below the LLJ nose (similar to the results of case $S1700$ presented in figure 9 a), which prevents instability growth inside the near-surface layer.