## Abstract

A spectral budget model is developed to describe the scaling behavior of the longitudinal turbulent velocity variance with the stability parameter and the normalized height in an idealized stably stratified atmospheric surface layer (ASL), where *z* is the height from the surface, *L* is the Obukhov length, and *δ* is the boundary layer height. The proposed framework employs Kolmogorov’s hypothesis for describing the shape of the longitudinal velocity spectra in the inertial subrange, Heisenberg’s eddy viscosity as a closure for the pressure redistribution and turbulent transfer terms, and the Monin–Obukhov similarity theory (MOST) scaling for linking the mean longitudinal velocity and temperature profiles to *ζ*. At a given friction velocity , reduces with increasing *ζ* as expected. The model is consistent with the disputed *z*-less stratification when the stability correction function for momentum increases with increasing *ζ* linearly or as a power law with the exponent exceeding unity. For the Businger–Dyer stability correction function for momentum, which varies linearly with *ζ*, the limit of the *z*-less onset is . The proposed framework explains why does not follow MOST scaling even when the mean velocity and temperature profiles may follow MOST in the ASL. It also explains how *δ* ceases to be a scaling variable in more strongly stable (although well-developed turbulent) ranges.

## 1. Introduction

The scaling laws of the root-mean-squared longitudinal turbulent velocity component normalized by the friction velocity with distance *z* from a solid boundary has received renewed attention in meteorology and turbulence (Alfredsson et al. 2011; Smits and Marusic 2013; Poggi et al. 2006; Cai et al. 2008; Hsieh and Katul 2009; Hansen et al. 2012; McKeon 2013; Stevens et al. 2014; Yang et al. 2014). While it has been known for some time that in the atmospheric surface layer (ASL) deviates from the expectations set by the Monin and Obukhov similarity theory (MOST) even for near-neutral conditions (Bradshaw 1978; Kader and Yaglom 1990; Kaimal 1978; Kaimal and Finnigan 1994; Charuchittipan and Wilson 2009), MOST scaling for remains the “workhorse” model for practical problems because alternative scaling laws that retain the simplicity of MOST remain lacking (Zilitinkevich and Calanca 2000). The goal of this work is to offer a potential alternative to MOST for the stable ASL for “idealized” conditions. These idealized conditions include stationary planar homogeneous flows in the absence of subsidence and well-developed turbulence so as to ensure an extensive inertial subrange spectrally separating production from viscous dissipation scales.

A universal logarithmic scaling with *z*, originally proposed by Townsend (1976) and expressed as for neutral turbulent boundary layers, has been recently confirmed by different high Reynolds number laboratory experiments and large-eddy simulations (Marusic et al. 2013; McKeon 2013; Smits and Marusic 2013; Stevens et al. 2014), where the atmospheric boundary layer (ABL) height *δ* is featured as a significant dynamical variable. However, a satisfactory phenomenological theory predicting how is altered by thermal stratification in idealized stably stratified ASL flows while recovering its well established near-neutral limit remains elusive and frames the scope here. By phenomenological theory, we mean a theory that relates a number of empirical observations on in a way that is consistent with the fundamental equations but the theory itself is not directly derived from these fundamental equations. Hence, as a starting point for developing such a phenomenological theory, a number of observations and findings must be reviewed.

First, under near-neutral conditions and for an idealized ASL, the mean velocity *U* exhibits a logarithmic scaling with *z* roughly for the same range of as . The logarithmic scaling for the mean velocity is expressed as , where *κ* is the von Kármán constant and is a surface roughness coefficient, and and are nondimensional velocity and length from the wall, respectively, and ν is the kinematic viscosity. Second, in the presence of thermal stratification, the mean velocity is reasonably described by incorporation of a stability correction function for momentum that is only a function of the stability parameter according to MOST, where *L* is the Obukhov length and indicates unstable stratification and indicates stable stratification. Links between the energy spectrum of turbulence and (or similar measures) have received renewed interest in the last decade (Sukoriansky et al. 2006; Katul et al. 2011; Li et al. 2012b; Salesky et al. 2013; Sukoriansky and Galperin 2013; Katul et al. 2013, 2014; Li et al. 2015b) but will not be elaborated on here. Third, several studies reported measurements of spectra, cospectra of momentum and heat fluxes as well as the corresponding integral length scales, and their variations with the stability (Kaimal et al. 1972; Kaimal 1973; Caughey 1977; Högström 1990; Launiainen 1995; Canuto et al. 2008; Basu et al. 2014). This dependence of integral scales on *ζ* will be explicitly accommodated in the proposed spectral budget model. Fourth, a balance between the turbulent kinetic energy dissipation rate, mechanical production, and buoyant destruction (Wyngaard and Coté 1971; Katul et al. 2014) is often assumed though several datasets suggest the existence of a finite transport term that can vary with *ζ* as reviewed elsewhere (Salesky et al. 2013). Implications and biases due to a stability-dependent nonnegligible transport term will be discussed later in an appendix but not explicitly included at this stage. Last, it was noted by Wyngaard and Coté (1971), Nieuwstadt (1984), and Olesen et al. (1984) among others that all dimensionless parameters resulting from local combinations of flow variables reach a uniform value at sufficiently high *ζ* (but still maintaining a sufficiently developed turbulent state), thereby losing their *z* dependence. This phenomenon is referred to as *z*-less stratification, and its limit is generally assumed to be close to 1 or 2 (Olesen et al. 1984). However, measurements of , compiled by Pahlow et al. (2001), suggest that a power-law scaling in *ζ* appears to be maintained even when , which was used to argue against the *z*-less scaling for . The analysis of Pahlow et al. (2001) was problematic because of the presence of the so-called submesomotions. The revised analysis (Basu et al. 2006) of the same dataset (augmented with others) did not show such spurious behavior. The proposed spectral model will be used to partly clarify some (but not all) of the conditions promoting *z*-less stratification, including their connection to the shape of .

