Parameterization of Sheared Entrainment in a Well-Developed CBL. Part I: Evaluation of the Scheme through Large-Eddy Simulations

The development of the ABL over land is dominated by surface heating during the daytime. Convective activities and turbulent mixing are common within the daytime ABL, i.e., the CBL. The top of the CBL is capped by an interface layer, where the stratified free atmospheric air aloft is entrained down into the mixed layer by overshooting thermals. The entrainment process can significantly influence CBL evolution and profiles of mean variables within the CBL (Hoxit, 1974; Arya and Wyngaard, 1975; Lemone et al., 1999). It is important in NWP and air pollution models.

There have been many studies on the entrainment process. (Kim et al., 2006) pointed out that much work had been done for the free CBL, whereas studies about the uncertainties in the entrainment process were quite limited compared to studies of the sheared CBL. In the past decade, research on the entrainment process has mainly focused on how to understand and parameterize the effect of wind shear (Kim et al., 2003, 2006; Pino et al., 2003, 2006; Conzemius and Fedorovich, 2006a, 2006b, 2007; Sun and Xu, 2009). Bulk models are often employed to describe CBL evolution. Two commonly used types of bulk models are the zeroth-order model (ZOM), which represents the inversion layer (IL) as an infinitesimally thin interface (Tennekes, 1973), and the first-order model (FOM), which assumes a certain IL thickness (Betts, 1974). The IL structure is important for CBL dynamics (Sorbjan, 1996a, 1996b, 2004; Lewellen and Lewellen, 1998; Sullivan et al., 1998; vanZanten et al., 1999; Otte and Wyngaard, 2001; Kim et al., 2003, 2006; Sun and Wang, 2008). (Conzemius and Fedorovich, 2007) reviewed previous studies of bulk models and suggested that at least the FOM is needed in order to adequately capture the entrainment process in a sheared CBL. (Huang et al., 2011) demonstrated that the FOM can adequately describe not only the entrainment heat flux but also the entrainment fluxes of water vapor and other conservative scalars such as carbon dioxide.

The deficiencies of FOMs were reviewed in (Gentine et al., 2015). In FOMs based on (Betts, 1974), both the potential temperature and heat flux profiles are assumed linear in the IL, and the mixed-layer top is located with the minimum heat flux height. (Deardorff, 1979) argued that the representation of the IL in such models is oversimplified. Firstly, the observed maximum vertical gradient of potential temperature is generally much higher in observations than in the FOM. Secondly the minimum heat flux level is located above the mixed-layer top. Thirdly, the assumption that the mixed-layer height is equal to the minimum heat flux height generates a singularity for the IL growth rate equation under strong inversions. Thereafter, the IL thickness has generally been parameterized (e.g., defined as a function of the convective Richardson number) to avoid this singularity. On the other hand, the linear profiles are incompatible: parabolic flux profiles should correspond to linear profiles of conserved variables. (Deardorff, 1979) proposed a more realistic representation of the IL——the so-called general structure model (GSM). However, the structure of the IL in the GSM needs to be parameterized (Fedorovich and Mironov, 1995; Fedorovich et al., 2004a), which limits its applicability. To overcome the limitations in previous FOMs, (Gentine et al., 2015) proposed a new IL model based on a second-order polynomial for the potential temperature profile, and a third-order polynomial for the heat flux profile. This model can accurately prognosticate the growth rate of the IL, and of the mixed layer, under purely convective conditions. However, our study focuses on the entrainment process under shear conditions. It is not clear how wind shear impacts the profile shapes of velocity and potential temperature and their fluxes. In order to simplify the derivations, we use an FOM with linear profiles in the IL. The relative errors between the linear and curving flux profiles are discussed in this study.

The bulk model consists of a set equations for the CBL. Parameterizations of the entrainment flux ratio A _e (defined as the ratio of entrainment heat flux at the top of the mixed layer to the surface heat flux) and IL thickness are needed for closure of the CBL equations (e.g., Kim et al., 2006; Pino et al., 2006). (Kim et al., 2006) developed a parameterization of A _e for the CBL under the condition of height-constant geostrophic velocity (GC case). (Conzemius and Fedorovich, 2007) developed a bulk CBL model under the condition of linearly sheared geostrophic velocity with a zero value at the surface (GS case), in which the A _e is not explicitly expressed and the IL thickness is parameterized assuming a constant gradient Richardson number. Note that large uncertainties exist in the parameterization of sheared entrainment. (Pino et al., 2006) suggested that about 70% of TKE produced by wind shear across the IL is available for entrainment, whereas (Conzemius and Fedorovich, 2006a) proposed a value of 40%. (Sun and Xu, 2009) argued that the fraction should be 30%. Such a large discrepancy among different studies indicates that further investigation of the entrainment process is necessary for a better understanding of the CBL.

The IL thickness is another key parameter in bulk models. (Pino and Vilà-Guerau De Arellano, 2008) suggested that the IL thickness is a natural length scale that characterizes the shear-produced turbulence in the TKE budget at the CBL top. (Kim et al., 2006) proposed three schemes to estimate the IL thickness based on different empirical considerations of the effect of wind shear. (Conzemius and Fedorovich, 2007) developed a scheme under the assumption that the flux Richardson number remains at 0.25 in the entrainment zone (the layer where vertical potential temperature flux is negative). However, the LES results of (Pino and Vilà-Guerau De Arellano, 2008) showed the flux Richardson number to be larger than 0.25. Therefore, the scheme of (Conzemius and Fedorovich, 2007) needs further validation and the parameterization of the IL thickness should be modified.

