## Abstract

Large-eddy simulation (LES) models are widely used to study atmospheric turbulence. The effects of small-scale motions that cannot be resolved need to be modeled by a subfilter-scale (SFS) model. The SFS contribution to the turbulent fluxes is typically significant in the surface layer. This study presents analytical solutions of the classical Smagorinsky SFS turbulent kinetic energy (TKE) model including a buoyancy flux contribution. Both a constant length scale and a stability-dependent one as proposed by Deardorff are considered. Analytical expressions for the mixing functions are derived and Monin–Obukhov similarity relations that are implicitly imposed by the SFS TKE model are diagnosed. For neutral and weakly stable conditions, observations indicate that the turbulent Prandtl number (Pr_{T}) is close to unity. However, based on observations in the convective boundary layer, a lower value for Pr_{T} is often applied in LES models. As a lower Prandtl number promotes a stronger mixing of heat, this may cause excessive mixing, which is quantified from a direct comparison of the mixing function as imposed by the SFS TKE model with empirical fits from field observations. For a strong stability, the diagnosed mixing functions for both momentum and heat are larger than observed. The problem of excessive mixing will be enhanced for anisotropic grids. The findings are also relevant for high-resolution numerical weather prediction models that use a Smagorinsky-type TKE closure.

## 1. Introduction

The large-eddy simulation (LES) modeling technique as proposed by Lilly (1967) is currently widely used to study atmospheric turbulence. In general, LES studies have greatly enhanced our understanding of the turbulence structure of the atmospheric boundary layer (Moeng 1984; Mason 1989; Kosović and Curry 2000; Fedorovich et al. 2004). LES results have provided the means to improve parameterizations of turbulent transport for use in large-scale weather forecast and climate models (e.g., Holtslag and Moeng 1991; Siebesma and Cuijpers 1995; Lock 1998). The inclusion of moist thermodynamics (Deardorff 1980) has allowed LES models to be applied to cloud-topped boundary layers (De Roode et al. 2016) and the ever-increasing computational power has enabled LES of deep convection (Böing et al. 2012; Khairoutdinov and Emanuel 2013).

An LES model solves the filtered budget equations for momentum, heat, moisture, and passive scalars down to scales as fine as the gridmesh size. The transport by eddies with sizes smaller than the grid size, the so-called subgrid- or subfilter-scale (SFS) eddies, needs to be parameterized. Nieuwstadt et al. (1993) compared results as obtained from four different LES models with observations collected in a buoyancy-driven convective boundary layer (CBL). Despite the fact that a rather coarse grid resolution was used, with 160 and 60 m in the horizontal and vertical directions, respectively, a good agreement between the modeling results and observations was found. The simulated boundary layer obtained a final depth of about 1600 m. Since the horizontal size of the dominant eddies is about the same as the boundary layer depth (Young 1987), the SFS contribution to the total fluxes was found to be less than about 10%, except for the surface layer and the entrainment zone near the top of the mixed layer. Furthermore, the results were found to be rather insensitive to details regarding the numerics, with most of the rather small differences between the models being attributable to differences in the SFS scheme.

The applicability of LES models to a wide range of turbulence regimes in the atmosphere enables the nesting of an LES model in a numerical weather prediction (NWP) model (Moeng et al. 2007; Neggers et al. 2012; Schalkwijk et al. 2015) or to operate the NWP as an LES model (Dipankar et al. 2015). Such an approach could arguably be beneficial for the stable boundary layer which remains rather poorly represented in NWP models (Storm et al. 2009; Banta et al. 2013). An accurate prediction of the boundary layer depth and the location of the low-level jet near its top are actually vital for the assessment of the expected wind power from wind turbines (Abkar et al. 2015; Lu and Porté-Agel 2011).

