## Abstract

In this study, shear-based parameterizations of turbulent mixing in the stable atmospheric boundary layer (SABL) are proposed. A relevant length-scale estimate for the mixing length of the turbulent momentum field is constructed from the turbulent kinetic energy and the mean shear rate *S* as . Using observational data from two field campaigns—the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment and the 1999 Cooperative Atmosphere–Surface Exchange Study (CASES-99)— is shown to have a strong correlation with . The relationship between and corresponds to the ratio of the magnitude of the tangential components of the turbulent momentum flux tensor to , known as stress intensity ratio, . The field data clearly show that is linked to stability. The stress intensity ratio also depends on the flow energetics that can be assessed using a shear-production Reynolds number, , where *P* is shear production of turbulent kinetic energy and is the kinematic viscosity. This analysis shows that high mixing rates can indeed persist at strong stability. On this basis, shear-based parameterizations are proposed for the eddy diffusivity for momentum, , and eddy diffusivity for heat, , showing remarkable agreement with the exact quantities. Furthermore, a broader assessment of the proposed parameterizations is given through an a priori evaluation of large-eddy simulation (LES) data from the first GEWEX Atmospheric Boundary Layer Study (GABLS). The shear-based parameterizations outperform many existing models in predicting turbulent mixing in the SABL. The results of this study provide a framework for improved representation of the SABL in operational models.

## 1. Introduction

In geophysical flows, turbulence (or turbulent mixing) is a primary constituent in the transport of momentum and scalars (active or passive) influenced by velocity shear, stratification, and boundaries. In the atmospheric boundary layer (ABL), turbulent mixing directly impacts small-scale (e.g., air quality and wind energy) and large-scale (e.g., circulation and climate change) processes. In the nighttime as Earth’s surface cools, stable stratification suppresses turbulence as restorative buoyancy and gravitational forces limit vertical fluctuations and transfer turbulent kinetic energy (TKE) to potential energy (PE). The turbulent kinetic energy is a fundamental quantity that is used to define the characteristic velocity scale of the turbulent field (e.g., ).

To describe the evolution of TKE, a prognostic equation for stratified homogeneous conditions, with the Boussinesq approximation and Einstein summation convention can be written as

where is the turbulent kinetic energy, is the turbulent momentum flux, *g* is the gravitational acceleration, is a reference potential temperature, is the potential temperature (or buoyancy) flux, is the Kronecker delta function, and is the dissipation rate of turbulent kinetic energy. For a statistically stationary flow, Eq. (1) yields the equilibrium equation , where is the shear production of TKE and is the consumption of TKE due to buoyancy forces.

To quantify stratification, we rely on stability parameters such as the flux Richardson number written as

For stationary, homogeneous conditions, it is widely accepted that tends to an upper limit, , in the range of 0.20–0.25 for atmospheric flows (see, e.g., Mellor and Durbin 1975). The quantity may be negative for nonstationary flows represented through imbalance in the production and dissipation terms in Eq. (1). Another closely related stability parameter is the gradient Richardson number,

where is the buoyancy (or Brunt–Väisälä) frequency, is the averaged potential temperature, and is the mean shear rate. The quantity quantifies stability through the ratio of buoyancy frequency to shear leading to the qualification of a flow as shear dominated, , or buoyancy dominated, . A transition between these regimes is presumed to be associated with turbulence “collapse,” or relaminarization of the flow, signified by a critical value of stability such as the critical gradient Richardson number, (e.g., Richardson 1920). The exact value of has been the focus of a vast body of research (e.g., Taylor 1931; Miles 1961; Businger et al. 1971; Turner 1973; Rohr et al. 1988; Itsweire et al. 1993; Armenio and Sakar 2002; Zilitinkevich and Baklanov 2002; Galperin et al. 2007; Ohya et al. 2008; Grachev et al. 2013) with a consensus on the range and a classic value of 0.25 from linear stability analysis (Howard 1961; Miles 1961). An alternative to this fixed prescription of is a quasi-dynamic definition where, as stability increases, turbulence is sustained up to , but for decreasing stability, the flow remains laminar until (e.g., Woods 1969) whereby there is not a deliminator for turbulence collapse (Ansorge and Mellado 2014).

