**3.1. Standard eddy-viscosity/diffusivity models**

Base on the Boussinesq hypothesis [14], a common class of SGS models, eddy-viscosity/diffusivity model, parameterizes the SGS stress' deviatoric part and the SGS flux as

$$
\pi\_{\rm ij} - \frac{1}{3} \delta\_{\rm ij} \pi\_{\rm kk} = -2 \nu\_{\rm sgs} \tilde{\mathbf{S}}\_{\rm ij} \qquad \text{and} \qquad q\_{\rm i} = -\frac{\nu\_{\rm sgs}}{\rm Sc\_{\rm sgs}} \frac{\partial \tilde{\theta}}{\partial \mathbf{x}\_{\rm i}^{\prime}} \tag{5}
$$

where *<sup>S</sup>ij* <sup>=</sup> <sup>1</sup> 2 *∂ui ∂xj* <sup>+</sup> *<sup>∂</sup><sup>u</sup><sup>j</sup> ∂xi* is the resolved (filtered) strain rate tensor, *νsgs* is the SGS eddy viscosity and *Scsgs* is the SGS Schmidt number. Several different models have been used to determine the eddy viscosity. The most common one was introduced by Smagorinsky [82] by assuming a local equilibrium between production and dissipation of SGS kinetic energy. The Smagorinsky model computes the eddy viscosity as

$$\nu\_{\rm sgs} = \left(\mathbb{C}\_{\rm S} \widetilde{\boldsymbol{\Delta}}\right)^2 \left|\widetilde{\boldsymbol{S}}\right| \,\tag{6}$$

where *S* = <sup>2</sup>*<sup>S</sup>ijSij*<sup>1</sup> 2 is the strain rate, and *CS* is a non-dimensional parameter called Smagorinsky coefficient. In isotropic turbulence, if a cut-off filter is used in the inertial subrange and the filter scale <sup>∆</sup> is equal to the grid size, then *Cs* <sup>≈</sup> 0.17 and *Scsgs* <sup>≈</sup> 0.5 [6, 53]. However, flow anisotropy, particularly the presence of a strong mean shear near have found the range of *Scsgs* is from 0.33 to 0.7.