The embedding of LES in an NWP (Schalkwijk et al. 2015) inevitably includes regime transitions associated with the diurnal forcing by solar radiation (Holtslag et al. 2013). In fact, if one aims to simulate the diurnal cycle two opposing constraints emerge. On the one hand, the domain size should be large enough to allow for the development of mesoscale fluctuations of quantities like water vapor that develop during daytime convective conditions (Jonker et al. 1999; De Roode et al. 2004; Honnert et al. 2011) or to permit the formation of deep convective clouds (Böing et al. 2012; Khairoutdinov and Emanuel 2013). On the other hand, the grid size should be sufficiently small to capture the small turbulent eddy sizes typically observed during nocturnal stably stratified conditions (Mahrt 2014). Despite the promising effects of integrating LES in NWP models, the reliability of LES itself in (strongly) stratified conditions is not without debate. The LES results of the stable boundary layer may depend on the specific SFS closure or the grid resolution used (van Stratum and Stevens 2015; Gibbs and Fedorovich 2016). Beare et al. (2006) performed sensitivity experiments with different grid sizes for the first Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS1) stable boundary layer (SBL) case and argued that a grid size of 3.125 m is required to obtain convergence of the modeling results. Sullivan et al. (2016) found that even for a further grid refinement down to a size of 39 cm the simulation results still showed sensitivity to resolution, with most notably the boundary layer depth becoming shallower for a higher resolution. Suppose that computational resources allow us to apply 1000 grid points in each horizontal direction. In that case a grid resolution of about 3 m as proposed by Beare et al. (2006) for the stable boundary layer will then enable a horizontal domain size of 3 × 3 km^{2}. Such a domain is clearly too small to allow for the development of mesoscale fluctuations in the convective boundary layer.

As an illustration, Fig. 1 shows results of the GABLS1 SBL case as obtained from the Dutch Atmospheric LES (DALES) model (Heus et al. 2010) for isotropic grids with varying mesh sizes in the range between 3.125 and 25 m. We find that the SFS turbulent kinetic energy (TKE) depends strongly on the grid size. The runs with a grid size ≥ 6.25 m exhibit a very different behavior of the resolved TKE as compared to the high-resolution run. In the run with the coarsest resolution there is no resolved TKE, and in the runs that used gridmesh sizes of 6.25 and 12.5 m, respectively, the resolved TKE is strongly time dependent, with intermittent periods of almost no resolved turbulence. Figure 2 shows time-mean results during the eighth hour of the simulation. During this period the flow in the simulation that used a gridmesh size of 6.25 m has laminarized, which means that the resolved turbulence has vanished and that the transport is mainly due to parameterized SFS eddy diffusion. By contrast, in the high-resolution run the SFS fluxes dominate the results only in the surface layer. Because the SFS TKE model controls the dissipation of resolved TKE, the disappearance of resolved motions in the coarse-resolution runs strongly hints at the SFS dissipation of resolved TKE becoming too strong.

The SFS TKE model captures the effect of relatively small eddies to the total vertical transport of heat, moisture, and momentum. For well-resolved turbulence this SFS contribution should ideally be smaller than about 10%. However, near the surface, and in the vicinity of thermal inversion layers where the temperature strongly increases across a relatively thin layer, the SFS fluxes may dominate. In general, LES models fail to capture observed flux–gradient relations for the horizontal wind speed in the near-surface layer correctly (Shao et al. 2013; Mason and Thomson 1992; Khanna and Brasseur 1997; Brasseur and Wei 2010; Maronga 2014). These flux–gradient relations are part of a family of functions that relate dimensionless flow properties in the surface layer to a dimensionless height parameter, a description that is widely known as the Monin–Obukhov (MO) similarity theory. Because the lowest few grid levels are always quite strongly affected by the SFS model we will diagnose MO similarity relations that are implicitly imposed by the SFS TKE model. To this end we will consider the special case in which the SFS fluxes dominate the total transport, a situation that closely mimics the flow conditions in the lowest LES model layers near the ground surface. By a comparison with mixing functions that have been derived on the basis of field observations, we will determine the conditions for which excessive SFS mixing will occur. Our philosophy follows Baas et al. (2008), who investigated the scaling behavior of a TKE closure model used in a regional weather forecast model for stably stratified conditions.

