A new turbulence scheme with two prognostic energies is presented. The scheme is an extension of a turbulence kinetic energy (TKE) scheme following the ideas of Zilitinkevich et al. but valid for the whole stability range and including the influence of moisture. The second turbulence prognostic energy is used only for a modification of the stability parameter. Thus, the scheme is downgradient, and the turbulent fluxes are proportional to the local gradients of the diffused variables. However, the stability parameter and consequently the turbulent exchange coefficients are not strictly local anymore and have a prognostic character. The authors believe that these characteristics enable the scheme to model both turbulence and clouds in the planetary boundary layer. The two-energy scheme was tested in three idealized single-column model (SCM) simulations, two in the convective boundary layer and one in the stable boundary layer. Overall, the scheme performs better than the standard TKE schemes. Compared to the TKE schemes, the two-energy scheme shows a more continuous behavior in time and space and mixes deeper in accordance with the LES results. A drawback of the scheme is that the modeled thermals tend to be too intense and too infrequent. This is due to the particular cutoff formulation of the chosen length-scale parameterization. Long-term three-dimensional global simulations show that the turbulence scheme behaves reasonable well in a full atmospheric model. In agreement with the SCM simulations, the scheme tends to overestimate cloud cover, especially at low levels.
The unified parameterization of turbulence and clouds in the planetary boundary layer (PBL) is one of the challenges in current numerical weather prediction (NWP) and general circulation models (GCMs; Siebesma et al. 2007; Golaz et al. 2002). The unified approach should ideally enable a better representation of interactions and transitions and eliminate compensating errors that arise from the separate modeling of turbulence and clouds.
One group of such unified parameterizations is based on the extension of the local downgradient parameterization with a nonlocal component that accounts for the transport caused by thermals. The nonlocal component can be introduced by a countergradient term (Deardorff 1972; Holtslag and Moeng 1991), which leads to upward turbulent heat transport in the presence of a slightly stable temperature gradient. A more sophisticated representation of thermals via mass flux is used in the eddy-diffusivity and mass-flux (EDMF) schemes (Siebesma et al. 2007; Neggers et al. 2009; Neggers 2009; Sušelj et al. 2012).
The other family of parameterization relies on a higher-order turbulent closure (HOC), where the higher-order moments enable the representation of nonlocal and countergradient transport. These parameterizations range from full third-order closure (e.g., Bougeault 1981) and schemes with prognostic equations for all second-order moments and a third-order moment for vertical velocity (e.g., Golaz et al. 2002) to closures for turbulence kinetic energy (TKE) together with prognostic scalar second-order moments (e.g., Nakanishi and Niino 2009) or scalar variances (e.g., Machulskaya and Mironov 2013).
The schemes with only prognostic TKE are able to simulate stratocumulus clouds (Duynkerke and Driedonks 1987; Bechtold et al. 1992) and, with proper treatment of cloudiness distribution and an advanced length-scale formulation, even trade wind cumulus (Bechtold et al. 1995; Bogenschutz and Krueger 2013).
In the present paper, we propose a new scheme that is based on prognostic equations for two turbulence energies. The additional prognostic turbulence energy is a compressed equivalent of the two scalar variances. In the closure hierarchy, this scheme is one step above schemes with only prognostic TKE.
Our scheme builds on the TKE scheme by Bašták Ďurán et al. (2014), which is used in the ALARO physical package of the ALADIN NWP model (Termonia et al. 2018; Wang et al. 2018) and is implemented in the Integrated Forecast System (IFS) model developed at the European Centre for Medium-Range Weather Forecasts (ECMWF). We added a second prognostic turbulence energy by following the work of Zilitinkevich et al. (2007, 2008, 2013), who proposed the usage of two turbulence energies in the stable boundary layer. We extend their derivation to all atmospheric stabilities and include the influence of moisture.
The additional information from the second turbulence energy is used solely for the computation of the stability parameter, which thus obtains a prognostic character and an independence from local gradients. The scheme uses the downgradient approach, where the turbulent fluxes are proportional to the local gradients but the turbulent exchange coefficients are not strictly local. This enables transport across locally stable layers and deeper mixing, which eliminates the major problems associated with traditional downgradient schemes. A similar behavior was achieved in Bretherton and Park (2009) by extending their convective layers into stable regions and merging them if they overlap and by using averaged stability functions across each of these layers. Still, our scheme is unable to parameterize countergradient fluxes and is therefore less accurate than higher-order closure schemes. On the other hand, the two-energy scheme is simpler in design, and because the scheme is easier to implement, it is computationally cheaper and numerically stable for longer time steps.
The influence of moisture and phase changes of water on vertical transport plays an important role in the parameterization of PBL clouds. In HOC and TKE schemes, this effect is realized through a modification of the buoyancy terms or the Brunt–Väisälä frequency. The cloud fraction or the buoyancy terms can be computed by assuming a Gaussian distribution of temperature and moisture (Sommeria and Deardorff 1977). However, the distributions are skewed when shallow convection is present (Lewellen and Lewellen 2004) and need to be adjusted accordingly (Bechtold et al. 1995). A more advanced approach is to use a trivariate (temperature, moisture, and vertical velocity) probability density function (PDF; Golaz et al. 2002; Larson et al. 2012; Bogenschutz and Krueger 2013) for the direct computation of the buoyancy terms.
In our scheme, the buoyancy terms are computed as a cloud-dependent linear combination of covariances for temperature and moisture. The influence can then be represented via a modification of the stability parameter used in the stability dependency functions. We use the Gaussian formulation of Sommeria and Deardorff (1977) to determine the cloud fraction and compute the buoyancy terms according to Marquet and Geleyn (2013) with a skewness-dependent heuristic adjustment.
This paper is organized as follows. Section 2 gives a general overview of our TKE scheme for the dry and the moist case. In section 3, we describe the addition of a second turbulence energy and the corresponding stability parameter. Details about the turbulence energy solver are presented in section 4. The scheme is evaluated in section 5 in the single-column model (SCM) setup of the OpenIFS model and in section 6 in a long-term three-dimensional simulation of the IFS model. We finish with a summary and conclusions in section 7.
2. The TKE turbulence scheme
a. The dry TKE turbulence scheme
To demonstrate the impact of moisture on the turbulence scheme, we will first present the TKE scheme for the dry case, following closely Bašták Ďurán et al. (2014).
The turbulence scheme is based on a prognostic equation for TKE :
where u, υ, and w are the wind components; I is the shear production term; is the buoyancy production/destruction term in the dry case; θ is the potential temperature; g is the gravitational acceleration; z is the height; t is the time; is the dissipation time scale for the TKE; is the turbulent diffusion coefficient; L is the length scale for molecular dissipation; is the share of vertical component in TKE computed according to Bašták Ďurán et al. (2014); and and are closure constants.
The starting point for the derivation of the expressions for the vertical turbulent fluxes is the Cheng et al. (2002) scheme. It is a Mellor–Yamada-type scheme with improved relationships for the pressure correlation terms. There, the turbulent fluxes are parameterized using downgradient closure:
where and are stability functions for momentum and heat, is the dry Brunt–Väisälä frequency (BVF), S is the wind shear, is the inverse of the Prandtl number at neutrality, and and (see Cheng et al. 2002) are closure constants.
The Cheng et al. (2002) scheme suppresses turbulence in stratification stronger than the critical gradient Richardson number. With the aim of removing this characteristics of the Cheng et al. (2002) scheme, Bašták Ďurán et al. (2014) assumed that the dissipation time scale for the heat flux is stability dependent and used equilibrium conditions to express the stability functions and as functions of flux Richardson number, :
where P and R are either constants or functions of , depending on the complexity of the chosen scheme. The resulting and stability functions are valid for whole stability range, and they include the influence of anisotropy (see Bašták Ďurán et al. 2014).
The exchange coefficients for momentum and heat are computed from an independent prognostic TKE value in the Bašták Ďurán et al. (2014) scheme, but the stability dependency functions are derived by assuming steady state in the TKE equation:
b. Length scale
The length scale L in the above expressions may be chosen independently from the stability functions and TKE. For our scheme, we will use a combination of two types of length scales.
The first one is a similarity law mixing length-scale after Cedilnik (2005):
where κ is the von Kármán constant; is the PBL height; and , , , and are constants. The length scale is then computed from according to Redelsperger et al. (2001) in order to match the similarity laws for the neutral surface layer:
The second one is derived after Bougeault and André (1986) and Bougeault and Lacarrere (1989), where is given by the balance between TKE and local stability. We use a modified combination of their upward free path and the downward free path :
Note that for consistency reasons, it is important that the expression for buoyant eddies in and shares the same discretization with the production term .
To ensure similarity law behavior in the lower levels via and limit the turbulent mixing according to in the upper levels, the final length scale is taken as the minimum of the above two formulations:
c. Density fluctuation in moist case
Compared to the dry case, two instead of one (potential temperature) conservative thermodynamic variables are required to uniquely define an air parcel in the moist case (Betts 1973; Marquet 2011). This adds five more Reynolds-averaged Navier–Stokes (RANS) equations for second-order moments to the system (see, e.g., Stull 1988). The presence of water and its phase changes influence the density fluctuations and thus modify the buoyancy. We will accommodate these extensions and, analogous to the dry case, derive the relationships for turbulent vertical fluxes in the moist case.
We have chosen to use static energy and total specific water content as our conservative variables and express the density fluctuation following the derivation of moist BVF in Marquet and Geleyn (2013):
where C is the cloud fraction; is the skewness-influenced turbulence cloud fraction (to be defined below); , , , and are the specific heats for dry air, water vapor, liquid water, and ice; and are the latent heats of vaporization and sublimation; T is the temperature; and are the specific contents for liquid and solid water; is the dry air gas constant; is the water vapor gas constant; is the specific content for water vapor; is the specific content for dry air; is the latent heat of vaporization (C) or sublimation (C); is the mixing ratio for saturating water vapor for vaporization (C) or sublimation (C); is the partial saturation pressure over liquid (C) or solid water (C); ρ is air density; and is the reference air density.
Compared to Marquet and Geleyn (2013), the dimensionless quantity now depends on instead of the mixing ratio for water vapor . This modification makes dependent only on the local variables T and p, causing to be invariant to changes of the atmospheric composition. Such a formulation does not change values of and at the extreme boundary conditions for and .
The cloud fractions C and are critical parameters in Eqs. (21) and (22). The determination of these quantities is beyond the scope of this paper. We will adopt here a practical approach and base our C computation on the well-established parameterization of nonprecipitating clouds after Sommeria and Deardorff (1977). Under the assumption that and the liquid water potential temperature are Gaussian distributed, the cloud fraction C can be computed as follows:
To represent the influence of skewness on the density fluctuations (see, e.g., Lewellen and Lewellen 2004), which is essential for the parameterization of shallow cumulus clouds in a downgradient scheme (see Bechtold et al. 1995), a skewness-dependent cloud fraction was used in the linear-type term in Eq. (22). The dependence of on skewness-linked parameter was heuristically determined using LES data (courtesy of D. Lewellen):
where can be interpreted in terms of buoyancy flux [see Eq. (35) below] as the effective cloud fraction for neutral stratification and shows a significant correlation with the skewness of vertical velocity.
d. The moist TKE turbulence scheme
Now that we have an expression for the density fluctuations [see Eq. (19)], we can directly express the relationship for the buoyancy terms (covariances with buoyancy) in the moist case:
where ψ stands for any variable from u, υ, w, , and .
As in the dry case (Bašták Ďurán et al. 2014), we will use the Cheng et al. (2002) paper as the baseline, but we first generalize the RANS equations for the moist case. To ensure convergence of the moist version of the equations to the dry one, we will not modify the values of the closure coefficients. We express the buoyancy terms via Eq. (35) and get the following set of diagnostic equations for the second-order moments:
where is the Kronecker delta.
With the help of the symbolic computation system wxMaxima, the system of Eqs. (36)–(49) can be solved. The resulting stability functions and have the same functional form as the stability functions in the dry case [Eqs. (8) and (9)], only that in the moist case, the functions depend on instead of .
As a consequence of the above, the derivation of the one-dimensional turbulent parameterization in Bašták Ďurán et al. (2014) can be applied also to moist case, where the flux Richardson number is substituted with its moist counterpart. The -dependent stability function and have the same shape and so also have the same properties as in the dry case [Eqs. (11) and (12)]. We will use these stability functions in our scheme for both the dry and the moist cases.
3. The two-energy turbulence scheme
Only a single stability parameter, namely, flux Richardson number , is necessary for the closure in the above TKE turbulence scheme. It can be computed in one of the following three ways:
We will investigate a third option and express as a function of prognostic second-order moments. This will result in the extension of the scheme with a second turbulence energy.
By definition, is closely related to the shear I and the buoyancy terms in the TKE equation:
In the spirit of our closure, where stability functions are derived for equilibrium conditions but TKE is prognostic, we can relate the two terms to TKE if we neglect the time tendency, the advection, and the vertical turbulent transport terms in the TKE equation [Eq. (55)]:
An additional quantity is required to obtain an independent expression for the shear and the buoyancy terms. As in Zilitinkevich et al. (2013) for the dry case, the buoyancy term in the TKE equation can be eliminated using the prognostic heat and moisture variance equations:
where TT represents a third-order moment term.
To do so, we multiply the variance equations by and , respectively, and combine them into
where is the turbulence energy, whose only source is the buoyancy destruction term in the TKE equation. The only sinks are dissipation and the buoyant generation of TKE. Adding the equation [Eq. (60)] to the TKE equation [Eq. (55)] results in
where is sum of TKE and , is the dissipation time scale for , and is the turbulent diffusion coefficient for , where we assume that we can parameterize the turbulent transport term in the equation using a downgradient closure. For the sake of simplicity, we will assume that = .
Neglecting the time tendency, the advection and the vertical turbulent transport terms in the TKE equation [Eq. (55)] and the equation [Eq. (61)], we can establish a link between the two turbulence energies and and the flux Richardson number :
and also express the dissipation time scale for in terms of :
The relationship Eq. (64) enables the computation of from the two prognostic turbulence energies.
The two-energy scheme can be easily implemented into the TKE scheme, where the prognostic equation for TKE is already implemented. The second prognostic variable can be either or , as they are linked via Eq. (62). We can assume that is always positive, since the shear term I under the downgradient assumption is always positive and the dissipation time-scale , computed according to Eq. (65), will be also always positive (, ). Thanks to the nonnegative property of , we can use the TKE solver ( is always nonnegative) also for the computation of . This simplified implementation favors the choice of as the second prognostic turbulence energy.
Please note that in the dry case in stable stratification, would be equal to the turbulence total energy, and would be equal to the turbulence potential energy (Zilitinkevich et al. 2013). This is not the case in the moist case because the turbulence potential energy is then equal to
which differs from .
In comparison to the downgradient prognostic model of Zilitinkevich et al. (2013), no assumption about stratification has been made in our scheme. This means our scheme can be used for both stable and unstable stratification.
4. Turbulence energy solver
An essential part of the turbulence scheme is the solver for the turbulence energies. The details of the TKE solver will be presented in the following subsection.
a. TKE solver
To ensure stability of vertical diffusion schemes, it is common practice to evaluate the turbulent exchange coefficients on half levels (i.e., between the levels where model variables are defined). Direct dependence of the exchange coefficients on TKE [see Eqs. (8) and (9)] implies that TKE would be best located on the half levels. Such a discretization however conflicts with the practical requirement to advect TKE consistently with the other prognostic variables located on the full levels. This issue can be solved either by putting TKE on the half levels and neglecting the advection of TKE or by using a full-level TKE and interpolating it to the half levels for the computation of the exchange coefficients. The former option works fine for scales coarser than 10 km and/or time steps longer than 5 min, where the advection of TKE contributes typically less than 1% to the total tendency. The latter option adds some implicit smoothing, which might not be desirable. Here, TKE is kept on full levels in order to ensure problem-free advection of TKE.
The numerical stability of the TKE solver is achieved by a careful discretization and by adopting methods that deal with stiffness typically occurring in the vertical diffusion schemes (Kalnay and Kanamitsu 1988; Bénard et al. 2000). To mitigate the stiffness problem, the prognostic equation for TKE is written as
with ADV representing the advection operator. The last term on the right-hand side can be regarded as representation of the TKE source and sink terms. Contrary to the standard formulation as the sum of the shear, buoyancy, and TKE dissipation term, it is expressed as a relaxation toward the static equilibrium value of TKE for which the sum of the above three terms is zero. By simple algebraic manipulation, it can be shown that
The simplified prognostic TKE equation [Eq. (67)] can be seen as a prognostic equation having only two sources: autodiffusion and relaxation. To prevent any potential stiffness issues caused by the source terms, an overimplicit treatment is adopted [following the conclusions of Kalnay and Kanamitsu (1988)]. The relaxation term is being discretized with the constant overimplicit parameter equal to 1.5. For the autodiffusion term, the overimplicit parameter β is set to the minimal value that ensures stability according to the method proposed by Bénard et al. (2000).
The discretization of the above TKE equation [Eq. (67)] without the advection term for a level i then becomes
with representing the thickness of the layer i, and δ and α are weights for interpolation from half levels to full levels. The index denotes half levels between the full levels i and . The levels are counted from the top to bottom. The is defined as , and the half-level value of the full-level quantity is obtained by standard averaging: .
The boundary conditions for the top and bottom are obtained from the flux boundary condition imposing . While the stationary TKE at the upper bound is set to zero (), at the lower bound, it is set to the value used for the surface layer:
with the friction velocity (see Redelsperger et al. 2001).
The numerical stability of the solver is further ensured by making the matrix of the solver [see Eq. (69)] diagonally dominant. This is achieved by controlling the relationship between the diffusion and the dissipation [see Eqs. (3) and (5)] in terms of
If the matrix is not diagonally dominant, the length scale L is locally adequately increased.
b. Modification for second turbulence energy
As mentioned in section 3, we have chosen as the second prognostic turbulence energy in our scheme because it is by definition always nonnegative like TKE. So we can use the TKE solver also for when we modify (i) the source terms by excluding the buoyancy term, (ii) the dissipation term by switching to the dissipation time-scale , and (iii) the turbulent transport term by using instead of .
5. Idealized single-column model simulations
To evaluate our turbulence scheme, we implemented it into the ECMWF OpenIFS SCM.
The OpenIFS (Carver and Váňa 2017) model is the portable version of the IFS at ECMWF. The model is based on the IFS model cy40r1 (ECMWF 2014). All SCM simulations are performed with a 900-s time step using the two-time-level semi-Lagrangian scheme. The SCM runs with 91 atmospheric vertical levels (ECMWF 2017d), having 17 levels in the lowest 2 km.
We tested the SCM for three idealized cases, all developed by the Global Energy and Water Cycle Experiment (GEWEX): the continental cumulus case based on measurements over the Atmospheric Radiation Measurement (ARM) program Cloud and Radiation Testbed (CART) site in Oklahoma (Brown et al. 2002; Lenderink et al. 2004), the trade wind cumulus case based on observations from the Barbados Oceanographic and Meteorological Experiment (BOMEX; Siebesma et al. 2003), and the third test case is based on the GEWEX Atmospheric Boundary Layer Study (GABLS) project (Holtslag 2006; Beare et al. 2006; Cuxart et al. 2006).
This particular set of cases was chosen to analyze the turbulence scheme without interaction with other physical parameterizations. Although this is a limitation, the situations tested do include stable stratification (GABLS) and unstable stratification. The latter are two shallow convection cases: the quasi-steady state over ocean (BOMEX) and the diurnally driven case over land (ARM). Especially the capability of our scheme to represent moisture effects (i.e., without the need for a separate shallow convection scheme) is tested with the latter two cases.
Thus, the OpenIFS runs include only the turbulence scheme, and the remaining IFS physics are turned off. As a reference, we also run OpenIFS with the EDMF scheme and shallow convection parameterization (SCP; Bechtold et al. 2014).
The SCM results for the above cases were compared with results from the large-eddy simulation (LES) model MicroHH (van Heerwaarden et al. 2017). In the ARM and BOMEX simulations, MicroHH runs with moist thermodynamics and in the GABLS simulation with dry thermodynamics. Microphysical and precipitation processes are turned off for all cases. For the ARM case, the domain size is 6.4 km × 6.4 km in the horizontal and 4.4 km in the vertical with horizontal and vertical grid spacings of 100 and 34.375 m, respectively. For the BOMEX case, the domain size is 6.4 km × 6.4 km in the horizontal and 3.0 km in the vertical with horizontal and vertical grid spacings of 100 and 46.875 m, respectively. For the GABLS case, the domain size is 400 m × 400 m in the horizontal and 400 m in the vertical with horizontal and vertical grid spacings of 3.125 and 3.125 m, respectively. All the MicroHH simulations reproduced the results of the respective LES intercomparison studies.
The LES statistics were averaged over the whole horizontal domain in space and over one hour in time. All figures were generated with the Python matplotlib package (Hunter 2007).
a. Continental shallow cumulus
To outline the prognostic properties of the two-energy turbulence scheme, we start the analyses with the ARM case, which has a pronounced diurnal cycle with onset and decay of shallow convection.
Figure 1 shows the vertical profiles of selected variables for the LES and the SCM in four different configurations. Three configurations of the turbulence scheme (see section 2d) differ in the formulation of the flux Richardson number stability parameter: the TKE model based on gradients [TKEM-GR; see Eq. (51)], the TKE model from fluxes [TKEM-FL; see Eq. (53)], and the TKE model from the two turbulence energies [TKEM-TE; see Eq. (64)]. The fourth configuration uses the IFS implementation of the EDMF and a shallow convection parameterization, the IFS-EDMF+SCP. Please note that the internal cloud fractions C [Eq. (27)] and [Eq. (32)] are not available for the IFS-EDMF+SCP configuration and differ from the IFS cloud-cover diagnostics. The LES is depicted for three different sampling criteria applied to the LES data: ALL refers to whole horizontal domain; cloud active (CA) refers to all grid points that have at least one active vertical level in the column, that is, a level that has nonzero liquid cloud water, is positively buoyant, and is in an updraft (positive vertical velocity); and noncloud active (NCA) refers to grid points that do not meet the criteria for CA.
In terms of liquid water potential temperature, total specific content of water, and wind, the TKEM-TE performs best of all SCM configurations. The TKEM-GR and the TKEM-FL have in comparison a too-low PBL height due to too-weak mixing across the cloud top, and as a consequence, their PBL is too moist and too cold. The IFS-EDMF+SCP configuration is better in this respect thanks to the activity of the mass flux (both in the EDMF and in the SCP) but generates jagged profiles. Please note that the PBL height of the IFS-EDMF+SCP configuration is improved if the prognostic cloud scheme of the IFS (see ECMWF 2014) is activated.
The percentage of cloud fraction C and turbulence cloud fraction [see Eqs. (22) and (32)] for the TKEM-TE is comparable to the LES, but the cloud layer is too thin. The TKEM-GR overestimates and the TKEM-FL underestimates the cloud fractions. The nonlinear relationship between cloud fraction and turbulence cloud fraction is visible mainly for the TKEM-TE configuration but still does not reach the level of the LES.
The evolution of the vertical turbulent fluxes—including the mass and convective fluxes for the IFS-EDMF+SCP—for heat, moisture, and wind components is presented in Figs. 2 and 3. The TKEM-GR and, even more, the TKEM-FL show an on–off behavior in the diurnal evolution of all displayed fluxes. This behavior can be also observed in the IFS configuration. Compared to that, the TKEM-TE configuration is more continuous in time, showing better agreement with the onset and decay of shallow convection. This positive feature can be attributed to the prognostic character of the stability parameter in the TKEM-TE approach.
The TKE and the flux Richardson directly influence turbulent fluxes via the turbulent exchange coefficient [see Eqs. (13) and (14)], so their characteristics can be used to assess the behavior of the turbulence scheme. The evolution of TKE and the flux Richardson number is depicted in Fig. 4. Additionally, profiles of the three types of turbulence energy (, , and ) and three types of stability parameters computation (, , and ) for the TKEM-TE simulation are displayed in Figs. 5 and 6.
We can see that the TKEM-GR has a sharp transition between the mixing and nonmixing states because of its dependence on local gradients. The sharp transition leads to the accumulation of heat and moisture under the stable layer, which at some point results in the formation of clouds at the top of the PBL. The associated latent heat release enables turbulent mixing across the top of the PBL but only up to the next vertical level, where the new sharp transition is formed. The process repeats itself then on the next higher level, which results in the on–off behavior of the scheme.
The oscillating behavior of the TKEM-GR and the TKEM-FL is also caused by the direct link of the stability parameter to the local gradients or fluxes. If the local buoyancy flux exhibits a stable intermediate layer (see right in Figs. 5 and 6), the turbulent mixing in the PBL is interrupted until accumulation of heat and moisture from lower levels restarts it.
In comparison, the TKEM-TE has a smoother transition between stable and unstable layers and enables mixing across and into the locally stable layers (see right in Figs. 5 and 6). This is due to its link to the turbulence energies, which are vertically transported by the turbulence scheme and thus do not depend necessarily on the local gradients or the local turbulent fluxes (see left in Figs. 5 and 6). Thanks to that, the TKEM-TE is able to sustain continuous growth and later decline of the PBL.
In LES, the transport is caused mainly by intense transport in the active cells, which are not inhibited by local stability, and only the nonactive cells have a capped subcloud layer. Since the transport does not linearly depend on local stability, the LES appears to have turbulent transport across a locally stable layer in the domain average. This means that turbulent transport in the SCM grid should not depend only on local stability. This feature is mimicked in the TKEM-TE approach, which is confirmed by the good agreement of the TKE budget terms between the TKEM-TE and the LES. The comparison of TKE budget terms for one profile is displayed in Fig. 7.
b. Marine shallow cumulus
Next, we turn to BOMEX, representing a quasi-steady case of marine shallow cumulus. The case is a test whether our prognostic turbulence scheme is able to reach and hold a balance between turbulent mixing and a constant large-scale forcing.
The vertical profiles of selected variables are shown in Fig. 8. The TKEM-TE has the best agreement with the LES in this respect. The IFS-EDMF+SCP tends to overestimate and the TKEM-GR and the TKEM-FL underestimate the turbulent mixing in the upper part of the convective boundary layer.
The evolution of fluxes, presented in Fig. 9, shows that all the SCMs are able to hold a balanced state, which is accompanied by periodically rising thermals. The thermals in the TKEM-TE are too intense and less frequent than in the LES. This indicates that the speed of growth of the thermals is limited because of the cutoff formulation of the length scale above the PBL [see Eq. (18)]. The frequency and intensity of the thermals improves for smaller time steps (not shown). As in the ARM case, the penetration of the thermals is not deep enough for the TKEM-GR and the TKEM-FL. The IFS-EDMF+SCP has a good representation of the thermals, which are however too frequent in comparison with LES.
c. GABLS case
To check the versatility of our turbulence scheme, we tested it also for the stable GABLS case. For the sake of brevity, we present only the vertical profiles of wind speed and potential temperature after 7 h 30 min of integration (see Fig. 10).
All of the SCMs overestimate the turbulent mixing. Nevertheless, the profiles of potential temperature indicate the development of an elevated inversion within the SBL and generate a supergeostrophic jet peaking near the top of the boundary layer. These pronounced features in the LES are best modeled by the TKEM-FL and the TKEM-TE. The upper inversion in the IFS-EDMF+SCP has a smoother transition than the other SCMs.
The overestimation of mixing can be attributed to our particular choice of length scale, which might not be optimal for the conditions in GABLS case. Still, the results of our scheme are somewhat better than the operational IFS-EDMF+SCP scheme.
d. Countergradient fluxes
To evaluate the importance of the countergradient fluxes of heat and moisture in the shallow cumulus cases, we marked the countergradient regions in the LES data (see Figs. 1, 2, 8, and 9). The main countergradient egion is located above the PBL. The relatively small amplitude of the fluxes in this region indicates that their lack in the downgradient schemes does not lead to large errors. Another countergradient region can be found in PBL layers with near-constant liquid water potential temperature. The downgradient mixing in this region produces also a near-constant profile but enforces opposite signs of the vertical gradient and the turbulent flux. Altogether, we conclude that the downgradient approach is sufficiently accurate for the ARM and BOMEX cases. Similar observations regarding the downgradient schemes were made by Yamada (1977) for the Wangara case, Yamada and Mellor (1979) for the BOMEX case, and Bechtold et al. (1995) for the Rain in Cumulus over the Ocean (RICO) case.
6. Long-term 3D simulation
The SCM simulations allows us to study the turbulence scheme in idealized conditions independent from other parameterizations and the dynamics of the model, though the real challenge is to ensure that the turbulence scheme works well together with the other parts of the model and that it improves the forecast. We tested our turbulence scheme additionally in long-term three-dimensional global simulations. The aim of these tests was not an overall evaluation of its performance but the identification of possible weaknesses of the scheme, which might not be visible in the idealized simulations. Therefore, scores for selected variables only are presented.
We used for this the diagnostic tool of the ECMWF (for details, see Váňa et al. 2017; ECMWF 2017b), where a four-member ensemble simulation of the uncoupled global IFS model is run and compared to the model reanalyses and measurements for the 2000/01 time period.
We run the cy43r3 (ECMWF 2017c) version of the IFS model (i) with the IFS EDMF and SCP scheme (in figures marked as go7n) and (ii) the TKEM-TE turbulence scheme with two prognostic energies but without SCP (in figures marked as gugn). Apart from switching off the SCP, the following additional changes were made for the implementation of the TKEM-TE configuration into the full IFS model:
Two additional prognostic variables TKE and with inter–time step memory were introduced and are being transported by the semi-Lagrangian advection scheme.
Three more variables with inter–time step memory, but no advection (mainly for cost reasons), were introduced: turbulent heat flux, turbulent moisture flux, and the shear term in the TKE equation.
The turbulence scheme is computed with the same time step as the remaining physics and dynamics (i.e., no substepping is applied).
The existing coupling of physics and dynamics (see Wedi 1999) treats fast processes (i.e., turbulence and microphysics) by a first-order method while benefiting from a second-order averaging method along the semi-Lagrangian trajectory for the other processes. Thanks to the improved time consistency of the prognostic two-energy turbulence scheme, the second-order coupling can be now applied also to the tendencies from vertical diffusion. The remaining part of the fast-evolving processes requiring the first-order interface is then reduced to the tendencies of cloud water, cloud ice, and the two new turbulent energies only.
To stay consistent with the sequential time stepping and the present physics dynamics interface (Beljaars et al. 2004), the model variables seen by the physics (including now the variables used by the turbulence scheme) are updated by explicit dynamics and tendencies from other parameterizations that were already evaluated.
Finally, to stay consistent also with the cloud scheme, all the prognostic advectable quantities (prognostic cloud fraction, cloud ice and water, moisture, TKE, , wind components, and temperature) are interpolated by the same Lagrangian cubic interpolation.
The two-energy turbulence scheme was implemented in the same configuration as in the SCM simulations and without further adjustment of the IFS physics.
The evaluation of 2-m temperature and 2-m dewpoint temperature compared to ERA5 data (ECMWF 2017a) is displayed in Fig. 11. We can see that the impact of the turbulence scheme on the surface variables is comparable to the IFS-EDMF+SCP. As expected, a bigger impact of the turbulence scheme can be seen in Fig. 12, where the liquid water path compared to the Special Sensor Microwave Imager (SSM/I), version 6, dataset (see, e.g., O’Dell et al. 2008) and the total cloud cover (diagnosed outside the turbulence scheme) compared with CALIPSO data (see, e.g., Chepfer et al. 2008) are depicted.
The TKEM-TE has a slightly better score than the IFS-EDMF+SCP for the liquid water path. The IFS-EDMF+SCP tends to underestimate the cloud cover and the TKEM-TE to overestimate the cloud cover. The low-level clouds especially are overestimated with the TKEM-TE. This agrees with the results from the SCM simulations, where the IFS-EDMF+SCP is drier than the LES and the TKEM-TE in the upper part of the PBL (see Figs. 1 and 8). Thus, a suitable adjustment of the IFS cloud diagnostics, which is calibrated for the IFS-EDMF+SCP, could partly compensate the overestimation of the cloud cover for the TKEM-TE scheme. This is confirmed by the good agreement of the internal cloud fraction of the turbulence scheme [Eq. (27)] with the LES data in Figs. 1 and 8. However, the two-energy turbulence scheme still seems to transport too much moisture into the upper part of the PBL for the shallow convection cases. As discussed in the previous section for BOMEX, we believe that an improvement of the length-scale formulation could mitigate this problem.
An increased influence of skewness in the relationship between the cloud fraction C and the turbulence cloud fraction would also reduce the moisture and increase the temperature in the upper part of the PBL in the presence of shallow convection. Nevertheless, such modifications should be done with the consideration of the stratocumulus cases, where the distribution of temperature and moisture is almost Gaussian. However, the current TKEM-TE setup already tends to underestimate the cloud fraction for stratocumulus (not shown).
7. Summary and conclusions
A turbulence scheme with two prognostic turbulence energies suitable for all atmospheric stabilities has been derived and presented. The influence of moisture and phase changes is represented via a modification of the buoyancy terms in the presence of clouds. Thanks to the assumed linear dependence of the buoyancy terms on the covariances of heat and moisture, this influence can be expressed via a modification of the stability parameter only. By using the “stability-dependent adjustment for turbulence energy modeling,” the stability parameter can be computed from the two prognostic turbulence energies. This way, the stability parameter gains independence from the local stability conditions and obtains a prognostic character. The scheme is still downgradient, and the turbulent fluxes are proportional to the local gradients of the diffused variables, but the exchange coefficients are not strictly local.
The performance of the two-energy scheme was tested and compared with a standard TKE scheme and the IFS EDMF scheme with shallow convection parameterization in three idealized SCM simulations. The test for the ARM case shows that our scheme is able to model shallow convection with a diurnal cycle. The turbulent fluxes are more continuous in space and time compared to the TKE scheme and result in deeper mixing in accordance with the LES results. The cloud fraction has acceptable strong amplitude, but the cloud layer is too thin. Similar conclusion can be drawn for the quasi-steady shallow cumulus case BOMEX. However, the thermals modeled by our scheme are too intense and too infrequent compared to the LES reference. This is due to the cutoff formulation of the length scale. The problem can be partly remedied with an increased temporal resolution. The scheme performs also satisfactorily for the stable stratification case GABLS but slightly overestimates turbulent mixing.
The long-term three-dimensional global simulations indicate that our turbulence scheme has reasonable behavior in a full atmospheric model even though no calibration of the full model has been undertaken. In fact, our scheme is able to outperform the default IFS configuration in some characteristics, for instance, in the 2-m temperature in the Arctic region. However, the two-energy scheme tends to overestimate the cloud cover, especially at low levels. This can be partly attributed to the specific adjustment of the cloud-cover diagnostics in the IFS model. Still, our turbulence scheme seems to transport too much moisture into the upper part of the PBL in the presence of shallow convection.
Besides the abovementioned tests made in the global model, the two-energy scheme has been implemented first in the ALARO canonical model configuration of the ALADIN system (Termonia et al. 2018) and is currently operationally exploited at several ALADIN NWP centers.
The aim of this paper was to present and evaluate a turbulence scheme with two prognostic turbulence energies, which is one step above the commonly used TKE schemes. We have shown that this approach has good potential. There are still possibilities for improvement, especially in the length-scale formulation and the cloud fraction/buoyancy terms computation. We plan to improve both of these aspects in future work. A promising approach for the length-scale formulation is usage of relaxation technique for dissipation rate of TKE proposed in Zilitinkevich et al. (2013) or the computation of the length scale according to Nakanishi (2001).
Furthermore, we are aware that the selected set of test cases is limited. We plan to extend our investigations to stratocumulus cases and cases with precipitation.
The authors thank D. Mironov, S. Zilitinkevich, P. Marquet, P. Bechtold, E. Machulskaya, J. Mašek, A. Schomburg, M. Hrastinski, and P. Smerkol, who helped to shape this paper with their constructive comments and discussions. Thanks also to M. Ahlgrimm and the OpenIFS team at ECMWF for their support related to the SCM experiments.
Ivan Bašták Ďurán and Juerg Schmidli were supported by the Hans Ertel Centre for Weather Research. This German research network of Universities, Research Institutes and the Deutscher Wetterdienst is funded by the Federal Ministry of Transport and Digital Infrastructure (BMVI).
Support of the ALADIN community, especially in its RC LACE component (ALARO physics project), was essential for the completion of the work. This work was supported by the Czech Ministry of Education, Youth and Sports within the National Sustainability Program I, Grant LO1415.
The LES data for estimation of cloud fraction were provided by D. Lewellen.