The proposed framework builds on a previous spectral budget method developed to explain the logarithmic scaling exhibited by when (Banerjee and Katul 2013). A later study expanded the spectral budget method to describe under unstable conditions (i.e., ) by the addition and appropriate modeling of a buoyant production term in the spectral budget (Banerjee et al. 2014). It has been found that *z*, *δ*, and *L* are all parameters required for describing the variation of . Analytical results for from the theory were also possible in two limiting conditions identified as and (Banerjee et al. 2014). While progress has been made toward understanding several characteristics of in the unstable boundary layer, there has been limited progress for the stable atmospheric boundary layer (SABL), which is the main focus here. This lack of progress is partly experimental because contamination from nonturbulent phenomena such as gravity waves; drainage forces from change in terrain slope, low-level jets, and mesoscale motions; and passage of clouds during nocturnal conditions, which are known to impact , are ubiquitous in field experiments under stable conditions (Pahlow et al. 2001; Cava et al. 2004). Thus, the SABL is characterized by different levels of balance between mechanical production of turbulent kinetic energy (TKE) and buoyant damping under different conditions ranging from almost nonturbulent to well-mixed conditions (Ohya et al. 1997). In low wavenumbers, the effects of gravity waves dominate close to the Brunt–Väisälä frequency *N* and can influence the turbulent spectra and cospectra of the velocity component (Einaudi and Finnigan 1981; Finnigan and Einaudi 1981; Leyi and Panofsky 1983; Finnigan et al. 1984; Hunt et al. 1985; Rohr et al. 1988; Edwards and Mobbs 1997; Mahrt et al. 2001; Cava et al. 2004; Campos et al. 2009; Mahrt et al. 2012). Heat and pollutants can be transferred vertically by waves and low-frequency turbulent motions and a vertical diffusion length scale of for heat and pollutant transfer has also been identified where is the vertical velocity variance (Hunt et al. 1985). Weinstock (1981) reported a strong correlation between the turbulent kinetic energy dissipation rate and . The proposed spectral budget here will not consider such nonstationarities, the role of gravity waves (and hence *N*), and intermittent switching between turbulent and nonturbulent states. The turbulence in the SABL is assumed to be sufficiently developed so that a sufficiently large separation between the integral length scale and the viscous dissipation length scale (or the Kolmogorov microscale) always exists and Kolmogorov’s theory for the inertial subrange (Kolmogorov 1941) approximates velocity and temperature spectra for all scales bounded between the integral and Kolmogorov scales. This is consistent with the observation by Huang and Bou-Zeid (2013) that this “continuously turbulent SABL” is different from the “intermittently turbulent SABL.” The continuously turbulent SABL is also described as “weakly stable” and is the subject of the present work while the intermittently turbulent SABL is often designated as a“very stable” boundary layer (Mahrt 1998; Huang and Bou-Zeid 2013). Approaches that analyzed similar problems where the spectrum of (vertical) velocity and temperature shift from waves at large scales (characterized by a −3 power law) to inertial at smaller scales (characterized by a −5/3 power law) are also not within the scope here and are reviewed elsewhere (Sukoriansky et al. 2006; Sukoriansky and Galperin 2013). The proposed derivation also does not require or assume any link between the turbulent kinetic energy dissipation rate and . However, the derivation here does assume that much of the turbulent kinetic energy resides in the longitudinal and lateral velocity components, which is plausible given that the damping effects of atmospheric stability are in the vertical direction.

## 2. Theory

### a. Background and definitions

A coordinate system is defined with *x*, *y*, and *z* forming the longitudinal, lateral, and vertical directions, respectively. Standard Reynolds decomposition notation is also employed so that the longitudinal *u*, lateral *υ*, and vertical *w* velocities are decomposed into a stationary mean and fluctuating or turbulent quantities using , , and with based on the choice of the coordinate system. Potential temperature and pressure are also decomposed into a mean (*T* and *P*) and turbulent fluctuations ( and ), respectively. As earlier noted, it is assumed that the flow is characterized by high Reynolds and Peclet numbers (negligible molecular viscosity and diffusivity compared to their turbulent counterpart), stationary [], planar homogeneous in the mean [] and without subsidence () or significant mean horizontal pressure gradient () as common to MOST. When these assumptions are incorporated into the equations for mean horizontal velocity and mean temperature, they result in and , where and are the momentum and sensible heat fluxes, respectively. Beyond ignoring molecular diffusion and viscosity, the high Reynolds and Peclet numbers are assumed to be sufficiently large so as to further suppress formation of gravity waves or laminarization of the flow. A large spectral separation is also assumed between the integral length scale of the flow and the Kolmogorov microscale so that scales bounded by these two limits maintain an approximate locally homogeneous and isotropic state as earlier noted.

### b. The TKE budget

The TKE budget for such an idealized ASL can be written as