We will derive analytical solutions for the Smagorinsky–Lilly-type SFS TKE model, which includes a buoyancy flux, mechanical shear production, and a viscous dissipation term (Lilly 1962; Smagorinsky 1963; Deardorff 1980; Mason 1989). Such SFS TKE models are widely used in LES models (Khairoutdinov and Randall 2003; Heus et al. 2010) and in some NWP models, such as the Weather Research and Forecasting Model (Skamarock et al. 2008; Talbot et al. 2012) or the Met Office Unified Model (Boutle et al. 2014; Efstathiou and Beare 2015). The SFS TKE models may differ in terms of the constants used or the length scale formulation. Our study includes analytical results as obtained for a constant length scale , in addition to a length scale that was designed by Mason and Thomson (1992) to match between the near-surface region and the interior, and a stability-dependent length scale as proposed by Deardorff (1980). The latter length scale is used by Raasch and Schröter (2001) and Sullivan et al. (2003) and has been critically discussed by Schumann (1991) and Gibbs and Fedorovich (2016). It will be explained that and do both depend on the mesh grid size . On the other hand, in the surface layer the relevant length scales of the turbulent eddies are the height above the surface *z* and the Monin–Obukhov stability length scale *L*. Our analysis will show that for a strong stability, the mixing functions for both heat and momentum as imposed by the SFS Smagorinsky TKE model tend to be larger than the ones derived from observations. It will also be argued that if the turbulent Prandtl number applied in the LES model is based on the CBL, then for stable conditions this will yield excessive mixing of heat.

The organization of the paper is as follows. Section 2 summarizes the LES equations, and section 3 presents solutions for the SFS TKE equation for stably stratified conditions. Section 4 gives a brief overview of MO relations and prepares the reader for section 5, which presents mixing functions and MO relations imposed by the SFS TKE model. From a comparison of these analytical results with the empirical MO relations we determine the conditions for which one can expect excessive mixing by the SFS TKE model. The last section summarizes the main findings.

## 2. Formulation of the LES model

Here we will briefly summarize the governing LES equations that apply to an atmosphere free of clouds. For a detailed description of an LES model like DALES we refer the reader to Heus et al. (2010) and Böing (2014), who explain the updated anelastic version for simulations of deep convection. We note that the model is freely available for downloading.

### a. Prognostic budget equations

LES models solve the budget equations for filtered variables including momentum and thermodynamic state variables, such as heat, entropy, or the total water specific humidity. After application of the LES filter the prognostic equation for an arbitrary scalar *φ* can be written as

with *t* the time, and the velocity vector components in the (*x*, *y*, *z*) direction, respectively. In the absence of clouds , with *θ* representing the potential temperature and *q* the water vapor specific humidity. A tilde indicates the filtered mean value and the SFS scalar flux is denoted by .

The Boussinesq form of the filtered momentum equation reads

where *g* is the gravitational acceleration; the reference state virtual potential temperature; the Kronecker-delta function; the virtual potential temperature is defined as

with a thermodynamic constant; *π* is the modified pressure (Deardorff 1973); and an overbar is used to indicate a horizontal slab-mean value. For compact notation, we have included the mean horizontal pressure gradient and the Coriolis force in the source function . The deviatoric part of the SFS momentum flux is computed from Deardorff (1980):

and

The factor that is subtracted in (4) does not arise from the filtering procedure. To compensate it has been added to the filtered pressure term to give the modified pressure. Here and represent the eddy viscosity for momentum and the eddy diffusivity for the thermodynamic scalars, respectively. In a TKE closure approach both are taken proportionally to the square root of the SFS TKE ,

with *λ* the characteristic length scale of the SFS turbulent eddies and and proportionality constants. By analogy with the molecular Prandtl number, which is defined as the ratio of the viscosity to the thermal diffusivity, the ratio can be interpreted as a turbulent SFS Prandtl number,

The budget equation for reads

with a reference density and *p* the pressure. The SFS flux is computed as following (4), and (5) is used to calculate the SFS fluxes of the *θ* and *q*, which in turn are used to calculate the SFS buoyancy flux. The total turbulent transport term is computed following a downgradient diffusion approach,

and the viscous dissipation of *e* by molecular viscosity *ε* is calculated as

with a proportionality constant.

In the remainder of the text we will omit the tildes. With this notation the parameterized equation for the SFS TKE can be written as

with

The classical Smagorinsky model assumes a balance between shear production and dissipation of TKE (Smagorinsky 1963). Stratification effects can be included by maintaining the buoyancy flux (Mason 1989),

