An Improved Surface Boundary Condition for Large-Eddy Simulations Based on Monin–Obukhov Similarity Theory: Evaluation and Consequences for Grid Convergence in Neutral and Stable Conditions

Monin–Obukhov similarity theory is used in large-eddy simulation (LES) models as a surface boundary condition to predict the surface shear stress and scalar fluxes based on the gradients between the surface and the first grid level above the surface. We outline deficiencies of this methodology, such as the systematical underestimation of the surface shear stress, and propose a modified boundary condition to correct for this issue. The proposed boundary condition is applied to a set of LES for both neutral and stable boundary layers with successively decreasing grid spacing. The results indicate that the proposed boundary condition reliably corrects the surface shear stress and the sensible heat flux, and improves grid convergence of these quantities. The LES data indicate improved grid convergence for the surface shear stress, more so than for the surface heat flux. This is either due to a limited performance of the Monin–Obukhov similarity functions or due to problems in the LES model in representing stable conditions. Furthermore, we find that the correction achieved using the proposed boundary condition does not lead to improved grid convergence of the wind-speed and temperature profiles. From this we conclude that the sensitivity of the wind-speed and temperature profiles in the LES model to the grid spacing is more likely related to under-resolved near-surface gradients and turbulent mixing at the boundary-layer top, to the SGS model formulation, and/or to numerical issues, and not to deficiencies due to the use of improper surface boundary conditions.


Introduction
One persisting problem in large-eddy simulation (LES) of the atmospheric boundary layer is the so-called logarithmic layer mismatch, the fact that the simulated wind-speed profile does not match the predicted logarithmic relation, a direct result of the inherent inability of LES models to resolve locally the dominant (small) eddies close to the surface. In this region, the subgrid-scale (SGS) turbulence parametrization dominates the flow. Also, the ability to resolve near-surface vertical gradients depends strongly on the grid spacing involved. It is commonly found that the mean wind shear is thus overestimated near the surface (Sullivan et al. 1994;Khanna and Brasseur 1997;Brasseur and Wei 2010).
The structure of the atmospheric surface layer can be essentially described through the turbulent exchange of momentum, heat, and moisture with the surface. Monin-Obukhov similarity theory (MOST, Monin and Obukhov 1954) provides a solid mathematical framework to describe this exchange in terms of turbulent fluxes and atmospheric stability. MOST includes the logarithmic law-of-the-wall in neutral conditions under the assumption of a constant-flux layer. In reality, however, the turbulent fluxes usually change with height from the surface value to zero at the top of the boundary layer. There is evidence from field and numerical experiments, though, that MOST provides reliable estimates of the surface fluxes and has been used for more than 50 years (Foken 2006). The vertical extent of the surface layer is commonly defined to be that region in which the turbulent fluxes do not vary more than ≈ 10% of their surface values. Due to the linear decrease with height (observed in steady-state conditions) the depth of the surface layer can be loosely estimated to be ≈ 10% of the boundary-layer depth.
To the authors' knowledge, MOST is nowadays used in most state-of-the-art LES models for atmospheric boundary-layer flows as a surface boundary condition to calculate the surface Reynolds stress and the turbulent surface sensible and latent heat fluxes (e.g., Heus et al. 2010;Maronga et al. 2015;van Heerwaaren et al. 2017). LES models typically use grid spacings at a metre scale, since in recent decades, due to continuously increasing computational power of state-of-the-art supercomputers, the grid spacing has decreased from about 100 m (Deardorff 1980) down to 1 m (e.g., Maronga and Bosveld 2017) or even less (e.g., Sullivan et al. 2016). Especially under stable conditions, the dominant eddies are often smaller than 10 m and thus demand very fine grids (Beare et al. 2006;Sullivan et al. 2016;Maronga and Bosveld 2017).
A problem arises from the systematic overestimation of the wind shear near the surface due to under-resolved flow (e.g., Nikitin et al. 2000;Kawai and Larsson 2012). As discussed above, atmospheric LES models apply MOST between the surface and the first grid level and thus precisely in the region of excessive wind shear, leading to a systematical underestimation of the surface shear stress. This can be explained by the fact that a region of excessive shear near the surface leads to too high wind speeds above this very region. As a direct consequence, the surface shear stress is too small to match the wind-speed profile in this upper region. This is illustrated in Fig. 1, where we see that the excessive shear near the surface compared to the MOST prediction leads to too high wind speeds above and that can only be compensated for by having a higher friction velocity in the MOST prediction. The same excessive vertical gradient was also observed for the near-surface temperature profile (e.g., Maronga and Bosveld 2017).
Other studies also showed that this excessive shear implies deviations from the expected MOST functions throughout the surface layer (Sullivan et al. 1994; Khanna and Brasseur 1997;Maronga and Bosveld 2017). On the one hand, one possibility to overcome this issue is to apply dynamic SGS schemes (e.g., Sullivan et al. 1994;Porte-Agel 2000;Porte-Agel et al. 2004;Basu and Porte-Agel 2006;Brasseur and Wei 2010;Lu and Porte-Agel 2013).  1 Schematic illustration of the excessive wind shear observed in LES models. Shown are exemplary mean profiles (denoted by <>; notation will be introduced later) of the horizontal wind speed u based on LES data and the respective MOST prediction based on the friction velocity calculated based on the wind shear between the surface and the first grid level z 1 . The region between the two profiles in which these do not run parallel is denoted as the excessive wind shear region On the other hand, Kawai and Larsson (2012) discussed that, theoretically, one can take any height level within the logarithmic region (i.e., the inertial sublayer) to evaluate the law-of-the-wall (or MOST to be more general). Based on this reasoning they proposed a correction for the neutral boundary layer by taking an elevated level for evaluating the lawof-the-wall in which the flow-and thus the wind shear-is fully resolved. They showed that the resulting wind-speed profiles displayed a significantly improved agreement with the logarithmic profile predicted by theory. However their method was essentially designed for engineering LES application and thus did not incorporate atmospheric stability. Also, they did not discuss whether and how their approach compromised the correlation between the near-surface turbulence structure and the local surface fluxes. Recently, Maronga and Bosveld (2017) employed a similar method for stable conditions in LES of a nocturnal radiation fog, where the simulated fog-formation time showed a decisive dependence on the grid spacing. By taking a fixed elevated level for evaluating the MOST functions they improved the grid convergence regarding fog formation significantly. Grid convergence here means that the simulation results in general (and the quantities to be studied) do not change if the grid spacing is further reduced. This is a general prerequisite for all kinds of LES applications.
The effect of grid resolution on simulation results has been studied by various authors (Sullivan et al. 1994;Khanna and Brasseur 1997;Beare et al. 2006;Brasseur and Wei 2010;Maronga and Bosveld 2017), and might be linked to the shortcomings outlined above (Brasseur and Wei 2010). 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.36 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. From our own experience, and in line with previous research, we found a non-convergence of the surface shear stress and heat flux for both neutral and stable boundary layers, suggesting that the non-convergence in LES models in the stable boundary layer might be related to the outlined issues of the surface boundary condition. The motivation for the present study thus was to, (a) develop a reliable methodology (i.e., a surface boundary condition for LES models) to correct the surface fluxes due to excessive wind-speed and temperature gradients close to the surface, and (b) assess the chosen surface boundary condition as a possible reason for the lack of grid convergence in neutral and particularly in stable conditions. It is known that the application of MOST as a boundary condition imposes further limitations and problems. First, when using grid spacings in the order of a few metres, the application of MOST can become problematic. Several researchers pointed out that MOST can only be considered valid within the inertial sublayer, the upper part of the surface layer, but fails in the roughness sublayer below, in which direct effects of single surface roughness elements are present (Garratt 1980;Raupach 1992;Physick and Garratt 1995;Harman and Finnigan 2007;Basu and Lacser 2017). Garratt (1980) found, for high vegetation and neutral and unstable conditions, that the lower boundary of the inertial sublayer (z * ), i.e., the top of the roughness sublayer, depends on the average horizontal distance between trees (δ), and estimated z * ≈ 3δ, which corresponded to z * = 35z 0 , where z 0 is the roughness length. Physick and Garratt (1995) noted that z * might be much smaller for stable conditions. For a similar experimental site reported by Physick and Garratt (1995), z * was found to be around 50z 0 instead. Note that these values are site and stratification specific and are taken as some first proxy at this point. Garratt (1980) and Harman and Finnigan (2007) suggested corrections for the MOST relationships below z * in order to include the roughness sublayer. This correction, however, was derived for a specific forest canopy and for neutral atmospheric stratification only. Referring to this previous work, Basu and Lacser (2017) suggested considering z * = 50z 0 as a general rule for LES models. As an example, following this suggestion and for typical roughness lengths for surfaces covered with low vegetation (heights in the order of 0.1 m), this imposes a minimum vertical grid spacing of 5 m for LES models with a common MOST boundary condition. Basu and Lacser (2017) reported correctly that, in practice, MOST also is applied for much smaller grid spacings, despite this violation of MOST assumptions and in the lack of an alternative (e.g., Beare et al. 2006;Basu et al. 2011;Maronga 2014;Sullivan et al. 2016;Udina et al. 2016). This is particularly true for LES of the stable boundary layer, where grid spacings < 5 m are required to resolve the small-scale turbulence (Beare and MacVean 2004). The violation of MOST assumptions for very fine grid spacings can thus be considered one possible reason for the insufficient grid convergence observed in LES of the stable boundary layer. Contemporary atmospheric LES codes nevertheless apply MOST as a boundary condition between the surface and the first computational grid level above the surface, providing a further motivation for reviewing the possible implications of this practice.
The starting hypothesis for the present paper is that the lack of grid convergence observed in LES of neutral and stable boundary layers might be attributed to the outlined problems when using MOST as the surface boundary condition. Consequently, we extend the work of Kawai and Larsson (2012) and Maronga and Bosveld (2017) through a proper description and a thorough evaluation of an improved surface boundary condition for LES models for neutral and stable stability regimes based on MOST. The improved boundary condition is designed in such a way that the MOST assumptions are not violated and the surface fluxes are corrected to fit to the resolved profiles of wind speed and temperature. For this purpose we employ the LES model PALM (Maronga et al. 2015) and conduct a set of idealized LES of neutral and stable boundary layers.
The paper is organized as follows: Sect. 2 gives an overview of the current and the proposed surface boundary conditions for LES models. Section 3 describes the LES set-up, while results are presented in Sect. 4. Finally, a summary and outlook is presented in Sect. 5.