where is the TKE; , , and ; is the buoyancy parameter; *g* is the gravitational acceleration; *ρ* is the mean air density; *t* denotes time and is the mean TKE dissipation rate; and *T* denotes potential temperature. The first, second, and third terms on the right-hand side (RHS) of Eq. (1) are the mechanical production and buoyant destruction of TKE (written as positive as per convention; note that would be negative for a stable atmosphere) and the transport of TKE by turbulence and pressure–velocity interactions, respectively. For several ASL applications, the production, buoyant, and dissipation terms can exist in a near balance (Townsend 1976; Pope 2000; Wilson 2008). However, some datasets demonstrated that the transport term cannot be neglected, especially in stable stratification, and the imbalance in the TKE budget is suggested to vary with the stability parameter (Salesky et al. 2013). A consequence of this imbalance on the scaling will be explored separately in appendix A but this imbalance is momentarily ignored. In case of a TKE balance (i.e., negligible transport), MOST (Monin and Obukhov 1954) yields and , where and is the stability correction function for momentum as mentioned earlier. In case the transport terms are significant, can be modified as following Salesky et al. (2013), where is a function of *ζ* that may be determined from experimental datasets.

Other common quantities are now defined. The flux Richardson number is given by the ratio of the buoyant destruction rate of TKE to the mechanical production rate of TKE :

and the gradient Richardson number is given as

where , is the mean velocity gradient, and is the mean potential temperature gradient. The is, as before, the friction velocity defined as and is defined as . Using these definitions, the flux Richardson number can be related to *ζ* using

Unless otherwise mentioned, can be described by the widely used empirical Businger–Dyer (Businger and Yaglom 1971; Dyer 1974) form: for , which is the range of stability conditions explored here. This is commensurate with the critical flux Richardson number (, same for the gradient Richardson number as well) for the applicability of local similarity theory since when , (Grachev et al. 2013). The concomitant temperature stability correction function is also assumed to be given by its MOST representation as . To a leading order, , though for , this approximation no longer holds (Zilitinkevich et al. 2013; Katul et al. 2014). With this background, the spectral budget for TKE and its connection to are presented.

## 3. A spectral budget formulation for TKE

At an arbitrary wavenumber *k* defined along the longitudinal direction, the spectral budget can be expressed as (Hinze 1959; Panchev 1971; Banerjee et al. 2014)

where the first, second, third, and fourth terms on the RHS represent the mechanical production and buoyant destruction of TKE in the range of , the transfer of TKE in the range , and the viscous dissipation in the range of , respectively. It is to be noted that the formulation here assumes that . The mechanical generation term is also positive and as buoyancy acts against mechanical generation of turbulence, the term would be negative for the SABL. The transfer and viscous dissipation terms are defined with a positive sign. If *k* tends to 0 without the presence of any buoyancy term, and . The sign of is negative so the production term is always positive here. Similarly, the normalizing property owing to in the SABL. Using the Heisenberg model, for all eddies between the wavenumbers 0 and *k*, the action of smaller eddies can be represented by an additional viscosity, since the way in which these smaller eddies transfer momentum is similar to the action of ordinary friction. This additional viscosity must depend upon the intensity of the small eddies—on the part of the spectrum with large wavenumbers by means of a certain integral that was formulated based on dimensional arguments by Heisenberg (1948). Using this argument, the transfer term can be “closed” as

where is a the wavenumber-dependent Heisenberg viscosity. With this simplification, the spectral budget Eq. (5) becomes (Banerjee et al. 2014)

### a. Modeling the mechanical production term

To link the spectral budget to MOST scaling, it is assumed that within the ASL, where varies with *ζ* and is specified later. The represents a stability-dependent crossover between the large scales and the inertial scales, which may be related to the integral length scale of the flow but is not necessarily identical to it. The reduction of *γ* with stability from a neutral limit magnitude of unity implies that the extent of the inertial subrange reduces with increasing *ζ*.

Thus for a wavenumber *k* that satisfies , the traditional inertial subrange scaling applies (Kolmogorov 1941) so that

where is the Kolmogorov constant when *k* is interpreted as a one-dimensional cut along the longitudinal direction and is the Kolmogorov constant associated with three-dimensional wavenumbers. The coefficient originates from a 1D interpretation (along *x*) of the TKE spectrum that is locally homogeneous and isotropic given by , with since the Kolmogorov spectral constants for the individual one-dimensional velocity spectra are , , and , respectively (Banerjee and Katul 2013). In conventional turbulence literature where index notation is prevalent, is equivalent to the as discussed by Pope (2000), where is the wavenumber along the longitudinal direction and , , and are component velocity spectra along the longitudinal, lateral, and vertical directions, respectively. Throughout, meteorological notation is employed unless otherwise stated, and *k* is interpreted as as noted before for notational simplicity.

Using the cospectral budget formulation provided by Katul et al. (2014), the following form can be derived for :

which is in agreement with the accepted “−7/3” scaling for the cospectra of momentum flux in the inertial subrange (Lumley 1967), where , and is the Rotta constant. Hence, the production term in the spectral budget [Eq. (7)] can be simplified by substituting in Eq. (9) with , as well as setting to yield

### b. Modeling the buoyant destruction term

Using the cospectral budget formulation for provided in Katul et al. (2014), the following form can be derived for :

where , is the Rotta constant for heat assumed to be the same as its momentum counterpart (equal to ) here, and varies with *ζ* as follows:

where is the Kolmogorov–Obukov–Corrsin constant for the temperature spectrum in the inertial subrange (along *x*), is a constant associated with the isotropization of the production term in the cospectral budget of temperature flux and can be predicted from rapid distortion theory (RDT) for isotropic turbulence, and is the thermal variance dissipation rate defined as

Hence, *Q* can be simplified as

It is interesting to note that is consistent with the findings from the Kansas experiment (Kaimal and Finnigan 1994). Because , the buoyant destruction term in the spectral budget [Eq. (7)] can be simplified as

### c. Spectral budget at

The spectral budget equation can now be used to find the spectral form of (which is later linked to ) in the range using the modeled forms of the mechanical production and buoyant destruction terms. It is necessary to recast the spectral budget equation [Eq. (7)] in the following form:

where it is assumed that , that is, the turbulent viscosity is much higher than the molecular viscosity for high Reynolds number flows. Using Eqs. (10) and (15), Eq. (16) is rewritten as

Assuming is inertial for *k* bounded by given by Kolmogorov’s theory (i.e., ignoring exponential cutoff due viscous effects at high *k*), can be determined as (Banerjee et al. 2014)

where is the Heisenberg constant defined as . Substituting the definition of from Eq. (18), Eq. (17) can be written as

which is consistent with Banerjee et al. (2014):

where here is defined as

In the neutral limit, , , , and , thus recovering to be

With appropriate values of the constants [, , , and ], yields a .

The function is now discussed. As earlier noted, the integral length scale decreases with increasing *ζ* (Kaimal et al. 1972; Kaimal 1973; Caughey 1977; Högström 1990; Launiainen 1995; Canuto et al. 2008). While is not identical to the integral length scale as earlier noted, its dependence on *ζ* is assumed to follow a similar pattern. For this purpose, so as to ensure a , where is a decay coefficient. Possible decay rates of *γ* with increasing *ζ* are discussed and are shown in Fig. 1.

The variations of the integral length scale of the longitudinal velocity normalized by *z* are shown as black circles on Fig. 1, which were computed by Salesky et al. (2013) using the Advection Horizontal Array Turbulence Study (AHATS) data collected in the ASL. A good fit to this data is . The dashed–dotted line refers to the integral length scale of the vertical velocity whose functional form has been computed from the Kansas experiment as (Kaimal and Finnigan 1994). It is to be noted that its neutral limit is greater than 1, which indicates a larger length scale at lower *k* than z. However, when normalized by this neutral limit, it can also be used as a guide to how the breakpoint wavenumber in the spectrum representing the crossover from large to inertial subrange scales shifts to higher wavenumbers with increasing stability. In addition, because the integral length scale of the vertical velocity is closely aligned with this breakpoint, it may be better suited for approximating than the integral length scale of the longitudinal velocity. As evident from Fig. 1, its variation is quite similar to the variation of the integral length scale of the longitudinal velocity. Interestingly, is also shown for reference and appears to be also similar to the one inferred from the integral length scale of the vertical velocity. For this reason, this aforementioned form is used for analytical tractability. Notwithstanding which precise approximation to use, it appears that a value of in the range from 3 to 4.7 is reasonable.

### d. Formulation for C_{st}: Numerical fitting

To proceed with an analytical derivation, the integral expressed in Eq. (20) must be simplified. Figure 2 shows the variation of the function with *ζ*. As observed, suggests the possible existence of a two-regime formulation somewhat similar to the result of Banerjee et al. (2014). In zone I, is roughly flat until , assuming the neutral limit of , which is designated as . In zone II, can be described by an approximate power law of the form , where *a* and *b* can be obtained by a numerical curve fitting as and with a coefficient of determination of 0.99. Although a smooth transition between the two regimes (zone I and zone II) can be possibly formulated by an interpolation scheme, it is assumed that is the location of this transition for all practical purposes. With this two-part formulation for , the shape of the can now be inferred for each regime in the limit .

### e. Formulation for zone I

Since the form of in the limit is unknown, it may be assumed to be a power law; that is,

Substituting this power-law expression in Eq. (20), the following identity can be derived:

Equating terms in the left-hand side (LHS) and RHS of Eq. (24), one can obtain

while the parameter is known a priori as discussed earlier.

Thus to summarize, for zone I,

where . Thus the spectrum recovers the much-discussed scaling as reviewed elsewhere (Banerjee et al. 2014). Furthermore, it can also be assumed that the scaling is valid until , where *H* is a measure of the size of largest turbulent eddies in the SABL and the spectrum is flat at lower wavenumbers; that is, form continuity requirement at .

### f. Formulation for zone II

Following the same approach presented in the previous subsection, is assumed to be a power law; that is,

Substituting this general power-law expression in Eq. (20), the following identity can be derived:

Equating terms in the LHS and RHS of Eq. (28), one can obtain

where *a* and *b* are known a priori from . It is interesting to note that , which indicates that the spectrum follows a scaling law analogous to its inertial value at higher *ζ*. Because of the similarity in the spectral exponents and to ensure spectral continuity, we simply extend the inertial subrange spectrum up to some scale , where is discussed later. However, it must be emphasized that the main mechanisms leading to the approximate −5/3 scaling here differ from those in the inertial subrange. For maintaining maximum simplicity here, the spectrum can be described as

Now the determination of is discussed. It is known that the turbulent (but not necessarily the total) kinetic energy in is quite suppressed and eddies are appreciably deformed by stratification at eddy sizes much larger than the Ozmidov length scale , where . At these large scales, nonturbulent phenomena such as gravity waves become energetically prevalent but the turbulent contribution to the energy is small (or negligible). For the idealized ASL considered here where , . To a first order, the declines at a rate commensurate with when setting and . An assumption of maximum simplicity again is that the turbulent component of the kinetic energy spectrum lacks significant energy at with declining with *ζ* at a rate commensurate with . However, instead of determining from , it will be determined so as to ensure continuity between zones I and II when deriving , as will be discussed later. Another interesting observation is that the spectrum exhibits an apparent inertial signature in its exponent (=−5/3) at higher stability. Given that the spectrum is also found to follow −5/3 scaling at highly unstable conditions (Banerjee et al. 2014), it can be stated that the signature of buoyancy is an apparent “inertialization” or transition to an approximate −5/3 scaling from a combination of −5/3 and −1 scalings dominating the spectrum for near-neutral conditions.

