## Abstract

This study seeks to quantitatively and qualitatively understand how stability affects transport in the continuously turbulent stably stratified atmospheric boundary layer, based on a suite of large-eddy simulations. The test cases are based on the one adopted by the Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS) project, but with a largely expanded stability range where the gradient Richardson number (Ri_{g}) reaches up to around 1. The analysis is mainly focused on understanding the modification of turbulent structures and dynamics with increasing stability in order to improve the modeling of the stable atmospheric boundary layer in weather and climate models, a topic addressed in Part II of this work. It is found that at quasi equilibrium, an increase in stability results in stronger vertical gradients of the mean temperature, a lowered low-level jet, a decrease in vertical momentum transport, an increase in vertical buoyancy flux, and a shallower boundary layer. Analysis of coherent turbulent structures using two-point autocorrelation reveals that the autocorrelation of the streamwise velocity is horizontally anisotropic while the autocorrelation of the vertical velocity is relatively isotropic in the horizontal plane and its integral length scale decreases as stability increases. The effects of stability on the overall turbulent kinetic energy (TKE) and its budget terms are also investigated, and it is shown that the authors' large-eddy simulation results are in good agreement with previous experimental findings across varied stabilities. Finally, Nieuwstadt's local-scaling theory is reexamined and it is concluded that the height *z* is not a relevant scaling parameter and should be replaced by a constant length scale away from the surface, indicating that the *z*-less range starts lower than previously assumed.

## 1. Introduction

The stable atmospheric boundary layer (SABL) forms when the underlying land surface is cooler than the air aloft. Typical SABLs, such as the nocturnal boundary layer and the polar boundary layer, significantly influence near-surface as well as large-scale atmospheric dynamics. Consequently, an in-depth understanding of the SABL is required for applications such as numerical weather prediction (NWP) and regional and global climate modeling. However, our current understanding of the SABL is largely hampered by the difficulties and limitations in field observations and numerical simulations of its dynamics. On the sensing side, the reduction in the scale of the turbulent eddies and in the magnitude of the fluxes, and the breakdown or absence of a constant flux layer, affects the quality and representativity of field data (Pahlow et al. 2001; Vickers and Mahrt 2003). On the modeling side, the main problem is that classic turbulence parameterizations and turbulent flow models have been found inadequate under stable conditions, especially under strong stability (Mahrt 1998). However, recent advances in the large-eddy simulation (LES) technique are rapidly making this technique a very valuable tool in the study of the SABL, offering access to unprecedented details of its dynamics.

In numerical modeling of the atmospheric boundary layer (ABL), LES is a major departure from the Reynolds-averaged Navier–Stokes (RANS) approach, which calculates only the mean velocities and the mean temperatures, while parameterizing the turbulent stresses and fluxes using the mean fields. In LES, the turbulent eddies of the size of the numerical mesh grid cells and larger are explicitly resolved, while the effects of the smaller subgrid-scale (SGS) eddies on the resolved ones and the mean fields are parameterized using SGS models. While RANS-type closures will remain a necessity in the near future since weather and climate models cannot afford a resolution that can capture turbulent scales, detailed simulations of the ABL increasingly rely on the LES approach. Unlike the unstable atmospheric boundary layer, which is dominated by eddies of the size of the whole boundary layer, the SABL is populated with small-scale eddies because of the damping of turbulent motions by thermal stratification. The SABL therefore requires higher grid resolutions and more accurate SGS models to achieve realistic simulation results. As such, LES studies of the SABL have been performed only recently, and problems in these studies have been attributed to poor SGS models, low grid resolutions, and numerical instabilities (Derbyshire 1999; Beare and MacVean 2004). Nevertheless, recent studies found that dynamic SGS models can adapt very well to stable conditions and produce correct SGS model coefficients (Kleissl et al. 2004) and that in the critical near-surface region, stability does not increase or alter the role of the subgrid scales (Bou-Zeid et al. 2010). These findings suggest that LES of the SABL can provide accurate results, though careful simulation setup and validation are critical for stable conditions.

Mahrt (1998) distinguished between different stable boundary layers as weakly stable or very stable. The weakly stable case is characterized by weak to moderate stratification and strong wind shear such that a continuously turbulent state can be maintained. However, the strongly stable case is much more complex with the ABL transiently and heterogeneously switching between turbulent and laminar flows. To avoid confusion, in this paper we use the terms “the continuously turbulent SABL” and “the intermittently turbulent SABL” to represent the weakly stable case and the very stable case, respectively. Meanwhile, we generally refer to “weakly” and “very” stable boundary layers to describe different cases and large-eddy simulations that we perform of the continuously turbulent SABL, which is the focus of our study. While the intermittently turbulent stable boundary layer is certainly of practical importance, direct numerical simulations (DNS) would be a more appropriate tool to investigate it. In addition, the continuously turbulent SABL is more ubiquitous and we focus on this scenario in this paper. In fact, recent LES studies of the weakly stable boundary layer have obtained reasonable profiles of turbulent stresses and fluxes that match field measurement results well (Kosović and Curry 2000; Beare et al. 2006; Stoll and Porté-Agel 2008; Zhou and Chow 2011). In addition to verifying the capability of LES in simulating the continuously turbulent SABL, these studies have also quantified turbulent structures and improved turbulence parameterizations in the SABL. These turbulence parameterizations relate gradients and fluxes and have to be reliable across the three vertical regions from bottom to top as described by Mahrt (1998): 1) the surface layer which is the lowest layer and where the turbulent stresses and fluxes are approximately constant with height and the Monin–Obukhov similarity theory (MOST) applies; 2) the variable flux layer where the local turbulent fluxes and dynamics are not in equilibrium with the surface fluxes, which nevertheless remain relevant, and MOST should be modified to form the local similarity theory (Nieuwstadt 1984; Sorbjan 1986b,a); and 3) a layer that is sufficiently constrained by buoyancy where the turbulence is not directly influenced by the fluxes at or distance to the ground; that is, the stratification becomes *z*-less.