Methodology
In the following we first outline the state-of-the-art methods applied in LES models to calculate the surface fluxes of heat and momentum (Reynolds stress), and confine ourselves to the case of a dry atmosphere without the presence of humidity. The methodology can, however, be easily extended to incorporate humid air and hence the surface latent heat flux. In a second step, we describe an improved method that accounts for the excessive wind shear and temperature gradient near the surface, while simultaneously conserving the correlation between local surface fluxes and the flow adjacent to the surface.

Traditional Surface Boundary Condition
In the MOST framework, the wind-speed profile can be expressed as with u h being the streamwise horizontal wind speed and z being the height above the ground; κ = 0.4 is the von Kármán constant. A tilde ( ) symbol indicates the ensemble average. The similarity function for momentum φ m depends on the stability parameter z/L, where L is the Obukhov length defined as where θ is potential temperature, g is the acceleration due to gravity, and a prime ( ) symbol indicates a turbulent fluctuation. A subscript 0 indicates a surface value. The friction velocity u * , is defined through the shear stress at the surface τ 0 (also referred to as the surface Reynolds stress), with u, v, and w being the velocity components in x, y, and z directions, respectively. This framework is commonly applied in LES models to predict u * by rearranging Eq. 1 and integration over z from the roughness length z 0 to the height of the first computational grid level above the surface (z 1 ). This results in a diagnostic equation for u * , where Here we have also introduced the dependence on x and y to account for the fact that a surface stress (or friction velocity) is needed at each surface element (x, y) of the computational domain. For convenience we omit this dependence hereafter. The function ψ m is defined as (see e.g., Panofsky 1963;Grachev et al. 2007) and follow the notation of Hultmark et al. (2013) and refer to this method as the "instantaneous logarithm" (IL) method. As discussed by Hultmark et al. (2013), this method leads to a systematical overprediction of the mean shear stress by (κ/ψ M ) 2 < u h > 2 (cf. Eq. 4). Note that hereafter we use the horizontal average (<>) over the entire model domain and additional time averaging, indicated by an overbar, instead of the ensemble average (˜), which is the LES variant of an ensemble average. Note that we use time averaging only during analysis of the results, but not in the realization of the used boundary conditions. The overestimation is caused by the fact that the equations are defined for mean quantities and the non-linearity in the logarithmic law. An analysis of LES data at hand suggest that this error is in the order of 1% as long as the mean wind speed does not tend to zero (i.e., in convective cases). In such cases, however, buoyancy usually dominates the flow and the surface shear stress becomes a minor contributor to the surface energy exchange. One common method that avoids this overestimation was developed by Schumann (1975) and was improved by Grötzbach (1987) (the so-called SG method, named after Schumann and Grötzbach, see Hultmark et al. 2013). This method is based on solving the averaged form of the equations above and imposing a local variation based on the ratio u h (z 1 )/< u h >(z 1 ). For a detailed description of the SG method and more elaborate methods, see Piomelli et al. (1989), Marusic et al. (2001), Stoll and Porte-Agel (2006), and Hultmark et al. (2013). While these methods might improve the boundary condition for purely neutral flows, they do not consider the effect of atmospheric stability. Also, their application is limited to very idealized cases as the surface is required (for most methods) to be entirely homogeneous. Many LES models, including the model PALM used in the present study, employ the IL method, which is purely local.
The derivations above can similarly be made for calculating the surface flux of sensible heat w θ 0 , which is defined through the temperature scale θ * , and which can be calculated as and the function ψ h given by, In order to solve Eqs. 4 and 8, it is required to obtain L, which is a function of both u * and θ * . In the PALM model this is realized using a Newton iteration method and the relationship with the bulk Richardson number Ri b being defined as Equation 11 is first solved for L with Ri b calculated via Eq. 12. Afterwards, u * , θ * , and the fluxes are calculated based on the obtained L.

