## Abstract

This paper presents a turbulence closure for neutral and stratified atmospheric conditions. The closure is based on the concept of the total turbulent energy. The total turbulent energy is the sum of the turbulent kinetic energy and turbulent potential energy, which is proportional to the potential temperature variance. The closure uses recent observational findings to take into account the mean flow stability. These observations indicate that turbulent transfer of heat and momentum behaves differently under very stable stratification. Whereas the turbulent heat flux tends toward zero beyond a certain stability limit, the turbulent stress stays finite. The suggested scheme avoids the problem of self-correlation. The latter is an improvement over the widely used Monin–Obukhov-based closures. Numerous large-eddy simulations, including a wide range of neutral and stably stratified cases, are used to estimate likely values of two free constants. In a benchmark case the new turbulence closure performs indistinguishably from independent large-eddy simulations.

## 1. Introduction

Atmospheric motion occurs from the planetary scales, limited from above by the size of the earth, to the turbulent scales, which are limited from below by the molecular viscosity of air. All scales are important for a complete understanding of the atmosphere. However, it is neither feasible nor desirable to resolve all motions in models of the atmosphere. Rather we prefer to parameterize the motions, which we are not explicitly resolving. Prime examples that require parameterization are convection, cumulus clouds, gravity waves, and turbulence. The latter is the subject of the present study.

The purpose of a turbulence closure model is to predict the tendencies of the small-scale turbulent disturbances on the large-scale mean flow. A classical turbulent boundary layer approach is obtained by assuming horizontal homogeneity and applying the Boussinesq approximation (e.g., Stull 1988). Prognostic equations for the mean-state wind and the mean potential temperature of the atmospheric turbulent boundary layer in their Reynolds-averaged form are

where and are the vertical momentum fluxes, is the vertical potential temperature flux, *f* is the Coriolis parameter, *D/Dt* denotes the total derivative, and *U _{g}* and

*V*are the zonal and meridional components of the background geostrophic wind vector. Uppercase letters denote means, while lowercase letters are turbulent departures from the mean. Obviously, numerical weather prediction and climate models also account for tendencies of, for instance, radiative and moist processes.

_{g}To solve the coupled equation system (1)–(3) it is necessary to parameterize the turbulent fluxes. This can be done in many different ways, but the most common is to define turbulent conductivity, = −*K*_{h}(∂Θ/∂*z*); and turbulent viscosity, *τ* = *K _{m}S,* where

**= −(, ) is the turbulent stress vector and**

*τ***S**= (∂

*U*/∂

*z*, ∂

*V*/∂

*z*) is the shear vector. The most simple case is to assume that

*K*and

_{m}*K*are equal and constants, in which case the problem decouples such that only (1) and (2) need to be solved simultaneously. Ekman (1902) found a spiraling flow steady-state analytical solution for the case of constant turbulent viscosity.

_{h}To properly account for the effects of the mean flow on turbulence a myriad of turbulence closure models have been developed. First-order closures are formulated in terms of the mean variables themselves, while higher-order closures add one or more prognostic variables to the problem; for instance, turbulent kinetic energy is a popular choice (see Holt and Raman 1988 for a general review). First-order closures tend to be used mostly in general circulation models, in operational weather forecasting, and in climate studies. They bear the advantage of being computationally efficient and most often quite insensitive to the rough vertical and temporal resolutions applied (Ayotte et al. 1996). Prime examples are based on the Monin–Obukhov scaling (Monin and Obukhov 1954), which relates the stability of the flow to a nondimensional ratio of the turbulent heat flux to the turbulent stress (e.g., Louis 1979). Unfortunately, Monin–Obukhov scaling suffers from the statistical problem known as self-correlation (Hicks 1978; Mahrt et al. 1998; Klipp and Mahrt 2004). Higher-order closures tend to be used in mesoscale to regional weather forecasting and in research, though there are exceptions. Most higher-order closures owe their legacy to the second-order turbulence closure models presented by Mellor and Yamada (1974), which treat all second moments of the flow.

Efforts have been made to understand the differences between the many different closures within the Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS). The first intercomparison case was an idealized setup with a constant background geostrophic wind and a slowly cooling surface. A weakly stable boundary layer grew against a background stratification, qualitatively resembling an Arctic stable case (Kosovic and Curry 2000). The intercomparison included both 11 large-eddy simulations (LESs; Beare et al. 2006) and 20 different turbulence closure schemes (1D models) from operational weather services and research-type models (Cuxart et al. 2006). Whereas the LES agreed on a boundary layer depth of 150–200 m, depending on resolution and subgrid-scale model, the 1D models varied between 120 and 483 m. The deepest boundary layer was found for the scheme presented by Viterbo et al. (1999), as implemented in the operational model at the European Centre for Medium-Range Weather Forecasts (ECMWF). The scheme was designed deliberately to be more diffusive in stably stratified conditions than what micrometeorological observations indicate in order to improve forecasts nighttime minimum temperatures and synoptic cyclones (Viterbo et al. 1999; A. Beljaars 2007, personal communication). This choice unfortunately leads to overprediction of the depth of the boundary layer.