This simplified form of the SFS TKE equation thus neglects the tendency, mean advection, and turbulent transport.

### b. Formulations of the length scales

#### 1) Constant length scale

Deardorff (1973) proposed to use the geometric mean of the filter mesh sizes , , and as a representative length scale for SFS eddies,

#### 2) Stability-dependent length scale

Deardorff (1980) argued that for a stable stratification the length scale of the eddies may become smaller than the grid size. The vertical stability can be expressed in terms of the Brunt–Väisälä frequency *N*,

Deardorff proposed the following stability-dependent length scale,

to be used only if its magnitude is smaller than ,

For , the quantity becomes dependent on the stability,

with and . This approach effectively lets the turbulent Prandtl number depend on the stability, with Pr_{T} approaching unity for a very strong stable stratification. The factor is also adapted according to

#### 3) Mason and Thomson length scale

Last we mention the length scale that was constructed by Mason and Thomson (1992) to let the resulting eddy viscosity better match observed MO similarity relations. Specifically, they proposed

with the roughness length. Brown et al. (1994) suggested to use *n* = 2.

## 3. Solutions of the subfilter-scale TKE equation for a stable stratification

Analytical solutions of the SFS TKE equation will be presented for the special case where the production and shear balances the buoyancy flux and the viscous dissipation. Solutions for different length scale formulations will be considered. In the remainder of the paper we will use subscripts , *D*, and *M* to indicate quantities that are derived with the constant length scale defined by (14), the stability-dependent length scale according to (16) and (17), and the length scale following (20), respectively.

### a. Constant length scale

with the gradient Richardson number defined by

SFS turbulence will vanish if exceeds a critical value ,

DALES has evolved from the LES code used by Nieuwstadt et al. (1993), and the original setting is still used; that is, .

From the definition of the eddy viscosity in (6) it follows that

where represents the Smagorinsky constant,

with the filter constant and the Kolmogorov constant (see Table 1 for their values used in DALES). We note that there is no general consensus on the optimum values of these quantities, causing differences in the value for the Smagorinsky constant. For example, Lesieur et al. (2005) uses , whereas Schumann (1975) and Meneveau and Katz (2000) use a value of 1.5 and 1.6, respectively. As compared to DALES, Mason (1989) uses a smaller filter constant of . Kleissl et al. (2003) used an array of sonic anemometers to measure SFS diffusion constants in the atmospheric surface layer. They actually found that is not constant but is reduced near the ground surface and also tends to become smaller with increasing stability.

### b. Stability-dependent length scale

Also for it is possible to determine a solution for *e* (VanZanten 2000). Substitution of (6) in (13), and using (16), (18), and (19), gives

In this case *e* becomes zero for a critical Richardson number of

DALES uses the constants as proposed by Deardorff (1980), which yields (see Tables 1 and 2 for a summary). This value is slightly smaller than the critical Richardson number for , .

With aid of (16), (17), and (26) we can assess the threshold Richardson number above which value is applied (VanZanten 2000),

where for the settings used in DALES.

With aid of the solution for *e* in (26), we can express in terms of the gradient Richardson number,

Interestingly, the stability-dependent length scale still depends on the grid size .

### c. Mason and Thomson length scale

Mason and Thomson (1992) proposed to replace the factor in (24) by to give

They recommended a turbulent Prandtl number of 0.7.

## 4. Eddy viscosity following Monin–Obukhov similarity relations

We will derive the respective eddy viscosities and diffusivities from empirical MO relations for comparison with the analytical SFS TKE results. For the convenience to the reader, some key definitions used in MO similarity theory are summarized in the appendix. In this study we will frequently use the following key relations:

Högström (1988) evaluated several formulas for and that were based on observations from field campaigns. In addition, Högström proposed modifications to correct for instrumental shortcomings (his Tables VI and VII). The parameterizations for a stable stratification as obtained by Businger et al. (1971) from the Kansas field experiment read