Proposed Surface Boundary Condition
At this point, we first make use of the reasoning of Kawai and Larsson (2012) that there is no requirement to apply MOST between the first grid level and the surface. It thus appears most logical to substitute z 1 with an elevated height z sl that fulfils the requirements to, (a) be located in the surface layer, and (b) be far enough from the surface where the bulk part of the turbulent transport is resolved so that the influence by the SGS model is negligible. However, as the IL method uses local information only, this would induce a spatial correlation between the wind-speed and temperature fields at height z sl and the surface fluxes. By the same token, however, the surface fluxes then would not be spatially correlated to the adjacent flow, which leads to an undesired and non-physical behaviour. In order to avoid this, we make use of the concept of the SG method and use only the horizontally-averaged quantities at z sl , and a modulation according to the local quantities at z 1 . This then yields and Note that the last bracket term acts as a correction to maintain the correlation between the quantities at the surface and at z 1 . At this point we need to stress that the turbulence at the first couple of grid points above the surface is usually not well resolved in LES models, and neither are the small structures. The number of grid points in question is model-dependent (e.g., seven in the present study, see below). Alternatively, we might use the local information at height z sl instead of its horizontal average; however, the turbulence spectrum at that height is neither a good approximation for the flow adjacent to or directly above the surface, because (a) the dominant eddies scale with distance to the surface so that the maximum of the spectrum shifts towards larger scales compared to the near-surface flow, and (b) eddies have already travelled according to the mean wind speed and are thus not correlated to the surface directly below. Note that the latter is the reason why (Piomelli et al. 1989;Stoll and Porte-Agel 2006) used a displacement factor and why (Bou-Zeid et al. 2005) used a filtering procedure in their studies. This is further supported by the fact that the viscous force near the surface is large and that we should not expect the momentum cospectra to follow those of the airflow. Experimental studies for neutral conditions have shown that the cospectra display no inertial subrange and thus implicate only a small separation between large and small scales. It is thus physically questionable to use the wind-speed and temperature data from an elevated height. Equations 13-14, on the other hand, maintain the correlation between the adjacent flow spectra and the surface and are thus in line with the traditional methods, but might be flawed by the under-resolved spatial turbulent fluctuations near the surface. In this framework, Ri b is redefined and approximated as As for the IL method, the above set of equations can be solved iteratively for L, and hereafter, we refer to this method as the "elevated SG" (ESG) method. For the special case, when z sl = z 1 is used, note that we obtain the SG method. We thus focus our analysis on a comparison of the IL and ESG methods, but evaluate the effect of using the horizontal average using selected large-eddy simulations with the SG method.

Requirements for the Evaluation Height
The ESG method requires prescribing the height for evaluating MOST (z sl ) explicitly. The choice is here limited by the following requirements: 1. This level must be within the atmospheric surface layer (also referred to as constant-flux layer). The surface-layer height H can be estimated to an extent over the lowest 10% of the boundary-layer height z i , As several definitions exist for defining z i , which also differ for different stability regimes, this remains a rather loose criterion. The chosen height should thus be as close to the surface as possible. 2. The level must be within the inertial sublayer and thus above the roughness sublayer.
Here, we follow (Basu and Lacser 2017) and demand that Note that the coefficient of 50 here is the recommendation of Basu and Lacser (2017), but that it might depend on stability and the characteristics of the roughness elements (see Garratt 1980;Physick and Garratt 1995). Also, note that this requirement could be eliminated by adding a correction term for the roughness sub-layer as discussed, e.g., by Harman and Finnigan (2007). 3. The flow at height z sl must be well-resolved by the model, a requirement that is not a physical constraint, but a numerical requirement. The height z sl is thus highly dependent on the model (including its SGS scheme). For the PALM model (static Deardorff SGS model), previous studies have found that the mean profiles follow MOST at height levels starting from the seventh grid level above the surface (Maronga 2014;Maronga and Reuder 2017). As the PALM model uses a Cartesian staggered grid with wind speeds and scalars defined at the vertical centre of the grid boxes, the requirement for the PALM model is where Δz is the grid spacing in the z-direction. Note that this height is model-specific and will also depend on the numerical schemes involved. The value reported here is valid for PALM's 5th-order advection scheme (Wicker and Skamarock 2002).
For given values of z 0 and z i , the grid spacing can no longer be chosen freely to meet these requirements and a maximum grid spacing is predefined. For example, for a neutral boundary layer with z 0 = 0.1 m and z i = 1000 m, z sl ∈ [5 m, 100 m] and hence the allowed grid spacing must be Δ z ≤ 15.4 m. However, this value will further decrease when choosing a level closer to the surface as suggested by requirement 1. Analogously, a stable boundary layer with identical roughness but a typical height of only 100 m implies that Δ z ≤ 1.5 m, without considering the dominant size of the eddies at all.

