Turbulent kinetic energy (TKE) is derived from a first-order planetary boundary layer (PBL) parameterization for convective boundary layers: the nonlocal K-profile Yonsei University (YSU) PBL. A parameterization for the TKE equation is developed to calculate TKE based on meteorological profiles given by the YSU PBL model. For this purpose buoyancy- and shear-generation terms are formulated consistently with the YSU scheme—that is, the combination of local, nonlocal, and explicit entrainment fluxes. The vertical transport term is also formulated in a similar fashion. A length scale consistent with the K profile is suggested for parameterization of dissipation.
Single-column model (SCM) simulations are conducted for a period in the second Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS2) intercomparison case. Results from the SCM simulations are compared with large-eddy simulation (LES) results. The daytime evolution of the vertical structure of TKE matches well with mixed-layer development. The TKE profile is shaped like a typical vertical velocity (w) variance, and its maximum is comparable to that from the LES. By varying the dissipation length from −23% to +13% the TKE maximum is changed from about −15% to +7%. After normalization, the change does not exceed the variability among previous studies. The location of TKE maximum is too low without the effects of the nonlocal TKE transport.
Mean turbulent kinetic energy (TKE) is one of the most important variables used to study turbulent boundary layers since it controls vertical mixing (Wyngaard and Coté 1971; Lenschow 1974; McBean and Elliott 1975; Yamada and Mellor 1975; André et al. 1978; Lenschow et al. 1980; Therry and Lacarrère 1983). Turbulence closure models can be classified by their order of closure. In large-scale and mesoscale atmospheric numerical models, the first-order and 1.5-order closure models are generally used to parameterize planetary boundary layer (PBL) processes (Cuxart et al. 2006; Svensson et al. 2011). Higher-order turbulence schemes, including 1.5-order schemes, use a TKE prognostic equation, and diffusion coefficients are calculated using TKE. Therefore, TKE can be directly evaluated and intercompared in these schemes (Holt and Raman 1988; Musson-Genon 1995; Cuxart et al. 2006; Svensson et al. 2011). On the other hand, TKE analysis in the first-order PBL models has not been conducted because TKE is not computed or used by these schemes.
In this paper, we introduce a derivation of TKE from a first-order nonlocal PBL model, the Yonsei University (YSU) PBL parameterization. One important goal of this research is to identify the characteristics of the K-profile YSU PBL scheme through the calculation of TKE, which has been computed for the TKE schemes. The YSU PBL parameterization and corresponding TKE prognostic equation are described in section 2. A description of single-column model (SCM) simulations and large-eddy simulations (LES) is provided in section 3. Section 4 shows the evaluation of the resultant TKE and TKE budget profiles compared to the LES results. Summary and concluding remarks follow in the final section.
2. YSU PBL parameterization and TKE prognostic equation
a. First-order, nonlocal YSU PBL (Hong et al. 2006)
In mesoscale and large-scale atmospheric models, impacts of vertical divergence (or convergence) of turbulent fluxes to prognostic mean variables are expressed by PBL parameterizations:
Here, the overbar and prime designate the mean and turbulent parts of variable c (e.g., u, υ, θ, qυ), respectively: . Basic assumptions in this formula are horizontal homogeneity and negligible mean subsidence.
In unstable conditions, the YSU PBL scheme parameterizes turbulent fluxes of heat, momentum, and moisture using the first-order (i.e., mean or grid) variables:
where Kc,KP is the K-profile diffusivity for variable c (the subscript KP is an abbreviation of K profile), γc is a nonlocal gradient-adjustment term applied for u, υ, and θ, and is the vertical flux of c at the PBL height h. The first term on the right-hand side of Eq. (2) refers to the vertical mixing proportional to the local gradient of , while the second one expresses the nonlocal mixing due to large convective eddies. The third term explicitly represents entrainment flux. The entrainment term is dependent on surface variables, and thus is also nonlocal. However, in this paper, the terminology “nonlocal” is used only for designating the term involved with γc. For more detailed explanations of each term, refer to Noh et al. (2003) and Hong et al. (2006). Note that this paper is limited to unstable conditions.
b. TKE prognostic equation for the YSU PBL
Under the assumptions of horizontal homogeneity and negligible mean vertical wind, the TKE (hereafter, the word “mean” is omitted for TKE) prognostic equation is written as (Stull 1988)
where ρ is the density, p the pressure, and g the gravitational acceleration. The term on the left-hand side represents the tendency. The two terms in the first set of parentheses on the right-hand side represent the vertical divergence of the vertical transport of TKE and pressure, respectively. The terms in the second parentheses and next term are shear production and buoyancy production, respectively. The viscous dissipation of TKE is ɛ.
The shear and buoyancy production/destruction terms are expressed using the flux parameterization of Eq. (2) for u, υ, θ, and qυ. Unlike the shear/buoyancy terms, the transport and dissipation terms require additional parameterizations to close the TKE equation. To describe the transport term, the parameterization method suggested by Therry and Lacarrère (1983, hereafter TL83) is combined with the explicit entrainment of TKE. From the third-order moment equation of TKE flux (i.e., ), TL83 derived a parameterization for TKE flux [cf. Eq. (26) of TL83]:
Here, αe, Ck, and C1 are numerical constants; lk and lɛ are a mixing length for momentum and a dissipation length used in TL83, respectively; and is the convective velocity scale. This equation is simplified as below:
Here, Ke is the diffusivity for TKE and γe represents nonlocal TKE transport involved in buoyancy effects. The variable γe is not a surface variable, but depends on z; it is expressed using in Eq. (2). We replace Ke,TL83 by KH,KP (the subscript H is for heat) to make consistency with the YSU PBL algorithm, and . TKE flux using Ke,TL83 is too small, since Ke,TL83 is about one-third of KH,KP in the interior of the PBL (not shown). Explicit entrainment flux of TKE is included to prevent TKE flux from being zero at h because of the K-profile diffusivity [cf. Eq. (8)]. The final formulation is as follows:
Here, Pr is the Prandtl number and the subscript M is an abbreviation of momentum. Equation (5c) expresses TKE flux at z = h [cf. Eq. (A11) in Hong et al. (2006)], which is proportional to the entrainment rate we and the jump in TKE at the inversion layer.
The dissipation is
Here, Λ1 and l are the dissipation length scale and master length scale, respectively; B1 is 16.6 (Mellor and Yamada 1982, hereafter MY82). Consequently, Cɛ ≈ 0.1704. A master length scale consistent with the K-profile model is computed; it is similar to that of MY82, but κz is multiplied by an amplification factor CKP that depends on surface-layer stability (more details in the next subsection):
Here, l0 is an asymptotic length scale, z is height above ground, κ is the von Kármán constant (0.4), and α1 is an empirical constant (0.3 in our simulations).
At each time step, vertical diffusion by the PBL scheme is conducted through Eqs. (1) and (2), and then TKE is computed [Eq. (3)]. We want to emphasize that the TKE calculation is a one-way method, and therefore does not affect the performance of the YSU PBL scheme. In association with the one-way feature the transport/dissipation terms are not exactly consistent with the YSU PBL scheme in a strict sense, because they are not necessarily in the forms above. The impacts of the formulations of the dissipation and transport will be discussed in sections 4b and 4c.
c. A master length scale for the K-profile model
where ws and are the mixed-layer velocity scale, and the friction velocity scale, respectively; ΦM ≡ (1 − 8z/L)−1/3 is the wind profile function, and L is the Obukhov length.
Our derivation of the master length scale for the YSU PBL scheme [cf. l1, Eq. (7)] focuses on allowing l1 in the surface layer to vary with surface-layer stability, instead of κz. In the surface layer, the K profile can be expressed as the product of a length scale and the frictional velocity:
Here, lKP is defined as a length scale of the K-profile model in the surface layer, which satisfies Eq. (10). The expression of the K profile using and lKP is not a standard formula, but it is used for an intermediate stage.
In the Mellor–Yamada (MY) model for the neutrally stratified boundary layer,
where SM = 0.55, and l2 is the master length scale in the MY model (MY82). The length scale is essentially for the neutral boundary layer. In the surface layer, KM,MY can be expressed similar to Eq. (10):
where lMY is defined as a length scale for the Mellor–Yamada model in the neutral surface layer, which satisfies the last equality in Eq. (13a). The approximation in Eq. (13a) is from surface-layer observations under near-neutral conditions (cf. Arya 2001, p. 208), , , and : ; σc is standard deviation of fluctuations of variable c.
When z ≪ h, and . Therefore, . This means that a length scale in the surface layer of the K-profile model defined under the convective/neutral conditions is approximately times larger a length scale of the MY model defined for the neutral conditions. Note that we replace z in ΦM [Eq. (9)] by ɛh. Here, ɛ is the ratio of surface-layer height to h, and it is set to 0.1. Therefore, the master length scale in the K-profile model (i.e., l1) is derived by allowing the length scale in the surface layer to vary with surface-layer stability as in Eq. (7):
Since the derivation of CKP is for z ≪ h, only κz is multiplied by the factor CKP.
Figure 1 compares CKP and an analogous amplification factor used by Nakanishi (2001) at z = ɛh: [cf. CN is a factor multiplied to κz in the master length scale suggested by Nakanishi (2001), as Eq. (20) in section 4b]. The term CKP is smaller than CN roughly for −h/L < 400, and larger for −h/L ≥ 400.
3. SCM simulations and LES
a. Model setup
We use the SCM version of the Weather Research and Forecasting (WRF) model (WRF-SCM), which is based on the WRF, version 3.2.1 [refer to Hacker and Snyder (2005) for WRF-SCM and to Skamarock et al. (2008) for WRFV3]. The WRF-SCM is set for a period in the second Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS2) intercomparison case (Svensson et al. 2011): 29 h from 1900 UTC 22 to 0000 UTC 24 October 1999. Initial vertical profiles of potential temperature θ, vapor mixing ratio qυ, and wind speed U are shown in Fig. 2a. Time evolution of surface temperature is prescribed as in the GABLS2 case, and surface heat fluxes and surface friction velocity are calculated using the Monin–Obukhov similarity theory with stability functions suggested by Jiménez et al. (2012). Figure 2b shows the calculated surface sensible heat flux H, latent heat flux (LH), and surface friction velocity . The daytime heat fluxes are smaller and decay more rapidly than those observed; is comparable to observations (not shown), but it is overestimated in stable conditions (cf. Kumar et al. 2010; Svensson et al. 2011).
It is noted that the direct comparison of the SCM results with the corresponding observations is less meaningful since the GABLS2 setup itself uses a simplified initial and boundary conditions, although they were derived from observations (Svensson et al. 2011). The SCM runs designed in this study follow the GABLS2 setup, but their results could diverge from the measurement. Thus, empirical reference curves (in the next subsection) and LES results will be used for the validation. A large-eddy simulation using the WRF-LES (e.g., Antonelli and Rotunno 2007; Moeng et al. 2007; Zhu 2008; etc.) is conducted over a 4 km × 4 km × 3 km domain; Δx = Δy = 25 m, and Δz increases with height between 20 and 30 m. Subgrid-scale (SGS) mixing is represented using the Smagorinsky-type SGS model [cf. chapter 4.2.3 in Skamarock et al. (2008)]. The LES run is initialized at 1800 UTC 23 October using the initial conditions based on the SCM simulation. The surface fluxes computed in the SCM and LES are almost identical to each other because the same surface temperature forcing and surface-layer scheme are used (not shown).
b. Reference (REF) TKE budget profiles
To evaluate TKE from the YSU PBL scheme, TKE budget profiles suggested by Lenschow (1974, hereafter L74) and Lenschow et al. (1980, hereafter LWP80) are used as references (REF) along with the LES results. BP, SP, TR, and DIS denote each TKE budget term (buoyant production/loss, mechanical production/loss, turbulent transport of TKE, and dissipation, respectively), normalized by surface buoyancy production .
In L74, BP is given as a function of as follows:
SP is given as
Lenschow used the assumption that momentum flux does not vary with height and has the value of . This works in situations with relatively small shear (L74). SP contributes less to TKE budget than BP does above surface layer, for −h/L > 10 (i.e., L74). This indicates that Eq. (15) does not cover situations of either strong shear across the top of the convective boundary layer (CBL; e.g., Conzemius and Fedorovich 2006) or relatively weak buoyancy forcing (e.g., −h/L < 3, from L74).
DIS is determined by LWP80 as shown below:
where a mixed-layer average is designated by angle brackets: a = 0.57/(〈SP〉 + 3.75) and b = 0.7.
For TR, two reference profiles are used. One includes the pressure-transport term and the other does not. L74 neglected the time rate of change of TKE when calculating TR. Using Eqs. (14)–(16), TR1 is
The transport of energy by pressure fluctuations is included in TR1.
On the other hand, LWP80 suggested a profile for normalized TKE flux and its vertical gradient is TR2. TR2 excludes the pressure effects:
It is known that the divergence of pressure flux works as a gain (loss) term at small (at mid-CBL) (Wyngaard and Coté 1971; McBean and Elliott 1975; Moeng and Wyngaard 1989) and it has the sign opposite the divergence of TKE flux. Therefore, the magnitude of TR2 is generally larger than that of TR1.
This model of the height variation of TKE budget is relatively simple but comparable to observations, as well as other higher-order turbulence model and large-eddy simulations [cf. Fig. 5.4 in Stull (1988)].
Our focus is on convective boundary layers. Hence, results will be presented for the last 12 h (from 1200 UTC 23 to 0000 UTC 24 October 1999) of the simulations. Note that the simulations begin at 1900 UTC 22 October according to the GABLS setup (cf. section 3a).
a. Overview of simulated boundary layer structures
Svensson et al. (2011) evaluated the performance of the YSU PBL scheme using the WRF-SCM for the GABLS2 case. Therefore, we restrict ourselves to briefly explaining some important characteristics of the simulated PBL. Vertical structures of the simulated virtual potential temperature and TKE of the control experiment (CTL; as described in sections 2b, 2c, and 3a) (cf. Table 1) are shown in Fig. 3. The YSU PBL scheme produces a weakly stable PBL during the daytime, as the LES does (Fig. 3a). This subadiabatic profile is a well-known characteristic of nonlocal PBL schemes (Holtslag and Boville 1993; Hong and Pan 1996), which also appears in the real atmosphere. In the CTL experiment, temporal evolution of the daytime boundary layer and TKE coincide (Fig. 3b).
Figure 4 shows TKE profiles averaged over 1 h centered at 1930 UTC (1330 LST). TKE from the LES slightly decreases with height because of the mean winds about 7.27 m s−1 (see Table 3) (cf. Stull 1988; Moeng and Wyngaard 1989). The vertical velocity w variance has its maximum at z/h = 0.4, together with the local maximum at 1.1h. The former is typical in CBL, and the latter occurs in the sheared CBL (Conzemius and Fedorovich 2006; Catalano and Moeng 2010). TKE for the YSU PBL is shaped similar to the LES w variance, and its maximum magnitude is comparable to the LES TKE (i.e., around 0.4h). Witek et al. (2011a) attributed the limited growth of near-surface TKE and the parabolic-like profile shape to the lower boundary conditions imposed in the TKE equation (cf. , where a and b are constants). The limitation of TKE production, in relation to the K profile, also results in the lack of TKE in the upper CBL (cf. section 4c).
b. Sensitivity to the master length scale
To estimate dissipation in the CTL experiment, the master length scale derived for the K-profile model l1 is used [cf. Eq. (6)]. The introduced TKE calculation is a one-way method such that calculated TKE is sensitive to the length scale, given the same mean meteorological fields. Therefore, two length scales are tested to assess uncertainty of the derived TKE owing to the length scale: the ML2 and ML3 experiments (cf. Table 1 and Table 2). In the ML2, the master length scale is computed as in MY82 [i.e. l2, Eq. (12)]. Comparing Eqs. (7) and (12), l2 applies to the neutrally stratified boundary layer. The amplification factor CKP > 1 for most unstable conditions (roughly for −h/L > 1: Fig. 1), so that l2 > l1. To reflect the large scale of convective eddies a larger l is often assumed in the CBL (Nakanishi 2001; Teixeira et al. 2004). Nakanishi (2001) modified the diagnostic equation for the master length scale of MY82 by introducing lS in the surface layer dependent on stability (i.e., dependent on ζ ≡ z/L) and lB as the length scale limited by stability (i.e., by N, the Brunt–Väisälä frequency) in the upper part of CBL:
where , and is a velocity scale with a reference temperature θ0 = 300 K; α2 = 1.0, α3 = 5.0, and α4 = 100.0 are empirical constants. This length scale l3 is used in the ML3 experiment.
In the lower CBL where ζ < 0 and lB → ∞,
Figure 5a shows the normalized dissipation length scale [cf. the dissipation length scale, Λ1 = B1l: Eq. (6)] in the CTL, ML2, and ML3 simulations. It is not surprising that the length scale for ML3 is larger than for CTL since CN = 2.95 > CKP = 2.05 for z = 0.1h and −h/L = 22.59. Dissipation decreases by increasing l (Fig. 5b), which results in a larger TKE (Fig. 5c). TKE maximum () increases from for ML2 to for CTL and to for ML3. The increased TKE strengthens TKE transport from the lower CBL (Fig. 5d). In summary, the predicted TKE varies between −15% and 7% by changing the dissipation length scale between −23% and 13%; the variation does not exceed the difference among the previous studies (Fig. 5c).
c. Role of nonlocal fluxes
To determine the role of nonlocal TKE flux in TR and heat flux in BP, the LTR (γe = 0) and LBP (γe = 0 and γθ = 0) experiments are conducted (Table 1). Dissipation is calculated using l1 as in the CTL [cf. the dissipation length scale, Λ1 = B1l: Eq. (6)].
TKE budget profiles for the CTL and LTR experiments are compared to the LES and REF profiles in Fig. 6. Note that in the TKE budget profiles, the TKE budgets at the lower boundary are excluded. This is because TKE is directly imposed by the lower boundary condition, rather than calculated by solving the TKE equation. Therefore, the underestimated near-surface TKE cannot be inferred from the TKE budgets (cf. Fig. 4). The simulated BP and SP resemble the LES and REF curves (Fig. 6a). This is because the heat and momentum fluxes of the YSU PBL scheme [cf. Eq. (2)] are designed to mimic the large-eddy simulations of convective boundary layers (Noh et al. 2003). The near-surface kink (that exceed 1) in BP for the LES is due to the fact that the WRF model uses a mass vertical coordinate (Moeng et al. 2007). For the CTL, the simulated TR (Fig. 6b) follows well the reference TR1 curve, which includes the pressure effects. However, it approaches zero near since TKE flux is parameterized using the prescribed K-profile diffusivity. The explicit entrainment also reduces TKE. In the CTL experiment, the calculated DIS is consistent with LWP80 (Fig. 6c), except for insufficient dissipation in the upper boundary layer above 0.6h. The underestimated DIS results from the limited TR at the PBL top. At a PBL top, a gain by turbulent and/or pressure transport balances losses by BP and DIS (Moeng and Wyngaard 1989).
Without the nonlocal term γe (i.e., the LTR experiment), the simulated TR approaches zero near the surface (Fig. 6b); upward TKE transport is not sufficient. Consequently, is located a little lower, at z ≈ 0.2h (Fig. 7a). Previous studies (e.g., TL83; Moeng and Wyngaard 1989; Witek et al. 2011a,b) argued that TKE flux is underpredicted when only an eddy-diffusion parameterization is used [cf. the gradient-diffusion term in Eq. (5a)], and an additional term is necessary. One exception is when an unphysically large diffusivity is used for the TKE flux (Moeng and Wyngaard 1989). This is because strong, organized updrafts themselves constitute a large part of TKE and they can redistribute TKE in a highly nonlocal way (Witek et al. 2011b). Our results (i.e., the impacts of γe) are consistent with the previous studies. Note that higher-order PBL models (e.g., the third-order models) explicitly predict TKE flux; therefore, the parameterization issue for the TKE flux is not relevant. Some second-order PBL models [e.g., the Mellor–Yamada–Nakanishi–Niino (MYNN) level 3.0 model (Nakanishi and Niino 2009)] also implicitly include the nonlocal TKE flux through stability functions used for calculating diffusivity.
In the present work, the inclusion of the nonlocal mixing in the buoyancy production term is essential, since it is used in the parameterization of temperature flux in the YSU PBL scheme [cf. Eq. (2)]. The simulated PBL is slightly stable in the upper part with the aid of the nonlocal heat transport (cf. Fig. 3a). In the case, the nonlocal buoyancy production is essential for preventing TKE damping above the middle of PBL (cf. Angevine et al. 2010). The LBP experiment (Fig. 7a) confirms the TKE limitation due to the negative buoyancy in the upper PBL. On the other hand, a well-known characteristic (or deficiency) of first-order or 1.5-order local PBL schemes is that they should maintain absolutely unstable profiles across the depth of CBL for upward heat transport (e.g., Holtslag and Boville 1993), while the observed upper boundary layer is slightly subadiabatic. Therefore, even though only the gradient-diffusion term is considered for buoyancy production, the buoyancy production term is positive in the upper part of boundary layer (i.e., because ). This results in the expected positive buoyancy production because of an inaccurate temperature profile.
We repeated five experiments shown in Table 1 without mean wind at the initial time, while other details of the experimental setup in section 3a were maintained. The Obukhov length L and −h/L at 1930 UTC are −2.54 and 367.9, respectively. Characteristics of TKE in the original GABLS2 case—impacts of l, role of nonlocal fluxes—are maintained (Fig. 7b). Since l1 and l3 have a very similar numerical value at that stability (cf. Fig. 1), the CTL and ML3 experiments produce similar TKE.
5. Summary and concluding remarks
We derive a turbulent kinetic energy profile from the first-order nonlocal YSU PBL parameterization. One important goal of this paper is to identify the characteristics of the YSU PBL scheme through the calculation of TKE-related variables that have been computed for standard TKE schemes (e.g., TKE, TKE budget, and length scale). A parameterization for TKE equation is developed to calculate TKE based on meteorological profiles of wind, temperature, and moisture from the YSU PBL model. For this purpose, buoyancy- and shear-generation terms are formulated by the combination of local, nonlocal, and explicit entrainment fluxes, consistent with the YSU scheme. The vertical transport term is also formulated in a similar fashion. To calculate dissipation length scale and parameterize dissipation, a master length scale consistent with the K profile is introduced. Note that this development is for convective boundary layers.
The results of the SCM simulations for a day from the GABLS2 case are compared with reference curves suggested by L74 and LWP80, as well as LES results. Several characteristics of the YSU TKE are discussed: 1) It is confirmed that the computed TKE develops well along with the vertical profile of the virtual potential temperature. In comparison with LES the vertical shape of TKE profile is similar to the shape of w variance, with maximum TKE comparable to the LES. 2) Changes in the dissipation length scale from around −23% to +13% result in the variation of estimated TKE from about −15% to +7%. The magnitude of the maximum TKE increases by enhancing the dissipation length scale as a direct effect of dissipation reduction. 3) A nonlocal TKE transport is necessary to reflect the TKE redistribution by strong organized updrafts.
Although there is no immediate benefit to incorporate a TKE budget into the YSU PBL scheme, the TKE, corresponding to turbulence in the YSU PBL is useful in evaluating the behavior of the scheme. This study introduces an example. Our results identify mixing characteristics of the widely used nonlocal K-profile PBL scheme through the TKE budgets, which can give us some knowledge of the scheme’s characteristics against prognostic TKE schemes. This approach can be applied to other non-TKE schemes. In a technical sense, this TKE calculation from a first-order scheme opens the way for coupling the scheme with shallow/deep cumulus convection modules that require information on the boundary layer TKE (e.g., Zhang and McFarlane 1995; Park and Bretherton 2009).
We appreciate Margaret A. LeMone and two anonymous reviewers, whose comments greatly improved our manuscript. This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (MEST) (2012-0000158). The first author would like to express her gratitude to Jinkyu Hong for his comments regarding the interpretation of TKE budget profiles.