Coherent structures in stably stratified wall-bounded turbulent flows

Abstract

To date, a growing body of literature has documented the existence and impacts of coherent structures known as large- and very-large-scale motions within wall-bounded turbulent flows under neutral and unstable thermal stratification. These coherent structures can account for a considerable fraction of the overall turbulent transport and have been found to modulate small-scale turbulent fluctuations near the wall. In the context of stably stratified flows, however, the examination of such coherent structures has garnered relatively little attention. Stable stratification limits vertical transport and turbulent mixing within flows, which makes it unclear the extent to which previous findings on coherent structures under unstable and neutral stratification are applicable to stably stratified flows. In this study, we investigate the existence and characteristics of coherent structures under stable stratification with a wide range of statistical and spectral analyses. Outer peaks in premultiplied spectrograms under weak stability indicate the presence of large-scale motions, but these peaks become weaker and eventually vanish with increasing stability. Quadrant analysis of turbulent transport efficiencies (the ratio of net fluxes to their respective downgradient components) demonstrates dependencies on both stability and height above ground, which is evidence of morphological differences in the coherent structures under increasing stability. Amplitude modulation by large-scale streamwise velocity was found to decrease with increasing gradient Richardson number, whereas modulation by large-scale vertical velocity was approximately zero across all stability ranges. For sufficiently stable stratification, large eddies are suppressed enough to limit any inner–outer scale interactions.

1. Introduction

In recent decades, researchers have come to recognise the importance of coherent, organised structures within neutrally stratified wall-bounded turbulent flows, which can account for a large fraction of turbulent fluxes close to the wall (Corino & Brodkey 1969; Wallace, Eckelmann & Brodkey 1972; Willmarth & Lu 1972; Guala, Hommema & Adrian 2006; Balakumar & Adrian 2007; Marusic et al. 2010b; Wallace 2016). The streamwise spatial extent of these structures can exceed the depth of the flow itself, $h$. Within the inner layer of inertia-dominated wall-bounded flows ($Re_\tau = u_\tau h/\nu \gtrsim O(10^3)$, where $u_\tau$ is the friction velocity and $\nu$ is the kinematic viscosity, Hutchins & Marusic 2007a), turbulence can be generated by hairpin vortices that eject low-momentum fluid upwards (i.e. $u' < 0$ and $w' > 0$). This process serves to generate additional hairpin vortices that continually bound regions of low-momentum fluid (Head & Bandyopadhyay 1981; Meinhart & Adrian 1995; Adrian 2007). Conversely, structures within the outer layer can affect turbulence within the inner layer by sweeping high-momentum fluid downwards towards the wall (i.e. $u' > 0$ and $w' < 0$). These collections of hairpin vortices, known as large-scale motions (LSMs), have been widely studied in the fluid mechanics community since the early 1970s (e.g. Kovasznay, Kibens & Blackwelder 1970; Brown & Thomas 1977; Nakagawa & Nezu 1981; Murlis, Tsai & Bradshaw 1982; Wark & Nagib 1991; Adrian, Meinhart & Tomkins 2000; Ganapathisubramani, Longmire & Marusic 2003; Tomkins & Adrian 2003; Del Álamo et al. 2004). Large-scale motions are characterised by regions of high- and low-momentum fluid in the log-layer region of high $Re_\tau$ flows (also known as the inertial sublayer, which can overlap between the inner and outer layers) that scale as $O(h)$ in the streamwise direction and are comprised of several successive hairpin vortices propagating at similar speeds (Adrian 2007).

It has also been found that LSMs can organise into superstructures known as very-large-scale motions (VLSMs) that can extend $O(10h)$ in the streamwise direction (Cantwell 1981; Kim & Adrian 1999; Guala et al. 2006; Balakumar & Adrian 2007; Hutchins & Marusic 2007a; Marusic & Hutchins 2008). Studies examining VLSMs have only been possible more recently due to limitations of $Re$ accessible by laboratory set-ups and direct numerical simulations (DNS), but since the early 2000s they have been identified in turbulent channel flows (Del Álamo et al. 2004; Chung & McKeon 2010), pipe flows (Guala et al. 2006) and atmospheric boundary layers (ABLs; Tomkins & Adrian 2003; Hutchins & Marusic 2007a,b; Lee & Sung 2011).

Although LSMs and VLSMs primarily exist in the logarithmic region of wall-bounded turbulence, they have been found to influence turbulent fluctuations in the near-wall region (Hutchins & Marusic 2007b; Mathis, Hutchins & Marusic 2009a; Marusic, Mathis & Hutchins 2010a). More specifically, Mathis et al. (2009a) explored these inner–outer interactions through the lens of amplitude modulation (AM). By leveraging the strong correlations that exist between the large-scale signal and the small-scale fluctuations, Marusic et al. (2010a) developed a predictive model of near-wall turbulence that requires only the large-scale velocity in the logarithmic region as an input. They found strong agreement between the predicted and observed streamwise velocity energy spectra close to the wall across $2800 \le Re_\tau \le 19\ 000$ based on wind tunnel data, thereby demonstrating the utility in the decoupling procedure to examine the effects of LSMs and VLSMs in wall-bounded turbulent flows.

While this study investigates the effects of stable stratification on the existence and characteristics of LSMs, relatively little literature exists compared with analogous studies under neutral and unstable stratification. Thus, in § 1.1 we review the impacts instability has on LSMs and VLSMs before discussing some more recent studies focused on stably stratified flows in § 1.2.

1.1. Unstable stratification

It is generally well understood that buoyancy can significantly affect the nature of turbulence under unstable stratification, for example, in terms of integral length scales (Khanna & Brasseur 1998; Sullivan et al. 2003; Salesky, Katul & Chamecki 2013), the turbulence kinetic energy (TKE) budget (Wyngaard & Coté 1971; Salesky, Chamecki & Bou-Zeid 2017), velocity and temperature spectra (Kaimal & Finnigan 1994; Khanna & Brasseur 1998), the statistics of uniform momentum and temperature zones (UMZs and UTZs, respectively, Salesky 2023), and the morphology of organised motions (Khanna & Brasseur 1998; Weckwerth, Horst & Wilson 1999; Smedman et al. 2007; Lotfy et al. 2019; Jayaraman & Brasseur 2021). Salesky et al. (2017) explored the role of instability on the organization of motions within the convective boundary layer (CBL) using a suite of large-eddy simulations (LES; Stoll et al. 2020) at varying levels of instability, and demonstrated a transition between modes from quasi-two-dimensional horizontal convective rolls (HCRs) under weak surface heat fluxes relative to large mean wind shear towards open cellular convection reminsicent of Rayleigh–Bénard convection as instability increases. These HCRs are typically aligned within 10–20$^\circ$ of the geostrophic wind vector (Weckwerth, Wilson & Wakimoto 1996; Weckwerth et al. 1997, 1999).

One common area of focus when studying the morphology of LSMs and VLSMs involves quantifying their inclination angle relative to the surface, $\gamma$, which for neutral stratification is typically found to be $\gamma = 15^\circ$ (Brown & Thomas 1977; Rajagopalan & Antonia 1979; Marusic & Heuer 2007). This angle is commonly defined as $\gamma = \arctan (\Delta z/\Delta x^*)$ by determining the streamwise lag $\varDelta _x^*$ of the maximum value of, e.g. two-point correlation for streamwise velocity at height $z$ (e.g. Kovasznay et al. 1970; Brown & Thomas 1977; Rajagopalan & Antonia 1979; Boppe, Neu & Shuai 1999; Ganapathisubramani et al. 2005; Marusic & Heuer 2007; Hutchins et al. 2012). Under increasing instability, the inclination angles of LSMs in the CBL have been found to increase past $50^\circ$ (Salesky & Anderson 2018, and references therein), which is consistent with the topological transition towards vertical buoyant plumes at high instability. Using sonic anemometer data from the ABL, Li et al. (2022) also examined the relationship between stability, inclination angle and aspect ratio of coherent structures in the context of self-similar wall-attached eddies after Townsend (1976) (also see Woodcock & Marusic 2015; Marusic & Monty 2019). They found that coherent structures have an aspect ratio close to unity under neutral stratification, and become progressively taller and wider under increasing unstable stratification. For unstable conditions, they also found structures to be inclined at greater angles at larger scales as compared with smaller scales.

The changes in LSM and VLSM structure under convective conditions also can be detected by examining how turbulent transport efficiencies (fraction of the net flux in the downgradient direction) change for momentum versus scalars such as heat and moisture. Using atmospheric surface data, Li & Bou-Zeid (2011) found that under near-neutral stratification, momentum and scalars are transported by the same updrafts with high correlations. With increasingly unstable conditions, they observed a reduction in the transport efficiency of momentum paired with an increase in scalar transport efficiency, indicating these processes are governed by differing mechanisms related to the structure of vertical plumes. Through quadrant analysis (also referred to as conditional sampling; e.g. Wallace et al. 1972; Willmarth & Lu 1972; Holland 1973; Grossman 1984; Finnigan 2000; Wallace 2016), they further identified that under higher instability, vertical motions preferentially organise into rapid, intense updrafts compensated by longer, weaker downdrafts. Salesky et al. (2017) later confirmed these findings, further noting that these differences are related to the spatial distribution of individual quadrant events that are in turn affected by global stability.

Following the procedure of Mathis et al. (2009a), Salesky & Anderson (2018) examined the influence of instability on AM coefficients in simulated CBLs. By repeating the signal decoupling procedure with virtual tower data at multiple heights within CBLs across varying stabilities, they found the strongest correlations for the least convective cases considered. They also noted that significant correlations existed for all cases as long as there existed sufficient separation between inner and outer peaks in the premultiplied spectrograms. Their results indicated that small-scale fluctuating velocity, temperature and instantaneous second-order moments can be modulated by the large-scale streamwise and vertical velocity components associated with LSMs. Salesky & Anderson (2018) conclude with a conceptual model illustrating the effects of buoyancy on LSM inclination angles and how LSMs at varying stabilities act to modulate surface-layer turbulence.