Experimental Set-Up
The LES model PALM 5.0 in revision 3230 was used for the present study, which solves the equations of conservation for momentum, heat, and moisture in Boussinesq-approximated form on a Cartesian staggered Arakawa-C grid. It has been widely applied for different flow regimes in the convective (e.g., Raasch and Franke 2011;Maronga and Reuder 2017), neutral (e.g., Knigge et al. 2015;Knigge and Raasch 2016) and under stable conditions (Beare et al. 2006;Maronga and Bosveld 2017). All simulations were carried out using cyclic lateral boundary conditions. The PALM model applies an SGS turbulence closure scheme after (Deardorff 1980) in the formulation of Saiki et al. (2000), which solves a prognostic equation for the SGS turbulence kinetic energy (TKE) with due regard for third-order TKE turbulent diffusion and applies simple down-gradient formulations for the components of the Reynolds stress and scalar fluxes. For selected runs, a dynamic SGS closure after (Heinz 2008;Mokhtarpoor and Heinz 2017) was employed instead. Discretization in space and time is achieved by a fifth-order advection scheme after (Wicker and Skamarock 2002) and a thirdorder Runge-Kutta time-stepping scheme (Williamson 1980). A geostrophic wind vector from the west-and hence along the x-direction of the model domain-was prescribed in all cases to maintain dynamically-generated turbulence. Note that the was no perfect alignment of the flow within the boundary layer as the Coriolis force led to the flow veering with height. Monin-Obukhov similarity theory is applied as a surface boundary condition using either the IL or the ESG method, both employing the similarity functions of Businger-Dyer (see e.g., Panofsky and Dutton 1984) for momentum, and heat: For details on the technical realization of the MOST boundary condition in the PALM model and the implementation of the SGS model, see Maronga et al. (2015). Two sets of simulations for a neutral and a stable boundary layer were performed and are outlined in the following.

Set-Up of Neutral Boundary Layer
First, an idealized neutral boundary layer with a geostrophic wind speed of 5 m s −1 and z 0 = 0.1 m was chosen to test the new boundary condition without taking into account effects of static stability. The horizontal domain size for all runs was about 2000 m × 2000 m, which is in the range of domain sizes found to be sufficient for simulating neutral boundary layers (see Andren et al. 1994;Drobinski and Foster 2003;Foster et al. 2006;Ludwig et al. 2009;Knigge and Raasch 2016). The Coriolis parameter was set to a latitude of 73 • N (to be consistent with the runs for stable conditions, see below). The wind speed was initialized by its geostrophic value and was constant with height. Initial perturbations were imposed on the velocity fields to trigger turbulence. Overall, 12 LES were performed, four using the IL (cases IL_dX ) and four the ESG method (cases ESG_dX , where X indicates the grid spacing used). Additionally, we conducted four runs with z sl = z 1 (i.e., SG method) in Cases not fulfilling all requirements are marked with * . Additionally, the computational demand is given in processor hours (CPUh) on the used Cray XC40 cluster equipped with Intel Haswell processors order to separate the effect of taking an elevated level and the outlined general differences between IL and SG methods (cases SG_dX ). For each method, the equidistant grid spacing (Δ = Δ x = Δ y = Δ z ) was varied from 16 to 8 m, 4 m, and 2 m. For the IL method, z sl was defined by the numerical grid and thus set to 0.5Δ z = z 1 , while for the ESG method z sl was chosen to be the height of the first grid level greater or equal 50 m. As z i varied between the runs from 880-1100 m, this was well within the inertial sublayer (equal to 0.05z i ). A known issue for simulations of the neutral boundary layer are inertial oscillations in the wind profile. By carefully looking at the time series of < u * >, we found that these oscillations are sufficiently damped throughout the boundary layer after 60 h of simulation time in all cases. In particular, we found that these oscillations are very weak in the lower boundary layer and thus they did not affect the results within the surface layer. The simulations all ran for 84 h and the data were averaged and analyzed during the last 24 h (14,400 samples). Note that the overbar symbol is used in the following to denote time-averaged quantities. Table 1 gives an overview of the simulations performed for neutral conditions and whether the individual runs fulfil the requirements defined in Sect. 2.2.1. Note that all cases with IL and SG methods by implication fail to fulfil these requirements, while for the ESG method only case ESG_d16, i.e. with the coarsest grid, does not fulfil all requirements. The height of the surface layer was estimated to be H = 0.1z i , where z i was calculated as in Beare et al. (2006), as the height at which the mean shear stress fell to 5% of its surface value.
For test purposes, we also repeated the four IL cases with a dynamic SGS model, which was implemented in the PALM model during the review process of the present paper. Thereby we wanted to check whether an advanced SGS model might already lead to improved grid convergence. Results from these runs showed neither a significant advantage over the Deardorff SGS scheme, nor improved convergence, and thus will not be discussed further (as a consequence they are also not listed in Table 1).

Set-Up of Stable Boundary Layer
A weakly-stable boundary layer was used to evaluate the performance of the proposed boundary condition, with the set-up similar to that used in the GABLS1 model inter-comparison outlined in Beare et al. (2006). The initial potential temperature was set constant with height throughout the model domain to 265 K in our simulations (i.e., a stable boundary layer (SBL) developed in a neutrally-stratified environment). Note that in GABLS1, a capping inversion starting at a height of 100 m was used instead. A model domain of about 500 m×500 m×500 m was used, with equidistant grid spacing, and which is sufficient for a SBL that reaches up to ≈ 200 m. Geographical latitude was set to 73 • N as in GABLS1. In order to stimulate turbulence, a random potential temperature perturbation of amplitude 0.1 K was applied for height levels of up to 50 m. The vertical velocity was set to zero at the surface and top boundary of the domain, with top boundary conditions set to free-slip conditions. A continuous surface cooling rate of 0.25 K h −1 was applied to create and maintain a stably-stratified flow. The geostrophic wind speed was set to 8 m s −1 and the roughness was set as in the neutral simulations to z 0 = z 0h = 0.1 m. As the simulated SBL was significantly shallower than for the neutral cases, i.e. z i = 150−170 m, we followed a more rigid concept and fixed z sl to exactly 13 m. Unlike for neutral conditions, where a height in the middle of the surface layer was employed, we decided to use a height value very close to the top of the surface layer. This was done in order to use as coarse grid spacings as possible while simultaneously fulfilling all requirements for the ESG method. The grid spacings were then varied in such a way that a prognostic level was exactly located at a height of 13 m. Five cases were hence performed, each using the IL and the ESG method with grid spacings of 5.2 m, 2.8 m, 1.37 m, and 1.04 m. Note that the choice of z sl is always case-specific and must be adjusted depending on simulation set-up.
The simulations ran for 18 h, which is somewhat longer than in the GABLS1 case. As for neutral conditions, we observed inertial oscillations in the neutral layer above the boundary layer, which did not penetrate into the boundary layer itself. As the surface was continuously cooled, no true steady temperature state could be achieved, but we observed convergence of time series (e.g., friction velocity) after about 16 h. As direct consequence, a shorter timeaveraging of 1 h was applied as in Beare et al. (2006) from hours 17 to 18 (600 samples). We also tested other averaging intervals and earlier averaging periods, but found no effects on the results. Table 2 summarizes the cases performed for stable conditions and their key parameters. In analogy to the neutral boundary-layer experiments, we performed two additional set of runs using the SG method and a dynamic SGS model. As the SG cases for the SBL did not provide any new insight, we will omit a discussion of their results (as a consequence they are also not listed in Table 2). Results for the cases with dynamic SGS model (cases DYN_dX _sbl) are briefly discussed in the Appendix.