### g. Linking the longitudinal velocity variance to TKE

For neutral and unstable atmospheric surface layers, it was argued elsewhere (Basu et al. 2010; Banerjee and Katul 2013; Banerjee et al. 2014) that and thus . It is assumed that this condition is also valid for stable conditions, especially at high Reynolds and Peclet numbers when turbulence is well separated from gravity waves. To test this assumption, experimental data (described in detail in section 4) collected over two distinct vegetation canopies at the Duke forest over 4 yr are used. Figure 3 demonstrates a one-to-one comparison between and *e* for two different seasons (i.e., summer and winter) and for two different vegetation canopies: a second growth oak–hickory hardwood (HW) forest and a loblolly (PI) plantation. In all cases it can be assumed that irrespective of surface cover and stability parameters (the slopes are 0.80, 0.84, 0.70, and 0.76 for winter pine, summer pine, winter HW, and summer HW, respectively). Thus the TKE spectrum can be integrated to obtain for the purposes here.

### h. The longitudinal velocity variance

Integrating the spectra from the formulation for zone I [Eq. (26)] to obtain the longitudinal velocity variance results in

which can also be written in a form consistent with Townsend’s attached eddy hypothesis as

where

and

where and thus the Townsend parameters now vary (mildly) with stability for near-neutral to slightly stable conditions.

If the second formulation for zone II is used, the can be integrated to derive the longitudinal velocity variance as

It is also interesting to note that the boundary layer height is not present in the formulation for zone II, which indicates that at higher stability is no longer sensitive to the boundary layer height (as expected). The parameter can be found out by imposing continuity condition on at (i.e., the transition between the two formulations). Equating Eqs. (32) and (36) at , can be calculated to have a near-neutral limit and can be assumed to reduce with stability somewhat similar to . Again for simplicity, we choose based on the scaling argument that the Ozmidov length scale decreases with increasing stability as . Unfortunately, this continuity constraint results in weak scale separation between all the key length scales here—, *z*, *L*, and *δ* (at least when compared to zone I). Given the number of assumptions used to separate the two formulations and the associated uncertainties, predictions based on the formulation for zone I and zone II are separately compared to data for all *ζ* and biases are then analyzed as will be shown in section 5 and appendix B. Both formulations predict the correct neutral limit of (between 2 and 3).

The onset of *z*-less scaling in can now be assessed by inserting the estimate of into Eq. (36) and rearranging to yield

Evidently, for , and *z*-less scaling is assured for the formulation derived in zone II. However, the onset of *z*-less scaling for stability regimes where is of interest here. As *ζ* increases, *z*-less scaling in appears to be recovered from Eq. (37) if increases at rate much faster than *ζ* so that . Conversely, absence of *z*-less scaling can arise when increases at rate slower than *ζ*. It can be surmised that *z*-less scaling in is dependent on whether becomes small compared to unity or simply independent of *ζ* with progressively increasing *ζ*. For the Kansas experiments, for the SABL so that , a constant at large *ζ*, and is suggestive of possible *z*-less scaling in .

## 4. Datasets

To evaluate the model, four datasets are used. One dataset is a published laboratory experiment (Ohya et al. 1997) where the data sources are published figures digitized by us, one is an experiment over a large lake, and two datasets are from long-term monitoring initiatives above a pine (PI) plantation and a HW forest. The forest experiments are separated into summer and winter seasons to reflect large and small leaf area index (LAI) values.

The laboratory experiment was conducted by Ohya et al. (1997) in a wind tunnel where stably stratified flows were generated by cooling the aluminum floor to 3°C and maintaining the ambient air at 50°C. The experimental runs covered a range of from 0.013 to 0.13 m s^{−1}, varying over an order of magnitude. The boundary layer height was measured and a range was reported in their results. Further details on the experiments and datasets can be found in Ohya et al. (1997).

The lake data include 20-Hz eddy covariance measurements of three-dimensional velocity and air temperature at four different heights (1.65, 2.30, 2.95, and 3.60 m) above an extensive lake surface. Details about the lake data and quality control measures can be found elsewhere (Vercauteren et al. 2008; Huwald et al. 2009; Li and Bou-Zeid 2011; Li et al. 2012a). In particular, the averaging interval is 30 min and data with fluxes measured at the four heights differing by more than 10% are excluded so as to maintain the constant flux assumption required by the derivation here and MOST scaling for mean velocity and temperature (Li et al. 2015a). The stability range covered in the lake measurements are mildly stable conditions (). As a result, although the boundary layer heights are not measured in the experiments, they are expected to be not very low (on the order of hundreds of meters). Figure 4 shows the variations of with *ζ* in the lake data. It is clear that the decreases with increasing stability, which is expected. Interestingly, the decrease of is captured by with (the black line) and (the red line), as shown in Fig. 4. Note that indicates the neutral m s^{−1}. This approach of multiplying the neutral with can also be used to capture the variation of for comparison with the laboratory data as well in absence of explicit for each run.