As mentioned earlier, the parameterization of the SABL is of crucial importance to large-scale atmospheric models. One of the efforts to improve our ability to model the SABL is the Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS) project launched in 2002 (Holtslag 2003). As a first step of this project, Beare et al. (2006) performed an intercomparison of LES to assess the reliability and sensitivity of the results to resolution and SGS models. The test case used in this intercomparison is the one studied by Kosović and Curry (2000) and in section 2b in this paper we simulate that case and show good agreement between our LES and other simulation and observation results. Building on this, we then modify the test case conditions by increasing the stability of the ABL in order to answer the following questions: 1) How does stability modify bulk dynamics and turbulent structures qualitatively and quantitatively in the SABL? 2) How are the different terms in the turbulent kinetic energy (TKE) budget equation influenced by stability? 3) What can we infer from these analyses about the parameterization of the SABL in climate and weather simulations and can the theoretical model proposed by Nieuwstadt (1984, 1985) be used as a basis for such parameterization across various stabilities? In section 2, we overview our LES code and the SGS model and the test cases we will use. In Section 3, we provide comparison of our LES results against the results in Beare et al. (2006). Then we proceed to investigate the effects of increasing stability on bulk and turbulent structures in section 4 and on TKE budgets in section 5. In section 6 we examine the applicability of Nieuwstadt's model across varied stabilities. Finally, we summarize the findings and their implications for SABL parameterization in section 7.

## 2. Methodology

### a. Large-eddy simulation code

Currently, LES is the state-of-the-art computational tool for the study of high-Reynolds-number flows in a wide variety of environmental settings. Its use for modeling the ABL traces back to the pioneering work of Deardorff (1972), and since then, developments in LES have been proliferating. The basic premise in LES is that large-scale turbulent eddies contribute the majority of turbulent fluxes, and the effect of small-scale turbulent eddies on the large scales can be parameterized with sufficient reliability, such that 3D high-Reynolds-number turbulent flows can be accurately simulated.

The LES code used in this study is a modified version of the code developed by Albertson and Parlange (1999a,b) and Bou-Zeid et al. (2005). This code solves the filtered Navier–Stokes equations written in rotational form (Orszag and Pao 1974). We use the scale-dependent Lagrangian dynamic SGS model as proposed and validated by Bou-Zeid et al. (2005); this model follows from the works of Germano et al. (1991), Meneveau et al. (1996), and Porté-Agel et al. (2000). As its name suggests, this SGS model has several advantageous features compared to the most widely used Smagorinsky–Lilly model (Smagorinsky 1963; Lilly 1967): 1) the SGS eddy viscosity is calculated dynamically, thus avoiding the need for ad hoc or tunable coefficients in the model; 2) the averaging of the model coefficients, which is required to stabilize dynamic SGS models numerically, is performed over fluid pathlines (the Lagrangian approach) to preserve local variability in the eddy viscosity, which allows the simulation of complex flows over heterogeneous boundary conditions; 3) the scale dependence of the model coefficient near solid boundaries is accounted for to better capture land–atmosphere fluxes and the dynamics of the surface layer. A similar approach was also recently implemented for scalar SGS fluxes (Stoll and Porté-Agel 2008) and we have adapted this approach for our LES. The approach mimics the dynamic computation of the Smagorinsky coefficient [as presented in Bou-Zeid et al. (2005)] to compute the equivalent coefficient for scalar SGS fluxes [see details in Shah and Bou-Zeid (2010)]. The code we use here thus uses the scale-dependent Lagrangian dynamic approach for both momentum and scalar SGS fluxes.

A pseudospectral approach is used to compute the horizontal derivatives such that the horizontal boundary conditions are periodic. To calculate the vertical derivatives, a staggered vertical grid is used to allow the implementation of a second-order accurate centered differences scheme. Full dealiasing of the nonlinear terms is done using the rule (Orszag 1970, 1971), and the Coriolis force is included. We also note that the horizontal pressure gradients driving the flow are represented using a geostrophic wind.

### b. Description of the simulation cases

In this study, we simulate a series of idealized SABL cases based on the one described in Kosović and Curry (2000) who imposed a steady surface cooling rate, but we also simulate cases with significantly higher cooling rates. The initial and boundary conditions imposed in Kosović and Curry (2000) were made consistent with the Beaufort Sea Arctic Stratus Experiment (BASE) dataset to resemble a clear-air SABL driven by a moderate surface cooling rate. With slight modifications, this case was also adopted by the GABLS LES intercomparison project (Beare et al. 2006) and other LES studies (e.g., Basu and Porté-Agel 2006; Stoll and Porté-Agel 2008) mainly because it is simple and well documented. The land surface is horizontally homogeneous with the roughness length of momentum set as *z*_{0m} = 0.1 m. Note that *z*_{0m} is enhanced from its typical value for sea ice in the Arctic ocean to avoid an underresolved surface layer (Kosović and Curry 2000). As such, the simulations are more representative of a nocturnal ABL over farmland. The roughness length of heat *z*_{0h} is set the same as *z*_{0m} for consistency with Beare et al. (2006) although we note that *z*_{0h} is typically an order of magnitude smaller than *z*_{0m} for such hydrodynamically rough walls (see Brutsaert 2005, 45–46). The initial mean temperature is 265 K up to 100 m with an overlying inversion of strength 0.01 K m^{−1}. A constant geostrophic wind of *U _{g}* = 8 m s

^{−1}is imposed, with a Coriolis parameter of 1.39 × 10

^{−4}s

^{−1}, which corresponds to latitude 73°N. The initial mean wind was geostrophic. Stress-free and no penetration conditions are imposed at the top of the computation domain; that is, ∂

_{3}

*u*

_{1,2}=

*u*

_{3}= 0, where 1, 2, and 3 (or

*x*,

*y*, and

*z*) refer to the streamwise, spanwise, and vertical directions, respectively. Monin–Obukhov similarity theory is applied at the surface as a wall model, with the coefficients consistent with Beare et al. (2006) (i.e.,

*β*= 4.8 and

_{m}*β*= 7.8). Other parameters include the von Kármán constant