with corrected coefficients , , , and . We have introduced the subscript “obs” to denote a result that is based on observations. The corrected fits for Dyer (1974) have slightly different slope factors: and . Högström used and a reciprocal turbulent Prandtl number for neutral conditions. With regard to the latter choice, Högström remarked that from a physical point of view it is difficult to find a reason why for neutral conditions should differ from . Besides, he noted that for the uncertainty in the observations could well accommodate . The large scatter in the observations for stable conditions and, particularly, the scarce amount of data for hindered to select the best fits, although Högström (1988) hinted at a slightly better performance for the modified of Businger et al. (1971) and the modified of Dyer (1974).

Grachev et al. (2007b) analyzed measurements collected during stable and very stable conditions over the Arctic sea ice pack during the Surface Heat Budget of the Arctic Ocean Experiment (SHEBA). They proposed the following nonlinear MO similarity functions:

The SHEBA functions therefore suggest a turbulent Prandtl number of unity for neutral conditions. Moreover, the three parameterizations presented above have in common that for a stable stratification , which implies that the turbulent Prandtl number depends on the stability. We note that direct numerical simulation results analyzed by Van de Wiel et al. (2008) suggest a value of about 5 both for and , whereas in a similar study by Ansorge and Mellado (2014) a slightly larger value of 5.7 was reported for .

With aid of (A10), it is possible to derive the corresponding eddy viscosity for MO functions that are linear in like (34). To this end we will first eliminate to give

In the surface layer we will approximate . Next we can rewrite (A6) as a function of ,

and can be obtained straightforwardly from (33),

## 5. Comparison of the SFS TKE eddy viscosity with Monin–Obukhov similarity relations

Here we will compare the eddy mixing functions as imposed by the SFS TKE scheme with those obeying MO similarity relations. Our analysis specifically applies to situations in which the SFS contribution dominates the total flux (i.e., in the surface layer). For simplicity, we will consider the situation in which the resolved vertical velocity is negligibly small, which allows to write .

To facilitate a direct comparison of the eddy mixing coefficient that follows from the SFS TKE model and MO functions, we express the eddy viscosity and diffusivity as (Duynkerke and De Roode 2001)

with and the mixing functions for momentum and heat, respectively. From a comparison of (39) with the definitions (A6) and (A11), it follows that

As schematically depicted in Fig. 3, this implies that for any , and for any arbitrary form of and , the mixing functions can be determined, while (A10) allows to compute the corresponding stability in terms of .

Table 3 summarizes the analytical expressions for the eddy viscosity, the mixing function, and the critical Richardson number for the length scales and , in addition to the ones obtained from the MO relations.

### a. Mixing function

For the special case of MO functions that are linear in [see (34)] the mixing function can be obtained from (37) with (39),

A comparison of (24) or (30) with (39) shows that the mixing functions for and can be written as, respectively,

We note that the maximum values for the mixing function are obtained at the lowest model level, , and that the magnitude of will decrease proportionally to the reciprocal of the height *z*.

Figure 4 shows that the mixing function evaluated at a height increases for increasing aspect ratio . Here denotes the horizontal grid size, which we take equal in the two horizontal directions. Except for a weak stable stratification, excessive SFS mixing will occur if . Grid size configurations of are commonly used to study convective cases. Our analysis suggests that one should be careful to apply an anisotropic grid, in particular if one aims to simulate the diurnal cycle including a nocturnal period of stable stratification. The results actually demonstrate that excessive mixing will occur for . It is possible to diminish the mixing function by reducing the value of the Smagorinsky constant. For example, Beare and Macvean (2004) reduce the Smagorinsky constant from 0.23 to 0.15 for coarser-grained simulations of grid length 10 m or more.

The stability-dependent length scale somewhat alleviates the problem of excessive mixing as compared to the use of , yet for the SFS TKE model also produces too much SFS mixing (see Fig. 4b).

An inspection of the imposed mixing function for the length scale in (20), which was proposed to closely match MO relations near the surface, shows that for an isotropic grid excessive mixing will occur if exceeds a value of about 0.11 (see Fig. 4c). We also can see that for a larger value of the horizontal grid size the range of where excessive mixing takes place expands toward smaller values.

### b. Mixing function and the turbulent Prandtl number