The boundary layer height (*H*) is an integral measure, usually defined as the level where the turbulent stress vanishes. Predicting the *H* may not seem practically important, for instance, compared to the 2-m temperature or the development of a synoptic storm. However, Cuxart et al. (2006) demonstrated that the turbulent surface stress and downward surface heat flux is nearly proportional to *H* within the models tested in the first GABLS case. Further, Svensson and Holtslag (2007, manuscript submitted to *Quart. J. Roy. Meteor. Soc.*) found that models predicting deeper boundary layers compared to LES had a too-small surface wind angle with the geostrophic wind and too-large cross-isobaric Ekman transport. Obviously, these deficiencies are bound to have severe consequences for weather and climate predictions, assuming LES results bear some resemblance with the behavior of the real atmosphere.

The first GABLS intercomparison considers a particular idealized case, whereas the real atmosphere experiences a wide range of conditions. Hence, successful emulation by a turbulence closure of a single case does not prove its worth as a generally applicable parameterization. Instead Esau and Zilitinkevich (2006) generated a database of idealized LES cases, both neutral and different types of stable boundary layers with *H* ranging from 12 to 1600 m. Given the same initial and boundary conditions as the LES we applied the turbulence closure suggested by Viterbo et al. (1999). The results are compared with the LES and plotted as gray symbols in Fig. 1; details of the figure will be given later. The figure shows that the Viterbo et al. (1999) closure performs quite acceptably for the deeper (>500 m) boundary layers. However, as stability increases and the boundary layer becomes shallower, this closure systematically overestimates *H*, in the worst case by an order of magnitude. A different version of the closure based on Monin–Obulhov scaling, intended to resemble micrometeorological observations, was also tested and found to yield overprediction of the stable boundary layer *H*, in the worst case only by a factor of 3. In the present study we present a new turbulence closure model, based on the prognostic total turbulent energy equation, with the aim of improving modeling of both the neutral and the stably stratified boundary layer.

## 2. Model

A common approach to the turbulence closure problem is to apply additional prognostic equations to the mean-state (1)–(3), so-called higher-order closures. Usually the prognostic turbulent kinetic energy (*E _{k}*) equation is applied [(A1)]. Richardson (1920) found, using the

*E*equation, that beyond a certain stability limit, which was later to be known as the critical Richardson number, turbulence would decay. Here, we consider instead the total turbulent energy (

_{k}*E*=

*E*+

_{k}*E*), which is the sum of both the turbulent kinetic and the turbulent potential energies. The latter is proportional to the density variations in the fluid, which can be expressed in the potential temperature variance (Zilitinkevich et al. 2007)

_{p}where *σ*^{2}_{θ} is the potential temperature variance, *β* = *g*/Θ is the buoyancy parameter, *N* ^{2} = *β*∂Θ/∂*z* is the squared Brunt–Väisälä frequency, and the vertical bars denote the absolute value. The two budget equations, for *E _{k}* and

*σ*

^{2}

_{θ}(

*E*), have been used in a turbulence closure model proposed by Zilitinkevich (2002). To derive the total turbulent energy equation one must multiply the temperature variance equation by

_{p}*β*

^{2}/

*N*

^{2}and add the result to the turbulent kinetic energy equation (appendix A). The total turbulent energy equation accordingly reads

where *γ* is the dissipation rate of *E* and *F _{E}* is the third-order flux. The terms on the right-hand side we refer to as shear production, dissipation, third-order flux divergence, and buoyancy production terms, respectively. The principle of (4) is illustrated in Fig. 2. Both the sources,

**·**

*τ***S**and 2

*β*, and the sink,

*γ*, owe their existence to the presence of turbulence such that reality is more complicated than it appears from Fig. 2. It is possible to solve a steady-state version of (4) instead, which may be beneficial in models that require long time steps.

In neutral and stratified, steady-state, horizontally homogeneous conditions the first three terms on the right-hand side of (4) should balance. In most situations the third-order flux divergence, which acts to redistribute *E* vertically, is much smaller than the shear production and dissipation terms. In neutral conditions *E _{p}* is actually zero, so here the prognostic equation for

*E*, (A1) applies in itself. However, the major step forward in applying (4) is in stably stratified conditions, where the term

_{k}*β*, which is negative, cancels out. The term is often referred to as buoyancy destruction in stratified conditions. From the present analysis it seems more appropriate to name it buoyancy redistribution, as the term merely serves to transfer energy from

*E*to

_{k}*E*in stratified conditions. In statically unstable conditions the buoyancy production term appears in (4). This means that even in the absence of shear, there will be production of