_{h}*κ*= 0.4 and gravitational acceleration

*g*= 9.81 m s

^{−2}.

The computational domain size is fixed at 800 m × 800 m × 400 m, consistent with Stoll and Porté-Agel (2008), but the shallow depth of the SABL ensures that our horizontal domain scale is larger than 4–5 times the largest eddy size (Moeng et al. 2007). The simulations are run with 80^{3} gridpoints for the first 6 h and then the resulting outputs are interpolated to 162 × 162 × 160 gridpoints and run for another 4 h. The simulations are first allowed to reach quasi-steady state and to develop the smaller turbulent scales at the higher resolution during the first hour, and only the last 3 h are used to compute the statistics of the SABL. The 162 × 162 × 160 gridpoint resolution was selected after numerical convergence tests showed that the results, for all the stabilities simulated here, did not vary much as this resolution was approached. It results in grid cell dimensions of approximately 5 m × 5 m × 2.5 m, which in our LES also represents the dimensions of the filter.

The full signal of a turbulent variable can be decomposed into a Reynolds average (represented by uppercase letters) and a turbulent part (represented by a prime) or into a resolved part (represented by a tilde) and an SGS part (represented by double primes) (e.g., ). Note that in LES we only solve for the resolved fields and hence we directly obtain only the resolved part of the statistical moments. For the first-order moments, the mean of the resolved part is in practice a very accurate estimates of the mean of the total (resolved + SGS). For second-order moments, we can model the SGS components, which are important in this case, to recover the totals for some important quantities (e.g., stresses and scalar fluxes). As such, in this paper, we use the term “mean” to refer to the mean of the total, including the modeled SGS part for the second-order moments only, unless otherwise stated.

Six steady surface cooling rates are prescribed to simulate increasingly stable ABLs, including the one studied by the GABLS project—namely, −0.25, −0.5, −1, −1.5, −2, and −2.5 K h^{−1}, denoted as cases A–F, respectively. The statistics are computed at quasi equilibrium, defined as the state of the SABL where the height of that layer and the surface fluxes of momentum and heat change relatively slowly with time, allowing turbulence to be close to equilibrium at any instant. A list of the primary mean characteristics at quasi equilibrium are tabulated in Table 1, which includes the ABL height *h*, friction velocity *u*_{*} = (−*τ _{s}*)

^{1/2}, surface temperature scale , the Obukhov length at the surface , and the Ozmidov length at the surface ;

*τ*represents the kinematic surface stress,

_{s}*q*the surface buoyancy flux, Θ

_{s}_{s}the surface potential temperature, ε the dissipation rate (we will detail how we compute it in section 5), and

*N*

_{BV}= [(

*g*/Θ)(

*d*Θ/

*dz*)]

^{1/2}the Brunt–Väisälä frequency. We calculate

*h*following Beare et al. (2006); that is, it is the height where the mean stress falls to 5% of its surface value divided by 0.95.

## 3. Validation of LES results

The LES code with the dynamic SGS model has been tested extensively for ABL flows over homogeneous and heterogeneous surfaces (Bou-Zeid et al. 2004, 2005), urban flows and wind tunnel flows over cubes (Tseng et al. 2006), and flow in plant canopies (Yue et al. 2007b). The code was also validated for stable and unstable ABL flows and for diurnal cycles (Kumar et al. 2006; Kleissl et al. 2006; Kumar et al. 2010) for values of the flux Richardson number at the surface up to about 4. These validations analyzed mean velocity profiles, Smagorinsky coefficient values, stress profiles, and other statistics that confirm that this LES code can realistically reproduce observed and theoretical ABL statistics. Yue et al. (2007a) compared the LES results to PIV data for flow in a plant canopy and concluded that the quadrant analysis results from the LES matched PIV results very well, suggesting that the coherent structures, which directly affect the quadrant analysis statistics, are well represented.

However, given the particular difficulties in modeling the SABL, we conduct in this section further tests comparing our LES results against those in Beare et al. (2006), which corresponds to case A and which were compared to observational data, before increasing stability and investigating how it modifies turbulence structures within the SABL. We will also later in the paper show direct comparison of the LES results with field experimental data for buoyant TKE destruction under all the simulated stabilities. The time step is set such that the maximum Courant–Friedrichs–Lewy (CFL) number is around 0.1, which is the low value needed for numerical accuracy and stability of the code. In Figs. 1 and 2, we compare the major results of case A at four different resolutions with those obtained by the LES models in Beare et al. (2006) (the horizontal resolutions of our simulations are increased proportionally to our vertical resolution, though only the vertical grid spacing Δ*z* is listed in the figures). Note that all the mean statistics at quasi equilibrium in this paper are obtained by averaging over horizontal planes and over the last 3 h except for the mean potential temperature, which we only average over hour 9 to be consistent with Beare et al. (2006). Beare et al. (2006) concluded that a resolution of Δ*z* = 3.125 m or less is ideal for simulating such a moderately stable ABL, while acceptable performance can still be achieved with a resolution of Δ*z* = 6.25 m. Basu and Porté-Agel (2006) showed that the scale-dependent Lagrangian dynamic SGS model has the advantage of producing relatively resolution-insensitive results for the moderately stable ABL. Here we show that all the vertical profiles are converging as the resolution increases for case A. For *τ _{s}* the difference between Δ

*z*= 3.125 m and Δ

*z*= 2.5 m is about 3%, and for

*q*it is about 4%. The convergence test has also been performed for the other cases and similar results have been obtained. As such, we will present LES results of varying stabilities at the highest resolution of Δ

_{s}*z*= 2.5 m for the rest of the paper.