In this paper, a parameterization scheme of A _e for a well-developed CBL is developed in an FOM framework. As in (Conzemius and Fedorovich, 2007), the CBL is assumed to develop under the condition of the GS case. The impacts of different factors on A _e are discussed. A new convective velocity scale for both buoyancy and shear effects is proposed to parameterize the IL thickness. LESs are conducted to evaluate the parameterization of A _e and the performance of parameterization for the IL thickness. In a companion paper, Part II, these parameterizations are further simplified according to the characteristics of entrainment derived from the LES output, and a simple model for predicting the growth rate of the well-developed CBL is proposed and evaluated.

Twenty-six CBL cases are simulated using an LES model to provide sufficient basic data in this study. The model used is DALES (the Dutch Atmospheric Large-Eddy Simulation model), which is based on the LES code of (Nieuwstadt and Brost, 1986) and developed by researchers from Delft University, the Royal Netherlands Meteorological Institute, Wageningen University and the Max Planck Institute for Meteorology (Heus et al., 2010). The domain size in this study is 10.0× 10.0× 4.0 km³, in the x, y and z directions respectively, with grid dimensions of 256× 256× 400. (Sullivan and Patton, 2011) pointed out that the lower-order moment statistics (means, variations and fluxes) become grid-independent when the ratio of CBL height to LES filter width is larger than 56. This value corresponds to their case D, in which the mesh resolution was thought to be fine enough to characterize the entrainment process. In the present study, the ratio of CBL height to LES filter width is about 31, which is slightly larger than that of case C in (Sullivan and Patton, 2011). Their results showed that the lower-order moment statistics change slightly while the entrainment rate is obviously overestimated in case C. However, their sensitivity experiments indicated that the finer vertical resolution can improve the LES estimates of entrainment rate efficiently. Our vertical mesh resolution is closer to that in case D than case C in (Sullivan and Patton, 2011). It is expected that the mesh resolution in this study will be able to describe the turbulence statistics and entrainment process reasonably. A third-order Runge-Kutta scheme with self-adaptive time stepping is used for time integration. The surface is treated as a semi-slip boundary at the bottom, and Monin-Obukhov similarity theory is applied at the lowest model level to calculate surface momentum flux. The top 1 km of the domain is a sponge layer and periodic boundary conditions are applied at the lateral boundaries. The closure scheme for the calculation of subgrid-scale fluxes is based on the TKE method (Deardorff, 1980).

For all cases in this study, the surface potential temperature flux is prescribed to be 0.1 K m s^-1, and the potential temperature at the surface is initially set to be 300 K. The Coriolis parameter f is set to be a constant value of 10^-4 s^-1. Half of the cases are conducted with a large gradient of potential temperature [γ_θ=0.006 K m^-1 (denoted as 6)], and the other half are conducted with a small gradient [γ_θ=0.003 K m^-1 (denoted as 3)]. Two cases are free-convection cases (denoted as NS00), and the others are divided into two groups, i.e., the GC and GS groups, with the geostrophic velocity along the x-direction. In the GC group, the geostrophic velocity is prescribed with three different values: 10 m s^-1, 15 m s^-1 and 20 m s^-1 (denoted as GC10, GC15 and GC20, respectively). In the GS group, the geostrophic velocity is zero at the surface and linearly increases with height at three different vertical gradients: 10 m s^-1, 15 m s^-1 and 20 m s^-1 per 2 km (denoted as GS10, GS15 and GS20, respectively). Two values of surface roughness lengths, z₀=0.1 m and z₀=0.01 m, are used to represent the rough surface (denoted as R) and the smooth surface (denoted as S). The case name GC15R3 means that the simulation is conducted under the conditions of a constant geostrophic velocity of 15 m s^-1, over a rough surface, with z₀=0.1 m, and an initial potential temperature gradient of 3 K km^-1. Results from the 26 cases are used to determine the empirical constants in the parameterization schemes.

The present study is based on linear equations of potential temperature and momentum for a horizontally homogeneous CBL. With Galilean transformation, CS cases can be easily transformed to GS cases. It is expected that the parameterizations derived from GS cases should be suitable for CS cases. However, this is not true for a nonlinear system such as the three-dimensional CBL. Furthermore, non-zero surface geostrophic velocity leads to changes in surface friction velocity and mixed-layer velocity, and consequently affects velocity and fluxes at the CBL top. The changes in these variables are not linear. For this reason, four additional CS CBL experiments (C5S10S3, C5S15S3, C5S15S6 and C5S15R3) are conducted to validate the parameterizations derived from GS cases. The case name C5S15S3 means that the surface value of geostrophic velocity is 5 m s^-1, the gradient of geostrophic velocity is 15 m s^-1 per 2 km, the surface is smooth (z₀=0.01 m), and the gradient of potential temperature is 3 K km^-1.

The integration for each simulation covers 28 800 s, and the results for the period from 4800 s to 28 800 s are output at a time interval of 100 s. The horizontally averaged idealized profiles of potential temperature, velocity and their fluxes can be obtained from the LES output. Figures 1 and 2 depict these LES cases and idealized profiles. In this paper, the mean and the fluctuating parts of a turbulent variable are denoted with uppercase and lowercase letters; for example, Θ and θ represent the domain averaged and fluctuating parts of potential temperature. A horizontally averaged turbulent flux is denoted by an overbar; for example, \(\overline{w\theta}\) means horizontally averaged kinematic heat flux. The idealized profile of Θ is assumed to vary linearly with height in all cases. (Fedorovich, 1995) gave the idealized profile of \(\overline{w\theta}\) in the mixed layer as \begin{equation} \label{eq1} \overline{w\theta}=\left(1-\dfrac{z}{h_1}\right)\overline{w\theta}_{\rm s}+\dfrac{z}{h_1}\overline{w\theta}_1 , (1)\end{equation}