_{p}*E*. It is, however, necessary to include a parameterization of the nonlocal effects of convection (e.g., Troen and Mahrt 1986). The main advantage of using the total turbulent energy is that Richardson’s (1920) critical limit in stable stratification disappears.

## 3. Turbulence closure assumptions

To solve the system of prognostic equations (1)–(4), we need five turbulence closure assumptions. This is because the system contains the unknown variables , , , *γ*, and *F _{E}*, besides the prognostic variables. The latter two unknown variables arose when we introduced (4). Adding additional prognostic equations to the problem introduces even more unknown variables, known as the turbulence closure problem. Here, we have decided to keep things as simple as possible while benefiting from the physics contained in (4). Further simplification can be obtained by assuming steady state in (4), making it a diagnostic equation.

For the second-order vertical fluxes, we utilize observed properties of turbulence at weak and strong stabilities by Mauritsen and Svensson (2007). They plotted several nondimensional entities, including normalized fluxes and variances, from six different atmospheric observational datasets as functions of the gradient Richardson number:

where *S*^{2} is the squared magnitude of the shear vector. Only positive values of Ri were considered, as the main focus is stable stratification. In short, Mauritsen and Svensson (2007) found that stably stratified turbulence is very different in weak (Ri < 0.1) and strong stability (Ri > 1). In weak stability the turbulent fluxes are proportional to the variance; that is, the turbulent eddies are actively transporting both momentum and heat. However, in strong stability the normalized turbulent stress is diminished to a fraction of the weakly stable level, while the normalized turbulent heat flux is statistically indistinguishable from zero. Note that this does not necessarily imply that the heat flux itself is zero. It may be explained by the temperature variance growing at a higher rate than the heat flux with increasing stability. Between the two regimes a rapid transition occurs in the range 0.1 < Ri < 1. Here, we use the nondimensional stress and heat flux

In Fig. 3 we plotted these quantities as functions of Ri; details of the analysis can be found in Mauritsen and Svensson (2007). Simple stability functions were subjectively fitted to the observations, *f _{τ}* = 0.17[0.25 + 0.75(1 + 4Ri)

^{−1}] and

*f*= −0.145(1 + 4Ri)

_{θ}^{−1}. These functions are shown in the plot. Given these functions and the turbulent variances (see appendix A), it is easy to calculate the fluxes of heat and momentum from (6). In the latter case, we assume the stress vector is aligned with the wind shear. It is common to formulate turbulence closures in terms of turbulent viscosity and conductivity, rather than fluxes, primarily for numerical reasons. These may be obtained analytically from the presented closure (see appendix B).

We assume the presence of a dissipation range for turbulent energy, following the ideas of Kolmogorov (1941; see, e.g., Frisch 1995), such that the dissipation rate can be approximated by

where *C _{γ}* is an empirical constant and

*l*is the dissipation length scale. It is, however, noteworthy that

*C*is not a free constant. In truly neutral, stationary, horizontally homogeneous conditions, neglecting vertical energy transport (

_{γ}*F*= 0), we have

_{E}*τ*=

*S*

^{2}

*l*

^{2}and

*E*=

*E*and using the budget (4) and the definition (6) we get

_{k}*C*=

_{γ}*f*(0)

_{τ}^{3/2}≈ 0.07. We approximate the dissipation length scale by a multilimit formulation as follows inspired by Blackadar (1962):

which takes the distance from the ground, the Coriolis effect, and static stability into account (Rossby and Montgomery 1935; Zilitinkevich 1972; Nieuwstadt 1984; Hunt et al. 1985; Zilitinkevich and Mironov 1996; Zilitinkevich and Baklanov 2002; Zilitinkevich and Esau 2005). Whereas the above-mentioned studies mostly considered the bulk of the boundary layer, we here intend to resolve the boundary layer. As a consequence we use local *τ* and *N*, rather than the surface stress and background stability, in the formulation of the length scale. This approach is more flexible and allows for decoupled, elevated turbulence in, for example, resolved low-level jets, breaking gravity waves, and turbulence in baroclinic flows. The von Kármán constant, *k*, is here taken to be 0.4. The free constants, *C _{f}* and

*C*, will be found by adjusting the model results to fit reasonably with multiple LES results.

_{N}Finally, for the turbulent energy flux we use *l* as an approximate mixing length for turbulent energy, such that

The effect of the third-order flux divergence in this parameterization is therefore local and only active where there is a nonzero curvature in the turbulent energy profile. The effect of the term is small in the steady state but may be important in developing boundary layers. A more advanced closure for *F _{E}* have been suggested by Zilitinkevich (2002) to account for the effect of gravity waves, which are emitted from the boundary layer to the free atmosphere. Further, in convective conditions nonlocal transport of

*E*occurs due to the vertical advection by thermals. Here we refrain ourselves from such further complications.

## 4. Numerical solution

