## Abstract

The nighttime high-latitude stably stratified atmospheric boundary layer (SBL) is computationally simulated using high–Reynolds number large-eddy simulation on meshes varying from 200^{3} to 1024^{3} over 9 physical hours for surface cooling rates *C*_{r} = [0.25, 1] K h^{−1}. Continuous weakly stratified turbulence is maintained for this range of cooling, and the SBL splits into two regions depending on the location of the low-level jet (LLJ) and . Above the LLJ, turbulence is very weak and the gradient Richardson number is nearly constant: . Below the LLJ, small scales are dynamically important as the shear and buoyancy frequencies vary with mesh resolution. The heights of the SBL and Ri noticeably decrease as the mesh is varied from 200^{3} to 1024^{3}. Vertical profiles of the Ozmidov scale show its rapid decrease with increasing , with over a large fraction of the SBL for high cooling. Flow visualization identifies ubiquitous warm–cool temperature fronts populating the SBL. The fronts span a large vertical extent, tilt forward more so as the surface cooling increases, and propagate coherently. In a height–time reference frame, an instantaneous vertical profile of temperature appears intermittent, exhibiting a staircase pattern with increasing distance from the surface. Observations from CASES-99 also display these features. Conditional sampling based on linear stochastic estimation is used to identify coherent structures. Vortical structures are found upstream and downstream of a temperature front, similar to those in neutrally stratified boundary layers, and their dynamics are central to the front formation.

## 1. Introduction

Despite ongoing intensive study, understanding and parameterizing stably stratified turbulence in geophysical flows and, in particular, in atmospheric and oceanic planetary boundary layers (PBLs) remains at the forefront of geophysical turbulence research. The dynamics of the stable boundary layer (SBL) are still in a state of discovery as the couplings between large-scale forcings, internal waves, and weak turbulence over undulating rough boundaries in the presence of strong stratification continues to resist a generic description; see Fernando and Weil (2010) and Mahrt (2014) for recent reviews. The importance of the SBL is well appreciated in climate and weather modeling as these forecasts show great sensitivity to the form of their SBL parameterization (Holtslag et al. 2013). Intermittent SBL turbulence is also important in numerous application areas—for example, electromagnetic wave propagation (Wyngaard et al. 2001), interpreting observations collected from wind profilers (Muschinskii and Sullivan 2013), and air quality (Weil 2012).

A crucial thread in SBL research for climate and weather applications is quantifying the connecting relationship between the mean wind and temperature fields and the average turbulent fluxes and variances used in large-scale modeling applications; for example, Brost and Wyngaard (1978), Large et al. (1994), Svensson and Holtslag (2009), Sorbjan (2010), and Huang et al. (2013) all propose single-column models for the SBL. However, to further improve these flux–gradient relationships, to better understand mixing in the SBL, and to guide the interpretation of observational data collected in field campaigns requires a more fundamental understanding of the building blocks—that is, the coherent structures, in turbulent SBLs. It is now widely appreciated that coherent structures, loosely defined as spatially organized entities long lived in a Lagrangian frame of reference, are the key flux-carrying structures in most geophysical boundary layers. Observational and numerical studies find that the properties of the organized flow structures vary with stratification and the external driving forces in the PBL. For example, thermal plumes dominate the daytime convective PBL (Deardorff 1970), large-scale rolls are pervasive in mixed shear–convective PBLs (Moeng and Sullivan 1994; Fedorovich et al. 2004), low-speed streaks dominate near-neutral surface layers (Marusic et al. 2010), large eddies in canopy turbulence are determined by an inflection-point instability centered near the canopy top (Finnigan et al. 2009), and surface-wave-generated Langmuir circulations populate the upper ocean (Sullivan and McWilliams 2010), to mention a few.

Because of the challenges in observing and simulating stably stratified weak turbulence, the coherent structures in the SBL are less studied, and likely more variable, than those in the convective and near-neutral PBL. In the very stable boundary layer characterized by large Richardson number (Ri), Mahrt (2014) shows that the morphology of structures in the surface layer is exceedingly diverse with intermittent turbulence mixed with wavelike motions and two-dimensional modes. These surface-layer features differ from the highly intermittent Kelvin–Helmholtz instabilities that appear to dominate the overlying residual turbulence above the SBL top (Balsley et al. 2003). At lower Ri in the weakly stable regime, the turbulence is near continuous and a new collection of turbulent structures emerge in the SBL. For example, in a slightly stratified wind tunnel flow, Chen and Blackwelder (1978) analyze time series from a vertical rake of instruments and find that “the most interesting observation was the existence of a sharp internal temperature front … that extended throughout the entire boundary layer.” Their temperature front appears to be part of the family of cliff–ramp structures observed by Thorpe and Hall (1980) in a lake under moderate wind conditions and by Gao et al. (1989) above a forest under slightly unstable conditions. Chung and Matheou (2012) also find cliff–ramp fronts in direct numerical simulations (DNS) of stably stratified shear flow with no solid boundaries; their flow visualization also reveals that the fronts tilt farther forward (or downstream) and can become intermittent with increasing stratification. The dynamics behind the cliff–ramp fronts in boundary layers is not completely explained but is possibly linked to hairpin packets found in neutral wall-bounded flows as discussed by Adrian (2007). Williams and Smits (2011) speculate that hairpin packets become elongated in the downstream direction in a thermally stratified laboratory flow. Cliff–ramp structures are also ubiquitous features of passive scalars in the turbulent flows described by Warhaft (2000).

The specific goals of this article are to identify and characterize coherent structures, and in particular temperature fronts, in large-eddy simulations (LESs) of a rough-wall weakly stable PBL. Previous pioneering work by Mason and Derbyshire (1990), Derbyshire (1999), and Saiki et al. (2000) used LES with some success to investigate the SBL under weak stratification. However, in these simulations the combination of small-scale turbulence and coarse mesh, , places a heavy burden on the subgrid-scale (SGS) model that masks the numerous small-scale features we wish to identify. Thus, we build on past work but use much finer-grid meshes than in previous LES of the SBL (e.g., Beare et al. 2006; Huang and Bou-Zeid 2013). Our fine mesh LES with utilizing 10^{9} grid points follows recent trends in the use of high-resolution LES (Bou-Zeid 2015) and provides an opportunity to examine in a systematic manner the sensitivity of the solutions to the mesh spacing similar to our work with the convective PBL (Sullivan and Patton 2011). The problem posed is a canonical SBL with a homogeneous lower boundary. We are aware that coupling with nonhomogeneous or time-varying surface conditions (Van de Wiel et al. 2002; Nieuwstadt 2005; Flores and Riley 2011; Ansorge and Mellado 2014; Mironov and Sullivan 2016) can lead to a different family of structures in the SBL and possibly global turbulence collapse at least in low–Reynolds number DNS. Van de Wiel et al. (2012) and Donda et al. (2015) using a combination of analytic models and DNS provide evidence that turbulence collapse is likely a transient state even with high surface cooling. A discussion of the turbulent-to-laminar transitional regime at large Ri is beyond the scope of the present article. The road map of the manuscript is as follows: a brief description of the LES equations and numerical algorithm is given in section 2, an outline of the numerical experiments is provided in section 3, results from the grid resolution tests and experiments with surface cooling variations are given in section 4, the identification of temperature fronts and coherent structures are presented in sections 5 and 6, and a summary of the findings is provided in section 7.