In contrast to the notion of turbulence collapse, numerous field campaigns have measured high mixing rates (e.g., strong turbulence) at strong stability (see, e.g., Kondo et al. 1978; Banta et al. 2002; Mahrt and Vickers 2006; Zilitinkevich et al. 2008), although the origin and sustenance of such events remains an area of intensive research. Strongly stable conditions are ideal for the formation and propagation of inertial phenomena such as low-level jets, meandering motions, and internal gravity waves wherein nonlinear interactions and decorrelation of layered motions can lead to influxes of shear-generated turbulence (e.g., Mahrt 1998). Increased shear from such events can lead to local values of remaining relatively constant in time despite strong gradients in potential temperature (or density), sometimes referred to as the Businger mechanism (Businger et al. 1971; Townsend 1976; Nieuwstadt 1984; Kim and Mahrt 1992; Derbyshire 1994; Mahrt 1998).

Numerical models are valuable tools for assessing and predicting atmospheric processes. While large-eddy simulation (LES) techniques have matured (see, e.g., Beare et al. 2006 and references therein), many large-scale operational models rely on turbulence parameterizations with the stable atmospheric boundary layer (SABL) represented rather poorly (e.g., Viterbo et al. 1999; Holtslag et al. 2013). For a Reynolds-averaged framework, the turbulent momentum and buoyancy fluxes for a horizontally homogeneous flow are parameterized through the flux–gradient relationships written as

where and are the averaged components of turbulent momentum flux, is the eddy diffusivity for momentum, and is the eddy diffusivity for heat.

From dimensional analysis, we can estimate from characteristic (relevant) velocity and length scales, , or length and time scales, . The specification of these scales relies on the sophistication of the prescribed closure scheme ranging from lower-order, algebraic formulations (e.g., similarity theory or *K* theory) to higher-order closures with transport equations for two or more turbulent parameters (e.g., – closure). Monin–Obukhov (MO) similarity theory (Monin and Obukhov 1954) is a prominent algebraic closure based on the nondimensional parameter, , where *L* is the Obukhov length given by

where is the surface friction (or shear) velocity, is the turbulent momentum flux, is the von Kármán constant, and is the vertical turbulent potential temperature flux at the surface (e.g., Obukhov 1946). MO similarity provides a direct link to the means gradients of velocity and potential temperature. To extend the conceptual framework of MO similarity theory beyond the surface layer (or lower ~10% of the SABL), Nieuwstadt (1984) developed a local similarity theory where the terms in Eq. (7) are defined locally. For weakly stable stratification under stationary, homogeneous conditions similarity theory performs very well (Kaimal and Finnigan 1994), but as stratification strengthens, turbulence statistics become independent beyond some height (or “*z*-less”) and traditional similarity theories begin to break down (see, e.g., Pahlow et al. 2001 and references therein). This breakdown warrants consideration of additional turbulence closure models.

One-equation closure continues to see usage in many operational models (see, e.g., Holt and Raman 1988; Cuxart et al. 2006; references therein). A common velocity scale in many models is given by ; thus, a prognostic equation for is solved [see Eq. (1)] and the mixing length is determined diagnostically. Prandtl’s mixing length hypothesis assumes the mixing length of the turbulent momentum field varies linearly with altitude, (Prandtl 1925). Extended surface-based mixing lengths have been proposed by Blackadar (1962) and Delage (1974) with many subsequent variations (see, e.g., Weng and Taylor 2003; Cuxart et al. 2006; Huang et al. 2013). Alternatively, the mixing length can be prescribed locally. A common local length-scale estimate is given by the buoyancy length scale, , where is the root mean square (rms) of the vertical velocity fluctuations (e.g., Brost and Wyngaard 1978; Nieuwstadt 1984; Hunt et al. 1985). Several additional local length scales have been investigated (e.g., André et al. 1978; Beljaars and Viterbo 1998; Weng and Taylor 2003) varying in accuracy and applicability within numerical models.

Recently, a renewed interest in model parameterizations of the SABL has led to new closure methods. Mauritsen et al. (2007) and Wilson (2012) developed one-equation closures with alternative velocity scales based on the total turbulent energy (TTE; ) and , respectively. Additionally, new developments in higher-order closures give explicit consideration to the existence of turbulence for high . Cheng et al. (2002) and Canuto et al. (2008) revisit the seminal – model of Mellor and Yamada (1982) providing a numerical framework that allows for turbulence to persist up to . Zilitinkevich et al. (2007) proposed the energy and flux budget (EFB) based on the transport of TTE without an explicit . While these recent closures represent significant advances in numerical modeling efforts of the SABL, we choose to approach this theoretical study maintaining as the velocity scale in developing local parameterizations of turbulent mixing for a broad range of stability.