Equations (1)–(4) are conveniently solved on a vertically staggered grid. The first level is a mass level, containing *U*, *V*, and Θ, at the roughness height, *z*_{0}. At this level we assume *U* = *V* = 0. The vertical resolution was usually on the order of a few meters, though coarser resolutions were tested. Between the mass levels, we placed the turbulence levels, where *E*, , , and were calculated. The computational grid is depicted in Fig. 4. Straightforward linear finite difference approximation was used everywhere, except at the first turbulence level, where we applied the logarithmic approximation (e.g., Arya 1991). We used explicit time stepping.

Close to the surface the velocity gradient is governed by the turbulent stress and the distance to the ground, since *l* ≈ *kz* as can be seen from (8); hence *S* ≈ *τ*/*kz* and it is easy to show that the wind profile is logarithmic, *τ* = *kU*(*z*)/ln(*z*/*z*_{0}) (e.g., Landau and Lifshitz 1987). In practice, however, the first level is often at a considerable fraction of the entire boundary layer, and therefore we need to compensate for flow stability in the same way as in the rest of the model. We assume steady state and replace *kz* → *l*_{k} → *l *[ *f*_{τ}(Ri)/*f*_{τ}(0)], where the turbulent mixing length with respect to momentum, *l _{k}*, is equal to the dissipation length in neutral conditions but smaller when Ri > 0, in accordance with observations by Tjernström (1993). Hence, the average stress between the first and second mass levels for the general case is

where the first fraction on the right-hand side is the logarithmic finite difference approximation of *S*^{2} taken at the first turbulence level, *z _{t}*(1), and the meaning of the symbols can be found in Fig. 4. Note how (9) can be rewritten into the form

*τ*(1) =

*C*[

_{m}*U*(2)

^{2}+

*V*(2)

^{2}], where

*C*is the well-known surface drag coefficient. As such, in neutral conditions, close to the surface the formula reverts to the logarithmic wind law. However, as

_{m}*l*<

*kz*and Ri > 0 the first-level stress decreases. By using (9) the results are less sensitive to resolution, as it provides a seamless transition from the logarithmic layer near the surface and the rest of the boundary layer. Note that

*τ*(1) is not the surface stress. Unlike many surface layer formulations, we do not assume the presence of a constant flux layer. If the first mass level is at some fraction of the entire boundary layer, there may be considerable flux divergence in the layer below. Therefore, a better way to diagnose the turbulent surface stress is to linearly extrapolate from the first two turbulence levels,

*z*(1) and

_{t}*z*(2), to the surface. This was not necessary, though, in the present study.

_{t}Next, we consider the heat transfer in the layer close to the ground. Again we seek the average flux between the lowest mass levels; for generality, we place the lowest level at the roughness height with respect to heat, *z _{T}*

_{0}. The heat flux, , is governed by the temperature gradient, ∂Θ/∂

*z*; the turbulent stress,

*τ*(1); and the turbulent length scale,

*l*≈

*lz*. In this case Landau and Lifshitz (1987) showed that = (∂Θ/∂

*z*)[

*kz*

*τ*/Pr(0)], where the additional parameter, Pr(0) =

*K*/

_{m}*K*, is the turbulent Prandtl number in neutral conditions. In terms of the present model Pr(0) ≈ 0.69 (see appendix A). The effective Prandtl number increases implicitly with increasing Ri within the model. Generalizing as above to some finite height above ground and compensating for the flow stability, we get

_{h}Again we must stress that (1) is not the surface heat flux. Rather the latter should be found by extrapolation. Note how the expression, by insertion of (9), again can be rewritten into the traditional form (1) = *C*_{h}[Θ(2) − Θ(1)]*U*(2)^{2} + *V*(2)^{2}, where *C _{h}* is the surface heat transfer coefficient. If the temperatures at the first two mass levels are given, it is straightforward to apply (10). Such a case is presented in section 6. In applications where the surface condition is the surface heat flux, it is necessary to linearly interpolate between the surface level,

*z*

_{T}_{0}, and the second turbulence level,

*z*(2), to obtain (1). In the latter case the surface temperature, Θ(1), can be diagnosed by rearranging (10). This technique is used in section 5.

_{t}## 5. Determining the free constants

The turbulent dissipation length scale contains two unknowns or free constants (8). One constant, *C _{f}* , accounts for the dampening effect of the earth’s rotation on the largest eddies and the other, C

*, for the suppression by the mean flow stratification. Ideally, data from the atmosphere should have been used to determine these constants. However, the atmosphere seldom obeys the criteria of being horizontally homogeneous and stationary on the synoptic scale as is convenient to assume in the model. To cover a large range of background conditions, a particular field experiment would not be enough. Furthermore, to be able to constrain the constants, information from the surface properties, through the turbulence profiles and background stability and flow in the free atmosphere, are needed. If not, the uncertainties in the experimental data will easily overpower the estimations. For the purpose of determining these constants, we have instead utilized a large database of LES cases presented by Esau and Zilitinkevich (2006). Here the background conditions were varied systematically in an attempt to span a major part of the parameter space for the horizontally homogeneous, neutral, and stably stratified dry atmospheric boundary layer. The LES cases consist of four types of boundary layers and are summarized in Table 1. Normalized profiles of the mean and turbulent variables are plotted in Fig. 5.*

