Effects of stability functions in a dynamic model convective boundary layer simulation

Dynamic subgrid models are increasingly being used in simulations of the atmospheric boundary layer. We have implemented several variant forms of dynamic models in the UK Met Office Large Eddy Model (MetLEM), including a state‐of‐the‐art Lagrangian‐Averaged‐Scale‐Dependent (LASD) model. The implementation includes optional use of stability functions in the specification of the eddy viscosity and diffusivity, as well as optional use within the dynamic calculation of the Smagorinsky parameter. This paper reports on the behaviour of the LASD model with different choices for the inclusion and treatment of stability functions in convective boundary layer simulations at different resolutions. Results are compared against a high‐resolution Large‐Eddy simulation (LES) and against simulations employing the Smagorinsky–Lilly subgrid model. We conclude that the use of stability functions improves the behaviour of the LASD model in the grey zone regime. Moreover, a careful treatment of the stability functions in the calculation of the dynamic parameters, while attractive theoretically, is found to be unnecessary in practical terms.


| INTRODUCTION
While a variety of sub-grid parameterizations for largeeddy models (LEMs) now exist, the simple Smagorinsky scheme and its variants are still probably the most widely-used approach (Yang, 2015).This model, first proposed by Smagorinsky (1963), is based on the application of the mixing length formulation in three dimensions (Mason, 1994).The eddy viscosity ν is obtained by assuming that the small scales are in equilibrium, so that energy production and dissipation are in balance.This yields an expression of the form ν c s Δ ð Þ 2 j S j where Δ is the filter width (which is proportional to the grid size), c s is the Smagorinsky constant, and a product of the two makes up the mixing length and S is the modulus of the rate of strain tensor S ij given by j S j = 2S ij S ij À Á 1=2 .
Abbreviations: LASD, Lagrangian-averaged scale-dependent dynamic model; NO-STABF, LASD model without stability functions; PAR-STABF, LASD model with stability functions partially implemented; SMAG, Smagorinsky model; STABF, LASD model with stability functions fully implemented.Mason (1994) showed that reducing the value of c s from about 0.8 results in finer scales being resolved, with c s ≈ 0.2 providing the most accurate simulation in the inertial-subrange for a neutral boundary layer; reducing the value of c s further results in larger finite-differencing errors.Mason and Callen (1986) used a centred difference scheme (Piacsek and Williams, 1970) in their model and found c s values of 0.1 and 0.07 to produce poor numerical representations of the flow, nonphysical structure, with significant features dominated by single grid-point values.The optimal value of c s is dependent on the flow being simulated (Yang, 2015) and a single, universal constant does not exist.Germano et al. (1991) proposed a dynamic Smagorinsky model that calculates a suitable, local, value of c s using test filtering at the smallest resolved scales in the flow.Such locally-derived values of c s have been found to be highly variable both in magnitude and sign (You and Moin, 2009).To reduce the variability Germano et al. (1991) took a horizontal plane average of the calculated values of c s .Plane-averaging is unsuitable for inhomogeneous flows over complex terrain.Meneveau et al. (1996) developed an alternative Lagrangian dynamic model that averages in time following fluid pathlines.
Both the Germano et al. (1991) and the Meneveau et al. (1996) models use 2Δ as a test filter scale and assume scale-similarity, so that c s,Δ = c s,2Δ .The scaleinvariance assumption is reasonable in the inertial subrange but is not expected to hold when Δ falls near a transition scale that separates different physical processes occurring in distinct ranges of scales (Honnert et al., 2011).Porte-Agel et al. (2000) proposed a scale-dependent version of the plane-averaged dynamic model in which a second test filter is applied to determine how the coefficient changes across scales.The plane-averaging implies that although the Porte-Agel et al. (2000) model is scaledependent, it cannot be applied to nonhomogeneous flow over complex topography.Bou-Zeid et al. (2005) combined the Lagrangian and scale-dependent methods to develop a model applicable for inhomogeneous flows over complex topography.
In various applications, the Smagorinsky model has also been extended to include stability functions in order to account for the effects of stratification (e.g., Beare and Macvean, 2004;Efstathiou and Beare, 2015).This is because when the ideal of the long inertial subrange is not achieved, buoyancy effects are expected to influence the subgrid model (Brown et al., 1994;Mason, 1994).Applications of dynamic models so far have typically excluded these stability functions (e.g., Kleissl et al., 2006;Huang and Bou-Zeid, 2013).An exception is Kirkpatrick et al. (2006) though the emphasis of that study was not on the effect of the stability functions per se.
This study is focused on the use of the scaledependent dynamic model in convective boundary-layer simulations at different resolutions from the near-LES through to the so-called grey zone regimes.The study shows the effects of stability functions on the CBL simulations with changing horizontal resolutions.The stability functions are employed fully in the calculation of c s , and the eddy viscosity and diffusivity, as well as partially by being used only in the calculation of the eddy viscosity and diffusivity.

