Addressing the Grid-Size Sensitivity Issue in Large-Eddy Simulations of Stable Boundary Layers

Abstract

We have identified certain fundamental limitations of a mixing-length parametrization used in a popular turbulent kinetic energy-based subgrid-scale model. Replacing this parametrization with a more physically realistic one significantly improves the overall quality of the large-eddy simulation (LES) of stable boundary layers. For the range of grid sizes considered here (specifically, 1 m–12.5 m), the revision dramatically reduces the grid-size sensitivity of the simulations. Most importantly, the revised scheme allows us to reliably estimate the first- and second-order statistics of a well-known LES intercomparison case, even with a coarse grid size of O(10 m).

1 Introduction

The first large-eddy simulation (LES) intercomparison study (Beare et al. 2006), organized under the auspices of the Global Energy and Water Exchanges (GEWEX) Atmospheric Boundary Layer Study (GABLS), has had a lasting impact on research on the stable boundary layer (SBL). In the past decade and a half, the key findings from this study (henceforth referred to as GABLS1–LES) were cited by numerous papers; a few examples are:

“Adequate WSBL [weakly stable boundary layer] resolution is attained with 2-m grids (Beare et al. 2006), but higher resolution is required for moderate and very stable stratification (e.g., SBL depths of less than 50 m or so).” (Fernando and Weil 2010)

“Even for weak to moderately stable conditions, LES of the NBL [nocturnal boundary layer] requires a grid spacing of O(1 m) (Beare et al. 2006), which greatly increases the computational burden.” (van Stratum and Stevens 2015)

“In particular, it is often observed that grid convergence for simulations of the stable boundary layer is lacking, see Beare et al. (2006) and Sullivan et al. (2016). The latter used fine grid spacings down to [0.39] m (pseudo-spectral code) and still reported a sensitivity of their results to the grid spacing. Until now, a convincing explanation for this behaviour has been lacking, creating a limitation for the application of LES models for simulating the stable boundary layer.” (Maronga et al. 2019)

Similar statements on the grid-size sensitivity can be found in other peer-reviewed publications and are often heard in any contemporary workshop or conference session on SBL. Such an overwhelming consensus among the SBL–LES community at large is somewhat disconcerting given the fact that a handful of papers established a while ago that certain dynamic (tuning-free) subgrid-scale (SGS) models perform rather well (in terms of first- and second-order statistics) with coarse resolutions. As a matter of fact, around the time of publication by Beare et al. (2006), Basu and Porté-Agel (2006) demonstrated that the GABLS1–LES case can be simulated reliably by a dynamic SGS parametrization called the locally-averaged scale-dependent dynamic (LASDD) model. They concluded:

“Moreover, the simulated statistics obtained with the LASDD model show relatively little resolution dependence for the range of grid sizes considered here. In essence, it is shown here that the new LASDD model is a robust subgrid-scale parametrization for reliable, tuning-free simulations of stable boundary layers, even with relatively coarse resolutions.” (Basu and Porté-Agel 2006)

Later on, Stoll and Porté-Agel (2008) and Lu and Porté-Agel (2013) compared different dynamic SGS schemes and also reported negligible sensitivity to grid sizes.

In this study, we revisit the GABLS1–LES case study and probe into the inherent cause of the grid-size sensitivity of a static (non-dynamic) SGS parametrization. This parametrization was originally proposed by Deardorff (1980) and appears in many LES codes (e.g., Moeng et al. 2007; Heus et al. 2010; Maronga et al. 2015; Gibbs and Fedorovich 2016). Henceforth, we refer to this as the D80 parametrization. Several past SBL–LES studies have reported grid-size sensitivity with the D80 SGS scheme (e.g., Jiménez and Cuxart 2005; Beare et al. 2006; de Roode et al. 2017; Maronga et al. 2019). After extensive numerical experiments, we have identified the SGS mixing-length ($\lambda$) parametrization in this scheme to be at the root of the grid-size sensitivity issue. We have found that a rather simple (yet physically realistic) modification of $\lambda$ alleviates the grid-size sensitivity substantially. In addition, the first- and second-order statistics from the LES runs utilizing this revised parametrization (named D80-R) agrees well with the ones produced by a pseudo-spectral LES code utilizing a dynamic SGS model. Most importantly, the D80-R scheme allows us to reliably simulate the GABLS1–LES case even with a coarse grid-size of O(10 m).

The organization of this paper is as follows: in the following section we describe the D80 parametrization and its fundamental limitations. Physical interpretation and analytical derivations pertaining to the D80-R parametrization are included in Sect. 3, and technical details of the simulations are provided in Sect. 4. Results of simulations with D80 and D80-R SGS schemes are documented in Sect. 5. A few concluding remarks are made in Sect. 6. In Appendix 1, we have documented a few SGS models which do not use the grid size as a mixing-length scale; instead, similar to the D80-R parametrization, they use flow-physics-dependent mixing-length-scale formulations. To further confirm the validity of the D80-R approach, results from simulations with an independent LES code are included in Appendix 1.

2 SGS Parametrization

The SGS parametrization (D80) by Deardorff (1980) utilizes the following mixing-length scale

$$\begin{aligned} \lambda = \min \left( \varDelta _g, L_b \right) , \end{aligned}$$

(1)

where $\varDelta _g = \left( \varDelta x \varDelta y \varDelta z\right) ^{1/3}$ and $L_b$ is the so-called buoyancy length scale that is typically represented as follows

$$\begin{aligned} L_b = c_n \frac{{\overline{e}}^{1/2}}{N}, \end{aligned}$$

(2)

where ${\overline{e}}$ denotes SGS turbulence kinetic energy (TKE) and N is the Brunt–Väisäla frequency. Deardorff (1980) assumed $c_n$ to be equal to 0.76. Many LES studies still assume this constant value (e.g., Jiménez and Cuxart 2005; Heus et al. 2010; Maronga et al. 2015) though there are a handful of exceptions. For example, Gibbs and Fedorovich (2016) assumed $c_n = 0.5$.