_{N}### a. LES database

In all cases an initial temperature profile, either neutral or with constant stratification, a constant background geostrophic wind, the surface roughness length, and the surface heat flux were prescribed. Each case was run for 15 h to achieve a quasi-steady state. The truly neutral class has uniform potential temperature and zero heat fluxes. While this class is somewhat academic, in the sense that truly neutral boundary layers hardly ever occur in the real atmosphere, it provides an important limit in parameter space. More realistic are the conventionally neutral and nocturnal classes. The former is neutral in the sense that the surface heat flux is zero, but it is growing against a stably stratified atmosphere. As a consequence the lower part of the boundary layer is well mixed, while at the same time it is capped by a stably stratified elevated inversion. These cases are representative of windy situations when the surface heat flux is negligible. The nocturnal boundary layers grow against a neutral atmosphere, while heat is lost at the surface. These boundary layers occur during nighttime over land where the near-neutral residual layer is remnant from the daytime convective boundary layer and the surface is cooling radiatively. Finally, the surface cooling and background stratification is combined in the long-lived stable boundary layer class. These boundary layers are frequently found at high latitudes, over land during wintertime, and over large glaciers. The bulk stability parameters, *L*, *L _{f}* , and

*L*, as described in Zilitinkevich (2002), along with the boundary layer heights, indicate the wide stability range covered by the database (see Table 1).

_{N}The LES used a dynamic subgrid-scale closure model and a resolution of 64^{3} grid points. This is moderate with present standards but necessary due to the large number of cases. Further, Esau and Zilitinkevich (2006) found only insignificantly different results at higher resolution for a limited number of cases. The grid spacing varied from case to case in order to well accommodate the boundary layer in its total depth. Quality of the included cases was ensured: cases where the LES domain was chosen smaller than 1.5 times *H* or when the turbulent fluxes did not tend to zero at the domain top were removed. The LES did not conserve heat exactly, unlike the presented model, and therefore cases when the integrated cooling did not correspond to within 20% of the imposed surface cooling were rejected for the nocturnal and long-lived classes. For the conventionally neutral class this objective criterion was not applicable, since the surface cooling was zero. Here cases, which did not conserve heat, were subjectively removed. In most cases these coincided with the ones that had too-shallow domains to accommodate the entire boundary layer. After quality control we were left with 90 cases in total.

### b. Strategy for determination of C_{f} and C_{N}

We decided to use the boundary layer height as a metric for determining the free constants. First, the truly neutral cases were used to determine *C _{f}* . Figure 6 shows the average

*H*normalized by the LES heights,

*H*

_{LES}, for varying

*C*. It appears that the truly neutral

_{f}*H*is proportional to the square root of

*C*and that the optimal value is about

_{f}*C*= 0.185. Also plotted are 95% confidence intervals indicating a possible range of 0.16–0.22, while the outliers indicated a range of 0.12–0.27 for

_{f}*C*. The use of confidence intervals is questionable because we cannot assume Gaussianity of the error distribution. Therefore the latter, more cautious, uncertainty range should be used. Setting

_{f}*C*= 0.185 we proceeded to determine

_{f}*C*with the other three classes. Figure 7 shows the normalized boundary layer heights for varying

_{N}*C*for each class. Class optimal values are 1.3 for the conventionally neutral, 1.8 for the long-lived, and 2.5 for the nocturnal class, respectively. The largest sensitivity to

_{N}*C*is for the long-lived class, which also includes the most shallow boundary layers. This means that perfect optimization across the classes cannot be achieved. However, zero bias against these LES cases is obtained for

_{N}*C*= 2.0. Here

_{N}*H*for the conventionally neutral and long-lived classes are overpredicted by 5% and the nocturnal class is underpredicted by 8%. Again, the uncertainty range is quite wide, roughly between 1 and 3. Also,

*C*have only been measured within the concept of mixing length (Nieuwstadt 1984; Hunt et al. 1985; Mahrt and Vickers 2003), finding values of

_{N}*C*

^{mix}

_{N}of 0.5–1. Remembering that the mixing length is shorter than the dissipation length by a factor of 2 to 4 for stratified cases (Tjernström 1993), these values are consistent with our findings.

Modeled boundary layer heights, using *C _{f}* = 0.185 and

*C*= 2.0, are shown in Fig. 1. It appears that the model predicts

