Addressing the Grid-size Sensitivity Issue in Large-eddy Simulations of Stable Boundary Layers

In this study, we have identified certain fundamental limitations of a mixing length parameterization used in a popular turbulent kinetic energy-based subgrid-scale model. Replacing this parameterization with a more physically realistic one significantly improves the overall quality of the large-eddy simulations (LESs) 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).


I. INTRODUCTION
The first large-eddy simulation (LES) intercomparison study [1], organized under the auspices of the Global Energy and Water Exchanges (GEWEX) Atmospheric Boundary Layer Study (GABLS), has had a lasting impact on the stable boundary layer (SBL) research field. 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 [1], but higher resolution is required for moderate and very stable stratification (e.g., SBL depths of less than 50 m or so)." [2] "Even for weak to moderately stable conditions, LES of the NBL [nocturnal boundary layer] requires a grid spacing of O(1 m) [1], which greatly increases the computational burden." [3] "In particular, it is often observed that grid convergence for simulations of the stable boundary layer is lacking, see [1] and [4]. 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." [5] Similar statements on the grid-size sensitivity can be found in other peer-reviewed publications and are often heard a few times in any contemporary workshop * sukanta.basu@gmail.com or conference session on SBLs. 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. [1], Basu and Porté-Agel [6] demonstrated that the GABLS1-LES case can be simulated reliably by a dynamic SGS parameterization 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 parameterization for reliable, tuning-free simulations of stable boundary layers, even with relatively coarse resolutions." [6] Later on, Stoll and Porté-Agel [7] and Lu and Porté-Agel [8] 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 parameterization. This parameterization was originally proposed by Deardorff [9] and appears in many LES codes [e.g., [10][11][12][13]. Henceforth, we refer to this parameterization as D80. Several past SBL-LES studies have reported grid-size sensitivity with the D80 SGS scheme [e.g., 1,5,14,15]. After extensive numerical experiments, we have identified the SGS mixing length (λ) parameterization 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 λ alleviates the grid-size sensitivity substantially. In addition, the first-and second-order statistics from the arXiv:2003.09463v1 [physics.ao-ph] 20 Mar 2020 LES runs utilizing this revised parameterization (named D80-R) agrees well with the ones produced by a pseudospectral 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 parameterization and its fundamental limitations. Technical details of the simulations are provided in Sect. 3. Results based on the D80 and D80-R SGS schemes are documented in Sect. 4. A few concluding remarks are made in Sect. 5. To further confirm the strength of the revised approach, simulated results from an independent LES code are included in Appendix 1.

II. SGS PARAMETERIZATIONS
The SGS parameterization (D80) by Deardorff [9] utilizes the following mixing length scale: where, ∆ = (∆x∆y∆z) 1/3 . L b is the so-called buoyancy length scale and is typically represented as follows: where, e denotes SGS turbulent kinetic energy (TKE). N is the Brunt-Väisäla frequency. Deardorff [9] assumed c n to be equal to 0.76. Many LES studies still assume this constant value [e.g., 10,12,14] though there are a handful of exceptions. For example, Gibbs and Fedorovich [13] assumed c n = 0.5. The eddy viscosity (K m ), eddy diffusivity (K h ), and energy dissipation rate (ε) are all assumed to be functions of λ as shown below: The unknown coefficients (c m , c h , and c ε ) are either prescribed or parameterized as follows: Please note that the values of the coefficients in Eq. 4 somewhat vary between different studies [e.g., 9,16,17]. For neutral condition, Eq. 1 reduces to: λ = ∆. As a direct consequence, K m , K h , and ε become stability independent, as would be physically expected. Furthermore, SGS Prandtl number (P r S = K m /K h = c m /c h ) equals to 0.33 for neutral condition.
For stably stratified condition, one expects that K m , K h , and ε should have a clear dependence on N . Deardorff [9] 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 λ = ∆, which fails to take account of the possibility that in a stably stratified region λ could become much smaller than the grid interval." [9] For the very stable case, in the limit of 1/N → 0, the D80 SGS model predicts λ = L b → 0, c ε → 0.19, and P r S → 1. Also, K m and K h are expected to approach negligibly small values for such conditions.
The behavior 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 SBLs. 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, λ becomes equal to ∆ 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 condition, K m , K h , and ε 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 parameterization. It is not influenced by the presence of a surface. Deardorff [9] recognized this problem and proposed a solution of increasing c ε near the surface [see also 18]. In this context, Gibbs and Fedorovich [13] 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., 10,12,15]. In other static SGS models (e.g., the Smagorinsky model and its variants), empirical wall functions are utilized to explicitly account for the near-surface shear effects [e.g., 19,20]. 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., 6,7].
In the present study, we replace the mixing length parameterization in D80 (i.e., Eq. 1) with the following formulation: where, κ is the von Kármán constant. This mixing length scale was originally proposed by Brost Wyngaard [21] for a Reynolds Averaged Navier-Stokes (RANS) model. It is still used in many present-day atmospheric models [e.g., 22]. Eq. 5 nicely captures both the stability and nearsurface shear effects. Please note that the grid-size (∆) does not directly impact the mixing length scale according to Eq. 5. Rather, e decreases with increasing resolution and impacts L b , λ, K m and K h , accordingly. In addition, the effect of ∆ is felt via c h , and c ε . Hereafter, we refer to Deardorff's SGS parameterization in conjunction with Eq. 5 as the D80-R parameterization.
In terms of SGS Prandtl number (P r S ), both the D80 and D80-R parameterizations suffer from unphysical prescription in different ways. This issue is discussed in detail in Sect. 4.