## 2. LES equations

The LES model equations for a dry atmospheric PBL under the Boussinesq approximation with system rotation and stable stratification with a flat bottom boundary are well documented (e.g., Moeng 1984; McWilliams et al. 1999; Moeng and Sullivan 2015):

The above equation set includes transport equations: (1a) for momentum ; (1b) for virtual potential temperature ; and (1c) for SGS turbulent kinetic energy *e*. In (1d) the divergence-free (incompressible) condition determines the elliptic pressure variable . Variables that appear in (1) are velocity components , geostrophic winds , rotation vector with Coriolis parameter *f*, unit vector in the vertical direction, and buoyancy parameter where *g* is gravity and is the still-air potential temperature. In the later discussion, we also make reference to the pressure and air density *ρ*, which do not appear explicitly in (1). The overbar notation denotes a spatially filtered quantity.

The LES equations are formally derived by applying a spatial filter term by term to the governing equations of motion. This operation introduces unknown SGS kinematic momentum and temperature fluxes (e.g., Sullivan et al. 2003):

An abundant number of prescriptions are possible for the parameterization of these SGS fluxes, [e.g., Kosović (1997), Bou-Zeid et al. (2005), Bhushan and Warsi (2005), Lévêque et al. (2007), and Ramachandran and Wyngaard (2011) to mention just a few]. For simplicity and because of the flow regime of interest, we adopt the two-part SGS model proposed by Sullivan et al. (1994) that utilizes the transport equation in (1c) and an eddy viscosity approach in the parameterization of the SGS fluxes given by (2). This parameterization is specifically tailored to high–Reynolds number LES that uses rough-wall surface boundary conditions based on Monin–Obukhov (MO) similarity theory. Past experience shows this model noticeably improves the agreement with empirical MO similarity functions for neutral flows in the surface layer. However, all SGS models invoke assumptions and thus we primarily rely on fine mesh resolution (see section 3) to provide adequate separation between the large resolved anisotropic and parameterized small-scale isotropic eddies and, thereby, minimize the reliance on the exact form of the SGS model. The right-hand side of (1c) for the SGS energy includes a standard suite of terms. In symbolic notation is the bidirectional energy transfer from resolved to SGS motions, is SGS buoyancy production/destruction, is a diffusion term, and is viscous dissipation. The specific formulas for these terms used in our LES implementation are not repeated here but are available in numerous references (Deardorff 1972; Moeng 1984; Moeng and Wyngaard 1989; Sullivan et al. 1994; McWilliams et al. 1999; Moeng and Sullivan 2015).

In our LES, the sidewall boundary conditions are periodic and a radiation boundary condition (Klemp and Durran 1983) is used at the top of the domain. As is common practice with geophysical flows, we impose rough-wall boundary conditions based on a drag rule where the surface transfer coefficients are determined from MO similarity functions (Moeng 1984). In the present application, the MO rules are applied point by point at the lower boundary as described by Mironov and Sullivan (2016). As a first approximation, in the limit of extremely fine computational meshes compared to the physical roughness the undulations can be resolved with point-by-point drag rules applied as a model of the even smaller unresolved roughness (Sullivan et al. 2014). The use of local surface exchange coefficients is an approximation but is supported by the analysis of Wyngaard et al. (1998), who performed an in-depth study of surface flux conservation equations. In the limit of large-grid aspect ratio , where is the first model level, they find the conservation equations are in a state of local equilibrium justifying the usual practice of applying local MO rules. For smaller aspect ratios, the conservation equations become stochastic but their coupling with the LES equations yields only marginal improvements in the prediction of the nondimensional mean shear and mean temperature gradient profiles. Wyngaard et al. (1998) conclude the form of the SGS model has a much greater impact on the overall LES predictive capabilities in the surface layer than the surface flux exchange rule.

We utilize well-established algorithms to integrate the LES equations in (1). The equations are advanced in time using an explicit fractional step method that enforces incompressibility at every stage of the third-order Runge–Kutta scheme. Dynamic time stepping with a fixed Courant–Fredrichs–Lewy (CFL) number is employed, which we have found naturally adapts to a wide spectrum of dynamical processes. The spatial discretization is second-order finite difference in the vertical direction and pseudospectral in horizontal planes. For the present application, the advective terms in the momentum equations are written in rotational form, while a flux-conserving form is used for the advective terms in the scalar (temperature) equation. The vertical velocity equation is solved for the deviation of *w* from its horizontal mean value at each height. The flow variables are explicitly filtered at each time step, or dealiased, using the 2/3 rule (Moeng and Wyngaard 1988). Further algorithmic details are given by Moeng (1984), Sullivan et al. (1994, 1996), McWilliams et al. (1999), Sullivan and Patton (2011), Moeng and Sullivan (2015), and the references cited therein.

To streamline the notation and text in the following discussion, we now drop the overbar symbol on all spatially filtered resolved variables and simply refer to virtual potential temperature *θ* as “temperature.”

## 3. Design of LES experiments

The SBL flow examined here is the first GEWEX Atmospheric Boundary Layer Study (GABLS1) described by Beare et al. (2006). This high-latitude SBL is a benchmark intercomparison case for canonical stable LES and, in addition, serves as a clean test case for the evaluation of single-column PBL schemes used in climate and weather models (Cuxart et al. 2006). GABLS1 does not include potentially important feedbacks from a land surface as investigated by Holtslag et al. (2007). The GABLS1 PBL is driven by steady geostrophic winds and a time-varying surface temperature; the surface is cooled at a constant rate starting from an initial surface temperature of 265 K. Other important input parameters are the Coriolis parameter , surface roughness for momentum and temperature, buoyancy parameter , and still-air potential temperature . Brost and Wyngaard (1978) and later Kosović and Curry (2000) found that the recipe of gradually cooling the surface at a fixed rate is advantageous as it does not introduce artificial long-lived transients but allows continuous stably stratified turbulence to develop in the SBL. For the GABLS1 set of inputs, multiple LES codes, including the present one, with different SGS models and numerics predict a well-developed near-equilibrium PBL is reached after 9 physical hours with stratified rotated winds featuring a low-level supergeostrophic wind maximum of about 10 m s^{−1} below the boundary layer top m (Beare et al. 2006); also see Fig. 1 from Huang and Bou-Zeid (2013).

We adopt the GABLS1 computational domain and impose the two-layer structure on the initial temperature sounding

where is the initial inversion height. The initial and final temperature profiles are compared in the right panel of Fig. 1. In terms of , the computational domain is , which is sufficient to allow turbulent flow fields to fully develop independent of the periodic sidewall boundary conditions. However, in order to assess possible impacts of the finite domain size on the solution results, and organized structures, a simulation Bw is also performed in a twice-as-wide domain . All simulations are started from small random perturbations in temperature near the surface with the initial profile of SGS energy for .