It is shown in Figs. 1 and 2 that our LES generally falls in or near the range of the corresponding results produced by the LES models employed in Beare et al. (2006). A close examination of this comparison reveals that our results are in better agreement with those produced by the Met Office (MO) model (Beare and MacVean 2004), the University of Hannover and Yonsei University (IMUK) model (Raasch and Schröher 2001), and the Universitat de les Illes Balears (UIB) model (Cuxart et al. 2000) than the other models employed in the LES intercomparison of Beare et al. (2006), probably because they used finer resolutions than the rest. Based on Fig. 2, we observe that a quasi-equilibrium state has been reached approximately after hour 7, as evidenced by the plateaus. In Fig. 3 we compare the results obtained using the new dynamic Pr_{sgs} model (Shah and Bou-Zeid 2010) and the static Pr_{sgs} (=0.6) model for cases A and F. An almost perfect match is found for case A between the dynamic and static Pr_{sgs} models for the profiles of the mean quantities as well as the turbulent fluxes (lines are indistinguishable on the figure); case F with higher stability shows some differences that remain however small. In the rest of this paper we will use the dynamic Pr_{sgs} since it should be able to better adapt to the varying flow stability for different cases and at different heights.

Finally, given the higher stabilities we simulate in this study, we plot the *u* spectra, scaled by *U _{g}* and

*z*, for case A and case E in Fig. 4. As stability is strengthened in the SABL, the highest frequencies, corresponding to the smallest resolved scales, still display a − slope. This indicates that they have inertial-subrange dynamics and that our cutoff grid scale reaches into the inertial subrange, except at the very first one or two grid points because of proximity to the wall rather than stability. This can also be deduced by comparing our grid cell size Δ to the Ozmidov-scale

*L*, which can be interpreted as the smallest scale that is affected by buoyancy (Bou-Zeid et al. 2010). Even for the highest stabilities

_{OZ}*L*≈ Δ, indicating that the SGS range and the smallest resolved scales are not strongly feeling the stabilizing effect of buoyancy. The inertial subrange however becomes narrower as the effect of stability reaches smaller and smaller scales, and the energy of the production range eddies (large scales) decreases markedly. The SGS model hence takes a more important role in determining the turbulent dynamics (in a later section we will analyze the fraction of resolved versus SGS energies). In addition, close to the wall, the SGS model still needs to represent the effects of unresolved production range eddies, which are now modified by buoyancy. Hence, despite the fact that the filter scale lies in the inertial subrange except near the wall, it is critical for the SGS model to drain the right amount of SGS fluxes and TKE from the resolved scales and to produce the correct fluxes. The scale-dependent Lagrangian dynamic SGS model that we use here has been shown to successfully accomplish this, even in the near-wall region where the cutoff filter scale is in the production range, and generally performs much better than the traditional Smagorinsky–Lilly model and the Lagrangian scale-invariant model (Bou-Zeid et al. 2005). Furthermore, Kleissl et al. (2004) have shown that the scale-dependent dynamic procedure is able to yield the optimal model coefficients near the wall (in the production range) as the stability of the surface layer is varied. Overall, based on the results presented in this section, we expect the LES to produce reliable simulations of the SABL for all the stabilities we simulate here.

_{OZ}## 4. Effects of stability on bulk dynamics and turbulence structures

We now proceed to investigate how stability affects bulk profiles and turbulence structures in the SABL. In Fig. 5 we contrast the vertical profiles of mean wind speed, mean potential temperature, mean stress, and mean heat flux of cases A–F at quasi equilibrium. A preliminary examination of Table 1 and Fig. 5 suggests, as expected, that an increase in stability results in a shallower boundary layer, a decrease in stress (due to a decrease in the momentum transfer from the top of the ABL), an increase in vertical temperature gradients and consequently in heat flux, as well as a lowered low-level jet (LLJ). Based on measurement data from the Microfronts project (Sun 1999), Mahrt (1998) proposed to use *h*/*L _{MO}* ≈ 1 as the criteria to distinguish the continuously turbulent SABL (

*h*/

*L*< 1) from the intermittently turbulent one, where in the latter regime turbulence tends to exhibit strong intermittency or to totally collapse. However, Flores and Riley (2011) reviewed other similar studies using experimental, laboratory, and numerical simulation data and concluded that such a parameter is not appropriate to characterize turbulence collapse. In our study

_{MO}*h*/

*L*ranges from 1.5 for case A to 9.1 for case F and no collapse or global intermittency of turbulence have been observed (global intermittency refers to the absence of turbulence everywhere in the domain), which is consistent with the perspective of Flores and Riley (2011). To illustrate this point, the time histories of

_{MO}*τ*,

_{s}*q*and

_{s}*h*are plotted for the six cases in Fig. 6. As with case A, cases B–F have all achieved quasi equilibrium after hour 7, if not earlier, with no sign of global intermittency. This is a critical aspect of the analysis since it has important consequences for the parameterization of the stable boundary layer, which according to our results can continue to rely on models of the continuously turbulent SABL even when the gradient Richardson number (Ri

_{g}) exceeds the critical value of 0.25 (cf. Fig. 13). In fact, a number of data sources from measurements (e.g., Poulos et al. 2002; Mahrt and Vickers 2005) and numerical simulations (e.g., Zilitinkevich et al. 2007, 2008) have shown that turbulent mixing persists even for Ri

_{g}> 1. However, one observation that we made in our analysis is that higher stabilities produce instantaneous turbulence fields that are more heterogeneous in space, with “hot spots” of turbulent activity and fluxes that are easier to delineate from the more quiescent regions of the flow.

After our analysis of the mean profile, we proceed to study the dynamics and scales of the large coherent turbulent structures and how they are affected by stability. An important statistic to quantitatively characterize the spatial features of coherent structures is the two-point spatial correlation function between turbulent velocities at an arbitrary reference location (*x*, *y*, *z*) and another roving location . For horizontally homogeneous turbulence, its definition can be expressed as

where the angle brackets represent horizontal and temporal averaging. Figure 7 depicts contour plots of *r*_{11}(*r _{x}*,

*r*, 0.25

_{y}*h*, 0.25

*h*) and

*r*