III. DESCRIPTION OF THE SIMULATIONS
In this work, we simulate the GABLS1-LES case study using the Dutch Atmospheric Large-Eddy Simulation (DALES) [10] and the PALM model system [12,23]. Since the configurations of the GABLS1-LES case study are well-known in the literature, we mention them in a succinct manner. The boundary layer is driven by an imposed, uniform geostrophic wind of 8 m s −1 , with a surface cooling rate of 0.25 K h −1 and attains a quasisteady state in about 8-9 h with a boundary layer 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 × 10 −4 s −1 , corresponding to latitude 73 • 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 × 400 m × 400 m. A wide range of isotropic ∆ 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 ∆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 time step of 0.075 s.
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 [24]. 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 [25], in order to apply MOST, the first model level in an LES model should not be at a height lower than ∼ 50z 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 λ depend on local N 2 values and are listed in Table I. 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.

IV. 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.
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 (λ) profiles. From Fig. 1 (top panel), it is evident that in the D80-based runs (left panel), λ equals to ∆ 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, λ 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 λ 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 λ.
In the D80-based runs, due to the underlying condition of λ = ∆, P r 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, λ is typically smaller than ∆ near the surface. Thus, P r S is much larger than 0.33. However, in Eq. 4b ∆ Eq. 4b Eq. 5 Eq. 4b ∆ PALM Eq. 4b Eq. 1* Eq. 4b ∆ c h = cm Eq. 5 Eq. 4b ∆ * PALM uses min(1.8z, ∆, L b ). the middle part of the SBL, depending on stratification, λ may become larger than ∆. For such cases, in D80-Rbased runs, P r S can be even smaller than 0.33 (bottomright panel of Fig. 1). In the upper part of the SBL, due to stronger stratification, λ becomes much smaller than ∆, and as a consequence, P r S monotonically approaches to unity for both the D80 and D80-R cases. In the case of the LASDD SGS model, the dynamically estimated P r S values remain more or less constant (around 0.5-0.6) inside the SBL and becomes greater than 1 in the inversion layer.
Gibbs and Fedorovich [13] recognized the problem with the P r S in D80 parameterization. They proposed to use P r S = 1 when N 2 is locally positive. For locally negative N 2 values they proposed an empirical formulation for P r S . In the present study, thus, the PALM model is employed with P r S = 1 for N 2 > 0 in all the D80-R-based runs. These results are shown in Appendix 1. Please note that for N 2 ≤ 0, both the DALES and the PALM models always use Eq. 4b (i.e., they effectively assume P r S = 0.33); we have not incorporated the empirical formulation by Gibbs and Fedorovich [13].
Simulated time series of surface friction velocity and sensible heat flux are documented in Fig. 2. In the D80based runs, no temporal fluctuations of surface fluxes are visible for ∆ ≥ 6.25 m. When ∆ 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 ∆ ≥ 6.25 m, λ (= ∆) 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 (∼ 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 ∆ ≥ 6.25 m. This unphysical behavior is solely due to the inherent limitations of the D80 SGS parameterization. Once ∆ 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 DALESbased simulated LLJ peak remains lower and weaker than the one simulated by the MATLES code. The positive impacts of the revised mixing length parameterization 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 ∆ = 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 [26] investigated the wind turning in SBLs in great details. 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 parameterization 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 strong agreement with one another. From the middle panel of Fig. 3, it is also clear that the original D80 parameterization performs quite poorly in capturing the turning of wind with height for coarse grid sizes.
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-sizes. Once again, incorporation of the revised SGS mixing length parameterization leads to better convergence. However, in this case, lower values of P r 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 MA-TLES code. In contrast, by using a higher value of P r S , the PALM model generates potential temperature profiles which are indistinguishable from the MATLES-based ones (refer to Fig. 11 in Appendix 1).
In the D80-based runs, with finer resolutions, the simulated mean potential temperature profiles also start to converge. However, for these cases, P r 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 (σ 2 θ ) do show the effect clearly (discussed later).
The u-and v-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 ∆ ≥ 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 in 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.
As depicted in Fig. 6, the overall behavior 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 ∆ = 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 and sensible heat flux profiles, the resolved TKE is non-existent in the bottom-half of the SBL for the D80-based simulations using ∆ ≥ 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-part (especially, in the lower part of the SBL). In other words, the flow is highly resolved  for all the simulations involving the D80-R parameterization.
Lastly, the resolved variances of vertical velocity and potential temperature are plotted in Fig. 8. The trends of the resolved σ 2 w profiles are similar to the resolved TKE plots. The resolved σ 2 θ profiles from the D80based runs are quite interesting. As discussed earlier, the coarse-resolution runs (i.e., ∆ ≥ 6.25 m) show the quasi-laminarization problem up to 75 m or so. Interestingly, for higher resolution runs with the original mixing length parameterization, the resolved σ 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 P r S (= 0.33) in the bottom part of the SBL. Similar decrease in P r S in the core of the SBL also causes resolved σ 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 behaviors [see 27, for some examples]. Spurious piling up of energy near the highwavenumber part of the spectrum is a telltale sign of under-dissipation. Since DALES and PALM are finitedifference 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.

V. DISCUSSION
In addition to Deardorff's SGS model, the Smagorinsky SGS model [28] is also quite popular in the boundary layer community. In this SGS model, the effective mixing length is C S ∆; where, C S is the so-called Smagorinsky coefficient. This coefficient is adjusted empirically [e.g.,  19,20,29] or dynamically [e.g., 6, 30, 31] to account for shear, stratification, and grid resolution. In contrast, the effective mixing length is c m λ in D80 SGS model. Most of the LES codes assume a constant value for c m and expects λ 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 parameterization by Brost and Wyngaard [21] 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 ∆-dependent and estimated via a dynamic approach [e.g., 32,33]. Our future work will be in that direction.
We will also need a better understanding of the SGS Prandtl number (P r S ) and its relationship to the turbulent Prandtl number (P r T ). Several field and laboratory experiments have reported that P r T should be on the order of one for stably stratified conditions [see 34,35, and the references therein]. However, within the SBL (excluding surface and inversion layers), the dynamic SGS models typically estimate P r S and P r T values to be around 0.5-0.6 and 0.7, respectively [6,7]. We agree with Li [35] that P r S is less important than P r T as the large-scale fluxes are resolved in LES. Nonetheless, we do believe the current P r S parameterizations in the D80 and D80-R SGS models are not acceptable and definitely need amendments. The formulation proposed by Gibbs and Fedorovich [13] is a good starting point and should be rigorously tested in future studies.

VI. CONCLUDING REMARKS
In the preface of the book titled "A History and Philosophy of Fluid Mechanics", G. A. Tokaty wrote: "The scientist should listen to every reasonable suggestion, but judge objectively. He should not be biased by appearances; have a favorite hypothesis; be of a fixed school of thought; or have a master in matters of knowledge. He should remember constantly that the progress of knowledge is often hampered by the tyrannical influence of dogma." In the present context, this quote seems rather appropriate. J. W. Deardorff was the father of large-eddy simulations of boundary layer flows. He made invaluable contributions to our field and we are all indebted to him to say the least. Perhaps, due to this adulation, it took us four decades to take a closer look at his SGS mixing length parameterization. Even though Deardorff proposed his parameterization for a stratocumulus-topped boundary layer, the boundary layer community decided to use it blindly for stable boundary layers without questioning its validity.
In this study, we have demonstrated that Deardorff's mixing length parameterization 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 SBLs, many open questions still remain: • Would this new mixing length parameterization work for very stable boundary layers? Would it allow us to simulate turbulent bursting events?
• Is the proposed parameterization under-dissipative near the surface? How can we (dis)prove this behavior?
• 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 [36][37][38]?
• How sensitive are the simulated results with respect to the chosen SGS coefficients (i.e., c n , c m , c ε )? Should they be dynamically determined instead of being prescribed? • How should we parameterize the SGS Prandtl number? Should it be a function of point-wise gradient Richardson number?
A few years ago, Basu and Lacser [25] cautioned against violating MOST in LES runs with very high resolutions. To overcome this issue, Maronga et al. [5] proposed certain innovative strategies; however, the results were somewhat inconclusive (possibly) due to the usage of the D80 mixing length parameterization 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.
We further recommend the SBL-LES community to revisit some of the past 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.

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 via zenodo.org. The simulated results from the PALM model are documented in this appendix. 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 ≤ ∆ ≤ 5 m, the results from ∆ = 1 m run stands out. We believe that λ is quite low for this particular case; even though it is small enough to sustain turbu-lence near the surface, it is not large enough to promote diffusion. The D80-R run with ∆ = 1 m does not portray such an unusual behavior.
The (positive) impacts of higher SGS Prandtl number (P r 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 σ 2 θ values are much larger than their DALES counterpart.