| NUMERICAL MODEL, DYNAMIC MODEL IMPLEMENTATION, AND SIMULATIONS
The simulations in this study were performed using the UK Met Office Large Eddy Model (MetLEM) described in Shutts and Gray (1994) and Brown et al. (1994), and briefly summarized in Section 2.1.The turbulent stress and subfilter-scale heat flux terms are parameterized using the Smagorinsky-Lilly scheme as well as different variations of the dynamic Smagorinsky scheme which were coded into the MetLEM specifically for the present study.Details of the schemes reported upon are given in Sections 2.2 and 2.3.
Here u is the flow velocity, θ is the potential temperature, p is the pressure, ρ is the density, B is the buoyancy, τ is the subgrid stress, h θ is the subgrid scalar flux of θ, δ ij is the Kroneker delta function, ϵ ijk is the alternating pseudo-tensor, and Ω is the Earth's angular velocity.The subscript s indicates a height-dependent reference state.The term on the left hand side of Equation ( 1) is the total time-derivative of momentum.The first term on the right is the pressure gradient; the second term, buoyancy, is nonzero only in the vertical; the third term is the divergence of the turbulent stress; the final term is the Coriolis acceleration.The left hand side of Equation (3) is the total derivative of potential temperature, while the first term on the right hand side is the divergence of the heat flux.The second and the final terms represent the effect of latent heating or cooling and the effect of radiation, respectively.These last two terms are zero in our study.The second-order centred difference scheme of Piacsek and Williams (1970) is used for advecting both momentum and scalars.The p 0 needed for Equation (1) is calculated using the Poisson-like elliptic equation given by. where which is solved by performing a Fourier transform in the horizontal.

| The Smagorinsky-Lilly Sub-Filter model
The sub-filter stress, τ ij , and heat flux, h θ i in Equations ( 1) and (3) as parameterized in the Smagorinsky-Lilly scheme (Smagorinsky, 1963;Lilly, 1965) are specified respectively by and where ν is the sub-filter eddy viscosity and ν h the corresponding diffusivity for scalars.The rate of strain tensor is defined by and the eddy viscosity and diffusivity are prescribed through where Ri p is the local Richardson number given by Ri p = ∂B=∂z S 2 : ð11Þ B = g θ s θ 0 is the buoyancy, and f m and f h are Richardson-number dependent functions applied to the viscosity and diffusivity, respectively as described by Shutts and Gray (1994).The stability functions are defined based on the Richardson number as described below.
For Ri p < 0 Default values of the constants shown in the above equations are: a = 1/Pr N , b = 40.0,c = 16.0,g = 1.2, h = 0.0, Pr N = 0.7, and Ri c = 0.25.The stability functions encourage mixing for smaller values of Ri p , and suppress mixing for larger values.Ri p < 0 is characterized by buoyant instability, 0 < Ri p < Ri c is Kelvin-Helmholtz instability, while Ri p > Ri c is associated with stable conditions.S is the modulus of S ij given by λ is a mixing length scale given by where κ is von Karman's constant with value 0.4, λ 0 = c s Δ, Δ = (Δx + Δy)/2 is the horizontal grid spacing, and c s is a constant diagnostic of the ratio of the subgrid filter scale to the grid scale known as the Smagorinsky constant, chosen to be 0.23 in this study.