The paper is organized as follows. In section 2, we briefly outline the datasets from field campaigns containing high-quality measurements of stably stratified atmospheric turbulence. We also introduce the LES data from the Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS). In section 3, we first present a discussion on relevant length scales and provide an assessment of the proposed length-scale estimate with the exact mixing length for momentum using observational data (e.g., Holtslag et al. 2013). This assignment reveals the shear-based length scale correlates remarkably well providing a basis for subsequent shear-based parameterization of and . Further evaluation of the proposed shear-based parameterizations is performed through an a priori analysis of LES data in section 4. Conclusions are given in section 5.

## 2. Description of datasets

### a. The SHEBA data

The Surface Heat Budget of the Arctic Ocean (SHEBA) experiment took place from October 1997 to September 1998 in the Beaufort Gyre north of Alaska drifting between approximately 74°N, 144°W and 81°N, 166°W. Full descriptions of the SHEBA campaign, measurement techniques, and data quality can be found in the works of Andreas et al. (1999), Persson et al. (2002), and Uttal et al. (2002). Primary mean and turbulent data were collected on a 20-m tower with five instrumented levels at 2.2, 3.2, 5.1, 8.9, and 18.9 m or 14 m depending on the season. The dataset contains hourly averaged (1 h) measurements. The gradients of mean velocity and potential temperature are calculated using a second-order polynomial fit and respective derivatives at the individual measurement levels (see, e.g., Grachev et al. 2005). The data from first level, 2.2 m, is excluded from our analysis because of prominent surface interactions with drifting snow leading to significant scatter (Sorbjan and Grachev 2010). Bin averaging is performed using to prescribe bins into which all other data are then sorted. The SHEBA data captured nearly 6000 h of measurements mostly for weakly stable conditions and absent of common contaminants (e.g., surface inhomogeneities, drainage flows).

### b. The CASES-99 data

The 1999 Cooperative Atmosphere-Surface Exchange Study (CASES-99) was a large collaborative field campaign lasting from 1 to 31 October 1999 in southeastern Kansas, located at approximately 38°N latitude (Poulos et al. 2002). Data were collected on a 60-m main tower with sonic anemometers at six levels and 34 thermocouples. Wind speeds were measured at 1, 5, 20, 30, 40, 45, and 50 m while temperatures were measured at 1, 5, 15, 25, 35, 45, and 55 m. The gradients of mean velocity and potential temperature are calculated using a sixth-order polynomial fit and respective derivatives at the individual measurement levels (Sorbjan and Grachev 2010). The original dataset of 5-min-averaged measurements is transformed to hourly averages (1 h) and bin-averaged based on following a similar procedure as for the SHEBA data. The CASES-99 dataset contains very stable measurements with weak and strong turbulent mixing events (Mahrt and Vickers 2006) as well as nocturnal phenomena such as internal gravity waves, Kelvin–Helmholtz instabilities, and drainage fronts (Poulos et al. 2002; Mahrt 2007).

### c. The GABLS LES data

The first GABLS initiative provided an intercomparison of LES (Beare et al. 2006) and single-column models (Cuxart et al. 2006) with an eye toward mixing-length specification for numerical weather prediction (NWP) and climate models (e.g., Holtslag et al. 2013). The LES study presents insights on the behavior of the SABL dynamics based on the simulations of Kosović and Curry (2000). We selected the 2.0-m-resolution LES reference dataset conducted by the National Center for Atmospheric Research (NCAR) for analysis. The LES is based on the subgrid-scale eddy viscosity model of Sullivan et al. (1994), which adheres to MO similarity theory in the surface layer region. For analysis, the turbulent quantities are determined from the summation of the resolved and SGS contributions (total = resolved + SGS).

## 3. Theoretical formulation

### a. Relevant length scales

From the flux–gradient relationships in Eqs. (4) and (5), the “exact” mixing length for the turbulent momentum field can be written as