the surface in high-Reynolds-number ABLs, makes the optimum values of those coefficients depart from their isotropic counterparts [e.g., 13, 46, 77]. [45] have also shown that the optimal value of the Smagorinsky coefficient decreases with increasing atmospheric stability in order to account for the reduction of characteristic length scales associated with thermal stratification. A common practice is to specify the coefficients in an *ad-hoc* fashion. The *ad-hoc* damping function proposed by Mason and Thomson [65] can be rewritten [76] as: *Cs* = *<sup>C</sup>*−*<sup>n</sup>* <sup>0</sup> + *κ z* <sup>∆</sup> <sup>+</sup> *<sup>z</sup>*<sup>0</sup> ∆ <sup>−</sup>*<sup>n</sup>* −1/*<sup>n</sup>* , where *n* is an adjustable parameter, and studies [e.g., 5, 65, 76] have reported that this formulation with values of *C*<sup>0</sup> ranging from 0.1 to 0.3, and *n* = 1, 2, or 3 can deliver a more realistic logarithmic velocity profile in the surface layer than 10.5772/57051

195

http://dx.doi.org/10.5772/57051

*(2)* For complex flows, it may not be possible to find a universal coefficient that is appropriate for the entire domain at all times. For instance, using the theoretical values in LES of boundary-layer turbulence, the Smagorinsky model yields too much transfer of energy and scalar variance from resolved to subgrid scales near the ground [64, 76]. This excessive SGS dissipation leads to unrealistic mean velocity and scalar profiles, with excessive non-dimensional mean shear and scalar gradients (Fig. 1(a)). The overestimation of the SGS dissipation leads also to velocity spectra (Fig. 1(b)) that decay too rapidly in the near-ground region [76]. Note that the normalized spectra do not show the collapse as found

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

*(3)* In anisotropic turbulence there can be a net transfer of kinetic energy from small to large scales [17, 84, 85]. The Smagorinsky model is by construction dissipative. For numerical stability, this is a desirable characteristic of the model. However, the actual SGS stress may also provide opportunities for inverse energy transfer. Studies of Khanna and Brasseur [43], Juneja and Brasseur [39], and Porté-Agel et al. [76] have showed that on coarse grids the Smagorinsky model may induce large errors because it is not able to account for the strong flow anisotropy in the near-wall region. Moreover, eddy-viscosity models do not have the same rotation transformation properties as the actual SGS stress tensor, which is not material frame indifferent (MFI). Recent studies [33, 47, 61, 62] revisited the importance of the MFI-consistency of the modeled SGS stresses. In LES of mesoscale and large-scale atmospheric turbulence including planetary rotation, the Smagorinsky model induces extra errors, yields excessive dissipation at large scales, and thus delivers unsatisfactory results, such as the failure of capturing large-scale cyclone/anti-cyclone asymmetry in favor of

To overcome some critical weaknesses of eddy-viscosity/diffusivity models, such as low correlations between the exact SGS terms and the modeling terms, modeling errors regarding the strong flow anisotropies, and a common issue of the early eddy-viscosity/diffusivity models that, in the context of ABL flows, the mean wind and temperature profiles in the surface layer differ from the Monin-Obukhov similarity forms [e.g., 15, 88] as shown in figure 1(a), a nonlinear SGS approach has been introduced and tested [57, 58, 60]. The new closure is based on gradient models, which, different than eddy-viscosity/diffusivity models, are derived from the Taylor series expansions of the SGS terms that appear in the filtered conservation equations, do not locally assume the same eddy viscosity/diffusivity for all directions, and make no use of prior knowledge of the interactions between the resolved motions and the SGS motions. At *a-priori* level, gradient models generally predict the structure of the exact SGS terms much more accurately than eddy-viscosity/diffusivity models (and therefore are better able to capture anisotropic effects and disequilibrium, e.g., 32, 56, 61, 62, 77). These features make gradient models attractive. However, when implemented in LESs, they are not able to produce the correct levels of the SGS energy production (energy transfer between resolved and SGS scales), and as a result, simulations

often become numerically unstable as reported in a variety of contexts [e.g., 80].

Many schemes have been proposed for resolving this insufficient-dissipation issue relative to gradient models. Bardina et al. [8] proposed a mixed procedure; and later on, Vreman et al. [94, 95] achieved the mixing of the gradient form with eddy-viscosities, and showed

in experimental studies [40, 73].

**3.2. A new nonlinear formulation**

cyclone [62].

does the standard Smagorinsky model using a constant coefficient. Studies [e.g., 5, 63, 65]

**Figure 1.** ((a) Vertical distribution of non-dimensional gradient of the mean streamwise velocity (Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup> u*∗ *d* <*u*> *d z* ) and scalar concentration (Φ*<sup>θ</sup>* = *<sup>κ</sup><sup>z</sup> θ*∗ *∂*<*θ* > *<sup>∂</sup><sup>z</sup>* ) obtained using standard eddy-viscosity/diffusivity models ( shows <sup>Φ</sup>*<sup>M</sup>* using *<sup>C</sup>*<sup>0</sup> <sup>=</sup> 0.17 and *n* = 1, △ shows Φ*<sup>θ</sup>* using *Scsgs* = 0.5, and � shows Φ*<sup>θ</sup>* using *Scsgs* = 0.7, figure is modified from[60]); the dashed line corresponds to the classical similarity profile; (b) averaged non-dimensional 1-D spectra of the streamwise velocity obtained using the traditional Smagorinsky model with the constant coefficient *Cs* = 0.17.

Despite the popularity of the eddy-viscosity/diffusivity models, they are known to have important limitations:

*(1)* The eddy-viscosity/diffusivity closure assumes a one-to-one correlation between the exact SGS terms and the eddy-viscosity/diffusivity term, and locally employs the same eddy-viscosity/diffusivity for all directions. *A-priori* studies using experimental data and DNS data have shown that eddy-viscosity SGS models cannot capture the local characteristics of the SGS fluxes [21, 29, 32, 56, 61, 68, 69, 81]. For example, the underlying assumption of strain rates being aligned with the SGS stress tensor is found to be unrealistic (see 32 and the references therein). Field studies have also shown that the correlation between measured (in the field) and modeled SGS stresses (using the eddy-viscosity model) is very low (about 20 %) [75, 77]. Moreover, the fully dissipative nature of the model precludes backscatter, or negative SGS dissipation (energy transfer from unresolved to resolved scales), which has been reported in *a-priori* studies.

*(2)* For complex flows, it may not be possible to find a universal coefficient that is appropriate for the entire domain at all times. For instance, using the theoretical values in LES of boundary-layer turbulence, the Smagorinsky model yields too much transfer of energy and scalar variance from resolved to subgrid scales near the ground [64, 76]. This excessive SGS dissipation leads to unrealistic mean velocity and scalar profiles, with excessive non-dimensional mean shear and scalar gradients (Fig. 1(a)). The overestimation of the SGS dissipation leads also to velocity spectra (Fig. 1(b)) that decay too rapidly in the near-ground region [76]. Note that the normalized spectra do not show the collapse as found in experimental studies [40, 73].

*(3)* In anisotropic turbulence there can be a net transfer of kinetic energy from small to large scales [17, 84, 85]. The Smagorinsky model is by construction dissipative. For numerical stability, this is a desirable characteristic of the model. However, the actual SGS stress may also provide opportunities for inverse energy transfer. Studies of Khanna and Brasseur [43], Juneja and Brasseur [39], and Porté-Agel et al. [76] have showed that on coarse grids the Smagorinsky model may induce large errors because it is not able to account for the strong flow anisotropy in the near-wall region. Moreover, eddy-viscosity models do not have the same rotation transformation properties as the actual SGS stress tensor, which is not material frame indifferent (MFI). Recent studies [33, 47, 61, 62] revisited the importance of the MFI-consistency of the modeled SGS stresses. In LES of mesoscale and large-scale atmospheric turbulence including planetary rotation, the Smagorinsky model induces extra errors, yields excessive dissipation at large scales, and thus delivers unsatisfactory results, such as the failure of capturing large-scale cyclone/anti-cyclone asymmetry in favor of cyclone [62].

### **3.2. A new nonlinear formulation**

4 Computational and Numerical Simulations

 *κ z* <sup>∆</sup> <sup>+</sup> *<sup>z</sup>*<sup>0</sup> ∆ <sup>−</sup>*<sup>n</sup>*

have found the range of *Scsgs* is from 0.33 to 0.7.

(a) (b)

using the traditional Smagorinsky model with the constant coefficient *Cs* = 0.17.

−1/*<sup>n</sup>*

θ

*Cs* =

 *<sup>C</sup>*−*<sup>n</sup>* <sup>0</sup> +

concentration (Φ*<sup>θ</sup>* = *<sup>κ</sup><sup>z</sup>*

important limitations:

*θ*∗ *∂*<*θ* >

been reported in *a-priori* studies.

the surface in high-Reynolds-number ABLs, makes the optimum values of those coefficients depart from their isotropic counterparts [e.g., 13, 46, 77]. [45] have also shown that the optimal value of the Smagorinsky coefficient decreases with increasing atmospheric stability in order to account for the reduction of characteristic length scales associated with thermal stratification. A common practice is to specify the coefficients in an *ad-hoc* fashion. The *ad-hoc* damping function proposed by Mason and Thomson [65] can be rewritten [76] as:

5, 65, 76] have reported that this formulation with values of *C*<sup>0</sup> ranging from 0.1 to 0.3, and *n* = 1, 2, or 3 can deliver a more realistic logarithmic velocity profile in the surface layer than does the standard Smagorinsky model using a constant coefficient. Studies [e.g., 5, 63, 65]

10−1

*<sup>∂</sup><sup>z</sup>* ) obtained using standard eddy-viscosity/diffusivity models ( shows <sup>Φ</sup>*<sup>M</sup>* using *<sup>C</sup>*<sup>0</sup> <sup>=</sup> 0.17 and

100

Eu(z)u−2

**Figure 1.** ((a) Vertical distribution of non-dimensional gradient of the mean streamwise velocity (Φ*<sup>M</sup>* = *<sup>κ</sup><sup>z</sup>*

*n* = 1, △ shows Φ*<sup>θ</sup>* using *Scsgs* = 0.5, and � shows Φ*<sup>θ</sup>* using *Scsgs* = 0.7, figure is modified from[60]); the dashed line corresponds to the classical similarity profile; (b) averaged non-dimensional 1-D spectra of the streamwise velocity obtained

Despite the popularity of the eddy-viscosity/diffusivity models, they are known to have

*(1)* The eddy-viscosity/diffusivity closure assumes a one-to-one correlation between the exact SGS terms and the eddy-viscosity/diffusivity term, and locally employs the same eddy-viscosity/diffusivity for all directions. *A-priori* studies using experimental data and DNS data have shown that eddy-viscosity SGS models cannot capture the local characteristics of the SGS fluxes [21, 29, 32, 56, 61, 68, 69, 81]. For example, the underlying assumption of strain rates being aligned with the SGS stress tensor is found to be unrealistic (see 32 and the references therein). Field studies have also shown that the correlation between measured (in the field) and modeled SGS stresses (using the eddy-viscosity model) is very low (about 20 %) [75, 77]. Moreover, the fully dissipative nature of the model precludes backscatter, or negative SGS dissipation (energy transfer from unresolved to resolved scales), which has

∗ z−1

101

, where *n* is an adjustable parameter, and studies [e.g.,

10−1 100

−1

k1 z

*u*∗ *d* <*u*>

−5/3

*d z* ) and scalar

To overcome some critical weaknesses of eddy-viscosity/diffusivity models, such as low correlations between the exact SGS terms and the modeling terms, modeling errors regarding the strong flow anisotropies, and a common issue of the early eddy-viscosity/diffusivity models that, in the context of ABL flows, the mean wind and temperature profiles in the surface layer differ from the Monin-Obukhov similarity forms [e.g., 15, 88] as shown in figure 1(a), a nonlinear SGS approach has been introduced and tested [57, 58, 60]. The new closure is based on gradient models, which, different than eddy-viscosity/diffusivity models, are derived from the Taylor series expansions of the SGS terms that appear in the filtered conservation equations, do not locally assume the same eddy viscosity/diffusivity for all directions, and make no use of prior knowledge of the interactions between the resolved motions and the SGS motions. At *a-priori* level, gradient models generally predict the structure of the exact SGS terms much more accurately than eddy-viscosity/diffusivity models (and therefore are better able to capture anisotropic effects and disequilibrium, e.g., 32, 56, 61, 62, 77). These features make gradient models attractive. However, when implemented in LESs, they are not able to produce the correct levels of the SGS energy production (energy transfer between resolved and SGS scales), and as a result, simulations often become numerically unstable as reported in a variety of contexts [e.g., 80].

Many schemes have been proposed for resolving this insufficient-dissipation issue relative to gradient models. Bardina et al. [8] proposed a mixed procedure; and later on, Vreman et al. [94, 95] achieved the mixing of the gradient form with eddy-viscosities, and showed that mixed gradient models can capture disequilibrium. In LES of rotating turbulence, to overcome the weakness that the traditional eddy viscosity (−2*νsgsSij*) is too dissipative at large scales, which may hinder kinetic energy from transferring to large scales, Lu et al. [62] considered a hyper-viscosity term (+*νu*∇2*<sup>S</sup>ij*) as suggested by previous researchers [9, 26]. The new mixed nonlinear model is defined as

$$\pi\_{\rm ij} = 2k\_{\rm sgs} \left( \frac{G\_{\rm ij}}{G\_{\rm kk}} \right) + \nu\_{\rm u} \nabla^2 \widetilde{S}\_{\rm ij} \, , \tag{7}$$

norm <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> <sup>=</sup>

 *G*2 *<sup>θ</sup>*,1 <sup>+</sup> *<sup>G</sup>*<sup>2</sup>

of the SGS kinetic energy, *usgs* = *C*

variance production, *<sup>P</sup><sup>θ</sup>* = −*qi*

<sup>−</sup> *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> *∂θ ∂xi*

the SGS kinetic energy production *<sup>P</sup>* (*<sup>P</sup>* = −*τij*

*k*3/2 *sgs*

kinetic energy, *ksgs* = <sup>1</sup>

dissipation is *ε* = *C<sup>ε</sup>*

this leads to

*<sup>θ</sup>sgs* <sup>=</sup> <sup>∆</sup>

*Cεθ*

*<sup>θ</sup>*,2 <sup>+</sup> *<sup>G</sup>*<sup>2</sup>

*ksgs* <sup>=</sup> *<sup>H</sup>* (*P*) <sup>4</sup><sup>∆</sup><sup>2</sup>

classical evaluation of the SGS scalar variance dissipation rate is *εθ* = *Cεθ*

*∂θ ∂xi*

and the dissipation rate are always non-negative *<sup>θ</sup>sgs* <sup>=</sup> *<sup>H</sup>* (*P<sup>θ</sup>* ) <sup>∆</sup>

time is proportional to the SGS turbulent characteristic time, *<sup>θ</sup>*<sup>2</sup>

following equation for the magnitude of the SGS flux

Also, the coefficient can be derived as *Cεθ* = <sup>1</sup>

*usgs* <sup>=</sup> (*uiui* <sup>−</sup> *<sup>u</sup>iu<sup>i</sup>*) <sup>=</sup> <sup>2</sup> *ksgs*), one obtains *<sup>C</sup>εθ* <sup>=</sup> 1.0.

<sup>|</sup>**q**<sup>|</sup> <sup>=</sup> *<sup>H</sup>* (*P<sup>θ</sup>* ) *<sup>H</sup>* (*P*) <sup>1</sup>

the scalar field) leads to satisfactory results. Using the model *εθ* = <sup>1</sup>

*Sc*

<sup>2</sup><sup>∆</sup><sup>2</sup> *C*2 *εθ* 10.5772/57051

197

http://dx.doi.org/10.5772/57051

<sup>+</sup> *<sup>∂</sup><sup>u</sup><sup>j</sup> ∂xi*

 is

*<sup>θ</sup>*,3. To close the approach, it is required to evaluate the SGS

*ksgs*. The value of *ksgs* is evaluated by using the resolved

2 *∂ui ∂xj*

, (10)

*θ*2 *sgsusgs*

*<sup>ε</sup>* . Tests have shown

*ksgs* , one obtains the

. (11)

<sup>∆</sup> . Using

. Although a

*<sup>∂</sup>xj* <sup>=</sup> <sup>−</sup>*τijSij*, where *<sup>S</sup>ij* <sup>=</sup> <sup>1</sup>

, and the SGS scalar variance dissipation rate *εθ* . A

*Cεθ*

*sgs εθ* <sup>∝</sup> *ksgs*

<sup>−</sup> *<sup>G</sup>ij Gkk* <sup>−</sup> *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> *∂θ ∂xi*

*Sc θ*2 *sgsε*

*Sij*

*<sup>C</sup>* . When adopting (as this study has

<sup>2</sup> *<sup>τ</sup>ii*, and the magnitude of the SGS flux vector, <sup>|</sup>**q**|. Even though a

Large-eddy simulation of turbulent flows with applications to atmospheric boundary layer research

<sup>∆</sup> . Simulations allow for no local energy transfer from unresolved to

*Sij*2

. The clipping procedure assumes the SGS scalar variance production

previous approach [20] places much emphasis on the scalar field, it is desirable, owing to the definition of the SGS flux vector, that the SGS flux magnitude encompasses both the velocity and the scalar fields. Therefore, it is proposed to model the flux magnitude as the multiplication of an SGS velocity scale and an SGS scalar concentration scale, |**q**| = *usgsθsgs*. It is straightforward to assume that the SGS velocity scale is proportional to the square root

velocities on the basis of the local equilibrium hypothesis, which assumes a balance between

the resolved strain rate tensor) and dissipation rate *ε*. A classical evaluation of kinetic energy

resolved scales, a step that is consistent with the non-negative value of the dissipation rate;

 <sup>−</sup> *<sup>G</sup>ij Gkk*

*C*2 *ε*

where *H*(*x*) is the Heaviside step function defined as *H*(*x*) = 0 if *x* < 0 and *H*(*x*) = 1 if *x* ≥ 0. The local equilibrium hypothesis assumes a balance between the SGS scalar

the proposed model formulation together with the local equilibrium hypothesis, one obtains

dynamic procedure [28, 54] might serve to determine the model coefficients *C<sup>ε</sup>* and *Cεθ* , for simplicity, a simple method for determining a constant value has been adopted. To model the SGS scalar variance dissipation, Jiménez et al. [38] assumed that the SGS scalar mixing

that the Schmidt number (or the Prandtl number depending on the physical significance of

<sup>−</sup> *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup>

> *Sc Cε*

done) *Sc* <sup>=</sup> 0.71 (the Prandtl number of air near 20◦C), *<sup>C</sup><sup>ε</sup>* <sup>=</sup> 1, and *<sup>C</sup>* <sup>=</sup> <sup>√</sup>2 (with

*∂θ ∂xi* *∂ui*

where *Gij* <sup>=</sup> <sup>∆</sup><sup>2</sup> 12 *∂ui ∂xk ∂uj <sup>∂</sup>xk* for the isotropic grid, and the hyper-viscosity magnitude *ν<sup>u</sup>* can be calculated by either *<sup>ν</sup><sup>u</sup>* <sup>=</sup> *<sup>C</sup>*′ *<sup>s</sup>*∆4|*<sup>S</sup>*<sup>|</sup> or *<sup>ν</sup><sup>u</sup>* <sup>=</sup> *<sup>C</sup>*′ *<sup>k</sup>*∆<sup>3</sup>*ksgs*. Note that *<sup>C</sup>*′ *<sup>s</sup>* and *<sup>C</sup>*′ *<sup>k</sup>* can be determined by dynamic procedures, for the sake of simplicity, an empirical constant *<sup>C</sup>*′ *<sup>k</sup>* = 0.008 has been adopted. The development of this new nonlinear model also shows the importance of anisotropy and MFI-consistent requirements for SGS models in rotating system. At last, an approximation for *ksgs* is obtained solving the transport equation

$$\frac{\partial k\_{\rm sgs}}{\partial t} + \tilde{u}\_{\rm j} \frac{\partial k\_{\rm sgs}}{\partial \mathbf{x}\_{\rm j}} = -\pi\_{\rm ij} \frac{\partial \tilde{u}\_{\rm i}}{\partial \mathbf{x}\_{\rm j}} - \mathbb{C}\_{\varepsilon} \frac{k\_{\rm sgs}^{3/2}}{\Delta} + \frac{\partial}{\partial \mathbf{x}\_{\rm j}} \left[ \left( \frac{\nu\_{\rm k}}{\sigma\_{\rm k}} + \nu \right) \frac{\partial k\_{\rm sgs}}{\partial \mathbf{x}\_{\rm j}} \right]. \tag{8}$$

Here, the three terms on the right-hand side represent, respectively, the production, the dissipation, and the diffusion. The constants are typically chosen as *Ck* = 0.05, *C<sup>ε</sup>* = 1.0 and *σ<sup>k</sup>* = 1.0 based on previous studies [e.g., 44, 100].

Liu et al. [56] revisited this energy cascade issue and recommended another meaningful choice, a clipping procedure, to control the amplitude of backward cascade induced by the model. Inconveniently, Vreman et al. [95] have reported that their coupling of the gradient model with the clipping procedure still does not provide sufficient SGS dissipation in simulations of mixing layer flows.

Integrate the clipping procedure and the model of Pomraning and Rutland [74] and Lu et al. [61, 62], a simple SGS formulation for the LES of ABL flows has been proposed. The new algebraic nonlinear models for the SGS stress tensor and for the SGS flux vector can be written as

$$\mathbf{r}\_{l\bar{j}} = 2k\_{\text{sgs}} \left( \frac{\tilde{\mathbf{G}}\_{l\bar{j}}}{\tilde{\mathbf{G}}\_{kl}} \right) \quad \text{and} \quad q\_{\bar{i}} = |\mathbf{q}| \left( \frac{\tilde{\mathbf{G}}\_{\theta,\bar{i}}}{|\tilde{\mathbf{G}}\_{\theta}|} \right) \,. \tag{9}$$

The method separates the modeling into two elements: the normalized gradient term serves to model the structure (relative magnitude of each component); and a separate model is needed for the magnitudes. To account for the grid anisotropy in the study (<sup>∆</sup>*<sup>x</sup>*, <sup>∆</sup>*<sup>y</sup>* and <sup>∆</sup>*<sup>z</sup>* are not equal), it is defined that *<sup>G</sup>ij* <sup>=</sup> <sup>∆</sup><sup>2</sup> *x* 12 *∂ui ∂x ∂uj <sup>∂</sup><sup>x</sup>* + ∆2 *y* 12 *∂ui ∂y ∂uj <sup>∂</sup><sup>y</sup>* <sup>+</sup> <sup>∆</sup><sup>2</sup> *z* 12 *∂ui ∂z ∂uj <sup>∂</sup><sup>z</sup>* , and *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>=</sup> <sup>∆</sup><sup>2</sup> *x* 12 *∂ui <sup>∂</sup><sup>x</sup> ∂θ <sup>∂</sup><sup>x</sup>* + ∆2 *y* 12 *∂ui ∂y ∂θ <sup>∂</sup><sup>y</sup>* <sup>+</sup> <sup>∆</sup><sup>2</sup> *z* 12 *∂ui <sup>∂</sup><sup>z</sup> ∂θ <sup>∂</sup><sup>z</sup>* , and the gradient vector's magnitude is computed using the Euclidean norm <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> <sup>=</sup> *G*2 *<sup>θ</sup>*,1 <sup>+</sup> *<sup>G</sup>*<sup>2</sup> *<sup>θ</sup>*,2 <sup>+</sup> *<sup>G</sup>*<sup>2</sup> *<sup>θ</sup>*,3. To close the approach, it is required to evaluate the SGS kinetic energy, *ksgs* = <sup>1</sup> <sup>2</sup> *<sup>τ</sup>ii*, and the magnitude of the SGS flux vector, <sup>|</sup>**q**|. Even though a previous approach [20] places much emphasis on the scalar field, it is desirable, owing to the definition of the SGS flux vector, that the SGS flux magnitude encompasses both the velocity and the scalar fields. Therefore, it is proposed to model the flux magnitude as the multiplication of an SGS velocity scale and an SGS scalar concentration scale, |**q**| = *usgsθsgs*. It is straightforward to assume that the SGS velocity scale is proportional to the square root of the SGS kinetic energy, *usgs* = *C ksgs*. The value of *ksgs* is evaluated by using the resolved velocities on the basis of the local equilibrium hypothesis, which assumes a balance between the SGS kinetic energy production *<sup>P</sup>* (*<sup>P</sup>* = −*τij ∂ui <sup>∂</sup>xj* <sup>=</sup> <sup>−</sup>*τijSij*, where *<sup>S</sup>ij* <sup>=</sup> <sup>1</sup> 2 *∂ui ∂xj* <sup>+</sup> *<sup>∂</sup><sup>u</sup><sup>j</sup> ∂xi* is the resolved strain rate tensor) and dissipation rate *ε*. A classical evaluation of kinetic energy dissipation is *ε* = *C<sup>ε</sup> k*3/2 *sgs* <sup>∆</sup> . Simulations allow for no local energy transfer from unresolved to resolved scales, a step that is consistent with the non-negative value of the dissipation rate; this leads to

6 Computational and Numerical Simulations

where *Gij* <sup>=</sup> <sup>∆</sup><sup>2</sup>

12 *∂ui ∂xk ∂uj*

> *∂ksgs <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup><sup>j</sup>*

in simulations of mixing layer flows.

written as

∆2 *y* 12 *∂ui ∂y ∂θ <sup>∂</sup><sup>y</sup>* <sup>+</sup> <sup>∆</sup><sup>2</sup> *z* 12 *∂ui <sup>∂</sup><sup>z</sup> ∂θ* 

calculated by either *<sup>ν</sup><sup>u</sup>* <sup>=</sup> *<sup>C</sup>*′

The new mixed nonlinear model is defined as

that mixed gradient models can capture disequilibrium. In LES of rotating turbulence, to overcome the weakness that the traditional eddy viscosity (−2*νsgsSij*) is too dissipative at large scales, which may hinder kinetic energy from transferring to large scales, Lu et al. [62] considered a hyper-viscosity term (+*νu*∇2*<sup>S</sup>ij*) as suggested by previous researchers [9, 26].

> *Gij Gkk*

been adopted. The development of this new nonlinear model also shows the importance of anisotropy and MFI-consistent requirements for SGS models in rotating system. At last, an

*<sup>∂</sup>xk* for the isotropic grid, and the hyper-viscosity magnitude *ν<sup>u</sup>* can be

*<sup>k</sup>*∆<sup>3</sup>*ksgs*. Note that *<sup>C</sup>*′

*∂ ∂xj*  *νk σk* + *ν*

<sup>+</sup> *<sup>ν</sup>u*∇2*<sup>S</sup>ij* , (7)

*<sup>s</sup>* and *<sup>C</sup>*′

 *<sup>∂</sup>ksgs ∂xj*

*<sup>k</sup>* can be determined

*<sup>k</sup>* = 0.008 has

. (8)

*τij* = 2*ksgs*

*<sup>s</sup>*∆4|*<sup>S</sup>*<sup>|</sup> or *<sup>ν</sup><sup>u</sup>* <sup>=</sup> *<sup>C</sup>*′

approximation for *ksgs* is obtained solving the transport equation

= −*τij*

*∂ksgs ∂xj*

*σ<sup>k</sup>* = 1.0 based on previous studies [e.g., 44, 100].

*τij* = 2*ksgs*

are not equal), it is defined that *<sup>G</sup>ij* <sup>=</sup> <sup>∆</sup><sup>2</sup>

 *Gij Gkk* *x* 12 *∂ui ∂x ∂uj <sup>∂</sup><sup>x</sup>* + ∆2 *y* 12 *∂ui ∂y ∂uj <sup>∂</sup><sup>y</sup>* <sup>+</sup> <sup>∆</sup><sup>2</sup> *z* 12 *∂ui ∂z ∂uj*

by dynamic procedures, for the sake of simplicity, an empirical constant *<sup>C</sup>*′

*∂ui ∂xj*

− *C<sup>ε</sup>*

Here, the three terms on the right-hand side represent, respectively, the production, the dissipation, and the diffusion. The constants are typically chosen as *Ck* = 0.05, *C<sup>ε</sup>* = 1.0 and

Liu et al. [56] revisited this energy cascade issue and recommended another meaningful choice, a clipping procedure, to control the amplitude of backward cascade induced by the model. Inconveniently, Vreman et al. [95] have reported that their coupling of the gradient model with the clipping procedure still does not provide sufficient SGS dissipation

Integrate the clipping procedure and the model of Pomraning and Rutland [74] and Lu et al. [61, 62], a simple SGS formulation for the LES of ABL flows has been proposed. The new algebraic nonlinear models for the SGS stress tensor and for the SGS flux vector can be

The method separates the modeling into two elements: the normalized gradient term serves to model the structure (relative magnitude of each component); and a separate model is needed for the magnitudes. To account for the grid anisotropy in the study (<sup>∆</sup>*<sup>x</sup>*, <sup>∆</sup>*<sup>y</sup>* and <sup>∆</sup>*<sup>z</sup>*

, and *qi* = |**q**|

*<sup>∂</sup><sup>z</sup>* , and the gradient vector's magnitude is computed using the Euclidean

 *Gθ*,*i* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> . (9)

*x* 12 *∂ui <sup>∂</sup><sup>x</sup> ∂θ <sup>∂</sup><sup>x</sup>* +

*<sup>∂</sup><sup>z</sup>* , and *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>=</sup> <sup>∆</sup><sup>2</sup>

*k*3/2 *sgs* <sup>∆</sup> +

$$k\_{\rm sgs} = H\left(P\right) \frac{4\tilde{\Delta}^2}{\mathbb{C}\_\varepsilon^2} \left(-\frac{\tilde{\mathcal{G}}\_{ij}}{\tilde{\mathcal{G}}\_{kk}} \tilde{\mathcal{S}}\_{ij}\right)^2,\tag{10}$$

where *H*(*x*) is the Heaviside step function defined as *H*(*x*) = 0 if *x* < 0 and *H*(*x*) = 1 if *x* ≥ 0. The local equilibrium hypothesis assumes a balance between the SGS scalar variance production, *<sup>P</sup><sup>θ</sup>* = −*qi ∂θ ∂xi* , and the SGS scalar variance dissipation rate *εθ* . A classical evaluation of the SGS scalar variance dissipation rate is *εθ* = *Cεθ θ*2 *sgsusgs* <sup>∆</sup> . Using the proposed model formulation together with the local equilibrium hypothesis, one obtains *<sup>θ</sup>sgs* <sup>=</sup> <sup>∆</sup> *Cεθ* <sup>−</sup> *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> *∂θ ∂xi* . The clipping procedure assumes the SGS scalar variance production and the dissipation rate are always non-negative *<sup>θ</sup>sgs* <sup>=</sup> *<sup>H</sup>* (*P<sup>θ</sup>* ) <sup>∆</sup> *Cεθ* <sup>−</sup> *<sup>G</sup><sup>θ</sup>*,*<sup>i</sup>* <sup>|</sup>**<sup>G</sup>** *<sup>θ</sup>* <sup>|</sup> *∂θ ∂xi* . Although a dynamic procedure [28, 54] might serve to determine the model coefficients *C<sup>ε</sup>* and *Cεθ* , for simplicity, a simple method for determining a constant value has been adopted. To model the SGS scalar variance dissipation, Jiménez et al. [38] assumed that the SGS scalar mixing time is proportional to the SGS turbulent characteristic time, *<sup>θ</sup>*<sup>2</sup> *sgs εθ* <sup>∝</sup> *ksgs <sup>ε</sup>* . Tests have shown that the Schmidt number (or the Prandtl number depending on the physical significance of the scalar field) leads to satisfactory results. Using the model *εθ* = <sup>1</sup> *Sc θ*2 *sgsε ksgs* , one obtains the following equation for the magnitude of the SGS flux

$$|\mathbf{q}| = H\left(P\_{\theta}\right)H\left(P\right)\frac{1}{Sc}\frac{2\tilde{\Delta}^2}{\mathcal{L}\_{\varepsilon\theta}^2}\left(-\frac{\tilde{\mathcal{G}}\_{\theta,i}}{|\tilde{\mathbf{G}}\_{\theta}|}\frac{\partial\tilde{\theta}}{\partial x\_i}\right)\left(-\frac{\tilde{\mathcal{G}}\_{ij}}{\tilde{\mathcal{G}}\_{kk}}\tilde{\mathcal{S}}\_{lj}\right)\,. \tag{11}$$

Also, the coefficient can be derived as *Cεθ* = <sup>1</sup> *Sc Cε <sup>C</sup>* . When adopting (as this study has done) *Sc* <sup>=</sup> 0.71 (the Prandtl number of air near 20◦C), *<sup>C</sup><sup>ε</sup>* <sup>=</sup> 1, and *<sup>C</sup>* <sup>=</sup> <sup>√</sup>2 (with *usgs* <sup>=</sup> (*uiui* <sup>−</sup> *<sup>u</sup>iu<sup>i</sup>*) <sup>=</sup> <sup>2</sup> *ksgs*), one obtains *<sup>C</sup>εθ* <sup>=</sup> 1.0.