| Implementation of the Lagrangian-Averaged Scale-Dependent (LASD) dynamic model in the MetLEM
We have implemented several variants of dynamic Smagorinsky models in the MetLEM.In this study, we focus on the Lagrangian Averaged Scale Dependent (LASD) model of Bou-Zeid et al. (2005) to dynamically compute a suitable value of c s from the smallest resolved scales.The Germano identity (Germano et al., 1991) relates the subgrid stresses at different scales as where the overbar and tilde indicate filtering on scales Δ and αΔ, respectively, and L ij represents the stress on scales intermediate between Δ and αΔ.
Using a Smagorinsky representation, the stresses are parameterized by: and the parameterization error, which is to be minimized by the dynamic procedure, is therefore Substituting the Smagorinsky formulation produces where . In the Germano model, a scaleinvariance assumption is made for c s so that β = 1.An optimum value for c s,Δ can then be obtained by minimizing the square error e ij e ij .
In practice, pointwise, instantaneous values of c s,Δ soderived can be noisy (Kleissl et al., 2005) and lead to numerical difficulties so it is common to introduce some averaging procedure (Meneveau and Katz, 2000).In the Lagrangian procedure, the coefficient is obtained by minimizing the weighted time-average of the local error contraction e ij e ij over pathlines.The average error metric is defined by where z(t 0 ) are the previous positions of the fluid elements and W(t) is a relaxation function that allocates larger weights to the more recent history of the coefficients (i.e., W(τ) is a decreasing function of τ).c s,Δ is then determined from the stationary point of E, resulting in where and Several weighting options could be used, but an attractive choice is to take an exponential weighting of the form W t−t 0 ð Þ= T − 1 e − t − t 0 ð Þ=T which allows the integral definitions to be replaced by transport equations: The generalization to a scale-dependent dynamic model can be motivated by consideration of model resolutions approaching the grey zone for which the simulation may not lie fully within the inertial subrange (e.g., Efstathiou et al., 2018).The idea is to extend the above model by dynamically computing β alongside c, and to do so by means of a second test-filter operation (Porte-Agel et al., 2000).The additional equations we have implemented are closely analogous to those presented above and can be found summarized in Efstathiou et al. (2018) or in a more detailed presentation and discussion by Bou-Zeid et al. (2005).It may be noted that the Smagorinsky representation used by Bou-Zeid et al. (2005) is given by which implies that f m = f h = 1 as opposed to Equation (20) where f m and f h are defined in Equations ( 12)-( 17).Hence, the c s,Δ obtained by Bou-Zeid et al. (2005) is not influenced by stability functions, while the one described here includes the stability functions.
Once the value of c s,Δ has been computed through the above procedure, the mixing length λ is determined as λ = c s,Δ Δ without the need for any special treatment close to the surface as in Equation ( 19) for the Smagorinsky-Lilly model.Finally, the mixing length so-derived is used within Equations ( 6) and ( 7) for the eddy viscosity and diffusivity.

| Simulations
Simulations were performed for a dry convective boundary layer with a strong inversion and a strong surface sensible heat flux, following Sullivan and Patton (2011).These simulations imposed a surface sensible heat flux of 241 Wm −2 and geostrophic winds of (U g , V g ) = (1, 0) ms −1 and (U g , V g ) = (10, 0) ms −1 were considered.Similar conclusions are drawn from these weakly and strongly sheared cases, and so only results from the 10 ms −1 simulations are presented here.The initial inversion height is z i = 1024 m and the surface roughness used is z 0 = 0.1 m.The initial sounding has a three layer structure with a strong inversion as follows: 300 K for 0 < z < 974 m, 300 K + (z − 974 m) 0.08 Km −1 for 974 < z < 1, 074 m, 308 K + (z − 1, 074 m) 0.003 Km −1 for z > 1, 074 m.This profile results in negative Richardson numbers up to an altitude of 270 m, less than 0.25 up to about 880 m and larger numbers beyond.
Horizontal grid lengths of Δ = 25, 50, 100, 200, and 400 m were used with a horizontal domain size of (9.6 km) 2 .These grid lengths correspond to 0.024z i , 0.049z i , 0.098z i , 0.195z i , and 0.39z i .The vertical grid length is constant for each simulation and was set to 0.4Δ for Δ = 25, 50, 100 m, and then kept at 40 m for Δ = 200 and 400 m.The domain height is 2 km and the horizontal boundaries are cyclic.The simulation time is 4 hours and the data is averaged from 12,600 to 14,400 s, which is the last 30 min of the simulation.This last 30 min represents a time when all simulations have long since reached equilibrium.
Separate runs were performed using the Smagorinksy-Lilly subgrid model (SMAG; Section 2.2) and with the LASD dynamic model (Section 2.3) using stability functions (STABF), without stability functions (NO-STABF) and with a partial inclusion (PAR-STABF) of stability functions.The NO-STABF formulation without stability functions corresponds to setting f m = f h = 1 throughout all of the calculations within the sub-grid model.In the PAR-STABF simulations, the Smagorinsky coefficient, and hence the mixing length, is calculated without the stability functions, but these functions are retained in Equations ( 6) and ( 7) for the final calculation of the viscosity and the heat diffusivity.The highest resolution run, an LES at 25 m grid spacing, was performed with SMAG only and will be used as a reference against which to compare the other runs.