In estimating , Prandtl’s mixing length implies a ubiquitous linear increase in eddy size with elevation. Blackadar (1962) observed that eddies tend to a maximum size and proposed the estimate , where is the maximum length, *G* is the geostrophic wind speed, and *f* is the Coriolis parameter. In practice, is typically assigned a constant value in the range 6–200 m (see, e.g., Huang et al. 2013 and references therein). Delage (1974) modified Blackadar’s proposition by adding a stability term such that , where is a model constant (Businger et al. 1971). Surface-based mixing lengths provide accurate model predictions when used in conjunction with model constants and stability functions (see, e.g., Cuxart et al. 2006) but lack generality in describing active scales of turbulent mixing.

In seeking a unified view of turbulent mixing in the SABL, pertinent length scales can be expressed locally for unforced (e.g., decaying), shear-dominated, or buoyancy-dominated turbulence. This approach is not meant to discount the contributions of MO or local similarity theories but rather to focus on fundamental aspects of stably stratified turbulence from which we can develop model parameterizations. For unforced or decaying turbulence, the pertinent length scale is given by (e.g., Pope 2000). While this regime is certainly important to a complete description, measurements in the SABL are difficult to obtain owing to wind speeds and flux magnitudes approaching instrument sensitivity. Since we are unable to quantify the unforced regime with atmospheric observations, we focus our efforts on the shear- and buoyancy-dominated regimes. In the shear-dominated regime, the Corrsin length scale represents the largest survivable eddy, (Corrsin 1958). Analogously, the Dougherty–Ozmidov length scale, , describes eddies in the buoyancy-dominated regime (Dougherty 1961; Ozmidov 1965). However, and more accurately describe small-scale turbulence in oceanic flows. In the atmosphere, direct measurements of are historically difficult to obtain (e.g., Sorbjan and Balsley 2008), which leads us to seek more suitable length scales from attainable quantities.

From the velocity scale , two additional length scales can be constructed from the shear and buoyancy time scales, and , respectively, written as

where corresponds to the average eddy size of active fluctuations for shear flow (Venayagamoorthy and Stretch 2010) and analogously for a buoyant flow (e.g., Deardorff 1976). We note that and have been previously studied in the context of numerical modeling application for LES (e.g., Kosović and Curry 2000) and single-column models (e.g., Grisogono and Belušić 2008). Using the SHEBA and CASES-99 datasets, we evaluate the correlation between the relevant scales, and , with the exact length scale, in Fig. 1. The quantity correlates quite well with with the data collapsing with a linear scaling of approximately 2 (e.g., ). On the other hand, a visual correlation between and is not as straightforward.

Continuing our assessment, we observe that the relationship between and is given by

Provided *S* and *N* are measured concurrently at the same resolution, and are interchangeable according to Eq. (11). In practice, vertical profile measurements (e.g., in the lower atmosphere or ocean) often yield higher resolution in potential temperature (or density ) corresponding to *N*, which favors over than in the velocity gradient measurements corresponding to *S*. Since measurements from fixed meteorological towers, such as the SHEBA and CASES-99 data, are taken at approximately the same intervals, we henceforth only use carrying forward in the definition of shear-based parameterizations without any loss in generality.

### b. Definition of constants

We now focus on defining the correlation, or model constant, between the estimate and the exact mixing length . The ratio of these two length scales is given by

where the ratio is simply the square root of the familiar stress intensity ratio, , given by

Thus, the shear-based mixing length estimate for is explicitly given by . To complete this assignment, we need a numerical value for *c* (or rather ).

Numerous experimental studies of unstratified turbulent wall-bounded flows dating back to the 1970s suggested a constant value of (see, e.g., Launder and Spalding 1974). Using measurements near neutral stability from the Cabauw experiment (e.g., Driedonks et al. 1978), Nieuwstadt (1984) arrived at a value of . Recent high-Reynolds-number (Re) direct numerical simulation (DNS) (e.g., Hoyas and Jiménez 2006) and experimental (e.g., Marušic and Perry 1995) studies suggest a possible constancy of for . While the value of has been extensively investigated for high- unstratified boundary layer flows, its dependency on stratification has not received much attention. Hence, a pertinent question is how does vary under stably stratified conditions (e.g., with )?

Mauritsen et al. (2007) studied this issue using an ensemble of SABL data and found a statistically significant decrease in with increasing stability as a function of [e.g., ], which we cast as