1.2. Stable stratification

While a majority of research on LSMs focus on neutrally and unstably stratified flows, analogous investigations of stably stratified flows are not as prominent. Turbulence within stably stratified flows is difficult to observe or simulate due to the buoyant suppression of vertical motions that results in turbulence that is increasingly weak and localised in space and time with increasing stratification (Mahrt 1999; Lan et al. 2018). With increasing stability, turbulent eddies become decoupled from the surface and scale with local stability (Nieuwstadt 1984; van de Wiel et al. 2008), and eventually $z$ loses relevance as a characteristic length scale (the so-called $z$-less stratification regime, e.g. Wyngaard & Coté 1972; Dias, Brutsaert & Wesely 1995; Grachev et al. 2013).

Due to the nature of scales involved and the way large eddies interact with small-scale flow features, most of the studies examining coherent structures in stably stratified flows leverage DNS of channel or free-shear flows (e.g. García-Villalba & del Álamo 2011; Watanabe et al. 2018, 2019; Atoufi, Scott & Waite 2021; Gibbs, Stoll & Salesky 2023). Watanabe et al. (2019) confirmed the existence of hairpin vortices within stably stratified free-shear layers, noting their strong similarity to those typically observed in wall-bounded turbulent flows. They observed that while these hairpin vortices could be found throughout the shear layer, so-called superstructures (collections of multiple individual hairpin vortices that can be up to 10 times larger than the depth of the shear layer) only exist in the centre of the layer. The authors also observed cospectral peaks at the wavelengths associated with the horizontal extent of individual hairpin vortices, and they determined the composite superstructures are responsible for large peaks in density and velocity spectra at wavelengths associated with the streamwise length of these structures. In a DNS investigation of stratified channel flow, García-Villalba & del Álamo (2011) considered a wide range of stability and noted several key findings. Two-dimensional spectral energy density analysis indicated that the primary effect of stratification is to damp the large-scale modulation of intensity of near-surface streaks caused by global stability modes. Close to the surface, vertical motions are largely unaffected by stability as the flow is dominated by wall effects and coherent structures within the outer layer of the flow are not tall enough to penetrate down to the surface due to the suppression of vertical motions by negative buoyancy. They argue that stratification prevents the formation of larger-scale structures by damping turbulent vertical fluxes at those scales.

Observational studies in the stable atmospheric boundary layer (SBL) largely agree with these findings, particularly when vertical wind shear is weak, which enables the development of strong vertical temperature gradients due to the lack of vertical mixing (Lan et al. 2018). In these cases, turbulence becomes highly localised into thin layers that are completely decoupled from the surface. In weakly stable boundary layers with high levels of coupling, Lan et al. (2019) found that large eddies can contribute equally to both turbulent production and transport, resulting in fluxes that were nearly constant with height. However, for increasing stability, such large eddies do not contribute evenly thereby resulting in non-zero vertical gradients of fluxes. Lan et al. (2022) found that sudden events of wind profile distortion can trigger large eddies that penetrate downwards and initiate a transition towards decreased stability as they induce enhanced regions of turbulent transport, increased fluxes and reduced TKE and flux gradients across layers. With such weak turbulent motions, these studies elucidate the importance of large eddies in the SBL when they are able to penetrate across scales in the vertical.

Extending the analysis of stability effects on inclination angle to stably stratified channel flow, Gibbs et al. (2022) recently found that structures become increasingly inclined with height above the lower boundary up to $z/h=0.15$, where $h$ is the boundary-layer depth. Above this height, the inclination angles level off, which they discuss is indicative of a region where local $z$-less scaling behaviour may no longer exist (Grachev et al. 2013). Moreover, they found that the inclination angle decreases with increasing stability at all heights, and that angles inferred from buoyancy structures are larger than those from momentum.

Although these studies have provided foundational context on the existence of turbulent coherent structures in stably stratified flows, they are limited in Reynolds number $Re_\tau$ by at least four orders of magnitude when compared with typical SBL flows. At these scales, the LES technique offers the ability to simulate the large scales of these high-$Re_\tau$ flows with relative computational efficiency at the expense of not being able to explicitly resolve the fine-scale dynamics. This tradeoff results in a statistical dependence on grid resolution that becomes especially important for SBL studies (Khani & Waite 2014; Sullivan et al. 2016; Khani 2018; Dai et al. 2021; Maronga & Li 2022; Greene & Salesky 2023). In one of the few studies on coherent structures in the SBL that employed LES, Sullivan et al. (2016) utilised a fine grid spacing of $\Delta = 0.39$ m to simulate the SBL with varying surface cooling rates to induce increasing levels of static stability. They focused on the nature of localised coherent boundaries in the temperature field, and how these so-called microfronts act upon the surrounding flow. Through conditional averaging, the authors identified ring and hairpin vortices along the frontal boundaries that also lie within the energy-containing range of the turbulent flow. These frontal boundaries were also present in the DNS experiments of Gibbs et al. (2022), who similarly noted how their inclination angles flattened with height above the surface. Huang & Bou-Zeid (2013) additionally presented two-point correlation statistics on horizontal planes at varying heights within the SBL along with profiles of integral length scales. They concluded that turbulence becomes increasingly local with stability and that coherent structures are buoyantly suppressed in the vertical, leading to elongated features in the streamwise direction. Heisel et al. (2023) recently explored the role of stability on the organization of turbulence in UMZs and UTZs using the suite of LES from Sullivan et al. (2016). These authors found that under weak stability, the vertical thickness of UMZs and UTZs scale with distance from the wall, but become thinner and less dependent on $z$ as stability increases. These results indicate that deviations from the canonical logarithmic mean profiles of wind speed and temperature are related to decreased eddy sizes with higher stratification.

1.3. This study

While the existence and general features of turbulent coherent structures within stably stratified wall-bounded flows have recently been explored, at present there is a relative dearth of studies examining their role in modulating turbulence within the SBL. This study aims to close this knowledge gap by addressing the following key questions.

1. (i) How does stability impact the properties of LSMs within the SBL?

2. (ii) How does stability affect transport efficiencies of momentum and temperature?

3. (iii) How do coherent structures with the SBL contribute to these differences?

The rest of this paper is organised as follows. In § 2 we provide an overview of the LES code and cases considered. We discuss the impact of LES grid resolution with respect to the relevant scales of motion in the SBL in § 3. In § 4 we present our results on mean profiles and instantaneous fields in § 4.1, spectral analysis including spectrograms and linear coherence spectra in §§ 4.2 and 4.3, transport efficiencies in § 4.4, AM in § 4.5 and conditionally averaged fields in § 4.6. We conclude with a general discussion and interpretation of results, along with a future outlook in § 5.

2. Large-eddy simulation and cases

2.1. Large-eddy simulation code

In this study we employ the LES code summarised in Greene & Salesky (2023), which is described in more detail in Albertson & Parlange (1999) and Kumar et al. (2006). Throughout this study, we employ the notation of a tilde denoting a resolved (filtered) value, such that $\tilde {u}_i$ represents the filtered velocity vector where $i = 1,2,3$ represents the streamwise ($x$), spanwise ($\kern0.06em y$) and wall-normal ($z$) components, respectively. The code solves the filtered rotational form of the incompressible Navier–Stokes equations for momentum and potential temperature. Spatial derivatives are calculated pseudospectrally in the horizontal plane and via second-order centred finite differencing in the vertical, and the second-order Adams–Bashforth method is utilised for time integration. We utilise the Lagrangian-averaged scale dependent (LASD) subgrid-scale (SGS) model (Bou-Zeid, Meneveau & Parlange 2005) along with a wall model based on Monin–Obukhov similarity theory applied locally with test filtering at a scale twice the grid spacing to improve average stress profiles (Bou-Zeid et al. 2005). The forms of the Monin–Obukhov universal functions in the wall model are consistent with those of Beare et al. (2006). To force stable thermal stratification when simulating the SBL, we apply a lower thermal boundary condition using a prescribed surface temperature with a constant cooling rate (after e.g. Beare et al. (2006), see table 1). The upper boundary is impenetrable and stress free, and we apply a sponge layer in the upper $25\,\%$ of the domain after Nieuwstadt et al. (1993) to suppress the reflection of gravity waves. The LES code is parallelised in horizontal ($xy$) slabs using the message passing interface (Aoyama & Nakano 1999).

Table 1. Mean simulation properties for cases A–E averaged over the last physical hour of the simulation, including: the surface cooling rate $C_r$, SBL height $h$, low-level jet (LLJ) height $z_j$, ratio of the LLJ height to the SBL height $z_j/h$, surface friction velocity $u_{\tau 0}$, surface potential temperature scale $\theta _{\tau 0}$, Obukhov length $L$, global stability $h/L$, bulk Richardson number $Ri_B$, bulk SBL inversion strength $\Delta \langle \theta \rangle / \Delta z$, large-eddy turnover time $T_L = h/u_{\tau 0}$ and number of large-eddy turnover times within the last hour $nT_L$.

2.2. Cases

The cases presented are designed to simulate the SBL, and are based off the benchmark simulations of the Beare et al. (2006) LES intercomparison project. These five simulations are run with a constant value of surface cooling rate $C_r = -\partial \langle \theta _0 \rangle / \partial t$ spanning from values of $0.10 \le C_r \le 1.00$ K hr$^{-1}$ (table 1). All other parameters are held constant during the simulations and are summarised in table 2. Notably, these include setting the Coriolis parameter $f = 1.39 \times 10^{-4}$ s$^{-1}$ valid at 71 $^\circ$N latitude, and a temperature inversion $\varGamma _{inv} = 0.01$ K m$^{-1}$ for $z \ge z_{inv} = 100$ m.