| RESULTS
Results from the 100 m grid spacing simulation are qualitatively similar to those obtained with 50 m spacing and so the discussions below focus on the relatively wellresolved 50 m case, the 200 m case towards the grey zone and the 400 m case within the grey zone.

| Temperature variance and mean temperature
The temperature variance profiles are shown in Figure 1 for 50, 200, and 400 m grid spacing, together with the associated LES simulation at 25 m coarse grained to 50, 200, and 400 m.Profiles for the different variants of the dynamic model are presented, as well as the original Smagorinsky.The peak in the upper boundary layer is associated with the occurrence of the inversion layer.For a 50 m grid spacing, the height of the peak is the same for all subgrid models, and occurs at a slightly higher altitude than for the LES.The results agree with previous studies (e.g., Sullivan and Patton, 2011) that showed that the boundary layer height somewhat increases as the resolution is reduced.There are relatively slight differences in the magnitudes of the peaks.Larger differences appear at coarser resolutions.When the grid spacing approaches the grey zone (i.e., at 200 m), the peak variance of the NO-STABF simulation occurs at a markedly greater height.The STABF and PAR-STABF versions of the LASD simulations produce a peak at a lower height with a relatively constant altitude for different horizontal resolutions.The peak for the NO-STABF run gets smaller with reduced resolution.These results are consistent with the temperature profiles of Figure 2, in which the inversion is seen to become smoother using the dynamic models.The peak of the variance with SMAG occurs at a lower height than the LES at a 400 m grid spacing.As the resolution is reduced, there is an increasing spread in the locations of the peaks for the different sub-grid models; the NO-STABF run also gives an increasingly wider peak.Larger differences in the variance are also observed close to the surface, where the NO-STABF simulation is associated with a larger variance.The SMAG run is associated with the least amount of variance close to the surface.
When a grid spacing of 50 m is used, the domainaveraged temperature profiles in Figure 2a are almost identical across all models, with slight variations in the inversion layer and close to the surface.For all the 50 m simulations, the inversion appears at a higher level than in the LES, consistent with earlier studies (e.g., Sullivan and Patton, 2011;Beare, 2014).As the resolution is made coarse, the temperature profiles from different subgrid averaged potential-temperature variance profiles obtained at different resolutions, using the Smagorinksy-Lilly (SMAG; blue dotted) and LASD subgrid models with (STABF; red), without (NO-STABF; green) and partially-implemented (PAR-STABF; cyan dashed) stability functions.Also shown is the corresponding profile obtained from the LES reference run (black dash-dotted), coarsegrained in both the vertical and the horizontal to the relevant grid length models start to show larger differences.For a grid spacing of 200 m, the Smagorinsky simulation is associated with a steeper inversion compared to all the versions of the LASD model.The NO-STABF version has the least steep inversion and remains cooler than all the other simulations between around 1,100 and 1,400 m.This simulation is also associated with higher temperatures within the boundary layer itself.The STABF and PAR-STABF versions of LASD are very similar to one another and lie between the NO-STABF and SMAG simulations.
As the resolution is made even coarser, for a grid spacing of 400 m in Figure 2c, the same trends continue.The inversion layer of the SMAG run moves to a lower height, to the extent that it becomes lower than the one associated with the LES.The inversion associated with the NO-STABF run is pushed further up, and is yet for the mean potential temperature profiles smoother such that it remains cooler than other simulations up to a height of 1,500 m.The STABF and PAR-STABF simulations remain similar.The boundary layer gets warmer with lower resolution, especially without the inclusion of stability functions.Table 1 shows the domainaveraged potential-temperature difference between different versions of the dynamic model with 400 m grid spacing and the corresponding coarse-grained LES for heights of 20, 500, 900, and 1,220 m.The selected heights represent the surface layer, the middle of the boundary layer, the entrainment region and the inversion layer.The difference between the NO-STABF version and the LES is much larger than with the STABF and PAR-STABF versions, confirming results shown by line plots.