where denotes the value of the stress intensity ratio at neutral stability. We take based on unstratified high-Re DNS and experimental data (see, e.g., Karimpour and Venayagamoorthy 2014), whereas Mauritsen et al. (2007) selected from an empirical fit to the SABL data at . Figure 2 presents a comparison of as a function of for the 1-h bin-averaged SHEBA and CASES-99 data. The modified stability function of Mauritsen et al. (2007), as defined here in Eq. (14), follows the general trend in data approaching high . If we assume that high mixing rates are possible at high , it is likely that is not an inclusive function. Thus, the energetic state of the flow also provides feedback on the assignment of .

To quantify this energetic influence, we introduce a parameter that we will term the shear-production Reynolds number, , defined by

where we can define the eddy diffusivity for momentum as , is a turbulent (production) Reynolds number, is a nondimensional shear parameter, and is the turbulence production time scale. The quantity quantifies the relative influence of the shear to the turbulent time scale which we discuss later in this article (see section 3c).

In Fig. 3, we evaluate as a function of for the field data. The coexistence of high mixing rates (e.g., high ) at high in Fig. 3, while likely intermittent in nature, confirms that a ubiquitous turbulence shut-off mechanism in geophysical flows may be physically unrealizable. While weak- and strong-mixing regimes at high have been previously discussed in literature (e.g., Mahrt and Vickers 2006), adds some clarity on this issue providing a measure of the turbulent mixing generated by shear mechanisms. From this analysis we observe that has a functional dependence on stratification (e.g., in Fig. 2) well described by the modified function in Eq. (14). The flow energetics (e.g., in Fig. 3) also significantly influences , although a definitive function remains elusive to this point. Shah and Bou-Zeid (2014) found an interdependent influence of Richardson and Reynolds numbers in wall-bounded flows. What remains uncertain are the exact mechanisms linking stratification to energetics and ultimately turbulence in environmental flows where external factors (e.g., topography) are significant. Without a definitive theory, we employ the modified stability function in Eq. (14) to estimate the actual mixing length according to .

### c. Eddy diffusivity for momentum

In the Reynolds-averaged framework, the eddy diffusivity for momentum is a fundamental parameter quantifying turbulent mixing processes. For predominantly horizontal flow conditions, can be written as

Maintaining , we estimate the eddy diffusivity of momentum from . The relationship between and scales as

where the quantity is the stress intensity ratio. From our previous discussion on the stability function , the shear-based estimate for is given by

In Fig. 4 we evaluate the estimated and actual eddy diffusivity for momentum for the SHEBA and CASES-99 datasets. The quantity performs remarkably well in predicting . Strong turbulent mixing rates at high are indicative of large-scale events such as the breaking of highly energetic internal waves, low-level jets, or Kelvin–Helmholtz shear instabilities perhaps leading to the agreement with .

We can also evaluate the turbulent time scale associated with our estimate . Applying the equilibrium assumption to Eq. (1) results in . Assuming the dissipation rate for TKE can be given a ratio of TKE to a turbulent time scale, , we can rearrange by substituting Eq. (2) to give

Equation (19) effectively describes the influences of shear and stratification on the turbulent time scale, provided that shear is maintained in the flow (e.g., within the boundary layer).

### d. The turbulent Prandtl number

One possible avenue toward solving for the eddy diffusivity for heat is the link provided by the turbulent Prandtl number, . The turbulent Prandtl number for neutral stability, , has an approximate range of 0.5–1.0 under simple flow conditions (Kays and Crawford 1993; Pope 2000). For more complex flows, the exact specification of remains an area of contention, although literature widely suggests that . For example, Businger et al. (1971) proposed from measurements of the ABL made during the 1968 Kansas experiment. Likewise, Venayagamoorthy and Stretch (2010) arrived at a value of from analysis of homogeneous shear flow DNS data.

Considering stably stratified flow conditions relevant to this study, results from DNS studies (e.g., Gerz et al. 1989; Holt et al. 1992; Shih et al. 2000), laboratory experiments (e.g., Rohr et al. 1988; Strang and Fernando 2001), analytical models (e.g., Schumann and Gerz 1995; Zilitinkevich et al. 2007; Venayagamoorthy and Stretch 2010), and atmospheric observations (e.g., Webster 1964; Kondo et al. 1978; Kim and Mahrt 1992) show a strong connection between and stability. By definition, is assumed to vary linearly with and inversely with . Numerous studies have concluded that can be described as a monotonically increasing function of (see, e.g., Townsend 1958; Mellor and Yamada 1982; Nakanish 2001) reaching an asymptotic limit for large (Venayagamoorthy and Stretch 2010). Thus, for homogeneous, stably stratified sheared conditions, is a weak linear function at low stability, , increasing to at strong stability. What this simple theoretical derivation does not address is the fact that many SABL flows exhibit nonstationarity and inhomogeneities not included in most stability-dependent prescriptions of .