The present suite of LES experiments extends the original GABLS1 problem design in the following ways. First, four different cooling rates are applied at the surface^{1} that increase the overall level of stratification in the SBL compared to GABLS1 and, second, three levels of grid refinement are considered, namely grid spacing , which corresponds to meshes with *M*^{3} = grid points. For each mesh, the grid spacing is constant (isotropic) across the three coordinate directions ; the wide-domain simulation Bw with mesh also has spacing of 0.78 m. It should be noted that previous LES investigations of the GABLS1 PBL by Beare et al. (2006) and Beare and Macvean (2004) used finest resolution meshes with spacing , respectively, Huang and Bou-Zeid (2013) use a fine mesh , and Matheou and Chung (2014) also use . Thus the finer spatial resolution and the grid mesh variation of more than 5 times in each coordinate direction used in the present simulations is a clear advance over past work.

The simulations are integrated for more than 9 physical hours, which equates to more than 40 nondimensional times for the simulation with the highest cooling rate. Significant computational resources are needed to simulate 9 physical hours on a mesh with grid points. Most often, the time step is limited by the CFL constraint based on the maximum horizontal wind , where . Thus, to increase the allowable , all simulations are performed on a horizontally translating mesh moving at a constant speed of (Sullivan et al. 1996). Further, the time step is determined adaptively based on a fixed CFL number equal to 0.5. Even with these enhancements, approximately 900 000 time steps are required on the finest mesh with an average time step over the last hour of the simulation. To complete each fine-mesh simulation running on 2048 processors requires ~2.5 × 10^{6} core hours or ~1130 wall clock hours on the NCAR peta-scale machine Yellowstone. The code parallelization and performance is described by Sullivan and Patton (2011).

The statistics presented here are generated in flight and also by further postprocessing of archived data volumes. As is customary practice, statistics are computed by averaging in horizontal *x*–*y* planes at each *z* and also over time; this is the LES approximation to an ensemble (mean) average. The averaging time window is from hours 8 to 9 as in the original GABLS1 LES comparisons. Averages are denoted by angle brackets with resolved turbulent fluctuations indicated by primes; for example, for variable *f* its mean is and its turbulent fluctuations are simply . Variables marked with subscript ⊥ lie in the horizontal *x*–*y* plane: . The average surface friction velocity and surface kinematic temperature flux are denoted *u*_{*} and *Q*_{*}, respectively, and the Monin–Obukhov length , where is the von Kármán constant. Table 1 provides a summary of the simulation details and bulk boundary layer variables.

## 4. Interpretation of low-order moments

In addition to the bulk quantities in Table 1, our statistical analysis of the SBL includes computation of typical vertical profiles of low-order moments, namely means, variances, and momentum and temperature fluxes. In the discussion of these statistics, we refer to the height of the low-level jet (LLJ) or wind maximum defined as the vertical location where reaches a maximum: . Then the instantaneous (mean plus fluctuation) horizontal wind aligned with the mean wind direction is given by the vector dot product . Also, we introduce the gradient Richardson number (Ri) defined in terms of the buoyancy frequency *N* and shear frequency *S*:

The boundary layer top is found using a modification of the local gradient method described in Sullivan et al. (1998). With this method, is defined as the average of the vertical locations where is a maximum along every vertical column. This is usually a robust recipe and because the horizontal averaging operates on grid points it yields a smooth variation of with *t*. However, in the present application with very fine meshes we find the local gradient method often finds a false low estimate of . The source of this bias is traced to multiple points in the middle to lower PBL with very large vertical temperature gradients, exceeding the mean gradient of the overlying temperature inversion, which spoils the area averaging used to estimate . The dynamical reasons for this are extensively discussed in section 5a. To avoid these false triggers in the stratified PBL, we instead search for the location of the maximum average temperature gradient—that is, here is defined as the location where is a maximum. This modification produces acceptable and consistent results for estimating the boundary layer top independent of the grid resolution. Beare et al. (2006) and Huang and Bou-Zeid (2013) use an alternate definition of based on extrapolation of the average vertical momentum flux, which gives lower estimates of compared to the gradient method, especially so for strong surface cooling. This minimum flux method estimates the boundary layer top near . We prefer the gradient method as our simulations show boundary layer turbulence is supported between and we view the LLJ as part of the SBL.

### a. Varying resolution with fixed cooling rate

Profiles of the mean wind speed, mean temperature, gradient Richardson number, and shear and buoyancy frequencies are displayed in Figs. 1 and 2 for varying mesh resolution at a fixed cooling rate . First, our wind and temperature profiles from simulation A with are indistinguishable from those we supplied to the GABLS1 intercomparison and, second, are in good agreement with results from other LES codes with similar resolution [see Figs. 2 and 3 from Beare et al. (2006)]. Our new simulations B and C with mesh spacing , which are identical to simulation A in all other respects, highlight the sensitivity to a change in mesh resolution. We observe a systematic variation where a coarse mesh produces a deeper SBL with an overall higher level of stratification. Also, *u*_{*} and *Q*_{*} are, respectively, 13% and 32% larger on the coarse mesh compared to the finest mesh (see Table 1) and in the SBL interior Ri is largest on the coarse grid (see Fig. 2). Further reexamination of the GABLS1 intercomparisons also shows the same overall trend for boundary layer depth as the mesh decreases from 12.5 to 2 m for all LES codes (Beare et al. 2006, p. 253). A broadly similar effect is also found in LES of a convective PBL where tended to be higher and the entrainment rate increased on coarser meshes (Sullivan and Patton 2011). On a fine mesh , the SBL is shallower, decreases, the magnitude of the LLJ slightly increases, and the wind turning with height is sharper. Figure 2 shows that the shear and buoyancy frequencies (or velocity and temperature gradients) and, thus, Ri are sensitive metrics for judging whether the LES solutions are mesh independent in this stably stratified PBL. The shear frequency in particular exhibits the greatest sensitivity to the mesh resolution as it tends to increase with decreasing , thus leading to a lower Ri; also increases but appears to be slightly less sensitive to . Subtle changes in near and above the LLJ lead to strong stratification with . Meanwhile, the LES carried out with a fine grid supports weak stratified turbulence as Ri is maintained at a constant level near the critical value of ~0.25 in the region . Finally, a comparison of the wind, temperature, and Richardson number profiles from simulations B and Bw are nearly identical; this demonstrates that the smaller computational domain 400^{3} m is adequate for capturing the largest scales of motion in our SBL and that the solution mesh sensitivity is due entirely to small scales. DNS of stratified homogeneous shear flow show a dependence on the computational domain size (Chung and Matheou 2012) and large scales are also reported in DNS of a stably stratified Ekman layer with zero buoyancy gradient aloft by Ansorge and Mellado (2014). We believe that their results are a consequence of the problem posing; these DNS do not contain a stably stratified capping inversion or a LLJ. In our SBL, and impose physical constraints on the biggest scales that can develop in the boundary layer. Based on the statistics, the LES solutions are still changing at but appear to nearly converge as is reduced from 0.79 to 0.39 m.

