Modelling of Stably Stratified Atmospheric Boundary Layers with Varying Stratifications

A recently developed explicit algebraic Reynolds-stress (EARS) model is validated for an idealized representation of the night-time high-latitude stably stratified atmospheric boundary layer. The simulations are made with four surface cooling rates that result in weakly to moderately stratified stable boundary layers. The predictions of the EARS model are compared to high-resolution large-eddy simulations (LES) of Sullivan et al. (J Atmos Sci 73(4):1815–1840, 2016). First- and second-order statistics are shown to be well predicted by the EARS model. The EARS model also predicts the horizontal turbulent fluxes and turbulence anisotropy and these compare well with the LES results. The sensitivity to the model coefficients is studied by comparing the EARS model results with LES results. Finally, we propose a new scaling for the production of turbulence kinetic energy and show that the EARS model captures the essential trends of the LES results for different cooling rates.


Introduction
A stable atmospheric boundary layer (SABL) generally occurs with longwave radiative cooling of the land surface when the underlying surface becomes colder than the air above it, or in the case of warm air advection over a colder surface. These processes are followed by the transfer of energy from the atmosphere to the surface and lead to stably stratified inversion layers.
An important challenge for SABL research is understanding the relationships between the mean wind and temperature fields, and the average turbulent fluxes and variances. The increase of computational power has made it possible to use techniques such as large-eddy simulation (LES) to investigate the structure and properties of the SABL, which can help improve turbulence models (Huang et al. 2013). A comparison of LES results is presented in Beare et al. (2006) for an idealized set-up of a SABL in dry conditions with prescribed B Velibor Želi velibor@mech.kth.se 1 Department of Mechanics, Linné FLOW Centre, KTH Royal Institute of Technology, 10044 Stockholm, Sweden temperature forcing on a flat homogeneous surface, the so-called first Global Energy and Watercycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS1). The LES results have been utilized in the intercomparison of turbulent closure models typically used both in operational weather forecasting and research (see Cuxart et al. 2006). The study shows large variations among the model results and indicates that the turbulence models have a tendency to overpredict turbulence mixing in the SABL. As pointed out by Holtslag et al. (2013), turbulence models have a significant effect on the overall accuracy of climate and weather numerical models.
Recently, a high-resolution LES study by Sullivan et al. (2016) extended the set-up of the GABLS1 case by modifying the rate at which the surface temperature is changing in time. Varying the cooling rates at the surface produces a set of different stratification levels in the SABL. Turbulence transport weakens with increased stratification and the height of the SABL decreases. These recent LES data make it possible to test different turbulence models to varying levels of stratification. A particularly relevant question is how stratification affects the production of turbulence kinetic energy since many of the turbulence models do not contain the relevant physical processes for capturing the interaction between the turbulence and density stratification. Simpler turbulence models rely on the eddy-viscosity hypothesis (also known as K-theory) in which turbulence kinetic energy is typically modelled using empirical stability correction functions that depend on the gradient Richardson number.
In the present study we investigate the ability of the recently proposed explicit algebraic Reynolds-stress (EARS) model to predict the effects of different levels of stratification in the SABL. A notable difference between the EARS model and simpler models that rely on the eddy-viscosity hypothesis is that, instead of directly modelling turbulent fluxes, these are found by solving algebraic approximations of the transport equations for turbulent momentum and heat fluxes, also known as the Reynolds-stress and scalar-flux equations. The model has been developed by Lazeroms et al. (2013), and since then, a new boundary-condition treatment has been implemented along with a simplified model formulation by Želi et al. (2019). Previously, the result of the EARS model was compared to LES data for the case of the original GABLS1 set-up (Lazeroms et al. 2016). However, simulations with varying cooling rates allow us to study how the EARS model responds to varying stratification levels by comparing the results against the new set of high-resolution LES data from Sullivan et al. (2016).
The following section gives a model description: Sect. 3 explains the set-up of the numerical simulations. In Sect. 4, results of the validation are presented for the EARS model in the model's standard form. Discussion of possible model modification and improvements are given in Sect. 5, and a summary is given in Sect. 6.
where D/Dt ≡ ∂/∂t + U m ∂/∂ m is the material derivative in the direction of the mean flow, and the remaining terms represent (from left to right) diffusion, shear production/destruction, pressure redistribution, dissipation, and buoyancy production/destruction. The large number of strongly coupled prognostic equations makes this class of models hard to evaluate in practice.
The system of equations in (1), (2) will be simplified as algebraic relations further below. Additional prognostic equations for scalar properties are then introduced, namely, the model for turbulent kinetic energy (TKE), K , and half the temperature variance, K θ , where P and G are the production or destruction of TKE due to mean-shear and buoyancy effects, respectively, and is the dissipation rate of TKE; P θ and θ are the production or destruction and the dissipation rate of half the temperature variance; D K and D K θ are the diffusion terms. It is worth mentioning that (3) is half of the trace of (1), meaning that it does not provide any additional information for the case when the transport equation (1) is solved directly, but will do so here since (1) is solved in an algebraic approximation. Compared to other atmospheric turbulence models, the EARS model does not explicitly model the turbulent length scale (mixing length). Instead, it adopts the approach of K − models where an additional prognostic equation for is added as where C 1 , C 2 , and C 3 are model coefficients discussed in Sect. 2.2; K / is the turbulent time scale and K 3/2 / is the turbulent length scale. The equation for the dissipation rate of TKE is to a lesser extent related to the exact equation compared to (3) but it has the same general form with production, destruction, and diffusion effects. The diffusion effects for the scalar properties are modelled as gradient diffusion and read where σ K = 1, σ K θ = 1, σ = 1.3, and ν T is the eddy viscosity obtained from the EARS model solution. Refer to Appendix 1 in Želi et al. (2019) for the full derivation and details regarding ν T in the EARS model. To simplify the system of equations in (1), (2) we use the approach of Rodi (1972Rodi ( , 1976) and introduce the weak-equilibrium assumption that states assuming similarity between the different components. The weak-equilibrium formulation is more general than simply neglecting the advection and diffusion of turbulent momentum and heat transport, and can be expressed with the Reynolds-stress anisotropy Equation 9 is equivalent to the advection and diffusion of a i j being zero, Correspondingly, the weak equilibrium for the heat flux is to neglect advection and diffusion of the normalized flux, and can be written as where Formulation of the weak-equilibrium assumption for heat flux is adopted from Wikström et al. (2000). In Rodi (1972Rodi ( , 1976) the heat flux is not considered. Applying the weak-equilibrium assumption in (1) and (2) results in These two algebraic relations, with modelling of i j , i j , θi , and θi , replace the eddyviscosity and eddy-diffusivity hypotheses.
The following source terms in (3), (4), (5), (14), (15) are given in closed form. This implies that these terms are exactly represented without further modelling where T 0 is the reference temperature, g i is the acceleration due to gravity, and S i j and i j are the mean strain and rotation rate tensors defined as For the remaining terms in (14), (15) we apply the following models, where c 1 , c 2 , c 3 , c θ 1 , c * θ 1 , c θ 2 , and c θ g are model coefficients discussed in Sect. 2.2 and r = 0.55 is the ratio of turbulence time scales for the temperature (τ θ = K θ / θ ) and the velocity (τ = K / ). The pressure terms are modelled according to Launder et al. (1975) and also adopted by Wallin and Johansson (2000), Wikström et al. (2000). Aside from the commonly found terms in (26) there is an additional non-linear term related to c * θ 1 used to avoid model singularity, see Wikström et al. (2000), where also the modelling of θ is discussed in more detail.
Equations 14, 15 are non-linear in the sense that the source terms depend on turbulent momentum and heat flux. Adopting the approach from Lazeroms et al. (2013) and expressed by Želi et al. (2019), it is possible to derive explicit expressions for solutions of the turbulent momentum and heat flux that have the following forms, where f i j and f θi are the functions for the turbulent momentum and heat flux, respectively. For the details of these functions, see Appendix 1 in Želi et al. (2019).
In the context of the Mellor and Yamada (1974) hierarchy of turbulence models, the EARS model belongs to the class of the Mellor and Yamada level-3 closure, see Lazeroms et al. (2016). Computational efficiency of the EARS model is higher than the Mellor and Yamada level-4 models and is comparable to that of the commonly used turbulence models, e.g., the Mellor and Yamada level-2.5 class of models.

Coefficients of the Explicit Algebraic Reynolds-stress Model
In (5) C 1 = 1.44, which is a commonly used value for this coefficient and is also found in Lazeroms et al. (2016). However, instead of C 2 = 1.92 as in Lazeroms et al. (2016), we set C 2 = 1.82. For the case of neutrally-stratified homogeneous shear flow (i.e. when G = 0), this parameter recalibration results in a change of the dimensionless production from P/ = 2.05 to P/ = 1.85. This is close to the experimental results in the non-stratified homogeneous shear flow reported by Tavoularis and Corrsin (1981) and gives a value of κ = 0.38 for the von Kármán constant in a neutrally-stratified atmospheric boundary layer.
We require that the EARS model with stratified effects returns to the baseline modelling in the limit of passive scalar flows where the turbulent velocity and density fields are decoupled. Therefore, in (25) and (26) the coefficients c 1 , c 2 , c θ 1 , c * θ 1 , and c θ 2 have the values from Wallin and Johansson (2000) and Wikström et al. (2000). Coefficients associated with the passive scalar are calibrated in Wikström et al. (2000) using DNS data. When calibrating the coefficients in the EARS model for the case of SABL it is important to retain the general behaviour of the EARS model in the limiting case of a passive scalar. Therefore, among the coefficients in (25), (26) only c 3 , c θ g , and C 3 are directly related to buoyancy effects in the EARS model and can be recalibrated. In stratified flows, buoyancy is believed to have a lesser influence on the dissipative scales and we set C 3 = 0. The remaining coefficients c 3 and c θ g are then calibrated by considering the critical gradient Richardson number Ri c = 0.2, and were changed from c 3 = 0.3 and c θ g = 0.75 to c 3 = 0.35 and c θ g = 1.0. Increasing c θ g decreased the values of K and K θ improving the overall agreement with LES. The current value of c θ g is still far from c θ g = 2 which corresponds to the two-component limit (compare Eqs. 21 and 26). The value of Ri c = 0.2 is chosen as a compromise between observed extreme conditions and the overall agreement with the present LES. However, it should be noted that the critical gradient Richardson number is an artificial asymptotic state in stratified homogeneous shear flow with a subtle balance between production and dissipation never exactly obtained in reality.

Explicit Algebraic Reynolds-stress Model in the Framework of Single-column Model
The single-column model approach is commonly used to model a vertical column of the atmospheric boundary layer when statistical variations in the horizontal directions can be neglected. The previously introduced formulation of the EARS model is applied in the framework of the single-column model by setting x 1 = x, x 2 = y, and x 3 = z. The first two coordinates are in the horizontal plane and the third is in the vertical direction. The main quantities in the single-column model are the mean horizontal wind velocity V = (U , V ) and the mean potential temperature . The main equations under the assumption of horizontal homogeneity and using the Boussinesq approximation are where V g = (U g , V g ) is the geostrophic velocity, f = 2 sin φ is the Coriolis parameter depending on the rotation rate of the Earth and φ is the latitude. The quantities uw and vw are the vertical components of the turbulent momentum flux and wθ is the vertical component of the turbulent heat flux.

Case Description
The modified GABLS1 set-up used herein is taken from Sullivan et al. (2016). The SABL is primarily driven by a steady geostrophic velocity (U g , V g ) = (8, 0) m s −1 and a time-varying surface temperature. The surface temperature is initially set to 265 K and thereafter decreases with a constant cooling rate C r = (0.25, 0.375, 0.5, 1) K h −1 . Increasing the cooling rate of the surface enhances the overall level of stratification in the SABL as compared to the original GABLS1 set-up  in which only C r = 0.25 K h −1 was studied.
The height of the vertical domain is H = 400 m with 231 grid points log-linearly distributed. The height of the first grid cell = 0.9 m. Other parameters are g = 9.81 m s −2 , f = 1.39 × 10 −4 s −1 , reference temperature T 0 = 265 K, and roughness length for momentum and potential temperature z 0 is 0.1 m.
The initial conditions are set in accordance to the LES in Sullivan et al. (2016). After an initial transient, the flow slowly develops according to the specified surface cooling rate. The initial potential temperature is constant at 265 K from the surface up to the initial inversion layer at 100 m and increases above it at a rate of 0.01 K m −1 . The initial wind profile is equal to the geostrophic velocity throughout the SABL. At the time of initialization, turbulence is described by the initial profiles of K , and K θ . The initial TKE profile is K = 0.4(1−z/250) 3 m 2 s −2 for z < 250 m. Above z = 250 m it is equal to 10 −5 m 2 s −2 . The initial profile is not specified in Sullivan et al. (2016). In the present study the initial profile is prescribed through specifying the turbulence time scale τ 0 , which is a diagnostic variable computed from K and . The turbulence time scale τ 0 = K / = 1 s at z 0 and increases linearly to τ 0 = 550 s at 250 m, above which it is constant. This implies that the turbulence eddies close to the surface are small and the size of the turbulence structures increases with height. The initial profile for K θ is set to 10 −4 K 2 , resulting in small levels of potential-temperature variance in accordance with a stably stratified atmospheric boundary layer. The effect of the initial conditions for K , , and K θ on the model results is negligible compared to the effect of the initial condition for .
In contrast to Sullivan et al. (2016), the boundary-condition treatment in the EARS model, with sufficiently small first cell height, does not use Monin-Obukhov functions for calculating friction velocity u * and the characteristic temperature scale θ * . The boundary-condition treatment is taken from Želi et al. (2019) where it was shown that the EARS model with more general boundary conditions gives the same Monin-Obukhov functions that are used in Sullivan et al. (2016) without a priori imposing such behaviour. For details on the boundarycondition treatment in the simulations, see the Appendix.

Results
In order to reach a quasi-steady state solution, the simulations are integrated with the EARS model for the duration of nine physical hours and a timestep of t = 60 s. The atmospheric quantities have been averaged over the last hour of the simulation, as in Sullivan et al. (2016). We explore the possibility of using two different sets of coefficients in the EARS model, and this section discusses the results of the EARS model with the coefficients defined as in Sect. 2. This choice of coefficients is referred to as the EARS model with the standard coefficient values. Section 5 discusses the results of what is referred to as the EARS model with the modified set of coefficients for Eq. 5. Turbulent fluxes in the LES are the total turbulent fluxes, i.e. the sum of the resolved and subgrid-scale contributions. Figures 1 and 2 show mean profiles of the horizontal wind speed V(z), where V = √ U 2 + V 2 , and potential temperature (z). The profiles for the EARS model show similar behaviour to those obtained with LES. Both the model and LES predict the existence of a low-level jet and the decrease of the boundary-layer height with increased surface cooling rates. The EARS model predicts a somewhat higher SABL and a higher low-level jet than in the LES model. The model predicts profile bending in the higher part of the SABL at higher cooling rates, i.e. the sharp change of the profile slope close to the inversion layer for the wind-speed and potential-temperature profiles (separating the inversion layer from the interior of the SABL).
A large number of grid points for the simulations of the SABL is used in order to obtain a smooth representation of the vertical profiles of atmospheric variables. However, a high resolution is not a requirement of the EARS model and the present second-order numerical method. Computationally efficient simulations were carried out with 29 log-linearly distributed grid points between z 0 and H , where the height of the first grid cell is = 1.5 m. Figure 1 shows that the solutions for the finer and the coarser grid resolutions, with 231 and 29 grid points respectively, are almost identical, indicating a grid independent solution. The same statement is valid for the second-order statistics (not shown in the figures). Wind turning is an important effect in the SABL. Figure 3 shows that the wind direction in the EARS model agrees well with the LES. Differences are only observed in the part of the SABL above the centre of the low-level jet where the wind direction in the EARS model approaches the geostrophic wind direction. In the near-surface region, wind direction and speed are strongly influenced by the boundary conditions. Good agreement between the model and LES indicates that the methodology for imposing the boundary conditions presented in Želi et al. (2019) gives an accurate treatment of the near-surface region. Interestingly, and in contrast to eddy-viscosity-based models, the EARS model predicts a non-zero uv component as a result of the wind turning. Figure 4 shows profiles of TKE. The EARS model with standard values of the coefficients, introduced in Sect. 2, tends to overestimate the amount of TKE in the SABL compared to the LES, and predicts a layer of constant TKE close to the surface. This behaviour is considered to be related to the non-optimal modelling of the -equation (5). The resulting K and K θ profiles are sensitive to the details of the modelling of this equation, as shown later in Sect. 5. This can be expected since generally an underestimated value of leads to overestimated TKE. However, in LES it is not possible to directly evaluate derivatives needed for , which means that this quantity will depend on the subgrid-scale model. For this reason profiles are not compared. Section 5 discusses a possible modification of the equation (5)  Profiles of K θ , shown in Fig. 5, are directly related to the level of stratification in the SABL and the temperature forcing at the surface. Similar to TKE, K θ is higher in the standard EARS model results than in the LES results. In Eq. 4, K θ is coupled to through θ . Therefore, changing the way reacts to the level of stratification in the SABL influences K θ (see Sect. 5 for details). Modifying the prediction of the K θ profile is also possible by changing r , as it represents the ratio between the two time scales K θ / θ and K / . However, in the current model, r is considered to be constant (different values of this constant could be considered).
The EARS model allows the individual stresses to be studied. Instead of plotting uw and vw profiles separately, Fig. 6 shows the total turbulent vertical momentum flux profile defined as uw T = − uw 2 + vw 2 .
Compared to the LES, the character of the profiles is well captured by the model, but the magnitude of the surface stress uw T (z = z 0 ) predicted by the model is greater, which is  Figure 7 shows profiles of the individual normal stresses scaled with TKE, for the lowest cooling rate. In the LES, the variance is computed by adding the isotropic (2/3)K sgs subgrid-scale contribution to the resolved part and the scaling is made with the total TKE. We note the horizontal turbulent momentum variances, uu and vv, in the EARS model are dominant and the vertical variance ww is damped due to stratification effects. Qualitatively, the results of the EARS model agree well with the results of LES also for the normal stress components. Additionally, we note that ww/K should have a finite value above the height of the roughness elements. Although this statement is valid for the LES results, close to the surface the profile of ww exhibits unphysical behaviour in the sense that it decreases near the surface, which is probably related to the surface boundary conditions. Figure 8 shows profiles of the vertical turbulent heat flux. The EARS model predicts that the vertical turbulent heat fluxes behave similarly to the vertical turbulent momentum flux, shown in Fig. 6. The predefined temperature forcing in the EARS and the LES is the same but surface fluxes are different, depending on the surface modelling and the way boundary conditions are imposed. As a consequence, the surface flux wθ ≈ −u * θ * in the standard formulation of the EARS model is greater than in the LES, which is related to the larger overall predicted mixing and in the deeper boundary layer.
Aside from the wθ component of the heat-flux vector, the EARS model predicts the turbulent heat fluxes in x-and y-directions as well, shown in Figs. 9 and 10, respectively. This is an improvement compared to the isotropic atmospheric turbulence models that use the eddy-viscosity concept and which are unable to predict (non-zero) turbulent heat fluxes in the horizontal plane. In the surface layer, the EARS model tends to predict |uθ/wθ| < 1, while Wyngaard et al. (1971) measured |uθ/wθ| > 2, which is closer to the LES result. Above this layer the EARS model performs qualitatively well, changing the sign of the solution in the same way as in the LES model. Results in Figs. 8,9,and 10 show that the turbulent heat flux in the horizontal plane in the EARS model is comparable to the turbulent heat flux in the vertical direction.  Figure 11 shows profiles of the gradient Richardson number, This quantity is sensitive to the turbulence modelling, i.e. the prediction of the mean profile for wind speed and potential temperature, which are results of the turbulence fluxes uw T and wθ. Note that in the EARS model there are no explicitly imposed dependencies between  Fig. 1 turbulent fluxes and the Richardson number. This is obviously not the case in the eddyviscosity models where such dependency usually is imposed through the stability functions. The coefficients in the model for are calibrated in the simulation of the homogeneousshear flow with a constant temperature gradient such that the critical Richardson number is Ri c = 0.20. In the SABL, however, there are turbulent-transport effects present and it is possible to have regions where Ri(z) > Ri c .
As long as there are no transport or transient effects, the change of TKE is ∂ K /∂t → 0 when Ri → Ri c . It is interesting to investigate how the total production of TKE changes as a function of Ri. To study this for different cooling rates, an alignment coefficient C is where II a = a i j a ji and II S = S i j S ji are the second invariants of turbulent anisotropy (10) and mean strain-rate tensor (23). Applying the Cauchy-Schwarz inequality on P gives |P| ≤ K √ II a II S where the equality is reached only when a i j and S i j are perfectly aligned. Misalignment can significantly reduce the production of TKE and even cause the production to become negative in the rare conditions when the TKE is transferred to the mean flow. Therefore, C is a dimensionless measure of local production that relates the alignment between the turbulence anisotropy and the strain-rate tensor with the total production of TKE. Figure 12  shows that C has close relation to the Richardson number inside the SABL in all four cases.
In the LES, C is calculated from the total turbulent quantities (resolved plus sub-grid scale) and reaches zero at Ri ≈ 0.26 implying that above Ri c = 0.25 no turbulence is produced in the LES. Results for the critical Richardson number in the LES agrees with the value reported in Howard (1961) and Miles (1961). One should be aware that C becomes singular at this limit and that the SABL is not in steady-state. Results of the EARS model scale and show a nearly linear dependence on the Richardson number. The alignment coefficient tends towards zero in the higher part of the SABL, close to the inversion layer. In the results of the EARS model, however, C never reaches zero because above the SABL TKE → 0 and C becomes singular. For this reason the profiles were cutoff at the height outside the SABL where the turbulence is weak and the condition TKE < 10 −5 m 2 s −2 is met. Close to the surface the stratification effects weaken so that G → 0 and Ri → 0. Using the data for a neutrally-stratified logarithmic layer from Wallin and Johansson (2000) in (35) results in C ≈ 0.76 near the surface. This agrees with the prediction of the EARS model C = 0.76, but the LES gives a significantly lower value. The reason for this could be that the LES, in the near-surface region, is strongly influenced by the boundary conditions and the value of ww/K in the surface layer is underpredicted (see Fig. 7) resulting in high values of II a . This leads to low values of C.

Discussion
The model presented so far is basically calibrated using generic flow cases such as homogeneous turbulence and passive scalars. The terms related to buoyancy were simply calibrated to give a realistic Ri c . Although the EARS model captures most of the features of the SABL, there are some quantitative differences compared with LES. In general, the heat flux was found to be higher in the model results, leading to a higher SABL. Hence, in this section we will discuss possible modifications and recalibration. However, we avoid any modifications that would influence predictions in the passive scalar limit. We propose two modifications that change the behaviour of the standard model by changing the coefficients in (5) and by introducing the Daly-Harlow transport model.

Modified C 2 and C 3
Equation 5 is modelled in analogy with the exact transport equation for and contains model parameters that need to be calibrated. Different values for C 2 and C 3 coefficients have been proposed in studies of stratified flows with K − model. Pereira and Rocha (2010) made C 2 a function of turbulence anisotropy and set C 3 = 0. In Goudsmit et al. (2002) C 3 < 0 was used for stable stratification, and by using data from a direct numerical simulation Venayagamoorthy et al. (2003) suggested that C 2 is a function of Ri, and C 3 a function of Froude number. Similar to the model suggested by Venayagamoorthy et al. (2003), we here propose to model C 2 as a function C 2 (Ri) with the following functional dependency, for 0.00 ≤ Ri < 0.13 1.82 + 2.57(Ri − 0.13), for 0.13 ≤ Ri ≤ 0.20 2.00, for 0.20 < Ri, and use C 3 = −0.70. The purpose of introducing the functional dependence of C 2 on Ri and modifying the value of C 3 is to study how such a modification in (5) affects the overall behaviour of the model. The seemingly higher values of the K and K θ , seen in Figs. 4 and 5, respectively, were discussed in Sect. 4. Functional dependency for C 2 is chosen such that the value of C 2 does not change in the case of neutral stratification and increases linearly up to 2.0 in the Ri range 0.13 -0.20. The critical Richardson number is thereby calibrated to be Ri c = 0.20, i.e. the same as in the standard EARS model. We note that the total change of C 2 is 10%. Figures 1-11 show the results of the modified EARS model for different stratification levels as dashed curves. Important differences as compared to the standard EARS model are in the levels of K and K θ , shown in Figs. 4 and 5. K and K θ are very sensitive to the modelling in (5) and responsible for intensity of the turbulent fluxes, shown in Figs. 6,8,9. Weaker turbulent fluxes (also surface fluxes, Figs. 6, 8) lead to the lower height of the SABL. In the simulation with the strongest cooling rate Cr = 1 K h −1 the model is predicting lower values of all mean atmospheric properties than in the LES. Interestingly, Fig. 11 indicates that the profiles of Ri also get significantly closer to the LES results, as compared to the standard EARS model. The new coefficients do not change the results shown in Fig. 12 and are hence not included there.

Implementing Daly-Harlow Transport Model
The standard EARS model, described in Sect. 2, models the diffusion in terms of an eddyviscosity-based diffusivity coefficient. However, this diffusion might be inaccurate because turbulence in the SABL is strongly anisotropic. Therefore, the Daly-Harlow diffusion model (Daly and Harlow 1970) can be used for the modelling of the diffusion terms D K , D K θ , and D in (3)-(5). According to the Daly-Harlow transport model: where c s K = c s K θ = 0.16 and c s = 0.21. Modelling of the diffusion terms will have an impact on Ri but will leave Ri c unchanged. It should be noted that the Daly-Harlow model does not impose any additional computational cost on the turbulence model. Figure 13 shows the comparison between the simple diffusion models (6) and Daly-Harlow model (36) used for transport of TKE. Particularly, in (6) we use ν T from the EARS model (see Želi et al. 2019) and the standard K − model approach with C μ = 0.09. The results are based on the TKE profile shown in Fig. 4. Using the same TKE profile ensures that the transport models are comparable. The results of TKE transport using ν T from the EARS model are closer to the Daly-Harlow model than the model where ν T is computed with the formulation from the K − model. This is probably related to the fact that ν T in the EARS model is a function of buoyancy in the EARS model (see Želi et al. 2019). The difference is consistent across all stratification rates and is similar for the transport of half the temperature variance (not shown). Figure 14 shows vertical profiles of the Richardson number, obtained using the Daly-Harlow model for diffusion and with the coefficients given in Sect. 2.2. The mean profiles are nearly the same as the ones described in Sect. 4, however, subtle differences in the mean profiles lead to qualitatively different profiles of Ri in the region around the low-level jet. Other atmospheric quantities are not significantly affected; this is consistent with the previous observation that there are only small differences between the Daly-Harlow and the current  Sullivan et al. (2016) diffusion model, so the results are not shown. Therefore, there is little value of using a more complex diffusion model.

Conclusion
The SABL with weak to moderate stratification is simulated with the EARS model. The results of the model are compared with high-resolution LES results. The mean profiles show that the EARS model is able to predict the general properties and features of stratified turbulence in such situations, including e.g. the low-level jet and wind turning. The mean profiles and vertical turbulent fluxes compare qualitatively well with the LES but the fluxes are greater in the EARS model, which leads to a higher SABL than in the LES. The weakest point in the EARS model is probably the modelling of the equation. It was shown that improvements can be made by modifying the coefficients, although a more general solution to this problem would be preferred. In order to improve the model for the equation it would be helpful to study dissipation in stratified boundary layers by direct numerical simulation and experimental atmospheric measurements.
Unlike many other models, the EARS model predicts all components of the turbulence momentum-flux tensor and heat-flux vector, including horizontal fluxes, without the need for any additional modelling. Apart from the results in the surface layer the horizontal-heat-flux profiles show qualitatively good agreement. This is an advantage compared to the simpler atmospheric turbulence models that are not able to predict horizontal fluxes in the singlecolumn model approach. The model predicts that the horizontal turbulence fluctuations carry more energy than the vertical fluctuations, which is in line with the observations of the SABL.
The study also presents a new type of scaling for the production of TKE as a function of Richardson number in the SABL. Results for different cooling rates in the case of the EARS and LES results collapse well, showing that the scaled production is dependent only on the local Ri and almost independent of the cooling rate. The collapse can be further investigated to see if the scaling has universal features.
Finally, the scaling properties for wind speed and potential temperature are defined from the neutrally-stable logarithmic law relations as where wθ(z 1 ) is the kinematic heat flux obtained from the solution of the EARS model at the first point.