The eddy viscosity ($K_m$), eddy diffusivity ($K_h$), and energy dissipation rate (${\overline{\varepsilon }}$) are all assumed to be functions of $\lambda$ as shown below

$$\begin{aligned} K_m= & {} c_m \lambda {\overline{e}}^{1/2}, \end{aligned}$$

(3a)

$$\begin{aligned} K_h= & {} c_h \lambda {\overline{e}}^{1/2}, \end{aligned}$$

(3b)

$$\begin{aligned} {\overline{\varepsilon }}= & {} \frac{c_\varepsilon }{\lambda } {\overline{e}}^{3/2}. \end{aligned}$$

(3c)

The unknown coefficients ($c_m$, $c_h$, and $c_\varepsilon$) are either prescribed or parametrized as follows

$$\begin{aligned} c_m= & {} 0.12, \end{aligned}$$

(4a)

$$\begin{aligned} c_h= & {} \left( 1 + 2\frac{\lambda }{\varDelta _g} \right) c_m, \end{aligned}$$

(4b)

$$\begin{aligned} c_\varepsilon= & {} 0.19 + 0.51 \frac{\lambda }{\varDelta _g}. \end{aligned}$$

(4c)

Please note that the values of the coefficients in Eq. 4 somewhat vary between different studies (e.g., Deardorff 1980; Moeng and Wyngaard 1988; Saiki et al. 2000).

For neutral conditions Eq. 1 reduces to $\lambda = \varDelta _g$. As a direct consequence, $K_m$, $K_h$, and ${\overline{\varepsilon }}$ become stability independent, as would be physically expected. Furthermore, SGS Prandtl number ($Pr_S = K_m/K_h = c_m/c_h$) is set equal to 0.33 for neutral conditions.

For stably stratified conditions, one expects that $K_m$, $K_h$, and ${\overline{\varepsilon }}$ should have a clear dependence on N. Deardorff (1980) accounted for such dependence by introducing $L_b$ in Eq. 1. Specifically, he stated (we have changed the variable notation for consistency):

“In past work it has been assumed that $\lambda = \varDelta _g$, which fails to take account of the possibility that in a stably stratified region $\lambda$ could become much smaller than the grid interval.” (Deardorff 1980)

For the very stable case, in the limit of $1/N \rightarrow 0$, the D80 SGS model predicts $\lambda = L_b \rightarrow 0$, $c_\varepsilon \rightarrow 0.19$, and $Pr_S \rightarrow 1$. Also, $K_m$ and $K_h$ are expected to approach negligibly small values for such conditions.

The behaviour of the D80 SGS model in the weakly and moderately stable boundary layer is rather problematic and is the focus of the present study. For such cases (including the GABLS1–LES case), $L_b$ is typically on the order of 10 m in the lower and middle parts of the SBL. As mentioned earlier, it has become customary to perform SBL–LES runs with grid sizes of 1–5 m or finer these days. In such simulations, by virtue of Eq. 1, $\lambda$ becomes equal to $\varDelta _g$ for lower and middle parts of the SBL. In other words, the ‘min’ operation in Eq. 1 is only utilized in the upper part of the SBL, and in the free atmosphere. Most importantly, akin to the neutral conditions, $K_m$, $K_h$, and ${\overline{\varepsilon }}$ become spuriously independent of stability in the lower and middle parts of the SBL.

There is another fundamental problem with Deardorff’s SGS mixing-length parametrization. It is not influenced by the presence of a surface. Deardorff (1980) recognized this problem and proposed a solution of increasing $c_\varepsilon$ near the surface (see also Moeng 1984). In this context, Gibbs and Fedorovich (2016) stated:

“To this end, it is unclear, however, whether the parameter adjustments incorporated in the original D80 scheme were based on some clear physical reasoning or were intended to merely produce more plausible effects close to the surface.”

To the best of our knowledge, such ad hoc solutions are not implemented in recent LES codes (e.g., Heus et al. 2010; Maronga et al. 2015; de Roode et al. 2017). In other static SGS models (e.g., the Smagorinsky–Lilly model and its variants), empirical wall functions are utilized to explicitly account for the near-surface shear effects (e.g., Mason and Derbyshire 1990; Brown et al. 1994). Interestingly, it is not a common practice to use wall functions with the D80 SGS scheme. In the case of dynamic SGS models, wall functions are not needed as the estimated SGS coefficient steadily decreases as one approaches any surface, and thus, reduces the SGS mixing length in a self-consistent manner (e.g., Basu and Porté-Agel (2006); Stoll and Porté-Agel (2008)).

In the present study, we replace the mixing-length parametrization in D80 (i.e., Eq. 1) with the following formulation

$$\begin{aligned} \frac{1}{\lambda } = \frac{1}{\kappa z} + \frac{1}{L_b}, \end{aligned}$$

(5)

where $\kappa$ is the von Kármán constant. Both the effects of stability and near-surface are nicely captured by this equation. In Sect. 3, we will show that the $\kappa z$ term of this equation can be derived from a well-known spectral scaling.

The origin of Eq. 5 can be traced back of Blackadar (1962) and Brost and Wyngaard (1978). Blackadar (1962) introduced the following length scale

$$\begin{aligned} \frac{1}{\lambda } = \frac{1}{\kappa z} + \frac{1}{\lambda _0}, \end{aligned}$$

(6)

where, $\lambda _0$ is an asymptotic length scale. Brost and Wyngaard (1978) and following studies (e.g., Baas et al. 2008) assumed

$$\begin{aligned} \lambda _0 = L_b^* \approx \frac{{\overline{E}}^{1/2}}{N}, \end{aligned}$$

(7)

where, ${\overline{E}}$ is the total TKE.