Our mesh resolution tests indicate that the small scales and smaller are dynamically active in LES of the SBL especially in the region near the LLJ, and it is likely that the SGS model viscosity tends to overly damp motions when . McWilliams (2004) shows that decreasing SGS viscosity, or effectively boosting the large-eddy Reynolds number, increases the number of overturning events at small scales, thus lowering the value of Ri. The analysis of stratified residual-layer turbulence by Balsley et al. (2008) uses 1D spatial filtering to show that Ri is scale dependent; small scales tend to have lower values of Ri compared to the average background stratification. This speculation is further clarified by considering the Doughtery–Ozmidov length scale where *ϵ* is the viscous dissipation; see Doughtery (1961) and Ozmidov (1965). Recall eddies with vertical scale are preferentially influenced by buoyancy and do not overturn while small eddies are conceptually free of stratification influences. In the SBL, vary with distance from the wall and thus it is illuminating to consider a vertical profile of . Figure 3 (left panel) shows how varies with *z* depending on the mesh resolution; the viscous dissipation is calculated from the SGS model with but varying with stratification (Moeng and Wyngaard 1988). We observe that over the middle of the SBL is indeed larger than typical LES mesh spacings, but near the surface and especially in the region above the LLJ, decreases markedly and suggested as adequate by GABLS1 and Huang and Bou-Zeid (2013). Typical SGS models are dissipative and not designed to operate in a regime with . Apparently, in the entrainment zone LES with coarse resolution is incapable of resolving important small-scale turbulent overturning events. The consequence is a reduced value of and a false sharp increase in Ri. Our mesh resolution tests show that and Ri both decrease with finer resolution, but for different causes; increases with decreasing , which lowers , but better resolution of small-scale turbulence enhances to an even greater extent with the combined effect of lowering Ri, as shown in Fig. 2. Underresolution of the top and bottom of the SBL further appears to couple with the interior flow generating increased values of Ri. Jonker et al. (2013) shows that adequately resolving wall layers and entrainment zones using DNS leads to solutions that are nearly independent of Re. The above remarks are specific to our pseudospectral code using the two-part SGS model described previously; further work is needed to establish grid resolution requirements for simulations of stratified turbulence using different SGS models (e.g., Khani and Waite 2014).

### b. Impacts of cooling rate with fixed resolution

The impacts of the surface cooling rate on the boundary layer statistics are next discussed for LES solutions obtained with a fine mesh ; coarser meshes are not considered because of the mesh sensitivity found in section 4a. The cooling rate in simulation F is 4 times larger than the value used in the original GABLS1 design, and at the temperature difference between the top of the SBL and the surface results in strong stratification. Overall, for our suite of simulations, the bulk stratification increases by more than a factor of 3 as varies from 1.7 to 6.0. Figure 3 adds a cautionary note on the use of coarse meshes in LES of the SBL. The maximum value of decreases by more than a factor of 6 over the range of stratification considered, and above can be considerably less than the Δ = 2-m mesh adopted as adequate resolution for the GABLS1. It is noteworthy that even when , increasing stratification reduces the TKE on a fixed grid and thus places an even greater reliance on the SGS model.

Table 1 and the sequence of Figs. 4–8 summarize the impacts of surface cooling on the bulk statistics and the vertical profiles of typical low-order moments in the SBL. First, we notice that with increasing surface cooling (or stratification) and both decrease, but increases, and then the region above the LLJ occupies a bigger fraction of the SBL. Thus, to encompass regions above and below the LLJ, most often we use the nondimensional vertical coordinate in presenting the vertical profiles. Huang and Bou-Zeid (2013) and others define a boundary layer top *h* based on a stress minimum, which tends to occur near , and plot their results over the restricted range . Thus, their results exclude the region above the LLJ.

Inspection of the wind profiles shows that, with increasing stratification, the SBL is shallower, the height of the LLJ descends, the winds turn more sharply with height, and the surface wind stress decreases; see Fig. 4. At the same time, the mean temperature profile develops sharper vertical gradients in the lower boundary layer and weaker gradients aloft—especially so for the simulation with the highest cooling rate. The SBL appears to divide into two regimes with increasing stratification, which is particularly apparent in the *θ* profile for simulation F and in the mean vertical gradients of wind and temperature shown in Fig. 5. Above the LLJ , and both decrease but at rates sufficient to maintain a constant Ri near or slightly below the critical value of ~0.25. Simulation F appears to be entering a more strongly stable regime as and decrease by more than a factor of 10 compared to their values at . The sharp increase and decrease in Ri in simulation F near occurs because of delicate transitions in the wind and temperature gradients with height; Fig. 5 shows that first starts to noticeably decrease at a vertical location where is near constant or slightly increasing inducing sharp changes in Ri. In all simulations, the minimum value of tends to occur where and . For all cooling rates, resolved stably stratified turbulence with small vertical momentum fluxes and low levels of turbulence kinetic energy (TKE) are maintained by the fine grid resolution above the LLJ; see Figs. 6 and 7. The inset of Fig. 7 indicates that the TKE above is less than 10% of the surface stress magnitude for simulation C and less than 2% for simulation F. However, in all simulations the resolved turbulence dominates, outside of a very thin layer near the surface, as the ratio even in the region above the LLJ as shown in Fig. 7.

Below the LLJ, , and above the surface layer, , the bulk shear is nearly constant with height but increases in magnitude with stratification. At the same time, the buoyancy frequency tends to increase with *z* and with stratification; and change at rates such that the gradient Ri increases with *z* and stratification but is always less than 0.25. For example, at , Ri smoothly increases from 0.12 to 0.22 as varies from 1.7 to 6.0. These changes are consistent with decreasing turbulence levels resulting from the increasing levels of stratification. The vertical momentum fluxes in Fig. 6 support the observations about the structure of the SBL mentioned earlier. There is a consistent trend for the streamwise momentum flux to decrease and spanwise momentum flux to increase with higher cooling rates inducing shallower boundary layers. A measure of the resolved nature of the flow fields is provided in the TKE plot shown in Fig. 7. For all cooling rates, the SGS energy computed from (1c) near the surface is less than 20% and at is less than 10% of the total TKE.

Profiles of the vertical and horizontal temperature fluxes for varying cooling rates are compared in the left and right panels of Fig. 8, respectively. These fluxes are normalized by the surface flux to preserve the sign and ease the interpretation of the fluxes. In the mid- to lower SBL, the vertical temperature fluxes are near-linear functions of as expected. In the upper region, the mean flux profile displays more curvature and above the vertical flux nearly collapses because of the increasing stratification; as shown in Fig. 5. We find that the horizontal scalar fluxes and are comparable in magnitude to the vertical scalar flux throughout the bulk of the SBL. Large net horizontal fluxes persist despite the horizontal homogeneity of the flow fields and surface boundary conditions; that is, the mean gradients . This is explained by the budgets for total horizontal scalar flux. For example, Wyngaard et al. (1971) shows assuming horizontal homogeneity and using the classic Reynolds decomposition that

