Parameterization of Sheared Entrainment in a Well-developed CBL. Part II: A Simple Model for Predicting the Growth Rate of the CBL

The depth of the CBL is an important parameter in air pollution and NWP models. In these models, the CBL height is often diagnosed from the critical bulk Richardson number or TKE profile (Skamarock et al., 2008). These diagnosed CBL heights exhibit large intra-scheme variances (Shin and Hong, 2011; Xie et al., 2012; Breuer et al., 2014). Thus, a proper method to estimate the CBL height is important for numerical models. The growth rate of a CBL is actually the entrainment rate when there is no background vertical velocity. In the 1960s, (Ball, 1960) and (Lilly, 1968) proposed a bulk model framework to describe the evolution of the CBL. Since then, the bulk model approach has been widely used to predict the CBL entrainment rate. For a sheared CBL, the results of LESs in (Kim et al., 2003) showed that the turbulence in the entrainment zone enhances due to the break of Kelvin-Helmholtz billows at the upper edge of the entrainment zone. As a result, the entrainment process accelerates and the IL deepens. In order to adequately capture the entrainment process in a sheared CBL, at least the first-order model (FOM) is needed (Kim et al., 2006; Conzemius and Fedorovich, 2007), and the set of equations for the CBL should be solved. Unfortunately, the bulk model is complex and difficult to apply in numerical models since it includes too many unknown variables. Therefore, it is imperative to develop a simple model for predicting the growth rate of a sheared CBL.

The entrainment rate is associated with the entrainment flux ratio A _e, which is defined as the ratio of heat flux at the CBL top to that at the surface. (Kim et al., 2006) proposed a parameterization of A _e for the sheared CBL in the FOM framework. They only considered the sheared CBL under the condition of height-constant geostrophic velocity (GC). However, the LES results in (Conzemius and Fedorovich, 2006) indicated that the entrainment process has different characteristics under the condition of sheared geostrophic velocity with a zero value at the surface (GS). Following the derivation in (Kim et al., 2006), (Liu et al., 2016) developed a parameterization of A _e for a well-developed and sheared CBL. This scheme takes into account the buoyancy effect and the shear effect in the surface layer, the mixed layer and the IL. The shear effect in the IL is represented by local momentum fluxes and velocity jumps at the CBL top, which is similar to that in (Pino and De Arellano, 2008). It still includes many variables and thereby cannot be applied directly. However, if the relations between the entrainment variables (e.g., the entrainment momentum fluxes, the potential temperature jump, and the velocity jumps) can be described by external parameters (such as the background stratification, the geostrophic wind gradient, and the surface friction velocity), the parameterizations of A _e and the entrainment rate can be simplified and become applicable. For this reason, the present study aims to develop a simple model appropriate for predicting the growth rate of the sheared CBL by using external parameters.

The parameterization of the entrainment rate suggested by (Sun and Xu, 2009) and the A _e scheme proposed in (Liu et al., 2016) are simplified with some assumptions in this study. The LES outputs obtained in (Liu et al., 2016) are used for analysis and evaluation. The paper is organized as follows: Section 2 discusses the theoretical considerations in the FOM framework, and a simplified parameterization for the entrainment rate is proposed. In section 3, the characteristics of the sheared entrainment are analyzed based on the LES results, and the proposed simple parameterization is verified using the LES outputs. Conclusions and discussion are given in the final section.

In this study, Θ, U and V represent horizontally averaged potential temperature and velocity components, while θ, u and v represent the fluctuation parts of potential temperature and the velocity components. \(\overline{w\theta}\), \(\overline{uw}\) and \(\overline{vw}\) represent horizontally averaged vertical fluxes of potential temperature and velocity components. γ_θ and γ_u represent the vertical gradients of the initial potential temperature and the geostrophic velocity in the x-direction. γ_θ and γ_u are external parameters, which are assumed to remain unchanged during CBL development. Figure 1 depicts the idealized profiles of Θ, U and V, and their fluxes, in a well-developed sheared CBL. The CBL height h₁ is defined as the level at which \(\overline{w\theta}\) reaches its minimum. h₀ is the first zero-crossing height of the \(\overline{w\theta}\) profile. h₂ is defined as the level toward which \(\overline{w\theta}\) becomes larger than 10% of its minimum. The layer between h₁ and h₂ is the inversion layer (IL), in which the idealized Θ increases with height. The IL thickness is ∆ h₂₁=h₂-h₁. The layer between h₀ and h₂ is the entrainment zone, in which \(\overline{w\theta}\) is negative. The entrainment zone thickness is ∆ h₂₀=h₂-h₀. Θ₁, U₁ and V₁ are values of Θ, U and V at h₁. Θ₂, U₂ and V₂ are values of Θ, U and V at h₂. The potential temperature jump is ∆Θ=Θ₂-Θ₁, and the two components of the velocity jump are ∆ U=U₂-U₁ and ∆ V=V₂-V₁. The A _e is defined as \(A_e=-\overline{w\theta}_1/\overline{w\theta}_s =h_1/h_0-1\). \(\overline{uw}_1\) and \(\overline{vw}_1\) are the momentum fluxes at h₁, which are obtained by integrating the momentum equations from the surface to h₂. For further details, please refer to (Liu et al., 2016).

Figure 1.  Idealized profiles of horizontally averaged potential temperature $\Theta$, velocity components $U$ and $V$, and their vertical fluxes $\overline{w\theta}$, $\overline{uw}$ and $\overline{vw}$ in the GS cases. Dash-dot lines represent $h_0$, $h_1$ and $h_2$, and dotted lines represent zero.