The forest experiments were conducted over two different but adjacent vegetation canopies: a loblolly pine plantation and a 100-yr-old second-growth oak hickory hardwood forest. The data include 10-Hz eddy covariance measurements of three-dimensional velocity and air temperature and the details of the experiments can be found elsewhere (Juang et al. 2008; Katul et al. 2012). The data span over 4 yr and are divided into summer and winter seasons. Since the experiments were conducted over the canopy, the is modified as , where *d* is the zero-plane displacement height. This *d* is set to be zero in winter for the HW canopy and set to be otherwise for both forests, where is the canopy height ( m for HW and m for PI). The height of measurement *z* is 39.5 m for the HW and 15.5 m for the PI sites. The PI forest had a LAI varying from 2.65 to 4.56 m^{2} m^{−2} and the HW forest had a maximum LAI about 6 m^{2} m^{−2}. The stability covered by this dataset ranges from mildly stable to highly stable conditions (). The also varies with stability with modal values of 0.06 m s^{−1} in the summer and 0.35 m s^{−1} in the winter for both hardwood and pine canopies. The two forests were also collocated so they were subjected to similar meteorological conditions. There are 11 449 and 14 977 30-min runs for the PI site in the summer and winter, respectively. For the HW, there are 13 049 and 15 045 30-min runs in the summer and winter, respectively. To assess the applicability of MOST in the canopy sub layer (CSL), which is roughly 2–3 times , some of the stored high-frequency data above the pine forest are used for illustration. The second-order structure function for *u* is calculated independently, which is regressed with where (Stull 1988) and *r* is the lag in meters in the longitudinal direction (computed by multiplying the sampling time interval with the mean velocity using Taylor’s frozen turbulence hypothesis). Within the inertial subrange, . It is found that the *r* range where this is true is approximately between 0.3 and 3 m. From this regression, the dissipation rate is computed. The mechanical [] and buoyant destruction terms [] in the TKE budget are then calculated using the dataset assuming MOST is valid in the CSL. If the independently computed production and dissipation rates are in balance, it suggests that MOST may be used in the canopy sublayer and that the TKE budget is, to reasonable approximation, represented by a balance between production and dissipation terms. Figure 5 shows a one-to-one comparison between production and dissipation terms for the high-frequency PI data collected over 12 different days for near-neutral and more stable conditions shown in black and red symbols, respectively. Overall, the terms are reasonably balanced suggesting that MOST may be applied to relevant mean-flow quantities in the CSL.

The boundary layer height was not measured in these experiments, either. However, since the dataset covers highly stable conditions, *δ* can be expected to drop significantly. On the other hand, *δ* was assumed not to appreciably exceed the average neutral ABL height of about 1000 m. This limit was determined from separate measurements in neutral conditions by Salesky et al. (2013). To have a realistic measure of *δ* at stable conditions, Zilitinkevich’s model (Zilitinkevich 1972) is invoked. It is estimated that , where is a constant (Zilitinkevich 1972), , , and rad s^{−1} is the Coriolis parameter. Using this approximation, it is estimated that m under highly stable conditions. As a more improved estimation over the model of Zilitinkevich (1972), another model by Arya (1981) is used and is given by , where . However, this modification is not found to impact the results significantly. Thus, these four datasets provide a wide range of conditions in terms of , stability, and roughness (smooth: wind tunnel and lake experiments; very rough: forest experiments) to assess the proposed model.

## 5. Results and comparisons

### a. Evaluation of the model for

Although two regimes have been identified for the formulation of , we first answer the following question: to what extent can the formulation for zone I (mildly stable regime) explain the experiments? The performance of the formulation for zone II will be explored later. Recall that the formulation for zone I recovers Townsend’s form and explains how the Townsend parameters evolve with *ζ*, which is why exploring this formulation for the entire dataset is of interest. Figure 6 shows the comparison with the wind tunnel experiment conducted by Ohya et al. (1997), where is plotted against . It is to be noted that the gradient Richardson number can be rewritten in terms of the stability parameter as (Kaimal and Finnigan 1994). Three different ratios, 0.1, 0.25, and 0.5, are used in order to cover the data and the comparisons are satisfactory. To cover the range of (0.013–0.13) reported in Ohya et al. (1997), is made dynamic by multiplying by guided by Fig. 4.

For the lake data, the range of *ζ* covered by the lake data spans from 0.001 to 0.025, as shown in Fig. 4, and hence the boundary layer height is estimated not to be small under such mildly stable conditions (recall that the boundary layer height is not measured). Figure 7 shows the comparison between modeled and measured for the lake data when *δ* is calculated based on the model by Arya (1981) as discussed in section 4. As can be seen, the calculated *δ* decreases with increasing from about 1000 m (under near-neutral conditions) to about 100 m under mildly stable conditions. Interestingly, the calculated *δ* from Arya’s model over the lake does not change significantly when . However, it is noted that using a variable *δ* is not appreciably better when compared to calculations with a constant and m.

Figure 8 shows the comparison with field data collected over the two forests. In Fig. 8, all panels show the comparison between measured and modeled (including MOST scaling) in terms of their variations with *ζ*. Ensemble averages of measured and modeled are shown and their variability is represented by black and red error bars, respectively. As expected, (in dimensional form) is found to reduce with stability, since the effects of mechanical turbulence is diminished by the effects of buoyancy. Figure 8a (summer, HW), Fig. 8d (winter, HW), Fig. 8g (summer, PI), and Fig. 8j (winter, PI) show the aforementioned comparisons with a preset boundary layer height m since the data were mostly collected at night time when the boundary layer height is reduced to great extent. For m, the comparisons are found to be satisfactory given that the ensemble averages almost follow each other. There is some unavoidable scatter in the data when they are observed at the low stability regime. However, instrument signal-to-noise ratio is also lower under very stable conditions and cannot be discarded. Figures 8b, 8e, 8h, and 8k in the middle column show these comparisons with a m, which is about the height of the ABL under mildly stable to near-neutral conditions (Stull 1988). With m, the model overestimates at low stabilities. Figures 8c, 8f, 8i, and 8l show the comparisons between measured and modeled from MOST, which is given as (Stull 1988). As shown, MOST predictions often overestimate the data. Figure 9 demonstrates this feature better in which all panels show one to one comparisons between the modeled (*x* axis) (proposed model and MOST) and measured (*y* axis) . As described before, the comparisons with m (Figs. 9a–d) are satisfactory but the proposed model here overestimates the data with a m (Figs. 9e, f). One important point to note for Figs. 7 and 8 is that the measured for each data point is used to drive the model since varies appreciably.