_{11}(

*r*,

_{x}*r*, 0.5

_{y}*h*, 0.5

*h*) for cases A and case F, respectively. Note that both

*x*and

*y*axes are scaled by

*h*such that the dimensional sizes of the contour patterns are not directly comparable between case A and case F. Also note that the

*y*symmetry is imposed since the boundary conditions of the simulations are periodic in the

*y*direction and the Navier–Stokes equations are invariant with respect to a reflection in

*y*. Isocontours of

*r*

_{11}appear as ellipses with their major axis aligned in the streamwise direction. For case A at

*z*= 0.25

*h*the coherent structure defined by

*r*

_{11}> 0.1 has a scale of about

*h*in

*x*and 0.3

*h*in

*y*. We also observe regions of weakly negative

*r*

_{11}centered at

*r*/

_{y}*h*≈ ±0.3, indicating the existence of vortical structures, aligned with the mean velocity, at this vertical level. These vortical structures will produce negative correlations between the streamwise velocities at their subsiding side (high

*u*for air parcels coming from aloft) and rising side (low

*u*for air parcels coming from below). For case A at

*z*= 0.5

*h*, the coherent structure defined by

*r*

_{11}> 0.1 grows to about 1.5

*h*in

*x*and 0.5

*h*in

*y*. Compared to case A, the size of the coherent structure defined by

*r*

_{11}> 0.1 increases significantly for case F at the same

*z*/

*h*—it extends to about 4.5

*h*in

*x*and

*h*in

*y*at

*z*/

*h*= 0.25, and even larger at

*z*/

*h*= 0.5. In Fig. 8, a similar analysis is depicted for

*r*

_{33}, where we compare case A and case B since the coherent structure revealed through

*r*

_{33}is too small for cases under higher stabilities. For these higher stabilities, the instantaneous turbulence fields are heterogeneous and the number of active structures present in a given instantaneous field is limited (as revealed by vertical velocity plots over an

*x*–

*y*plane not shown here). Hence, the correlation decreases quickly because of the presence of large quiescent areas with few active eddies and almost zero correlation. Compared with the results of

*r*

_{11}in Fig. 7, the size of the coherent structure defined by

*r*

_{33}> 0.1 is much smaller (note that the axes in Figs. 7 and 8 are not equal), but its scale also increases as

*z*/

*h*increases from 0.25 to 0.5. However, unlike

*r*

_{11},

*r*

_{33}is less anisotropic in the

*x*–

*y*plane. In addition, the size of the coherent structure defined by

*r*

_{33}decreases with stability.

The analysis of the characteristic size of a typical coherent structure, at a given height, using two-point correlation data in Figs. 7 and 8 is relevant for our aim of improving the parameterization of the ABL since it can be used to calculate integral length scales and how they vary with height. Following Shaw et al. (1995) and Huang et al. (2009), the longitudinal integral length scale is obtained by integrating the correlation tensor of zero lag in *y* and *z* over *r _{x}*:

where no summation is implied by the repeated index *ii*. It is clear from this definition that *L _{i}*(

*z*) characterizes the longitudinal extent of the coherent structures at

*z*for any velocity component

*i*. In Fig. 9, the vertical profiles of

*L*

_{1}and

*L*

_{3}for cases A–F are shown. Note that instead of monotonically increasing with stability,

*L*

_{1}tends to decrease slowly with stability from case A to case C in the lower SABL, and then increases sharply from case C to cases D–F. It is interesting to note that

*L*

_{1}for case D–F appears to be unimodal with the peak ranging between

*z*/

*h*= 0.4 and

*z*/

*h*= 0.5. We also see that

*L*

_{3}tends to monotonically decrease with stability and increase with height, and the increase with height is more rapid in the upper SABL. These nonmonotonic behaviors result from the range of complex dynamics that are influencing variance and autocorrelation in the stable ABL. From cases A–C, the decrease in both length scales is most likely related to an increase in the buoyant destruction of TKE (see Fig. 11); this affects all variance components and the scale of the large eddies. For higher stabilities, it is likely that a decrease in pressure redistribution (main production term, −〈

*u*′

*w*′〉∂

*U*/∂

*z*, is in the

*u*-variance budget) and/or a weakened coupling between different layers are allowing an increase in the integral scale based on

*u*, while the one based on

*w*is damped by buoyancy and reduced redistribution. The results indicate that there is a critical point, as stability increases from case C to case D, beyond which coherent turbulent structures become rather very flat in the middle SABL.

The sharp decrease of *L*_{3} for cases D–F is primarily a consequence of the reduction in the number of active eddies, and also implies a reduction in the size of these eddies. For these cases, active eddies occupy a small fraction of the domain and can still be observed in LES fields (and thus are resolved by the computational grid). A computation of *L*_{3} over areas with active fluxes only, or based on some conditional averaging technique, might be more suited since it would exclude the contribution of the large nonturbulent regions. However, we do not perform this computation here and we rely on the classic definition of the integral-scale *L*_{3} since it still captures the variability of the length scale of active turbulence with stability, and is more appropriate than *L*_{1} because *w* is directly related to the active turbulence that results in vertical transfer of momentum and scalars. In contrast, *u* is associated with the larger-scale quasi-horizontal motions that do not make significant contributions to vertical transfer (e.g., Townsend 1976; Perry et al. 1986; Raupach et al. 1996). The variability trends of *L*_{3} illustrate that the vertical mixing length for an ABL model needs to increase nonlinearly with height and decrease with stability.

## 5. Effects of stability on TKE budgets

We now proceed to investigate the effects of stability on total TKE budgets. Since in LES of incompressible flows the SGS TKE (*E*_{sgs}) is not directly calculated, this term has to be modeled. By using the − energy spectrum relationship in the inertial subrange, Knaepen et al. (2002) were able to relate the trace of the Leonard stress tensor, , to *E*_{sgs}:

where the tilde and the overbar denote the grid filter scale and the test filter scale (here taken as 2 × the grid filter scale), respectively, and Δ represents the width of the filters. We use this model because of its simplicity and relatively good performance and because our SGS range is not strongly affected by buoyancy as previously discussed [see discussion and more accurate models for stable and unstable turbulence in Salesky and Chamecki (2011)]. In Fig. 10 the vertical profiles of the resolved TKE (*E*_{res}), *E*_{sgs}, the total TKE (*E* = *E*_{res} + *E*_{sgs}), and the fraction of *E*_{sgs} in *E* are shown. Note that the results for the lowest vertical levels are not shown in Fig. 10 since near the surface and will lie in the production range of the energy spectrum (cf. Fig. 4), and consequently, Eq. (3) will give erroneous results. We should not expect *E*_{res} and *E*_{sgs} to be exactly zero at *z* = *h* because *h* is defined based on stress rather than TKE. It is clear that *E*_{res} generally decreases with stability as expected; this decrease is more significant from case A to case C than from case C to case F. However, *E*_{sgs} is rather insensitive to stability, which agrees with the fact that buoyancy primarily affects the resolved scales and their TKE in our SABL simulations, rather than the subgrid scales; this is in agreement with the findings of Bou-Zeid et al. (2010). Consequently, the trend of *E* with stability is in accord with that of *E*_{res}. Taking the values at *z* = 0.5*h* as an example, *E* decreases by approximately 42% from case A to case C and 14% from case C to case F. This decrease is associated with the shrinking of the inertial subrange that is now squeezed by the buoyancy subrange, as evidenced by the decrease of *L _{OZ}* with stability illustrated in Table 1 and by the spectra plotted in Fig. 4. Finally, the decrease of

*E*

_{res}with, and insensitivity of

*E*

_{sgs}to, stability lead to the general increase of

*E*

_{sgs}/

*E*with stability, although from case C to case F this ratio tends to converge. At

*z*= 0.5

*h*,

*E*

_{sgs}is only about 24% of

*E*for case A and this value increases to around 40% for case F.

For the stably stratified surface layer over horizontally homogeneous surface at quasi equilibrium, the pressure and turbulent transport of TKE are typically negligible such that the mean TKE budget equation roughly reduces to a simple balance between the mean mechanical production 〈*P*〉 as a source, and the mean buoyant destruction 〈*B*〉 and the mean viscous dissipation 〈ε〉 as sinks (Wyngaard 1992), which gives

where 〈*ξ*〉 is the total mean TKE destruction by buoyant and viscous effects and 〈*B*〉 consists of the SGS part and the resolved part *B*_{res} = −*g*Θ^{−1}〈*w*′*θ*′〉; that is,

Following Sullivan et al. (1994) and Bou-Zeid et al. (2010), and since production of TKE by interaction with the mean flow and TKE transport are negligible for the SGS range, 〈ε〉 is expressed as a balance between the energy flux across the filter-scale and the SGS buoyant destruction *B*_{sgs}:

We calculate 〈*P*〉 as

where *τ _{uw}* is the total stress between

*u*and

*w*and

*τ*is the total stress between

_{υw}*υ*and

*w*. In Fig. 11 the vertical profiles of 〈

*P*〉, 〈

*B*〉, 〈ε〉, and 〈

*B*〉/〈

*ξ*〉 are plotted. Note that 〈

*B*〉 increases with increasing stability across the entire SABL. However, stability affects 〈

*P*〉 and 〈ε〉 more noticeably in the surface layer (roughly

*z*< 0.1

*h*) than in the outer layer: an increasing stability generally results in a decreasing 〈

*P*〉 and 〈ε〉 in the surface layer; however, farther away from the surface, this trend continues for 〈ε〉 but reverses for 〈

*P*〉. An increase in 〈

*P*〉 with stability can be noted around the location of the LLJ as depicted in Fig. 5 owing to narrower LLJ and higher resulting shears in the more stable cases. This increase in the middle of the SABL of 〈

*P*〉, which is first produced mainly as

*u*variance, is likely contributing the strong increase in the integral scale based on

*u*depicted in Fig. 9 for cases D–F. The fraction of 〈

*B*〉 in 〈

*ξ*〉 increases with stability across the entire SABL as well. For

*z*≈ 0.1

*h*, 〈

*B*〉/〈

*ξ*〉 ranges from less than 10% for case A to over 30% for case F. We see that 〈

*B*〉 decreases approximately linearly with height as expected since it is proportional to the buoyancy flux (this will be revisited in section 6), but clearly, 〈

*P*〉 and 〈ε〉 decrease more rapidly with height; this results in an increase of 〈

*B*〉/〈

*ξ*〉 with height. At the top of the SABL, 〈

*B*〉 contributes about 32% of 〈

*ξ*〉 for case A and 84% for case F. Bou-Zeid et al. (2010) studied the effects of stability on the relative importance of 〈ε〉, 〈

*B*

_{res}〉, and 〈

*B*

_{sgs}〉 in contributing to the total mean TKE dissipation 〈

*ξ*〉 using field experimental data collected over an extensive glacier in Switzerland and found that at high stabilities in the stable surface layer 〈ε〉 and 〈

*B*

_{res}〉 become equally important, while 〈

*B*

_{sgs}〉 also grows to be a significant fraction of total mean TKE destruction. In Fig. 12, we reproduce Fig. 9 of Bou-Zeid et al. (2010) and compare to our LES results. Note that instead of using a scatterplot as in their original figure, here we plot the bin-averaged values of Bou-Zeid et al. (2010) with the error bars representing the standard deviations in the corresponding bin. We use the data at the first full vertical node (

*z*= 2.5 m) to ensure that the field data and the LES data share approximately the same Δ/

*z*(≈2), because this ratio controls the fraction of turbulent fluxes falling into the SGS range. The match between the field data and the LES data is generally good with the LES data falling in or near the range of the error bar of the field data; some discrepancy is expected since Δ/

*z*is not exactly matched. This is an important validation result for our code that supplements the previous LES comparison. It is very encouraging that LES can reproduce these fractions rather well.

## 6. Examination of Nieuwstadt's model