For the purposes of discussion, we have changed our notation in (5) and use angle brackets to denote Reynolds averaging with total fluctuations indicated by primes. The right-hand side of the above budgets contains the usual suite of terms; in order, they are mean shear and mean stratification production/destruction, turbulent transport, and pressure gradient interaction (e.g., Hanjalić 2002). Scalar flux budgets derived from the LES equations are identical to (5) when the SGS contributions are negligible (e.g., Mironov and Sullivan 2016). In (5) it is important to notice that the vertical scalar flux is coupled to horizontal fluxes through mean shear tilting and . In the lower SBL, , the mean velocity gradients and are greater than zero and hence mean shear tilting and also mean stratification production are significant sources of positive scalar flux. We note that SGS eddy viscosity models, including the dynamic approach, do not contain these tilting terms (Hatlee and Wyngaard 2007). Figure 8 shows that the horizontal fluxes are of similar magnitude as the vertical scalar flux and near the surface are 2 times larger; also see section 6a. Our LES results for vertical and horizontal scalar fluxes are in good agreement with surface-layer observations obtained under stable conditions (Wyngaard et al. 1971, their Fig. 4). In stable boundary layers with heterogeneous surface boundary conditions, one expects the horizontal temperature fluxes to play an even larger role in the budgets of temperature variance and flux.

## 5. Temperature fronts

### a. Flow visualization

Extensive visualization of the LES flow fields is carried out using selected 2D planes (*x*–*y*, *x*–*z*, *y*–*z*) and 3D volumes for all mesh resolutions. From this large database we concentrate on the temperature field *θ* from simulations using a mesh of 1024^{3} grid points. One of the most ubiquitous features common to all simulations is displayed in Figs. 9 and 10. The top panel of Fig. 9 shows a grayscale image of the temperature difference *θ* minus *θ*_{o} in an *x*–*z* plane late in the simulation *t* ≈ 9 h from simulation C; the top panel of Fig. 10 shows the same field but displayed as a family of contour lines. The bottom panel of Fig. 10 displays results for higher stratification from simulation F. Inspection of these figures reveals an impressive array of locally compacted temperature lines; that is, sharp temperature fronts tilted in the streamwise (or downstream) direction. Each front marks a sharp boundary between warm upstream and cool downstream air with more uniformly well-mixed air between fronts. In the lower boundary layer, the fronts tilt upward often with an apparent origin near the surface. Near the low-level jet, the fronts are weaker with values of tilt angle that transition from positive, to zero, to negative (downward tilt) as *z* moves across the LLJ. Animations show that the images in Fig. 9 are not isolated special instances and that the fronts meander and evolve in time but maintain clear coherence as they propagate with the winds. The temperature fronts are mainly tilted into the mean wind direction with finite spatial extent in the crosswind direction—that is, the fronts are not 2D sheets as shown in the bottom panel of Fig. 9. Increasing stratification reduces their vertical tilt—the fronts are tipped more toward the downstream direction—and also narrows the spacing between individual fronts; see the bottom panel of Fig. 10. We emphasize that in the current problem posing the surface boundary condition for temperature varies in time but is spatially homogeneous, constant across the *x*–*y* domain, and thus the observed fronts are internally generated by the dynamical interaction between well-resolved turbulence and a stably stratified temperature field.

Closer inspection of Figs. 9 (top panel) and 10 shows that the front tilt varies depending on the *x*–*z* location in the domain. For example, consider the two fronts marked by dashed white lines in the top panel of Fig. 9. The front that starts and ends at points (35, 50) and (155, 93) m, respectively, has an average tilt angle of , while the front that starts and ends at (290, 20) and (320, 50) m, respectively, has a steeper tilt with angle ; is measured positive counterclockwise from the *x* axis. It is illuminating to estimate the tilt of a constant-*θ* surface in the presence of turbulence. The gradient vector points in the direction normal to the surface and is utilized in vector cross products and , where and are unit vectors, to define vectors lying in the plane of a temperature front. Then a constant-*θ* surface has instantaneous tilt angles in the directions:

In (6), we invoke and split the vertical temperature gradient into its mean and fluctuation to expose the turbulence effect. Equation (6) is interesting: it has limits for low-amplitude temperature fluctuations acting on a strong stably stratified mean vertical temperature gradient. At the other extreme, for high-amplitude fluctuations superimposed on a weak mean vertical gradient, depend only on the ratio of horizontal to vertical gradients of . Recall and thus for a forward-tilting warm–cool front. In section 6a, based on conditional sampling, we find the horizontal gradient of the *θ* fluctuations aligned with the mean wind is comparable to or larger than for our weakly stratified SBL. Thus if the turbulent temperature fluctuations are large compared to the background vertical gradient of *θ*. Evidently the fronts at the top of the boundary layer in Fig. 10 are nearly flat, level surfaces because of weak horizontal temperature fluctuations while those in the lower boundary layer are tilted upward when the temperature fluctuations are vigorous with strong local horizontal and vertical gradients.

A horizontal *x*–*y* slice of taken at is provided in Fig. 11. In this image we spot multiple fronts sprinkled throughout the horizontal domain. Again, the characteristic signature is warm (cool) air upstream (downstream) of the frontal boundary. Recall that the spacing and thus the warm–cool regions are highly resolved features in our LES. Based on the contours in Fig. 11, a bulk estimate of the temperature jump is typically more than 0.2°. The four-panel image in Fig. 12 illustrates the *x*–*y* spatial coherence of a front with increasing distance from the surface. The frontal boundary, although irregular in *x*–*y*, remains sharp and clearly identifiable. It appears to rotate slightly in a clockwise direction and is found at a larger downstream *x* distance with increasing *z*. This apparent propagation with increasing *z* reflects the downstream tilt observed in Fig. 9, while the front’s rotation appears to track the rotation of the mean wind vector with height. The positive correlation between the wind and temperature fluctuations upstream and downstream of a frontal boundary in *x*–*y* planes leads to positive horizontal temperature fluxes and in the lower boundary layer as shown in the vertical profiles of Fig. 8.

Figure 13 displays a 2D image of the fluctuating vertical temperature gradient normalized by the local mean temperature gradient computed from the *θ* field shown in the top panel of Fig. 9. The most intense positive values of this statistic nicely overlap with the front locations found previously, and the compression of the vertical gradients into thin streaky filaments provides a sense of the front sharpness. Notice that the largest fluctuating vertical gradients are often more than 10 times larger than the local mean gradient. Overall the image shows that at a fixed *x*–*y* location the temperature progression from the surface (cool) to the SBL top (warm) occurs in a series of jumps or in staircase fashion, and between jumps the temperature is relatively well mixed despite the overall bulk stable stratification; see Fig. 4. The probability density functions (PDFs) shown in Fig. 14 quantify the skewed character of the vertical and horizontal temperature gradients. Near the surface , the vertical gradient is positively skewed while at the same location the horizontal temperature gradient aligned with the mean wind is negatively skewed consistent with a forward-tilting front. Above the LLJ, the PDF of the horizontal gradient appears to be slightly skewed toward positive values in agreement with the flow visualization. We note that the thickness of a typical front is multiple grid points wide but the thickness tends to vary with the LES grid resolution—that is, it depends on the SGS eddy viscosity . In our simulations, decreases by more than a factor of 8.5 and the standard deviation of the spatial temperature gradients increases by more than a factor of 5 as the mesh varies from Δ = 2 to 0.39 m. A sensitivity of the frontal width to viscosity is supported to a limited extent by the measurements of decaying turbulence by Tong and Warhaft (1994), who find that their passive scalar fronts become sharper and more intense as the molecular Reynolds number increases (Warhaft 2000, p. 215).