_{N}*H*in individual cases to within 20%. Closer inspection reveals that the conventionally neutral and long-lived cases are overpredicted, while the nocturnal cases are consistently underpredicted. Figure 5 graphs modeled normalized profiles of mean wind, potential temperature, and turbulent fluxes for each class. The normalization is done using the LES values as described in the figure caption in order to emphasize the biases between the model and the LES. In general the agreement is good, as the model is able to capture the overall vertical structure in the different classes. The heat flux profiles most clearly show the biases. A deeper boundary layer yields a stronger downward heat flux and, consequently, a warmer surface temperature. This is the case for the conventionally neutral and long-lived classes, while the opposite is the case for the nocturnal class. Finally, we plotted the surface wind angle to the geostrophic wind vector in Fig. 8. Large surface wind angles are usually associated with shallow boundary layers (Csanady 1974). Here performance is good up to about 40°, while the model underpredicts compared to LES at larger angles. It is likely that even the LES is underpredicting the extreme angles by a few degrees, as some wind turning may occur below the lowest computational level. The more diffusive Viterbo et al. (1999) closure seems to be unable to produce surface wind angles larger than about 30°.

## 6. The first GABLS case

Within GABLS, an idealized case was set up to resemble a typical Arctic stably stratified dry boundary layer. A constant geostrophic wind of 8 m s^{−1} allowed an initially 100-m-deep neutral boundary layer to grow against a background stratification of 0.01 K m^{−1} aloft. At the same time the surface temperature was prescribed to decrease by 0.25 K h^{−1}. The experiment was run for 9 h, achieving a quasi-steady state. In principle this is a long-lived case as described earlier, except that the surface temperature is prescribed rather than the surface heat flux. Eight LESs were run with a resolution of 3.125 m (Fig. 9). The vertical potential temperature distribution exhibits a surface inversion, followed by a less stratified mixed layer, which in turn is capped by an elevated inversion. The wind profile shows a typical Ekman spiral with a jet in the vicinity of the capping inversion.

Along with the LES results, we also plotted the results from the presented model and our implementation of the Viterbo et al. (1999) turbulence closure model. The latter overpredicts the boundary layer depth compared to LES roughly by a factor of 2, in good agreement with the earlier results of Cuxart et al. (2006), despite the differences in implementation. Hence, the deficiency of the scheme is neither related to the resolution nor the numerical method but rather the physical parameterization itself. The high-resolution version of the presented model (4), on the other hand, performs indistinguishably from the LES.

We also performed tests with sparse vertical grids. An example is shown in Fig. 9 using a stretched grid with mass levels at 30, 78, 155, 278, and 474 m, respectively. This resolution is typical for operational models, though not specific to any particular model. Even though the model now hardly resolves the inversion and the jet, there is a good resemblance with the high-resolution version. Essential for this good agreement is the boundary conditions (9) and (10), which extend the turbulence closure scheme properly all the way to the surface. Obviously, the results are bound to degrade for more shallow boundary layers or with poorer resolutions. Note how a linear extrapolation of the first- and second-level turbulent fluxes down to the surface is a good approximation of the surface fluxes.

The approach we used to determine the free constants provides us with uncertainty intervals rather than definite values. As an illustrative example, we may investigate the sensitivity of the model to the choice of *C _{N}* in the first GABLS case. The potential temperature profile is shown in Fig. 10. The three values of

*C*correspond to the three points of intersection with unity for the individual classes in Fig. 7. The constant

_{N}*C*only impacts the inversion layer because the surface temperature is fixed. The smaller

_{N}*C*provides a sharper inversion, while the larger

_{N}*C*allows for a deeper inversion layer. The spread among the different realizations of the model is about the same as the spread among the LES, though there is a tendency toward a better qualitative resemblance for the smaller values. While the sensitivity to the choice of

_{N}*C*may seem to be small, it may have a large impact on more stably stratified cases and, for instance, the evolution of shallow cloud layers as it controls the amount of entrainment.

_{N}## 7. Discussion

A turbulence closure model is founded on closure assumptions or hypotheses. These hypotheses can take many different forms. Hence the problem leaves room for several subjective choices, with respect to both formulation and constants. We have presented a new turbulence closure model, based on (4), and determined optimal values of the free constants with respect to the average bias on the boundary layer height. Obviously, this does not ensure reasonable physical behavior of the model. For example, a closure model that gives a constant boundary layer depth could have yielded optimal values of the free constants on a finite dataset. Success in a single case, such as the first GABLS case, also does not provide a sufficient justification for further interest in yet another turbulence closure. Investigating the detailed results of the cases in the LES database does not provide independent evidence for its value nor does comparison with the complex reality give good guidance to the worthiness of the model due to the uncertain initial and background conditions. So can we be sure the new model is of any value?

