## Abstract

A new turbulent kinetic energy (TKE)-based moist eddy-diffusivity mass-flux (EDMF) vertical turbulence mixing scheme (EDMF-TKE) is developed, where the nonlocal transport by large turbulent eddies is represented by a mass-flux (MF) scheme while the local transport by small turbulent eddies is represented by an eddy-diffusivity (ED) scheme, which is given by a function of a prognostic TKE. In the scheme, an MF approach is employed for the stratocumulus-top-driven downdrafts as well as for the thermals in the daytime unstable boundary layer. The scheme includes parameterizations for enhanced buoyancy due to moist adiabatic processes in condensation and for TKE interaction with cumulus convection. A scale-aware parameterization is proposed for the grid sizes where the large turbulent eddies are partially resolved. The single-column model (SCM) tests show that both the EDMF-TKE and the current operational NCEP GFS hybrid EDMF scheme (EDMF-CTL) simulate a daytime dry-convective boundary layer well and agree with a benchmark large-eddy simulation (LES). For a marine stratocumulus-topped boundary layer case, the EDMF-TKE better simulates the liquid water and wind speed profiles than the EDMF-CTL compared to the LES. For a stable boundary layer (SBL) case, the EDMF-TKE also agrees better with the LES than the EDMF-CTL, although it tends to produce a deeper SBL compared to the LES. On the other hand, three-dimensional medium-range forecast experiments show that the EDMF-TKE slightly improves forecast skill in the 500-hPa height anomaly correlation and wind vector, while it has a neutral impact on precipitation forecasts over the continental United States.

## 1. Introduction

Recently, the National Centers for Environmental Prediction’s (NCEP) Global Forecast System (GFS) has implemented a hybrid eddy-diffusivity mass-flux (EDMF) planetary boundary layer (PBL) scheme replacing an eddy-diffusivity countergradient (EDCG) mixing scheme that underpredicted the daytime convective boundary layer development (Han et al. 2016). In the hybrid scheme, the EDMF approach was applied to strongly unstable PBLs and the EDCG approach to weakly unstable PBLs. Use of the EDMF approach for weakly unstable PBLs resulted in overproduction of the amount of low clouds over the tropical ocean where strongly unstable PBLs are rarely found. This could be due to neglect of the enhancement of vertical turbulence mixing due to moist adiabatic processes such as condensational heating.

In the EDMF parameterization, the nonlocal transport by large turbulent eddies is represented by a mass-flux (MF) scheme and the local transport by small eddies is represented by an eddy-diffusivity (ED) scheme. The EDMF parameterization was introduced by Siebesma and Teixeira (2000) and since then this approach has been widely used in convective PBL modeling (Soares et al. 2004; Hurley 2007; Siebesma et al. 2007; Angevine et al. 2010; Witek et al. 2011; Sušelj et al. 2012) and successfully incorporated into weather and climate prediction models (Köhler et al. 2011; Hourdin et al. 2013; Sušelj et al. 2014; Han et al. 2016). A comprehensive and detailed description for the EDMF approach can be found in Siebesma et al. (2007).

Most of the EDMF schemes developed until now use a first-order (e.g., prescribed ED profile) or 1.5-order (e.g., turbulent kinetic energy closure) turbulence closure model for the ED scheme. The first-order closure models generally use a prescribed ED profile for the ED scheme (e.g., Siebesma et al. 2007; Köhler et al. 2011; Han et al. 2016). Most of the 1.5-order closure models are based on a prognostic turbulent kinetic energy (TKE) for the ED model (e.g., Soares et al. 2004; Witek et al. 2011; Sušelj et al. 2012). For higher-order accuracy, in this study we develop a new TKE-based moist EDMF parameterization for vertical turbulence mixing (hereafter EDMF-TKE), aiming to replace the current operational GFS EDMF PBL scheme using a prescribed ED profile (hereafter EDMF-CTL). Unlike the EDMF-CTL using a hybrid approach (i.e., EDMF for strongly unstable PBL and EDCG for weakly unstable PBL), the EDMF-TKE is applied to all unstable PBLs. The EDMF-TKE includes parameterizations for TKE interaction with cumulus convection as well as for enhanced buoyancy due to moist adiabatic processes in condensation. For the grid sizes where the large turbulent eddies (thermals) are partially resolved (e.g., subkilometer resolution), the EDMF-TKE assumes that the mass flux for the updraft thermals decreases with decreasing grid size (so-called, scale-aware parameterization).

## 2. Development of the TKE-based moist EDMF parameterization