We are aware that warm–cool temperature fronts are also found in idealized low-Re direct numerical simulations of homogeneous shear flows (no walls); for example, see Gerz et al. (1994), Warhaft (2000), and Chung and Matheou (2012). The warm–cool fronts found here are cousins to those in homogeneous shear flows but differ because of the presence of a rough wall, the vertically varying stratification , and the presence of a stably stratified capping inversion. In an SBL, a noteworthy difference is the variability in the frontal tilt angle with increasing *z* especially in the region of the LLJ.

### b. Instantaneous vertical profiles of θ

We attempt to use our LES results to help interpret outdoor observations by simply asking the following: What is the signature of a coherent warm–cool temperature front in a fixed-tower-based measurement? To expand on this idea, we position a “virtual tower” in the LES domain at a fixed horizontal location and monitor the instantaneous temperature profile as function of time. Typical results are shown in Fig. 15 for a virtual tower placed at in Fig. 9 over a time period of 25 s. Assuming a horizontal advection speed of 5 m s^{−1}, the system of fronts in Fig. 9 propagates more than 100 m over this time period. In this example, three temperature fronts are found in profile 1, near matching the temperature contours in Fig. 10. As time advances, we observe that the jump in *θ* at these locations stiffens, relaxes, and stiffens, but generally a vertically well-mixed region is maintained between the fronts. Note that the horizontal advection speed increases with height and thus the fronts appear at different times in our tower measurement. Because of the front’s forward tilt, the location of the sharp vertical gradient in *θ* descends with increasing time when observed at a fixed location. We emphasize that this apparent descent is not a consequence of downward vertical advection as we are sampling the same frontal structure as time advances but at different locations along its tilted boundary because of streamwise advection. Mahrt (2014) discusses the very stable boundary layer where downward advection of turbulence from the LLJ occurs above the surface layer. Inspection of numerous profiles shows that the process sketched in Fig. 15 is generic with the vertical temperature profile often taking on a staircase shape with increasing *z*.

The LES findings shown in Fig. 15 are tantalizing targets for observations. Thus, we next search for SBL temperature fronts in the observational database collected during the 1999 Cooperative Atmosphere–Surface Exchange Study (CASES-99) field campaign [for an overview, see Poulos et al. (2002)]. Typical temperature profiles collected from the 55-m-tall tower are displayed in Fig. 16 under conditions broadly similar to the LES. At 0730:00 UTC 21 October 1999, the observations reported (*u*_{*}, *Q*_{*}) ~ (0.16 m s^{−1}, −0.025 m s^{−1} K), a LLJ with winds 12–13 m s^{−1} at a height *z*_{j} ~ 100–150 m (Balsley et al. 2003), and winds from 5 to 6 m s^{−1} at *z* = 45 m. Although qualitative, inspection of the profiles clearly shows the main features predicted by the LES; that is, multiple vertically well-mixed regions bounded by sharp temperature gradients aloft. Also, the maximum temperature gradient shows a steady descent with increasing time—a key feature of a streamwise advecting tilted temperature front shown in Fig. 15. Note the temperature jump across the front can approach 2 K because of stronger surface cooling during the observation period. However, because of the limited tower height, we are only able to observe fronts that propagate into and across the surface layer; Fig. 15 shows multiple fronts distributed vertically for *z* > 60 m. Finescale measurements collected from a slowly ascending kite during CASES-99 (Sorbjan and Balsely 2008; see their Fig. 1a) also support the LES predictions of intermittent temperature jumps in the middle and upper regions of the SBL.

## 6. Linear stochastic estimation

The flow visualization described in section 5 paints an image of many randomly distributed but locally organized warm–cool temperature fronts populating the weakly stable boundary layer. These sharp fronts are tilted in the downstream direction, exhibit spatial spanwise and vertical coherence, and propagate in time as organized entities. We seek to identify in an average sense the flow structures and dynamics that couple with the formation of these many warm–cool temperature fronts.

To establish a connection between flow structures and temperature fronts in our simulations, we use event-based conditional sampling—specifically, linear stochastic estimation (LSE); for an excellent essay on the topic of stochastic estimation, see Adrian (1996). LSE has several advantages over other eduction techniques. Conditional averages are obtained using unconditional variances, covariances, and two-point spatial correlations, and thus LSE conditional estimates are robust since all available data are used in forming the averages. Also, with LSE it is straightforward to use complicated detection events at multiple locations in the estimation procedure. As LSE conditional estimates depend linearly on the event data once the underlying statistical functions are computed and stored, parametric variations in the event data are easily formed. Finally, numerous tests using experimental and simulation data for a variety of turbulent flows show that LSE is an excellent approximation of a conditional average; for example, see Adrian et al. (1987), Guezennec (1989), Adrian et al. (1989), Adrian (1996), Christensen and Adrian (2001), and Richter and Sullivan (2014).

In the terminology of turbulent structure identification, we are interested in the average state of the SBL flow fields subject to a particular set of events; that is, we wish to compute the conditional average where the events are chosen specific to our application; are the velocity, pressure, and temperature fields; and angle brackets denote averaging over many realizations. Linear stochastic estimation is a technique for approximating the conditional average .

To make the ideas concrete, a short description of the technical details of our LSE implementation as applied to the SBL are outlined below. Based on the flow visualization of the temperature and velocity fields (Fig. 11) and the vertical profiles of the scalar flux statistics (Fig. 8), we choose a horizontal temperature front aligned with the mean wind direction at a particular *z* as our detection event . More complicated events utilizing multiple events are possible, but we believe that a temperature jump with higher (lower) amplitude on the left (right) side of its frontal boundary is a generic feature of the SBL flow as shown previously and, because of its simplicity, does not preordain the type of structure that the conditional sampling will extract from the turbulent flow fields. The properties of the front event are its amplitude , spatial scale *d*, and orientation in an *x*–*y* plane; it has zero mean. To capture a horizontal front we specify two temperature events with equal but opposite amplitudes and ; for example, *E* = +0.1 and −0.1 K on the warm and cool sides, respectively, of the fronts in Fig. 11. Nondimensionalizing this event amplitude in terms of surface-layer variables, , shows that it equates to a scalar flux about 5 times the surface flux using the values for for simulation C in Table 1. The locations of the events define the scale of the front . Although our algorithm allows arbitrary location of the points , these points are initially picked to lie in an *x*–*y* plane at the same *z*. The first linear term in the stochastic estimate of given these two events is then