Venayagamoorthy and Stretch (2010, hereafter VS10) developed a novel proposition from the time scales of momentum and scalar mixing inclusive of nonstationary and irreversible mixing contributions given by

where , , and is the mixing efficiency nearing strong stability. In practice, the formula of VS10 does not differ significantly from other stability-dependent formulations that exhibit similar behavior for weak and strong stability (see, e.g., Kim and Mahrt 1992, Schumann and Gerz 1995, Zilitinkevich et al. 2007, and references therein). Figure 5 reveals that is indeed a function of from the mean bin-averaged SHEBA and CASES-99 datasets despite some scatter near-neutral stability for the SHEBA dataset. The data presented in Fig. 5 clearly illustrates that the increases with for strongly stable flow conditions.

### e. Eddy diffusivity for heat

We can predict the diffusivity of heat from . Initially, we maintain the exact formulation for [Eq. (16)] and evaluate a constant turbulent Prandtl number, , and the stability-dependent formula of VS10, . Figure 6 compares the estimated and actual eddy diffusivity of heat where noticeably improves the prediction of as compared to . We see significant improvement for data at (green and red symbols) consistent with our analysis of in Fig. 5.

Using the shear-based parameterization , we compare the estimate, , with for the SHEBA and CASES-99 data in Fig. 7. Overall, the shear-based parameterization in conjunction with predict quite well with the exception of weakly stable SHEBA data (blue diamonds) attributed to the lows values of nearing neutral stability (see Fig. 5).

## 4. A priori analysis of the SABL vertical profile

Field campaigns, such as SHEBA and CASES-99, provide measurements from fixed meteorological towers and are limited to the lower portion of the ABL. Measuring beyond this region leads to many difficulties in precise measurement techniques. LES data provides additional information on upper-SABL vertical structure, which we can use to further evaluate the proposed shear-based parameterizations.

### a. Revisit of model constants

The SHEBA and CASES-99 datasets clearly indicate that the stress intensity ratio is well described by the modified stability function in Eq. (14) (Fig. 2). We revisit as a function of and in Fig. 8 to include the 2.0-m-resolution NCAR LES dataset from the first GABLS intercomparison. In the LES analysis, the turbulent momentum flux is taken to be the sum of the resolved and SGS terms. We observe that the boundary layer in the LES is less energetic than the field data exhibiting a slightly sharper decay in with than the function . The also approaches the asymptotic upper limit of 0.25 for increasing but exhibits much smaller values of at low mixing rates. In contrast to the field observations, we do not observe the coexistence of high mixing rates (e.g., high ) and strong stratification in the LES data.

### b. Eddy diffusivity for momentum

To provide a more complete evaluation of our proposed parameterization, , we include existing zero-, one-, and two-equation closures. The local similarity model of Nieuwstadt (1984, hereafter N84) is given by

where is the boundary layer height using the approximation of Zilitinkevich (1972), where (Garratt 1982) and are empirical model constants. Additionally, the one-equation *k*– closure model of André et al. (1978, hereafter A78) can be written as

where is an empirical constant, is the mixing length with , and is an empirical model constant (Duynkerke and Driedonks 1987). This implementation represents the modifications proposed by Weng and Taylor (2003) to the original one-equation model of A78. Finally, the two-equation *k*–*ε* model (Launder and Spalding 1972) modified for stable stratification based on the equilibrium assumption is given by

where is an empirical model constant. It is noted that the empirical model constants and assume a stress intensity ratio of .

Figure 9a shows a comparison of the estimates with the actual eddy diffusivity for momentum as functions of nondimensionalized height , where m is the maximum height of the computational domain. While the estimate of N84 preforms remarkably well, it relies strongly on an accurate prediction of the boundary layer depth *h*. A78 misses the shape and magnitude of and the –*ε* model significantly overestimates . This result is a noted difficulty with the –*ε* model in model applications to the SABL (e.g., Detering and Etling 1985; Apsley and Castro 1997). Overall, estimates remarkably well and does not require any foreknowledge of the boundary layer. These results present a potentially very significant advantage under very stable conditions where the layers in the SABL tend to decorrelate and the applicability of similarity theory tends to break down.