Results
In order to evaluate the two different applied boundary conditions, we first compare the results of IL and ESG methods for purely neutral conditions and thus neglect the possible complications imposed by stratification. In the second part of this section, we evaluate the performance and implications for stable conditions. Cases not fulfilling all requirements are marked with *

Neutral Boundary Layer
In neutral conditions, the boundary condition is reduced to the prediction of the surface shear stress and thus provides an ideal prerequisite for analyzing the behaviour of LES using different methods to determine this stress.

Logarithmic-Layer Mismatch
Dimensionless vertical profiles of the time-and horizontally-averaged horizontal velocity component from the cases using the IL method are shown in Fig. 2 (black lines). The focus at this point will not be on grid convergence. Nevertheless, it makes sense to look at the effect of the chosen boundary condition method for different grid spacings and thus differently wellrepresented surface layers in the simulation. In Fig. 2 there is clear evidence that the wind speed is higher in the bulk of the surface layer than MOST predicts for the given value of u * . Also, it is visible that the wind speed profiles strictly follow the predicted value by MOST near the surface, a direct implication by the IL method. Moreover, the difference between the actual wind-speed profiles and the predicted wind-speed profiles remains almost constant throughout the surface layer. This means that, despite the fact that the actual wind speed is too high, the vertical gradient of the wind speed agrees well with the theoretical gradient predicted by MOST. We can conclude from these observations that the excessive wind shear near the surface leads to a too low u * (with respect to the resolved wind-speed profiles above) when calculated with the IL method, independent of the grid spacing involved. Furthermore, we can identify a tendency of the discrepancy between the theoretical profiles and the actual LES profiles to increase for decreasing grid spacing. In order to evaluate the effect of using the horizontally-averaged wind speed instead of the local value (IL method) separately, we additionally plotted the SG method data in Fig. 2 (red lines). As expected, the behaviour is similar to the IL method close to the surface (first grid point). At higher levels, the SG method is closer to the theoretical MOST profile than the IL method, which is indicating higher u * values compared to the IL method runs. We might ascribe this to the fact that the horizontal mean gradient is used for calculating the surface shear stress instead of the local gradient in the IL method so that the excessive wind shear near the surface is smaller. However, note that the overestimation is still well pronounced (see Fig. 2), so that the SG method alone does not resolve the logarithmic layer mismatch. Figure 2 (blue lines) shows the same data for the ESG method. Here we immediately recognize that the wind speed profiles agree remarkably well to the predicted ones and are also much better than for the SG method. This result is in agreement with Kawai and Larsson (2012), who showed similar plots in their study of a neutral boundary layer configuration. The main reason for this is the higher value of u * (and thus < u * >) which was calculated to match the wind speed at a height around 50 m instead of the first grid level above surface. Furthermore, it is noteworthy that for coarser grids (cases ESG_d16 and ESG_d8) we observe the expected overestimation of the wind speed near the surface, being a direct implication of the correction involved in the ESG method, which provokes a shift of the dimensionless wind-speed profiles towards smaller values (as the scaling parameter < u * > is higher). Rather surprising is the fact that this overestimation is not visible for finer grid spacings (cases ESG_d4 and ESG_d2). Here, the ESG method apparently is able to also correct the excessive wind shear near the surface. Also, it should be noted that the ESG method also seems to provide a good correction in cases when the surface layer is not well-resolved and the criteria for the application of MOST as a boundary condition are not all met (i.e., case ESG_d16).