Please note that despite of apparent similarity, the length scales $L_b$ and $L_b^*$ are quite different. Scale $L_b$ is proportional to SGS TKE (${\overline{e}}$), and thus implicitly depends on the filter size $\varDelta _f$ (see Sect. 3). In contrast, $L_b^*$ is used in Reynolds-averaged Navier–Stokes (RANS) models, and by definition, it does not depend on $\varDelta _f$.

Hereafter, we refer to Deardorff’s SGS parametrization in conjunction with Eq. 5 as the D80-R parametrization. In this parametrization, with increasing resolution, both ${\overline{e}}$ and $L_b$ decrease. As a result, $\lambda$, $K_m$, and $K_h$ also decrease. Additional effect of grid-size is also felt via $c_h$, and $c_\varepsilon$ coefficients. In terms of SGS Prandtl number ($Pr_S$), both the D80 and D80-R parametrizations suffer from unphysical prescription in different ways. This issue is discussed in detail in Sect. 5.

3 Physical Interpretation and Analytical Derivation

In this section, we make direct associations between the energy spectra of turbulence and several key elements of the D80-R parametrization. In addition, we introduce a generalized form of Eq. 5 which has the potential to further extend the regime of application of the proposed D80-R parametrization.

3.1 Dependence of ${\overline{e}}$ on Filter Size ($\varDelta _f$)

A simple model of longitudinal velocity spectrum, spanning all the scales of turbulence, can be written as (see page 232 of Pope 2000)

$$\begin{aligned} E_u\left( k\right) \sim {\overline{\varepsilon }}^{2/3} k^{-5/3} \varPhi _\eta \left( k\eta \right) \varPhi _{L_I}\left( k L_I \right) , \end{aligned}$$

(8)

where k is wavenumber, $E_u$ is the energy spectrum for longitudinal velocity component, $L_I$ and $\eta$ denote integral length scale and Kolmogorov’s scale, respectively. The non-dimensional functions $\varPhi _{L_I}$ and $\varPhi _{\eta }$ represent the buoyancy range and dissipation range, respectively. In the inertial range, both these functions are close to unity. For small values of $k L_I$ and large values of $k \eta$, they deviate from unity.

If the filter size $\varDelta _f$ is within the inertial-range (as inherently assumed in LES), the subgrid-scale variance of longitudinal velocity component $\sigma _{us}^2$ can be estimated from Eq. 8 as

$$\begin{aligned} \sigma _{us}^2 \sim {\overline{\varepsilon }}^{2/3} \int _{k_\varDelta }^{k_\eta } k^{-5/3} \varPhi _\eta \left( k\eta \right) dk, \end{aligned}$$

(9)

where, $k_\eta$ is the dissipation wavenumber. The wavenumber associated with the filter is $k_\varDelta = \frac{\pi }{\varDelta _f}$. Within the range of $k_\varDelta$ to $k_\eta$, $\varPhi _{L_I}$ is unity. Since the contribution of dissipation range is typically small towards $\sigma _{us}^2$, one can assume $\varPhi _{\eta } \approx 1$ in Eq. 9.

Owing to isotropy in the inertial range, Eq. 9 can be integrated as follows

$$\begin{aligned} {\overline{e}} = \frac{3}{2} \sigma _{us}^2 \sim {\overline{\varepsilon }}^{2/3} \left[ \frac{1}{k_\varDelta ^{2/3}} - \frac{1}{k_\eta ^{2/3}} \right] . \end{aligned}$$

(10)

Since $k_\eta \gg k_\varDelta$, we can simplify this equation as

$$\begin{aligned} {\overline{e}} \sim {\overline{\varepsilon }}^{2/3} \varDelta _f^{2/3}. \end{aligned}$$

(11)

Eq. 11 can be re-written as

$$\begin{aligned} {\overline{\varepsilon }} \sim \frac{{\overline{e}}^{3/2}}{\varDelta _f}. \end{aligned}$$

(12)

This equation is analogous to Eq. 3c, if one replaces $\lambda$ with $\varDelta _f$.

3.2 Estimation of ${\overline{e}}$ in Surface Layer

The derivations in Sect. 3.1 assume $\varDelta _f$ falls within the inertial range. However, in the surface layer, the inertial range is rather limited. For an extensive range of scales, the longitudinal velocity spectrum follows a $k^{-1}$ power law. Thus, in most LES studies, it is likely that $\varDelta _f$ falls within the $k^{-1}$ range in the surface layer and not within the inertial range.

Following Tchen (1953), numerous studies have reported the $k^{-1}$ scaling in the literature; please refer to a comprehensive list in Table 1 of Katul and Chu (1998). By combining the $k^{-1}$ scaling with the inertial-range scaling (i.e., $k^{-5/3}$), the energy spectrum for the longitudinal velocity spectrum can be written as

$$\begin{aligned}&E_u \sim u_*^2 k^{-1} \text{ for } k_\varDelta \le k \le k_o{,} \end{aligned}$$

(13a)

$$\begin{aligned}&E_u \sim {\overline{\varepsilon }}^{2/3} k^{-5/3} \varPhi _\eta \left( k\eta \right) \text{ for } k_o \le k \le k_\eta {,} \end{aligned}$$

(13b)

where the friction velocity is denoted by $u_*$. The crossover wavenumber $k_o$ equals to $\frac{1}{\gamma z}$. For unstable conditions, $\gamma$ was found to be equal to unity by Kader and Yaglom (1991) and others. For stable conditions, the value of $\gamma$ decreases from unity with increasing stability (Banerjee et al. 2016).

By integrating Eq. 13, we can estimate the SGS variance of longitudinal velocity component as follows

$$\begin{aligned} \sigma _{us}^2 = c_1 u_*^2 \int _{k_\varDelta }^{k_o} k^{-1} dk + c_2 {\overline{\varepsilon }}^{2/3} \int _{k_o}^{k_\eta } k^{-5/3} \varPhi _\eta \left( k\eta \right) dk, \end{aligned}$$