### c. Eddy diffusivity for heat

In section 3, we observed that the proposition of VS10 fits the general trend of a linear increase with (Fig. 5) and accurately predicts the eddy diffusivity for heat (Figs. 6 and 7). On this basis, we use only to estimate the eddy diffusivity for heat using the aforementioned models for in addition to the proposed shear-based parameterization. We compare the estimated and actual eddy diffusivity for heat for the LES data in Fig. 9b. The comparison in Fig. 9b reveals the same general trend for predicted eddy diffusivity for momentum and heat because of the consistent use of . Once again the model of N84 and shear-based parameterization, , predict closely.

## 5. Conclusions

In this research, we sought to characterize the active turbulent eddies with pertinent length scales using local (or *z*-less) formulations. Using datasets from the SHEBA and CASES-99 field campaigns, two length-scale estimates, and , are evaluated for efficacy in predicting the mixing length of the turbulent momentum field, . The quantity is shown to be closely correlated with from near-neutral to strongly stable conditions as measured by . Also, given that and are interchangeable through a factor of and data are measured at the same resolution, we have developed a shear-based parameterizations based on the length-scale .

We observed that the ratio of the shear-based parameters for turbulent mixing and exact Reynolds-averaged quantities corresponds to the stress intensity ratio, . We confirmed that depends both on stability and flow energetics. A modified form of the stability function proposed by Mauritsen et al. (2007) with at neutral stability correlates well with the SHEBA and CASES-99 datasets. Using as a measure of flow energetics, we observed that the stress intensity ratio exhibits an upper limit, , at high for low and moderate values of but decreases for low . Further studies are needed before any definitive conclusions can be drawn on the functional dependence of with owing in large part to the complex interactions between stability and flow energetics in geophysical flows.

From this assessment of length scales and model constants, we proposed a shear-based parameterization for the eddy diffusivity for momentum, , which provides efficacy in predicting the actual . Notably, the CASES-99 dataset includes observations of nocturnal phenomena (e.g., low-level jets and Kelvin–Helmholtz shear instabilities) that may inject shear-generated three-dimensional turbulence. Such events may lead to estimating quite well even in the strongly stable regime. We further evaluated from atmospheric observations and the formulation of VS10 yielding a more accurate prediction of as opposed to a constant value for especially for moderately and strongly stable flow conditions.

These turbulence parameterizations are further evaluated with an a priori analysis of the GABLS LES intercomparison data for the SABL vertical structure. While the LES data exhibit only moderately stable stratification, it is worth noting that the shear-based parameterization predicts very well, both in magnitude and shape over the boundary layer depth. Using with the stability-dependent of VS10, , provides an accurate estimation of the eddy diffusivity for heat, .

Continuing the trajectory of this current work, a natural extension is to implement the proposed turbulence parameters in an operational model for the stable atmospheric boundary layer within the larger context of turbulence modeling in NWP, global circulation, and climate models (e.g., Holtslag et al. 2013).

## Acknowledgments

We thank the three referees for their helpful comments and recommendations. We would like to acknowledge financial support of the National Science Foundation under Grant OCE-1151838. The authors also gratefully acknowledge the SHEBA Atmospheric Surface Flux Group and the individuals and institutions of the CASES-99 initiative for their skillful collection and analysis of the valuable SABL datasets. We also thank the organizers of the GABLS LES Intercomparison and Dr. Peter Sullivan for providing access to the detailed postprocessed data.

## REFERENCES

*13th Symp. on Boundary Layers and Turbulence,*Dallas, TX, Amer. Meteor. Soc., 550–555.

*Clear and Cloudy Boundary Layers,*A. A. M. Holtslag and P. G. Duynkerke, Eds., Academy of Arts and Sciences, 287–304.

*Seminars on the Treatment of the Boundary Layer in Numerical Weather Prediction,*European Centre for Medium Range Weather Forecasts, 234–284.

*Fourth Symp. on Meteorological Observations and Instrumentation,*Denver, CO, Amer. Meteor. Soc.,

**18,**011702, doi:.

**83,**

**153,**355–387, doi:.

*Proc. Roy. Soc. London,*

**132A,**499–523, doi:.