Effect on Grid Convergence
In order to evaluate whether the ESG boundary conditions lead to an improved convergence of u * and thus to the surface forcing via shear stress, statistics for u * were calculated and are shown for both IL and ESG methods in Fig. 3. For the IL method, we see that < u * > decreases continuously with decreasing grid spacings from values around 0.24 m s −1 to 0.23 m s −1 (decrease by 4.5%). This trend holds for the mean as well as the minimum and maximum values, while the standard deviation has a minimum at a grid spacing of 8 m (case IL_d8), probably due to changes in the flow structures related to the different resolution (e.g., streaks and roll-like patterns related to dynamic instabilities). However, there seems to be a trend to reach grid convergence for the mean value at grid spacings between 4 and 2 m, while the minimum and maximum values suggest a continuous trend towards smaller values.
The SG method, though showing consistently slightly larger values than the IL method, apparently does not improve grid convergence of < u * >. For the cases using the ESG method, however, we note that < u * > is generally higher than for the IL and SG methods, an expected consequence of the excessive wind shear leading to smaller values in the IL method. The data also reveal a weak decrease in the mean value with decreasing grid spacing from 0.26 m s −1 in case ESG_d16 to 0.25 m s −1 in case ESG_d4. This corresponds to a decrease by 2.3%, which is thus only half of the decrease observed for the IL method, indicating better convergence for the ESG method. However, we also observe a slight increase in the friction velocity from case ESG_d4 to ESG_d2, which might be suggestive to approach grid convergence. Note that the minimum and maximum values have converged, though. We thus ascribe the small increase in < u * > from case ESG_d4 to ESG_d2 to dynamic instabilities in the flow that temporarily alter the mean state of the boundary layer (streaks or roll-like structure, see above), but which are probably not resolved at coarser grid spacings. In general, we can conclude from Fig. 3 that the grid convergence of the surface shear stress is improved when using the ESG method.
In addition to Fig. 3, the spatial probability density functions (p.d.f.) of u * − < u * > are shown for an instantaneous point in time in Fig. 4, while Table 3 provides mean, standard deviation, skewness, and kurtosis for u * of the same data set. First of all, Table 3 reflects what   we learned from Fig. 3. The mean u * for the ESG method is higher than for the IL method, while the changes in the mean between individual runs are significantly smaller. Furthermore, the p.d.f. reveal that u * displays a narrower distribution for the ESG method than for the IL method. All cases display a positive skewness (towards larger values of u * ), which is in line with the results of Stoll and Porte-Agel (2006). The skewness is higher, though, for the IL method than for the ESG method, where u * appears to be almost of Gaussian distribution. This can be explained by the fact that u * is proportional to u 2 h and u h for the IL and ESG methods, respectively (see also Stoll and Porte-Agel 2006;Hultmark et al. 2013). Following the reasoning of Stoll and Porte-Agel (2006) this might lead to more damping of velocity fluctuations for LES using the ESG method. It is also visible that the maximum p.d.f. values for the ESG method coincide with their mean, while for the IL method we see that the peak p.d.f. values are smaller than their mean. The kurtosis shows values around 3 for all cases, indicating that the distributions are prone to outliers similar to Gaussian distributions. However, it is apparent that for the ESG method the kurtosis is consistently below 3, suggesting that the distribution here has a tendency to have higher flatness than a Gaussian distribution. Stoll and Porte-Agel (2006) calculated the kurtosis of the surface shear stress for several boundary condition methods and reported values greater than 3 (values around 3.5 were reported for the SG method), while our experiment suggests smaller values of at most 3.23.
In analogy to Marusic et al. (2001), we also calculated the spectra of u * in the x-direction, which can be considered to be the streamwise component (Fig. 5). Note however, that the nearsurface mean flow is affected by the Coriolis force so that it also has a small y-component (the wind direction near the surface was found to be around 250 • ). Here, we note that all spectra are qualitatively similar and do not show modifications in their shape between the two methods. The spectra reveal that the large scales follow the proposed k −1 behaviour for large scales under neutral conditions (see Stoll and Porte-Agel 2006) and a negligible inertial subrange between large and small scales, suggesting a small-scale separation, which is expected close to the wall where the viscous force is large (Hultmark et al. 2013). Qualitatively, these spectra compare remarkably well to measurements of the surface shear stress shown by Marusic et al. (2001) and Hultmark et al. (2013). The latter also performed LES, but their spectra displayed a much more pronounced inertial subrange, leading to a significant underestimation of the intermediate wavenumber range compared to observations. We can thus conclude that both methods can represent the expected spectral characteristics of u * very well and the ESG is able to retain these characteristics. Besides, we note a general difference in the energy level, where the runs with ESG method yield smaller energy levels than for IL method. This is a consequence of the narrower distribution and smaller standard deviation for the ESG method (see Table 3). Figure 6 shows the vertical profiles of the horizontal wind speed for both methods and all grid spacings employed throughout the boundary layer. First, we note that the wind speed profiles appear to be converged for the IL method at a grid spacing of 4 m, while no convergence is reached when using the ESG method. This appears initially as rather surprising given the better convergence observed in u * for the ESG cases. Assuming that this behaviour is solely an effect of the differences in u * , one would expect hardly any difference between cases ESG_d4 and ESG_d2, but a pronounced difference between IL_d4 and IL_d2. Figure 6 reveals, however, that there is no correlation between the behaviour of u * and the convergence of the mean wind-speed profile. These findings let us conclude that the ESG method, despite ensuring a much better near-surface wind shear, does not improve grid convergence of the mean wind-speed profiles. This might be ascribed to the fact that even with very fine grid spacings the LES technique is prone to erroneous flow near the surface, which is dominated by the SGS model. From Fig. 2 (ESG method, red lines) we learned, on the one hand, that the near-surface wind-speed profile did show deviations from the theoretical profile near the surface for cases ESG_d16 and ESG_d8, which indicates that these grid spacings are too coarse and must be disqualified for looking at grid convergence because the near-surface profiles are not well-resolved. On the other hand, both cases ESG_d4 and ESG_d2 no longer displayed this deficiency. Furthermore, note that we generally observed a slightly deeper Fig. 6 Mean profiles of the horizontal wind speed for, a cases with IL method, and b cases with ESG method. The profiles were averaged over 24 h of simulation time boundary layer for the ESG method, which is indicated by a vertical shift of the maximum wind speeds.
Profiles of the horizontal variances are shown in Fig. 7. While we find that grid convergence has not improved when using the ESG method, the profiles reveal that there is consistently more variance than for the IL method, indicating that turbulent mixing is stronger (variance of vertical velocity and the momentum fluxes, not shown, revealed the same behaviour). This agrees with the observations of a deeper boundary layer induced by increased surface friction that was found in the mean profiles (see Fig. 6). Also, we see that the variance is increasing for decreasing grid spacings, which is opposite of what we would expect from the decrease of < u * > with decreasing grid spacings. This suggests that the effect of finer grid spacing and therewith generally better resolved temperature gradients and turbulence has a dominant effect on the mean profiles and variances over the effect of variations in the surface friction induced by different grid spacings.
To summarize the results for neutral conditions, it was shown that the ESG method corrects the calculated surface shear stress efficiently for the exceeding wind shear near the surface. For the cases studied, however, grid convergence of the wind-speed profiles was not improved, which suggests that the finest grid spacing of 2 m was possibly not sufficient to resolve the near-surface wind-speed profile good enough. Also, deficiencies of the SGS model might play an important role here. Nevertheless, the overall differences in < u * > between the runs with different grid spacings were not found to be larger than 2.3%, which is smaller than the 4.5% observed for the IL method, and which appears tolerable for LES applications.

Stable Boundary Layer
The SBL provides a much more complex test bed for the evaluation of the ESG method, mainly because the Obukhov length now is an additional scaling parameter that accounts for the combined effect of shear production and buoyancy. The Obukhov length is used in the Fig. 7 Variance profiles of the horizontal wind components for, a cases with the IL method, and b cases with the ESG method. The profiles were averaged over 24 h of simulation time similarity functions to correct the logarithmic profile. Furthermore, not only is a boundary condition required for u * , but one is also required for θ * . In the next sections we will generally follow the same strategy as for the neutral boundary layer set-up and first evaluate the IL against the ESG method for linking u * and θ * to the surface layer profiles of wind speed and temperature. In a second step we will then discuss the issue of grid convergence. Figure 8 shows, in analogy to Fig. 2, the dimensionless profiles of the horizontal wind speed for the stable cases. The results reveal the same behaviour as observed for neutral conditions (cf. Fig. 2) where the wind speed is underestimated by the theoretical profile using the calculated u * that was used as a boundary condition. Also, the tendency of larger differences for decreasing grid spacings is recovered. Figure 8 (blue lines), however, reveals slightly different results than for neutral conditions. For stable conditions, we first note that the ESG method is (as for neutral conditions) effective in closing the gap between theoretical and LES profiles (see differences between solid and dashed lines in Fig. 8) and thus provides a good correction for u * . However, here we also note that there is a continuous tendency in the wind-speed profiles. For cases ESG_d5_sbl and ESG_d3_sbl, the LES profiles suggest lower wind speeds than predicted by MOST, while for cases ESG_d15_sbl and ESG_d1_sbl, the LES profiles more and more overestimates the theoretical profile. Before we try to find an explanation for this sensitivity to the grid spacing, we first analyze the temperature profiles.