where z is height and h₁ is the CBL height, which is defined as the height at which \(\overline{w\theta}\) from the LES output reaches its minimum; \(\overline{w\theta}_\rm s\) and \(\overline{w\theta}_1\) are the kinematic heat flux at the surface and at h₁ respectively. Based on this equation, A _e is expressed as \begin{equation} \label{eq2} A_{\rm e}=-\frac{\overline{w\theta}_1}{\overline{w\theta}_{\rm s}}=\frac{h_1-h_0}{h_0} , (2)\end{equation} where h₀ is the first zero-crossing height of the \(\overline{w\theta}\) profile. h₂ is defined as the level at which \(\overline{w\theta}\) is larger than 10% of its minimum value, which is the same as in (Pino et al., 2006) and (Conzemius and Fedorovich, 2006a). The layer from h₁ to h₂ is the so-called IL, and its thickness is ∆ h₂₁=h₂-h₁. The layer between h₀ and h₂ is the so-called entrainment zone, and its thickness is ∆ h₂₀=h₂ -h₀. It should be noted that the definition of the IL is based on the idealized profile of Θ , while the entrainment zone is based on the profile of \(\overline{w\theta}\). Following (Pino et al., 2006), Θ₁ (the potential temperature in the mixed layer) is determined from the LES Θ profile at the center of the CBL (h₁/2), Θ₂ is determined from the LES Θ profile towards h₂, and the potential temperature jump across the IL is ∆Θ=Θ₂-Θ₁. Following the approach of (Fedorovich, 1995), the profile of \(\overline{w\theta}\) in the IL is defined as a quadratic function of z if Θ linearly increases with height. Calculations show that the time-averaged relative error of the integral \(\overline{w\theta}\) in the IL has a maximum value of 7.8% in GC20S6. Because the errors introduced by the linear assumption are small in all of the cases, we prefer to use the linear \(\overline{w\theta}\) profile in this study. (Kim et al., 2006) pointed out that the linear profile of \(\overline{w\theta}\) gives a larger A _e than the curving profile. However, the errors of A _e are associated with the empirical constants in the parameterization scheme. They are obtained by a least squares fit to the LES outputs. The results show that the derived A _e parameterization can perform very well (details given in section 3).

Figure 1.  Idealized profiles (solid lines) of a CBL with constant geostrophic wind (GC). Top: horizontally averaged potential temperature $\Theta$ and its vertical flux $\overline{w\theta}$; middle: horizontally averaged $x$-component velocity $U$ and its vertical flux $\overline{uw}$; bottom: horizontally averaged $y$-component velocity $V$ and its vertical flux $\overline{vw}$. Thick dashed lines represent LES profiles，dash-dot lines represent $h_0$, $h_1$ and $h_2$. The vertical axis represents height above the surface.

Figure 2.  As in Fig. 1, but with linearly increasing geostrophic wind (GS).

The idealized profiles of the velocity components U and V are also assumed to be a linear function of z (Figs. 1 and 2). For the GC cases, U and V are constant in the mixed layer; thereby, U₁ and V₁ (velocity componets at h₁) are determined from the LES velocity profiles at the center of the CBL (h₁/2). For the GS and CS cases, the idealized U is constant in the mixed layer, whereas the idealized V linearly increases with height in the mixed layer. Thus, the determination of U₁ in the GS and CS cases is the same as in the GC cases. Values of V at the surface and in the middle of the CBL (h₁/2), which are denoted as V _s and V_1/2, are determined from the LES profile of V at 0.1h₁ and 0.5h₁, respectively. Then, V₁ is given as 2V_1/2-V _s. For all sheared cases, the idealized U and V are assumed to be linear functions of height in the IL. U₂ and V₂ are determined from the LES profiles of U and V towards h₂. The wind jumps across the IL are ∆ U=U₂-U₁ and ∆ V=V₂-V₁.

Looking at the idealized profiles of momentum fluxes, it is found that \(\overline{uw}\) and \(\overline{vw}\) below h₁ vary linearly with height in the GC cases. In the GS and CS cases, \(\overline{uw}\) below h₁ is a linear function of height, whereas \(\overline{vw}\) below h₁ is a quadratic function of z (Fedorovich, 1995). According to the idealized profiles, the momentum fluxes at h₁ are written as \begin{eqnarray} \overline{uw}_1&=&\dfrac{1}{1+{\Delta h_{21}}/{h_1}}\Bigg[\dfrac{\Delta h_{21}}{2h_1}\overline{uw}_{\rm s}-\left(\Delta U- \frac{1}{2}\gamma_u\Delta h_{21}\right)w_{\rm e}+ \frac{1}{4}f\Delta h_{21}(V_{\rm s}-V_1)\Bigg]\ ({\rm for\ all\ cases}) ,\nonumber\\[1.2mm] \overline{vw}_1&=&\dfrac{1}{1+{\Delta h_{21}}/{h_1}}\left(\dfrac{\Delta h_{21}}{2h_1}\overline{vw}_{\rm s}-w_{\rm e}\Delta V\right)\ ({\rm for\ GC\ cases}) ,\quad (3)\\[-2mm]\nonumber \end{eqnarray} and \begin{eqnarray} \overline{vw}_1&=&\dfrac{1}{1+{\Delta h_{21}}/{h_1}}\Bigg\{ \dfrac{\Delta h_{21}}{h_1}\overline{vw}_{\rm s}-\dfrac{1}{2}\Delta h_{21}\dfrac{\partial V_{\rm s}}{\partial t}-\nonumber\\[1.2mm] &&\left(\Delta V+\dfrac{V_{\rm s}-V_1}{2h_1}\Delta h_{21}\right)w_{\rm e}- \dfrac{1}{2}f\Delta h_{21}(U_{\rm s}-U_{{\rm g,s}}) \Bigg\}\ ({\rm for\ GS\ and\ CS\ cases}) . (4)\end{eqnarray} The above equations are derived by vertically integrating the momentum equations (the derivations are given in Appendix A), where U _g is geostrophic velocity; U _g,s is the surface geostrophic velocity; γ_u is the gradient of U _g, w _e is the entrainment rate, which is defined as \(w_\rm e=\partial\langle h_1\rangle/\partial t\), \(\langle h_1\rangle\) is the least squares fit of h₁, according to the relation \(h_1\propto \sqrt t\), to avoid a negative value of w _e. (Fedorovich et al., 2004a) showed that h₁ follows this relation in the NS case. Our LES outputs show that this relation is also effective in the sheared CBL cases (figures presented in the companion paper, Part II). Similarly, the departure of the idealized \(\overline{uw}\) and \(\overline{vw}\) profiles from the curving ones in the IL also introduces some errors. Calculations show that the maximum relative error of \(\overline{uw}\) is 8% in the sheared CBL cases. The maximum relative errors of \(\overline{vw}\) are 10.6%, 658% and 371% in the GC, GS and CS cases, respectively. However, the LES outputs show that in the GS and CS cases, the contribution of IL wind shear to the entrainment is negligibly small in the y-direction when compared with that in the x-direction (results reported in the companion paper, Part II). Therefore, the idealized profiles can characterize the IL shear effect on the entrainment reasonably.