(14)

where, $c_1$ and $c_2$ are unknown constants. Due to strong anisotropy in the surface layer, we cannot estimate SGS TKE (${\overline{e}}$) from $\sigma _{us}^2$. However, it is expected that ${\overline{e}}$ will be proportional to $\sigma _{us}^2$.

The first integration term of Eq. 14 can be simplified as

$$\begin{aligned} \sigma _{us1}^2 \sim u_*^2 \int _{k_\varDelta }^{k_o} k^{-1} dk = u_*^2 \log \left( \frac{k_o}{k_\varDelta }\right) = u_*^2 \log \left( \frac{\varDelta _f}{\pi \gamma z}\right) . \end{aligned}$$

(15)

Whereas the second term of Eq. 14 amounts to

$$\begin{aligned} \sigma _{us2}^2 \sim {\overline{\varepsilon }}^{2/3} \int _{k_o}^{k_\eta } k^{-5/3} dk \sim {\overline{\varepsilon }}^{2/3} \left[ \frac{1}{k_o^{2/3}} - \frac{1}{k_\eta ^{2/3}} \right] \sim {\overline{\varepsilon }}^{2/3} \left( \gamma z \right) ^{2/3}. \end{aligned}$$

(16)

Please note that, according to Eq. 16, $\sigma _{us2}^2$ is not dependent on $\varDelta _f$, rather it is dependent on z. Due to logarithmic operation in Eq. 15, $\sigma _{us1}^2$ weakly depends on $\varDelta _f$. Thus, ${\overline{e}}$ is expected to be weakly dependent on $\varDelta _f$ in the surface layer. More importantly, Eq. 12 cannot be used to estimate ${\overline{\varepsilon }}$ in the surface layer.

In general, both the terms $\sigma _{us1}^2$ and $\sigma _{us2}^2$ contribute to ${\overline{e}}$. It is not possible to neglect either of them for further simplification. However, if we assume that both these variances are proportional to ${\overline{e}}$, we get qualitatively similar results regarding the mixing length.

Based on Eq. 16, we can write

$$\begin{aligned} {\overline{e}} \propto \sigma _{us2}^2 \sim {\overline{\varepsilon }}^{2/3} \left( \gamma z \right) ^{2/3}, \end{aligned}$$

(17a)

and

$$\begin{aligned} {\overline{\varepsilon }} = \frac{c_\gamma }{z} {\overline{e}}^{3/2}, \end{aligned}$$

(17b)

where, $c_\gamma$ is an unknown constant. By comparing Eq. 3c with Eq. 17b, we can assert that the mixing length scale $\lambda$ in D80 and D80-R should be proportional to z in the surface layer instead of $\varDelta _f$. In the following section, we arrive at the same conclusion via a different route.

3.3 Shear-Based Mixing Length Scale