where the estimation coefficients and are obtained by minimizing the mean-square error of the estimate (Adrian 1996). For our application, and are obtained by solving the 2 × 2 matrix

at all separation points for all variables of interest. Inspection of (8) shows that the most computationally expensive part of the algorithm is calculation of its right-hand side—that is, the unconditional two-point correlations and between the temperature *θ* and the flow variables . These correlations are computed at 256 *z* locations using fast Fourier transforms at 1024^{3 }*x*–*y* points with further averaging over 64 3D volumes collected over the last hour of the simulation. The correlations are archived for later use in (7).

### a. Conditional fields

Conditional fields , velocity gradient tensor , and vorticity fields are estimated for a range of event amplitudes, spatial scales, and vertical locations for all simulations. For discussion, we only present results from simulation C where the reference vertical location . As shown in section 5a, this simulation exhibits numerous temperature fronts at this location above the surface layer. The choice of the event spatial scale *d*, needed by our LSE algorithm, is guided by the flow visualization and the 2D energy spectra for and 2D cospectra between and shown in Fig. 17. These 2D spectra, obtained by averaging in angular rings at constant values of the horizontal wavenumber , are advantageous as they exhibit a peak in the energy and temperature flux (Wyngaard 2010, p. 351). Figure 17 shows that for the SBL the main contribution to the energy and temperature fluxes occurs in the wavenumber range at , with a peak near . This corresponds to spatial scales roughly in the interval .

Figure 18 shows typical 3D vortical structures extracted from the turbulent flow fields using our LSE algorithm with (total jump 0.4 K) and at . These vortices are identified by first estimating followed by computation of the eigenvalues of the velocity gradient tensor . Zhou et al. (1999) show that the imaginary part of the complex eigenvalue; that is, is a measure of the local swirling strength due to a vortex and can be used to identify regions with true flow rotation. The method is well tested and has proven to be an excellent vortex identification technique. In Fig. 18, the isosurface displayed corresponds to a low value of swirling strength of the swirl’s maximum value. To provide a sense of rotational direction, the isosurface of is further colored using the vertical component of the vorticity vector .

Inspection of the 3D image displayed in Fig. 18 from many viewing angles shows a familiar but also unexpected pattern of vortices at this vertical location in the SBL. The downstream vortical structure is identified as a head-up hairpin vortex similar to the structures found in neutral flat wall boundary layers. The strong vertical vorticity in its vertical legs generates a Reynolds stress event (i.e., when the stress is decomposed into four quadrants, this is a Q2 event). Changes in the isosurface contouring, however, show no evidence of our head-up hairpin developing the vigorous streamwise-oriented legs characteristic of near-wall-turbulent boundary layers (e.g., Adrian 2007). At the top of the structure, an arch with spanwise-oriented vorticity connects the vertical legs while at the bottom the vertical legs nearly reconnect, which differs from the vortical structures found near the viscous sublayer in turbulent flow over smooth walls. The primary flow from the head-up hairpin in the SBL is familiar—strong upward vertical velocity combined with significant upstream (negative) horizontal velocity through the legs of the structure resulting in negative momentum flux . Also, the vortex induces horizontal and vertical scalar fluxes and .

The upstream vortical system is slightly weaker in vortical strength compared to its downstream counterpart and even with lengthy averaging its shape is less defined. Based on LSE, we cannot unambiguously claim that it is a head-down or head-up hairpin. Head-down hairpins tend to appear away from boundaries as found by Finnigan et al. (2009) in neutral shear flow over a canopy and also by Gerz et al. (1994) in homogeneous shear flow. Our structure more closely resembles a ring vortex with weak connecting arches at the top and bottom of the structure; Figs. 19 and 20 further support this interpretation. Adrian (1996, see p. 183) found in low–Reynolds number DNS isolated hairpins tend to evolve into ring structures away from the wall. Notice the vertical legs of the upstream structure are rotating in opposite directions compared to its downstream counterpart. Then the induced flow is reversed with strong downward-vertical and positive downstream-horizontal velocities through the center of the vortical structure generating a Q4 Reynolds stress event ; the scalar fluxes are again and .

We emphasize these vortical structures are outputs from our conditional sampling; that is, they are the average flow fields extracted from fully developed turbulent flow fields for a generic temperature front in the SBL. In this aspect our LSE differs from previous applications where Reynolds stress Q2 and Q4 events are simultaneously used as events to identify strong shear layers (e.g., Guezennec 1989). Also recall that the front scale *d* used in the conditioning event lies in the energy containing range of the turbulence and thus the coherent structures are characteristic of the main energy and flux-carrying turbulent eddies.

The upstream and downstream vortices combine to create complex 3D velocity patterns that result in a positive pressure maximum (stagnation point) near the sharp jump in temperature. At the reference height , the horizontal winds shown in Fig. 19 tend to rotate and align with the mean flow direction in the far field but deflect sideways near the stagnation zone. These flow vectors also illustrate the clockwise and counterclockwise rotation of the vertically tilted vortices and provide a sense of the vigorous flow through the legs of the vortical structures. In this *x*–*y* view, one also observes a temperature front rotated in the direction similar to the instantaneous images in Fig. 11. Next, to illustrate the vertical tilt and spatial distribution of the perturbation temperature field relative to the vortical structures, we show and fields in a rotated vertical *x*_{h}–*z* plane; see Fig. 20. This plane is first centered at and then rotated to align with the mean wind direction. The plane cuts midway between the pairs of vortices shown in Fig. 19 and, above and below , slices through the connecting arches bridging the vertical legs of the upstream and downstream vortical structures. Inspection of the figure shows that the upstream vortical structure is indeed a circular ring vortex while the downstream structure is a head-up hairpin vortex. The positive and negative extremes in are found at the centers between the vortex legs. In the interior region between the upstream and downstream vortices, a tilted temperature front develops that spans the same vertical extent as the vortices. For our combination of *d* and *E*, is tilted upward from the axis at an angle of approximately 42° at .

### b. Coherent structure and temperature front variability

Based on the results in section 6a, we conclude that temperature fronts are formed by dynamical coupling with vortical structures located upstream and downstream of a frontal boundary. To further investigate this connection, and also determine the sensitivity to the parameters used in the conditional sampling, we sweep across a broad 2D parameter space spanned by at using LSE. First, because of the linearity of (7), the conditional velocity field and hence the vortex swirling strength both simply increase or decrease in proportion to the event amplitude *E* for a given value of *d*. Thus, an isosurface of normalized by its maximum value does not change. Swirling strength does increase (decrease) with decreasing (increasing) separation distance *d*. Somewhat surprisingly, however, our flow visualization shows that the overall vortical patterns displayed in Fig. 18 remain remarkably robust to varying scale *d*. In other words, even at small and large separations, the upstream and downstream structures are roughly ring and hairpin vortices, respectively. The vortices are tilted at angle of 45° with respect to the horizontal at .