The current operational NCEP GFS is a global spectral model for weather and climate prediction. As a part of the next generation global prediction system development, the hydrostatic spectral GFS dynamic core is being replaced by Geophysical Fluid Dynamics Laboratory’s nonhydrostatic Finite Volume Model, version 3 (FV3) (https://www.gfdl.noaa.gov/fv3), which will become NCEP’s next global forecast model in 2019. For the three-dimensional medium-range forecast experiments discussed later, we use the nonhydrostatic FV3 dynamic core.

For the model physics, the radiation parameterization uses the Rapid Radiative Transfer Model (RRTM) adapted from AER Inc. (e.g., Mlawer et al. 1997; Iacono et al. 2000; Clough et al. 2005). The orographic gravity wave drag parameterization is adopted from Kim and Arakawa (1995). The cloud condensate is a prognostic quantity with a simple cloud microphysics parameterization (Zhao and Carr 1997; Sundqvist et al. 1989; Moorthi et al. 2001). The fractional cloud cover used for radiation purposes is diagnostically determined by the predicted cloud condensate, following Xu and Randall (1996). The Noah land surface model (Ek et al. 2003) is used, and the surface layer similarity formulations are based on Long (1986, 1989).

The mass-flux parameterizations for deep and shallow cumulus convection (Han and Pan 2011) have been recently updated (Han et al. 2017) to have scale-awareness capability. Convective triggering was made more restrictive in order to reduce excessive drizzle, convective precipitation, and unrealistically noisy rainfall over high terrain. As noted earlier, the current GFS PBL model uses a hybrid EDMF parameterization for the convective PBL (Han et al. 2016), where the EDMF scheme is applied only for the strongly unstable PBL, while the EDCG scheme is used for the weakly unstable PBL. In this study, we update the PBL model using a TKE-based EDMF approach.

In the current GFS hybrid EDMF PBL scheme (Han et al. 2016), the vertical turbulent flux of a field *ϕ* is given by

for weakly unstable PBL,

for strongly unstable PBL, and

for the atmospheric layers above the mixed layer and nighttime stable boundary layer (SBL). Here overbars indicate horizontal averages across a grid cell; primes represent turbulent fluctuations; *K*_{sfc} and *K*_{Sc} are the surface and stratocumulus-top driven eddy diffusivity profiles, respectively; $K\u2061(Ri)$ is the eddy diffusivity given as a function of the gradient Richardson number (Ri); *γ*_{ϕ} is the nonlocal countergradient mixing term due to large nonlocal convective eddies; *M* denotes the mass flux; the subscript “*u*” refers to the updraft properties and the subscripts “sfc” and “Sc” imply “surface driven” and “stratocumulus top driven,” respectively.

In the TKE-based EDMF scheme,

where the subscript “*d*” refers to the downdraft properties. The eddy diffusivity *K*_{ϕ} is now given as a function of TKE, that is,

where *e* is the TKE, *l*_{k} is a turbulent mixing length (described later), and *c*_{ϕ} is a proportionality coefficient. For the unstable PBL, the coefficient for heat *c*_{h} is obtained from that for momentum *c*_{m} (i.e., *c*_{h} = Pr^{−1}*c*_{m}, where Pr is the Prandtl number). Similar to Bougeault and Lacarrere (1989), *c*_{m} is set to 0.4. For the stable conditions, *c*_{m} is obtained from *c*_{h} (i.e., *c*_{m} = Pr*c*_{h}), assuming that *c*_{h} = 0.4 for the SBL and *c*_{h} = 0.2 for the stable layers above the PBL as optimal values. The coefficients for the TKE and the other scalar variables such as moisture and tracers are assumed to be same as *c*_{h}. Within PBL, Pr is given as Pr = Φ_{h}/Φ_{m}, where Φ_{h} and Φ_{m} are the surface layer (*z* ≤ 0.1*h*, where *h* is the PBL height) nondimensional gradient functions for heat and momentum, respectively, which can be found in any micrometeorology text book (e.g., see Arya, 1988). We assume that Pr = 0.67 in the stratocumulus layers and the unstable layers above PBL. For the stable layers above PBL, Pr is given as Pr = 1 + 2.1Ri (Kennedy and Shapiro, 1980).

Note that the stratocumulus-top driven eddy diffusion in Eq. (1) is replaced by the downdraft mass-flux mixing in Eq. (2), and Eq. (2) is applied to all stability conditions [i.e., there is no separate equation depending on PBL stability like Eq. (1)]. TKE is predicted to obtain *K*_{ϕ}. This procedure is the so-called 1.5-order turbulence closure. The prognostic TKE can be written as

where *u* and *υ* are the horizontal velocities; *w* is the vertical velocity; *ρ* is the air density; *p* is the pressure; *θ*_{υ} is the virtual potential temperature; *g* is the gravity. The turbulent transport term is parameterized as

Here TKE is treated as a tracer, allowing a nonlocal MF transport. The turbulent momentum and heat fluxes in Eq. (4) are parameterized using Eq. (2), implying that the shear and buoyancy production of TKE is strongly influenced by the MF transport. The eddy diffusivity for TKE *K*_{e} is assumed to be same as those of heat and moisture. The TKE dissipation rate *D* is parameterized as

where *l*_{d} is the turbulent dissipation length scale (described later) and the constant *c*_{d} = 0.7. To obtain the TKE dissipative heating, on the other hand, Han et al. (2016) parameterized the TKE dissipation based on the assumption that shear and buoyancy production of TKE is balanced by TKE dissipation in Eq. (4). The TKE dissipative heating is now obtained directly from Eq. (6) as

where *T* is the temperature, *c*_{p} is the specific heat at constant pressure, and the coefficient *d*_{1} is set to 0.5.

The turbulent mixing length *l*_{k} is given by a combination of the surface layer length scale *l*_{1} and a characteristic length scale *l*_{2}:

Following Nakanish (2001), *l*_{1} is given by

where *κ* is the von Kármán constant; *a*_{1}, *a*_{2}, and *a*_{3} are stability-dependent coefficients; and *L* is the Monin–Obukhov length. Following Bougeault and Lacarrere (1989), the *l*_{2} and dissipation length scale *l*_{d} are given by

where *l*_{up} and *l*_{down} are the distances that a parcel having an initial TKE can travel upward and downward before being stopped by buoyancy effects, obtained from the equations:

Similar to Han et al. (2016), the mass flux for the updrafts is given by

where *a*_{u} (assumed to be *a*_{u} = 0.13) is a mean updraft area fraction in a large grid box where many large turbulent eddies reside. The updraft velocity *w*_{u} is computed as (Simpson and Wiggert 1969; Soares et al. 2004)

where *ε*_{u} is the lateral entrainment rate. The coefficients *b*_{1} = 2.0 and *b*_{2} = 4.0. Following Siebesma et al. (2007), *ε*_{u} is given by

where Δ*z* represents the vertical grid size and *c*_{ε} = 0.4 is an empirical coefficient.

The updraft properties for scalar fields are obtained using a simple entraining plume model as

The updraft property for momentum is given by

where **V** is the horizontal wind vector and *d*_{e} = 0.55 is an empirical constant representing the effect of the updraft-induced pressure gradient force (Zhang and Wu 2003; Han and Pan 2006). A more detailed discussion of Eqs. (12)–(16) is given in Han et al. (2016).

For parcel properties in condensation, moist adiabatic processes are considered and the liquid water potential temperature *θ*_{l} and total water *q*_{T }*= q* + *q*_{l} (where *q* and *q*_{l} are the specific humidity and the liquid water specific humidity, respectively) are used, which are conserved during both dry and moist adiabatic processes. The *θ*_{l} is defined by

where *L*_{υ} is the latent heat of vaporization of water, *P* is the layer pressure, and *P*_{0} is the surface pressure. The *θ*_{l} for the updraft (i.e., *θ*_{l,u}) is obtained using Eq. (15). The liquid water temperature for the updraft, *T*_{l,u}, is derived from the following relation:

When the rising parcel becomes saturated, the liquid water specific humidity of the parcel is obtained as (Sommeria and Deardorff 1977)

where *q*_{s} is the saturation specific humidity. The virtual potential temperature of the parcel in cloudy layers is computed as

where *θ*_{u} is obtained using Eq. (17). For the dry layers, *q*_{s} = *q*_{u} and *q*_{l,u} = 0 in Eq. (20).

The buoyancy at the lowest level is given by

where $(w\u2032\theta \upsilon \u2032\xaf)0$ is the surface virtual kinematic heat flux and the coefficient *c*_{1} = 1.0. The velocity scale *w*_{s} (Troen and Mahrt 1986) is represented by the value scaled at the top of the surface layer:

where *u*_{*} is the surface friction velocity, *α* is the ratio of the surface layer height to the PBL height (specified as 0.1), *κ* = 0.4 is the von Kármán constant, and *w*_{*} is the convective velocity scale defined as

The PBL height *h* in Eq. (14) is taken as the smaller of the height of *w*_{u} = 0 from Eq. (13) and a PBL height from Troen and Mahrt (1986), given by

where *U*(*h*) is the horizontal wind speed at *z = h*, *θ*_{υa} is the virtual potential temperature at the lowest model level, *θ*_{υ}(*h*) is the virtual potential temperature at *z = h*, and *θ*_{s} is the temperature scale near the surface defined as

where *θ*_{T} is a thermal excess proportional to the right-hand-side term in Eq. (21). The critical Bulk Richardson number (Rb_{cr}) is a constant for the unstable PBL as set to be Rb_{cr} = 0.25. For the SBL where *θ*_{s} is the virtual potential temperature at the surface without the thermal excess, however, Rb_{cr} varies with the surface Rossby number *R*_{0} as given by (Vickers and Mahrt 2004)

where *U*_{10} is the wind speed at 10 m above the ground surface, *f*_{0} is the Coriolis parameter, and *z*_{0} is the surface roughness length. To avoid too much variation, we restrict Rb_{cr} to vary within the range of 0.15–0.35.

Similar to the updrafts, the mass flux for the stratocumulus-top-driven downdrafts is given by

where *a*_{d} (assumed to be *a*_{d} = 0.12) is a mean downdraft area fraction. When the condition for cloud top entrainment instability (CTEI) is met, however, *a*_{d} is enhanced to 0.5. The condition for CTEI is given by (Randall1980; Deardorff 1980)

where $\Delta \theta e$ and $\Delta qt$ are the jumps in equivalent potential temperature and total water content across the cloud top. A constant *f*_{1} = 0.7 (MacVean and Mason 1990) is used, which is much more restrictive than that (*f*_{1} = 0.23) derived by Randall (1980) and Deardorff (1980) but supported by observational data.

The downdraft velocity *w*_{d} is computed as

where *ε*_{d} is the lateral entrainment rate for the downdrafts, given by

where *z*_{t} and *z*_{b} are the heights at the cloud top and the downdraft bottom, respectively. The downdraft properties for scalar fields are obtained as

The downdraft property for momentum is given by

A parcel descent from the cloud top *z*_{t} to determine *z*_{b} is made by perturbing the cloud top the liquid water potential temperature *θ*_{υl} by an amount equal to the cloud-top radiative cooling rate, multiplied by an assumed cloud-top residence time scale of 500 s (Lock et al. 2000). The final *z*_{b} is taken as the larger one between the level of *w*_{d} = 0 from Eq. (29) and the grid level at which this parcel’s *θ*_{υl} exceeds that of the environment. Parcel properties in condensation for downdraft are computed similar to those for the updraft using Eqs. (17)–(19).

Buoyancy flux by eddy diffusion is enhanced in the cloudy layers. Following Bretherton and Park (2009), the buoyancy is diagnosed from vertical gradients of the conserved variables in both dry and moist adiabatic process, *s*_{l }*= c*_{p}*T* + *gz* − *Lq*_{l} (liquid water static energy) and *q*_{T}, and the cloud fraction *c*_{f}. That is,

where *T*_{υ} is the virtual temperature. The thermodynamic coefficients *b*_{hs} and *b*_{qs} describe the contribution of the two conserved variable gradients to *N*^{2} in saturated air and similarly for *b*_{hu} and *b*_{qu} in unsaturated air, given as *b*_{hs} = *α*_{1}*β*, *b*_{hu} = *α*_{1}, *b*_{qs} = *α*_{1}*L(β* − *ϵ)*, and *b*_{qu} = *δg*. Here *α*_{1} = *g*/(*c*_{p}*T*_{υ}), *δ* = 0.608, and *β* = [1 + *γϵ*(1 + *δ*)]/(1 + *γ*) ≈ 0.5, where *ϵ* = *c*_{p}*T*/*L* and *γ* = (*L*/*c*_{p})∂*q*_{s}/∂*T.* Following Xu and Randall (1996), *c*_{f} is given as