| Turbulent kinetic energy
In comparison with the Smagorinsky model, the LASD model is associated with larger values of the resolved as well as subgrid TKE close to the surface (Figure 3).At 400 m grid spacing, the SMAG run is associated with the least amount of resolved TKE throughout most of the boundary layer.The NO-STABF run is associated with the most resolved TKE close to the surface and the least in the vicinity of the boundary layer inversion as compared to other versions of LASD.The LES and SMAG runs, both of which use constant values of c s , are associated with the least amount of resolved TKE in those regions.
The LES is associated with the least amount of subgrid TKE as expected because of a smaller reliance on subgrid models at higher resolution (Figure 3).Close to the surface, the LASD runs are associated with more subgrid TKE than the SMAG run.The shapes of the TKE profiles close to the surface in the LASD model are more similar than SMAG to the shapes at higher resolution.While all the others increase from the surface to a peak slightly above the surface and start to reduce immediately, the SMAG profile does not exhibit a distinct peak at 200 and 400 m grid spacings (Figure 3b,c).A near-surface spike in the TKE was found in other dynamic-model studies (Porte-Agel et al., 2000;Efstathiou et al., 2018) due to poorly resolved turbulence, and Efstathiou et al. (2018) who used the same model as in this study found this spike to be a result of horizontal velocity fluxes.From just over 100 m to the upper boundary layer, the SMAG run is associated with the largest values of the subgrid TKE.The NO-STABF run is associated with the least amount of subgrid TKE compared to other runs throughout the boundary layer.Efstathiou and Beare (2015) showed that the energy close to the boundary layer top is larger in the grey zone especially when low resolution is used in the vertical.