Since the forest dataset spans over such a wide range of stability conditions, the comparisons with the formulation for zone II are also shown in Fig. 10. Figure 10a (summer, HW), Fig. 10c (winter, HW), Fig. 10e (summer, PI), and Fig. 10g (winter, PI) show one-to-one comparisons between the modeled and measured . Figure 10b (summer, hardwood), Fig. 10d (winter, hardwood), Fig. 10f (summer, pine), and Fig. 10h (winter, pine) show variations of with *ζ*. Both one-to-one comparisons and variations with *ζ* show that is also reasonably captured by the formulation for zone II. However, it is again stressed that the boundary layer height is not present in the formulation for zone II, which is different from the formulation for zone I. The errors and biases in the predictions using both formulations for the complete range of stability are discussed in the appendix B. The main finding is that both model formulations are not particularly biased when comparisons are conducted separately for and , presumably because of the fact that the value is derived by matching the two formulations at . Hence, in the zone-II formulation, for instead of an asymptote to 0 (because of a zero additive parameter). Also, the fact that becomes on the order of 0.5 means that the term becomes dominant (instead of ) for in the zone-I formulation. A consequence of this outcome is that the effects of *δ* also become muted in the zone-I formulation.

### b. Scaling variables

An important point to be stressed in Fig. 8 is that the measured and modeled dimensionless quantity is not compared. Rather, is compared directly to avoid the known problem of self-correlation (Hicks 1978, 1981; De Bruin 1982; Pahlow et al. 2001; Andreas and Hicks 2002; Hartogensis and De Bruin 2005; Cava et al. 2008; Banerjee et al. 2014). The reason for such self-correlation is the common presence of in both variables—the ordinate and the the abscissa , where dictates the Obukhov length *L*.

Some previous studies (Pahlow et al. 2001; Juang et al. 2008) reported an increase of with *ζ* and this increase exhibits an approximate 1/3 scaling, which is indicative of possible artificial self-correlation. The forest experiments used here also show a similar 1/3 scaling behavior when is plotted against *ζ*, as illustrated in Fig. 11. Note again that itself reduces with increasing stability as shown in Fig. 8. Hence, to separate out the effect of stability on , the variation of with *ζ* from the model is shown for three different but prefixed (i.e., no artificial self-correlation) in Fig. 12a. Figure 12b shows the variation of nondimensional with *ζ* and since the is fixed (i.e., only the sensible heat flux varies), the curves show a self-similar behavior and decrease with increasing stability. Note the boundary layer height is also fixed in Fig. 12 for illustration though in reality *δ* will vary with increasing stability.

### c. Connection between shape of the stability functions and z-less stratification

In Fig. 12 and also in Fig. 8, it can be observed that becomes independent of *ζ* after , indicating that is no longer dependent on *z* when using the Kansas derived as noted earlier. This observation indicates the possible onset of *z*-less stratification and is commensurate with the same limit reported by other studies such as Olesen et al. (1984) and Kaimal and Finnigan (1994). Figure 13 shows the variation of with *ζ* for the second formulation (zone II) under high stability conditions. While the range of variation of is quite small, it also reaches *z*-less stratification reasonably at (the model results are extended to only for illustrating the near-leveling off when no other constraint is employed to suppress turbulence). As a result, the spectral budget model here predicts the *z*-less scaling for under high stability conditions. The “1/3” scaling of reported in some previous studies (Pahlow et al. 2001) may be the outcome of self-correlation and cannot be used to argue against the existence of the *z*-less stratification in isolation of other metrics. However, it is to be noted that while the behavior of cannot be predicted from and alone, the variation of is heavily dependent on and . Still, the choices of and are not ad hoc and are standard in the literature. MOST predicts an unbounded increase of and with *ζ* so they are not used beyond . Nevertheless, a recent work by Chenge and Brutsaert (2005) provided alternative functional forms for and that level off after . The effect of this variation on the onset of *z*-less scaling is discussed in appendix C. As foreshadowed by Eq. (37), the leveling off in precludes a *z*-less scaling, which is confirmed by detailed calculations shown in appendix C. Hence, the model results here establish a connection between the shape of the stability correction functions for mean velocity and temperature and *z*-less scaling in .

## 6. Conclusions

A spectral budget framework has been developed to describe the longitudinal turbulent intensity under stable atmospheric stratification. In presence of stable stratification, although turbulence is generated by mechanical processes, buoyancy effects counter mechanical effects and consequently reduce the turbulent intensity. Thus, in the spectral budget framework, the mechanical production term and the buoyant destruction term compete with each other. Two stability regimes are identified and separated at . In the first regime (), the resulting normalized longitudinal velocity variance is found to be consistent with the log law proposed by Townsend although the Townsend’s coefficients are found to vary mildly with the atmospheric stability parameters. The second regime () results in a truncated spectrum exhibiting a −5/3 scaling even at low wavenumbers instead of a −1. Interestingly, the second regime formulation does not contain the boundary layer height. Thus, it can be concluded that at , is not sensitive to the boundary layer height.