Table 2. Values of various simulation parameters.

To balance computational efficiency with the numerical demands of accurately representing SBL processes within LES, these simulations were run on a grid consisting of $N = n_x n_y n_z = 192 \times 192 \times 384$ total nodes for nine physical hours. This grid aspect ratio of four was chosen to optimize the ability of this model to capture the buoyant suppression of eddies in the vertical and improves upon the resolution presented by Greene & Salesky (2023) by following a similar strategy as Kimura & Sullivan (2024). As will be discussed in § 4, the domain size ($800 \times 800 \times 400$ m$^3$) is large enough to resolve LSMs but not VLSMs for near-neutral stratification, and was chosen in order to resolve the small-scale structure of the SBL as finely as possible. All simulations herein were performed on the National Center for Atmospheric Research (NCAR) peta-scale supercomputer Derecho (Computational and Information Systems Laboratory 2023).

As in Greene & Salesky (2023) and unless otherwise specified with subscripts, we use angle brackets $\langle {\cdot } \rangle = \langle {\cdot } \rangle _{xyt}$ to denote averaging in horizontal planes and over the final hour of simulation (consistent with Beare et al. 2006), which corresponds to 4.92–7.27 large-eddy turnover times (table 1). Quantities with a prime indicate fluctuations about the mean, e.g. $\tilde {u}' = \tilde {u} - \langle \tilde {u} \rangle$. Additional analysis of these periods indicate they are quasi-stationary (not shown), and to improve statistical robustness, we implement linear detrending when calculating second-order turbulent parameters.

We determine the SBL depth $h$ as in Beare et al. (2006). Namely, $h$ is the normal distance from the wall where the mean stress profile $u_{\tau }(z)$ falls to a value equal to 5 % of its surface value $u_{\tau 0}$ and then linearly interpolated to zero by dividing by a factor of 0.95. Although many definitions of SBL depth exist in the literature, this one is most consistent with other LES studies of the SBL and still retains the majority of the low-level jet (LLJ) peak within the SBL for all but case C (table 1).

Here we define the mean stress profile by combining the streamwise and spanwise components as

(2.1)\begin{equation} u_{\tau}^2(z) = ( \langle \tilde{u}'\tilde{w}' + \tau_{xz} \rangle^2 + \langle \tilde{v}'\tilde{w}' + \tau_{yz} \rangle^2 )^{1/2}, \end{equation}

where $\tau _{xz}$ and $\tau _{yz}$ are the SGS contributions to the kinematic momentum flux such that $\tau _{xz} = \widetilde {uw} - \tilde {u}\tilde {w}$ and $\tau _{yz} = \widetilde {vw} - \tilde {v}\tilde {w}$. Other parameters defined in table 1 include the potential temperature scale $\theta _{\tau 0} = -\langle \tilde {\theta }' \tilde {w}' + {\rm \pi}_3^{\theta } \rangle _0 / u_{\tau 0}$, where ${\rm \pi} _3^{\theta }$ is the SGS heat flux ($\pi _3^{\theta } = \overline {\theta w} - \tilde {\theta }\tilde {w}$), the Obukhov length $L = u_{\tau 0}^2 \langle \theta _0 \rangle /\kappa g \theta _{\tau 0}$ which depends on the von Kármán constant $\kappa = 0.4$, the mean lapse rate $\Delta \langle \theta \rangle / \Delta z$ between the top of the SBL and lowest grid point, and the bulk Richardson number $Ri_B$ defined as

(2.2)\begin{equation} Ri_B = \frac{\dfrac{g}{\langle \theta_0 \rangle} \dfrac{\Delta \theta}{\Delta z}}{\left(\dfrac{\Delta \langle u \rangle}{\Delta z}\right)^2 + \left(\dfrac{\Delta \langle v \rangle}{\Delta z}\right)^2}, \end{equation}

where the differences are likewise calculated between $z=h$ and $z=\varDelta _z/2$.

For the results presented in § 4 and unless otherwise stated, we utilise the full volumetric fields for analysis. In § 4.5 we additionally employ output from a virtual tower centred in the domain at $(x_0,y_0,z) = (L_x/2, L_y/2, z)$ that emulates measurements from eddy-covariance systems at each domain height $z$ sampling at 50 Hz for time-series analysis (see Salesky & Anderson 2018; Greene & Salesky 2023).

To account for the flow rotation due to the Coriolis force, unless otherwise specified all statistics are computed with rotated volumetric fields in the coordinate system $(x',y',z) = (x'(z), y'(z), z)$ such that $x'$ is rotated about the $z$ axis to where $\langle \tilde {v} \rangle =0$ at each height. We accomplish this via a bilinear interpolation of each rotated coordinate system in horizontal slabs. To preserve the high-frequency representation of the flow that may be lost due to the bilinear interpolation, we first interpolate the unrotated volumetric fields in spectral space to increase the horizontal wavenumber resolution by a factor of two. As will be seen in § 4, this rotated coordinate system is closely aligned with coherent features in the simulations and ensures more physically representative interpretations of statistical and spectral analyses.

3. Resolution impacts on simulated SBL dynamics

Before presenting the results regarding coherent structures, it is worth discussing the performance of our LES model with regards to resolving the fine-scale dynamics under increasingly stable stratification. In this context it is useful to consider the Ozmidov scale $L_O$ (Dougherty 1961; Ozmidov 1965), which is defined based on the mean TKE dissipation rate $\epsilon$ and the Brunt–Väisälä frequency $N$ as

(3.1)\begin{equation} L_O = \sqrt{ \frac{\epsilon}{N^3} }. \end{equation}

The Ozmidov scale has the physical interpretation of being the largest eddy size unaffected by buoyancy, and has been shown to be the characteristic size of momentum transporting eddies within the SBL (Bou-Zeid et al. 2010; Li, Salesky & Banerjee 2016). Li et al. (2016) found that $L_O$ actually constrains the energy-production end of the inertial subrange under increasing stability. Sullivan et al. (2016) discuss the importance of explicitly resolving this scale when using LES, as typical SGS models are dissipative and are not designed to effectively emulate the small-scale overturnings when $L_O < \varDelta _f$. Because they utilise a fine mesh grid of $\varDelta _f = 0.39$ m, this is only an issue close to the surface and near the top of the SBL in their LES. Huang & Bou-Zeid (2013) address this explicitly, noting that when the LES filter scale lies within the inertial subrange, it is imperative for the SGS model to correctly drain TKE from the resolved to SGS scales to produce correct fluxes. They further argue that the LASD SGS model has proven capable of this difficult task, and has superior performance to traditional Smagorinsky–Lilly or scale-invariant models (Bou-Zeid et al. 2005). We therefore can maintain reasonable confidence in the ability of our LES model set-up to reproduce accurate total fluxes in the SBL. As will be discussed further, some caution will be necessary when considering spectral analysis of only the resolved velocity and temperature fields, as it is not always possible to include the SGS contributions with these analyses.

Profiles of $L_O$ relative to the effective filter width $\varDelta _f$ are presented in figure 1(a). For these cases, we define TKE dissipation rate for LES as the energy flux across the filter scale (Pope 2000) based on the subgrid stress tensor $\tau _{ij}$ and the filtered strain rate tensor $\tilde {S}_{ij}$ such that $\epsilon = -\langle \tau _{ij} \tilde {S}_{ij} \rangle$. We obtain $\tau _{ij}$ directly from the LASD SGS model and compute its product with $\tilde {S}_{ij}$ during the simulation to be output with the rest of the standard variables. From figure 1(a), it is apparent that $L_O$ is explicitly resolved for all cases A–D throughout the entire SBL, and for $z/h < 0.6$ in case E. When considering the ratio of $L_O$ relative to the vertical grid spacing $\varDelta _z$, all cases are explicitly resolved for the entire SBL depth, and by a factor of 10 in cases A–D in the lower half of the SBL (not shown).

Figure 1. Profiles of the (a) ratio of the Ozmidov scale $L_O$ to the LES characteristic filter size $\varDelta _f$ (in log coordinates), and ratios of subgrid (b) momentum and (c) heat flux to the total (resolved plus SGS) fluxes for all simulations A–E.

With such a reliance on the SGS model when simulating the SBL, it is also useful to consider the ratios of the SGS fluxes to the total (resolved plus SGS) fluxes. Profiles of these ratios for momentum and heat fluxes are included in figures 1(b) and 1(c), respectively. As expected, this ratio is large across all cases close to the wall, but in the middle of the SBL the ratios increase monotonically with stability. The ratios of SGS momentum and heat fluxes for all cases A–E remain at or below 30 % for $0.05 < z/h < 0.9$, with these ratios below 10 % in the middle of the SBL for cases A–C. Note that the SGS momentum flux ratio grows large in proximity to the LLJ at $z/z_j \approx 1$ as the denominator crosses zero for most cases at this height and is not necessarily indicative of an increased stress on the subgrid model.

We recognise that the resolution of $L_O$ by an order of unity is not a perfect global test on our model's resolution. We have found, however, that one-dimensional spectra of streamwise and vertical velocities close to the ground in all cases produce at least a decade of an inertial subrange (not shown), and the model filter width is located within the inertial subrange for each simulation. These results, in combination with the LASD SGS model employed in this study, yields confidence that the suite of simulations herein are of sufficiently high resolution to represent the turbulent processes of interest as described in § 1.3.

4. Results

4.1. Mean profiles and instantaneous fields