In the model intercomparison study of the dry convective boundary layer reported by Nieuwstadt et al. (1993) the four participating LES models applied turbulent Prandtl numbers in the range between 0.33 and 0.46. Other LES models use a constant value of 0.7 (Mirocha et al. 2010), 0.4 (Savre et al. 2014), or 1/3 (Matheou et al. 2011), respectively. However, these turbulent Prandtl numbers are clearly not very representative for stably stratified conditions, as is illustrated in Fig. 5. The “default value” appears to be much smaller than . Interestingly, as can be seen from the definition of for (18), for a strong stability , and consequently in this limit .

For the mixing function for heat can be simply expressed as

This result, for , is shown in Fig. 6a. It is striking to find that at the used turbulent Prandtl number pushes toward excessive mixing for any positive value of . Because this is highly undesired, this result provides another strong argument to use for the stable regime. Excessive mixing is also found for the stability regime where the stability-dependent length scale is applied (see Fig. 6b).

Mason and Thomson (1992) proposed a turbulent Prandtl number of 0.7. It can be directly seen from Fig. 4c that this value will yield excessive mixing for heat for even smaller values of as compared to the mixing of momentum.

### c. MO similarity functions imposed by the subfilter-scale TKE model

Given the expressions for the mixing functions it is straightforward to diagnose the relation imposed by the SFS TKE equation in terms of . From (A6) and (39) we can write

With use of (A10) the dimensionless height can also be expressed as a function of ,

For any we can therefore diagnose both and , which provides another means to compare the SFS TKE scaling behavior with MO relations.

Figure 7 shows the MO function imposed by the SFS TKE model. For the lowest model level we find that the slope of is much smaller than the observed value for . For a neutral stratification the observed relation is generally not obeyed. For an isotropic grid, the value is too large and becomes too small toward a stronger anisotropy of the grid. In the high range of the use of gives a slope that is much closer to the observations.

We notice that for the results look satisfactorily with regard to the observations; however, since it is only applied for a sufficiently strong stable stratification, we find that the stability regime in which is actually applied gives too high values for . The approximately linear behavior for the high-stability regime can be understood by taking the limit , which according to (18) and (19) yields constant values for and , respectively. In that case the SFS TKE equation [(13)] simply becomes

Using the relation (A10) we can obtain the following analytical solution:

This solution is similar to the one found by for the Reynolds-averaged TKE equation with stability-dependent length scale as applied in the Regional Atmospheric Climate Model (RACMO) (Baas et al. 2008). They discussed in detail how the coefficients of the TKE model can be chosen in order to match the observed relation. According to the numbers presented in Table 1 the proportionality factor .

## 6. Large-eddy simulations of the GABLS1 case

As an illustration of our analysis we will discuss four simulations of the GABLS1 case (Beare et al. 2006), which all used an isotropic grid but with different mesh sizes of 3.125, 6.25, 12.5, and 25 m, respectively. The simulations shown all used and the full SFS TKE model including the tendency and turbulent transport terms. Following the GABLS1 case description, each simulation lasted 9 h and was run with 128 grid points in each direction. Here we note that DALES applies the MO relations as proposed by Beljaars and Holtslag [1991, their (28) and (32), respectively] to compute the fluxes at the ground surface from the local gradients of momentum and heat.

### Sensitivity to grid resolution and length scale

Figure 8a shows that for the height of the maximum wind speed is located at about 175 m, in agreement with the high-resolution simulation results presented by Beare et al. (2006). Figure 9 shows the SFS TKE results as well as its analytical solution in (21) for the simulations in which there was hardly any resolved turbulence, *l*_{Δ} = 12.5 and 25 m, respectively. For the latter two cases, we find a very good correspondence between the LES results and the analytic solution. Since we performed the runs with the full prognostic equation for *e* in (8), including the SFS turbulent transport as well as the tendency of *e*, we conclude that for SFS-dominated SBL flow the prognostic SFS TKE model equation [(8)] almost behaves identically to the Smagorinsky model with a buoyancy flux contribution. According to Nieuwstadt (1984), the total transport of TKE for the quasi-steady stable boundary layer can be assumed to be negligibly small with respect to the other terms. For the laminarized cases, Fig. 8c shows that is rather small in the lower part of the boundary layer, which prevents the application of because is smaller than the critical value of 0.075.

