Direct numerical simulation of intermittent turbulence under stably stratified conditions

In this paper, we simulate intermittent turbulence (also known as bursting events) in stably stratified openchannel flows using direct numerical simulation. Clear signatures of this intriguing phenomenon are observed for a range of stabilities. However, the spatiotemporal characteristics of intermittency are found to be strongly stability dependent. In general, the bursting events are much more frequent near the bottom wall than in the upper-channel region. A steady coexistence of laminar and turbulent flows is detected at various horizontal planes in very stable cases. This spatially intermittent pattern is found to propagate downstream and strongly correlate with the temporal evolution of intermittency. Lastly, a long standing hypothesis by Blackadar (1979), i.e., the strong connection between local stability and intermittent turbulence, is corroborated by this modeling study.


Introduction
The term "intermittency" has different connotations in various scientific communities.Even among the turbulence researchers, it does not have a unique definition.Tsinober (2001) attempted to distinguish between the two most common usages of intermittency in a very succinct manner: The term intermittency is used in two distinct (but not independent) aspects of turbulent flows.The first one is the so-called external intermittency.It is associated with what is called here partially turbulent flows, specifically with the strongly irregular and convoluted structure and random movement of the "boundary" between the turbulent and nonturbulent fluid.The second aspect is the so-called small scale, internal or intrinsic intermittency.It is usually associated with the tendency to spatial and temporal localization of the "fine" or small scale structure(s) of turbulent flows.
There is extensive literature on small-scale intermittency.A wide range of conceptual models (e.g., multifractal model) have been proposed to explain various traits of intermittency.Chapter 8 of Frisch (1995) provides an excellent review on this topic.We would like to note that the present study does not contribute to this research arena.Instead, it focuses on the global (a.k.a.external) intermittency.Townsend (1948) was one of the first to report global intermittency in a laboratory experiment.Its existence was also well documented in the atmospheric boundary layer (ABL) literature (e.g., Mahrt, 1989;Sun et al., 2002Sun et al., , 2004;;Nakamura and Mahrt, 2005).Global intermittent turbulence (a.k.a.turbulent bursting events) is usually characterized as alternately quiescent and bursty portions of an observed time series, representing laminar and turbulent states of the ABL flow, respectively.This type of non-stationary time series has been widely observed in geographically and climatologically diverse regions around the world.For example, they were observed in field experiments in Kansas, USA (Nakamura and Mahrt, 2005), Cabauw, the Netherlands (Holtslag et al., 2013), and Summit, Greenland (Cullen et al., 2007) (see Fig. 1a, b, c).In addition, intermittent turbulence was generated in wind-tunnel experiments under idealized settings (Ohya et al., 2008) (see Fig. 1d).From Fig. 1, one appreciates the fact that intermittent turbulence is a truly multiscale (bursting duration ranging from seconds to hours) and a dynamically complicated (and perhaps complex) phenomenon.It is also known in the literature that this phenomenon por-  (Nakamura and Mahrt, 2005).(b) Cabauw, the Netherlands.(c) Summit, Greenland (Cullen et al., 2007).(d) Wind-tunnel experiment (Ohya et al., 2008).
Quantification of spatiotemporal characteristics of intermittent turbulence is of considerable importance from practical standpoints.For example, it has been shown that the turbulent bursts can cause an unusually high concentration of surface layer ozone by transporting it from higher altitudes in the ABL (Salmond and McKendry, 2005).It has also been speculated that these events might be causing excessive fatigue damage to wind turbines located in the US Great Plains (Park et al., 2013).Furthermore, these episodic events have been associated with the degradation of the performance of astronomical telescopes as they tend to generate excessive optical turbulence (Lawrence et al., 2004).
Given its ubiquitous occurrence and practical importance, much effort has been devoted to investigate the characteristics and generation mechanism of intermittent turbulence.In the recent years, various mesoscale atmospheric phenomena, including (but not limited to) low-level jets (LLJ) (Sun et al., 2002), density currents (Sun et al., 2002), and gravity waves (Sun et al., 2004), were identified as possible triggers for turbulent bursts.Even though these physical processes are extremely relevant for atmospheric intermittent turbulence generation, in this study we would like to address a more fundamental question: is intermittent turbulence an inherent property of stably stratified boundary layer flows?That is to say, we would like to find out if this phenomenon can be simulated in the absence of any external forcings (e.g., topography, Coriolis force and associated inertial oscillation, mesoscale atmospheric phenomena).In some way, our investigation can be perceived as a numerical testing of an age-old hypothesis by Blackadar (1979).
Several decades ago, Blackadar (1979) associated the occurrence of intermittent turbulence with the modulation of a gradient Richardson number (a stability parameter reflecting the balance between buoyancy and shear).He hypothesized that under stably stratified conditions, a turbulent flow becomes laminarized when the buoyancy forcing strongly dominates over shear causing the gradient Richardson number to exceed a critical value.The resulting reduction in turbulent mixing leads to flow acceleration and, in turn, triggers a turbulent burst.In nature, this decaying and bursting process can occur repeatedly generating turbulent intermittency (Businger, 1973;Stull, 1988;Mahrt, 1999).To the best of our knowledge, this hypothesis by Blackadar (1979) has not been confirmed using high-fidelity modeling approaches such as direct numerical simulation (DNS) or large-eddy simulation (LES).However, the wind-tunnel experiments by Ohya et al. (2008) seem to corroborate this hypothesis.The signatures of bursting events are very clear in their data (see Fig. 1d).However, Ohya et al. ( 2008) utilized a LLJ-type velocity profile to mimic typical nocturnal ABL flows over land.It is possible that this velocity profile played a significant role in triggering the intermittent events.To avoid any ambiguity, in the present study we will refrain from using any artificial profiles or boundary/forcing conditions in our numerical simulations.
At present, LES is one of the most efficient computational techniques available for high Reynolds number turbulent flow simulations.However, its applicability to simulate strongly stratified flows has remained an unsettled issue (Basu and Porté-Agel, 2006;Flores and Riley, 2011).Furthermore, it is well known in the literature that the LESbased flow characteristics are strongly dependent on the employed subgrid-scale (SGS) model, spatial resolution, filter type, etc.To circumvent these ad hoc numerical influences, in this study we opt to use a viable alternative: the DNS approach.In this approach, the Navier-Stokes (N-S) equations are solved without any averaging, any filtering, or any other approximations for turbulence.Owing to its precise representation of the N-S equations, the DNS-generated data sets are considered to be of extremely high fidelity and are utilized to provide a better physical understanding of various types of turbulent flows (see the review by Moin and Mahesh, 1998).Unfortunately, due to exorbitant computational costs, the DNS approach is limited to low Reynolds number flows.Despite this unavoidable limitation, we are hopeful that the DNS approach will enable us to answer our posed question in an unequivocal manner.
During the past decades, several DNS studies focused on the turbulence laminarization problem under stably stratified conditions (Garg et al., 2000;Barnard, 2000;Iida et al., 2002;Nieuwstadt, 2005;Flores and Riley, 2011;García-Villalba and del Álamo, 2011;Brethouwer et al., 2012;Ansorge and Mellado, 2014).Some of these studies discussed the transient nature of the laminarization phenomenon; a few of them documented steady coexistence of laminar and turbulent flows; and a handful of these studies proposed certain criteria for turbulence decay.Below, we briefly summarize some of these relevant works as the present study borrows scientific ideas and numerical strategies from them.Nieuwstadt (2005) investigated the maximum level of stratification that a turbulent flow can survive in an open channel.He proposed a stability parameter (h/L N ) to probe this problem; it was found that turbulence intensity decayed rapidly with h/L N > 1.25.Iida et al. (2002) studied the turbulence suppression in stratified channel flows, and it was indicated that, with increasing friction Richardson number (Ri τ = ρgh/ρ 0 u 2 τ with u τ the friction velocity), the flow became laminar on one wall while it was still turbulent on the other wall.Partial laminarization phenomenon was also observed in Flores and Riley (2011); however, the laminarization was transient due to the limited computational domain size.Domain size sensitivity was later discussed by García-Villalba and del Álamo (2011).The authors found that the turbulence was non-stationary with a domain size of 4π h × 2h × 2π h and Ri τ > 60.With increasing the domain size to 18π h × 2h × 8π h, the turbulence became stationary even with a relatively high friction Richardson number (Ri τ = 480), and turbulent-laminar patches were also observed in the channel.The coexistence of laminar and turbulent flows in a stably stratified open channel was further studied by Brethouwer et al. (2012).The authors found that the laminar and turbulent flows can stably coexist within a wide range of Reynolds number and Richardson number.Intriguing spatial intermittency was also observed by Ansorge and Mellado (2014) in their DNS study of stably stratified Ekman layer.Similar to Nieuwstadt (2005) and Flores and Riley (2011), a transient process of turbulence collapse and recovery was generated by imposing a suddenly cooled surface to a neutrally stratified flow field.Spatial characteristic of intermittent turbulence was then analyzed in detail during this transient process.However, whether the spatial pattern of intermittency can persist in the long term was not discussed in the paper.Barnard (2000) performed a DNS study of intermittent turbulence in a stratified Ekman layer, and intermittent time series of vertical velocity and heat flux were documented.He proposed that the intermittent turbulence originated from the Ekman instability that lifted colder air over warmer air.Since this proposed mechanism is modulated by inertial oscillation, it is not effective in the absence of the Coriolis force.
Since the coexistence of laminar and turbulent flows has been observed in the papers by Flores and Riley (2011), García-Villalba and del Álamo (2011), Brethouwer et al. (2012), and Ansorge and Mellado (2014), we speculate that the time variation of velocity, temperature, and fluxes in a channel should manifest temporal intermittency.Unfortunately, no such results were documented in those studies.We would like to emphasize that the field experimental studies (Sun et al., 2002(Sun et al., , 2004;;Nakamura and Mahrt, 2005;Cullen et al., 2007;Coulter and Doran, 2002) almost always elaborate on temporal intermittency.Thus, there is a need for documenting simulated temporal intermittency results in order to make the DNS studies relevant to the ABL research.In the present study, we fill that void by presenting time series and time height plots of several variables.In addition, we quantify the spatiotemporal characteristics of intermittent turbulence in stably stratified channel flows.Most importantly, we investigate (i) if intermittent turbulence can be generated in the absence of any external forcings, (ii) if our DNS results can corroborate Blackadar's hypothesis, i.e., the strong correlation between local stability and turbulence intermittency.
The structure of this paper is as follows: in Sect.2, we introduce the numerical method and summarize the computational configurations and modeling details.Spatiotemporal characteristics of intermittent turbulence are documented in Sect.3. In this section, we also discuss the relationship between local stability and intermittent turbulence, as well as the sensitivity to Reynolds number.In Sect.4, the conclusions and future directions are stated.Finally, verification for the DNS solver is performed in Appendix A.

Method
In this paper, the incompressible turbulent flow in a plane open channel is studied under stably stratified conditions with three friction Reynolds numbers (Re τ = 180, 395, and 550).The friction Reynolds number, Re τ = u τ h/ν, is based on the open-channel height (h) and the friction velocity u τ = √ τ w /ρ 0 , where ν is kinematic viscosity, τ w is wall shear stress, and ρ 0 is reference density.Several stratification levels are considered by imposing different friction Richardson numbers Ri τ = ρgh/ρ 0 u 2 τ , where ρ is the density difference between the top and the bottom walls, and g denotes gravitational acceleration.In this study, seven simulations are performed with a wide range of friction Richardson numbers (cases S0-S2880 in Table 1).It is important to note that the open-channel flow is equivalent to a y symmetrized (closed-channel) plane Poiseuille flow; however, the definitions of the friction Richardson number in these two cases are slightly different.The conventional definition for Ri τ in closed-channel flows (Armenio and Sarkar, 2002;García-Villalba and del Álamo, 2011) is 2 times larger than the definition used here for open-channel flows.In other words, the stratification level with Ri τ = 240 in a closed channel is similar to that in open-channel flows with Ri τ = 120.

P. He and S. Basu: DNS of intermittent turbulence
According to García-Villalba and del Álamo (2011), a large domain is needed to simulate stably stratified flows with high friction Richardson numbers.For the current study, the domain size is chosen to be L x × L y × L z = 8πh × h × 4π h.Here, L x , L y , L z are the domain lengths in the streamwise (x), vertical (y), and spanwise (z) directions, respectively (see Fig. 2).The domain size is shown to be sufficient to accommodate the spatiotemporal evolution of intermittent turbulence (please refer to Appendix B).The grid distributions in the streamwise and spanwise directions are uniform.To fully resolve turbulence in near-wall regions, the first grid level near the bottom wall is located at y + = 0.1, and 15 grids are set within y + < 10, where y + is the normalized grid distance from the bottom wall, y + = yu τ /ν.A summary of the simulation configurations can be found in Table 1.
Under the Boussinesq approximation, the governing Navier-Stokes and the scalar transport equations are presented in Eqs. ( 1 The DNS solver used here is based on an open-source CFD package called OpenFOAM.OpenFOAM has been well verified for a number of scientific and engineering applications (Flores et al., 2013;Hoang et al., 2013).The governing equations are discretized using the finite-volume method.The spatial and temporal derivatives are discretized with the second-order central scheme and implicit second-order backward scheme, respectively (OpenFOAM, 2015).The pressure equation is solved using pressure implicit with splitting of operators (PISO) algorithm (Issa, 1986).Periodic boundary conditions are prescribed along the streamwise and spanwise boundaries.A non-slip boundary condition (u * , v * , w * = 0) is specified at the bottom wall (y * = 0), whereas at the top wall (y * = 1), a free-slip boundary condition is imposed (∂u * /∂y * = 0, v * = 0, ∂w * /∂y * = 0).Here u * , v * , and w * are streamwise, vertical, and spanwise normalized velocity components, respectively.Fixed normalized temperatures (T * top = 0, T * bot = −1) are also prescribed along the top and the bottom walls.For the neutrally stratified case (Ri τ = 0), the buoyancy term (the third term on the right- hand side of Eq. 2) is removed from the momentum equation.Under this scenario, the temperature equation actually represents the transport of a passive scalar.
In order to obtain initial conditions for the simulations, coarse-resolution cases are run for about 200 nondimensional time units (i.e., t u τ / h ≈ 200) to obtain quasistationary flow fields.Then the flow variables are interpolated to fine-resolution initial fields.For the cases with Re τ = 180, the simulations are conducted for 50 non-dimensional time to obtain the mean flow and turbulence statistics profiles, as well as time series results.For Re τ = 395 and 550, the simulations are run for 10 non-dimensional time.
According to Kim et al. (1987) and Spalart (1988), in order to fully resolve the wall-bounded turbulence using spectral methods, the streamwise and spanwise grid resolutions should satisfy x + < 15 and z + < 8.One would expect the corresponding requirement for a finite-volume solver to be more stringent.In the current study, our analyses mainly focus on the cases with Re τ = 180.In lieu of any established criteria for DNS grid resolutions, we wanted to get some guidance from the literature.So, we summarized different grid resolutions used in the past DNS studies of neutral and stratified channel flows in Table 2.Then, to be on the conservative side, we chose a relatively fine grid resolution reported in Table 2 for the present study.Several verification cases (cases V1-V5 in Table 1) are conducted which show that the simulation results from the current solver agree well with the benchmark results for both neutral and stratified cases (see Appendix A for details).In addition, we increase the horizontal grid resolution by 50 % and no essential change is observed in the simulation results (see Appendix B).These results indicate that the grid resolution, employed in this study for Re τ = 180, is sufficient for the DNS of stably stratified flows.In order to evaluate the impact of Reynolds number on our key findings, we also conduct two additional simulations with Re τ = 395 and 550.The grid resolutions for Re τ = 395 are similar to those used for Re τ = 180 (see Table 1); however, for Re τ = 550, a relatively coarse resolution is used.For the case with Re τ = 550, an estimation of the ratio of Nu is Nusselt number; L x and L z are the domain lengths in the streamwise (x) and spanwise (z) directions, respectively; h is the openchannel height in the vertical (y) direction (h is set to be 1.0 for the current study); N x , N y , N z are the cell numbers in x, y, z directions, respectively; and x + , y + , z + are the normalized grid resolutions (e.g., x + = xu τ /ν).The smaller and larger y + represent the vertical grid resolution at the bottom and the top walls, respectively.

Case
Re denote the open-channel and closed-channel flows, respectively; and FD, FV, and SP represent finite-difference, finite-volume, and spectral methods, respectively.The numbers in the FD and FV methods represent the order of accuracy for the adopted spatial discretization schemes.grid size to the Kolmogorov length scale following Nieuwstadt (2005) shows that x ≈ 7η and z ≈ 3.5η, where η is the Kolmogorov length scale.Although these grid sizes are near the lower limit of the DNS resolution (Moin and Mahesh, 1998), it is nevertheless acceptable for a qualitative estimation of the impact of Reynolds numbers in the current study.

Results
According to the linear stability theory in Gage and Reid (1968), the Poiseuille channel flow is stable for approximately Ri τ > 900 at Re τ = 180.As was mentioned in Sect.order to obtain comprehensive characteristics of intermittent turbulence at Re τ = 180.Specifically, time series and time height plots are presented to delineate the temporal variations of turbulence with various stratification levels.In parallel, different horizontal (x-z) and spanwise (x-y) planes are shown to describe various patterns of spatial intermittency.In addition to these visual documentations, several statistics (e.g., turbulence enhancement index, turbulent fraction) are provided for quantification of turbulence intermittency.Furthermore, a strong emphasis is placed on understanding the relationship between gradient Richardson number and turbulence intermittency.The impact of Reynolds number on the aforementioned characteristics of intermittent turbulence is discussed at the end of this section.

Stratification effects on mean flow and turbulence statistics profiles
Before investigating the spatiotemporal features of intermittency, it is informative to document the influence of stable stratification on the overall flow structure in the open channel by analyzing the mean flow and turbulence statistics profiles.Vertical distributions of mean streamwise velocity, temperature, and turbulent momentum and heat fluxes with different Ri τ are shown in Fig. 3 (Re τ = 180).One can see that the flow accelerates with increasing Ri τ in the channel.The slope of the u profiles remains invariant between different cases in the near-wall region (y/ h < 0.1), implying that a constant viscous wall shear is enforced by the streamwise pressure gradient (see P * in Eq. 2).In the central-channel region, however, the u profiles change from well-mixed (Ri τ = 0) to curved (Ri τ > 0); these profiles qualitatively agree with those observed in the closed-channel simulations (Armenio and Sarkar, 2002;García-Villalba and del Álamo, 2011).In addition, for the case with Ri τ = 480, a parabolic profile of streamwise velocity is observed, indicating that the flow is laminarized; our simulation agrees with the linear stability theory in Gage and Reid (1968).The acceleration of mean velocity is caused by the reduced vertical mixing due to stable stratification, which can be clearly seen in Fig. 3b.The maximal and overall values of turbulent momentum flux (−u v ) decrease with increasing Ri τ , and the peak −u v moves from y/ h ≈ 0.18 (Ri τ = 0) to y/ h ≈ 0.25 (Ri τ = 240).For the case with Ri τ = 480, the flow becomes laminar and the turbulent momentum flux decreases to zero at all channel heights as would be expected.
The temperature profiles are also modified by the stratification.With increasing Ri τ , the gradient of mean temperature decreases near the wall, implying that the surface heat flux is significantly reduced by the stratification.This can be clearly indicated by the Nusselt number in Table 1.The Nusselt number, defined as the heat transfer due to turbulence with respect to its laminar value (pure conduction), can be used to describe how close the flow is to a laminar state: One can see that the Nusselt number equals the nondimensional temperature gradient at the bottom wall.For the neutral stratification case (Ri τ = 0), the Nusselt number is relatively high (Nu = 4.64), whereas for the very stable case (Ri τ = 240), Nu = 1.72, implying that the flow is becoming closer to laminar.Eventually, for the case with Ri τ = 480, the flow is fully laminarized and Nu = 1.0.This is further confirmed in Fig. 3c by the fact that the mean temperature profile for Ri τ = 480 is linear.
Under the influence of the stable stratification, the turbulent heat flux is significantly suppressed, which can be seen in Fig. 3d.As was noted by Armenio and Sarkar (2002), the gradient of mean temperature is related to the turbulent heat flux.This can be illustrated by integrating the Reynoldsaveraged temperature equation in channel flows: In the near-wall (y/ h < 0.1) and top-wall (y/ h > 0.9) regions, the turbulent heat fluxes (−T v ) change sharply in the vertical (y) direction for Ri τ = 0 − 240.By virtue of Eq. ( 5), sharp changes in temperature profiles near the walls are expected as shown in Fig. 3c.Please note that with increasing stratification, both the gradients in turbulent heat flux and mean temperature decrease rapidly in a monotonic fashion for any stratification level.In the central-channel region, the −T v remains almost unchanged resulting in a nearly constant gradient of mean temperature.Similar to that observed for the turbulent momentum flux, the turbulent heat flux decreases to zero for Ri τ = 480.In summary, the stable stratification significantly reduces the transport and mixing of turbulent flows in the vertical direction as would be expected.Therefore, in the following subsections, the vertical velocity field will be primarily utilized to probe if the simulated flow becomes laminarized or intermittent with high stratification levels.

Temporal characteristics of intermittent turbulence
To analyze the temporal characteristics of intermittent turbulence, several probe points are placed in the channel to monitor the time variations of the stratified turbulent flows.Time series of vertical velocity in the channel with different Ri τ are shown in Fig. 4. For the weakly stratified case (Ri τ = 30), the turbulent flow is continuous and no intermittency can be observed.However, with increasing the stratification level (i.e., Ri τ = 120), intermittent turbulence can be clearly seen, especially at the near-wall location (around y + = 10).Interestingly, the bursting and decaying turbulence follows quasiperiodic cycles; approximately 6-7 bursting events are observed in every 10 non-dimensional time (t u τ / h ≈ 10) at y + = 10.In contrast to the near-wall region, in the upperchannel region (i.e., around y + = 150) the intermittency can be hardly observed.We would like point out that this heightdependency behavior was earlier observed in a wind-tunnel experiment by Ohya et al. (2008).By further increasing the friction Richardson number to Ri τ = 240, intermittency becomes more pronounced.For example, a large segment of decaying turbulence is found in approximately 10 < t u τ / h < 15.In contrast to the Ri τ = 120 case, the time interval between the turbulent burst and decay is much less rhythmic and instead more irregular.Eventually, when the stratification level increases to Ri τ = 480, the turbulence is fully laminarized and no turbulent fluctuation is found.
Time series of streamwise velocity and temperature in the channel with Ri τ = 120 are shown in Fig. 5. Similar to what was found for the vertical velocity, turbulent bursts of streamwise velocity and temperature can be easily identified in the near-wall region.However, compared to the "symmetric" distributions of bursting vertical velocity (i.e., the vertical ve- locity is equally fluctuating up and down around zero mean), ramp-like variations are found for the streamwise velocity and temperature time series; the u and T increase sharply just after the bursting events and then gradually decrease to their laminar value.This phenomenon is especially evident in the near-wall region.It is quite exciting to see that the ramplike time series of streamwise velocity and temperature, simulated in the current DNS study, remarkably resemble those measured in stratified wind-tunnel flows in Ohya et al. ( 2008) (please refer to Figs. 12 and 13 in their work).
As were depicted in Figs. 4 and 5, the intermittent characteristics vary at different vertical locations.In order to comprehensively illustrate the temporal variability of intermittent turbulence across the entire channel, time height plots of square of instantaneous vertical velocity, namely v 2 = (v − v) 2 , are shown in Fig. 6.Here, v denotes the long-term temporal average of v (essentially, v ≈ 0).v 2 can be used to measure the bursting intensity, and high (low) v 2 represents turbulent (laminar) flow states, respectively.In the vertical direction, v 2 is generally higher in the central-channel regions than near the walls.At the bottom and the top walls, v 2 decreases to zero as a result of the constraint of the boundary conditions (v = 0).Alternating distributions of high and low v 2 can be clearly observed at most of the channel heights (except in close proximity of the walls) which indicates the burst and decay of turbulent flows.Moreover, an interesting temporal feature can be observed that the turbulent bursts are not synchronous along the vertical direction but instead are detected earlier at higher vertical locations.For example, for Ri τ = 120 at t * ≈ 39, the flow becomes turbulent at the channel center (y/ h ≈ 0.5), whereas at the same time, the flow remains laminar in the near-wall region (y/ h ≈ 0.1).There is a time lag ( tu τ / h ≈ 0.5) for the busting events between the central-channel and the near-wall regions.This time lag is related to the angular orientation of the near-wall turbulent flow structures which are elaborated on in Sect.3.3.
In order to quantitatively analyze the vertical distribution of turbulent intermittency, the turbulence enhancement index (TEI) is used to detect the turbulent bursting events.The TEI is defined as the ratio of increase of vertical velocity variance (v 2 ) between two consecutive windows with respect to the variance in the first window (Eq.6): where the subscripts 1 and 2 represent the current and next windows, respectively (Nakamura and Mahrt, 2005).The window is the prescribed time span for the calculations of v 2 .
The choice of window size should be small enough to eliminate the influence of wave-like motions of vertical velocity and large enough to capture the turbulent bursts (Nakamura and Mahrt, 2005).A window range of t + w = tu 2 τ /ν = 90 is used here and the interval between two consecutive windows is set to be t + i = 0.4t + w .The threshold of TEI is set to be 3, above which the rapid increase of velocity is considered to be a turbulent burst.In order to avoid an insufficiently small (spurious) interval between two consecutive bursting events, the minimal bursting interval is set to be t + min = 0.8t + w , below which multiple bursting events are considered to be a single event.In addition, we would like to emphasize that the turbulent bursting events obtained from the TEI approach can be very sensitive to the calculation settings mentioned above, e.g., the threshold.By utilizing a relatively high TEI threshold, we exclude classical bursting events in neutral channel flows, e.g., the burst-sweep cycles (Johnson, 1998) generated by the coherent vortices in the near-wall region.In other words, we only consider the bursting events under the impact of stable stratification.A detailed discussion on the determination of the TEI threshold can be seen in Appendix C.
Vertical distribution of turbulent bursting events in the channel is shown in Fig. 7.One can see that, no bursting event can be detected for Ri τ = 30 or 480, which coincides with the previous analyses.For Ri τ = 120 and 240, the most frequent bursting events can be found in the near-wall region, whereas in the central-channel region, the bursting frequency decreases drastically.The peak bursting frequency for Ri τ = 120 (32 events) is larger than that for Ri τ = 240 (26 events); however, for the case with Ri τ = 120, the bursting frequency decreases more rapidly especially in y/ h > 0.3.Interestingly, it is found that the vertical distribution of turbulent bursting events has a similar trend compared to that of turbulent momentum flux (see Fig. 3b).We speculate that this trend is attributed to the competition between the shear generation and negative buoyancy (due to stable stratifications).Due to the high turbulent momentum flux and therefore the stronger turbulent mixing in the near-wall region, the turbulent bursts are much more frequent near the wall than in the upper-channel region where the dominant buoyant force significantly suppresses the turbulent motion.
Using observational data based on the CASES-99 experiment, Nakamura and Mahrt (2005) analyzed naturally occurring intermittent turbulence in the atmosphere.The authors calculated TEI using the velocity time series collected from the eight levels of sonic anemometers on a 60 m tall meteorological tower.It was found that the frequency of turbulent bursting events decreased with height; the bursting events happened 5-6 times a night at and under 20 m, and 3-4 times a night above 20 m.It is quite encouraging to find that the DNS-based trend of turbulent bursting events in y/ h > 0.1 qualitatively agrees with that observed in atmospheric boundary layer flows by Nakamura and Mahrt (2005).

Spatial characteristics of intermittent turbulence
As was mentioned in Sect. 1, intermittent turbulence varies in both space and time.Although intermittent time series signals have been widely observed, little is known about the spatial distribution of these intermittent structures in the atmosphere.In recent years, DNS studies of stably stratified channel flows, conducted by Flores and Riley (2011), García-Villalba and del Álamo (2011), Brethouwer et al. (2012), and Ansorge and Mellado (2014), have presented intriguing patterns of spatial intermittency.In light of these findings, in this subsection, we further investigate the spatial characteristics of intermittent turbulence in detail.Instantaneous contours of v 2 at a horizontal (x-z) plane are shown in Fig. 8 for four simulations of varying Ri τ .The x-z plane is located at y + = 45, and the definition of the x-z plane can be seen in Fig. 2. For the case with Ri τ = 30, the entire flow is turbulent at the x-z plane.With increasing the friction Richardson number to Ri τ = 120, oblique "bands", with an angle ∼ 30 • with respect to the streamwise direction, can be seen at the x-z plane.These bands alternate between high and low v 2 (i.e., corresponding to turbulent and laminar flows, respectively), and the space between two turbulent bands is ∼ 12h.These oblique patterns qualitatively coincide with previous studies of stratified channel flow in García-Villalba and del Álamo (2011), and Brethouwer et al. (2012).For the case with Ri τ = 240, a coexistence of laminar and turbulent flows is also observed; however, in this case the patterns are not organized, rather randomly distributed.Eventually, with increasing the friction Richardson number to Ri τ = 480, the entire flow at the x-z plane becomes laminar, which concurs with the time series results in Fig. 4d.
The aforementioned results show that, for Ri τ = 120 and 240, the coexistence of laminar and turbulent flows can be observed at a given time.As was stated in Sect. 1, previous studies have shown complex transient phenomena (Flores and Riley, 2011) in stably stratified channel flows.In order to verify if laminar and turbulent flows can coexist continually, instantaneous contours of v 2 at different times (t * = 36, 38, 40) are shown in Fig. 9.One can see that, with increasing time, continuous coexistence of laminar-turbulent flows can be observed, which coincides with those reported in García-Villalba and del Álamo (2011) and Brethouwer et al. (2012).Moreover, a persistent pattern of oblique bands of laminarturbulent flow can be found for the case with Ri τ = 120; the oblique bands seem to propagate downstream with similar  widths and spacings.Once again, no particular temporal pattern can be found for the Ri τ = 240 case.
In the vertical direction, the spatially intermittent patterns can vary under different flow conditions.For example, in Couette flow, the spatial patterns can extend from the bottom to the top walls at low Reynolds numbers, whereas at high Reynolds numbers, these patterns can be confined to each wall, and uncorrelated (Brethouwer et al., 2012).Addition- ally, in stably stratified closed-channel flows, laminarization has been observed near the top wall, whereas near the bottom wall, the flow remained turbulent (Iida et al., 2002).In order to study the vertical distribution of the spatial patterns, instantaneous contours of v 2 at different horizontal planes are shown in Fig. 10.One can see that the oblique bands can be found at all of the channel heights; however, the intensity of these turbulent bands varies.In the near-wall region, the turbulent band is very weak, whereas at y + = 45, the intensity of v 2 in the turbulent band becomes much higher and the interface between laminar and turbulent flow can be easily identified.However, in the central-channel region, the turbulent band expands and the laminar-turbulent interface becomes diffused.At a higher channel location, the turbulent intensity in the oblique bands decreases again.In conjunction with the results in Fig. 9, it is indicated that an oblique pattern of laminar-turbulent coexistence is found at Re τ = 180 and Ri τ = 120 in the open-channel flow, and this pattern can steadily persist through time at all of the channel heights.
Additionally, by increasing the x-z plane from y + = 10 to 135, one can observe that the center of the turbulent band moves downstream ∼ 7h (see the red-dashed line for the approximate center of the turbulent bands in Fig. 10).This indicates that the turbulent bands are not upright in the vertical direction.This can be clearly seen in Fig. 11, where the contours of v 2 at the mid-spanwise plane are shown (see Fig. 2 for the definition of the spanwise plane).In Fig. 11a, the black-dashed lines denote the approximate boundary of the turbulent bands at different vertical locations and the reddashed lines represent the width of the turbulent bands observed in the x-z plane.One can see that, in the vertical direction, the turbulent bands lean downstream, and the inclination angle with respect to the streamwise (x) direction is approximately 5 • (Ri τ = 120).For the case with Ri τ = 240, similar inclination can also be found.The inclined structure conforms well with our previous observations in the context of the time height plots of v 2 in Fig. 6.When these forward leaning turbulent bands move downstream, the turbulent bursts are detected earlier in the higher vertical locations.
The aforementioned results clearly show that a coexistence of laminar and turbulent flows can be found in the stratified open-channel flows, and the spatial characteristic may vary with respect to vertical location.In order to quantify the turbulence coverage with respect to different x-z planes, turbulent fraction (F T ) is used to analyze the spatial distribution of intermittent turbulence.F T is defined as the spatial area, over which the flow is considered to be turbulent, normalized by the area of the entire x-z plane (Duguet et al., 2010).The calculation of F T consists of two steps.First, the square of instantaneous vertical velocity (v 2 ) is filtered using a Gaussian filter at a certain x-z plane.Then, based on the filtered v 2 (i.e., v 2 f ), a threshold is set to separate the turbulent and laminar flows for the F T calculation.For the Gaussian filter function G(r) = (2π σ 2 ) −1 e −0.5r 2 /σ 2 , the variance is chosen to be σ 2 = 2 /12 in order to match the second moments of the Gaussian and box filters (Leonard, 1975;Manneville, 2011;Pope, 2000).The filter width x,z is chosen to be 100 normalized wall distance ( x,z u τ /ν = 100) and the threshold is set to be a 10 % value of the averaged v 2 over the whole plane.
Vertical distribution of F T with different Ri τ is shown in Fig. 12.One can see that, for all the Ri τ , F T drop drastically to zero at both of the bottom and the top walls owing to the constraint of the boundary conditions (v = 0).For Ri τ = 30, as was mentioned above, there is no intermittency observed and F T = 1.0.On the other hand, for Ri τ = 480, the entire flow in the channel is fully laminarized so F T = 0 at all of the vertical locations.For Ri τ = 120, however, F T is observed to vary vertically.F T is relatively low (∼ 0.63) at y/ h = 0.06, whereas at y/ h = 0.75, F T increases to ∼ 0.9.For Ri τ = 240, a similar trend can be found except that the overall value of F T decreases for ∼ 0.15.In conjunction with the aforementioned results from Fig. 7, one can see that the higher frequency of turbulent bursts is generally associated with a lower turbulent fraction at each vertical location.
An intriguing phenomenon in stratified flows is the internal gravity wave.As reported in Komori et al. (1983), wave-like turbulence motions were measured in the core region of strongly stratified open-channel flows.These wave motions were then corroborated by the LES simulations in Kosović and Curry (2000), Saiki et al. (2000), and Armenio and Sarkar (2002).As shown in Armenio and Sarkar (2002), in the strongly stratified case, the prominent internal waves can even extend to the near-wall region.Based on the DNS results, here we further investigate the behavior of internal waves in intermittent turbulence.Instantaneous contours of temperature at the mid-spanwise plane is shown in Fig. 13.In the upper-channel region (y/ h ≈ 0.8), internal waves can be identified for both cases (Ri τ = 120 and 240).Compared with the prominent internal waves observed in Armenio and Sarkar ( 2002), here the wave motions seem to constraint themselves in y/ h > 0.7.In the meantime, intermittent turbulence can be observed in the lower-channel region (please also refer to Figs. 4 and 5).Our DNS results indicate that there is a coexistence of significant wave activity and intermittent turbulence in stably stratified flows, which coincides with the field observations in the atmospheric boundary layer (e.g., Sun et al., 2004).This coexistence implies that the internal gravity wave can be a possible trigger to generate intermittent turbulence; however, this topic is beyond the scope of the current study, and will be investigated in our future work.

Relationship between spatial and temporal intermittency
As reported in Coulter and Doran (2002), the temporal intermittency was found to be spatially correlated in the atmospheric boundary layer.The authors analyzed the time series of turbulent heat flux measured at different site locations of the CASES-99 field campaign.It was shown that the temporal distributions of intermittent turbulence were strongly correlated between different sites, separated by at least 1 km.For example, the episode of turbulent bursts observed at one site location was found to have similar features with that at another site with a time delay of half an hour (please refer to the discussion in Sect.4.2 in Coulter and Doran, 2002).
These results indicate that the turbulent bursting event may not be a fully localized phenomenon in stably stratified flows.
In other words, it is possible that the disturbance of intermittent turbulence can be advected by mean flows.In light of these findings, in this subsection we evaluate whether our DNS results can corroborate the strong correlation between spatial and temporal intermittency as observed in Coulter and Doran (2002).Profiles of mean velocity (u) and convection velocity (u c ) of intermittent turbulence are shown in Fig. 14.The convection velocity is calculated based on the space-time correlation approach reported in Kim and Hussain (1993).Specifically, the convection velocity is determined by u c = x max / t, where x max is the streamwise distance for which the space-time correlation of v 2 , i.e., R( x, y, t), is maximum at a given time delay ( t).In the current study, the time delay is chosen to be 0.5 non-dimensional time ( tu τ / h = 0.5).Following Brethouwer et al. (2012), in order to consider only the large-scale motion of laminarturbulent flows, a Gaussian filter is applied to the v 2 fields before calculating the space-time correlation.Both u and u c are normalized by the bulk/average velocity (u b ) in the channel.One can see that for both Ri τ = 120 and 240 cases, the convection velocity is lower than the bulk velocity.Moreover, the profile of u c seems to be essentially invariant along the vertical direction, which coincides well with that observed in Brethouwer et al. (2012).In Brethouwer et al. (2012), the averaged u c in the channel was found to be ∼ 0.8u b for a case with Re τ = 192.In the current study, a slightly lower averaged u c is observed, i.e., approximately 0.76u b and 0.73u b for Ri τ = 120 and Ri τ = 240, respectively.Now, using Figs.4, 9, and 14, we can analyze the relationship between spatial and temporal intermittency.As mentioned in Fig. 9, the laminar-turbulent bands seem to propagate downstream.When the turbulent and laminar bands pass through the probe points in the channel, the time series of velocity or temperature should alternate between turbulent and laminar states.In order to test this concept, we analyze the case with Ri τ = 120.In Fig. 9, one can see that the spacing between two turbulent bands is approximately 0.5L x (4π h).With a constant convection speed of 0.76u b (∼ 19u τ ), the interval between two bursting events ( t burst ) should be about 0.66 non-dimensional time.As discussed in Fig. 4b, approximately 6-7 bursting events can be observed in every 10 nondimensional time.This corresponds to a bursting interval of 0.6-0.7 non-dimensional time, which coincides well with the above analysis ( t burst ≈ 0.66).In other words, a strong correlation between the spatial and temporal intermittency is observed in the simulated stratified flows.This finding corroborates the analysis in Coulter and Doran (2002) that intermittent turbulence may not be a fully localized phenomenon, and it can be advected by mean flows.

Relationship between intermittency and stability parameter
In stratified shear flows, the gradient Richardson number (Ri g ) has been identified as the key parameter which determines the influence of stratifications.The gradient Richardson number is defined as the ratio between the buoyancy frequency (N) and the shear rate (S): Compared with the bulk quantities, such as Ri τ and h/L, the gradient Richardson number is a very useful parameter which reflects the local stability of the stratified flow.Linear stability analysis performed by Miles (1961) and Howard (1961) shows that Ri g ≥ 0.25 is a sufficient condition for laminarization of inviscid stratified shear flows.As mentioned in Sect. 1, the Blackadar's hypothesis (Blackadar, 1979) (see also , Stull, 1988;Mahrt, 1999) associated the intermittent turbulence with the modulation of Ri g .In this subsection, we would like to explore if Blackadar's hypothesis, i.e., a strong relationship between Ri g evolution and intermittent turbulence, holds for the DNS results.
As a first step, we analyze the mean profiles of Ri g in the channel, which are shown in Fig. 15.All the three Reynolds numbers (Re τ = 180, 395 and 550) are considered, and more analyses of Reynolds number sensitivity will be presented in Sect.3.6.One can see that, for all the stratification levels (Ri τ ), the local stability (Ri g ) monotonically increases along the vertical direction as a result of decreasing mean shear.Ri g becomes singular at the top wall where the mean shear vanishes.At a given vertical location, Ri g generally increases with increasing overall stratification level (Ri τ ).In addition, one notices that the slope of the Ri g profiles increases sharply at Ri g ≈ 0.2 for all the Re τ cases.The mean profiles of Ri g , as well as the sharp increase in Ri g slope, qualitatively coincides with those observed in the stratified closed-channel flows in Armenio and Sarkar (2002) and García-Villalba and del Álamo (2011).
Next, we discuss the overall impact of Ri g on turbulence statistics.Distribution of Ri g with respect to correlation coefficient of momentum and heat fluxes is shown in Fig. 16.One can see that in the low Ri g region, the correlation coefficients are relatively invariant.However, in 0.15 < Ri g < 0.2, there is a sharp drop of correlation coefficients of both u v and T v .The agreement of this sharp drop trend is relatively good between different Re τ cases, especially for u v .These results corroborate the linear stability analysis by Miles (1961) and Howard (1961), as well as the LES results reported in Armenio and Sarkar (2002), implying that the vertical motion of turbulence activity is significantly suppressed in the region with Ri g > 0.25.
In addition to the impact of Ri g on the overall stability of stratified flows, in this study, we also explore the relation between Ri g and local turbulence.This is done by correlat-  ing the instantaneous gradient Richardson number (henceforth Ri gi ) and v 2 .The definition of Ri gi is similar to that in Eq. ( 7), except that the instantaneous gradients of temperature and velocity are utilized: The concept of instantaneous gradient Richardson number is similar to the pointwise gradient Richardson number proposed in Mason (1994).Instantaneous contours of Ri gi and v 2 at the x-z plane with y + = 90 are shown in Fig. 17.One can observe that the Ri gi pattern (Fig. 17a) remarkably resembles the laminar-turbulent pattern in Fig. 17b.Similar to that observed in Sect.3.3, the Ri gi pattern is also found to convect downstream, closely connected to the laminarturbulent pattern (not shown).By using the approach introduced in Sect.3.3, we separate the turbulent and the laminar regions in the x-z plane.The plane-averaged Ri gi in the laminar region and the turbulent region are 0.34 and 0.19, respectively.The strong spatial-correlation between the high (low) Ri gi and the laminar (turbulent) regions, observed at the central-channel, well corroborates the Blackadar's hypothesis.In addition, the plane-averaged Ri gi in the laminar region coincides with the linear stability analysis that Ri g > 0.25 is sufficient for laminarization.
Similarly, we also plot the profiles of plane-averaged Ri gi for the turbulent (Ri T gi ) and the laminar (Ri L gi ) regions (see Fig. 18).As expected, Ri L gi is generally higher than Ri T gi .In the bottom and the top wall regions, the magnitudes of Ri T gi and Ri L gi are nearly identical.While in the central-channel region, Ri gi becomes quite distinguishable between the laminar and the turbulent regions, especially for the case with Ri τ = 120.The profiles of Ri T gi and Ri L gi qualitatively resemble the profiles of Ri g , calculated at the locations inside and outside of the oblique turbulent-bands, in Brethouwer et al. (2012).
In addition to investigate the spatial-correlation between Ri gi and intermittent turbulence, we also explore the temporal evolution of Ri gi associated with turbulent bursts.In the wind-tunnel experiments performed by Ohya et al. (2008), it was observed that the local Ri g was closely correlated to turbulent bursting events; for Ri g < 0.25, the velocity time series exhibited turbulent behavior, whereas for Ri g > 0.25, the velocity time series became quiescent.Here we would like to see whether our DNS results can corroborate the above dynamical characteristic as reported in Ohya et al. (2008).Similar to Ohya et al. (2008), a moving-average filter is applied on the simulated time series of Ri gi (hereinafter called Ri gi ) with a filter width of tu τ / h = 0.2.Time series of v 2 and Ri gi at the central-channel are shown in Fig. 19 (Re τ = 180, Ri τ = 120).One can see that, during the period with quiescent portions of v 2 (i.e., laminar state of flow; see Fig. 19b for reference), Ri gi is generally larger than the critical value (0.25).While in the turbulent portions of v 2 , Ri gi drops sharply (< 0.25).The temporal correlation between Ri gi and intermittent turbulence qualitatively coincides with those reported in Ohya et al. ( 2008) (see Fig. 16 in their work).Furthermore, it is suggested that Ri gi can be a good localized parameter to predict spatial and temporal evolution of intermittent turbulence in stably stratified flows.It is important to note that, in this study, the prominent dynamical correlation between Ri gi and turbulent burst is most distinguishable in the central-channel region (0.4 < y/ h < 0.6).However, in the near-wall region, the above dynamical relation becomes less evident (not shown).We speculate that this incoherent correlation is due to the significant impact of viscosity in the near-wall region, as discussed in Turner (1979) and Defina et al. (1999).The detailed impact of viscosity on the correla- tion between local stability and intermittent turbulence will be studied in the future.

Sensitivity of intermittency with respect to Reynolds numbers
The above analyses mainly focused on a relatively low Reynolds number case (Re τ = 180).In this subsection, we discuss whether the spatial and temporal characteristics of intermittent turbulence, observed in Sect.3.2 to 3.5, hold for higher Reynolds numbers.To do this, we perform two additional simulations with Re τ = 395 (Ri τ = 960) and Re τ = 550 (Ri τ = 2880).Time series of vertical and streamwise velocity for these higher Reynolds number cases are shown in Figs.20 and 21, respectively.One can see that the temporal evolutions of both vertical and streamwise velocities qualitatively resemble those shown in Figs. 4 and 5. Intermittent turbulence is much more distinguishable in the near-wall region (y + = 10) than in the central-channel (y + = 150).Moreover, the ramp-like variations of streamwise velocity can be clearly observed for both Re τ = 395 and 550, especially near the wall.Similar ramp-like variation is also found in the temperature time series (not shown).Next, we analyze the vertical distribution of turbulent bursting events for the higher Reynolds number cases (please refer to Fig. 22).Similar to what was observed for Re τ = 180 (see Fig. 7), the turbulent bursting events near the bottom wall are more frequent than that in the upper-channel region.Moreover, compared with Re τ = 180, the bursting events decrease much more sharply along the vertical direction.Lastly, we evaluate whether the close connection between Ri gi and intermittent turbulence can be observed for the higher Reynolds number cases.Time series of Ri gi and v 2 for Re τ = 395 are shown in Fig. 23.Again, the dynamical correlation between the temporal variation of Ri gi and intermittent turbulence is clearly observed, coinciding with what was observed in Fig. 19 for Re τ = 180.
A similar correlation can also be observed for the case with Re τ = 550 (not shown).Above all, the spatial and temporal characteristics of intermittent turbulence, observed for the higher Reynolds numbers, are qualitatively similar to those shown in Re τ = 180, although there are some differences in detail (e.g., the sharper drop of turbulent bursting events in vertical direction).

Conclusions
In this paper, direct numerical simulation is utilized to simulate open-channel flows under stably stratified conditions.The spatiotemporal characteristics of intermittent turbulence are discussed, and several conclusions are made as summarized below.Under stably stratified conditions, the vertical transport of turbulent momentum and heat fluxes is suppressed and the mean flow and turbulence statistics profiles are significantly modified.For the weakly stratified case with Ri τ = 30 (Re τ = 180), the turbulence is continuous and no intermittency is found.With increasing the stratification level to Ri τ = 120 and 240, intermittent turbulence is clearly ob- served from the time series of velocity and temperature, especially in the near-wall regions.Eventually, the entire flow is laminarized in the channel for Ri τ = 480.This result is in agreement with linear stability analysis by Gage and Reid (1968).Turbulence enhancement index is utilized to detect turbulent bursting events.It is found that the most frequent bursting event occurs in the near-wall region; in the upperchannel region where the buoyant force is dominant, the bursting event frequency decreases significantly.The vertical variation of bursting events qualitatively agrees with those observed in laboratory and field experiments (Ohya et al., 2008;Nakamura and Mahrt, 2005).
Stable coexistence of laminar and turbulent flows is found at various horizontal planes.For the case with Ri τ = 120 (Re τ = 180), alternating laminar and turbulent bands are observed at an angle ∼ 30 • to the streamwise direction.This oblique pattern persists through time and the bands propagate downstream with similar widths and spacings.With increasing Ri τ to 240, however, the coexistent pattern of laminar and turbulent flows becomes less organized.Turbulent fraction is utilized to analyze the vertical variation of intermittent turbulence.It is found that, for Ri τ = 120 and 240, the turbulent fraction is relatively low in the near-wall region, whereas it increases rapidly in the upper-channel region.
Strong correlation between the spatial and temporal intermittency is observed in the channel.The laminar-turbulent bands are found to propagate downstream with a uniform convection velocity, and the time series signal alternates between turbulent and laminar states when the corresponding laminar-turbulent bands pass though the probe points.This prominent correlation, observed in the direct numerical simulation (DNS) results, corroborates the observation-based analyses in Coulter and Doran (2002) that intermittent turbulence may not be a fully localized phenomenon, and it can be advected by mean flows in the atmosphere.
The instantaneous gradient Richardson number (Ri gi ), a local stability parameter reflecting the balance between wind shear and buoyancy, is found to be a quite relevant quantity to predict the spatial and temporal characteristics of www.nonlin-processes-geophys.net/22/447/2015/ Nonlin.Processes Geophys., 22, 447-471, 2015 intermittent turbulence.It is found that the alternating highlow Ri gi pattern remarkably correlates with the laminarturbulent bands at the horizontal plane in the central-channel region.Furthermore, the temporal evolution of the filtered Ri gi is strongly connected to the bursting and quiescent state of v 2 time series.Indirectly, our DNS results corroborate an age-old hypothesis by Blackadar (1979) that the intermittent turbulence is strongly correlated with the modulation of Ri g in the atmosphere.Additionally, the Reynolds number sensitivity study, performed at the end of the paper, suggests that the spatial and temporal characteristics of intermittent turbulence, observed for the cases with Re τ = 180, qualitatively holds for higher Reynolds numbers (Re τ = 395 and 550).
It is evident that some of the findings of the present work have striking similarities with results gleaned from atmospheric boundary layer (ABL) flows.However, given the Reynolds number disparity, these similarities should be considered with caution.In the near future, we will explore the strengths and weaknesses of a suite of large-eddy simulation subgrid-scale (LES-SGS) models in capturing low Reynolds number intermittency by using our DNS-based results as benchmark.Several contemporary LES-SGS models assume local balance of turbulent kinetic energy (TKE) production with its dissipation.We speculate that these models might not fare well in this exercise.On the other hand, we expect the SGS models with dynamic TKE formulations and energy backscatter options (e.g., Ghosal et al., 1995;Davidson, 1997;Kim and Menon, 1999) to be highly competitive.We envisage that a subset of these SGS models might be (not guaranteed) successful in mimicking the intermittent ABL flows.
Our future work will also involve characterization of intermittent turbulence from a complex systems perspective.Over the years, a number of simple dynamical systems have been reported to exhibit intermittency (Heagy et al., 1994;Marthaler et al., 2001;Toniolo et al., 2002).Several mechanisms were proposed, such as on-off intermittency, Pomeau-Manneville mechanism, etc. (Pomeau and Manneville, 1980;Platt et al., 1993), to explain such intrinsic dynamical behaviors.It was noticed that bursting activity can only be triggered within certain regions of a parameter space (Toniolo et al., 2002).More interestingly, the duration of quiescent period was found to follow certain power laws (Toniolo et al., 2002).We are interested to find out if analogous results can be discovered for intermittency in stably stratified flows.

Appendix C: Determination of TEI threshold
As mentioned in Sect.3.2, we chose a relatively high TEI threshold to exclude classical bursting events in neutral or near-neutral flows, e.g., burst-sweep cycles (Johnson, 1998).In this Appendix, determination of the TEI threshold is discussed.Time series of v 2 and TEI at y + = 10 is shown in Fig. C1.For Ri τ = 30, sharp spikes are observed in the v 2 time series (see Fig. C1a), which can be caused by the burst-sweep cycles in the near-wall region.However, the ratio between the magnitudes of bursting and quiescent v 2 is much less significant compared with that observed in the very stable case (see Fig. C1b).As a result, the TEI value for Ri τ = 30 is much less than that for Ri τ = 120 (see Fig. C1b  and d).In this study, we only consider the bursting events under the impacts of stable stratification, the TEI value is set to be 3.0 (see the red dash line in Fig. C1b and d) to exclude the classical bursting events in neutral or near-neutral flows.

Figure 2 .
Figure 2. Illustration of the computational domain, where x, y, and z are streamwise, vertical, and spanwise directions, respectively.The red short-dashed and blue long-dashed lines define a horizontal (x-z) plane and a spanwise (x-y) plane, respectively.

Figure 7 .
Figure7.Vertical distribution of turbulent bursting events with different Ri τ in the channel (Re τ = 180).The turbulent bursting events are calculated within 50 non-dimensional time.

Figure 11 .
Figure 11.Instantaneous contours of v 2 at the mid-spanwise plane (not to scale, Re τ = 180).Please refer to Fig. 2 for the definition of the spanwise plane.The black-dashed lines denote the approximate boundary of the turbulent bands at different vertical locations and the red-dashed lines represent the width of the turbulent bands observed in the x-z plane in Fig. 10.(a) Ri τ = 120, z/ h = 2π, (b) Ri τ = 240, z/ h = 2π.

Figure 15 .
Figure 15.Vertical distribution of gradient Richardson number (Ri g ) in the channel.

Figure 16 .
Figure 16.Distribution of Ri g with respect to correlation coefficient of (a) momentum and (b) heat fluxes.

Figure A2 .
Figure A2.Comparisons of variables for neutral flows, Ri τ = 0.The lines are the results from the current solver (OF), and the symbols are results from Moser et al. (1999) (MKM99).The production of TKE are normalized by ν/u 4 τ .(a) Mean velocity, (b) rms velocity, (c) turbulent momentum flux, and (d) production of turbulence kinetic energy.
− T top )/(T top − T bot ) and ρ * = (ρ − ρ top )/(ρ bot − ρ top ), respectively, where the top and bot subscripts represent the top and bottom walls, respectively.P * is the streamwise pressure gradient driving the flow and p * is the dynamic pressure.

Table 1 .
Summary of the simulation configurations used in this study.Re τ is friction Reynolds number; Ri τ is friction Richardson number;

Table 2 .
Summary of different grid resolutions used for DNS studies of neutral and stratified channel flows.Re τ is friction Reynolds number; Re g is Reynolds number based on geostrophic wind; Ri τ is friction Richardson number; Ri b is bulk Richardson number; h/L is stability parameter; x, y, z are streamwise, vertical, and spanwise directions, respectively; x + , y + , z + are the normalized grid resolution (e.g., x + = xu τ /ν); O and C