Quantitative, controlled, and detailed observations of turbulence in the stably stratified boundary layer are very difficult to obtain (e.g., Mahrt 1999). Of particular interest for the purpose of the model validation is the vertical structure of turbulent fluxes. Models for the stable boundary layer are commonly validated by normalized, observed profiles by Caughey et al. (1979) and Nieuwstadt (1984) as in, for example, Andrén (1990). In the former case the heat flux profile was strongly concave, while in the latter it was slightly convex. It is not possible to tell whether steady state has been achieved nor to which extent horizontal heterogeneity and topography influenced the results. In fact, quasi-steady state is seldom achieved in the real world nocturnal boundary due to the passing nature of night and the wealth of natural phenomena and variability (e.g., Mahrt 1999; Poulos et al. 2002). The presented model is able to provide both concave heat flux profiles early in the nocturnal cases, while convex profiles similar to Nieuwstadt’s are found in some of the long-lived cases. The latter can also be obtained in the nocturnal case if the surface reaches an equilibrium temperature. This only demonstrates the importance of knowing the actual initial conditions, forcings, and transient behavior when comparing with observations.

Our first hypothesis is that the turbulent fluxes of momentum and heat are related to the turbulent variances through unique functions, which depend on the local gradient Ri only. As a hypothesis it has already failed, as several of the datasets shown in Fig. 3 disagree significantly with the empirical curves (Popper 1934). Other parameters than Ri may be of importance. However, it is a matter of debate if these differences are due to physical differences or whether they reflect different measuring and dataset preparation techniques (Mauritsen and Svensson 2007). Our second hypothesis is that the turbulent energy dissipation can be described by a dissipation length scale according to (8). As a hypothesis this choice is less fortunate, as testing this hypothesis is difficult if not impossible with present-day atmospheric observational systems.

The process of determining the model constants after the model has been run is often referred to as “tuning” (Randall and Wielicki 1997). The use of tuning is disliked in the scientific community, for example, because tuning of models could mask fundamental errors in our understanding of the problem. Hence, tuning can be detrimental to the scientific development. However, given the importance of the processes, clearly seen in Figs. 6 and 7, we are left with no choice but to tune the constants. While tuning the performance of a complex model to some dataset, one should be aware of the dataset limitations and biases. For instance, the LES cases contain forcings and boundary conditions, which were chosen rather subjectively. Further, the approach allows for a number of subjective choices, such as the metric with which we assess the performance. The model that we are tuning could also be formulated in a slightly different manner with the sole purpose of resembling the LES results. However, given the large number of cases that has been used, it would not have been possible to tune the model using only two constants in the absence of fundamentally sound physics.

## 8. Conclusions

An essentially new turbulence closure scheme for neutral and stably stratified atmospheric boundary layers is presented. The scheme is based on the prognostic turbulent energy equation, which is the sum of the turbulent kinetic and turbulent potential equations. The use of this equation is preferable over the turbulent kinetic energy equation, because it predicts turbulence in the presence of wind shear, regardless of stratification.

The effect of stable stratification is to reduce turbulence and diminish the vertical turbulent fluxes. The scheme is founded on simple relations between the turbulent variances and fluxes and the idea of a characteristic turbulent length scale. Together, these relations account for the effects of stratification, within the framework of the presented model. Observations of the ratios, plotted as functions of the gradient Richardson number, avoid the problem of self-correlation. This is an improvement over the widely used Monin–Obukhov-based turbulence closure models. Surface conditions are derived, which are essentially an extension of the closure to the layer close to ground. Here we avoid the assumption of the presence of a constant flux layer.

Two free constants were undetermined by micrometeorological observations. Therefore the scheme was tuned, or calibrated, against a large number of LES cases. The method may be rightfully questioned, because in principle the model is forced to give the result we desire. However, without fundamentally sound physics it would probably not have been possible to mimic 90 LES cases only by varying two constants. The procedure provides us with guidance in the form of intervals with likely values of the constants.

The performance of the turbulence closure model is tested in an independent case against other LES. The resulting profiles of high-resolution integrations are indistinguishable from the LES results. Using a coarse vertical resolution, typical of operational models, does not degrade the performance significantly for this particular case.

## Acknowledgments

The authors wish to thank Bob Beare, Ye Cheng, Boris Galperin, Bert Holtslag, Nathan Kleeorin, Peter Lundberg, Anne McCabe, Igor Rogachevskii, Gert-Jan Steeneveld, Semion Sukoriansky, Michael Tjernström, and the members of the GABLS working group for contributing to the present study. Support for this work came from The Swedish Research Council, The Norwegian Research Council, and the activity in frameworks of the Nordic research network on boundary layer climates.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

### APPENDIX A

#### Total Turbulent Energy

The primary aim of this appendix is to derive (4) for the general horizontally homogeneous case. Second, we approximately diagnose the relative contributions of *E _{k}* and

*E*given

_{p}*E*. Among the equations for the second-order moments (e.g., Stull 1988), we study two equations, those for turbulent kinetic energy,