In the FOM framework and under the GC condition, the parameterization of the entrainment rate (w _e) and A _e proposed by (Kim et al., 2006) and evaluated by (Pino et al., 2006) (hereafter KP06) can be expressed as \begin{equation} \label{eq1} w_{{e,KP}}=\dfrac{\partial h_1}{\partial t}=\dfrac{\Delta h_{21}+(2h_1+\Delta h_{21})A_{{e,KP}}}{h_1(2\Delta\Theta-\gamma_\theta\Delta h_{21})} \overline{w\theta}_{s} (1)\end{equation} and \begin{eqnarray} \label{eq2} A_{{e,KP}}&=&\dfrac{-\overline{w\theta}_1}{\overline{w\theta}_{s}}=A_{1{,KP}}\dfrac{1}{1+\frac{\Delta h_{21}}{h_1}}+ A_{2,{KP}}\dfrac{1}{1+\frac{\Delta h_{21}}{h_1}}\dfrac{u_\ast^3}{w_\ast^3}+\nonumber\\ &&A_{3{,KP}}\dfrac{-\frac{1}{2}\overline{uw}_1\Delta U-\frac{1}{2}\overline{vw}_1\Delta V}{\left(1+\frac{\Delta h_{21}}{h_1}\right)w_\ast^3} , (2)\end{eqnarray} where \(\overline{w\theta}_s\) is the surface heat flux; w_* is the convective velocity scale in the shear-free CBL, defined as \(w_\ast^3=(g/\Theta_0)\overline{w\theta}_sh_1\); u_* is the surface friction velocity; and A_{1 ,KP}=0.20, A _2,KP=0.26, and A _3,KP=1.44. In order to conveniently compare with the parameterization scheme proposed in (Liu et al., 2016), here we do not express Eq. (3) in the same form as in KP06, in which \(-\overline{uw}_1\Delta U/2-\overline{vw}_1\Delta V/2\) is replaced by a derived relationship [(i.e., Kim et al., 2006, Eq. (5)].

Following the work of (vanZanten et al., 1999), (Sun and Wang, 2008) derived a parameterization of the entrainment rate for a shear-free CBL, in which w _e is normalized by w_* and is proportional to the inverse of the convective Richardson number. They argued that their scheme is still valid for the sheared entrainment process because the scheme is derived from the profiles of potential temperature and its flux, which have the same shape in both sheared and shear-free CBLs. This argument is supported by the LES outputs for a sheared CBL (Sun and Xu, 2009). The alternative form of (Sun and Wang, 2008) given in (Sun and Xu, 2009) is expressed as \begin{equation} \label{eq3} w_{{e,SW}}=\dfrac{\Delta h_{21}+(h_1+\Delta h_{21})A_{{e,SW}}}{h_1\Delta\Theta}\overline{w\theta}_{s} , (3)\end{equation} and \begin{equation} \label{eq4} A_{{e,SW}}=\dfrac{-\overline{w\theta}_1}{\overline{w\theta}_{s}}=\dfrac{\Delta h_{10}}{h_0} , (4)\end{equation} where ∆ h₁₀=h₁-h₀ is the thickness of the lower part of the entrainment zone. As pointed out in (Sun and Xu, 2009), KP06 and (Sun and Wang, 2008) are equivalent in describing sheared entrainment, although the expressions are different. KP06 explicitly includes the effects of shear-produced turbulence, whereas (Sun and Wang, 2008) uses a geometric relation to implicitly represent the wind shear effects.

In both KP06 and (Sun and Wang, 2008), the effect of wind shear on the entrainment rate is represented by the A _e. It should be noted that Eq. (3) is derived for GC CBLs. A parameterization of A _e in a well-developed and sheared CBL is derived in (Liu et al., 2016). The expression reads: $$ A_e=A_1 \dfrac{1}{1+\frac{\Delta h_{21}}{h_1}}+A_2 \dfrac{C^{-1/2}_D}{1+\frac{\Delta h_{21}}{h_1}} \frac{u^3_*}{w^3_*}+$$ $$ Term I \ \ \ \ \ \ \ \ \ \ Term II$$ $$\ \ \ \ A_3\dfrac{(-\frac{1}{2} \overline{uw}_1 \Delta U-\frac{1}{2} \overline{vw}_1 \Delta V)}{(1+\frac{\Delta h_{21}}{h_1})w^3_*}+$$ $$\ \ \ \ \ \ \ \ \ \ Term III$$ $$\ \ \ \ \ \ \ \ \ \ A_4 \dfrac{(V_1-V_S)(-frac{1}{2} \overline{vw}_S-\frac{1}{2} \overline{vw}_1+\frac{1}{12}f\gamma_u h_1^2}{(1+\frac{\Delta h_{21}}{h_1})w^3_*}, (5)$$ $$\ \ \ \ \ \ \ \ \ \ Term IV$$ where V _s is the velocity in the y-direction at 0.1h₁ (the top of the surface layer), and C _D=u_{*²/(U s²+V s²)} is the surface drag coefficient. Note that the geostrophic velocity in the y-direction is zero in this study. The coefficients have been determined in (Liu et al., 2016), i.e., A₁=0.21, A₂=0.01, A₃=0.86, and A₄=0.70.

The surface drag coefficient in Eq. (5) is not a constant. It is worth noting that an additional term that represents the effect of wind shear in the mixed layer is included in Eq. (5). In the GC case, Eq. (5) reduces to Eq. (3). On the other hand, Eq. (5) is the explicit form of Eq. (5) for a sheared CBL. The combination of Eq. (4) and Eq. (5) can be used to predict the evolution of the sheared CBL depth. However, with so many variables in the two equations, they are inconvenient for application. In this study, we attempt to develop a simplified scheme based on characteristics of sheared entrainment obtained from the LES outputs.

The relation between ∆Θ/∆ h₂₀ and γ_θ is often used as a parameter to characterize the thermal structure in the entrainment zone (Deardorff, 1979; Fedorovich, 1995; Gentine et al., 2015). It is called the relative stratification parameter and is defined as \begin{equation} \label{eq6} G=\dfrac{\gamma_\theta}{\Delta\Theta/\Delta h_{20}} . (6)\end{equation}

The LES results in (Fedorovich et al., 2004) indicate that the dimensionless parameter G is a constant of around 1.2 for shear-free entrainment. Whether it is a constant for sheared entrainment has not been discussed in previous studies. A relationship between G and A _e is derived from the equation of temperature by order analysis using our LES results (details of the derivation given in Appendix A). It is expressed as \begin{equation} \label{eq7} G=\frac{1+K_1A_{e}}{1+A_{e}} , (7a)\end{equation} where the parameter K₁ is defined as K₁=2-∆ h₂₁/2∆ h₂₀. In order to further simplify the problem, we assume ∆ h₂₁≈ ∆ h₂₀/2 (details about G and the validation of the approximation are discussed in the next section). Thus, K₁ can be treated as a constant of 7/4, and Eq. (8) can be written as \begin{equation} \label{eq8} G=\frac{1+\frac{7}{4}A_{e}}{1+A_{e}} . (7b)\end{equation} Equation (7b) indicates that G is a function of A _e rather than a constant. Substituting Eqs. (6) and (7b) into Eq. (3) and applying the relation ∆ h₂₀=∆ h₂₁+∆ h₁₀, we get \begin{equation} \label{eq9} w_{e}=\left(1+\dfrac{7}{4}A_{e}\right)\frac{\overline{w\theta}_{s}}{\gamma_\theta h_1} . (8)\end{equation} In (Liu et al., 2016), a new convective velocity scale (w _m), which includes the contributions of shear-produced TKE in the whole CBL, is proposed. The results in (Liu et al., 2016) indicate that w _m is suitable for characterizing the vertical turbulent motion at the mixed layer top. It is expressed as \begin{eqnarray} \label{eq10} w_{m}^3&=&w_\ast^3+\frac{A_2}{A_1}C_{D}^{-1/2}u_\ast^3+\frac{A_3}{A_1}\left(-\dfrac{1}{2}\overline{uw}_1\Delta U-\dfrac{1}{2} \overline{vw}_1\Delta V\right)+\nonumber\\ &&\dfrac{A_4}{A_1}(V_1-V_{s})\left(-\dfrac{1}{2}\overline{vw}_{s}-\dfrac{1}{2}\overline{vw}_1+\dfrac{1}{12}f\gamma_uh_1^2\right) . (9)\end{eqnarray} Then, Eq. (5) can be rewritten as \begin{equation} \label{eq11} A_{e}=\dfrac{A_1}{1+\Delta h_{21}/h_1}\dfrac{w_{m}^3}{w_\ast^3} . (10)\end{equation} The variation of ∆ h₂₁/h₁ is relatively large during the development of a CBL. However, the variation of 1+∆ h₂₁/h₁ is so small (Table 1) that it can be approximately treated as a constant. The LES outputs show that the average value of 1+∆ h₂₁/h₁ is 1.19. Since the value of A₁ is set to 0.21 (Liu et al., 2016), Eq. (10) is further simplified to be \begin{equation} \label{eq12} A_{e}=0.18\frac{w_{m}^3}{w_\ast^3} . (11)\end{equation} It is worth noting that Eq. (11) is not the scheme in the zeroth-order jump model (ZOM), although the expression has the same form as that in the ZOM. As shown in Eq. (9), the parameterization is still complex and includes many unknown variables. Thus, how to simplify Eq. (9) is the problem we need to solve.

(Liu et al., 2016) shows that the effects of wind shear in the surface layer and mixed layer [the second and fourth terms on the right-hand side of Eq. (5)] on the entrainment flux are quite small. Since the terms on the right-hand side of Eq. (9) are proportional to the terms on the right-hand side of Eq. (5), it can be concluded that the second and fourth terms on the right-hand side of Eq. (9) are small when compared with the first and third terms. The ratio of the sum of the second term and the fourth term to w _m³ is smaller than 5% in most simulated cases (Table 1). For the purpose of simplicity, the two small terms are neglected in the simplified parameterization. Therefore, the simplification of Eq. (9) begins with the third term, hereafter expressed as \(X_3=-\overline{uw}_1\Delta U/2-\overline{vw}_1\Delta V/2\). In order to derive ∆ U, ∆ V, \(\overline{uw}_1\) and \(\overline{vw}_1\), the velocity equation is vertically integrated from the surface to the level just above the entrainment zone top. Rearrangement of the integrations yields (see details of derivation in Appendix B) \begin{eqnarray} \label{eq13} h_1\dfrac{\partial U_1}{\partial t}&=&\overline{uw}_{s}-\overline{uw}_1\!+\!\frac{1}{2}fh_1[(V_{s}+V_1)-(V_{{g,s}}\!+\!V_{g,1})] ,(12)\nonumber\\\\ \left(1+\dfrac{\Delta h_{21}}{2h_1}\right)\overline{uw}_1&=&\dfrac{\Delta h_{21}}{2h_1}\overline{uw}_{s}- \left(\Delta U-\dfrac{1}{2}\gamma_u\Delta h_{21}\right)\dfrac{\partial h_1}{\partial t}+\nonumber\\ &&\dfrac{1}{4}f\Delta h_{21}[(V_{s}-V_1)-(V_{{g,s}}-V_{{g},1})] ,(13)\\ \dfrac{\partial \Delta U}{\partial t}&=&\gamma_u\dfrac{\partial h_1}{\partial t}-\dfrac{\partial U_1}{\partial t} . (14)\end{eqnarray} By using the differential transition, such as \(\partial/\partial t=\partial/\partial h_1\cdot\partial h_1/\partial t=(S/h_1)\partial/\partial h_1\), where \(S=\overline{w\theta}_s[1+(7/4)A_e]/\gamma_\theta\) comes from Eq. (8), Eq. (14) becomes \begin{equation} \label{eq14} \frac{S}{h_1}\frac{\partial \Delta U}{\partial h_1}=\gamma_u\frac{\partial h_1}{\partial t}-\frac{\partial U_1}{\partial t} . (15)\end{equation} Substituting Eqs. (12) and (13) into Eq. (15), we rearrange the equation to be: \begin{eqnarray} \frac{\partial \Delta U}{\partial h_1}&=&-\frac{1}{2h_1+\Delta h_{21}}(2\Delta U-2\gamma_u\Delta h_{21}-2\gamma_uh_1)-\nonumber\\ &&\frac{2h_1}{2h_1+\Delta h_{21}}\frac{\overline{uw}_{s}}{S}-\frac{1}{S}\frac{h_1}{2h_1+\Delta h_{21}}f\Delta h_{21}(V_1-V_{{g},1})-\nonumber\\ &&\frac{1}{S}\frac{h_1}{2h_1+\Delta h_{21}}fh_1(V_{s}-V_{{g,s}}+V_1-V_{{g,1}}) . \label{eq15} (16)\end{eqnarray} The solution of Eq. (16) can be obtained, in which the coefficients are substituted in for the purpose of simplicity (see details in Appendix C). It reads: \begin{eqnarray} \label{eq16} \Delta U&=&0.57\gamma_uh_1-0.48\frac{h_1}{S}\overline{uw}_{s}-\nonumber\\ &&\frac{fh_1^2(0.19V_1-0.19V_{{g},1}+0.16V_{s}-0.16V_{{g,s}})}{S} .\qquad (17)\end{eqnarray} To solve Eq. (17), we assume that ∆ U depends only on h₁ and 1+∆ h₂₁/h₁=1.19. Substituting Eq. (17) into Eq. (13) gives a simple expression of \(\overline{uw}_1\). It is written as \begin{eqnarray} \label{eq17} \overline{uw}_1&=&0.57\overline{uw}_{s}-0.47\gamma_uS+\nonumber\\ &&fh_1(0.14V_1-0.14V_{{g},1}+0.20V_{s}-0.20V_{{g,s}}) ,\quad (18)\end{eqnarray} where V _g,s and V _g,1 are the geostrophic velocity at the surface and CBL top, respectively.

The expressions of ∆ V and \(\overline{vw}_1\) are derived in the same way as ∆ U and \(\overline{uw}_1\). For the CBL under the GC condition, V is constant while \(\overline{vw}\) varies linearly with height in the mixed layer. Their profiles present the same characteristics as that of U and \(\overline{uw}\). Thus, the expressions of ∆ V and \(\overline{vw}_1\) have the same form as ∆ U and \(\overline{uw}_1\). According to Eq. (17) and Eq. (18), they can be written as (the terms associated with the vertical gradient of the geostrophic velocity in the y-direction γ_v diminish because γ_v=0) \begin{equation} \label{eq18} \Delta V\!\!=\!-0.48\frac{h_1}{S}\overline{vw}_{s}\!+\!\frac{fh_1^2(0.19U_1\!\!-\!0.19U_{{g},1}\!\!+\!0.16U_{s}\!-\!0.16U_{{g,s}})}{S} (19)\end{equation} and \begin{equation} \label{eq19} \overline{vw}_1=0.57\overline{vw}_{s}-fh_1(0.14U_1-0.14U_{{g},1}+0.20U_{s}-0.20U_{{g,s}}). (20)\end{equation} Therefore, in the GC case, the third term in Eq. (9) can be written as (the terms associated with γ_u diminish because γ_u=0 in the GC case) \begin{equation} \label{eq20} X_{3,{GC}}=-\frac{1}{2}(\overline{uw}_1\Delta U+\overline{vw}_1\Delta V)=0.14\frac{h_1}{S}u_\ast^4+F_{f,{GC}} , (21a)\end{equation} where F_{f, GC} represents the sum of the terms, including the Coriolis parameter (expression not shown). In order to estimate X_{3, GC}, Eq. (9) with the exact coefficients is written here (the second and fourth terms have been omitted): \begin{equation} \label{eq21} w_{m}^3=w_\ast^3+4.10X_3 . (9')\end{equation} The LES results show that w _m³=w_*³+5u_*³, the formula proposed by (Moeng and Sullivan, 1994), is a good approximation for describing w _m under the GC condition. Comparing Eq. (9') with Eq. (21a), it is found that 0.14(h₁/S)u_*⁴+Ff, GC≈ u_*³, which can lead to w _m³≈ w_*³ +5u_*³.

However, for a CBL under the GS condition, V is not constant in the mixed layer, while the shape of the \(\overline{vw}\) profile is different from that in the GC case. The expressions of ∆ V and \(\overline{vw}_1\) should be different to those in the GC case. The LES results show that, in the GS case, \(\overline{vw}_1\Delta V\) is far smaller than \(\overline{uw}_1\Delta U\), and thus can be neglected (Table 1). That is, \(\overline{uw}_1\Delta U+\overline{vw}_1\Delta V\approx\overline{uw}_1\Delta U\). The LES results also show that, in the GS case, \(\overline{vw}_s\approx 0\) and \(\overline{uw}_s\approx-u_\ast^2\). Thereby, X_{3, GS} can be obtained from Eq. (17) and Eq. (18). It is written as \begin{eqnarray} \label{eq22} X_{3,{GS}}&=&-\frac{1}{2}\overline{uw}_1\Delta U\nonumber\\ &=&0.28u_\ast^2\gamma_u h_1+0.14\gamma_u^2Sh_1+0.14\frac{h_1}{S}u_\ast^4+F_{f,{GS}} ,\qquad\ (21b)\end{eqnarray} where F_{f, GS} represents the sum of the terms, including the Coriolis parameter (expression not shown). When the geostrophic velocity gradient becomes zero (i.e., γ_u=0), X_{3, GS} should reduce to X_{3, GC}. Based on this consideration, we assume that the approximation 0.14(h₁/S)u_*⁴+Ff, GS≈ u_*³ is also appropriate for the GS case. Thus, the expression of X₃ for both the GC and GS cases is \begin{equation} \label{eq23} X_3=\eta u_\ast^3+0.28\gamma_u u_\ast^2h_1+0.14\gamma_u^2Sh_1 , (22)\end{equation} where η is an empirical constant introduced by the above approximations. Replacing S in Eq. (22) by the relation \(S=\overline w\theta_s[1+(7/4)A_e]/\gamma_\theta\), and further replacing A _e by using Eq. (11), then substituting Eq. (22) into Eq. (9'), we get the following equation: \begin{equation} \label{eq24} \left(1-0.37a_3\frac{\Theta_0\gamma_u^2}{g\gamma_\theta}\right)w_m^3=\left(1+a_3\frac{\Theta_0\gamma_u^2}{g\gamma_\theta}\right)w_\ast^3+a_1 u_\ast^3+a_2u_\ast^2\gamma_u h_1 , (23)\end{equation} where a₁=4.10η, a₂=1.15 and a₃=0.57. Since some approximations and assumptions are introduced in the derivation, these coefficients should be adjusted. Multiple linear regression of the LES outputs in the GC and GS cases gives a₁ =6.02, a₂ =0.24 and a₃=0.86.

(Pino and De Arellano, 2008) used local momentum fluxes and velocity jumps to represent the shear production rate of TKE. Equations (17) and (18) indicate that these local quantities are related to surface fluxes, entrainment rate, CBL height and geostrophic velocity gradient. Equation (22) also indicates that the net shear production rate of TKE in the IL comprises the dynamic effect (second term), the interaction between mean shear and environmental stratification (third term), and one other term that includes the Coriolis effects (first term).

Equations (8), (11) and (23) constitute a simple model to predict the height of a sheared CBL. Compared with the KP06 and Sun and Wang (2008) methods, this model does not need those variables at the CBL top, such as ∆Θ , ∆ h₂₁, and so on. Owing to the coarse resolution of NWP and air pollution models, these variables cannot be resolved explicitly (Hong et al., 2006). Thereby, they must be obtained from a bulk model. However, the input parameters in our simplified model can be easily derived in NWP and air pollution models. \(\overline{w\theta}_s\) and u_* can be estimated by a land surface model. γ_θ can be derived as the mean gradient of potential temperature within a certain thickness (for example, 1000 m) above h₂. γ_u can be treated as the mean gradient of velocity above h₂ and determined by the same method for the calculation of γ_θ. In the real atmosphere, the treatment of γ_u may introduce some errors. However, the CBL's growth is mainly controlled by surface heating, which explains 70%-90% of this process (Canut et al., 2010; Sühring et al., 2014). This implies that the contribution of entrainment to CBL growth is about 10%-30%. It is expected that the errors induced by γ_u cannot significantly influence the prediction accuracy. Therefore, this model is more convenient to apply in numerical models.

The outputs of the 30 LES runs, described in (Liu et al., 2016), are used to verify the parameterizations under the conditions of varying geostrophic wind or wind shear, external stratification, and surface roughness length. The main features of these simulations are described in detail in (Liu et al., 2016). There are two simulations for the shear-free CBL (NS00S3 and NS00S6), while all other simulations are for the sheared CBL and divided into three groups. One group is for CBLs under the GC condition, with vertically uniform geostrophic velocities of 10 m s^-1, 15 m s^-1 and 20 m s^-1, respectively. The second group is for CBLs under the GS condition, with geostrophic wind gradients of 10 m s^-1 (2 km)^-1, 15 m s^-1 (2 km)^-1 and 20 m s^-1 (2 km)^-1, respectively, and zero surface geostrophic velocity. In each group under the GC and GS conditions, the simulations are conducted with external temperature gradients of 3 K km^-1 and 6 K km^-1, and surface roughness length values of 0.01 m and 0.1 m, respectively. The name of a simulation case is given according to the simulation conditions. For example, GC20R6 means that the simulation is conducted under conditions of U _g=20 m s^-1, a rough surface with z₀=0.1 m, and ∂Θ/∂ z=6 K km^-1. In section 2, the derivation of Eq. (23) is based on assumptions and approximations that are obtained from results of the GC and GS cases, and the coefficients are fitted from simulations of these 24 cases. In order to confirm their validity, a third group of simulations are conducted under the CS condition (C5S10S3, C5S15S3, C5S15S6 and C5S15R3). The CS condition can be regarded as a combination of the GC and GS conditions, in which the geostrophic wind shear exists while the surface geostrophic velocity is not zero. In the four CS cases, the surface geostrophic velocity is 5 m s^-1 (denoted as C5), while its vertical gradient is 10 m s-1 (denoted as S10) or 15 m s^-1 (denoted as S15) per 2 km, and the meanings of the last two letters in the case names are the same as for the GC and GS cases. The integration covers 28 800 s, and the results from 4800 s to 28 800 s are output at an interval of 100 s for further calculations and analyses. The methods to determine the variables used in calculations and analyses are introduced in (Liu et al., 2016).

According to the definition expressed in Eq. (6), the LES outputs are used to calculate the relative stratification parameter G, and the results are shown in Fig. 2 (blue dots). In the GC cases, G has a slightly decreasing trend during CBL development. However, the decrease is negligibly small (only about 0.1 in a long period of 24 000 s). In the GS cases, G almost does not vary with time. As shown in Table 1, the average value of G varies slightly in different cases. It increases with increasing geostrophic velocity in the GC cases and increasing geostrophic velocity gradient in the GS cases. However, the difference in G among the simulated cases is very small. The mean value of G in all of the simulated cases is 1.17 (Table 1), which is very close to the result for shear-free CBLs in Fedorovich et al. (2004, therein G≈ 1.2).

Figure 2.  Comparison of $G$ (the relative stratification parameter) calculated from the definition and from the parameterization scheme in the GC and GS cases. The blue dots represent the results from Eq. (6) and the red dots represent the results from Eq. (7b).

The parameterization of the relative stratification parameter, i.e., Eq. (7a), is derived in the FOM framework. As shown in Appendix A, an approximation is applied to simplify Eq. (A6). The LES outputs are used to calculate every term on the left-hand side of Eq. (A6), and results indicate that the second and third terms are one order smaller than the other two terms. The LES outputs provide a solid basis for Eq. (7a). However, Eq. (7a) still includes variables at the CBL top (∆ h₂₁ and ∆ h₂₀). A further simplified parameterization of G, Eq. (7b), is obtained by using the assumption of ∆ h₂₁≈∆ h₂₀/2, which may not be true, but the LES outputs show that K₁ is almost a constant and approximately equals 7/4 (Table 1). In order to evaluate the performance of Eq. (9), the LES outputs are used to calculate G according to Eq. (7b). The results are also plotted in Fig. 2 (red dots) to compare to those calculated according to Eq. (6) [Fig. S1 in electronic supplementary material (ESM) for CS cases]. The spread of LES outputs (especially for GC20 cases when the integration time is longer than 20 000 s) is mainly attributed to the determination method of the upper edge of the entrainment zone, h₂, where the instantaneous potential temperature flux profile is not equal to zero but fluctuates around zero. Significant fluctuations of ∆ h₂₁ can lead to large fluctuations of ∆Θ. However, the fluctuation of ∆ h₂₀=∆ h₂₁ +∆ h₁₀ is relatively small. A small value of ∆Θ always corresponds to a large value of G calculated from Eq. (6). Figure 2 shows that the results from Eq. (6) and Eq. (7b) agree very well in all of the simulation cases, indicating that Eq. (7b) can accurately describe the behavior of G.

The parameterization of entrainment rate, Eq. (8), is derived from Eq. (3) and Eq. (7b). (Sun and Xu, 2009) demonstrated that Eq. (3) works in sheared CBLs. Equation (7b) is validated in the previous section. Thus, we expect that Eq. (8) can predict the CBL height correctly. In order to prove this point, Eq. (10) is solved numerically using the Euler predictor-corrector method over the period from 5700 s to 28 800 s under the CS condition. The time step is 100 s, the same as the time interval of LES outputs. The parameterization schemes of A _e and the convective velocity scale, Eqs. (11) and (23), are used to close Eq. (8). The LES outputs for the period from 4800 s to 6600 s are averaged to provide the initial condition at 5700 s. Results are illustrated in Figs. 3-5.

Figure 3.  Convective velocity scales calculated from the CBL height prediction model [red line, Eq. (23)] and from LES outputs [blue dots, Eq. (9)] in CS cases.

Figure 4.  $A_e$ values calculated from the CBL height prediction model [red line, Eq. (11)] and from LES outputs [blue dots, Eq. (4)] in CS cases.

Figure 5.  Values of CBL height predicted by the model (red line) and obtained from LES outputs (blue dots) in CS cases. Green squares show the nonlinear fit of the result of LES CBL height based on the relation that the CBL height is proportional to the square root of time.

Figure 3 shows that Eq. (23) slightly underestimates and overestimates w _m at the beginning and end, respectively. This difference is mainly due to the approximation of 0.14(h₁/S)u_*⁴+F_f ≈ u_*³, and a larger u_* gives a larger error. However, the biases of estimated w _m are very small, making it reasonable to conclude that the simplified form of the convective velocity scale Eq. (23) agrees well with its original form, Eq. (9). Figure 4 shows that the A _e estimated by Eq. (11) is in good agreement with that derived from the LESs. As presented in previous studies (e.g. Conzemius and Fedorovich, 2006; Pino et al., 2006), the LES A _e spreads widely because it is determined from instantaneous LES profiles (calculations show that the spread of LES A _e is reduced significantly when the LES heat flux profiles are averaged over 500 s). It is satisfactory that the values of the parameterized A _e are contained within the fluctuations of the LES outputs. Figure 5 indicates that Eq. (8) can correctly predict the CBL height.

(Fedorovich et al., 2004) suggested that, in a shear-free CBL, the CBL height is proportional to the square root of time, i.e., \(h_1\propto\sqrt t\). The direct numerical simulations of the CBL driven by a constant momentum flux in (Jonker et al., 2013) yielded the same result, i.e., \(h_1 \propto \sqrt t\). Our LES results show that this relation also exists in the GC, GS and CS cases (Fig. 5 only displays the results in the CS case; see Fig. S2 in ESM for the GC and GS cases).

The entrainment rate is often parameterized as w _e/w_*=ARi*-1, where Ri_*=(g/Θ₀)h₁∆Θ/w_*² is the bulk convective Richardson number. In the ZOM, the coefficient A is just the A _e. In the FOM, this expression still applies, but the convective Richardson number is different to that in the ZOM because of the different definition of ∆Θ. Determination of the coefficient A has been explored (e.g., Lewellen and Lewellen, 1998; Sullivan et al., 1998; vanZanten et al., 1999; Sun and Wang, 2008). (Sun and Wang, 2008) gave a relation of A=(1+A _e)∆ h₂₀/h₁, indicating that A is associated with not only the A _e but also the entrainment zone thickness and CBL depth. The parameterization scheme for the entrainment rate in (Sun and Wang, 2008) can be written as Eq. (3). KP06 provided an equivalent scheme, expressed as Eq. (1). Both Eq. (1) and Eq. (3) are valid for a sheared CBL, and both imply that the effect of wind shear can be represented by the A _e. On the other hand, the two schemes both include the variable ∆Θ, which is usually unknown. This problem is solved by introducing G. A simple relationship between G and A _e, expressed as Eq. (7b), is derived in this study; and the parameterization of the entrainment rate turns out to be Eq. (8).

For a sheared CBL, the influence of wind shear is included in the A _e. Equation (11) further indicates that the effects of wind shear can be represented by the convective turbulent velocity scale in a sheared CBL, which is always larger than that in a shear-free CBL. Equation (23) gives a parameterization of the convective velocity scale, in which the CBL bulk variables, such as ∆ U, ∆ V, ∆ Θ and ∆ h₂₁, which cannot be resolved explicitly in most numerical models, are not needed. This is the fundamental difference to the parameterization of KP06.

Actually, Eq. (3) is not the final form of the parameterization of A _e in KP06. The final form is expressed as

\begin{equation} \label{eq25} A_{{e,KP}}=\dfrac{A_{1,{KP}}\frac{h_1}{h_2}+A_{2,{KP}}\frac{u_\ast^3}{(w'_\ast)^3} +A_{3,{KP}}\frac{\Delta h_{21}}{2h_1+2h_2}\left[\frac{u_\ast^2 \Delta V_{e}}{(w'_\ast)^3}+ \frac{\Theta_{m}(\Delta V_{e})^2}{gh_2(\Delta\Theta-0.5\gamma_\theta\Delta h_{21})}\right]} {\left[1-A_{3,{KP}}\frac{\Theta_{m}(\Delta V_{e})^2}{2gh_2(\Delta\Theta-0.5\gamma_\theta\Delta h_{21})}\right]} , (24)\end{equation}

where Θ _m is the mean potential temperature in the mixed layer, \((w'_\ast)^3=g\overline{w\theta_s}h_2/\Theta_m\), and ∆ V _e is the velocity jump across the IL. (Liu et al., 2016) shows that the value of 1.44 for A_{3, KP} overestimates the contribution of shear-produced TKE in the IL to entrainment. Our LES outputs are used to optimize the parameters in Eq. (26). In (Kim et al., 2006), ∆ V _e=0.5(|∆ U|+|∆ V|), and the linear regression yields A_{1, KP}=0.22, A_{2, KP}=-0.54 and A_{3, KP}=2.39. In (Pino et al., 2006), $(\Delta V_e=\sqrt{(\Delta U)^2+(\Delta V)^2})$, and the linear regression yields A_{1, KP}=0.22, A_{2, KP}=0.23 and A_{3, KP}=0.74. Because A_{2, KP}=2C _D^-1/2(1-α₂), A_{3, KP}=2(1-α₃), where α₂ and α₃ represent the proportions of the dissipation rate to the production rate, and C _D is the surface drag coefficient; A_{2, KP} and A_{3, KP} must be positive and A_{3, KP} must be less than 2. The negative value of A_{2, KP} and the larger than 2 value of A_{3, KP} imply that the treatment of ∆ V _e in (Kim et al., 2006) is not reasonable. In the following calculation, only the ∆ V _e in (Pino et al., 2006) is adopted. In order to identify the differences among (Sun and Wang, 2008), LS (LS is the scheme in this study) and the optimized KP06 schemes, the LES outputs are used to calculate w _e in the simulated cases, and the relative errors are compared in Fig. 6. The relative error is defined as \begin{equation} \label{eq26} {Err}=\frac{1}{n}\sum\left|\dfrac{w_{{e,p}}}{w_{{e,LES}}}-1\right| , (25)\end{equation} where w _e,p is the predicted entrainment rate, \(w_e,LES=\partial\langle h_1\rangle/\partial t\), and \(\langle h_1\rangle\) is the least squares nonlinear fit of LES h₁ based upon the relation \(\langle h_1\rangle\propto \sqrt t\) (Fig. 5). In order to understand whether the error of KP06 comes from the parameterization scheme of the entrainment rate [namely, Eq. (2)] or the parameterization scheme of the A _e [i.e., Eq. (26)], the errors of the entrainment rate predicted by Eq. (2) and the simplified A _e parameterization scheme [i.e., Eq. (13)] (denoted as KP _LS) are also calculated and shown in Fig. 6. Among all the three schemes, the LS scheme performs best, implying that the approximations applied in the derivation of Eq. (25) are reasonable and the derived parameterization can successfully capture the characteristics of the entrainment rate in a sheared CBL. Note that the LS scheme performs better than the (Sun and Wang, 2008) scheme, although the former is developed on the basis of the latter by using some approximations. The reason is because the spread of instantaneous LES variables (∆Θ, ∆ h₂₁, and so on) used in the (Sun and Wang, 2008) scheme is larger than that of the instantaneous LES variable (namely, h₁) used in the LS scheme. KP _LS performs better than the KP06 scheme with optimized parameters but slightly worse than the (Sun and Wang, 2008) scheme. This suggests that Eq. (2) performs slightly worse than Eq. (4) and the A _e estimated by Eq. (26) has large errors in some cases, even though the parameters have been optimized. The errors of Eq. (26) are partially due to the wide spread of instantaneous LES variables used in Eq. (26). Based on the above discussion, it is concluded that the assumptions used in the derivations are reasonable and the simplified parameterizations proposed in this study can correctly predict the A _e and entrainment rate. Meanwhile, the simplified parameterizations proposed in this study do not include entrainment variables, which may introduce large calculation uncertainties.

Figure 6.  Relative errors of the entrainment rate predicted by different schemes against the LES outputs in the simulated cases.

This study aims to simplify the parameterization of the entrainment rate so that it can be conveniently applied to predict the growth rate of a well-developed and sheared CBL. To achieve this goal, G (the relative stratification parameter) is introduced into the parameterization scheme, and the A _e and convective velocity scale in a sheared CBL are simplified according to the characteristics of sheared entrainment, which are derived from the LES outputs. The major findings can be summarized as follows:

(1) G, which is defined in Eq. (6), can be used to characterize the thermal structure in the entrainment zone. (Fedorovich et al., 2004) suggested that it is a constant of around 1.2. Theoretical analysis indicates it is a function of the A _e as expressed in Eq. (7b), rather than a constant. This result is supported by the LES outputs.

(2) When the relationship between G and the A _e is introduced, the parameterization scheme of the entrainment rate proposed by (Sun and Wang, 2008) can be rewritten as a simple function of A _e, surface heat flux, and the potential temperature gradient in the free atmosphere and the CBL height, as expressed in Eq. (8).

(3) Shear-produced TKE at the CBL top enhances entrainment. (Pino and De Arellano, 2008) used the local momentum fluxes and velocity jumps to parameterize the shear production rate of TKE. Our results show that these local quantities are related to surface momentum fluxes, the entrainment rate, the geostrophic velocity gradient, the stratification in the free atmosphere, and the CBL height. The net shear production rate of TKE in the IL is expressed by Eq. (24), which comprises the dynamic effect, the interaction between mean shear and environmental stratification, and one other term that includes the Coriolis effect, and can be approximately characterized by u_*³. Based on Eq. (22), the convective velocity scale proposed in (Liu et al., 2016) can be further expressed by Eq. (23).

(4) In the framework of FOM, the IL thickness affects the A _e. The LES outputs show that this effect can be described by a constant, and thus the parameterization of A _e proposed in (Liu et al., 2016) is simplified by Eq. (11).

(5) The parameterizations of the entrainment rate, A _e and convective velocity scale constitute a prediction model for CBL height. The LES outputs show that it is an appropriate model for a well-developed and sheared CBL. Compared with bulk models, the parameters needed by this model can be easily derived in an NWP or air pollution model. Therefore, this model is more convenient for application. However, the performance in NWP and air pollution simulations needs further validation.

As pointed out in (Liu et al., 2016), the parameterizations obtained in this study may only be suitable for well-developed CBLs under idealized conditions. The storage term in the TKE budget is ignored in the derivations. However, this term is not negligibly small in the early stage of CBL development. Under this condition, the performance of the simple model needs to be verified. In this study, the potential temperature gradient in the free atmosphere and the geostrophic velocity gradient are assumed to be constant. For a CBL that is growing through a pre-existing inversion or a residual layer, application of Eq. (8) becomes problematic, especially when there is a residual layer, because the relative stratification parameter may not exist. How to parameterize the entrainment rate under these conditions needs further investigation. In the real atmosphere, the geostrophic wind may not vary linearly with height. Thereby, the applicability of the simplified parameterizations under this condition needs more evaluation. As discussed in section 2.3, the simple model developed in this study is convenient to apply in numerical models. We will incorporate this scheme into numerical models and evaluate its performance under real conditions in future work.