Logarithmic-Layer Mismatch
In Fig. 9 (black lines) the dimensionless profiles of potential temperature are shown (IL method). As expected, we also note a systematical discrepancy between the theoretical profile and the LES data. However, unlike for the wind-speed profile, we observe generally too low temperatures (or too high θ * ). Also, Fig. 9 reveals that this discrepancy becomes slightly smaller for smaller grid spacings. The ESG method also seems to be able to correct reliably for this discrepancy. For the bulk part of the surface layer, the temperature profiles are now close to the predicted profiles by MOST. However, we must also recognize that the ESG method is not able to completely remove the systematical underestimation of temperature. Here we also see a different trend than for the IL method. The discrepancy for the ESG method apparently becomes slightly, but continuously, larger with decreasing grid spacing. In searching for an explanation for the comparably worse (compared to neutral conditions) representation of the surface-layer wind speed and temperature profiles under stable conditions, we varied z sl and L in post-processing in order to produce better agreement with the LES-based profiles. This experiment indicated, though, that it was not possible to find better-matching theoretical profiles with a single set of parameters. From this, we must conclude that either the LES profiles are deficient (see discussion below), or that the used MOST formulation for stable conditions is not able to represent these profiles correctly. In the next section we will discuss 316 B. Maronga et al. Fig. 8, but for temperature profiles for the SBL cases a possible reason for the non-convergence of θ * and its implication on grid convergence of the mean profiles.

Effect on Grid Convergence
In the previous section, we identified two sensitivities when using the ESG method. First, the wind-speed profiles showed an increasing overestimation of the theoretical profile (or in other words, a too high u * ) with decreasing grid spacing. Second, the temperature profiles showed the opposite tendency, underestimating the theoretical ones with decreasing grid spacing (i.e., too low θ * ). The respective statistics for u * and θ * are shown in Fig. 10a, b. As for neutral stratification, we note that u * is about 23% higher for the ESG method when compared to the IL method, an implication from the correction of the exceeding wind speed near the surface. Besides, it is evident that there is much better grid convergence for the ESG method compared to the IL method. In contrast, θ * reveals that for both methods grid convergence is lacking and θ * is monotonically decreasing with decreasing grid spacing, showing an almost linear trend for grid spacings of 2 m and less. For the ESG method we observe generally lower values of θ * than for the IL method, which is consistent with the findings from our discussion of Fig. 9 (see above). Mathematically, a smaller θ * results in larger values of L as can be inferred from Eq. 8 when entered in Eq. 2. Figure 10c (see also Table 2) shows the statistics of L for both methods. Indeed, we observe a continuous increase in L for the ESG method. However, the opposite trend is observed for the IL method. As L is defined through both u * and θ * , the trend of decreasing surface shear stress with decreasing grid spacing also comes into play. As u * has a higher weight in the calculation of L (as its number enters squared), the resulting L reflects basically the trend of u * for the IL method. In contrast, u * for the ESG method was found to be rather constant (i.e., converged) so that the trend in θ * then dominates the behaviour of L, leading to different trends of L for IL and ESG methods.
The finding that L increases with decreasing grid spacing for the ESG method while it remains rather constant or slightly decreasing for the IL method would indicate a different effect on the wind-speed and temperature profiles. For the ESG method we might expect stronger and larger eddies with increasing values of L (by increasing wind shear and/or smaller surface heat flux) and hence possibly a greater depth of the boundary layer, while we would suspect that the opposite is the case for the IL method. Interestingly, Fig. 11 reveals that for both methods, the boundary-layer depth decreases to the same amount (and that apparently there is neither a correlation of the boundary-layer depth changes to L nor to u * ). This rather surprising finding leads us to the conclusion that the lack of grid convergence here is not linked to the surface fluxes and stresses. The variance profiles shown in Fig. 12 indicate a monotonically decreasing mixing with decreasing grid spacing throughout the boundary layer, giving further support to the fact that no grid convergence is reached (the same is found for the heat-flux profiles, not shown). Our findings are in line with Sullivan et al. (2016) who also investigated the GABLS1 case with grid spacings down to 0.39 m and showed similar results for mean profiles (thus suggesting a similar decrease in the turbulent mixing). We thus ascribe this to generally under-resolved turbulence in the LES runs. It remains unclear, however, why the variance profiles (or the total TKE, not shown) decrease with decreasing grid spacing, as more and more turbulence should be resolved. It goes beyond the scope of this paper to further investigate this. While Fig. 11 suggests that differences are largest in the upper boundary layer, these still might be the result of inadequate resolution of the turbulence near the surface. This might lead to an erroneous surface friction and result in different boundary-layer depths, which in turn are most prominently visible at the top of the boundary layer. However, we can also suspect that mixing processes near the top of the boundary layer might contribute to the grid dependence. Moreover, numerical dissipation of the advection scheme which is a known feature for large wavenumbers might contribute to the grid dependence. Furthermore, the SGS model formulation might play an important role as suggested by previous studies (e.g., Porte-Agel 2000;Porte-Agel et al. 2004;Basu and Porte-Agel 2006;Lu and Porte-Agel 2013). Previous studies that simulated the GABLS1 case could achieve much better grid convergence when using dynamic SGS models (Beare et al. 2006;Basu and Porte-Agel 2006;Lu and Porte-Agel 2013). Note, however, that the GABLS1 case differs from the present set-up by a capping inversion, which might inhibit the growth of the boundary layer and might lead to under-resolved turbulence in the entrainment zone. The capping inversion thus potentially has an effect on grid convergence. A more detailed analysis of these processes goes beyond the scope of the present study. For a more rigorous discussion on under-resolved turbulence in stable conditions, see Sullivan et al. (2016). In Fig. 10 As Fig. 3, but for the SBL cases. a-c Show data for u * , θ * , and L, respectively summary, however, we must note that the observed grid convergence effects are likely LES model dependent.
Finally, comparing the resulting boundary-layer depths between the IL and ESG method, we observe a tendency for deeper boundary layers with the ESG method (cf. surface-layer height in Fig. 2) by about 6%, which can be ascribed to the higher surface friction, but which was less clearly visible for neutral stratification, possibly because the changes due to different grid spacing dominated over those imposed by the changes in u * . The variance profiles for stable conditions confirm the relationship between boundary-layer depth and u * and reveal consistently increased values for the ESG method (Fig.12).
Besides the behaviour of the mean values of u * and θ * that was discussed above, Fig. 10 interestingly also shows that the variation of these parameters in time (standard deviation and minimum/maximum values) has a strong dependence on the grid spacing with much less variability for smaller grid spacings. This holds for both methods applied and is related to the fact that the first computational grid level is placed much higher in the atmosphere (e.g., 1.4 m and 0.52 m in cases IL_d3_sbl and IL_d1_sbl, respectively). This results in much more fine-scale turbulence structures for finer grid spacings (see Fig. 13) and has direct implications for calculating the horizontal average of u * . As can be seen in Fig. 13, only few large turbulence structures occupy the computational domain so that the sampling error for a horizontal average is larger than for finer grid spacings where more structures are sampled. Also, due to the higher elevation for coarser grids, higher wind speeds are suggestive to broaden the frequency distribution. This finding supports our reasoning above that very fine grid spacings are required to resolve the turbulence under stable stratification reliably.