*E*, and potential temperature variance,

_{k}*σ*

^{2}

_{θ}:

where *ε* and *ϕ* are the dissipation rates, and *F*_{k} = + /*ρ* and *F*_{θ} = /2 are the third-order vertical fluxes of *E _{k}* and

*σ*

^{2}

_{θ}, respectively.

We note that in stable stratification < 0 and ∂Θ/∂*z* > 0. Hence, in this case the second term on the right-hand side of (A1) is negative, while the first term in (A2) is positive. Next, in unstable stratification > 0 while ∂Θ/∂*z* < 0. Thus, both buoyancy terms in (A1) and (A2) are positive. Keeping this in mind we may multiply (A2) by *β*^{2}/|*N* ^{2}| and use the definition of *N* ^{2} to get

We then add (A1) to (A3) in order to obtain (4). Here we have defined the dissipation rate *γ* = *ε* + *ϕβ*^{2}/|*N* ^{2}| and the third-order energy flux *F _{k}* =

*F*+

_{E}*F*

_{θ}β^{2}/|

*N*

^{2}|. We parameterize

*ε*=

*C*

_{ε}

*E*

_{k}

*E*/

*l*and

*ϕ*=

*C*

_{ϕ}½

*σ*

^{2}

_{θ}

*E*/

*l*, assuming the same dissipation length scale applies to

*E*,

*E*, and

_{k}*E*. In neutral conditions it is clear that

_{p}*C*

_{ε}=

*C*, where

_{γ}*C*was the dissipation constant with respect to

_{γ}*E*. We shall also assume

*C*, which simplifies the following derivations slightly.

_{ϕ}= C_{γ}In the present study, we only solve (4) not (A1) and (A3) separately. Thus, we need to diagnose *E _{k}* and

*σ*

^{2}

_{θ}to calculate the fluxes of momentum and heat. Assuming steady state and neglecting vertical energy transport, dividing (A3) by (A1) and using the above definitions we obtain

When Ri → −∞, that is, |** τ** ·

**S**| ≪

*β*||, we get

*E*/

_{p}*E*→ 1. In the neutral limit, when Ri → ±0, we have |

_{k}**·**

*τ***S**| ≫

*β*||. Using the definitions of

*K*and

_{m}*K*and that the turbulent Prandtl number is Pr =

_{h}*K*/

_{m}*K*, we obtain

_{h}where Pr(0) = ( *f* ^{2}_{τ}/2*f* ^{2}_{θ}) is the neutral limit Prandtl number in terms of the present model (see appendix B). If *C _{ϕ}* ≠

*C*

_{ε}a factor of

*C*/

_{ϕ}*C*

_{ε}appears in the expression of Pr(0). Finally, we consider the limit Ri → +∞. Recalling Fig. 2, the sole source of

*E*in stably stratified conditions is the shear production. Hence it is not possible that the buoyancy redistribution term exceeds the shear production. Recognizing that (A4) is equal to

*E*

_{p}/

*E*

_{k}=

*R*

_{f}(1 −

*R*

_{f})

^{−1}, where

*R*

_{f}is the flux Richardson number, and taking the limit

*R*

_{f}= 0.3 ± 0.1 at Ri → ∞ after Zilitinkevich et al. (2007) and Esau and Grachev (2007), we found an upper bound of (A4) for very stably stratified conditions as

This limit is in good agreement with observations presented by Mauritsen and Svensson (2007) and somewhat higher than Schumann and Gerz (1995). For a field of linear internal gravity waves, in a nonrotating fluid, the ratio is unity (e.g., Nappo 2002). The difference is that in the latter case the energy source is external to the flow. In the present study we diagnose the relative contribution of *E _{p}* to

*E*by simply interpolating between the limits (A5) and (A6).

### APPENDIX B

#### Eddy Viscosity and Conductivity

Every so often models are formulated in terms of turbulent viscosity and conductivity, *K _{m}* and

*K*, rather than directly using the turbulent fluxes. Approximate values of

_{h}*K*and

_{m}*K*can be diagnosed from (A1) and (A2). Assuming steady state again and neglecting vertical energy transport, we multiply (A1) by

_{h}*K*and rewrite

_{m}Note how the expression reverts to the form *K*_{m} ≈ ( *f*_{r}(0)^{2}/*C*_{ε})*l**E*_{k} in near-neutral conditions. Also, the redistribution term reduces *K _{m}* in stably stratified conditions, while it enhances the turbulent viscosity in unstable stratification. We proceed in the same way with (A2) to get

It is interesting to note that the turbulent conductivity does not include temperature variance; thus, even a neutrally stratified fluid is conductive to heat.

## Footnotes

*Corresponding author address:* Dr. Thorsten Mauritsen, Department of Meteorology, Stockholm University, S-106 91 Stockholm, Sweden. Email: thorsten@misu.su.se