Study of Transitions in the Atmospheric Boundary Layer Using Explicit Algebraic Turbulence Models

We test a recently developed engineering turbulence model, a so-called explicit algebraic Reynolds-stress (EARS) model, in the context of the atmospheric boundary layer. First of all, we consider a stable boundary layer used as the well-known first test case from the Global Energy and Water Cycle Experiment Atmospheric Boundary Layer Study (GABLS1). The model is shown to agree well with data from large-eddy simulations (LES), and this agreement is significantly better than for a standard operational scheme with a prognostic equation for turbulent kinetic energy. Furthermore, we apply the model to a case with a (idealized) diurnal cycle and make a qualitative comparison with a simpler first-order model. Some interesting features of the model are highlighted, pertaining to its stronger foundation on physical principles. In particular, the use of more prognostic equations in the model is shown to give a more realistic dynamical behaviour. This qualitative study is the first step towards a more detailed comparison, for which additional LES data are needed.


Introduction
Turbulence modelling has been an important subject for the atmospheric sciences in order to describe the significant turbulent transport of heat and momentum in the atmospheric boundary layer (ABL, Holtslag et al. 2013). Furthermore, turbulence models have important applications in engineering, e.g. for calculating the drag on ground vehicles or aircraft. The most basic assumption present in many models is the eddy-viscosity hypothesis (first formulated by Boussinesq 1877), which states that the turbulent fluxes are related to the mean gradients of wind (and temperature) through the eddy viscosity (or eddy diffusivity), in analogy to molecular diffusion. Further insights into the structure of turbulent flows, such as the mixing-length model of Prandtl (1925) and energy considerations of Kolmogorov (1941), have further contributed to the development of these models. Here one differentiates between models that aim to capture the average effect of all turbulence scales (Reynolds-averaged Navier-Stokes or RANS models) and large-eddy simulation (LES) models where only the subgrid scales are modelled.
In particular, many models used in climate simulations or weather prediction are based on the eddy-viscosity concept, and buoyancy effects are often included in the expressions for the turbulent fluxes through empirical functions based on the similarity theory of Monin and Obukhov (1954), though strictly valid only for the surface layer. An important quantity in these models is the turbulence length scale, which is often based on the empirical model proposed by Blackadar (1962). However, the eddy-viscosity hypothesis is not sufficient for correctly describing the turbulent fluxes. In particular, it does not give a correct description of the anisotropy of the Reynolds stresses u i u j . A more general approach was described by Yamada (1974, 1982), who derived a well-known hierarchy of turbulence models based on transport equations for the Reynolds stresses and turbulent heat flux u i θ in the ABL. Later developed models are also based on this model hierarchy (e.g. Andrén 1990). However, due to practical limitations in climate and weather-prediction simulations, only the simplest models in this hierarchy are considered in practice, and a minimal number of prognostic equations is used, mostly for the turbulent kinetic energy (TKE). Recent insights have shown the need for additional prognostic equations, in particular for the turbulent potential energy (TPE, Mauritsen et al. 2007;Zilitinkevich et al. 2007Zilitinkevich et al. , 2013. The effects of buoyancy and density stratification, which are crucial in the ABL, are generally speaking less important on the scale of engineering applications. Not taking into account buoyancy amounts to a significant simplification of the equations of turbulence [compare e.g. Wallin and Johansson (2000) and Lazeroms et al. (2013)]. Popular examples of turbulence models used in engineering applications are the K − ε model (Jones and Launder 1972) and the K − ω model (Wilcox 1993), which use two prognostic equations for TKE and, respectively, the dissipation rate or the inverse time scale of turbulence. 1 Similar models have been used in atmospheric contexts as well (e.g. Duynkerke 1988; Apsley and Castro 1997;Sogachev et al. 2012), but to a far lesser extent. The standard formulations of these two-equation models are still based on the eddy-viscosity hypothesis, but substantial progress has been made in recent years to derive more general models from the transport equations for the Reynolds stresses and the turbulent heat flux.
A key concept here is the so-called weak-equilibrium assumption (Rodi 1972(Rodi , 1976, which neglects advection and diffusion of certain dimensionless quantities. This assumption yields algebraic equations for the turbulent fluxes (similar to the Mellor and Yamada hier-archy), and mathematical techniques have been developed to obtain an explicit model for u i u j and u i θ from these equations in a coordinate-free form (i.e. in terms of tensors rather than components with respect to a certain coordinate system). Notable examples for nonbuoyant flows that are well-established are the explicit algebraic Reynolds-stress (EARS) model developed by Wallin and Johansson (2000) and the passive-scalar model by Wikström et al. (2000). The non-buoyant EARS model has also been developed independently by Ying and Canuto (1996) and applied to the atmosphere.
Recently, attempts have been made to extend this more general modelling approach to turbulent flows with buoyancy, i.e. with temperature or density as an active scalar, thus bringing the approach back to the domain of atmospheric models. The inclusion of an active scalar further complicates the mathematical derivation of the model due to the coupling between the Reynolds stresses and the heat fluxes that the buoyancy terms create. Examples of models that were developed through this approach can be found in So et al. (2002So et al. ( , 2004 and Vanpouille et al. (2013Vanpouille et al. ( , 2015. In our study, however, we focus on the model framework developed in Lazeroms et al. (2013Lazeroms et al. ( , 2015, which according to the current authors has led to the first model of this type that is self-consistent and fully explicit (i.e. completely expressed in terms of the prognostic variables). In the previous work, this model has been thoroughly tested in canonical test cases such as homogeneous shear flow and turbulent channel flow, both with stable stratification. It should be noted that similar approaches have been used in LES in both the atmospheric and engineering communities (Findikakis and Street 1979;Wyngaard 2004;Rasam et al. 2013Rasam et al. , 2014Enriquez and Street 2014).
Although EARS models have proven themselves in engineering applications, deriving a model from more fundamental principles and testing it in canonical test cases does not imply that it is more accurate in more complex contexts such as the atmosphere. Such cases might contain more physical phenomena that have not yet been considered in the model. Therefore, the question arises how these newly developed models based on more general principles perform when applied to the ABL. The current study is a first step in answering this question. The aim is to test the new explicit algebraic model presented in Lazeroms et al. (2013Lazeroms et al. ( , 2015 in two test cases, the first of which is the GABLS1 case 2 Beare et al. 2006), where a purely stable boundary layer is formed. This case has particular attention because the stable boundary layer is currently one of the difficult issues in atmospheric turbulence models (Holtslag et al. 2013;Svensson and Lindvall 2015). In the second test case, we consider an ideal diurnal cycle forced by a fully periodic potential surface temperature.
In the following section, we give an overview of the hierarchy of single-column RANS models in order to show the (often detailed) differences and similarities that we would like to address. Section 3 gives a more in-depth description of the model by Lazeroms et al. (2013Lazeroms et al. ( , 2015. In Sect. 4 we present the results for the GABLS1 test case, which are compared with data from LES and two other turbulence closures implemented in an operational weather-prediction model, namely the ARPEGE model at Météo-France (Pailleux et al. 2000). Section 5 covers the case of the (idealized) diurnal cycle, in which we compare our model with a simpler first-order model. This comparison is mainly qualitative due to the lack of LES data for this specific case.

Overview of Turbulence Models
We consider models derived from the Reynolds-averaged Navier-Stokes (RANS) equations, in which the variables have been decomposed into mean and fluctuating parts. The starting point takes the equations for the mean velocity U = (U, V, W ) and the mean potential temperature Θ. In the case of the atmospheric boundary layer, the mean flow equations can be written in the following simplified form, where f = 2Ω sin ϕ is the Coriolis parameter depending on the rotation rate Ω of the Earth and the latitude ϕ, U g = (U g , V g ) is the geostrophic wind (i.e. the large-scale pressuregradient force), and D/Dt ≡ ∂/∂t + U k ∂/∂ x k is the material derivative in the direction of the mean velocity. The quantities u i w and wθ are components of the turbulent momentum flux (or Reynolds stresses) and the turbulent heat flux, respectively, which contain unknown correlations of the velocity fluctuations (u, v, w) and the temperature fluctuations θ . Note that (1) only accounts for net turbulence transport in the vertical direction and assumes horizontally homogeneous conditions for the turbulence statistics. When (1) is solved on a vertical grid only, it is called a single-column model (SCM). The aim is to find a model for the turbulent fluxes in order to close the system of Eq. (1), i.e. the turbulent fluxes need to be expressed in terms of the mean velocity and mean potential temperature.

The Eddy-Viscosity and Eddy-Diffusivity Approach
The first assumption that is often made in turbulence models uses the concept of eddy viscosity or eddy diffusivity to directly relate the turbulent fluxes to the mean gradients, where ν t is the eddy viscosity and κ t the eddy diffusivity of heat. They are the exchange coefficients of momentum and heat due to turbulent eddies, in analogy to molecular diffusion. The eddy-viscosity hypothesis was first formulated by Boussinesq (1877) and is still widely used, although it is not a generally valid assumption, as we discuss in the next sub-section. Finding a closure for Eq. 1 now amounts to modelling the eddy coefficients ν t and κ t ; a simple scaling analysis suggests where q and l are characteristic velocity and length scales of the turbulent eddies, respectively. The fact that we actually need to model the scales q and l leads to the following hierarchy of turbulence models. The simplest class of models are the so-called first-order models (or zero-equation models), where purely algebraic expressions for q and l are used. The eddy coefficients can then be written as, where l m , l h are (possibly different) length scales for momentum and heat, and the velocity scale has been modelled by assuming q = l m |∂U/∂z|. The coefficients f m , f h are stability functions of e.g. the gradient Richardson number, where g is the acceleration due to gravity and T 0 is a reference temperature. These functions are mostly modelled through empirical considerations and differ a lot between different models. The same holds for the length-scale models; one basic form for the length scale that is often used in atmospheric turbulence models is due to Blackadar (1962), which tends to κz (where κ = 0.4 is the von Kármán constant) in the surface layer and to the asymptotic value λ 0 higher up in the ABL. Modifications of this basic form exist, as well as completely empirical functions (see e.g. Yamada 1974, 1982;Andrén 1990). First-order models as described here are still used in recently developed climate models (e.g. Neale et al. 2010Neale et al. , 2012Hazeleger et al. 2010) and numerical weather prediction (NWP) models (Sandu et al. 2013). The next step in the hierarchy is the one-equation model where one of the scales q or l is determined from a prognostic equation. In atmospheric models, the most commonly used equation is that for TKE, denoted by K = (u 2 + v 2 + w 2 )/2, where P is shear production, G is buoyancy production, ε is the dissipation rate and D is a diffusion term. 3 This equation is used to determine the velocity scale q ∼ √ K . In order to close Eq. 7, the dissipation rate ε is modelled through where l ε is the dissipation length scale, for which empirical expressions similar to l m and l h are used. The diffusion term D can be modelled by e.g. gradient diffusion, where σ K is the so-called Schmidt number for TKE. The TKE equation together with algebraic models for l m , l h , and l ε (e.g. Eq. 6) have recently been adopted for the first time in a global model for operational NWP, namely the ARPEGE model at Météo-France , and will also be used in the climate version of ARPEGE for the next Coupled Model Intercomparison Project (CMIP6, see Meehl et al. 2014). Adding another prognostic equation, one obtains the so-called two-equation models, where the best known is the K − ε model, which uses Eq. 7 for TKE and a prognostic equation for the dissipation rate ε instead of the algebraic assumption (8), where C ε1 , C ε2 , C ε3 are constants and D ε is a diffusion term modelled in a similar way as D.
To determine the eddy coefficients in (3), the velocity scale is again taken to be q ∼ √ K and the length scale depends on ε by means of l ∼ K 3/2 /ε, which yields, where C μ is a constant and Pr t is the turbulent Prandtl number, often assumed constant. Since the 1970s (see e.g. Jones and Launder 1972), the K − ε model has been used extensively in many computational fluid dynamical calculations. Another example of a popular twoequation model in engineering flows is the so-called K − ω model, where a prognostic equation for an inverse time scale ω is used instead of ε (Wilcox 1993). Both two-equation models have also been studied in the context of the ABL (e.g. Duynkerke 1988; Apsley and Castro 1997;Sogachev et al. 2012), but their use in atmospheric applications has been problematic. As pointed out by Sogachev et al. (2012), the incorporation of buoyancy effects into these models is not trivial. Other two-equation models for environmental flows make use of a prognostic equation for the product of K and l (e.g. Mellor and Yamada 1982;Blumberg et al. 1992).

Reynolds-Stress Modelling
As noted above, the eddy-viscosity approach is not generally valid and does not yield an appropriate model for complex flow cases. In particular, it does not correctly take into account the anisotropy of the Reynolds-stress tensor u i u j and inter-component exchange of kinetic energy. Therefore, a more sophisticated class of models can be found by considering the transport equations for all components u i u j , as well as the turbulent heat flux u i θ , where in each of these equations, we have the following terms (from left to right): advection, diffusion, shear production, pressure redistribution, dissipation and buoyancy production. The production terms have the following form, with δ i j the Kronecker delta and K θ = θ 2 /2 is one half of the temperature variance. These production terms are in closed form, while the other terms in (12), most notably the pressureredistribution terms, contain further unknown correlations that need to be modelled. For the dissipation term in (12a) one often chooses an isotropic expression: ε i j = 2 3 εδ i j . The pressure-redistribution terms Π i j and Π θi describe inter-component exchange due to pressure fluctuations. Different models exist, though Π i j is usually decomposed into a slow part (or return-to-isotropy term, Rotta 1951) and a rapid part directly dependent on the mean gradients. A general linear model for Π i j was formulated by Launder et al. (1975) and extended with buoyancy-dependent terms by Launder (1975). A similar model for Π θi − ε θi can be found in e.g. Wikström et al. (2000). For the diffusion terms D i j and D θi (generalized) gradient diffusion models can be used (e.g. Daly and Harlow 1970). To completely close Eq. 12, additional prognostic equations or models for ε and the temperature variance θ 2 (as well as its dissipation rate ε θ ) are needed. Note that Eq. 7 for the TKE is obtained by taking half the trace of (12a).
Despite the high level of sophistication of Eq. 12, these so-called differential Reynoldsstress models (DRSM) are difficult to evaluate in practice due to the large number of prognostic equations. It is therefore desirable to simplify these equations to a purely algebraic form. Various approaches exist to accomplish this, which all amount to modelling the advection and diffusion terms on the left-hand sides of (12). The simplest way to accomplish this is to neglect these terms altogether. Although this assumes a steady-state for the turbulent fluxes that is not generally valid, this approach is still used in recently developed models (e.g. Zilitinkevich et al. 2013). A slightly different approach is that of Cheng et al. (2002), which neglects the advection and diffusion terms of the anisotropic part of the Reynolds stresses, b i j = u i u j − 2 3 K. One can argue that neglecting any of these dimensional quantities can lead to problems in more complex situations. Instead, an analysis based on dimensionless fluxes is preferred. A well-known method in atmospheric turbulence models is that presented by Yamada (1974, 1982), who consider Eq. 12 in terms of the following dimensionless quantities, where the tensor a i j is called the (dimensionless) Reynolds-stress anisotropy and ξ i is a dimensionless heat flux normalized with TKE and the temperature variance. Mellor and Yamada perform a scaling analysis of the equations for a i j , ξ i , K and K θ , and neglect terms based on powers of the anisotropy. Neglecting terms of different order leads to a hierarchy of models of different complexity, leading from the full prognostic equations (level 4) back to the standard first-order eddy-diffusivity models (level 1). This hierarchy serves as the basis of many later atmospheric turbulence models (e.g. Andrén 1990). Note that the dissipation of TKE is still modelled by means of a length-scale parametrization. The common approach in engineering, however, is based on the work by Rodi (1972Rodi ( , 1976, who postulated that the advection and diffusion terms of dimensionless quantities can be neglected. This is called the weak-equilibrium assumption. In practice, this neglects only advection and diffusion terms of a i j and ξ i as defined in (14). Going back, for example, to Eq. 12a for u i u j , this assumption amounts to, i.e. temporal and spatial variations of u i u j are directly determined by temporal and spatial variations of TKE. This assumption is clearly more general than neglecting advection and diffusion of u i u j completely. Data from direct numerical simulations, in which all turbulence is resolved, suggest that the weak-equilibrium assumption is valid in many contexts (e.g. Vanpouille et al. 2013), though it fails most notably in near-wall turbulent regions, which are usually not resolved in atmospheric models. The weak-equilibrium assumption results in a model consisting of algebraic equations for u i u j and u i θ together with prognostic equations for K , ε and K θ , which would be most similar to the Mellor and Yamada level-3 scheme (e.g. Niino 2004, 2006). However, there is a fundamental difference between the weak-equilibrium assumption and the Mellor and Yamada approach, leading to models of somewhat different complexity, as explained in Appendix 2.
Much work now focusses on finding the explicit solutions to the aforementioned algebraic equations for u i u j and u i θ . For practical reasons, it is important to solve these equations algebraically and not use any numerical iterations at this level in the model. In many atmospheric models, this is less of an issue because they are formulated with respect to a certain coordinate system (i.e. the geometry of the ABL, see Eq. 1), so that one can solve the algebraic equations for each component of u i u j and u i θ separately (see e.g. Mellor and Yamada 1974). The more general models, however, aim for a coordinate-free formulation without reference to a certain coordinate system and suitable for use in arbitrary geometries found in engineering applications. This requires more advanced mathematical considerations that are beyond the scope of our study. The resulting models for u i u j and u i θ are called explicit algebraic Reynolds-stress models, examples including the model of Wallin and Johansson (2000) for neutral turbulent flow and of Wikström et al. (2000) for passive scalar transport. Extensions to buoyant flows/active scalars were made by e.g. So et al. (2002So et al. ( , 2004, Lazeroms et al. (2013Lazeroms et al. ( , 2015 and Vanpouille et al. (2015). The differences between these models mainly consist of details in the mathematical framework that are not discussed here. Rather, we make use of the model derived in Lazeroms et al. (2013) and of which some aspects were extended and improved in Lazeroms et al. (2015)-for a more detailed description see Sect. 3 and Appendix 1. We refer to this model as the current model.

Turbulent Potential Energy and Counter-Gradient Heat Flux
The EARS model discussed in Sect. 2.2 constitute one of the most advanced classes of turbulence models used in engineering. On the other hand, many atmospheric models still make use of a standard TKE scheme (i.e. Eq. 7). However, a recent development in atmospheric turbulence modelling Zilitinkevich et al. 2007Zilitinkevich et al. , 2013 is the addition of a prognostic equation for the turbulent potential energy (TPE), defined as, The main idea behind these models is the conversion of kinetic to potential energy by buoyancy, whereas in standard TKE schemes the buoyancy term G would just act as a source or sink of energy. The aim of this approach is to predict the presence of turbulence in the strongly stable regime, the possible existence of which has recently been pointed out in observations Tjernström et al. 2009). In standard models, on the other hand, there exists a critical Richardson number above which turbulence is completely damped. These considerations have recently led to turbulence models using a single prognostic equation for the total turbulent energy Pithan et al. 2015), as well as the so-called Energy-and Flux-Budget (EFB) closure that uses separate prognostic equations for TKE and TPE (Zilitinkevich et al. , 2013. The latter is also derived from the momentum and heat-flux budgets in Eq. 12, albeit with more empirical considerations than the EARS model. On the other hand, the Reynolds-stress models described in Sect. 2.2 naturally depend on the temperature variance θ 2 = 2K θ , a quantity that is essentially equivalent to the TPE in Eq. 16. This dependence comes from the buoyancy term G θi in Eq. 12b and the scaling used for the turbulent heat flux (Eq. 14). An additional prognostic equation for K θ is easily added to models like EARS, which can be written as, where P θ is the production term, ε θ is the dissipation rate and D θ is a diffusion term. Interestingly, the possible existence of turbulence above the critical Richardson number has not been an issue in engineering turbulence models. Indeed, the EARS models described above appear to have a critical Richardson number close to the theoretical value of 0.25 known from the analysis of Miles (1961) and Howard (1961) (also discussed in Lazeroms et al. (2013) and in Sect. 4). At the same time, these models are derived from rather basic principles such as general expressions for the pressure terms Π i j and Π θi and the weakequilibrium assumption, while models such as the EFB closure make use of more empirical relations and ad hoc corrections in e.g. these pressure terms. It is therefore likely that such corrections are always necessary in order to obtain a turbulence model without a critical Richardson number. Another interesting feature of using an additional prognostic equation for K θ or TPE is the appearance of a counter-gradient heat-flux term, whereas standard eddy-diffusivity models have a heat flux that is always proportional to the mean temperature gradient. In the current model, for example, the turbulent heat flux can be expressed as, where κ t is an effective eddy diffusivity and Φ cg is the counter-gradient term proportional to K θ . The appearance of the counter-gradient term can be of high importance in regions with strong temperature fluctuations (i.e. large K θ ) and a positive temperature gradient, for example in the inversion layer. On the other hand, in models using a steady-state expression for K θ (see the Mellor and Yamada level 2, or Andrén 1990), the term Φ cg becomes proportional to the mean temperature gradient, so that it can be incorporated into the κ t term. A prognostic equation for K θ such as (17) would therefore be necessary for obtaining the correct dynamics in the counter-gradient term.

Description of the Current Model and The Numerical Solver
In the following, we focus on the EARS model presented in Lazeroms et al. (2013Lazeroms et al. ( , 2015 and apply it to a number of test cases in the ABL, in order to compare its performance with some of the other modelling approaches considered in the previous section. Here we summarize the main aspects of the model specifically for these ABL cases. The ABL is simulated by solving Eq. 1 for the mean wind and mean temperature. The models for the Reynolds-stress components in these equations can be expressed by (2a) and (2b) and the model for the turbulent heat flux by (18). The effective eddy viscosity ν t , effective eddy diffusivity κ t and the counter-gradient flux term Φ cg can be written symbolically as, In contrast to the empirical expressions used in standard eddy-viscosity models (e.g. Eq. 4), the functions f m , f h , f Φ now follow directly from the derivation outlined in Sect. 2.2. These expressions depend on a number of invariants of the flow and are given in full detail in Lazeroms et al. (2013Lazeroms et al. ( , 2015, though some aspects of their structure are discussed in Appendix 1. The model parameters in this part of the model are contained in the expressions f m , f h , f Φ and are entirely due to the (general) models used for the pressure-redistribution terms Π i j and Π θi in (12). These parameters have been calibrated for various flow cases in previous work Wikström et al. 2000;Lazeroms et al. 2013). Note that apart from the turbulent fluxes in (1), our general derivation also yields expressions for the other components of the (anisotropic) Reynolds stresses u i u j and turbulent heat flux u i θ , as explained in Appendix 1.
In order to close these expressions, we use the K −ε model discussed in Sect. 2.1, together with a prognostic equation for K θ as shown in (17). This leads to the following, with C ε1 = 1.44, C ε2 = 1.92, C ε3 = −0.8, σ K = σ K θ = 1.0, σ ε = 1.3 and r = 0.55. Note that the diffusion terms have been modelled using gradient diffusion as in Eq. 9. Furthermore, the dissipation rate of K θ (second term on the right-hand side in (20c)) is parametrized by ε θ = K θ /(r τ ) with relaxation time scale τ = K /ε, equal to the time scale of the velocity fluctuations. The value of r is based on empirical considerations from direct numerical simulations of Johansson and Wikström (1999). The rest of the parameter values in (20) are standard values for the K − ε models found in the literature, except for C ε3 , for which no standard value has been found. Here we treat this constant as a free parameter that we tune 4 for an optimal performance based on the case described in Sect. 4. This approach is simpler than e.g. that of Sogachev et al. (2012), who developed a consistent way of incorporating buoyancy (as well as vegetation effects) into standard two-equation models for use in ABL applications. However, a detailed investigation of the equations for the turbulence scales is beyond the scope of this work, as the main feature of the EARS model is the explicit algebraic formulation of u i u j and u i θ . Here we regard the K − ε − K θ equations as a secondary, yet necessary, part of the model. The constant C ε3 in the ε equation is not expected to be universal, and further investigations with more advanced models for the turbulence scales (such as that of Sogachev et al. 2012) are desirable.
Note that the current model is slightly different from that presented in Lazeroms et al. (2013Lazeroms et al. ( , 2015, where a so-called K −ω model was used for reasons connected to the near-wall treatment in engineering flows (also explained in Lazeroms et al. 2015). More precisely, the standard K − ω model is known to behave more realistically when the model is evaluated in the viscous sublayer and exactly at the surface. On the other hand, the standard K − ω model is known to have issues in boundary-layer flows in the region of transition to freestream conditions. Since we use logarithmic layer conditions instead of resolving the viscous sublayer (see below), the K − ε model seems to be the best choice for the ABL cases considered here. One should note, however, that these properties depend on the exact form of the two-equation model, since Sogachev et al. (2012) found that their K − ω model behaves more appropriately than their K − ε model. The system of differential equations (1) and (20) is solved on a domain [z 0 , H ], where z 0 is the aerodynamic roughness length at the surface and H is the domain height. By taking the lower boundary of the domain at z = z 0 , we can easily define the following boundary conditions for the mean flow equations (1): where Θ surf is a problem-dependent function of time. For the K − ε − K θ model in (20), we assume so-called logarithmic layer conditions, which can be derived from a balance of (shear) production and dissipation terms in the logarithmic layer, using a standard eddy-viscosity model for simplicity. These boundary conditions read as follows, with C μ = 0.09, Pr t = 0.9 and k = 0.4. The friction velocity u * and friction temperature θ * are determined by a linear fit of the following logarithmic profiles near the surface, Note that these boundary conditions are somewhat simpler than the often used empirical relations based on Monin-Obukhov similarity, which also take into account buoyancy production (Monin and Obukhov 1954). The present approach is justified if the grid size used for fitting the logarithmic functions near the surface is sufficiently small for neglecting buoyancy effects. This means that the first grid points at the surface have to be within the logarithmic layer in order to obtain the correct slope in the velocity profile near the surface (cf. Fig. 3) and, subsequently, appropriate values for u * and θ * from Eq. 22. At the top of the domain, the velocity equals the geostrophic wind velocity (U = U g and V = V g ) and a fixed value of the potential temperature gradient (corresponding to the initial profile). For the quantities in (20), we take ∂/∂z = 0 at the top. In this study, the same numerical solver as in previous work is used 5 (Lazeroms et al. 2013(Lazeroms et al. , 2015. Customized grid points, timesteps and initial conditions are added to complete this method. In the test cases below, the grid resolution near the surface is chosen fine enough to obtain a numerical solution that is independent of the grid, i.e. yielding velocity and temperature profiles consistent with the logarithmic layer boundary conditions discussed above (see also Fig. 3). In such a way, we can focus on the physical correctness of the model rather than discrepancies due to the numerics.

Model Comparison for a Stable Boundary Layer (GABLS1)
The first test case considered here is the GABLS1 case used in the intercomparison study by Cuxart et al. (2006). This represents a purely stably stratified ABL, which is formed by applying a constant geostrophic wind speed and constant cooling at the surface. The main point of interest is to investigate how well the model performs in predicting stably stratified turbulence in the atmosphere by comparing with LES data available from the GABLS1 study .
The constant geostrophic wind speed in this case is U g = 8 m s −1 with the wind vector directed in the x-direction. The case is set in the Arctic at a latitude of 73 • N, corresponding to f = 1.39 × 10 −4 s −1 , and the potential temperature at the surface Θ surf = Θ(z = z 0 ) is given by, with Θ 0 = 265.1 K (including a small correction with respect to the original data 6 ) and b = 0.25 K h −1 . The aerodynamic roughness length z 0 is taken equal to 0.1 m. For the initial profiles and other parameters, we refer to Cuxart et al. (2006). In addition, the initial ε profile is chosen such that the turbulence time scale τ = K /ε initially equals 1 s, while the initial K θ profile is put equal to zero. The model is now evaluated on a domain with height H = 400 m and 127 grid points, with the grid size ranging from 0.1 m near the surface to 10 m at the top. As in the original GABLS1 case, we run the simulation for 9 h (model time), after which a quasi-steady state is reached; the timestep in the simulation is 60 s. At the end of the simulation, mean values over the last hour are taken. Similar computations have been performed with the ARPEGE model at Météo-France, using two different turbulence schemes: the standard (operational) TKE scheme ) and a prototype of the EFB closure described in Sect. 2.3 (coded by E. Bazile, internal communication), both with a vertical grid size of 6.3 m. For the LES data, we have datasets from five different simulations with a 2-m resolution . Figure 1 compares the results of all these computations for the mean wind speed and mean potential temperature results for both the momentum flux and heat flux, as well as the TKE profile. Again we conclude that the current model produces results closer to the LES than the standard TKE scheme, and once more the EFB closure gives very similar results.
To demonstrate the validity of the boundary conditions used in the model, we plot the mean wind-speed and mean potential temperature profiles on a logarithmic scale, shown in Fig. 3. The model results are compared with the EFB closure implemented in the ARPEGE model and the LES dataset that is closest to the current model in Figs. 1 and 2 (i.e. the NCAR dataset, see Beare et al. 2006). Clearly, both the models and the LES obey the same logarithmic law near the surface (up to ≈5 m). This also shows that the logarithmic layer boundary conditions shown in Eqs. 21 and 22 are consistent and would work for the present case as long as the grid size near the wall is smaller than 5 m. Higher up in the boundary layer, the effects of buoyancy and a correction to the logarithmic profile based on Monin-Obukhov similarity are needed (Monin and Obukhov 1954). This correction is a linear function of the dimensionless parameter z/L, where L = −u 3 * T 0 /(kgwθ surf ) = u 2 * T 0 /(kgθ * ) is the Obukhov length (Obukhov 1946) depending only on surface fluxes. The wind profile agrees well with the theoretical curve if the well-known coefficient 4.8 proposed by Högström (1988) is used in the linear term, which also appears to work fairly well for the temperature profile. 7 Another indication of the correct treatment of the near-surface region can be seen in Fig. 4, where we have plotted the hodograph of the mean velocity vector for both the current model and the NCAR LES dataset. Clearly, the wind direction near the surface (making an angle of ≈36 • with the x-direction) resulting from the model corresponds very well to the LES data, and also the wind turning in the rest of the boundary layer is well predicted.
In the context of this stably stratified test case, it is relevant to investigate the issue of the critical Richardson number discussed in Sect. 2.3. In Fig. 5a we show the current model results for the velocity variances plotted against the gradient Richardson number. From Fig. 5a it is clear that all velocity fluctuations decay at a specific critical Richardson number. In this case, the critical Ri value in the current model appears to be close to 0.2, which is slightly lower than the critical value of 0.25 found in Lazeroms et al. (2013) with the K − ω model. , v 2 (green ×) and w 2 (red ). The solid lines in panel a are extracted from the NCAR LES data (same colour coding), which use the subgrid-scale model specified in Sullivan et al. (1994) Interestingly, the NCAR LES data shown in the same figure appear to decay at similar values of Ri, but residual turbulence remains slightly above this critical value. The LES results appear to be in accordance with observations that suggest the existence of turbulence above the critical Ri Tjernström et al. 2009). The present model predicts no sustained turbulence for conditions above the critical Ri, but it allows the existence of turbulence that decays very slowly if the turbulence length scales are large (see Sun et al. (2015) for a review of these issues). Note that there are clear discrepancies between the model and the LES for low Ri. However, these points are very close to the surface where neither the model nor the LES resolves the small scales in the viscous sublayer, so it is hard to draw any conclusions for this region (see also Rasam et al. (2011) for a discussion on the overprediction of near-wall peaks in the velocity fluctuations for different subgrid-scale models). Figure 5b shows the same variances scaled with TKE, which represents the contribution of each component to the total TKE. Generally speaking, one can conclude that the vertical contribution decreases slightly while the horizontal contributions remain around the same value. Interestingly, the horizontal contributions appear to oscillate at the critical Richardson number, most likely due to a slight turning of the wind in the upper part of the ABL, causing different levels of shear production in the x-and y-directions. Not shown in the figure are the data points above z ≈ 225 m where Ri becomes very large and where u 2 = v 2 and w 2 /2 < 0.1K , with K very close to zero. Also note that Fig. 5 reveals the ability of the model to predict anisotropic turbulence with clearly different components of TKE, in contrast to standard models based on the eddy-viscosity hypothesis.

Qualitative Study of an Idealized Diurnal Cycle
We now consider a test case that includes the dynamics of a diurnal cycle; this is partly based on the GABLS2 intercomparison study (Svensson et al. 2011), in which the models were forced by a prescribed diurnal variation of the surface temperature taken from observations in Kansas, USA (CASES-99, Poulos et al. 2002). Here, we investigate a more idealized version of the GABLS2 case that is forced by a completely periodic surface temperature and without the effects of humidity and subsidence. This case was also considered in Lazeroms et al. (2015) in order to validate certain improvements to the model (i.e. the approximation of Eq. 32 in Appendix 1). Such an idealized case is easier to handle in the current model framework and allows us to focus primarily on the physics of turbulent heat transfer, which is the main effect included in the current model. It also allows for simulations that are longer than the three-day period considered in GABLS2. The disadvantage of this approach is that no LES data are available yet for a quantitative comparison, but in the following we focus on qualitative results. We consider a constant geostrophic wind velocity with U g = 3 m s −1 and V g = −9 m s −1 ; the latitude of the GABLS2 case is 37.6 • N, corresponding to f = 8.9 × 10 −5 s −1 . The periodically varying potential temperature at the surface is given by, with Θ 0 = 283 K, ΔΘ = 11 K and t is in h. This function was chosen in order to obtain a temperature variation similar to the final day of the GABLS2 case. For the aerodynamic roughness length we take z 0 = 0.03 m. The initial velocity profile is put equal to the geostrophic wind velocity, and for the initial potential temperature we take a piecewise linear profile defined by the values given in Table 1. The initial TKE is given by K (z) = 0.5(1 − z/200) for z ≤ 200 m and a near-zero value for z > 200 m. The initial ε profile is again chosen such that the initial turbulence time scale is 1 s, and the initial K θ = 0. Other parameters can be found in Svensson et al. (2011) (note again that no effects of humidity and subsidence are used in the current case). The model is evaluated on a domain with height H = 4000 m and 213 grid points, the grid size is 0.5 m near the surface and increases to 100 m at the top of the domain, and the model is run for 480 h (20 days) with a timestep of 60 s. An interesting feature of the current case is the constant amplitude of the periodic surface temperature variation. Initially, the mixed layer (located between z ≈ 100 m and z ≈ 1 km to z ≈ 1.2 km) is heated up rapidly because the maximum surface potential temperature (294 K) is significantly larger than the initial potential temperature in the mixed layer (283 K), as seen in Fig. 6a. Since the amplitude of the potential temperature at the surface is constant throughout the simulation, a quasi-steady state is expected to appear when there is a (a) (b) Fig. 6 Diurnal variation of a the mean potential temperature, and b the turbulent kinetic energy at z = 503 m as calculated by the current model. The equilibrium value of daytime TKE is indicated by the dashed green line balance between daytime heating, nighttime cooling and the effects of entrainment at the top of the boundary layer. This can be seen in Fig. 6a, where the daytime increase of potential temperature in the mixed layer appears to gradually slow down, and in Fig. 6b, where the maximum value of TKE at a fixed height appears to reach an equilibrium value. Further test runs with a length of up to 40 days have shown that this quasi-steady state remains.
We now investigate several qualitative features of the current model. For this purpose, we show the diurnal variation of TKE over three consecutive days in the quasi-steady state in Fig. 7a. For comparison, we also show in Fig. 7b the corresponding results obtained by a typical first-order model, namely the level-2 version of the model framework derived in Yamada (1974, 1982). This model is based on a steady-state expression for the TKE, which is not directly obtained from any prognostic equation. The equivalent TKE can be calculated in the following way, which corresponds to Eq. 59 in Mellor and Yamada (1974). The length scale l in this model is determined by a modified version of the expression proposed by Blackadar (1962) (see Eq. 6), where the asymptotic value l 0 = 100 m. The first-order model has been evaluated with the same solver as used for the EARS model. 8 A comparison of Fig. 7a, b reveals that the two models give quite different levels of TKE and different boundary-layer heights. However, the intensity of the turbulence in the Mellor and Yamada model depends heavily on the chosen expression for the length scale l, making it difficult to compare the models quantitatively without any LES data. Nevertheless, some interesting qualitative features can be discerned. The most apparent is the fact that, compared to the first-order model, the current model gives a much slower decay of the TKE present  Mellor and Yamada (1974). Note that the colour scale is logarithmic. Also shown is c, the variation of the surface potential temperature  Mellor and Yamada (1974) (dashed line) in the mixed layer during nightfall. In other words, there is still weak turbulence left in the residual layer predicted by the current model. The phenomenon can be seen more clearly in Fig. 8, which shows the evolution of TKE in the two models at a fixed height. This effect is  Mellor and Yamada (1974). Note that the colour scale is logarithmic. Also shown is c the variation of the surface potential temperature interesting in the light of recent studies of observational data that suggest the presence of weak turbulence in the residual layer, even above the critical Richardson number (Tjernström et al. 2009). Whether this residual turbulence in the model is truly related to the effects described by Tjernström et al. (2009) remains to be investigated.
The other aspects that we discuss are apparent near the beginning of the simulation before the quasi-steady state has been reached. The evolution of the TKE in both models near the start of the simulation is shown in Fig. 9. The second qualitative feature that requires attention is the relatively smooth transition between the stable and convective boundary layers during each morning. In the current model, the smooth transition is mainly caused by the weak stability of the residual layer (also shown in Lazeroms et al. (2015)). These features have also been found in observations (e.g. Angevine et al. 2001). On the other hand, many firstorder models are known to exhibit a near-instantaneous growth of the boundary layer in the morning because the residual layer remains neutral (Beare 2008;Svensson et al. 2011). Also in this case one can see that the first-order model in Fig. 9b has a somewhat more rapid morning transition than the current model in Fig. 9a.
Moreover, Fig. 9a shows a more dynamical behaviour at the top of the daytime boundary layer (inversion layer) for the current model as compared to the first-order model shown in Fig. 10 Diurnal variation near the start of the simulation of K θ = θ 2 /2 obtained from Eq. 20c in the current model. Note that the colour scale is logarithmic Fig. 11 Diurnal variation near the start of the simulation of the equivalent length scale l eq (Eq. 26) in the current model. Note that the colour scale is logarithmic Fig. 9b, for which the contour plot of TKE appears rather "flat". This effect is likely due to the inclusion of a prognostic equation for the temperature variance θ 2 = 2K θ in the current model and the counter-gradient heat-flux term that it generates, whereas the first-order model is derived from a steady-state expression for this quantity. The importance of the temperature variance in the inversion layer can be seen in Fig. 10. Determining this quantity from a prognostic equation allows for a more dynamically adapting behaviour and might improve the description of entrainment processes in this part of the boundary layer.
As mentioned earlier, many turbulence models currently used in weather prediction and climatology depend heavily on the chosen formulation for the length scale l. Therefore, we investigate how this quantity is represented in the current model. Since the model makes use of the K − ε equations instead of a separate formulation for l, an equivalent expression can be constructed from K and ε. This can be done as follows (cf. Eqs. 8 and 11a), with C μ = 0.09. Figure 11 shows the diurnal variation of this equivalent length scale for three consecutive days at the start of the simulation. One can see that this length scale has a very dynamic behaviour, in contrast to the constant length scale (Eq. 6) used for the Mellor and Yamada model. In practice, however, some atmospheric models tend to use different (empirical) formulations for the length scale in stable and unstable situations, as well as for the surface layer and the upper part of the boundary layer (e.g. Andrén 1990). These different empirical relations are not necessarily continuous when transitions from stable to In such cases, a dynamically adapting length scale as provided by the K − ε model can have significant computational advantages, as well as a more realistic description of the turbulence physics. This also holds for the equations for TKE and the temperature variance. The effect of using more prognostic equations in the current model can also be seen in Fig. 9a, where the results are clearly different each day (i.e. TKE is dynamically adapting during the first half of the simulation, before the quasi-steady state has been reached), while the first-order model results in Fig. 9b are nearly identical each consecutive day. The advantage of using dynamical equations for the turbulence scales has also been noted by Sogachev et al. (2012). Of course, the K −ε framework is not the only way to produce a dynamical length scale, and several dynamical equations for the length scale have been introduced in the past in the atmospheric community (e.g. Mellor and Yamada 1982). In a more quantitative study, it would be helpful to test the EARS model in conjunction with other length-scale formulations as well. See also Sun et al. (2015) for a discussion on length-scale formulations.
To complete the discussion of this test case, we show the time evolution of the turbulent momentum and heat fluxes calculated with the current model in Figs. 12 and 13. These results are taken at the end of the simulation when the system has reached the quasi-steady state. One can observe in Fig. 12a that the momentum flux in the residual layer is zero, while the TKE is non-zero here (cf. Fig. 7), implying that these residual motions are uncorrelated and would naturally decay. Interestingly, Fig. 13a reveals patches of positive but small ( 10 −4 K m s −1 ) turbulent heat-flux in this part of the residual layer, which is weakly stably stratified. We therefore obtain a small positive contribution of the counter-gradient heat-flux term (Eq. 18) in this region. Such a contribution would of course be absent in the first-order model. Apart from this, the first-order model was found to give qualitatively similar results for the momentum and heat fluxes, albeit with lower intensities compared to the current model, as already noted for the TKE in Figs. 7, 8, and 9.

Conclusions
We have applied a newly developed explicit algebraic Reynolds-stress (EARS) model to the atmospheric boundary layer. The EARS model was derived from fundamental principles as described in Sect. 2.2. An essential feature of the model is its ability to capture an anisotropic distribution among the components of TKE, in contrast to standard eddy-viscosity models. This is crucial for the ABL, and stratified flows in general, which can exhibit a high degree of anisotropy. Apart from the evolution description for K , ε and K θ given by Eq. 20, the model framework contains only a few model parameters that appear through the general linear model for the pressure-redistribution terms, and which have been calibrated in Lazeroms et al. (2013Lazeroms et al. ( , 2015 and earlier work. Also the boundary conditions are based on rather general relations for the logarithmic layer. We were able to make a thorough comparison of model simulations with LES data in the GABLS1 case . The results in Figs. 1 and 2 show very good agreement with available LES data (especially the dataset from NCAR) and the agreement appears significantly better than with a standard TKE scheme. The good results of the model for this stably stratified case are consistent with those presented in Lazeroms et al. (2013) for turbulent channel flow. A very interesting result in Figs. 1 and 2 is the close agreement of the current model and the EFB closure as implemented in the ARPEGE model, even though they are based on somewhat different theoretical considerations. Most notably, the EFB closure as presented in Zilitinkevich et al. (2013) appears to have more empirical assumptions where our EARS model is based on more basic principles. A thorough comparison of the two models in different test cases would be an interesting subject of future work, in order to grasp the differences and similarities between them.
Furthermore, we have shown some qualitative results for an idealized diurnal cycle based on the GABLS2 case (Svensson et al. 2011). A comparison with a simpler first-order model has shown some interesting features of the current model, which shed light on the aspects where the model might give an improvement over other models (e.g. the smooth morning transition and the dynamics in the inversion layer). However, the big drawback of this test case is the lack of LES data for a thorough comparison. It would be interesting to conduct a LES study of the currently presented diurnal cycle with sinusoidal surface temperature variation to systematically assess the EARS model and other turbulence closures. A similar idealized case based on strongly stable conditions in the Antarctic is currently being considered as part of the GABLS4 intercomparison project. Furthermore, there are previous LES studies that have focussed only on certain aspects of the diurnal cycle (Beare 2008;Kumar et al. 2010).
One specific result that requires attention is explained by Fig. 5a, which shows that, although the current model is derived from quite fundamental considerations, it predicts a critical Richardson number above which turbulence will not be sustained indefinitely, an important difference with the EFB closure specifically designed for predicting sustained turbulence for high Richardson numbers. The question remains as to whether the regime of weak turbulence found in observations is caused by a phenomenon that in principle cannot be captured by RANS models (cf. Sun et al. 2015), so that further (empirical) corrections are the only means of removing the critical Richardson number. The residual TKE predicted by the current model (Fig. 7a) in the residual layer is also interesting in this respect, since it is a form of (decaying) turbulence present in very stable conditions (cf. Balsley et al. 2008). A crucial aspect of the model might be its ability to correctly predict anisotropic flows, so that the vertical motions are damped much more strongly in stable conditions than are horizontal motions. These anisotropic effects can even be extended by using an anisotropic diffusion length scale for TKE (e.g. Daly and Harlow 1970) instead of the simpler formulation in Eq. 9. Another effect that may be of interest in the issue of the critical Richardson number is the interaction of gravity waves and turbulence, which is normally not captured by turbulence models (Nappo et al. 2014).
The current approach may go beyond the computational limitations of climate models and NWP models. Despite its more fundamental derivation, the current model is not necessarily more accurate than simpler models, which are specifically designed and tuned for use in practical codes. However, a derivation from more general principles should make the model applicable in a wide variety of conditions, in contrast to being tuned for very specific cases. From this point of view, a more practical study of the model's behaviour in a full climate model is also necessary. The current model might prove to be more numerically stable than simpler models because it describes the physics of turbulent transport in the ABL in a more dynamical and continuous way, as discussed in the previous section.
coefficients also depend on the aforementioned invariants, as well as the model parameters appearing in (27). For example, the coefficient β 1 in (29a) can be written as, where C 2 = 4/9 is a model parameter and the factor N is defined as, where c 1 = 1.8 and P, G and ε are terms in the TKE equation (7). As discussed in Lazeroms et al. (2013Lazeroms et al. ( , 2015, the exact value of N obeys a complicated non-linear equation, but an efficient and robust method to approximate this quantity, as well as the factor N θ in (27), has been presented in Lazeroms et al. (2015). Note that we use the expressions for parallel flows, even though the geometry of the (single-column) ABL as given by (1) is technically not a parallel flow. However, as already mentioned in Lazeroms et al. (2015), one can show that the geometric properties of the ABL cases considered here are equivalent to those of parallel flows because the relations between the invariants of the tensor groups in the model are the same. Furthermore, we have assumed that effects due to the Earth's rotation can be neglected in the turbulence model, i.e. the Earth's rotation only affects the mean flow equations through the Coriolis term. This can be justified by considering the Rossby number for the turbulent motions in the ABL, Ro = q/( f l), which is large ( 10) for a velocity scale q ≈ 1 m s −1 , length scale l 10 3 m and Coriolis parameter f ≈ 10 −4 s −1 , as also noted in Lazeroms et al. (2015).

Appendix 2: Distinction Between the EARS Model and the Mellor and Yamada Level-3 Closure
Although the weak-equilibrium assumption underlying the EARS model leads to a set of equations similar to the level-3 closure models in Yamada (1974, 1982), there are some important fundamental differences between the two types of models. Here we explain these differences in more detail focussing on the equations for the momentum flux (similar arguments apply to the heat flux). As explained in Lazeroms et al. (2013), the transport equation for the dimensionless anisotropy a i j (see Eq. 14) can be derived from (7) and (12a) and has the form, in which the left-hand side denotes the advection and diffusion terms of a i j . In the weakequilibrium assumption, each term on the left-hand side is neglected, leading to the approximations in Eq. 15 for the advection and diffusion of u i u j , and the following balance of terms, Since the shear production P and buoyancy production G depend on u i u j and u i θ , the terms on the left-hand side of Eq. 34 are non-linear in the fluxes. This is the reason why, in order to obtain a consistent EARS model, we have to solve the non-linear equation for the factor N (Eq. 32) as explained in Appendix 1 and Lazeroms et al. (2015).
On the other hand, the simplifications made in the level-3 model of Yamada (1974, 1982) are based on a scaling analysis and subsequently neglecting terms of order a 2 , where a is the magnitude of a i j , rather than considering the variations in this quantity. This leads to a different approximation of advection and diffusion of u i u j , viz.
which is an isotropic formulation of these terms, in contrast to the (possibly more general) anisotropic formulation of Eq. 15. In fact, the same approximation would be obtained if one neglects advection and diffusion of the dimensional anisotropy, b i j = u i u j − 2 3 δ i j , as was done by Cheng et al. (2002). Neglecting these terms seems to be less plausible than neglecting advection and diffusion of a i j because the assumption of slow variations in space and time (i.e. the weak-equilibrium assumption) does not readily apply to dimensional quantities. Hence, the Mellor and Yamada approach is strictly valid only in strong equilibrium where also the turbulence scales are in equilibrium. The consequence is that rapidly developing flows cannot be captured, see Wallin and Johansson (2000). Furthermore, the Mellor and Yamada approach leads to the following balance of terms, 2 3 δ i j (P + G ) which is fully linear in the turbulent fluxes (assuming a linear model for Π i j is used). Therefore, Mellor and Yamada level-3 models do not contain a non-linear equation for the production-to-dissipation ratio (i.e. the factor N in (32)). This is a major difference in complexity between these types of models and EARS models.