The amplitude, spatial distribution, and tilt of the temperature fronts, however, do depend on the strength and separation of the vortical structures, essentially the *d*–*E* combination. To show this, consider the analog to (6) for the tilt angle of a temperature front based on conditional fields

The prediction from (9) for varying scale *d* is shown in Fig. 21 with event amplitude as a parameter. The tilt angle is shown at the center of the front . Because depends linearly on *E* in (7), the ratio of its horizontal to vertical gradients [i.e., ] depends only on the vortex-scale separation. However, the total tilt of a temperature front, as shown in Fig. 10, includes a mean background temperature gradient that does introduce a dependence on the front amplitude in (9). Examination of the curves in Fig. 21 shows that the tilt angle varies from low values less than 10° for widely spaced vortex systems with weak amplitudes to nearly 60° for closely spaced vortex systems with strong amplitudes. At the extreme end of small *d*, the tilt angle tends to become independent of *E* as the large vertical gradients in overwhelm ; also see Fig. 13. Notice also in this small *d* limit the horizontal gradient of exceeds its vertical gradient. Based on these results, we conclude that the temperature fronts and their variability in tilt observed in the top panel of Fig. 10 are largely created by the dynamical coupling between upstream–downstream vortices with scales in the energy and flux containing range of turbulence.

## 7. Summary

A canonical stably stratified atmospheric boundary layer (SBL) is simulated using high–Reynolds number large-eddy simulation (LES). The problem design is modeled after the first GEWEX Atmospheric Boundary Layer Study (GABLS1) described by Beare et al. (2006). The present set of LES experiments extend the original GABLS1 problem setup by using four different surface cooling rates and three decreasing levels of mesh spacing corresponding to grid meshes with grid points. For each mesh, the spacing is constant across the three coordinate directions. The simulations are integrated for more than 9 physical hours and require more than 900 000 time steps on the finest mesh. The bulk stratification varies from , where is the SBL top and *L* is the Monin–Obukhov length.

The major findings from the study are as follows:

For the posed problem, over the range of stratification considered continuous weakly stratified turbulence is maintained in the SBL both above and below the low-level jet (LLJ) with no global turbulence collapse on the finest LES mesh.

The SBL splits into two regions depending on the height of the LLJ and the surface cooling. Above the LLJ, the turbulence is very weak and the gradient Richardson number is nearly constant at . Below the LLJ, small scales are found to be dynamically important as the shear and buoyancy frequencies squared change with mesh resolution. Both and Ri decrease with increasing mesh resolution. The largest changes are found when the mesh spacing decreases from 2 to 0.78 m.

With increasing stratification the SBL is shallower, the height of the LLJ descends, the winds turn more sharply with height, and the surface wind stress decreases. Also, the mean temperature profile develops sharper vertical gradients in the lower boundary layer and weaker gradients aloft especially so for the simulation with .

Vertical profiles of the Ozmidov scale show a rapid decrease for the strongest cooling rate and over a large fraction of the SBL at high cooling. As a consequence, LES with meshes are likely too coarse to capture the turbulent flow dynamics near the LLJ and near the surface as the cooling increases.

Flow visualization identifies numerous warm–cool temperature fronts in the SBL. The fronts occupy a large vertical fraction of the SBL and tilt forward with increasing stratification. They propagate as coherent entities in space and time. The horizontal and vertical gradients of fluctuating

*θ*are often 10 times larger than the local mean gradient.The LES results can be used to interpret measurements collected from a fixed observational tower. In a height–time (

*z*–*t*) frame of reference, instantaneous vertical profiles of*θ*appear intermittent. Temperature increases in a staircase pattern from the surface up to the location of the LLJ, and between the fronts the temperature is nearly well mixed. Analysis of observations shows these patterns are also found in the temperature profiles in CASES-99 (Poulos et al. 2002). Although is highly intermittent, under heavy space and time averaging, takes on a very uniform vertical structure. Similarly, the average buoyancy frequency*N*^{2}is also smooth and the companion Ri profile varies continuously over the bulk of the SBL. Because of the intermittent character of*θ*, a local Ri(*z*) is inadequate to describe the flow dynamics at short time scales.Conditional sampling based on linear stochastic estimation (Adrian 1996) is used to identify coherent structures in the SBL. The conditioning event is a two-point model of a horizontal temperature front with varying amplitude and scale. For a weakly stratified SBL at

*z*/*z*_{i}= 0.2, we find that the coherent structures are ring and head-up hairpin vortices upstream and downstream, respectively, of the frontal boundary similar to those in neutrally stratified boundary layers. The scale of these structures lies in the energy and flux containing range of the turbulence. Although the vortical structures are oriented at an angle of 45° from the horizontal, they generate temperature fronts of varying amplitude, spatial distribution, and tilt.

## Acknowledgments

PPS and EGP were supported by the National Science Foundation through the National Center for Atmospheric Research (NCAR). PPS also recognizes support from the Office of Naval Research through the Physical Oceanography Program. HJJJ and DVM acknowledge past support from the NCAR Geophysical Turbulence Program. This research benefited greatly from computer resources provided by the NCAR Strategic Capability program managed by the NCAR Computational Information Systems Laboratory (http://n2t.net/ark:/85065/d7wd3xhc). We thank the three anonymous reviewers for their constructive comments.

## REFERENCES

*Eddy Structure Identification*, J. P. Bonnet, Ed., Springer-Verlag, 145–196.

*Studying Turbulence Using Numerical Simulation Databases: Proceedings of the Summer Program 1987*, Stanford University/NASA Ames Research Center, 7–19.

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

*16th Symp. on Boundary Layer and Turbulence*, Portland, ME, Amer. Meteor. Soc., P4.7. [Available online at https://ams.confex.com/ams/BLTAIRSE/techprogram/paper_78656.htm.]

*Atmospheric Turbulence and Mesoscale Meteorology*, E. Federovich, R. Rotunno, and B. Stevens, Eds., Cambridge University Press, 35–49.

*Air-Sea Exchange: Physics, Chemistry, Dynamics, and Statistics*, G. Geernaert, Ed., Kluwer, 507–538.

*Encyclopedia of Atmospheric Sciences*, 2nd ed. G. R. North, F. Zhang, and J. Pyle, Eds., Vol. 4, Academic Press, 232–240.

*IEEE Int. Symp. on Antennas and Propagation*, Orlando, FL, IEEE, 2321–2322, doi:.

*Air Pollution Modeling and Its Application XXI*, S. T. Castelli and D. Steyn, Eds., NATO Science for Peace and Security Series C: Environmental Security, Springer, 57–61, doi:.

*Proc. Seventh Int. Symp. on Turbulence and Shear Flow Phenomena*, Ottawa, ON, Canada, Organizing Committee of the TSFP, 7 pp. [Available online at http://www.tsfp-conference.org/proceedings/2011/6a5p.pdf.]

*Turbulence in the Atmosphere*. Cambridge University Press, 393 pp.