Let us assume that S represents the magnitude of velocity shear. Thus, gradient Richardson number ($Ri_g$) equals to $N^2/S^2$. For shear-dominated flows, Hunt et al. (1988, (1989) proposed a length scale

$$\begin{aligned} L_H^* \sim \frac{{\overline{E}}^{1/2}}{S}, \end{aligned}$$

(18)

where ${\overline{E}}$ is the total TKE. It is a common knowledge that the effects of shear are more prevalent than the buoyancy effects near the surface. For such a situation, $L_H^*$ is a more relevant length scale than $L_b^*$ (defined earlier in Eq. 7).

For SGS modelling, if ${\overline{E}}$ is replaced with ${\overline{e}}$, an analogous shear-based length scale can be defined

$$\begin{aligned} L_H \sim \frac{{\overline{e}}^{1/2}}{S}. \end{aligned}$$

(19)

In contrast to $L_H^*$, the length scale $L_H$ implicitly depends on $\varDelta _f$ due to its explicit relationship to ${\overline{e}}$.

Earlier we showed that Eq. 15 holds near the surface. Now, if we assume ${\overline{e}}$ is proportional to $\sigma _{us1}^2$, we can combine Eq. 15 with Eq. 19 and obtain

$$\begin{aligned} \lambda \sim L_H \sim \frac{u_*}{S} \left[ \log \left( \frac{\varDelta _f}{\pi \gamma z} \right) \right] ^{1/2} = \frac{\kappa z}{\phi _m} \left[ \log \left( \frac{\varDelta _f}{\pi \gamma z} \right) \right] ^{1/2}, \end{aligned}$$

(20)

where $\phi _m$ is the non-dimensional velocity gradient.

For finite-difference-based LES codes, $\varDelta _f$ is typically 4–6 times larger than $\varDelta _g$. For isotropic grids, the height of vertical levels (z) equal to $m \varDelta _g$; where, m are half-integers. For vertically stretched grids, z could be a small fraction of $\varDelta _g$ for the first few levels. For stably stratified conditions, $\phi _m \ge 1$. As mentioned earlier, in such situation $\gamma \le 1$. Based on these estimates, it is reasonable to state that $\lambda$ is proportional to z in the surface layer; the exact value of the proportionality constant is unknown. In this study, based on other usage described in the literature (specifically in RANS modelling), we assumed the proportionality constant to be equal to $\kappa$.

3.4 A Generalized Mixing-Length Parametrization

The proposed D80-R mixing-length parametrization (i.e., Eq. 5) is valid for stable conditions. However, for near-neutral conditions, as N approaches zero, the value of $\lambda$ could become very large. To account for such stability regimes, one can adopt a more generic parametrization for the mixing length. Essentially, one can combine $L_H$ and $L_b$ in the following manner

$$\begin{aligned} \frac{1}{\lambda } = \frac{c_H}{L_H} + \frac{c_b}{L_b}, \end{aligned}$$

(21)

where the unknown coefficients $c_H$ and $c_b$ should be prescribed. For near-neutral condition, the term involving $L_H$ will dominate and will lead to realistic $\lambda$ values. As discussed in the previous section, this term will also perform well in the surface layer.

Unfortunately, we do not know how to estimate optimal values of $c_H$ and $c_b$ in a meaningful way. Running numerous large-eddy simulations for various case studies with different combinations of $c_H$ and $c_b$ is not a computationally viable option. Hopefully, an efficient strategy will emerge in the near future. For the time-being, we utilize an approximation of Eq. 21, i.e., Eq. 5, as a working substitute for SBL simulations.

4 Description of the Simulations

In this work, we simulate the GABLS1–LES case study using the Dutch Atmospheric Large-Eddy Simulation (DALES; Heus et al. 2010) and the PALM model system (Maronga et al. 2015, 2020). Since the configurations of the GABLS1–LES case study are well described in the literature, we mention them in a succinct manner. The boundary layer is driven by an imposed uniform geostrophic wind speed of 8 m s$^{-1}$ and a surface cooling of 0.25 K h$^{-1}$. It attains a quasi-steady state in about 8–9 h with a depth of approximately 200 m. The initial mean potential temperature is 265 K up to 100 m with an overlying inversion of 0.01 K m$^{-1}$. The Coriolis parameter is set to $1.39 \times 10^{-4}$ s$^{-1}$, corresponding to latitude 73$^{\circ }$ N. Both the aerodynamic roughness length $z_0$ and the scalar roughness length $z_{0h}$ are assumed to be equal to 0.1 m.

For all runs, the computational domain is fixed at 400 m $\times$ 400 m $\times$ 400 m. A wide range of isotropic $\varDelta _g$ values are used to investigate the aforementioned grid-size sensitivity issue. For the DALES code, in order to avoid any temporal discretization error, the time step $\varDelta t$ is kept at a constant value of 0.1 s for all the simulations. In contrast, an adaptive time-stepping approach is used by the PALM model system. In addition to the results from the finite-difference-based DALES and PALM codes, we also report results from a pseudo-spectral code (called MATLES) utilizing the LASDD SGS model along with a grid size of 3.5 m and a fixed timestep of 0.075 s.

Table 1 SGS coefficients and mixing length values in the DALES and PALM codes

In terms of the numerical schemes in the DALES and PALM codes, a third-order Runge–Kutta scheme is used for time integration in conjunction with a fifth-order advection scheme in the horizontal direction (Wicker and Skamarock 2002). In the vertical direction, a second-order and a fifth-order scheme (which reduces to a second-order scheme near the surface) are used by the DALES and PALM codes, respectively. The MATLES code uses a second-order Adams–Bashforth scheme for time advancement.

In all the simulations, the lower boundary condition is based on the Monin–Obukhov similarity theory (MOST). As discussed by Basu and Lacser (2017), in order to apply MOST, the first model level in an LES model should not be at a height lower than $\sim 50 z_0$. In this study, for simplicity, this requirement has been violated for all the runs involving high-resolution. At this point, the impact of this violation on the simulated statistics is not noticeable and a thorough investigation on surface-layer physics will be conducted in a follow-up study.

Before discussing the results, we point out that some of the prescribed SGS constants differ between the DALES and the PALM codes. For example, $c_m$ is assumed to be equal to 0.12 in DALES; whereas, PALM sets it at 0.1. The other specifications related to $c_h$ and $\lambda$ depend on local $N^2$ values and are listed in Table 1. In addition, the coefficients in Eq. 4c are also slightly different in DALES and PALM. In spite of these differences, the results from DALES and PALM codes follow the same trends as depicted in the following section and in Appendix 1.

5 Results

The left and right panels of all the following figures represent the statistics from the D80- and D80-R-based runs, respectively. As prescribed in the GABLS1–LES study, all the statistics are averaged over the last one hour (i.e., 8–9 h) of the simulations.

Fig. 1

Vertical profiles of mixing length ($\lambda$) (top panel) and turbulent Prandtl number (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 2

Time series of surface friction velocity (top panel) and sensible heat flux (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Prior to discussing the first- and second-order statistics, it is important to highlight the differences between the original and the revised runs in terms of the SGS mixing-length ($\lambda$) profiles. From Fig. 1 (top panel), it is evident that in the D80-based runs (left panel), $\lambda$ equals to $\varDelta _g$ for the lower and middle parts of the SBL. In contrast, the D80-R-based runs show clear dependence on the distance from the surface. In both types of runs, $\lambda$ values monotonically decrease to zero in the upper part of the SBL and in the free atmosphere due to increasing stratification (as quantified by N). The $\lambda$ profiles from the D80-R-based simulations also show clear dependence on grid size in the core region of the SBL. This is due to the fact that SGS TKE decreases with increasing resolution, and in turn, decreases $\lambda$.

In the D80-based runs, due to the underlying condition of $\lambda = \varDelta _g$, $Pr_S$ becomes exactly equal to 0.33 for the lower part of the SBL (bottom-left panel of Fig. 1). Whereas, in the case of D80-R, $\lambda$ is typically smaller than $\varDelta _g$ near the surface. Thus, $Pr_S$ is much larger than 0.33. However, in the middle part of the SBL, depending on stratification, $\lambda$ may become larger than $\varDelta _g$. For such cases, in D80-R-based runs, $Pr_S$ can be even smaller than 0.33 (bottom-right panel of Fig. 1). In the upper part of the SBL, due to stronger stratification, $\lambda$ becomes much smaller than $\varDelta _g$, and as a consequence, $Pr_S$ monotonically approaches to unity for both the D80 and D80-R cases. In the case of the LASDD SGS model, the dynamically estimated $Pr_S$ values remain more or less constant (around 0.5–0.6) inside the SBL and becomes greater than 1 in the inversion layer.

Fig. 3

Vertical profiles of mean wind speed (top panel), wind direction (middle panel), and potential temperature (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Gibbs and Fedorovich (2016) recognized the problem with the $Pr_S$ in D80 parametrization. They proposed to use $Pr_S = 1$ when $N^2$ is locally positive. For locally negative $N^2$ values they proposed an empirical formulation for $Pr_S$. In the present study, thus, the PALM model is employed with $Pr_S = 1$ for $N^2 > 0$ in all the D80-R-based runs. These results are shown in Appendix 2. Please note that for $N^2 \le 0$, both the DALES and the PALM models always use Eq. 4b (i.e., they effectively assume $Pr_S$ = 0.33); we have not incorporated the empirical formulation by Gibbs and Fedorovich (2016).

Fig. 4

Vertical profiles of total (top panel) and resolved (bottom panel) momentum-flux (x-component) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Simulated time series of surface friction velocity and sensible heat flux are documented in Fig. 2. In the D80-based runs, no temporal fluctuations of surface fluxes are visible for $\varDelta _g \ge$ 6.25 m. When $\varDelta _g$ is finer than 6.25 m, temporal fluctuations do appear approximately 1–2 hours into the simulation. Increasing grid-resolution helps in systematically reducing this turbulence kick-off time. In the D80-R-based runs, similar artifacts are not visible in the dynamical evolution of surface fluxes. In the D80-based runs for $\varDelta _g \ge$ 6.25 m, $\lambda$ (= $\varDelta _g$) values are excessively large near the surface and simply do not allow for the sustenance of turbulence. The inclusion of the surface dependence term in Eq. 5 reduces the mixing length in a meaningful way and promotes the production and transport of turbulence.

In Fig. 2, the simulated time series from the MATLES code are also overlaid. In the D80-R-based runs, the agreement between the DALES- and MATLES-based results is remarkable. All the time series of surface friction velocity reach quasi-equilibrium stage after approximately 5 h of simulation. The surface sensible flux time series reach quasi-equilibrium stage at a later time ($\sim$ 6 h).

The vertical profiles of mean wind speeds are included in the top panel of Fig. 3. The presence of the low-level jet (LLJ) is clearly visible in both the plots. However, the height of the LLJ peak shows strong sensitivity with respect to grid-size in the D80-based runs when $\varDelta _g \ge$ 6.25 m. This unphysical behaviour is solely due to the inherent limitations of the D80 SGS parametrization. Once $\varDelta _g$ becomes smaller than 6.25 m, turbulence is sustained near the surface, and as a result of adequate diffusion, the simulated LLJ peak heights are elevated. With further increase in spatial resolution, the mean wind-speed profiles reach convergence; albeit, the height of the DALES-based simulated LLJ peak remains lower and weaker than the one simulated by the MATLES code. The positive impacts of the revised mixing-length parametrization in D80-R can be seen in the top-right panel of Fig. 3. Almost all the DALES-based runs (with the exception of the one with $\varDelta _g =$ 12.5 m) converge on a single line and also agree strongly with the MATLES-based results.

Utilizing the results from the GABLS1 intercomparison study, Svensson and Holtslag (2009) investigated the wind turning in SBL in great detail. With certain assumptions, they analytically derived the relationship between the vertically integrated cross-isobaric flow and surface friction. Earlier, we have shown that the runs using the D80-R parametrization lead to $u_*$ series which are insensitive to grid sizes. Thus, it is not surprising that those simulations also lead to wind direction profiles which are in good agreement with one another. From the middle panel of Fig. 3, it is also clear that the original D80 parametrization performs quite poorly in capturing the turning of wind with height for coarse grids.

Fig. 5

Vertical profiles of total (top panel) and resolved (bottom panel) momentum-flux (y-component) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

The profiles of mean potential temperature are shown in the bottom panel of Fig. 3. All the profiles do show similar convex curvature as reported by earlier studies on GABLS1–LES. However, as with the mean wind-speed profiles, the D80-based runs exhibit strong dependence on grid size. Once again, incorporation of the revised SGS mixing length parametrization leads to better convergence. However, in this case, lower values of $Pr_S$ in the middle part of the SBL (see bottom-right panel of Fig. 1) cause more heat diffusion; not surprisingly, the potential temperature profiles from the D80-R-based runs are more convex than the one simulated by the MATLES code. In contrast, by using a higher value of $Pr_S$, the PALM model generates potential temperature profiles which are indistinguishable from the MATLES-based ones (refer to Fig. 11 in Appendix 2).

Fig. 6

Vertical profiles of total (top panel) and resolved (bottom panel) sensible heat flux from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 7

Vertical profiles of resolved (top panel) and subgrid-scale (bottom panel) TKE from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison in the resolved TKE plots

In the D80-based runs with finer resolutions, the simulated mean potential temperature profiles also start to converge. However, for these cases, $Pr_S$ values are equal to 0.33 in the lower part of the SBL. The (negative) impact of such high eddy diffusivity values is hard to notice in the mean potential temperature profiles; however, the vertical profiles of variance of potential temperature, $\sigma _\theta ^2$, do show the effect clearly (discussed later).

The x- and y-components of momentum-flux profiles are shown in Figs. 4 and  5, respectively. Several observations can be made from these plots. First of all, the resolved momentum fluxes are virtually non-existent for the D80-based runs with $\varDelta _g \ge 6.25$ m. This result is not surprising given the other statistics shown in the previous plots. The total fluxes are strongly grid-size dependent for the D80-based runs; however, they are not for the D80-R-based ones. Also, the fluxes are well resolved in the revised case; near the surface, the resolved fluxes increase with increasing resolution as would be expected. The revised results are more or less in line with the fluxes generated by the MATLES code.

Fig. 8

Vertical profiles of resolved variances of vertical velocity (top panel) and potential temperature (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the DALES code. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

As depicted in Fig. 6, the overall behaviour of the total and resolved sensible heat flux profiles are qualitatively similar to those of the momentum-flux profiles. In the case of the D80-R-based runs, the grid convergence is slightly less than satisfactory for the total sensible heat-flux profiles. Especially, the simulation with $\varDelta _g = 12.5$ m consistently overestimates the magnitude of heat flux across the SBL.

The resolved and SGS TKE plots are shown in Fig. 7. The MATLES code does not solve for the SGS TKE equation; thus, only the resolved TKE values are overlaid for comparison. As with the momentum-flux and sensible heat flux profiles, the resolved TKE is non-existent in the bottom-half of the SBL for the D80-based simulations using $\varDelta _g \ge$ 6.25 m. In the lower part of the SBL, the resolved TKE values are larger in the D80-R-based runs in comparison to the D80-based ones. In that region, with increasing resolution, resolved TKE values increase as would be physically expected. Most importantly, the SGS TKE values are much smaller in magnitude than their resolved counter parts (especially, in the lower part of the SBL). In other words, the flow is highly resolved for all the simulations involving the D80-R parametrization.

Lastly, the resolved variances of vertical velocity and potential temperature are plotted in Fig. 8. The trends of the resolved $\sigma _w^2$ profiles are similar to the resolved TKE plots. The resolved $\sigma _\theta ^2$ profiles from the D80-based runs are quite interesting. As discussed earlier, the coarse-resolution runs (i.e., $\varDelta _g \ge$ 6.25 m) show the quasi-laminarization problem up to 75 m or so. Interestingly, for higher resolution runs with the original mixing-length parametrization, the resolved $\sigma _\theta ^2$ values decrease significantly near the surface. The opposite trend is seen in the D80-R-based simulations. We believe this decrease in the D80-based runs is due to the low value of $Pr_S$ (= 0.33) in the bottom part of the SBL. Similar decrease in $Pr_S$ in the core of the SBL also causes resolved $\sigma _\theta ^2$ to significantly decrease in the D80-R-based cases.

Even though most vertical profiles from the DALES code look physically realistic and agree quite well with the corresponding results from the MATLES code, some discrepancies are noticeable in the case of resolved variances and fluxes. It is possible that the LASDD SGS model is slightly over-dissipative near the surface; the alternative scenario is that the proposed D80-R SGS model is slightly under-dissipative near the surface. Typically, in pseudo-spectral codes, spectral analysis can shed light into such undesirable behaviours (see Anderson et al. 2007 for some examples). Spurious piling up of energy near the high-wavenumber part of the spectrum is a telltale sign of under-dissipation. Since DALES and PALM are finite-difference codes, detection of such a subtle feature in the energy spectrum is a challenging task due to strong numerical dissipation. As such, we will be implementing the D80-R SGS model in the MATLES code for in-depth spectral analysis.

6 Discussion

In addition to Deardorff’s SGS model, the Smagorinsky–Lilly SGS model (Smagorinsky 1963; Lilly 1966a, b) is also quite popular in the boundary-layer community. In this SGS model, the effective mixing length is $C_S \varDelta _f$; where, $C_S$ is the so-called Smagorinsky coefficient. This coefficient is adjusted empirically (e.g., Mason and Derbyshire 1990; Brown et al. 1994; Kleissl et al. 2004) or dynamically (e.g., Porté-Agel et al. 2000; Bou-Zeid et al. 2005; Basu and Porté-Agel 2006) to account for shear, stratification, and grid resolution.

In contrast, the effective mixing length is $c_m \lambda$ in D80 SGS model. Most of the LES codes assume a constant value for $c_m$ and expect $\lambda$ to account for shear, stratification, and grid resolution. As elaborated in Sect. 2, the original mixing-length formulation does not account for shear or stratification properly. The parametrization by Brost and Wyngaard (1978) provides a major improvement as both the shear and stratification effects are now explicitly included. The effects of grid size is somewhat indirect. The performance of the D80-R SGS model may be improved if $c_m$ is assumed to be $\varDelta _g$-dependent and estimated via a dynamic approach (e.g., Ghosal et al. 1995; Krajnović and Davidson 2002). Our future work will be in that direction.

We will also need a better understanding of the SGS Prandtl number ($Pr_S$) and its relationship to the turbulent Prandtl number ($Pr_T$). Several field and laboratory experiments have reported that $Pr_T$ should be on the order of one for stably stratified conditions (see Sukoriansky et al. 2006; Li 2019, and the references therein). However, within the SBL (excluding surface and inversion layers), the dynamic SGS models typically estimate $Pr_S$ and $Pr_T$ values to be around 0.5–0.6 and 0.7, respectively (Basu and Porté-Agel 2006; Stoll and Porté-Agel 2008). We agree with Li (2019) that $Pr_S$ is less important than $Pr_T$ as the large-scale fluxes are resolved in LES. Nonetheless, we do believe the current $Pr_S$ parametrizations in the D80 and D80-R SGS models are not acceptable and definitely need amendments. The formulation proposed by Gibbs and Fedorovich (2016) is a good starting point and should be rigorously tested in future studies.

7 Concluding Remarks

In this study, we have demonstrated that Deardorff’s mixing-length parametrization is not suitable for SBL simulations. Instead an older scheme, proposed by Brost and Wyngaard for RANS, gives promising results. Even though we have made some progress in the arena of LES of SBL, many open questions still remain:

• Would this new mixing-length parametrization work for very stable boundary layers? Would it allow us to simulate turbulent bursting events?

• Is the proposed parametrization under-dissipative near the surface? How can we (dis)prove this behaviour?

• Is a buoyancy-based length scale really appropriate for weakly or moderately stable boundary layer? Or, should we opt for a shear-based length scale (Hunt et al. 1988, 1989; Basu et al. 2020)?

• How sensitive are the simulated results with respect to the choice of SGS coefficients (i.e., $c_n$, $c_m$, $c_\varepsilon$)? Should they be dynamically determined instead of being prescribed?

• How should we parametrize the SGS Prandtl number? Should it be a function of point-wise gradient Richardson number?

A few years ago, Basu and Lacser (2017) cautioned against violating MOST in LES runs with very high resolutions. To overcome this issue, Maronga et al. (2019) proposed certain innovative strategies; however, the results were somewhat inconclusive—(possibly) due to the usage of the D80 mixing-length parametrization in all their simulations. In light of the findings from the present work, we will revisit the MOST issue in very high resolution LES in conjunction with the D80-R SGS model. In addition, we will investigate the interaction of the D80-R parametrization with a coupled land-surface model.

We further recommend the SBL–LES community to revisit LES intercomparison studies organized under GABLS with revised or newly proposed SGS models. We speculate that some of the conclusions from the previous studies will no longer be valid.

8 Data and Code Availability

The DALES code is available from: https://github.com/dalesteam/dales. The PALM model system is available from: https://palm.muk.uni-hannover.de/trac. The MATLES code is available from S. Basu upon request. Upon acceptance of the manuscript, all the analysis codes and processed data will be made publicly available at http://doi.org/10.5281/zenodo.3972345.

Appendices

Appendix 1: SGS Models with Implicit $\varDelta _f$-Dependence

Traditional eddy-viscosity models, commonly used by the boundary-layer meteorology community, explicitly involve the filter-size ($\varDelta _f$) as the mixing-length scale. For example, in the case of Smagorinsky–Lilly SGS model, one uses

$$\begin{aligned} K_m= & {} \left( C_S \varDelta _f\right) ^2 |{\tilde{S}}|, \end{aligned}$$

(22a)

$$\begin{aligned} \lambda= & {} \left( C_S \varDelta _f\right) , \end{aligned}$$

(22b)

where $|{\tilde{S}}|$ is the magnitude of the resolved strain rate tensor and $C_S$ is known as the Smagorinsky coefficient.

In contrast to these popular SGS models, the proposed D80-R SGS model depends on $\varDelta _f$ in an implicit manner. Similar implicit $\varDelta _f$-dependence can be found in the engineering LES literature. Recently, Piomelli and Geurts (2011) proposed an eddy-viscosity model

$$\begin{aligned}&K_m \propto \frac{{\overline{e}}_r}{\omega _r^2}|{\tilde{S}}|, \end{aligned}$$

(23a)

$$\begin{aligned}&\lambda \propto \frac{{\overline{e}}_r^{1/2}}{\omega _r}, \end{aligned}$$

(23b)

where, ${\overline{e}}_r$ and $\omega _r$ denote resolved TKE and resolved magnitude of vorticity, respectively. In a more recent work, Piomelli et al. (2015) formulated the integral length scale approximation (ILSA) model

$$\begin{aligned}&K_m \propto \frac{{\overline{e}}_r^3}{{\overline{\varepsilon }}_T^2}|{\tilde{S}}|, \end{aligned}$$

(24a)

$$\begin{aligned}&\lambda \propto \frac{{\overline{e}}_r^{3/2}}{{\overline{\varepsilon }}_T}, \end{aligned}$$

(24b)

where ${\overline{\varepsilon }}_T$ represents the total energy dissipation rate. A dynamic version of the ILSA model was developed by Rouhi et al. (2016) and utilized for simulations in complex geometries by Lehmkuhl et al. (2019). The advantages of these flow-physics-dependent length scales (i.e., Eq. 23b and Eq. 24b) over $\varDelta _f$ as a mixing-length scale have been extensively discussed by Piomelli (2014) and Geurts et al. (2019) and will not be repeated here for brevity.

In addition to the aforementioned eddy-viscosity type SGS models, certain non-eddy-viscosity type SGS models also do not include explicit dependency on $\varDelta _f$. A case in point is the similarity model and its variants (Bardina et al. 1980; Liu et al. 1994).

Appendix 2: Results from the PALM Model System

The simulated results from the PALM model are documented in this appendix (see Figs. 9, 10, 11, 12, 13, 14, 15 and 16). The trends are very similar to those reported in Sect. 4. Thus, we do not discuss most of these figures for brevity. However, we would like to point out a noticeable feature in the D80-based runs. Even though most of the simulated profiles converge for 2 m $\le \varDelta _g \le$ 5 m, the results from $\varDelta _g = 1$ m run stands out. We believe that $\lambda$ is quite low for this particular case; even though it is small enough to sustain turbulence near the surface, it is not large enough to promote diffusion. The D80-R run with $\varDelta _g = 1$ m does not portray such an unusual behaviour.

The (positive) impacts of higher SGS Prandtl number ($Pr_S$) utilized in the D80-R-based runs can be seen in some of the figures. First of all, Fig. 11 (bottom-right panel) shows that the vertical profiles of potential temperature are less convex than the corresponding profiles from the DALES code (i.e., bottom-right panel of Fig. 3). This is due to less heat diffusion. For the very same reason, the PALM model-generated resolved $\sigma _\theta ^2$ values are much larger than their DALES counterpart (see Fig. 16).

Fig. 9

Vertical profiles of mixing-length $\lambda$ from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 10

Time series of surface friction velocity (top panel) and sensible heat flux (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 11

Vertical profiles of mean wind speed (top panel), wind direction (middle panel), and potential temperature (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 12

Vertical profiles of total (top panel) and resolved (bottom panel) momentum-flux (x-component) from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 13

Vertical profiles of total (top panel) and resolved (bottom panel) momentum-flux (y-component) from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 14

Vertical profiles of total (top panel) and resolved (bottom panel) sensible heat flux from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison

Fig. 15

Vertical profiles of resolved (top panel) and subgrid-scale (bottom panel) turbulence kinetic energy from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison in the resolved TKE plots

Fig. 16

Vertical profiles of resolved variances of vertical velocity (top panel) and potential temperature (bottom panel) from the D80 (left panel) and D80-R (right panel) based simulations using the PALM model system. Different coloured lines correspond to different grid sizes ($\varDelta _g$). Results from the MATLES code are overlaid (dashed black lines) for comparison