The derivation begins with the TKE (E) budget in Boussinesq approximation. It is expressed as \begin{eqnarray} \label{eq3} \dfrac{\partial E}{\partial t}&=&\dfrac{g}{\Theta_0}\overline{w\theta} -\left(\overline{uw}\dfrac{\partial U}{\partial z}+\overline{vw}\dfrac{\partial V}{\partial z}\right)- \left(\dfrac{\partial\overline{wE}}{\partial z}+\dfrac{1}{\rho_0}\dfrac{\partial\overline{wp}}{\partial z}\right)-\varepsilon ,\quad\ (5)\end{eqnarray} where ρ₀ is the air density (Moeng and Sullivan, 1994), p is the fluctuating part of pressure, g represents the acceleration of gravity, ε is the viscous dissipation rate of TKE. The TKE storage term on the left-hand side is small compared to other terms, except in the early stage of CBL evolution (Driedonks, 1982; Randall, 1984). The LES outputs indicate that after 4800s of integration, when the CBL is well-developed, this term is small and can be neglected. The first and second terms on the right-hand side are the production rates of TKE by buoyancy and wind shear, respectively. The third term is the vertical transport rate of TKE. This is a divergence term, and thus its integration from the surface to h₂ should be zero (Moeng and Sullivan, 1994). ε is usually assumed to be proportional to its production rate (Flamant et al., 1999; Conzemius and Fedorovich, 2006b; Kim et al., 2006). Applying idealized profiles of potential temperature, velocity and their fluxes, as shown in Fig. 2, the vertical integration of the TKE budget can be written as (see derivation in Appendix B) \begin{eqnarray} &&-\dfrac{1}{2}\dfrac{g}{\Theta_0}\overline{w\theta}_1(h_1+\Delta h_{21})\nonumber\\ &&=\dfrac{1}{2}(1\!-\!\alpha_1)\dfrac{g}{\Theta_0}\overline{w\theta}_{\rm s}h_1+(1-\alpha_2)C_{\rm D}^{-1/2}u_\ast^3+\nonumber\\ &&\quad(1-\alpha_3)\left(-\dfrac{1}{2}\overline{uw}_1\Delta U-\dfrac{1}{2}\overline{vw}_1\Delta V\right)+\nonumber\\ &&\quad(1-\alpha_4)(V_1-V_{\rm s})\Bigg(-\dfrac{1}{2}\overline{vw}_{\rm s}-\dfrac{1}{2}\overline{vw}_1+\dfrac{1}{12}f\gamma_uh_1^2\Bigg) .\qquad \label{eq4} (6)\end{eqnarray} where the coefficients α₁,α₂,α₃, and α₄ are the proportions of the dissipation rate to the corresponding production rate and do not vary with height, $(u_\ast=\sqrt[4]{\overline{uw}_\rm s^2+\overline{vw}_\rm s^2})$ is the friction velocity, and C _D=u_{*²/(U s²+V s²)} is the surface drag coefficient. Together with the definition of the convective velocity scale, i.e., \(w_\ast^3=(g/\Theta_0)\overline{w\theta}_\rm sh_1\), the above equation yields the parameterization of A _e, which is expressed as $$ 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_*}, (7)$$ $$\ \ \ \ \ \ \ \ \ \ Term IV$$ where A₁=1-α₁, A₂=2(1-α₂), A₃=2(1-α₃), and A₄=2(1-α₄). The terms on the right-hand side represent the contributions of the buoyancy (Term I), the surface layer shear (Term II), the IL shear (Term III) and the mixed layer shear (Term IV), respectively. Compared with the parameterization scheme for a GC case in (Kim et al., 2006), the only difference is that Eq. (7) has an additional term, i.e., Term IV. If the geostrophic velocity gradient vanishes, the GS case will become a GC case (that is, V _s=V₁), and Eq. (7) will turn out to be the same as the parameterization scheme for the GC case described in (Kim et al., 2006) [see their Eq. (22); the term \(-\overline{uw}_1\Delta U/2-\overline{vw}_1\Delta V/2\) is equivalent to that expressed in their Eq. (5B)].