The model outcomes using the formulation for regime I are compared to data from a controlled laboratory experiment, a field campaign over a smooth lake surface, and two extensive datasets collected over tall vegetation canopies. The formulation for regime II is also compared with the field data over canopies. It is also found that decreases with increasing stability as expected. However, some previous studies had predicted an increase of with increasing stability parameter . This may be partly explained by the artificial self-correlation as a consequence of choosing common scaling variables in abscissa and ordinate. Hence, to separate out the effect of increasing stability on turbulent intensity, is deemed to be a more robust variable compared to . With this choice, the same data that show increasing trend of with display a reduction of with increasing , which is satisfactorily predicted by the model with zone I formulation. The modeled with a prefixed flattens as approaches 0 and recovers the correct neutral limit. The modeled also flattens and becomes asymptotic after . This result may signify that the model supports the existence of *z*-less stratification for a preset and also correctly captures the limit as the onset of *z*-less stratification when the Kansas stability correction functions for momentum and mean air temperature gradients are used. The second formulation as a limiting scenario for high stability regimes predicts an upper limit for the onset of *z*-less stratification after . From an operational perspective, if the boundary layer height is not measured precisely, it can be calculated from a simple model and the predicted would not suffer much, because of the low sensitivity to *δ* (the order of magnitude is important, not the precise value). The proposed framework here brings together seemingly unrelated theories of turbulence such as Kolmogorov’s hypothesis, Heisenberg’s eddy viscosity, Townsend’s hypothesis, and MOST scaling under a common framework. It also provides a physical basis for illustrating links between a number of well-established constants.

## Acknowledgments

T.B. and G.K. acknowledge support from the National Science Foundation (NSF-EAR-134470, NSF-AGS-1102227) and the U.S. Department of Energy (DOE) through the office of Biological and Environmental Research (BER) Terrestrial Ecosystem Science (TES) Program (DE-SC0006967 and DE-SC0011461). DL acknowledges support from the NOAA (U.S. Department of Commerce) Grant NA08OAR4320752 and the Carbon Mitigation Initiative at Princeton University, sponsored by BP. The statements, findings, and conclusions are those of the authors and do not necessarily reflect the views of the NOAA, the U.S. Department of Commerce, or BP.

### APPENDIX A

#### Effect of Imbalance in the TKE Budget

The first-order effect of a stability-dependent imbalance in the TKE budget is now explored. As discussed before, in stable conditions, there might be an imbalance between the production and dissipation terms indicating that the transport terms might be significant. In that case, the spectral budget should also contain terms arising out of Fourier transformation of those terms but it becomes overwhelmingly complicated. Instead, a simple first-order correction accounting for the imbalance is sought that can be perhaps encapsulated by the imbalance term in . The stability variation of this term from the AHATS campaign is compiled by Salesky et al. (2013) and the data are fitted with a power law for the present work (; see inset in Fig. A1). It is to be noted that the neutral limit of this term is zero and thus the formulation I is still valid. The only modification is in the Townsend coefficient :

Figure A1 shows the effect of this term. The top panel shows computed from the first formulation using a that is set to zero (original formulation, solid line) and is set to be a power-law function of *ζ* (dashed line). The bottom panel shows the variation of the associated Townsend parameter with stability for the two approaches. As can be observed, the formulation with the correction is slightly higher compared to the first formulation, but the difference is not highly significant.

### APPENDIX B

#### Biases in the Comparison to Forest Data

The biases arising from using both formulation I and formulation II for the entire range of are now explored. Errors are calculated as the absolute value of the difference between the modeled and measured . Probability density functions (pdf) of errors are shown in Fig. B1. Black lines indicate the pdf of errors in the comparison with the formulation for zone I. Dashed black lines indicate the pdf of errors in the comparison with the formulation for zone II. Dotted red lines indicate the pdf of errors in comparison with MOST. The top row indicates the domain of formulation I () and the bottom row indicates the domain for the formulation II (). As observed, most errors are concentrated on the lower side for our formulation, whereas for MOST, errors are more concentrated on the middle ranges, especially for lower stabilities. Note that for the forest data with , MOST () performs well and compares reasonably with the proposed models. This finding perhaps proves the point that boundary layer height is not so significant at such stable conditions. Also, the stability dependence of is weak. Both of those results are consistent with the formulations for zone II.

### APPENDIX C

#### Effect of Alternative Forms of and

The effect of alternative forms of and are now discussed. The following forms were provided by Chenge and Brutsaert (2005):

where , , , and . Figure C1 shows these functions and in the top panel (black lines) along with traditional MOST forms (blue lines). The bottom panel of Fig. C1 shows a sample variation of computed using and both with MOST (solid line) and Eqs. (C1) and (C2) (dashed line). It is noted that does not flatten after (like MOST) but reduces further with these modified functions. Thus the onset of *z*-less stratification is further “delayed.”

## REFERENCES

*Phys. Fluids*,

**23**, 041702, doi:.

*Phys. Fluids*,

**25**, 125106, doi:.

*Quart. J. Roy. Meteor. Soc.*,

**141**, 1699–1711, doi:.

*J. Geophys. Res.*,

**114**, D08124, doi:.

*Turbulence*. McGraw-Hill, 586 pp.

*Water Resour. Res.*,

**45**, W08431, doi:.

*Atmospheric Boundary Layer Flows: Their Structure and Measurement.*Oxford University Press, 289 pp.

*Phys. Rev. Lett.*,

**107**, 268502, doi:.

*Phys. Rev.*,

**86E**, 066311, doi:.

*Phys. Rev.*,

**87E**, 023004, doi:.

*Phys. Rev.*,

**89E**, 023007, doi:.

*Phys. Fluids*,

**24**, 105105, doi:.

**157**, 1–21, doi:.

*J. Atmos. Sci.*,

**72**, 2394–2410, doi:.

*Phys. Fluids*,

**10**, 855–858, doi:.