Explicit Algebraic Reynolds-stress Modelling of a Convective Atmospheric Boundary Layer Including Counter-Gradient Fluxes

In a recent study (Želi et al. in Bound Layer Meteorol 176:229–249, 2020), we have shown that the explicit algebraic Reynolds-stress (EARS) model, implemented in a single-column context, is able to capture the main features of a stable atmospheric boundary layer (ABL) for a range of stratification levels. We here extend the previous study and show that the same formulation and calibration of the EARS model also can be applied to a dry convective ABL. Five different simulations with moderate convective intensities are studied by prescribing surface heat flux and geostrophic forcing. The results of the EARS model are compared to large-eddy simulations of Salesky and Anderson (J Fluid Mech 856:135–168, 2018). It is shown that the EARS model performs well and is able to capture the counter-gradient heat flux in the upper part of the ABL due to the presence of the non-gradient term in the relation for vertical turbulent heat flux. The model predicts the full Reynolds-stress tensor and heat-flux vector and allows us to compare other important aspects of a convective ABL such as the profiles of vertical momentum variance. Together with the previous studies, we show that the EARS model is able to predict the essential features of the ABL. It also shows that the EARS model with the same model formulation and coefficients is applicable over a wide range of stable and moderately unstable stratifications.


Introduction
Turbulence parametrization schemes that are used for physical parametrization of the atmospheric boundary layer (ABL) are usually tailored with disparate descriptions for stable and unstable stratification, or sometimes even for different levels of stratification (see He et al. 2019). An attractive alternative is a turbulence model with sufficiently general formulation B Velibor Želi velibor@mech.kth.se 1 Department of Engineering Mechanics, FLOW Centre, KTH Royal Institute of Technology, 10044 Stockholm, Sweden that enables accurate predictions of both the stably stratified and convective ABL. One challenge with the latter is the existence of counter-gradient heat fluxes in the upper part of the convective ABL. These fluxes are usually modelled with additional non-gradient terms (see Holtslag and Moeng 1991;Siebesma et al. 2007) in addition to an eddy-diffusivity type of terms. Lazeroms et al. (2013) derived the so-called explicit algebraic Reynolds-stress (EARS) model for the general case of stratified flows and applied it in the context of the ABL (Lazeroms et al. 2015(Lazeroms et al. , 2016. The EARS model is based on transport equations for Reynolds-stress tensor and heat-flux vector, i.e., turbulence momentum and heat flux, respectively. The full transport equations are simplified by applying the weak-equilibrium assumption (Rodi 1972(Rodi , 1976) into a set of algebraic equations where the unknown correlations require modelling. To close these equations, three additional transport equations for turbulence kinetic energy (TKE), dissipation rate of TKE, and half the potential temperature variance are solved. Recently, atmospheric turbulence models have been developed that share similarities with the EARS model in the sense that they also include additional transport equations for turbulence scalar quantities (aside from the TKE). Zilitinkevich et al. (2013) added a transport equation for turbulence potential energy and Machulskaya and Mironov (2020) added a transport equation for potential-temperature variance. Both scalar quantities are directly related to half the potential temperature variance. The coefficients in the EARS model are not calibrated for a particular case study, but have instead been calibrated by use of a set of generic turbulent flows (see Želi et al. 2020). The complete description and detailed derivation of the EARS model, together with the boundary-condition treatment, are given in Želi et al. (2019).
In Želi et al. (2020) it was shown that the EARS model gives a good prediction of first-order statistics and other details in a stably stratified ABL. In the present study, we extend this work by applying the same EARS model formulation in an idealized case of a convective ABL and comparing the results with the large-eddy simulation (LES), results in Salesky and Anderson (2018). The present results illustrate the predictive capability of the EARS model for an unstably stratified ABL, and together with the previous results, illustrate the applicability of the model for the ABL under different types of stability conditions. This can be regarded as a step towards implementation in an operational ABL model.
The process of turbulent heat transport in the vertical direction governs the amount of TKE and drives the convective ABL. The complete transport equation for a horizontally uniform, dry ABL without subsidence takes the following form where the terms on the right-hand side (r.h.s.) are (from left to right) mean-field production, buoyancy production, pressure scrambling, dissipation, and diffusive transport; where w denotes the velocity fluctuation along the vertical direction z (x 3 ), and θ the mean and fluctuating potential temperature, g the acceleration due to gravity, and T 0 the reference temperature. We note that the production terms contain both a term proportional to the mean potential temperature gradient and one that is independent of this gradient, and the relative magnitude of these depends on the type and level of stratification. This can also be said about the model expression for the pressure scrambling term. Correspondingly, turbulence models for the convective ABL reflect this behaviour and split vertical heat flux into gradient and non-gradient parts (see Holtslag and Moeng 1991;Siebesma et al. 2007). A critical feature of such models is that the direction of heat flux is not only governed by potential-temperature gradient and can predict counter-gradient fluxes. In the EARS model, turbulent heat flux in the vertical direction has the following representation where K denotes is the TKE, is the dissipation rate of TKE, K θ is the half-temperature variance (i.e., θ 2 /2), f h and f are the so-called core functions, which are key elements of the EARS model. The core functions are non-constant and dimensionless coefficients that follow directly from the model's derivation and are functions of the state of flow in the ABL. They are implicitly dependent on the complete set of transport equations, not only the vertical flux equation given in (1). For more information about the core functions see Appendix 1 in Želi et al. (2019). The terms in (2) represent the gradient and non-gradient parts of the turbulent heat flux, respectively. Similar non-gradient terms exist in the turbulent heat, as well as momentum, fluxes. Hence, the EARS model predicts non-zero horizontal turbulent fluxes of heat and momentum even in the single-column formulation of the EARS model (see Želi et al. 2020). Although the expression for vertical heat flux in the EARS model has the same form and function as proposed in Holtslag and Moeng (1991), there are some significant differences between these models. Namely, in the EARS model the mean-field and buoyancy production terms are expressed in closed form, without the need for modelling. This ensures a high level of fidelity of the modelling in general and non-gradient part of the heat flux in particular. Another consequence and prominent difference between the EARS model and the models such as those in Holtslag and Moeng (1991) and Hong et al. (2006) is that the non-gradient term does not depend explicitly on the height of the convective ABL. In the EARS model this property is an integrated part of the resulting solution that is carried in the core functions. In contrast, the height of the convective ABL in the models such as Holtslag and Moeng (1991) and Hong et al. (2006) is a critical parameter that is usually modelled in an ad hoc manner.

Case Description
Simulation results for convective ABLs are produced using the EARS model implemented in the context of the single-column model from Želi et al. (2020). The detailed model formulation, given in symbolic form, is automatically discretized by the use of symbolic manipulation software in order to minimize the risk of human-induced errors. The resulting machine-generated Fortran code is then used for obtaining the numerical solution of the time-dependent one-dimensional problem by second-order central difference in space and Crank-Nicolson time discretization.
Reference data and simulation set-up are taken from the LES study of Salesky and Anderson (2018). A dry convective ABL is driven by a time-constant kinematic heat flux at the surface q (m K s −1 ) and geostrophic wind velocity V g = (U g , V g ) (m s −1 ) at the top of the boundary layer oriented along the x-axis, i.e., V g = 0. Five simulations with the following combinations of the forcing parameters are studied (q, U g ) = [(0.07, 15), (0.10, 15), (0.24, 15), (0.24, 11), (0. 24,9)]. The wind forcing and surface forcing directly influence the production of TKE in the ABL due to the buoyancy and shear effects. The choice of the forcing parameters corresponds to simulations of ABL conditions ranging from weakly to moderately convective in which three simulations share the same geostrophic wind velocity and vary surface forcing and three simulation have the same surface forcing while geostrophic wind velocity varies.
The numerical domain in this study has a stretched grid with 106 grid points. The height of the domain is H = 2000 m. Near the surface and at the height of z = 1100 (entrainment zone) the grid size z is refined in order to correctly represent rapid spatial change in wind speed and potential-temperature profiles, z ≈ 2 m and z ≈ 7 m, respectively. Between these regions the grid is smoothly stretched up to z ≈ 30 m. This resolution is fine enough to obtain grid independence. In more general implementations it can be of interest to use adaptive grid techniques, such as that of van Hooft et al. (2018a), which has shown promising results for both quasi-steady and transitional ABL cases (see van Hooft et al. 2018bvan Hooft et al. , 2019. Other parameters that define the simulations are g = 9.81 m s −2 , Coriolis parameter f = 10 −4 s −1 , reference temperature T 0 = 300 K, and roughness lengths for momentum and temperature z 0 = 0.1 m. As in Salesky and Anderson (2018), the mean potential-temperature profile was initialized using the three-layer profile where Γ 1 = 0.08 K m −1 and Γ 2 = 0.003 K m −1 determine the gradients of the potential temperature. The initial wind is equal to the geostrophic throughout the numerical domain. At the time of the initialization, the ABL is stably stratified. The initial TKE profile is In the present study the initial profile is prescribed through specifying the turbulence time scale τ 0 = K / = 1 s at z 0 that increases linearly to τ 0 = 550 s at 250 m, above which it is constant. This means 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. 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 .

Results
Simulations of the idealized the convective ABL layer with varying forcing parameters q and U g , described in Sect. 2, are compared to the LES of Salesky and Anderson (2018). The overview of bulk properties from the LES and present simulations are given in Table 1, where the boundary-layer depth z i is found at the location of the minimum of the vertical turbulent heat flux. Results of the simulations with the combination of forcing parameters (0.10 K m s −1 , 15 m s −1 ) and (0.24 K m s −1 , 9 m s −1 ) are not shown in the following figures because the results of these simulations can be interpolated from the remaining three simulations. The convective ABL develops in time with the constant boundary forcing. The EARS model and LES model with the same initial conditions have different initial transients. To make a clean comparison between the results, we ran the simulations of the EARS model until the thermal energy in the mixed layer of ABL become the same as in the LES model. In practice, this means that the value of potential temperature in the mixed layer of the EARS model and LES model are identical. In the LES model the statistics are averaged in time between the fourth and fifth simulation hour, while in the model simulations the Table 1 Bulk properties including kinematic heat flux (q), geostrophic wind speed (U g ), Obukhov length (L), friction velocity (u * ), Deardorff convective velocity scale (w * ), and boundary-layer height (z i ) corresponding to the simulation with the EARS model and LES model of Salesky and Anderson (2018) model instantaneous results are taken at 5.45 h so that the amount of heat related to the potential temperature in the mixed layer of the ABL is the same. Figure 1 shows vertical profiles of turbulent heat flux scaled with the value of the heat flux at the surface. Heat flux in the LES model varies linearly with height up to the entrainment zone where it changes sign and gives a negative peak that agrees reasonably well with the value reported by Conzemius and Fedorovich (2006). The EARS model captures the linear behaviour and the existence of a negative peak in the entrainment zone, although the model gives a significantly smaller amplitude of the peak. Note that the linearity of heat flux in the mixed layer is a consequence of the mean equation for the potential temperature (and is obtained with any turbulence model). Thus, the performance of the EARS model will be evaluated according to the correctness of the potential-temperature profiles.
Near the surface the vertical heat flux is equal to the value of the surface heat flux u * θ * , where u * is the friction velocity and θ * is the characteristic temperature scale. This shows that the methodology for imposing the boundary conditions presented in Želi et al. (2019) gives consistent behaviour of the model in the near-surface region.
Although heat-flux linearity is a mere direct result of the quasi-steadiness of the system, Fig. 2 shows that the contribution of the gradient and non-gradient terms are non-trivial, i.e., far from linear themselves. Those shapes reflect the actual physics being at play (see Eq. 2). Figure 2 shows both parts of the vertical heat flux scaled in the same way as in Fig. 1. The gradient and non-gradient parts of the model individually strongly depart from a linear variation. The former dominates near the surface while the latter dominates in the major part of the layer up to the entrainment zone. The LES study from Zhou et al. (2018) suggests that both parts have a positive value for z < 0.8z i . However, both the local and non-local fluxes are implicitly part of the LES solution, making the separation of the two parts non-trivial. Comparison of the non-gradient part can only be made at the particular location where the potential-temperature gradient in the LES model is zero. This result is marked in Fig. 2 with a black point and agrees with the result of the EARS model. The results in the EARS model are qualitatively similar to the ones reported in Holtslag and Moeng (1991) and Hong et al. (2006) with a maximum of the non-gradient part found at z ≈ 0.2z i .
In the model simulations with stronger surface forcing the mixed layer constitutes a larger part of the ABL. The gradient part of (2) is smaller in the regions where the potential tem-  Salesky and Anderson (2018). Large-eddy simulation data are not available for the upper part of the domain and therefore is not plotted perature gradient approaches zero and therefore the non-gradient term becomes dominant in the larger part of the ABL. The sign (direction) of the heat flux is not only governed by the gradient of the potential temperature but also by the magnitude of the non-gradient term. The gradient term shown in Fig. 2 changes sign around z ≈ 400 − 500 m, while the total heat flux is positive up to the entrainment zone (z ≈ 1000 m). Figure 3 shows that the results of the EARS model for potential temperature (solid lines) in a large part of the mixed layer are reasonably close to the LES (dashed lines), ensuring that the amount of heat related to the potential temperature in the convective ABL is comparable. The potential-temperature profile in the entrainment zone in the LES has a weaker gradient, and the boundary layer is somewhat higher (see Table 1). This could partly be associated with the averaging used in the LES model where profiles are smeared out as the top of the the convective ABL develops over time. The gradient of potential temperature in the results of the EARS model changes sign at z ≈ 400−500 m at the same point where the non-gradient part of the heat flux from Fig. 2 changes sign. This is similar with the values in the LES. Figure 3 also shows the results of the EARS model in which only the gradient part of the vertical heat flux is active and the non-gradient part is set to zero (dotted line). The vertical heat flux retains the linear dependence with respect to height and deviates slightly in the entrainment zone from the results of the original EARS model with the complete vertical heat flux (not shown). More profound changes are observed in the profile of potential temperature. In the simulations without modelling of the non-gradient term in the heat flux, the gradient of the potential temperature shows a physically incorrect behaviour with the same sign as the vertical heat flux throughout the domain. This agrees with findings of Siebesma et al. (2007) and Wang et al. (2016), who reported that different local and nonlocal turbulence models have quantitatively similar vertical heat-flux profiles but significantly different potential-temperature profiles. As in Siebesma et al. (2007) and Wang et al. (2016),  Salesky and Anderson (2018) this study indicates that the turbulence models that carry only the gradient term are not able to accurately predict the sign of a gradient of potential temperature and turbulent heat flux at the same time.  Figure 4 shows that the horizontal wind-speed profiles are well predicted by the EARS model. The agreement is better in the simulations with a larger Obukhov length, where the surface heat flux is smaller and the wind shear is larger. Similar to the potential-temperature profile, the horizontal wind speed in the LES shows a smooth gradient in the entrainment zone. The high wind shear close to the top and the bottom of the entrainment zone generates Kelvin-Helmholtz waves that mix and enhance transport in the vertical direction (Rayment and Readings 1974) The EARS model, however, cannot predict this instability-related contribution.
The EARS model gives a solution for the full Reynolds-stress tensor and heat-flux vector. As an example, Fig. 5 shows the variance of vertical velocity fluctuation ww for the EARS and LES models. The EARS model predicts a significantly larger ww than the LES in all cases. The large value of ww is likely related to imperfect modelling of the dissipation rate of TKE, which leads to high TKE and half the potential-temperature variance. A similar observation is reported in Želi et al. (2020) for the stably stratified case. In the model simulations with stronger surface forcing the agreement is better. Apart from the intensity amplitude, the results show qualitatively the same behaviour as in the LES model. Due to internal waves, the results from the LES model have a non-zero value for the variance of vertical velocity fluctuation at the top of the boundary layer, which cannot be captured by the EARS model.  (2018) with a well-resolved LES. The comparison shows that the EARS model is able to predict bulk properties and first-order statistics with fairly good accuracy. The gradient of potential temperature in the EARS model is quantitatively similar to the LES model in the sense that it can be both positive and negative in the mixed layer. The horizontal wind speed in the mixed layer is close to the LES data. The model is local and unable to replicate waves at the top of the the convective ABL, i.e., Kelvin-Helmholtz and gravity waves, therefore the boundary-layer height predicted by the model is smaller than in the LES model. This also affects the overall behaviour of the results in the entrainment zone. Particularly, the normalized turbulent heat flux in the entrainment zone is weaker than what other studies report. Girimaji and Balachandar (1998) reported that the EARS model has a stable solution for the purely convective case. Our preliminary results also show that the EARS model gives realistic results in strongly convective conditions. However, detailed analysis for the case of strong convection should be carried out in the future, together with reimplementation of boundary conditions which are more suited for this propose. The EARS model contains a non-gradient term in the expression for the vertical turbulent heat flux. This term results from the weak-equilibrium approximation and modelling of the terms in the full Reynolds-stress turbulence model, and is not an explicit function of the height of the ABL. This makes the EARS model less ad hoc than commonly used models. Also the boundary-layer height is computed as an integral part of the model solution carried in the core functions. The non-gradient term of the heat flux has an important role and strongly influences the predicted temperature profiles. It allows the gradient of potential temperature to be positive (indicating stable stratification) under conditions where the turbulent heat flux is positive and heat is being transported from a level with lower potential temperature to a level with higher, the so-called counter-gradient heat flux. Commonly, atmospheric turbulence models are designed for only one region of stability and have to use different turbulence parametrizations in order to handle cases with transitioning atmospheric stability. The EARS model is a generalized model with the same formulation and calibration for stable, neutrally stable, and convective ABLs (see applications to diurnal cycle in citealtvzeli2019consistent) making it more appealing for a wider range of applications.