By applying stepwise regression to outputs from all of the NS, GC and GS cases, the coefficients in Eq. (7) are given as A₁=0.21, A₂=0.01, A₃=0.86 and A₄=0.70. The value of A₁ is very close to the classical value of 0.2 (Stull, 1976; Fedorovich et al., 2004a, 2004b). From the definition of A₂, it can be easily shown that α₂=99.5%, which means that surface shear-produced TKE dissipates locally (Conzemius and Fedorovich, 2006a; Pino and Vilà-Guerau De Arellano, 2008). A₃=0.86 means that the fraction of the shear-generated TKE used for the entrainment process is 43%, which is approximately the same as in (Conzemius and Fedorovich, 2006a) and supports the argument of (Sun and Xu, 2009) that the value of 1.44 for A₃ proposed in (Pino et al., 2006) is overestimated. The stepwise regression also shows relatively large uncertainties in the determination of A₄. However, the fourth term is too small to significantly influence the accuracy of Eq. (7) (see the results in the next section).

Figure 3.  Horizontally and 30-min averaged vertical profiles of the total (upper panels) and the subgrid (lower panels) TKE budget for the GC15S3, GS15S3 and C5S10S3 cases. S, B, T $\varepsilon $ and $\partial E/\partial t$ represent shear production, buoyancy production, transport and dissipation rates of TKE, respectively. $\varepsilon_\rm p$ represents a linear combination of the shear and buoyancy production. Thin grey lines from bottom to top in each panel represent $h_0/h_1$, 1 and $h_2/h_1$, respectively.

In the above derivations, ε is treated as the sum of the dissipation rates of buoyancy-and shear-produced TKE. The regression of Eq. (6) to the LES outputs shows that the dissipation rate of shear-produced TKE varies in different parts of the CBL. That is, the parameterized dissipation rate ε _p can be calculated as