Following Beare et al. (2006), Fig. 10 compares the MO similarity functions as obtained from the LES results with the observed MO relations. For the high-resolution run, the slopes for both and are about , resulting in a turbulent Prandtl number of about unity. However, the effect of the prescribed ratio is clearly visible from the diagnosed ratio , which tends toward this value for small values of . This result clearly demonstrates the effect of a SFS TKE-dominated regime, which for the high-resolution run applies to the surface layer, and is in accord with a Prandtl number as predicted by (7). For the SFS-dominated solutions the slope of is too large and too small for as compared to the observations.

## 7. Conclusions

We have studied the Smagorinsky SFS TKE equation including shear, a buoyancy flux, and a dissipation term. More specifically, we derive analytical solutions for the special case for which there is no resolved turbulence. Three different length-scale formulations are considered. The simplest one is equal to the geometric mean value of the gridmesh sizes in all three directions . The length scale as proposed by Deardorff (1980) depends on the magnitude of the SFS TKE and the static stability. A third length scale involves the size of the turbulent eddies near the surface for neutral conditions (Mason and Thomson 1992). The results are presented in terms of analytical expressions for the mixing functions, which in turn are used to diagnose MO relations that are implicitly imposed by the SFS TKE model. For each length scale, their imposed MO relations are found to be strongly dependent on the gridmesh size. The length scales that are actually used in MO flux–gradient relations are the height above the surface and the MO stability length. For this reason the SFS fluxes will generally not match the observed flux–gradient relations that are applied as a lower boundary condition at the ground surface. In any case, the SFS TKE model should not produce excessive mixing. From a direct comparison of the diagnosed SFS mixing functions with empirical fits from field experiments, we have assessed the conditions that give rise to this undesirable situation. We will summarize our findings for the neutral or weakly stable regime and the very stable regime, respectively.

### a. Neutral or weakly stable stratification,

The widely used value of 1/3 for the turbulent Prandtl number will inevitably yield excessive mixing of heat at the lowest model levels. This finding applies to any of the three length scales studied, although it should be kept in mind that for , is not used. The mixing function for heat can be reduced if is increased. This is actually supported by observations which show that the value of is close to unity. In fact, if one sets in the SFS TKE model, the problem of excessive mixing of heat will vanish. For an isotropic grid, the LES SFS imposed mixing function for momentum will always be smaller than the MO-based mixing function. To match the turbulent Prandtl number for different stabilities, Gibbs and Fedorovich (2016) present a parameterization that allows to vary smoothly between 1/3 for unstable conditions and 1 for stable conditions.

### b. Very stable stratification,

The spread in the observations of the turbulent Prandtl number does not permit to recommend a turbulent Prandtl number (Sukoriansky et al. 2005). The studies by Ohya (2001), Zilitinkevich et al. (2007), and Anderson (2009) report that the turbulent Prandtl number increases with increasing stability, which is hypothesized to be due to internal wave activity that mix the momentum but do not mix a scalar. By contrast, Grachev et al. (2007a) and Sorbjan (2010) suggest that if outliers are rejected, decreases from 0.9 in nearly neutral conditions to 0.7 for increasing stability. In any case, in this regime excessive SFS mixing of momentum and heat will occur both for and . Because for increasing stability diminishes and increases, the mixing functions for this length scale are much closer to the ones based on MO relations.

Field observations analyzed by Kleissl et al. (2003) suggest that the Smagorinsky constant is reduced near the ground surface and also tends to become smaller with increasing stability. A part of the excessive mixing in the surface layer might therefore be attributed to the prescribed constant value of .

### c. Anisotropic grids

Because depends on the horizontal grid sizes, the application of an anisotropic grid with will yield enhanced mixing functions for both heat and momentum. This is true for each of the length scales considered here because they all include some dependency on . To avoid excessive SFS mixing of both momentum and heat, the use of an anisotropic grid should therefore be strongly discouraged for a stable stratification.

### d. NWP models

For we derived an analytical expression for its imposed function in the limit of a very strong stable stratification. This solution can actually be extended toward cases of weak stability by setting the constants and to zero. In this way, might be calibrated to match an observed fit. As discussed by Baas et al. (2008), such an approach might be useful for NWPs that use a Smagorinsky type of TKE closure.