A major departure, of particular relevance for turbulence parameterization, of the SABL from the unstable ABL is associated with the size of turbulent eddies: in the SABL turbulent eddies are relatively small and localized because of the suppression of vertical motion by buoyancy while in the unstable ABL turbulent eddies are of the scale of the entire ABL. This scale reduction leads to an increase in the time scale the SABL requires to reach equilibrium with surface fluxes, such that these surface fluxes are no longer the main determinants of ASL dynamics as is assumed in the Monin–Obukhov similarity theory. Inspired by this, Nieuwstadt (1984, 1985) proposed a local-scaling theory for the SABL. This theory states that for stably stratified turbulent flows, dimensionless combinations of variables (e.g., variances, covariances, and eddy viscosity coefficients) measured at a specific vertical level are functions of a sole parameter *z*/Λ at this level, where Λ is the local Obukhov length (rather than the surface one) defined as −*τ*^{3/2}Θ(*κgq*)^{−1}, where *τ* is the local total stress, *q* is the local total buoyancy flux, *κ* = 0.4 is the von Kármán constant, *g* is the gravitational acceleration, and Θ is the mean potential temperature at that height. Surface fluxes are still relevant in the local approach since they influence the local fluxes, and the height *z* thus remains a relevant parameter. In addition, in the limit of *z*/Λ → ∞, all dimensionless forms of the variables derived from the local-scaling assumption approach a constant value; that limit is called the *z*-less stratification since variables no longer depend on the distance to the surface, or on surface fluxes, under such high stabilities. Since some SABL parameterizations are based on this local-scaling theory and the *z*-less limit, in this section, we examine the effects of stability on the applicability of Nieuwstadt's theory. An assumption for this theory to hold is that both Ri_{g} and the flux Richardson number (Ri_{f}) are constant with height (≈0.2), where

Physically, Ri_{f} is interpreted as the ratio between 〈*B*〉 and 〈*P*〉 in the TKE budget equation, and Ri_{g} is related to Ri_{f} through the gradient-diffusion hypothesis. In Fig. 13, we show the vertical profiles of Ri_{g}, Ri_{f}, and turbulent Prandtl number (Pr), which can be conveniently calculated as Pr = *K _{m}*/

*K*= Ri

_{h}_{g}/Ri

_{f}, where

*K*and

_{m}*K*are turbulent diffusivities for momentum and heat flux, respectively. Unfortunately, the assumption of constant Ri

_{h}_{g}and Ri

_{f}is violated for all cases studied here—they increase monotonically with height and increasing stability. Around

*z*= 0.2

*h*, Ri

_{g}ranges from 0.07 for case A to 0.23 for case F while it is 0.25 for case A and 0.91 for case F around the top of the SABL. As such, no asymptotic Richardson number has been observed as far as our study cases are concerned. Stability has pronounced effects on Pr:

Pr generally ranges between 0.6 and 0.8 except close to the surface.

Pr becomes more dependent on height as stability increases. For cases A–C,

*Pr*is more or less constant in the middle and upper SABL. However, Pr starts to vary with*z*for case D, with its minimum located around*z*= 0.2*h*, and this variability increases from case D to case F. This suggests that, under very stable conditions, the transport efficiency of momentum flux relative to heat flux first decreases and then increases with height (as the effect of the proximity to the wall decreases and the effect of stability starts to dominate at higher elevations); this is consistent with the finding in Howell and Sun (1999), based on experimental data, that Pr at 3 m is generally higher than that at 10-m level.Pr tends to increase with increasing stability, which is in agreement with the modeling results of Venayagamoorthy and Stretch (2010) and with field-observed trends for the SGS Pr (Bou-Zeid et al. 2010).

A parameter that is more directly relevant for SABL parameterizations related to the local-scaling theory is the turbulent diffusivity. Following Nieuwstadt (1984), *K _{m}* and

*K*are nondimensionalized as

_{h}and, similarly,

The left panels in Fig. 14 depict and as functions of *z*/Λ for cases A–F, along with the LES results at a 2-m resolution in Beare et al. (2006), which correspond only to case A in our simulation. Beare et al. (2006) claim a good agreement arising from the comparison between the LES results they collected and the local-scaling theory. Although our results of case A do generally fall in the range of the LES results reported in Beare et al. (2006), it is clear from Fig. 14 that the agreement of Beare et al. (2006) with the local-scaling theory reflects the limited range of *z*/Λ and stabilities tested in that study. Instead of approaching a constant value, both and monotonically increase with increasing *z*/Λ and stability; and the collapse of the curves when plotted against *z*/Λ is only observed near the surface where *z*/Λ is small. In fact, the variability of and is directly associated with the violation of the constant Richardson number assumption as = *κ*Ri_{f} and = *κ*Ri_{f}/Pr. In contrast, we define a constant length-scale *λ* and plot and against *λ*/Λ in the right panels of Fig. 14 [assuming *λ* = 1 m for now; in Part II of this study (Huang et al. 2013) an appropriate length scale is determined], and we find that and collapse well with *λ*/Λ except close to the ground where *λ*/Λ is relatively small. This suggests that far from the land surface in the SABL, the *z*-less limit of the local-scaling theory applies well for the stability range we impose here, *z* is not a relevant length scale anymore, and its inclusion in eddy diffusivity models away from the surface can be problematic since it does not result in universal functions. Consequently, *z*/Λ is an appropriate dimensionless parameter for the SABL near the surface. Farther aloft, *z*/Λ is not appropriate for the observed *z*-less scaling behavior, and instead, our LES data indicate that *z* should be replaced by a constant length-scale *λ* to explain it. This issue is further explored in Huang et al. (2013), where the value of *λ* is determined and these findings are used to provide better parameterizations of the SABL.

In addition to the local-scaling argument, Nieuwstadt (1984, 1985) also derived an expression for the vertical profiles of stress and temperature flux at quasi equilibrium:

and