\begin{eqnarray} \label{eq6} \varepsilon_{\rm p}&=&-\alpha_1\dfrac{g}{\Theta_0}\overline{w\theta}_{\rm s}\left(1-\dfrac{z}{h_1}\right)-\alpha_xS ,(8)\\ \label{eq7} \alpha_1&=&0.79 ,\nonumber\\ \alpha_x&=&\left\{ \begin{array}{l@{\quad}l} \alpha_2=0.995 & (z\le 0.1h_1)\\[1mm] \alpha_3=0.57 & (z\ge h_1)\\[1mm] \alpha_4=0.65 & (0.1h_1 \le z\le h_1) \end{array} \right., (9)\end{eqnarray} where S is the shear production rate of TKE. Figure 3 depicts the profiles of ε _p and the forcing terms on the right-hand side of the TKE budget from the LES cases with weak inversion. ε _p is very close to the dissipation rate ε calculated from the LES outputs, suggesting that the coefficients in Eq. (7) are reasonable. The subgrid TKE budgets in these cases are also illustrated in Fig. 3. The results indicate that the subgrid TKE is negligibly small in the IL. The cases with strong inversion have the same situation [Fig. S1 in electronic supplementary material (ESM)]. Therefore, the resolved motions dominate the TKE budget in the IL and the derived parameterizations based on the LES outputs are reasonable. Figure 4 shows that the A _e estimated by Eq. (7) agrees well with that derived from the LES outputs. As presented in previous studies, the value of A _e calculated from LES outputs fluctuates significantly because of the fluctuation of 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 value of the parameterized A _e is contained within the fluctuations of the LES outputs. Therefore, the parameterization expressed as Eq. (7) can capture the characteristics of entrainment flux for a well-developed sheared CBL.

Figure 4.  Entrainment heat flux ratios in CS cases from LES outputs (blue dots) and calculated by Eq. (6) (red dots).

Figure 5.  Each term on the right-hand side of Eq. (7) for the parameterization of $A_\rm e$ in the GC cases. The blue dots represent Term I, the green dots represent Term II, and the red dots represent Term III.

Figure 6.  As in Fig. 5, but for the results in the GS cases. The cyan dots represent Term IV. Note that green dots and cyan dots overlap in some cases.

The evolution of each term in Eq. (7) during CBL development is illustrated in Fig. 5 for the GC case, and in Fig. 6 for the GS case. It is clear that Term I and Term III are the dominant terms for CBL development. The stratification and wind shear have little influence on Term I. Its value is about 0.18 and remains almost unchanged throughout CBL development in all of the simulation cases. Term II is always very small, albeit its value differs among cases. This result suggests that Term II has little influence on entrainment. The behavior of Term III is quite different in the GC and GS cases. Term III decreases with time and increases with stratification in the GC case, whereas in the GS case it remains almost constant throughout CBL development and decreases with stratification. Term IV exists in both the GS and CS cases. Figure 6 shows that this term is negligibly small in the early stage of a developing CBL. However, it increases during CBL development. Under the condition of a weak geostrophic velocity gradient, such as in GS10, Term IV can be neglected, since it remains very small throughout the entire CBL development process. On the other hand, if the geostrophic velocity gradient and h₁ are sufficiently large, Term IV becomes relatively large. For example, in GS20R3, at the end of the simulation, when h₁ is about 2000 m, the contribution of Term IV to the A _e is about 17%. However, this situation seldom happens in the real atmosphere because it is difficult for a geostrophic velocity gradient as large as that shown in case GS20R3 to occur.

In the GC case, the velocity jump across the IL increases slightly with time, but the momentum flux at h₁ decreases with time. Thus, their total shear production rate of TKE has a slight decreasing trend (Fig. S2 in ESM). Meanwhile, the denominator of Term III (i.e., w_*³) increases remarkably during this process. This is why Term III decreases with increasing CBL depth. The LES outputs indicate that a larger gradient of potential temperature can significantly enlarge the velocity jump across the IL and slightly decrease the momentum flux at h₁ (Fig. S3 in ESM). Thereby, the shear-produced TKE and A _e enlarge under stronger stratification. However, this does not mean that the growth rate of the CBL under stronger background stratification increases, since the capping inversion strength also enhances, which suppresses the CBL's development (Sun and Xu, 2009). The effect of the rough surface is to enlarge the value of Term III. This is because, under such a condition, the velocity in the mixed layer is smaller and the velocity jump at the CBL top is larger, as compared to under a smooth surface condition.

In the GS case, the value of the velocity jump across the IL increases while the momentum fluxes remain almost constant with time; thus, the shear production of TKE in the IL enhances during the CBL's development (Fig. S4 in ESM). Meanwhile, the denominator of Term III increases steadily during this process. Thus, the value of Term III does not change significantly. This implies that the shear production of TKE (i.e. \(-(\overline{uw}_1\Delta U+\overline{vw}_1\Delta V)/2\)) is approximately proportional to w_*³. The reduction effect of strong stratification on Term III means that the shear production of TKE may be proportional to the inverse of γ_θ. A larger geostrophic velocity gradient leads to a larger momentum flux at h₁ and velocity jump across the IL, and consequently a larger value of Term III. Figure 6 also shows that Term III is almost not influenced by surface roughness.

Based on the above results, it can be deduced that the shear production rate of TKE at the CBL top can be divided into two parts. One part is proportional to w_*³, γ_u and 1/γ_θ; the other part is insensitive to CBL development and γ_θ but sensitive to surface roughness. The former part dominates in the GS case, whereas the latter works only in the GC case. In the CS case, these two parts cooperate, but the former still dominates. Thus, Term III shows a slight decreasing trend and becomes weak with larger γ_θ (Fig. S5 in ESM). The expressions and meaning of these two parts are discussed in detail in the companion paper, Part II.

Equation (7) can be rewritten as \begin{equation} \label{eq8} A_{\rm e}=\dfrac{A_1}{1+{\Delta h_{21}}/{h_1}}\frac{w_{\rm m}^3}{w_\ast^3} ,(10) \end{equation} and \begin{eqnarray} w_{\rm m}^3&=&w_\ast^3+\frac{A_2}{A_1}C_{\rm D}^{-1/2}u_\ast^3+\dfrac{A_3}{A_1}\left(-\dfrac{1}{2}\overline{uw}_1\Delta U-\dfrac{1}{2}\overline{vw}_1\Delta V\right)+\nonumber \end{eqnarray} \begin{eqnarray} &&\dfrac{A_4}{A_1}(V_1-V_{\rm s})\left(-\dfrac{1}{2}\overline{vw}_{\rm s}-\dfrac{1}{2}\overline{vw}_1+\dfrac{1}{12}f\gamma_uh_1^2\right) , \label{eq9} (11)\end{eqnarray} where w _m can be interpreted as a new characteristic convective velocity scale that includes the contributions from both the buoyancy and the wind shears in a CBL. The results in Figs. 5 and 6 suggest that the characteristic convective velocity is mainly enhanced by the IL shear. In the ZOM, Eq. (9) reduces to A _e=A₁w _m³/w_*³, which agrees with the result of (Tennekes, 1973), (Driedonks, 1982), and (Moeng and Sullivan, 1994) that the A _e can be expressed as 0.2w _m³/w_*³.

For the GC case, the simplified form of Eq. (11) is often used to characterize the convective velocity scale in a sheared CBL, which includes only w_* and u_* on the right-hand side of the equation (Tennekes, 1973; Zeman and Tennekes, 1977; Tennekes and Driedonks, 1981; Driedonks, 1982; Boers et al., 1984; Batchvarova and Gryning, 1994; Moeng and Sullivan, 1994; Pino et al., 2003). For example, (Tennekes, 1973) suggested that w _m³=w_{*³ +12.5u*³}, while (Moeng and Sullivan, 1994) proposed that w _m³=w_*³+5u*³. Note that the equation w _m³=w_*³+Bu*³ only includes the contribution of shear-produced TKE in the surface layer. Actually, it can be regarded as the simplified form of Eq. (11) by assuming that \(-\overline{uw}_1\Delta U/2-\overline{vw}_1\Delta V/2\) is approximately proportional to u_*³ [the last term of Eq. (11) is zero under the GC condition). As mentioned in the previous section, this term is insensitive to CBL development and stratification strength but sensitive to surface roughness. Our LES outputs also show that the result of (Moeng and Sullivan, 1994) is a good estimate of w _m. However, for the GS and CS cases, the last term on the right-hand side of Eq. (11) is relatively small and can be neglected (although it is not zero), but the third term on the right-hand side of Eq. (11) is closely related to w_*³ , γ_u and 1/γ_θ. In this situation, the simplified form w _m³=w_*³+Bu*³ is not a good approximation of Eq. (11).

In the FOM, the IL thickness (∆ h₂₁=h₂-h₁) is a key parameter that is often used in the mixed-layer model, as described in (Pino et al., 2006) and (Conzemius and Fedorovich, 2007). According to parcel theory, after the overshooting thermal rises across the IL, its kinetic energy is transformed to potential energy. That is to say, w _m²∝(g/Θ₀)∆Θ∆ h₂₁. Based on this assumption, (Kim et al., 2006) gave the parameterization of the IL thickness in the form of \begin{equation} \label{eq10} \frac{\Delta h_{21}}{h_1}=aRi^{-1}+b , (12)\end{equation} and \begin{equation} \label{eq11} Ri=\dfrac{g}{\Theta_0}\Delta\Theta h_1/{w_{\rm m}^2} , (13)\end{equation} where a and b are empirical constants, and Ri is the convective Richardson number. (Kim et al., 2006) proposed an empirical formula to characterize the turbulence velocity scale under the GC condition, expressed as \begin{equation} \label{eq12} w_{\rm m}^2=w_\ast^2+cu_\ast^2+d(\Delta U^2+\Delta V^2) , (14)\end{equation} where c and d are empirical constants. (Kim et al., 2006) provided three groups of these empirical constants. (Pino et al., 2006) used this scheme in a mixed-layer model to evaluate their parameterization of A _e. By applying stepwise regression to the outputs of all the NS, GC and GS cases, cu_*² in Eq. (14) is excluded (because the existence of this term makes the significance of regression reduced), and the constants are a=0.37, b=0.13 and d=0.19. When the constant b is constrained to be zero, the regression also excludes cu_*² in Eq. (14), and the constants are a=2.46 and d=0.056. The exclusion of cu_*² suggests that the wind shear in the surface layer has little effect on entrainment, which agrees with the result in section 3.1 that the surface shear-produced TKE dissipates locally. This scheme with different constants is denoted as KP1 and KP2, respectively (Table 1).

In the previous section, a new convective velocity scale is proposed. We expect that it is appropriate for the estimation of IL thickness, and thus we use Eq. (12) but replace Eq. (14) with Eq. (11). By a least squares fit to the LES outputs of NS, GC and GS cases, the empirical constants a and b are determined to be 0.70 and 0.14, respectively. This scheme is denoted as LS1 (Table 1).

The parameterization of the IL thickness is usually evaluated by field observations (e.g., Boers et al., 1984), experiments (e.g., Deardorff et al., 1980) and LESs (e.g., Fedorovich et al., 2004a). Uncertainties in h₂ and Θ₂ determined from LES outputs will subsequently result in biases in the calculation of ∆ h₂₁ and ∆Θ. A positive bias of ∆ h₂₁ (as well as ∆Θ) causes a negative bias of Ri^-1, resulting in low correlation between ∆ h₂₁/h₁ and Ri^-1 (Sun et al., 2005). For this reason, another Richardson number, Ri_N, which is based on the buoyancy frequency \(\sqrt{(g/\Theta_0)\gamma_\theta}\) in the free atmosphere, is used in (Fedorovich et al., 2004a). As proposed by (Stull, 1973), the time taken for the rising thermal to penetrate into the free atmosphere should be related to the buoyancy frequency. That is, \(\Delta h_21/w_\rm m\propto 1/\sqrt{(g/\Theta_0)\gamma_\theta}\). The parameterization scheme can be written as \begin{equation} \label{eq13} \frac{\Delta h_{21}}{h_1}=a_N Ri_N^{-\frac{1}{2}}+b_N , (15)\end{equation} and \begin{equation} \label{eq14} Ri_{N}=\frac{g}{\Theta_0}\frac{\gamma_\theta h_1^2}{w_{\rm m}^2} , (16)\end{equation} where a_N and b_N are empirical constants. For comparison purposes, we use Eq. (15) to parameterize the IL thickness and employ Eq. (11) to characterize the turbulent velocity in Eq. (16). The least squares fit to our LES outputs of NS, GC and GS cases yields a_N=1.30 and b_N=0.00. This scheme is denoted by LS2 (Table 1). It should be pointed out that if ∆Θ∝γ_{θ∆ h21}, Eqs. (12) and (15) should be identical. The LES outputs indicate that ∆Θ≈1.88γ_θ∆ h₂₁ (Fig. S6 in ESM). The LS1 and LS2 schemes are actually equivalent; we denote them as LS1 and LS2 simply because they use different variables and have different constants.

(Mahrt and Lenschow, 1976) and (Conzemius and Fedorovich, 2006a) suggested that a balance exists between the shear production and buoyancy destruction of TKE in the entrainment zone, which can be described by the flux Richardson number or gradient Richardson number. (Conzemius and Fedorovich, 2007) set the bulk gradient Richardson number (Ri _b) in the IL to be a critical value of 0.15. The IL thickness is then given by \begin{equation} \label{eq15} \Delta h_{21}=Ri_{\rm b}\frac{\Delta U^2+\Delta V^2}{\frac{g}{\Theta _0}\Delta \Theta} . (17)\end{equation} The constant Ri _b is found to be 0.23 by a least squares fit to our LES outputs of GC and GS cases. This scheme is denoted as CF (Table 1).

In order to evaluate the performance of the above five parameterization schemes for the IL thickness, the relative errors are calculated and illustrated in Fig. 7. The relative error (Err) is defined as \begin{equation} \label{eq16} {\rm Err}=\dfrac{1}{n}\sum\left|\dfrac{\Delta h_{21{\rm ,p}}}{\langle{\Delta h_{{\rm 21,LES}}}\rangle}-1\right| , (18)\end{equation} where ∆ h_{21 ,p} is the inversion layer thickness predicted by each parameterization scheme. As mentioned above, ∆ h₂₁ is determined from the instantaneous LES profile of \(\overline{w\theta}\). This method can result in large errors that completely conceal differences between different parameterizations. In order to reduce errors, an equal weighted nine-point moving average operator is applied to ∆ h₂₁, and the result is denoted by \(\langle\Delta h_\rm 21,LES\rangle\). The CF scheme has the largest error in most cases. Further analysis shows that the bulk gradient Richardson number varies from 0.17 to 2.31 in different cases, whereas only in a few cases is the bulk gradient Richardson number very close to 0.23. This is why the CF scheme has relatively large errors in most cases. The two KP schemes apply well, although KP2 has slightly lager errors than KP1. The new empirical constants significantly improve the performance of the KP schemes (Fig. S7 in ESM), implying that the original ones proposed by (Kim et al., 2006) and (Pino et al., 2006) are not very representative because they were derived from only a few LES cases. LS1 is not improved in comparison with the two KP schemes; and the reason is probably that they all use the convective Richardson number (Ri), as discussed previously. However, LS2 has the best performance, and the errors are less than 20% in all cases, suggesting that the Ri_N scheme is more suitable for characterizing the IL thickness.

Figure 7.  Relative errors of the predicted capping IL thickness against the LES results. The parameterization schemes \small \parbox[t]14cmare listed in Table 1.

The most recent study on shear-free entrainment by means of direct numerical simulation (Garcia and Mellado, 2014) suggests a two-layer model might be appropriate for studying the entrainment zone. The upper sub-layer thickness of the entrainment zone (δ) is defined based on the maximum potential temperature gradient in (Garcia and Mellado, 2014) [see Fig. 5 and Eq. (20) in their paper]. It is different to the IL thickness (∆ h₂₁) defined in this study. The bottom of δ is located at the level of the maximum gradient of potential temperature that is higher than the bottom of ∆ h₂₁, while the top of δ is lower than the top of ∆ h₂₁. Thus, the upper sub-layer of the entrainment zone defined in (Garcia and Mellado, 2014) is part of the IL defined in this study. Their results show that the upper sub-layer thickness (δ) of the entrainment zone is actually the mean penetration depth of an overshooting thermal, which is directly affected by the background stratification N² \(N=\sqrt{(g/\Theta_0)\gamma_\theta}\). The following relation is obtained [their Eq. (24)]: \begin{equation} \label{eq17} \delta=c_\delta({w_\ast}/N) , (19)\end{equation} where c_δ=0.55 is the coefficient obtained from the direct numerical simulation results. In fact, the LS2 scheme is equivalent to ∆ h₂₁=a_N(w _m/N). Following the method in (Garcia and Mellado, 2014), δ is determined from our LES Θ profile. The LES results show that ∆ h₂₁/δ=2.44, which is quite close to 2.36, the ratio of a_N to c_δ. This implies that the LS2 scheme is similar to Eq. (18) in a shear-free CBL, because both ∆ h₂₁ and δ are the overshooting distances of thermals rising in the stably stratified environment. The difference between the LS2 scheme and Eq. (18) is attributed to different definations of ∆ h₂₁ and δ. Thus, in a shear-free CBL, the IL thickness is dominated by overshooting thermals. However, in a sheared CBL, the effect of wind shear on the IL thickness is also important. Our results suggest that w _m is suitable for characterizing the joint effects of thermal overshooting and wind shear on IL thickness.

In an FOM framework, the parameterization of A _e at the top of a well-developed CBL under the GS condition is derived by vertically integrating the TKE budget. Compared to the parameterization scheme under the GC condition proposed by (Kim et al., 2006) and (Pino et al., 2006), our scheme includes an additional term that represents the contribution of shear-produced TKE in the mixed layer. When the geostrophic velocity gradient becomes zero, the parameterization scheme turns out to be the one under the GC condition. This scheme is also valid for the CS case. Thus, the new parameterization developed in the present study is appropriate for entrainment approximation in a well-developed CBL under different linearly sheared geostrophic velocity conditions.

The new parameterization contains four terms representing the effects of the buoyancy, surface layer shear, IL shear and mixed layer shear, respectively. The buoyancy and IL shear are the dominant terms among these four terms. The LES outputs indicate that the shear-produced TKE in the surface layer dissipates locally, and 43% of the shear-produced TKE at the CBL top contributes to the entrainment, which is approximately the same as the results in (Conzemius and Fedorovich, 2006a).

A new convective velocity scale in the sheared CBL is proposed. It includes the contributions of buoyancy and wind shears. In the GC cases, the convective velocity scale is equivalent to the simplified form proposed by (Moeng and Sullivan, 1994), in which the effect of wind shear in the entire CBL can be approximately represented by the friction velocity. LES outputs show that the direct contribution of surface shear to the entrainment is relatively small. However, as pointed out by (Conzemius and Fedorovich, 2006a), the surface shear has an indirect effect on entrainment by slowing the flow in the CBL interior and inducing shear at the CBL top, and it is the IL shear that enhances the entrainment. Apparently, the contribution of IL shear to the entrainment process has been considered in the simplified formula by the friction velocity. However, note that in the GS and CS cases, the simplified form of the convective velocity scale is not valid because the shear-produced TKE in the IL is mainly related to w_*, the geostrophic velocity gradient and stratification strength, rather than the friction velocity.

The parameterization schemes of the IL thickness proposed in previous studies are evaluated by the LES outputs. The schemes suggested by (Kim et al., 2006) apply well when the new empirical constants are used. The empirical constants are derived by the stepwise regression to our LES outputs, which excludes the term representing the surface shear. This result supports that buoyancy and IL shear are the dominant factors of sheared entrainment. The parameterization scheme proposed by (Conzemius and Fedorovich, 2007) can only perform well in a few cases, because the bulk Richardson number varies widely in different cases. However, the buoyancy Richardson number approach (the Ri_N scheme), combined with the new convective velocity scale, can characterize the IL thickness well in all cases.

The A _e and IL thickness are important parameters in the mixed layer model. Our aim is to obtain a simplified scheme that can predict the developing process of a sheared CBL well. The parameterization scheme developed in this study represents our initial efforts to achieve this goal. The simplified model is further explored and discussed in the companion paper, Part II.

Finally, it is worth noting that the parameterizations proposed in this study may only be applicable for CBLs under special conditions. In the derivations we neglect the storage term in the TKE budget and only consider the linearly sheared geostrophic velocity and stable background stratification. However, when the CBL is in its early developing stage, the storage term in the TKE budget is not negligibly small, and the entrainment process may exhibit different characteristics. There often exists a residual layer in the real atmosphere. When the CBL is growing through this layer, the entrainment process is different to that of the stratified free atmosphere above the CBL. In addition, in the real atmosphere the geostrophic wind may not vary linearly with height. The applicability of the parameterization schemes under these conditions is not well investigated and needs further evaluation. These problems will be investigated in future work.