### e. Variants of the Smagorinsky-type SFS TKE model

The central question addressed in this study is under which conditions, in terms of atmospheric stability and grid configuration, excessive SFS mixing will take place. Although our analysis is dedicated to the traditional Smagorinsky SFS TKE model, this question is actually relevant to any arbitrary SFS TKE model. In this context we would like to mention two variants of the Smagorinsky model. In the SFS model by Sullivan et al. (1994), the SFS TKE production is strongly diminished by taking out the mean shear from the shear production term. Their LES results are found to be insensitive to grid anisotropy. Basu and Porté-Agel (2006) applied a dynamic SFS modeling approach in which the Smagorinsky constant and the turbulent Prandtl number are diagnosed from the resolved flow at each grid point and for each time step. They simulated the GABLS1 case and concluded that their results showed relatively little resolution dependency, even for relatively coarse resolutions. Last, in the context of the present study we would like to mention the work of Shao et al. (2013), who presented a method that aims to give SFS fluxes of heat and moisture in the surface layer that are consistent with the surface fluxes.

## Acknowledgments

We would like to kindly thank Erwin de Beus for his technical assistance with DALES and Dr. Arnold Moene for improving the implementation of the lower boundary condition for the subfilter TKE. The comments and suggestions by three anonymous are very much appreciated and were helpful to improve the manuscript. Dr. van de Wiel received financial support from a Consolidator Grant of the European Research Council (648666), which is gratefully acknowledged.

### APPENDIX

#### Summary of Definitions Used in Monin–Obukhov Similarity Theory

The MO functions and describe the relations between the mean vertical gradients of the horizontal wind speed and the virtual potential temperature, respectively, in terms of their vertical fluxes according to

with *κ* the von Kármán constant, the local friction velocity,

and defines the absolute value of the mean shear of the horizontal wind velocity,

We note that in practice the turbulent fluxes in these relations are usually obtained from point measurements. Following Reynolds decomposition, an overbar indicates a time-mean value and a prime a fluctuation with respect to the mean. The dependencies of and on the stability need to be determined empirically from observations. They are generally expressed in terms of the nondimensional group , with the Obukhov length *L* defined by

Here the flux values of momentum and heat are taken at the surface. However, Nieuwstadt (1984) analyzed observations collected during stable conditions and hypothesized that generally the local, height-dependent, MO length is a more appropriate length scale. The difference between *L* and is that for the latter quantity the fluxes observed at a height *z* above the surface are used. Formally, one can show that MO similarity is an asymptotic case of local similarity (Van de Wiel et al. 2012), which is used in the present study. In the surface layer, which is defined as the lower 10% of the boundary layer, the difference between the local and the surface-based MO lengths is less than 10%. The benefit of using is that it allows to extend MO similarity to the whole stable boundary layer (Nieuwstadt 1984; Wyngaard 2010)—a strategy that was also utilized in the GABLS1 model intercomparison study (Beare et al. 2006).

With some simple manipulations, it will be possible to directly compare this expression with the analytic solutions for the eddy viscosity as found for the SFS TKE model. More specifically, as schematically depicted in Fig. 3, we wish to present the MO-based eddy viscosity as a function of . To this end we will first substitute the definition of in (A1) to obtain

to eliminate and from (A7),

Last, with aid of (A2), (A6), (A8), and (A9), we can express the eddy diffusivity in terms of MO relations

which in turn allows us to write the turbulent Prandtl number [see (7)] as Grachev et al. (2007b),

## REFERENCES

*Phys. Fluids*,

**22**, 021303, doi:.

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

*Topics in Micrometeorology: A Festschrift for Arch Dyer*, B. B. Hicks, Ed., Springer, 55–78.

*Large-Eddy Simulations of Turbulence*. Cambridge University Press, 232 pp.

*Proc. IBM Scientific Computing Symp. on Environmental Sciences*, Yorktown Heights, NY, IBM, 195–209.

*Phys. Fluids*,

**23**, 065101, doi:.

*Turbulent Shear Flows 8*, F. Durst et al., Eds., Springer, 343–367, doi:.

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

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).