Figure 15 examines these predictions and compares them to our LES results for different stabilities. Note that an exact match between Nieuwstadt's theory and the LES at the top of the SABL should not be expected because *h* is not defined through zero stress or zero temperature flux. Other than that, the agreement between the prediction and the LES results is generally good for both stress and temperature flux at all stabilities, in agreement with what has been reported in Beare et al. (2006). It is not surprising that the violation of the assumption of constant Richardson number does not affect the applicability of Eqs. (12) and (13) significantly as it does for the *z*-less scaling. A careful examination of the derivation process in Nieuwstadt (1984) reveals that in addition to the boundary conditions, the only necessary assumption for Eqs. (12) and (13) to hold is that ∂Θ/∂*z* and ∂*U*/∂*z* are time insensitive, which is a much weaker condition than the constant-Richardson-number assumption and is generally valid at quasi equilibrium.

## 7. Summary

In this study, fine-resolution LES of the SABL are performed for a wide range of stabilities using the scale-dependent Lagrangian dynamic SGS model. The LES is first validated using an established weakly stable case. Then, six cases with increasing stabilities are studied by prescribing steady surface cooling rates ranging from −0.25 to −2.5 K h^{−1}. After approximately 7 h of simulation, all of the cases reach a quasi-equilibrium state in which a continuous turbulent state is maintained. The results are validated and analyzed with the aim of quantitatively and qualitatively investigating the effects of stability on bulk profiles and turbulent structures, and consequently on energy, momentum, and heat transfer in this continuously turbulent SABL. The focus is on the properties of the SABL that are critical for successful parameterization in weather and climate models, specifically the turbulent length scale, the variation of the TKE with stability, and the validity of the widely used local similarity proposed by Nieuwstadt.

The results reveal that an increasing stability in the SABL generally gives rise to 1) a mean vertical temperature profile with a stronger gradient, 2) a mean velocity profile with a lowered low-level jet, 3) a decreased momentum transfer, 4) an increased heat flux, and 5) a lowered boundary layer height at quasi equilibrium. Furthermore, as stability increases, its effects on the mean flow velocity and the magnitude of the turbulent fluxes become less significant. Two-point correlation analysis is used to illustrate coherent turbulent structures in the horizontal plane and calculate the integral length scale. Stability affects the two-point correlation of *u* in a different way from that of *w*; the effect on *w* is more important for SABL parameterization since *w* is associated with active turbulence that results in vertical transport, while *u* is not. The pattern of the *u* correlation is anisotropic in the horizontal plane, with the streamwise extent typically over 3 times larger than the spanwise extent. An increasing stability tends to largely increase the size of the pattern of the *u* correlation. The integral length scale based on *u* appears rather uniform with height for case A, but changes with height at stronger stabilities with a peak around *z*/*h* = 0.4. However, the pattern of the *w* correlation is approximately isotropic in the horizontal plane and, contrary to its *u* counterpart, the *w* correlation is reduced as stability is strengthened. The integral length scale based on *w* generally increases with height monotonically. These results suggest that while increasing stability normally decreases the size of coherent turbulent structures, there is a critical point in stability beyond which the correlations based on *u* increase while those based on *w* decrease; coherent turbulent structures become rather very flat, and pancake-like, beyond that critical point.

We also investigated how stability influences TKE and the different terms in its budget equation. Although LES does not yield the SGS TKE, we use a concise model proposed by Knaepen et al. (2002) that utilizes the − power-law spectral relationship in the inertial subrange to approximate it. As stability increases, the inertial subrange shrinks, with the lower wavenumber of the range shifting to a higher value. Consequently, the resolved TKE in LES decreases. With the SGS part of TKE approximately unchanged, the percentage that the SGS TKE contributes to the total TKE increases significantly. The effect of buoyancy dissipation generally increases while that of the viscous dissipation generally decreases. In very strong SABLs, it is shown that the buoyancy dissipation effects dominate over the viscous dissipation effects. The distribution of the total TKE destruction into viscous dissipation, resolved buoyant destruction, and SGS buoyant destruction is also compared with observational data analyzed in Bou-Zeid et al. (2010) and a generally good agreement is observed. This is an important validation of our LES code under the higher stabilities simulated here.

Finally, the applicability of Nieuwstadt's model under the whole range of stabilities investigated here is also examined. Nieuwstadt's model requires that the Richardson number approaches a constant across the SABL. However, it is shown that this basic assumption does not generally hold for all of the stability cases studied here. Rather, the Richardson number increases with height and stability, and no asymptotic value is observed. The Prandtl number, which characterizes the ratio of the transport efficiencies of momentum and heat fluxes, varies with height as well and generally increases with stability. Using turbulent diffusivities as an example, we found that the local-scaling assumption applies for different stabilities near the surface, but *z* should be replaced by a constant length scale for the *z*-less local scaling that we observe further aloft. This calls for a reexamination of previous analysis using *z*/Λ as a scaling parameter throughout the SABL. However, the power-law vertical profiles of turbulent stress and buoyancy flux predicted by Nieuwstadt are found to be generally acceptable across the stability range we examined.

## Acknowledgments

This work is supported by NSF Physical and Dynamic Meteorology Program under AGS-1026636 and by The Siebel Energy Challenge of Princeton University. The simulations were performed on the supercomputing clusters of the National Center for Atmospheric Research and those of Princeton University.

## REFERENCES

*Hydrology: An Introduction.*Cambridge University Press, 618 pp.

*GEWEX News,*No. 13, International GEWEX Project Office, Silver Spring, MD, 7–8.

*Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences,*H. H. Goldstine, Ed., Thomas J. Watson Research Center, 195–210.

*Turbulence and Diffusion in Stable Environments,*J. C. R. Hunt, Ed., Clarendon Press, 149–179.

*Advances in Geophysics,*Vol. 18, Academic Press, 224–236.

*Proc. 63rd Annual Meeting of the APS Division of Fluid Dynamics,*Vol. 55, Long Beach, CA, QD.00003.

*The Structure of Turbulent Shear Flow.*Cambridge University Press, 429 pp.