Profiles of mean quantities from cases A–E over the final hour of simulation (between hours 8–9, following Beare et al. 2006) are included in figure 2. The mean wind speed profiles $U_h/G = \sqrt {\langle \tilde {u} \rangle ^2 + \langle \tilde {v} \rangle ^2} / \sqrt {U_g^2 + V_g^2}$ (figure 2a) generally increase beyond the background geostrophic wind in the middle of the SBL, with the maximum speed at $z = z_j$ increasing with stability. The mean potential temperature profiles $\varTheta = \langle \tilde {\theta } \rangle$ (figure 2b) display a strong sensitivity to the surface cooling rate, as mean lapse rates increase monotonically from cases A–E for $0.2 \lesssim z/h \lesssim 0.8$. The normalised root-mean-square velocity, $u_{rms} = \sqrt {0.5 (\langle \tilde {u}'^2 \rangle + \langle \tilde {v}'^2 \rangle + \langle \tilde {w}'^2 \rangle )} = e^{1/2}$, where $e$ is the TKE (figure 2c), generally decreases with stability throughout the depth of the SBL likely due in part to the buoyant suppression of TKE with increasing stability.

Figure 2. Mean profiles from all simulations A–E of (a) horizontal wind speed $U_h = \sqrt {\langle \tilde {u} \rangle ^2 + \langle \tilde {v} \rangle ^2}$ as a fraction of the geostrophic wind vector magnitude $G$ (see table 2), (b) potential temperature $\varTheta = \langle \tilde {\theta } \rangle$ differences from the lowest grid point, (c) root-mean-square resolved velocity $u_{rms} = \sqrt {0.5 (\langle \tilde {u}'^2 \rangle + \langle \tilde {v}'^2 \rangle + \langle \tilde {w}'^2 \rangle )}$, (d) total momentum flux $u_\tau ^2$ (2.1), (e) total potential temperature flux and ( f) gradient Richardson number $Ri_g$, with a reference line at $Ri_g = 0.2$. Statistics are calculated using the final hour of each simulation.

Profiles of the gradient Richardson number (figure 2g),

(4.1)\begin{equation} Ri_g = \frac{g}{\varTheta_0} \frac{\partial \varTheta}{\partial z} \left[ \left(\frac{ \partial \langle \tilde{u} \rangle}{\partial z} \right)^2 + \left(\frac{ \partial \langle \tilde{v} \rangle}{\partial z} \right)^2 \right]^{{-}1}, \end{equation}

in general increase monotonically with surface cooling rate for a given height. The weakly stable cases are largely within the subcritical regime associated with Kolmogorov turbulence, $Ri_g < 0.2$, as identified by Grachev et al. (2013), whereas simulation E lies above this threshold for $z/h > 0.6$. Finally, the mean profiles of non-dimensional total (resolved plus SGS) momentum and heat flux (figure 2d,e) generally collapse, although the weakly stable cases A and B are more linear than the rest as they are closer to neutral stratification. The irregularities in the lowest grid points for these profiles can be attributed to the wall model. These mean flux profiles are in reasonable agreement with the semi-empirical formulations presented by Nieuwstadt (1984) and Sorbjan (1986) that were validated against data from the 1973 Minnesota experiment (Caughey, Wyngaard & Kaimal 1979).

To visually highlight coherent structures within these SBL flows, in figure 3 we present horizontal and vertical cross-sections of the instantaneous fluctuating streamwise velocity and potential temperature fields from simulations A, C and E. Inspection of the horizontal ($x$–$y$) cross-section of $\tilde {u}'/u_{\tau 0}$ located at $z/h=0.05$ (figure 3a–c) indicate elongated streaks of high and low momentum that decrease in size and magnitude with stability. In simulation A (figure 3a) these streaks are $\approx 0.5$–$1.5h$ in length and are rotated with respect to the geostrophic wind (table 2), which is consistent with the streamwise extent of LSMs under neutral and unstable stratification. This is analogous to how HCRs in the CBL are typically rotated $\approx$15–20$^\circ$ to the left of the geostrophic wind due to surface drag and momentum flux divergence (e.g. Salesky et al. (2017), and references therein). Ansorge & Mellado (2014) include discussion of these features within stably stratified turbulent flows, but is otherwise beyond the scope of this study. With increasing stability, in cases C and E (figure 3b,c) the velocity field is organised into fine ribbons of weaker fluctuations than in case A, and areas of locally similar magnitudes are less well defined.

Figure 3. Instantaneous cross-sections from simulations A (left column), C (middle column) and E (right column) including: (a–c) streamwise velocity perturbations and (g–i) potential temperature perturbations in the unrotated $x$–$y$ plane at $z/h=0.05$, and (d–f) streamwise velocity perturbations and (j–l) potential temperature perturbations in the rotated $x'-z$ plane as indicated by the superimposed axes in (a–c) and (g–i). The horizontal lines in (d–f) and (j–l) denote the heights above which $Ri_g \ge 0.2$. The thin angled lines in (d) and (j) annotate the presence of a wall-attached coherent structure.

Unlike streamwise velocity, the horizontal cross-sections of potential temperature fluctuations (figure 3g–i) do not demonstrate corresponding elongated streaks. This result is strikingly different than what is found in the CBL (Salesky et al. 2017), where the vertical velocity and potential temperature fields are visually analogous. Instead, the potential temperature field is composed of a patchy network of warm and cold pockets whose magnitudes depend on stability. In case A there are a few locations where seemingly organised patches of temperature overlap with the long streaks in velocity, e.g. around $(x/h \approx 1, y/h \approx 1)$. It is clear, however, that the frequency of these overlapping patterns is lower in cases C and E than for case A. Here we note that in addition to $\tilde {u}'$ and $\tilde {\theta }'$, cross-sections of $\tilde {w}'$ (not shown) demonstrate generally weak vertical motions without significant organization and, thus, are omitted.

It is apparent from the vertical cross-sections in streamwise coordinates (figure 3d–f) that the elongated velocity streaks extend into the vertical under weak stability. For example, in case A between $0.1 < x'/h < 1.1$, a region of $\tilde {u}'/u_\tau > 0$ extends from the surface up to $z/h \approx 0.25$. These dimensions ($\Delta x'/h \approx 1$, $\Delta z/h \approx 0.25$) correspond to an inclination angle of $\gamma = \arctan (0.25/1) = 14.0^\circ$ with respect to the surface, which agrees well with the values based on two-point correlation statistics as presented by Gibbs et al. (2022) under weakly stable stratification. With increasing stability, however, analogous structures become decreasingly prominent within the flow (figure 3i). The vertical cross-sections of potential temperature fluctuations (figure 3j–l) highlight similar features as those from the streamwise velocity fluctuations. In case A there are elongated regions of high and low perturbations with sharp boundaries in between, which is highly reminiscent of the temperature microfronts presented by Sullivan et al. (2016). It is apparent by examining the warm anomaly attached to the surface around $x'/h \approx 0.1$ (corresponding to the one discussed in figure 3d) that the temperature structures do not necessarily incline at the same relative angles as for momentum. As will be discussed further, this eludes to differing mechanisms for vertical transport of heat and momentum under stable stratification. The overall spatial correlation between the momentum and temperature fluctuations in this region align with a sweep of relatively warmer, high-momentum fluid moving towards the surface.

From an analysis of the instantaneous fields presented in figure 3, it is clear that buoyancy acts to suppress vertical organization more strongly than in the horizontal. This has been previously documented (e.g. Huang & Bou-Zeid 2013; Chinita, Matheou & Miranda 2022) and will be important to consider when analysing the results in the following sections.

4.2. Spectrograms

One common method for identifying the presence of coherent structures within turbulent flows is through the analysis of spectrograms (Hutchins & Marusic 2007a; Mathis et al. 2009a; Anderson 2016; Baars, Hutchins & Marusic 2017; Salesky & Anderson 2018). Spectrograms are premultiplied power spectra presented as functions of both wavelength and height above the surface, and evidence for LSMs exists when an outer peak is present at large wavelengths. Included in figure 4 are spectrograms for cases A–E of streamwise and vertical velocities, potential temperature and cospectra of $\langle {\tilde {u}'\tilde {w}'} \rangle$ as well as $\langle {\tilde {\theta }'\tilde {w}'} \rangle$.

Figure 4. Premultiplied spectrograms from simulations A–E (columns) for (a–e) streamwise velocity,( f–j) vertical velocity, (k–o) potential temperature, as well as cospectra of (p–t) $\langle {\tilde {u}'\tilde {w}'} \rangle$ and (u–y) $\langle {\tilde {\theta }'\tilde {w}'} \rangle$. Each is plotted versus streamwise wavelength $\lambda _x$ and wall-normal height $z$ normalised by the SBL depth $h$. Horizontal lines at $\lambda _x = h/4$ indicate the cutoff frequency utilised in the decoupling procedure outlined in § 4.5 that roughly separates the inner and outer peaks (where they exist). Vertical lines are plotted at heights where $Ri_g \ge 0.2$.

Inspection of figure 4 reveals distinct inner and outer peaks in case A across all parameters except for vertical velocity. In case A (figure 4a, f,k,p,u), the outer peak is located approximately at $z/h \approx 0.05$ and $\lambda _x / h \approx 1$, which is as expected within the logarithmic region of the wall-bounded flow at streamwise wavelengths approximately scaling with the flow depth. The wavelength $\lambda _x / h \approx 1$ associated with these outer peaks in case A is also consistent in scale with the velocity streaks in figure 3(a). With this evidence we can reasonably conclude that LSMs are present under weak stability. The outer peak scales are also slightly smaller than those reported by, e.g. Baars et al. (2017) for neutrally stratified channel flow, but roughly an order of magnitude smaller than those in the CBL reported by Salesky & Anderson (2018). These differences may likely be due to the lack of VLSMs in the domain considered, so energy peaks at the scale of individual coherent structures instead of a collective superstructure.

For increasing stability, the outer peaks in all of the spectrograms attenuate until they disappear entirely. This behaviour, specifically in streamwise velocity (figure 4a–e), is also in contrast with the findings of Salesky & Anderson (2018), who found that within the CBL, the peak distinctly moved toward smaller wavelengths until it merged with the inner peak with increasing instability. This is undoubtedly the signature of buoyant suppression of vertical motions, which has been shown to act at the large scales (García-Villalba & del Álamo 2011) so that large eddies do not traverse the full depth of the SBL. This behaviour is also consistent with the work of Kaimal et al. (1972), who noted decreasing energy at large scales under increasing stability. Recall from the discussion in § 3 and from the results of Li et al. (2016) that the Ozmidov scale is a characteristic eddy size within the SBL that denotes the beginning of the inertial subrange. Since $L_O$ increases with stability, this implies that the beginning of the inertial subrange shifts towards higher wavenumbers as energy carried by large eddies is damped by buoyancy. This notion is further supported by the lack of an outer peak in the vertical velocity across all cases, although a noticeable ridge does extend from the surface towards larger scales and heights that decreases with increasing stability. The combined effect of these results is evident in the $\langle {\tilde {u}'\tilde {w}'} \rangle$ cospectra (figure 4p–t), indicating a decreasing correlation between $u$ and $w$ with increasing stability. The outer peak in the potential temperature spectrogram $k_x \varPhi _{\theta \theta } / \theta _{\tau 0}^2$ (figure 4k–o) notably is higher in the SBL and occurs at longer wavelengths than those for $u$ for each case. For example, in case B (figure 4l) this peak is centred on $z/h \approx 0.3, \lambda _x/h \approx 1$. Moreover, an outer peak in the potential temperature spectrogram persists through at least case C in similar fashion to the $u$ spectrograms. There also appears to be an outer peak in the $\langle {\tilde {\theta }'\tilde {w}'} \rangle$ cospectra $k_x \varPhi _{\theta w}/\theta _{\tau 0} u_{\tau 0}$ (figure 4u–y) through cases A–C. These differences in momentum and scalar spectrograms indicate underlying differences in their respective transports, which will be discussed further in § 4.4.

Recall from figure 2( f) that the majority of the SBL for cases A–D are in the range $Ri_g < 0.2$, whereas the upper third of the SBL in case E becomes supercritical. Turbulence in these regions is not necessarily persistent in time and is heavily suppressed by buoyancy, thereby inhibiting the development of elongated coherent structures as observed close to the wall. Studies on intermittent turbulence in these regimes typically rely on DNS (e.g. Deusebio et al. 2011; Ansorge & Mellado 2014, 2016; Deusebio, Augier & Lindborg 2014a; Deusebio et al. 2014b; Deusebio, Caulfield & Taylor 2015), as the LES technique can struggle to adequately represent the local nature of turbulence production, transition and dissipation under strong stratification. Herein we therefore take extra caution when interpreting further results from these regions of supercritical stability.

From these spectrograms, it is apparent that buoyancy acts strongly to attenuate vertical motions at large streamwise wavelengths, resulting in turbulence that is increasingly local with increasing stability. There is evidence that velocity and potential temperature organise into large-scale coherent structures through at least case C (recall from table 1 that $C_r = 0.33$ K h$^{-1}$, $h/L=1.83$) based on the presence of outer peaks. In § 4.3 we explore further how these fields are affected by stability across scales.

4.3. Linear coherence spectra

In addition to spectrograms, another method of diagnosing the relevant scales affected by coherent structures is through computation of the linear coherence spectrum (LCS; Baars et al. 2017). The LCS is a measure of the linear coupling between variables across scales, and is defined as

(4.2)\begin{equation} \gamma_{uu}^2(z,z_R;\lambda_x) = \frac{| \langle \hat{u}(z;\lambda_x) \hat{u}^*(z_R;\lambda_x) \rangle_{yt} |^2}{\langle |\hat{u}(z;\lambda_x)|^2 \rangle_{yt} \langle |\hat{u}(z_R;\lambda_x)|^2 \rangle_{yt}}, \end{equation}

where $\hat {u}(z;\lambda _x) = \mathcal {F} \{ u(x,z) \}$ is the Fourier transform of $u(x,z)$ along the streamwise dimension with the asterisk $*$ denoting its complex conjugate, $z_R$ represents a constant reference height above ground level for comparison against all heights $z$ and $|{\cdot }|$ refers to the modulus of a complex value. In this context, $\gamma _{uu}^2$ falls within the range $\gamma _{uu}^2 \in [0, 1]$, and can be interpreted as the squared value of the correlation coefficient at a specific scale $\lambda _x$ between fluctuating values of $u$ at two different heights $z$ and $z_R$. An example of this is included in figure 5, where we calculate $\gamma _{uu}^2$ and $\gamma _{\theta \theta }^2$ using a reference height $z_R = \varDelta _z/2$ as the lowest grid point in each simulation. The strongest coupling across all simulations and parameters is noticeably found at horizontal wavelengths of $O(h)$, as was found with the outer peaks in the premultiplied spectrograms in the previous section. Moreover, the vertical extent of the LCS peaks diminishes with increasing stability. For example, $\gamma _{uu}^2$ for case A extends beyond $z/h=0.1$ whereas by case C the contour for $\gamma _{uu}^2=0.1$ only extends to $z/h \approx 0.04$. The $\gamma _{uu}^2$ and $\gamma _{\theta \theta }^2$ peaks for case A (figure 5a, f) can also be attributed to the coherent features identified in the instantaneous fields that extend from the surface up into the outer region of the flow (figure 3a,e). We note here that due to vertical resolution limitations using a wall-modelled LES, the contours near the surface do not provide significant amounts of information at higher stabilities. Baars et al. (2017) argue that a 1 : 1 slope of these peaks in log-log coordinates, specifically in streamwise velocity under near-neutral stratification (case A), is consistent with the attached-eddy hypothesis (Townsend 1976) across a self-similar hierarchy of scales. Analysis of these cases in the framework of the attached-eddy hypothesis is beyond the scope of this study, but certainly warrants further investigation, ideally with higher resolutions close to the wall.

Figure 5. Linear coherence spectra for (a–e) $u$, ( f–j) $w$ and (k–o) $\theta$ for cases A–E (columns) calculated with the reference point $z_R$ as the lowest grid point and plotted against non-dimensional wavelength and wall-normal distance. The horizontal line in each panel is the same as in figure 4.

It is also possible to define (4.2) for two independent variables at the same height, which provides information on the coupling of two parameters across scales. For example, the LCS between $u$ and $w$ at wall-normal distance $z$ would be determined as $\gamma _{uw}^2 = \gamma _{uw}^2(z,z;\lambda _x)$. The resulting values of $\gamma _{uw}^2$ and $\gamma _{\theta w}^2$ are included in figure 6 for cross-sections at constant heights $z/h$ within the SBL. Near the top of the surface layer ($z/h=0.1$, figure 6a) it is apparent that the coupling between $u$ and $w$ is strongest for $\lambda _x / h \ge 0.5$ in all cases. With increasing stability, the maximum in $\gamma _{uw}^2$ monotonically shifts towards smaller wavelengths and decreases in magnitude from a value of $\approx 0.4$ in case A to $\approx 0.3$ in case E. The magnitude of $\gamma _{uw}^2$ for wavelengths larger than their respective peaks also decreases monotonically with increasing stability. For vertical transport of potential temperature at this height (figure 6c), the impacts from buoyancy are noticeable in the attenuation of $\gamma _{\theta w}^2$ with increasing stability at larger scales ($\lambda _x/h > 0.5$). In general, we observe that the peaks at this height for $\gamma _{uw}^2$ are larger than those for $\gamma _{\theta w}^2$, but at small wavelengths $\lambda _x/h < 0.1$ this trend is reversed.

Figure 6. Cross-sections of LCS from cases A–H at constant heights: (a,c) $z/h=0.1$ and (b,d) $z/h=0.25$ for (a,b) $\gamma _{uw}^2$ and (c,d) $\gamma _{\theta w}^2$. Vertical lines at $\lambda _x / h = 0.25$ are included for reference.

Higher in the SBL at $z/h = 0.25$ (figure 6b,d), we largely observe the same patterns in $\gamma _{uw}^2$ and $\gamma _{\theta w}^2$ as those for $z/h = 0.1$ but with slightly smaller peak values at high stability. These results again reflect how buoyancy suppresses vertical transport of both momentum and heat at large scales, but the differences in $\gamma _{uw}^2$ and $\gamma _{\theta w}^2$ at fine scales allude to differing transport mechanisms. This topic is discussed in further detail in § 4.4.

4.4. Transport efficiency

In § 4.2 we identified the existence of an outer peak in the premultiplied spectrograms in cases A–C, and in § 4.3 found enhanced linear coupling at the scales of these outer peaks. Previous studies of the CBL have shown that turbulent transports of momentum and scalars become increasingly dissimilar with increasing instability (Li & Bou-Zeid 2011; Dupont & Patton 2012; Patton et al. 2016; Salesky 2023), and Salesky et al. (2017) were able to connect this breakdown of Reynolds’ analogy to varying modes of convective organization. However, the relationship between momentum and heat transport in stably stratified turbulent shear flows remains relatively unexplored. To study these effects, it is useful to consider the partitioning of turbulent fluxes into contributions by individual positive and negative fluctuations in either term. This technique is known as quadrant analysis (also conditional sampling; see Wallace (2016), and references therein), and is outlined as follows. Using the resolved vertical momentum flux $\langle {\tilde {u}'\tilde {w}'} \rangle$ as an example, we define the four quadrants as

• (i) Quadrant I: $\tilde {u}' > 0, \tilde {w}' > 0$,

• (ii) Quadrant II: $\tilde {u}' < 0, \tilde {w}' > 0$,

• (iii) Quadrant III: $\tilde {u}' < 0, \tilde {w}' < 0$,

• (iv) Quadrant IV: $\tilde {u}' > 0, \tilde {w}' < 0$.

With this definition, quadrants II and IV are respectively referred to as ejections and sweeps. The quadrants for potential temperature flux are defined likewise by replacing $u$ with $\theta$, and for stable thermal stratification, quadrants II and IV also refer to the downgradient direction. The turbulent transport efficiencies based on these quadrants are defined based on the fraction of the total flux occurring in the downgradient direction (Wyngaard & Moeng 1992; Li & Bou-Zeid 2011; Salesky et al. 2017). For the SBL where $\partial U_h / \partial z > 0$ and $\partial \varTheta / \partial z > 0$, the transport efficiencies for momentum and heat are defined (adopting the notation of the present paper) as

(4.3)\begin{equation} \eta_{uw} = \frac{\langle {\tilde{u}'\tilde{w}'} \rangle}{\langle {\tilde{u}'\tilde{w}'} \rangle^{II} + \langle {\tilde{u}'\tilde{w}'} \rangle^{IV}} \end{equation}

and

(4.4)\begin{equation} \eta_{\theta w} = \frac{\langle {\tilde{\theta}'\tilde{w}'} \rangle}{\langle {\tilde{\theta}'\tilde{w}'} \rangle^{II} + \langle {\tilde{\theta}'\tilde{w}'} \rangle^{IV}}. \end{equation}

Here the superscripts II and IV refer to the individual quadrant contributions to the total flux from specifically quadrants II and IV. Note that we intentionally neglect the SGS flux contributions in (4.3) and (4.4) for use with LES output, as the quadrant assignment for, e.g. $\tau _{xz}$ is not directly discernible from the signs of $\tilde {u}'$ and $\tilde {w}'$.

Profiles of these transport efficiencies along with their ratio are displayed in figure 7. There is a modest dependence on stability for both $\eta _{uw}$ and $\eta _{\theta w}$ that is most apparent in the profile of their ratio (figure 7c). The ratio $\eta _{uw}/\eta _{\theta w}$ is close to unity for $0.2 < z/h < 0.6$ for all cases A–E, and generally does not vary with stability in the layer $z/h < 0.2$. For $z/h > 0.6$, this ratio has a strong dependence on stability as $\eta _{\theta w}$ does not decrease as sharply with height for cases C–E as compared with cases A and B. In general, this decrease with height in $\eta _{uw}$ is mostly consistent across cases for this region of the SBL. In the middle of the SBL, there is a minor dependence on both $\eta _{uw}$ and $\eta _{\theta w}$ with stability, for example, at $z/h=0.4$ both transport efficiencies are highest for case A and lowest for case E. Closer to the surface ($z/h=0.05$), this trend is somewhat reversed. Regardless of these differences, the transport efficiencies in both momentum and heat are markedly lower than those reported in both observed and simulated CBLs at weak instability (e.g.  Li & Bou-Zeid 2011; Salesky et al. 2017). Therefore, it is apparent that even weak stratification plays a strong role in inhibiting the flow's ability to vertically redistribute momentum or heat.

Figure 7. Profiles of transport efficiencies (a) $\eta _{uw}$ (4.3), (b) $\eta _{\theta w}$ (4.4) and (c) their ratio $\eta _{uw}/\eta _{\theta w}$ for cases A–E.

These differences in transport efficiencies can be traced to changes in turbulent motions from each quadrant I–IV as displayed in figure 8. Plotted here are the individual quadrant fractions $Q_{uw}^k$ and $Q_{\theta w}^k$, which are defined as

(4.5)\begin{equation} Q_{uw}^k = \frac{| \langle {\tilde{u}'\tilde{w}'} \rangle^k |}{\sum | \langle {\tilde{u}'\tilde{w}'} \rangle^k |} \end{equation}

and

(4.6)\begin{equation} Q_{\theta w}^k = \frac{| \langle {\tilde{\theta}'\tilde{w}'} \rangle^k |}{\sum | \langle {\tilde{\theta}'\tilde{w}'} \rangle^k |}, \end{equation}

where $k \in \{ \mathrm {I, II, III, IV}\}$ represents the individual quadrant contributions to the absolute sum.

Figure 8. Individual quadrant fractions (a–d) $Q_{uw}^k$ (4.5) and (e–h) $Q_{\theta w}^k$ (4.6) for cases A–E. Vertical lines at $Q^k=0.25$ are included for reference.

Recalling that quadrants II and IV in the momentum flux (figure 8b,d) denote ejections and sweeps, respectively, it is apparent that motions in these quadrants dominate the total flux profile with quadrant fractions $Q_{uw}^{II}, Q_{uw}^{IV} > 0.25$ for all cases A–E. The fraction of ejections remains roughly constant with height for each case for $0.1 < z/h < 0.8$, whereas the sweeps decrease with height in this range. In the range $0.1 < z/h < 0.9$, the fraction of ejections also tends to decrease gradually with stability, and the upper bound on this range also decreases with stability as the profiles break towards $Q_{uw}^{II}=0.25$ at progressively lower heights. In the upper half of the SBL, the decreasing ejections with stability are primarily compensated for in the countergradient motions of quadrant I and by sweeps (figure 8a,d). The differences in quadrant III (figure 8c) are comparatively smaller with changes in stability, indicating that changes in transport efficiency largely depend on how positive vertical motions interact with relatively high or low streamwise momentum parcels within the SBL. At the top of the SBL ($z/h \approx 1$), all four quadrants reach equilibrium with an even distribution of $Q_{uw}^k = 0.25$.

The heat flux profiles display a somewhat different pattern, however, with quadrants II and IV dominating the contributions at all levels (figure 8f,h). Otherwise, the general trends in the heat flux quadrant fractions are largely similar to those of momentum fluxes: ejections (upwelling relatively cold parcels) are relatively constant with height whereas sweeps decrease with height, reaching a minimum around $z/h \approx 0.85$. In the middle of the SBL around $z/h \approx 0.5$, the thermal ejections decrease in fraction with stability, which is accounted for by an increase in countergradient motions in quadrants I and IV at this level whereas quadrant III is relatively invariant (similar to those for momentum).

We explore the contributions to each quadrant further by examining joint probability density functions (JPDFs) and their corresponding covariance integrands. These two are related by considering two random variables $a$ and $b$ such that their covariance $\langle a'b' \rangle = \int _{-\infty }^{\infty } abP(a,b) \,\mathrm {d}a \,\mathrm {d}b$ depends on their JPDF, $P(a,b)$ (Wallace 2016). By considering both of these quantities we can identify the pairs of, e.g. $u'$ and $w'$ that contribute most strongly to each quadrant fraction $Q_{uw}^k$.

The JPDFs and their covariance integrands of $uw$ and $\theta w$ are included in figure 9 for cases A, C and E at heights of $z/h = 0.1$ and $0.5$. Corresponding to the decrease in $uw$ ejections (figure 8b), it is apparent from figure 9(a,e,i) that the quadrant II peak in $uwP(u,w)$ shifts towards smaller values of $u'/u_{\tau 0}$ and $w'/u_{\tau 0}$ with increasing stability. Additionally, the extent of $uwP(u,w)$ in quadrant IV reaches towards larger values of $u'/u_{\tau 0}$ with increasing stability, which accounts for the increase in $Q_{uw}^{IV}$ at $z/h=0.1$. These changes are even more drastic at $z/h=0.5$ (figure 9b, f,j), where the quadrant I peak of $uwP(u,w)$ increases in magnitude with stability as the quadrant II peak weakens and the overall distribution of $u'$ and $w'$ becomes more evenly distributed in all four quadrants.

Figure 9. Joint PDFs and their premultiplied covariance integrands selected at (first and third columns) ${z/h=0.1}$ and (second and fourth columns) $z/h=0.5$ for cases (row one) A, (row two) C and (row three) E. Joint PDFS are shaded in greyscale, and the covariance integrands are overlaid in red and blue contours. Included are (first two columns) $uw$ and (last two columns) $\theta w$.

The joint distribution of $\theta$ and $w$ at $z/h=0.1$ does not change drastically with stability (figure 9c,g,k), which is expected given their quadrant fractions observed in figure 8. In the middle of the SBL, however, these distributions are extremely sensitive to increasing stability (figure 9d,h,l). The spread in values of $w'/u_{\tau 0}$ decreases markedly with stability, most notably in quadrants II and IV in $\theta wPDF(\theta,w)$, and there is not a corresponding decrease in the spread of $\theta '/\theta _{\tau 0}$. Additionally, the peaks of $\theta wPDF(\theta,w)$ tend towards smaller values of $|\theta '/\theta _{\tau 0}|$ and $|w'/u_{\tau 0}|$ in all quadrants, which corresponds to overall weaker and less efficient transport of heat in the middle of the SBL (figure 7b).

It is apparent from the turbulent transfer efficiencies (figure 7), individual quadrant fractions (figure 8) and JPDFs (figure 9) that momentum and heat are transported differently as stability increases. Under weakly stable stratification (case A), our results generally match those from Li & Bou-Zeid (2011) and Salesky et al. (2017) under weakly unstable conditions. Therefore, it is likely that coherent turbulent structures such as hairpin vortices exist in case A (Adrian 2007), but increasing stratification flattens motions into largely horizontal features. Vertical motions become more localised and contributions from countergradient fluxes reduce the overall efficiencies in turbulent transport of both momentum and heat. This is also consistent with the small-scale circulations around microfronts in the SBL as observed by Sullivan et al. (2016) from their high-resolution LES.

4.5. Amplitude modulation

To further examine how LSMs affect the smaller scales within the SBL, in this section we perform the decoupling procedure outlined by Mathis et al. (2009a) and more recently by Salesky & Anderson (2018). The decoupling procedure as implemented with single-point correlations is summarised as follows.

First we consider two random variables $a=a(z;t)$ and $b=b(z;t)$. We are interested in computing the extent to which the large scales of signal $b$ at height $z$ modulate the small-scale amplitude of signal $a$ also at height $z$. To extract the large-scale components of these signals $a_l$ and $b_l$, we lowpass filter each such that $a_l(z;t) = G * a(z;t)$, where $G$ is the impulse response function of a sharp spectral filter that is convolved with $a$. For the present study, we define the filter function to have a cutoff wavelength equal to half the height of the LLJ, $\lambda _c = h/4$, which generally is in the range separating the inner and outer peaks in the premultiplied spectrograms (figure 4). We further extract the small-scale component of each signal as $a_s(z;t) = a(z;t) - a_l(z;t)$.

The next step involves a Hilbert transform $\mathcal {H}$, which for the small-scale signal $a_s$, is defined as

(4.7)\begin{equation} \mathcal{A}_s(t) = \mathcal{H}\{ a_s(t)\} = \frac{1}{\rm \pi}\mathcal{P} \int_{-\infty}^{+\infty} \frac{a_s(\tau)}{t - \tau}\, {\rm d}\tau, \end{equation}

where $\mathcal {P}$ is the Cauchy principal value of the integral for time shift $\tau$. Mathematically, (4.7) is the convolution integral between $a_s(t)$ and the quantity $1/{\rm \pi} t$ such that $\mathcal {A}_s(t) = a_s(t) * (1/{\rm \pi} t)$. From the fundamental properties of the Hilbert transform (Mathis et al. 2009a; Bendat & Piersol 2010), $a_s(t)$ and $\mathcal {A}_s(t)$ form a complex analytic signal $Z(t)$ such that

(4.8)\begin{equation} Z(t) = a_s(t) + i\,\mathcal{A}_s(t) = A_s(t)\,\mathrm{e}^{{\rm i}\phi_s(t)}, \end{equation}

where $A_s(t)$ and $\phi _s(t)$ are the instantaneous modulus and phase of $Z(t)$ (Sreenivasan 1985; Tardu 2008; Mathis et al. 2009a). The modulus $A_s (t)$ of the analytic signal,

(4.9)\begin{equation} A_s(t) = \sqrt{a_s^2(t) + \mathcal{A}_s^2(t)}, \end{equation}

represents the envelope of the original signal, $E(a_s)$. Next, we lowpass filter the envelope of $a_s$ such that $E_l(a_s) = G * E(a_s)$, which is the final element required to determine the AM coefficients.

The AM coefficient in this example is given as the correlation coefficient between the large-scale component of $b$ and the large-scale component of the envelope of small-scale $a$, i.e.

(4.10)\begin{equation} R_{b_l,a_s}(z) = \frac{\langle b_l'(z;t) E_l'(a_s(z;t)) \rangle_{t}}{\sqrt{\langle b_l'^2(z;t) \rangle_{t}} \sqrt{\langle E_l'^2(a_s(z;t)) \rangle_{t}}}. \end{equation}

We note that (4.10) differs from that presented by Salesky & Anderson (2018) in generality since they considered both two- and one-point statistics, whereas in this present study we consider only one-point AM coefficients (i.e. at the same height). Their results indicate that the one-point statistics provide a stronger signal in terms of correlations when compared with the two-point AM coefficients. Since increasing stability further limits vertical turbulent transport, we expect two-point AM coefficients in the present study to be small.

For the decoupling procedure to be implemented appropriately, there needs to be adequate scale separation between the inner and outer peaks (Mathis et al. 2009a). Investigation of figure 4 indicates this condition is only met for cases A–C, with cases D and E ($C_r = 0.50,\ 1.00$ K h$^{-1}$, $h/L = 2.49, 4.17$, respectively) on the fringe. Herein, we present results using virtual tower output at 50 Hz frequency from cases A–E with the added caveat of marginal scale separation existing in cases D and E. In figure 10 we include the AM coefficients between large-scale $u_l$ and $w_l$ with small-scale $u_s$, $w_s$, $\theta _s$, $(uw)_s$ and $(\theta w)_s$. The AM coefficients are presented for each case as functions of wall-normal distance $z/h$ to identify the role of global stability on coupling between the large and small scales. As one can discern from the AM by large-scale $u_l$ (figure 10a–e), the largest correlations are found in the upper portion of the SBL, namely for $z/h > 0.4$. In this region, the values of $R$ are mostly negative and decrease moderately in magnitude with stability for a given height. Using the example of $R_{u_l,u_s}$ (figure 10a), a positive correlation near the surface can physically be interpreted as follows: the small-scale velocity $u_s$ increases due to modulation by a high-momentum LSM with $u_l > 0$, or decreases due to modulation by a low-momentum LSM with $u_l < 0$. Conversely, a negative correlation implies that on average, a high-momentum LSM will act to suppress small-scale perturbations, and a low-momentum LSM will excite small-scale perturbations. For the weakly stable case A, $R_{u_l,u_s}$ is negligible near the surface, decreases towards small negative values around $z/h \approx 0.2$ and further decreases to $R_{u_l,u_s} \approx -0.2$ in the upper half of the SBL. This behaviour is consistent with both the weakly convective case presented by Salesky & Anderson (2018) as well as the neutrally stratified case by Mathis et al. (2009b). This similarity also holds for all the other AM coefficients with modulation by $u_l$ except for $R_{u_l,\theta _s}$ (figure 10c), which is weakly positive throughout most of the SBL across all cases. In terms of overall magnitude for modulations by $u_l$, the largest impact is observed near the top of the SBL for $R_{u_l,(uw)_s}$, which reaches values as low as $-0.4$ in cases A and B (figure 10d). Negative coupling between large- and small-scale streamwise velocity at these heights is most likely associated with turbulence production by the LLJ.

Figure 10. The AM coefficients $R$ from cases A–E bin averaged versus $z/h$ for correlations with (a–e) $u_l$ and ( f–j) $w_l$. Small-scale envelopes include (a, f) $u_s$, (b,g) $w_s$, (c,h) $\theta _s$, (d,i) $(uw)_s$ and (e,j) $(\theta w)_s$.

By contrast, the AM coefficients for large-scale vertical velocity $w_l$ (figure 10f–j) are markedly smaller than those for $u_l$ across all cases. The only non-negligible coefficients occur under weak stability (cases A and B) for coupling between $w_l$ and the instantaneous second-order moments $(uw)_s$ and $(\theta w)_s$ (figure 10i,j) for the near-neutral case A. Even these values are modest, however, and again may likely be associated with turbulent transport below the LLJ. These weak $w_l$ AM coefficients may also stem from to the lack of an appreciable outer peak in the vertical velocity spectrograms of cases D and E (figure 4f–j).

There are a few core similarities and differences between the AM coefficients displayed in figure 10 versus those presented by Salesky & Anderson (2018) under varying convective stratifications. First, the coupling with small-scale instantaneous momentum flux $(uw)_s$ is the largest observed among all considered combinations of parameters in both the CBL and SBL. Moreover, the correlations with large-scale $w_l$ were relatively unaffected by global stability in both the CBL and SBL. However, in the SBL this is because the correlations were negligible across all simulations, whereas they were substantial in the CBL. With increasing stability, the results from figure 10 indicate that AM may occur due to large-scale $u_l$ but not necessarily for $w_l$, which is also consistent with the notion that buoyancy suppresses large-scale vertical motions (e.g. García-Villalba & del Álamo 2011).

In an attempt to better characterise the effects of local stability on AM, included in figure 11 are the AM coefficients plotted against $Ri_g$. Recall from figure 2( f) that $Ri_g$ nearly monotonically increases with height and stability for cases A–E throughout the SBL. We composited all five simulations A–E and bin averaged the AM coefficients based on evenly logarithmically spaced bins in $Ri_g$ with vertical error bars denoting the standard deviation of each bin in figure 11. It is apparent that on average, $u_l$ weakly correlates positively with small-scale parameters under weak local stability ($Ri_g < 0.05$, figure 11a–e), and these correlations decrease and eventually become negative under higher stability ($Ri_g > 0.1$). The largest spread in the correlations with $u_l$ occur around the critical Richardson number (Grachev et al. 2013), $Ri_g \approx 0.2$. All of these AM coefficients tend to level off with further increasing stability at relatively small values $-0.1 < R < 0.1$. The correlations with large-scale $w_l$ (figure 11f–j) are virtually zero for $Ri_g > 0.01$ across all parameters.

Figure 11. As in figure 10 but composited across all cases and plotted against $Ri_g$. Bin medians are plotted in blue and means in black with error bars denoting $\pm 1$ standard deviation.

The results following from the decoupling procedure discussed in this section are consistent with those throughout this study and in the literature: negative buoyancy in the SBL suppresses vertical motions at large scales, forcing coherent structures to become increasingly confined to horizontal planes and at increasingly local scales (recall figure 1a).

4.6. Conditional averaging

A common feature of LSMs in wall-bounded flows is their association with low- and high-speed streaks in the logarithmic layer (Adrian 2007), so it is therefore advantageous to composite snapshots of the flow when these conditions are present. This process is known as conditional sampling (Antonia 1981), which we can define for the streamwise velocity (in notation following Salesky & Anderson (2020), and with adapted conventions) as

(4.11)\begin{equation} \frac{\tilde{u}'^{{\dagger}}(\boldsymbol{x},t)}{u_{\tau 0}} = \left\langle \left.\frac{\tilde{u}'(\boldsymbol{x},t)}{u_{\tau 0}} \right| \frac{\tilde{u}'(\boldsymbol{x}_c,t)}{u_{\tau 0}} <{-}2 \frac{\sigma_{\tilde{u}'(\boldsymbol{x}_c)}}{u_{\tau 0}} \right\rangle_{N_{\alpha^-}}, \end{equation}

where $\tilde {u}'^{{\dagger}} (\boldsymbol {x},t)/u_{\tau 0}$ is the streamwise velocity averaged over $N_{\alpha ^-}$ instances where the flow is below the threshold $\alpha ^- = -2 \sigma _{\tilde {u}'(\boldsymbol {x}_c)}/u_{\tau 0}$ at the coordinate $\boldsymbol {x}_c = (x', y', z=0.05h)$ and $\sigma _{\tilde {u}'(\boldsymbol {x}_c)}$ is the standard deviation of velocity fluctuations. Here, $({\cdot })^{{\dagger}}$ is used to denote a conditionally averaged variable.

Conditionally averaged streamwise velocity, vertical velocity and potential temperature fields based on $\alpha ^-$ in (4.11) are included in figure 12 for cases A–E. The effects of stability are immediately apparent in all three averaged fields, with the extent of the conditionally averaged coherent structures diminishing in spatial extent (both horizontally and vertically) with increasing stability. In case A the streamwise extent of the central $\tilde {u}'^{{\dagger}} /u_{\tau 0}$ feature is $\approx 2h$, which corresponds well to the wavelength associated with the outer peak in the streamwise velocity spectrogram (figure 4a) of $\lambda _x / h \approx 1$–2. When conditioning on low-speed streaks near the surface, the $\tilde {u}'^{{\dagger}} /u_{\tau 0}$ minimum extends vertically up to $z/h \approx 0.3$ in case A, and appears to tilt upwards downstream from the central peak (figure 12a). This $\tilde {u}'^{{\dagger}} /u_{\tau 0}$ minimum diminishes under increasing stratification until the conditional low-speed streak becomes confined vertically and longitudinally (follow figure 12a,d,g,j,m sequentially). Combined with the evidence from the spectrograms, these results largely agree with the conceptual model of LSMs by, e.g. Marusic et al. (2010a). Interestingly, the corresponding field of $\tilde {w}'^{{\dagger}} /u_{\tau 0}$ in case A (figure 12b) features a maximum in vertical velocity directly overlaid with the low-speed streak, which highlights the dynamics of an ejection (recall § 4.4) in both momentum and temperature (figure 12c). These $\tilde {w}'^{{\dagger}} /u_{\tau 0}$ features include a similar downstream inclination as those in the $\tilde {u}'^{{\dagger}} /u_{\tau 0}$ minima, but also exhibit another vertically extending lobe upstream. The vertical and horizontal extents of these upstream and downstream lobes, respectively, decrease rapidly with increasing stratification. The correlation between $u$ and $\theta$ throughout the SBL is highlighted by how similarly the $\tilde {u}'^{{\dagger}} /u_{\tau 0}$ and $\tilde {\theta }'^{{\dagger}} /\theta _{\tau 0}$ fields evolve under increasing stability, as they are nearly identical qualitatively. In combination with the asymmetrical distribution of $\tilde {w}'^{{\dagger}} /u_{\tau 0}$, this again may be related to the presence of temperature microfronts (Sullivan et al. 2016) within the SBL that concentrate gradients in velocity and temperature.

Figure 12. Average fields conditioned on $\tilde {u}'/u_{\tau 0} < \alpha ^-$ as in (4.11) from simulations (a–c) A, (d–f) B, (g–i) C, (j–l) D and (m–o) E. Conditional fields include (a,d,g,j,m) $\tilde {u}'^{{\dagger}} /u_{\tau 0}$, (b,e,h,k,n) $\tilde {w}'^{{\dagger}} /u_{\tau 0}$ and (c, f,i,l,o) $\tilde {\theta }'^{{\dagger}} /\theta _{\tau 0}$. Here $\Delta x'$ represents the distance upstream and downstream from each instance meeting the conditioning criteria of a low-speed streak.

The conditionally averaged fields in figure 12 are a visual representation of the statistical results presented in §§ 4.2–4.5: buoyancy suppresses large-scale vertical circulations within SBL flows. Even under weak stability, the updrafts associated with ejections do not penetrate far above the surface and are roughly 70 % as wide as their corresponding low-speed streaks and cold air parcels. These fields are also reminiscent of the two-point correlations of $u$ and $w$ presented by Huang & Bou-Zeid (2013), which identified elongated coherence in the streamwise direction in the lower SBL. A similar analysis of our simulations (not shown) demonstrated the same pattern, which also decreased in extent with increasing stability.

5. Discussion and conclusions

Over the past half-century, investigations of turbulent wall-bounded flows have increasingly focused on the existence and dynamics of coherent structures (e.g. Kovasznay et al. 1970; Brown & Thomas 1977; Nakagawa & Nezu 1981; Murlis et al. 1982; Wark & Nagib 1991; Adrian et al. 2000; Ganapathisubramani et al. 2003; Tomkins & Adrian 2003; Del Álamo et al. 2004; Hutchins & Marusic 2007a; Marusic et al. 2010a; Salesky & Anderson 2018). A majority of the investigations into turbulent coherent structures have focused on flows under neutral and unstable stratification, but recent advances in computational resources and observational techniques have enabled further studies of stably stratified flows (e.g. García-Villalba & del Álamo 2011; Lan et al. 2018, 2019, 2022; Watanabe et al. 2018, 2019; Atoufi et al. 2021; Gibbs et al. 2022). This study builds upon previous research by simulating a suite of five SBLs using LES (Stoll et al. 2020) to examine the existence of turbulent coherent structures along with their role in governing SBL dynamics. We analyse these SBL simulations through a synergistic combination of mean profiles, instantaneous cross-sections, premultiplied spectrograms, linear coherence spectra, turbulent transport efficiencies, JPDFs, AM coefficients and conditionally averaged fields. Our key findings as related to the questions posed in § 1.3 are as follows.

1. (i) The outer peak in premultiplied spectrograms at weak stability diminishes with increasing stability until only an inner peak remains. This is notably different than what occurs under unstable stratification, for which the outer and inner peaks actually merge at intermediate wavelengths for increasing instability (e.g. Salesky & Anderson 2018).

2. (ii) For weak stability, the ratio between turbulent transport efficiencies of momentum and heat $\eta _{uw}/\eta _{\theta w}$ is nearly unity and is constant with height in the middle of the SBL, which is consistent with observed and simulated CBLs under weak instability (e.g. Li & Bou-Zeid 2011; Salesky et al. 2017). In the upper third of the SBL, $\eta _{uw}$ decrease towards zero at the height of the LLJ, whereas $\eta _{\theta w}$ decreases less rapidly and as a function of stability. The individual transport efficiencies are also appreciably smaller in the SBL than those reported in the CBL.

3. (iii) In § 4.1 we observe the existence of low- and high-speed streaks at weak stability, a telltale signature of canonical LSMs. These features decrease in coherence with increasing stability in conjunction with the attenuation of outer peaks in the spectrograms (§ 4.2). Analysis of linear coupling between flow parameters across scales in § 4.3 indicates that increasing stratification limits the vertical extent of coherent structures in the SBL. Without the added flux contributions by LSMs, vertical turbulent transport efficiencies decay for both momentum and heat with increasing stability and proximity to the LLJ (§ 4.4). By decomposing the simulated flows into large and small scales in § 4.5, we find that under increasing stability, horizontal motions remain correlated across scales whereas vertical motions are buoyantly suppressed throughout the SBL. Finally, by conditionally averaging on the presence of low-speed streaks near the surface in § 4.6, the resulting vertical cross-sections of $u$, $w$ and $\theta$ indicate the clear presence of LSMs under weak stability that are largely consistent with the conceptual models proposed in the literature (e.g. Marusic et al. 2010a; Baars et al. 2017; Salesky & Anderson 2018). Under increasing stratification, however, these coherent structures decrease in streamwise and vertical extent, and their intensities are attenuated.

Results from this study elucidate how vertical motions are unable to penetrate far beyond their initial levels, resulting in turbulence that is disproportionately horizontal (García-Villalba & del Álamo 2011; Huang & Bou-Zeid 2013). Further increasing stability acts to suppress turbulence in all directions, with characteristic motions becoming increasingly local in scale, weak in magnitude and decoupled from the surface (Grachev et al. 2013; Lan et al. 2018, 2019). These results share some similarities with those from free-shear flows under stable stratification (e.g. Watanabe & Nagata 2021; Akao, Watanabe & Nagata 2022), namely that inner–outer interactions are greatest in regions of strong wind shear. However, increasing stable stratification on free-shear layers has been found to increase the prominence of streamwise superstructures, whereas this was the opposite case found in our study. In wall-bounded flows, shear stress is maximized near the surface resulting in a relative minimum in the gradient Richardson number, thereby providing a more hospitable environment for the maintenance of turbulence.

While these conclusions are consistent with literature, further studies with high-resolution LES or DNS that better capture the spectrum of turbulent motions under moderate to high stability are certainly warranted (see discussion in Maronga & Li 2021). Additionally, it appears that LES of stably stratified wall-bounded flows can be rather sensitive to the wall model employed. Further studies should evaluate the performance of additional formulations of the Monin–Obukhov universal functions such as those of Grachev et al. (2007) that perform better under high stability, or with a gradient-based scaling framework (e.g. those evaluated by Sorbjan 2010; Greene et al. 2022).