Summary and Outlook
In the present study we introduced an improved boundary condition based on MOST for use in LES models. The main concept behind this boundary condition is the use of an elevated level for evaluating the MOST relationships in order to calculate the surface shear stress and scalar fluxes, instead of being limited to the first grid level above the surface. This concept is based on the previous work for neutral conditions by Kawai and Larsson (2012) and was here extended to non-neutral conditions. By taking into account three criteria for the choice of this elevated level, violation of the assumptions of MOST as outlined by Basu and Lacser (2017) is consistently avoided. Also, the improved method ensures the correlation between the near-surface turbulent structures and the surface fluxes, and thus does not corrupt the flow features near the surface. For evaluation of the improved boundary condition, we conducted a Fig. 11 Mean profiles of the horizontal (upper panels) wind speed and potential temperature (lower panels) for cases with the IL method (left) and cases with the ESG method (right). The profiles were averaged over 1 h of simulation time comprehensive set of LES for both a neutral and a weakly stable boundary-layer configuration using both conventional boundary conditions based on evaluating MOST at the first grid level and with the proposed new boundary condition at an elevated level.
The results of the performed LES indicate that the improved boundary condition is a solid method to avoid discrepancies between the surface fluxes (shear stress and surface sensible heat) and the resolved surface-layer wind-speed and temperature profiles by the LES model under neutral and stable conditions. For neutral conditions, it could be shown that the surface shear stress is not only corrected towards larger values that are in line with the resolved wind-speed profile, but it also effectively corrects the excessive wind shear observed near the surface. We also found that the correction implied by the improved boundary condition performed better for the surface shear stress than for the surface sensible heat flux in stable conditions. We could ascribe this to the fact that under stable conditions, both u * and θ * must be calculated. These are linked via a single parameter, L. Here, we found that the calculated L represented a compromise to achieve best agreement with both the simulated wind-speed and temperature profiles, but has an inherent bias to better represent the windspeed profile. Moreover, we could show that the simulated boundary layers with the proposed boundary condition are deeper (under stable conditions by 6%), which is an expected result as the correction of the surface shear stress involves higher values of u * and thus increased turbulent mixing.
The starting hypothesis for the present study was that the well-known lack of grid convergence in LES of the SBL might be attributed to the lack of grid convergence of the surface shear stress and the surface sensible heat flux. Our results gave evidence that this is not the case and that there is no link between the differences in the mean profiles and differences in the surface fluxes, particularly for stable conditions. Instead we suppose that LES models with current resolution are still not able to capture strong curvature very close to the surface, even at grid spacings in the order of 1 m, which might explain the observed tendency of shallower boundary layers for finer grid spacings. Previous studies also suggested that the formulation of the SGS model is a key aspect for grid convergence. Furthermore, numerical aspects like numerical dissipation, the reduction of higher-order advection schemes near the surface, or effects at the top of the boundary layer might be possible causes for the observed lack of grid convergence. LES runs using very sophisticated SGS models and a vertical nesting approach with very fine grid spacings near the surface would be a logical continuation to investigate at what grid spacing true grid convergence can be achieved and what the responsible processes are. We are currently working on the implementation of such a nesting approach in the PALM model and plan to tackle this issue in a follow-up study.
Our results clearly demonstrate that the proposed boundary condition eliminates key issues in the application of MOST as a boundary condition in LES models and can thus be considered as an addition or alternative to advanced SGS schemes in order to avoid the excessive nearsurface wind shear and generate more realistic boundary-layer heights in LES models. Also, the MOST-based boundary conditions proposed by Marusic et al. (2001) and Hultmark et al. (2013) can probably be improved based on our findings. However, our proposed method also has some limitations. First, as horizontally-averaged values are used to calculate the surface fluxes, the method is not applicable in the presence of surface heterogeneity, imposed for instance by heterogeneous land use or roughness as in such cases the atmospheric state in an elevated level might no longer be representative for the underlying surface (e.g., due to an internal boundary layer). Furthermore, the method is difficult to apply in complex terrain or in the presence of buildings, so usage might be restricted to academic idealized studies. Second, Eq. 13 in its current form is only valid for non-zero mean wind speeds and needs to be modified to be used under free convective conditions, i.e, by falling back to the IL method in such cases to preserve a minimum friction velocity (see also Schumann 1987). Furthermore, note that we did not vary the parameters z sl and z 0 in our simulated cases. Hultmark et al. (2013) and Stoll and Porte-Agel (2006), however, showed that there might be sensitivities on the results. More rigorous testing of the robustness of the methodology and the dependence of the results on surface roughness are desirable. Due to the computational demands involved, particularly for the finest grid spacings used, we could only do a few test runs where we varied z sl and therein we did not see a strong sensitivity of our general results and findings. Note, however, that a height-dependence of the results is to be expected as the constant-flux layer is a theoretical construct and, in reality (or the model), fluxes do vary with height within the surface layer.
Finally, in order to relax the strict requirements for the choice of the grid spacing in LES models, we plan to incorporate the roughness sublayer in the boundary condition formulated as outlined by Harman and Finnigan (2007).
Two different methods to improve grid convergence are proposed in literature: usage of a dynamic and thus sophisticated SGS model, or the correction of the surface boundary condition proposed in the present study. At the time when this study was conducted, the PALM model did not offer a dynamic SGS model. However, during the revision of the manuscript, based on the reviewers' comments, we decided to repeat some LES runs with a dynamic SGS model after (Heinz 2008;Mokhtarpoor and Heinz 2017), which had been implemented in the PALM model in the meantime. Results from these runs for stable stratification using the IL method as a boundary condition are shown in Figs. 14 and 15 (in comparison with the data using the default Deardorff SGS model together with the IL method). The mean profiles of wind speed and temperature shown in Fig. 14 reveal that no improvement regarding grid convergence could be achieved compared to the Deardorff scheme (cf. Fig. 11) and that the simulated boundary-layer depth decreases for decreasing grid spacings as observed for both IL and ESG methods in conjunction with the Deardorff scheme. Figure 15 shows that the same holds for convergence of u * , θ * , and L. For fine grid spacings, we only note marginal differences between Deardorff and dynamic SGS model results. Significant differences are only visible for coarser grid spacings of 5.2 m and 2.8 m. This finding lets us conclude that the importance of the SGS model weakens as the turbulent flow is more and more resolved, while it might affect the results most in such cases where the grid spacing is too coarse to resolve the turbulent flow sufficiently. These findings are of course based on one flavour of dynamic SGS models and, more sophisticated models might give a different answer. The data at hand, however, give no indication that an improved SGS model potentially solves the issue of grid convergence for the stable boundary layer.