| Viscosity and Smagorinsky coefficient
The shape of the viscosity profile is similar to the shape of the subgrid TKE (Figure 4).The viscosity increases as the grid spacing is made larger, which indicates a greater reliance on the subgrid models as the resolution is lowered.The shape of the viscosity profile with the LASD model is similar across resolutions, in contrast to SMAG which has a peak at increasing height for increasing grid spacing.Larger values of viscosity are associated with larger values of subgrid TKE.The largest values of viscosity in the LASD runs occur close to the surface, while in the SMAG runs they occur above a height of 100 m with the 200 m grid spacing and above 200 m with 400 m grid spacing.The viscosity in the SMAG simulations is the largest through most of the boundary layer.The NO-STABF simulation is associated with the smallest amount of viscosity in the lower half of the boundary layer.
The differences in viscosity for the LASD simulations with and without stability functions cannot be simply explained as a result of the mean value of the calculated Smagorinsky coefficient.However, for the LASD model as used here, it is important to recall that the calculated Smagorinsky coefficient is applied locally: that is, different values are applied for different grid points in the horizontal as well as in the vertical.To assess whether there are differences in the spread of the derived values, cumulative probability density functions (PDFs) were computed.The cumulative PDFs close to the surface (at 40 m) are presented in Figure 5 and show that the distributions are rather similar amongst the different versions of the LASD model.The same is true for PDFs at other levels investigated.
It may be noted that the values of the square of the Smagorinsky coefficient, c 2 s , are not allowed to become  clipped; this number increases to about 25% with a grid spacing of 200 m and about 30% with a grid spacing of 400 m.Thus, as the resolution is coarsened the amount of clipping needed increases.
Clipping is also applied to excessively large values, so that c 2 s is not allowed to be greater than 0.6.The simulations did not produce large values of c 2 s throughout much of the boundary layer.Most of the clipping is needed far above the boundary layer and this is associated with very large values of the Richardson number.The STABF simulation is associated with less clipping of the smaller numbers and more clipping of the larger numbers because of the influence of the Richardson number in the calculations of c s .The viscosity profile in Figure 4 however shows that this does not have an impact on the simulations produced.functions improve the dynamic model simulations, making their predictions closer to the LES results.The stability functions also increase the horizontal winds in the vicinity of the inversion layer, while reducing the vertical winds (not shown).
The dynamic model reacts naturally to the surface boundary, avoiding any need to prescribe a reduction of the mixing length scale close to the surface.Additional advantages are that profiles of viscosity and TKE in the LASD model are more similar in shape to those of the high-resolution simulations.As the resolution is reduced into the grey zone, the level of the inversion associated with the Smagorinsky model decreases and so do the resolved velocities.Using a constant value of c s dampens the resolved flow and reduces the resolved-scale transports of momentum and heat.The dynamic model maintains values of TKE and transports that are more comparable to the high-resolution simulations, with the stability functions serving to prevent excessive turbulent motions, particularly close to the surface.The height of the maximum variance remains close to the LES value, whereas it is either too low with Smagorinsky, or too high without stability functions.
Using the LASD model in the grey zone (e.g., at 400 m grid spacing) results in larger values of TKE close to the surface and has more mixing than the original Smagorinsky model.However, this increased vertical turbulence may be excessive, causing a marked reduction in the steepness of the inversion, with higher temperatures in the boundary layer, and a larger negative heat flux in the entrainment zone (not shown).The use of stability functions as part of the calculation of the diffusion and viscosity improves the dynamic model simulations.For example, the near-surface TKE is reduced somewhat, and the height of the peak in the potential-temperature variance increases only slightly compared to the LES.Moreover, the boundary-layer warming is reduced.
Accounting for these functions in the derivation of a dynamic Smagorinsky coefficient as was done for STABF has very little impact.Their inclusion in Equation (24), for example, is certainly more satisfactory from a theoretical perspective in providing a self-consistent dynamical treatment of the subgrid stresses based on Equations ( 6) and ( 9), but it does not seem to be necessary.The introduction of stability functions only in the calculation of viscosity and diffusion as done for PAR-STABF provides results that are as good as including the functions in calculation of c s .
Overall, these results show that the use of stability functions in a quasi-equilibrium convective boundary layer simulation results improves in the performance of the LASD model in the grey zone regime.The work presented in this paper represents an idealized daytime scenario which is characterized by a quasi-steady CBL.It may be noted that the same model including stability functions has also been tested for a morning transition case and showed improvements over the standard Smagorinsky model (Efstathiou et al., 2018).Turbulent mixing in the presence of clouds and latent heating remains an open subject in atmospheric convection, with more observations needed to compare directly with model results and establish the performance of sub-grid parametrizations (Feist et al., 2018).
As in Figure 1, for the profiles of (a-c) resolved and (d-f) subgrid contributions to the TKE.The LES profiles on their native 25 m grid are shown for reference in (a) only As in Figure1, for viscosity.The LES profiles on their native 25 m grid are shown for reference in (a) only CONCLUSIONS Introducing stability functions in the calculation of the viscosity and diffusivity has a substantial impact on the behaviour of the LASD model at near grey-zone resolutions.The inversion layer remains steeper than for the case without stability functions, the temperature bias is reduced in the boundary layer, and excessive near-surface resolved TKE is also ameliorated.Higher temperatures are associated with a deeper boundary layer, and so the stability functions also help to better control the boundary-layer height, and likewise the associated peak in the temperature variance.Overall, the stability As in Figure1, but here for the LASD runs only, showing the cumulative PDFs of the Smagorinsky coefficient at the height z = 40 m