where RH is the grid mean relative humidity.

The TKE equation [i.e., Eq. (4)] is integrated using time splitting to avoid a numerical instability as

TKE is defined on model half levels same as the other variables in the GFS and its fluxes at model bottom and top are assumed to be zero. To avoid an abrupt increase of TKE with a large model physics time step *∆t* due to strong wind shear, Eq. (37) is integrated with a smaller time step *∆t*′ (i.e., *∆t’ = ∆t*/*n*, where *n* is the number of integration, and *ê*^{t} is the *e* value after *n* integration) and with TKE and TKE dissipative rate *D* updated over every *∆t*′.

To represent TKE influence on cumulus convection, the entrainment rates in the cumulus updrafts and downdrafts are modified to be proportional to subcloud mean TKE. For example, the entrainment rate increases (decreases) by 30% when the subcloud mean TKE is larger than 0.65 m^{2} s^{−2} (less than 0.05 m^{2} s^{−2}). For the subcloud mean TKE of 0.05–0.65 m^{2} s^{−2}, the entrainment rate linearly increases with increasing subcloud mean TKE. To include the impact of cumulus convection on TKE, on the other hand, TKE is allowed to be transported by cumulus convection, similar to tracer transport by cumulus convection. In addition, TKE production by cumulus convection deduced from the cumulus mass flux *M*_{c} is added to TKE:

where $lc\sigma c$ is the cumulus area fraction.

Figure 1 displays the effects of the interaction of TKE and cumulus convection in the tropical regions where the cumulus convection is the most active. The increased TKE due to cumulus convection in upper troposphere (Fig. 1a) enhances the turbulent diffusion, which would have caused a reduced cloud condensate (Fig. 1b) and wind speed (Fig. 1d), but an increased relative humidity in the 300–500-hPa layers (Fig. 1c) where the relative humidity has a minimum value.

Conventional PBL models do not have a grid-size dependency because they assume that all the turbulent eddies are subgrid. For the grid sizes where the large turbulent eddies are partially resolved (e.g., subkilometer resolution), however, the subgrid-scale (SGS) turbulent fluxes decrease with increasing grid resolution (Honnert et al. 2011; Shin and Hong, 2013) and thus, a grid-size-dependent (i.e., scale-aware) parameterization will be necessary. Recently, Shin and Hong (2015) developed a scale-aware parameterization using spatially filtered data from their benchmark 25-m-resolution LES.

Since the present model uses an MF scheme for the nonlocal turbulent transport in the unstable PBL, in this study we propose a scale-aware parameterization based on the scale-aware cumulus convection parameterization, which also uses a mass-flux approach for the convective updrafts and downdrafts. Following Arakawa and Wu (2013) and Han et al. (2017), the mass flux for the updraft is assumed to decrease with increasing the updraft area fraction *σ*_{u} as

and

where *M*_{u} is the mass flux given by Eq. (12), *S*(*σ*_{u}) is a grid-size-dependent (scale-aware) function, and *M*′_{u} is a mass flux reduced with a finite *σ*_{u}, which is larger than *a*_{u} in Eq. (12). Following Grell and Freitas (2014) and Han et al. (2017), *σ*_{u} is parameterized as

where *A*_{grid} is the gridbox area, *R*_{u} is the radius of the updraft, and the lateral entrainment rate is averaged over the whole updraft from Eq. (14). Equations (40) and (41) imply that the nonlocal mass flux decreases with decreasing grid size for the subkilometer resolution. The similar scale-aware parameterization is also applied to the stratocumulus-top-driven downdraft.

The grid-size-dependent function using Eqs. (40b) and (41) for a typical PBL height of 1 km is displayed in Fig. 2 compared with that for the SGS nonlocal turbulent transport from Shin and Hong (2015). For the grid sizes less than *∆*/*h* = 0.5, the present scale-aware parameterization shows much faster reduction of the nonlocal turbulent fluxes with decreasing grid size than that of Shin and Hong (2015). Note that the grid-size-dependent function from Shin and Hong (2015) is obtained by fitting their own SGS turbulent transport model using an eddy diffusivity profile to the LES output and would significantly differ for different SGS models as indicated in Shin and Dudhia (2016).

The present scale-aware parameterization above has almost no impact on the current operational NCEP GFS of which grid size is about 13 km and would have a negligible effect until the grid size becomes less than *∆*/*h* ~ 0.7 as indicated in Fig. 2. A test or evaluation of the scale-aware parameterization requires a subkilometer resolution run, which is outside of the scope of this study and remains as a topic for future study.

## 3. Single-column model tests

To show the performance of the new TKE-based moist EDMF PBL scheme (EDMF-TKE) compared with the results from large-eddy simulations (LES) as well as the current hybrid EDMF PBL scheme in GFS (EDMF-CTL), we conduct the single-column model (SCM) simulations for three diverse boundary layer cases that have been used in past intercomparison studies. These include a convective boundary layer (CBL) with no mean wind, a marine stratocumulus-topped boundary layer, and a moderately stable boundary layer (SBL). These cases were also adopted in the study by Bretherton and Park (2009) for tests of their moist TKE parameterization.

### a. Convective boundary layer

Figure 3 compares the vertical profile of potential temperature simulated after 8 h by the LES with the GFS SCM results from EDMF-TKE, EDMF-CTL, and a TKE model without an MF term. For the simulations, an initial potential temperature profile *θ* = 288 K + (3 K km^{−1})*z* is used and a constant surface buoyancy flux (8 × 10^{−3} m^{2} s^{−3}) is specified. The LES model used here is the version 6 of the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall, 2003). The vertical grid size Δ*z* for the SCM is 50 m, which is the same as that in the LES. Figure 3a displays a typical potential temperature profile of a TKE closure scheme without the nonlocal MF term, showing a lack of a well-mixed CBL feature (i.e., unstable profile throughout the whole CBL) as well as an underprediction of the CBL growth compared to LES. Figure 3b shows that EDMF-TKE successfully simulates a daytime well-mixed CBL, giving excellent agreement with the LES as well as EDMF-CTL. Figure 4 compares the TKE profiles from the EDMF-TKE and LES at 8 h. Compared to the LES, the EDMF-TKE produces much weaker TKE (e.g., about 3 times weaker TKE at about 800-m height). This is because while the EDMF-TKE predicts only vertical component of TKE, the LES has extensive TKE in horizontal and vertical eddy motions. Despite much smaller TKE, the vertical mixing in the EDMF-TKE is very similar to that in the LES.

The SCM experiments were also conducted with the current coarse operational GFS vertical resolution, which is about Δ*z =* 170 m at *z =* 1 km and Δ*z =* 260 m at *z =* 2 km with 64 vertical levels. For the well-mixed CBL, both the EDMF-CTL and EDMF-TKE show good agreement with the LES even for the coarse resolution (Fig. 5), although they tend to slightly overestimate the CBL temperature.

### b. Marine stratocumulus-topped boundary layer

The mean LES result from the Second Dynamics and Chemistry of Marine Stratocumulus Experiment (DYCOMS2) Research Flight 1 (RF01) intercomparison (Stevens et al. 2005) is chosen as our second comparison with the SCM. For the simulations, the initial thermodynamic profile consists of a mixed layer up to a height of 840 m with a cloud layer in 600–840-m heights, topped by a sharp 7-K inversion and above that a stable layer. Surface fluxes, a mean horizontal divergence, and a vertically uniform geostrophic wind are specified. An idealized radiative flux profile is used. The detailed case specifications can be found in Stevens et al. (2005).

Figure 6 shows the SCM test results (averaged over 3–4 h) from the EDMF-TKE compared to the LES data and EDMF-CTL. The vertical grid size Δ*z* for the SCM is 10 m, which is the same as that in the LES. The liquid water profile from EDMF-TKE (Fig. 6a) agrees well with that of the LES, whereas the EDMF-CTL largely underestimates the LES liquid water profile, indicating that the EDMF-CTL is more diffusive. The wind speed profile from EDMF-TKE (Fig. 6b) also agrees well with that of the LES showing a well-mixed feature of momentum, whereas the EDMF-CTL fails to simulate the well-mixed momentum due to the lack of nonlocal momentum mixing in weakly unstable PBL. Both the EDMF-CTL and EDMF-TKE present well-mixed features of liquid water potential temperature but with slight overprediction compared to the LES (Fig. 6c). Compared to the LES, for the total water (Fig. 6d) the vertical mixing is somewhat too strong for the EDMF-TKE, whereas it is too weak for the EDMF-CTL especially in the lower PBL. Note that for the weakly unstable conditions of the present case, the EDMF-CTL does not have the nonlocal MF mixing for the scalar variables except temperature, whereas the EDMF-TKE has it for all the variables.

Figure 7 shows the SCM results with the coarse GFS vertical resolution (i.e., about Δ*z* = 170 m at *z* = 1 km and Δ*z* = 260 m at *z* = 2 km). EDMF-TKE agrees well with LES for the liquid water profile (Fig. 7a), although it tends to underestimate the LES liquid water more than the fine resolution SCM (Fig. 6a). EDMF-CTL is even more diffusive with much smaller peak value compared to the LES. For the wind speed (Fig. 7b), the EDMF-TKE agrees better with the LES than the EDMF-CTL, which is similar to the feature in the fine resolution (Fig. 6b). The liquid water potential temperature (Fig. 7c) and the total water (Fig. 7d) in the coarse resolution also display similar features to those in the fine resolution (i.e., overprediction of *θ*_{l} in both the EDMF-TKE and EDMF-CTL, overmixing in the EDMF-TKE and undermixing in the EDMF-CTL for the *q*_{T}). Further modification for the EDMF-TKE to reduce the overmixing in the weakly unstable condition is under way.

### c. Stable boundary layer

The Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study intercomparison case (GABLS1; Beare et al. 2006; Cuxart et al. 2006) is used to test our parameterization for a stable boundary layer. This is an idealized case for a moderately stable boundary layer. For comparison, we use ensemble-mean profiles from several LES model runs (Beare et al. 2006) at sufficiently high resolution that were quite similar. For the SCM simulations, the initial potential temperature equals to 265 K up to 100 m, and then it increases at a rate of 0.01 K m^{−1}. A vertically uniform geostrophic wind of 8 m s^{−1} is specified, and the latitude for the case is 73°N. The surface roughness is set to 0.1 m and there is no moisture. The surface is cooled at a constant rate of 0.25 K h^{−1} from the initial potential temperature of 265 K. The detailed case specifications are referred to Cuxart et al. (2006).

Figure 8 shows the SCM test results averaged over 8–9-h simulations for the EDMF-TKE and EDMF-CTL compared to the LES. The vertical grid size Δ*z* for the SCM is 6.25 m, which is the same as that in the LES. Both the EDMF-TKE and EDMF-CTL display the SBL depth larger than the LES (Fig. 8a) with wind speed peaks higher than the LES (Fig. 8b), while the EDMF-TKE has smaller SBL depth and lower height of the peak than the EDMF-CTL, showing a better agreement with the LES. The overestimated SBL depth is similar to the SCM results from the operational models in the GABLS1 (Cuxart et al. 2006). The abrupt changes of potential temperature gradients in both EDMF-TKE and EDMF-CTL just below inversion layers (Fig. 8a) are because the background diffusivities are set to be much smaller in the inversion layers and play dominant role in vertical diffusion for the strongly stable condition.

The background diffusivity *K*_{0} in the current operational GFS exponentially decreases with height from *K*_{0} = 1.0 m^{2} s^{−1} at the surface, given as

where *d*_{k} = 1.0. To better agree with the LES, we initially reduced *c*_{m} and *c*_{h} from the default value of 0.4, but that was not sufficient to reduce the SBL depth due to too large *K*_{0}. We find that profiles of potential temperature and wind speed similar to those of the LES can be obtained only with a substantially reduced *K*_{0} as well as reduced *c*_{m} and *c*_{h}, as shown in Fig. 9. For the three-dimensional forecasts such as medium-range forecasts in the following section, *d*_{k} in the new scheme is assumed to decreases with increasing grid resolution as

where *d*_{k0 }*=* 1.0 m^{2} s^{−1} and Δ*x*_{ref }*=* 25 000 m. Equation (43) yields *d*_{k} ~ 0.5 for the current global model grid size of Δ*x =* 13 km.

Figure 10 shows the SCM results with the coarse GFS vertical resolution (i.e., about Δ*z* = 170 m at *z* = 1 km and Δ*z* = 260 m at *z* = 2 km). Even with zero *K*_{0}, the SBL depth is overestimated in the coarse grid resolution. This indicates that a finer grid resolution is essential for the models to produce a more realistic SBL.

## 4. Medium-range forecast results

To assess the impacts of the new TKE-based moist vertical turbulent mixing scheme on forecast skill, 6-day forecasts every 5 days with the FV3 GFS were conducted for the period of 1 December 2016–6 December 2017. The initial forecast times were at 0000 UTC for all the forecasts. The GFS used in this test has 64 vertical sigma-pressure hybrid layers and the horizontal grid size is about 13 km. Since the forecasts were performed with no data assimilation, the analysis data from the operational GFS were used as the initial conditions. Although tests with data assimilation would be desirable, they are computationally expensive, and Park and Hong’s (2013) study seems to indicate that initial conditions are not always essential when evaluating the impact of model physics changes.

TKE usually displays a maximum near the surface due to large wind shear production of TKE (Fig. 11a), having a value of *O*(1) m^{2} s^{−2} or less. For a strong wind case such as a hurricane, however, the magnitude of the TKE near the surface can be on the same order of magnitude as the wind speed (Fig. 11b). Figure 12 displays a diurnal variation of PBL height, 2-m temperature, 2-m specific humidity, and 10-m wind speed above the ground surface with 6-hourly outputs for 5-day forecasts over an area of the Great Plains. Compared to the EDMF-CTL, the EDMF-TKE presents quite similar diurnal variation, although the differences between the two are become larger in later times. Figure 13 shows that although the EDMF approach is used for all the unstable PBLs, the EDMF-TKE only slightly increases low cloud amounts over the tropical regions compared to the EDMF-CTL. Note that the EDMF-CTL does not take into account moist processes, uses the EDCG approach for the weakly unstable PBL, and produces too much low clouds if the EDMF approach is used for all the unstable PBLs including the weakly unstable PBL (see Fig. 4b in Han et al. 2016). Figure 13 indicates that inclusion of moist processes in the EDMF-TKE may enhance the vertical turbulence mixing due to the enhanced buoyancy and consequently, prohibited an excessive accumulation of cloud condensate near the PBL top, yielding much less increase of low clouds compared to Fig. 4b in Han et al. (2016).

A comparison of anomaly correlations for the 500-hPa height, which illustrates how well synoptic-scale systems are represented over the globe, is shown in Fig. 14 as a function of forecast length. In both the Northern (20°–80°N) and Southern (20°–80°S) Hemispheres, the new scheme displays better anomaly correlations than the control forecast for the most forecast hours except for the forecast hour 24, although the improvement is not statistically significant especially for later forecast hours of 120 and 144 h.

Comparisons of the mean equitable threat and bias scores (Gandin and Murphy, 1992) for the 12–36-, 36–60-, and 60–84-h precipitation forecasts over the continental United States are shown in Fig. 15. Compared to the control forecasts, the new scheme does not display a significant difference for both equitable threat and bias scores. Although there is a slight improvement at the very high thresholds (i.e., thresholds of 50 and 75 mm day^{−1}), it is statistically not significant (not shown).

Figure 16 displays the wind vector root-mean-square errors (RMSE) for the tropics (20°S–20°N) and Northern (20°–80°N) and Southern (20°–80°S) Hemispheres as a function of height and forecast hour. Compared to the control forecasts, the new scheme generally produces smaller RMSEs for both the Northern and Southern Hemispheres. Over the tropics, the new scheme shows smaller RMSE than the control for the troposphere, whereas it displays larger RMSE than the control for the stratosphere.

## 5. Summary and conclusions

A new TKE-based moist EDMF vertical turbulence mixing scheme (EDMF-TKE) has been developed. Its advanced features over the current operational GFS vertical turbulence mixing scheme (EDMF-CTL) are as follows: 1) higher-order accuracy in turbulence representation, 2) advection of turbulence by the grid-mean flows, 3) inclusion of moist processes, 4) EDMF parameterization for the stratocumulus-top-driven turbulence mixing, 5) interaction of TKE with cumulus convection, and 6) scale awareness.

The SCM tests show that the EDMF-TKE well simulates daytime well-mixed PBL similar to the EDMF-CTL and agrees well with the LES result. For the marine stratocumulus-topped boundary layer, the EDMF-TKE better simulates the liquid water and wind speed profiles than the EDMF-CTL compared to the LES. In particular, the EDMF-TKE displays a well-mixed feature of momentum similar to the LES, whereas the EDMF-CTL fails to simulate the well-mixed momentum due to the lack of nonlocal momentum mixing. For the SBL, the EDMF-TKE agrees better with the LES than the EDMF-CTL, while it tends to produce deeper SBL compared to the LES. The LES results can be reproduced by largely reducing the empirical coefficient for diffusivity and background diffusivity *K*_{0} in the EDMF-TKE. The current GFS *K*_{0} appears to be too large and has been modified to decrease with increasing grid resolution in the new scheme.

Three-dimensional medium-range forecast experiments with the new scheme shows that the EDMF-TKE slightly improves forecast skill in 500-hPa height anomaly correlation, while it has a neutral impact on precipitation forecasts over the continental United States. For the wind vector forecasts, the new scheme generally improves the forecast skill in both the Northern and Southern Hemispheres as well as in the troposphere over the tropics, while it slightly degrades the skill in the stratosphere over the tropics.

As shown in this study, the new scheme not only agrees better with the LES compared to the current operational GFS scheme, but it also displays a positive impact on three-dimensional medium-range forecasts. In addition, since the new scheme has more information on turbulence with the prognostic TKE, it would have more potential for further improvement. The new scheme has been ported to the NCEP FV3 GFS system for possible operational use in 2021, replacing the current EDMF-CTL.

## Acknowledgments

This work was supported by the NOAA MAPP/CPO program as part of the Sc-to-Cu Transition Climate Process Team, and contribution by the second author, Chris Bretherton, was funded by NOAA grant NA13OAR4310104. The authors highly appreciate the help from Fanglin Yang and Ruiyu Sun at NCEP/EMC for the GFS forecast and postprocessor runs. Internal reviews from Shrinivas Moorthi and Ruiyu Sun at NCEP/EMC are also greatly appreciated. We are also grateful to Joao Teixeira at NASA/JPL for his encouragement on this work.

## REFERENCES

*Introduction to Micrometeorology*. Academic Press, 303 pp.

*14th Symp. on Boundary Layers and Turbulence,*Aspen, CO, Amer. Meteor. Soc., 4.16, https://ams.confex.com/ams/AugAspen/techprogram/paper_14840.htm, .

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).