Bulk Scaling Model of Entrainment Zone Thickness in a Convective Boundary Layer, with a Shear Effect Promoted by Velocity Difference

Abstract

Studying the thickness of the convective boundary layer (CBL) is helpful for understanding atmospheric structure and the diffusion of air pollutants. When there is velocity shear in CBL, the flow field structure is very different from that of shear-free CBL, which makes the thickness model of the entrainment zone deviate. A large-eddy simulation (LES) approach is carried out for a horizontally homogeneous, atmospheric CBL, with a shear effect promoted by velocity difference to explore the bulk scaling model of the entrainment zone thickness. The post-processed data indicate that the existing bulk scaling models cannot synthetically represent the effects of shear and buoyancy on entrainment, resulting in reduced accuracy or limited applicability. Based on the fraction of turbulent kinetic energy (TKE) used for entrainment, a different form of the characteristic velocity scale, which includes the shear effect, is proposed, and a modified bulk scaling model that uses a potential temperature gradient to replace the potential temperature jump across the entrainment zone is constructed with the numerical results. The new model is found to provide an improved prediction of the entrainment zone thickness in a sheared CBL.

1. Introduction

Due to the compound influence of surface radiation, potential temperature gradient (also named inversion strength), wind shear, and other factors, the flow field in a sheared convective boundary layer (CBL) is complex; this is especially true in the entrainment zone formed by shear and buoyancy between the mixed zone and the free atmosphere, which restricts the exchange of matter and heat between the layers. Accurately describing and quantifying the thickness of the entrainment zone in a sheared CBL is of significance for meteorological prediction [1,2], as well as numerous environmental and engineering applications [3,4].

The current bulk scaling models for the entrainment zone are primarily based on the CBL framework of a “zero-order model” (ZOM) [5] and a “first-order model” (FOM) [6]. An FOM assumes a linear variation of parameters with height in the CBL, which is more consistent with the observed data. The height of the entrainment zone, $z_i$, as sketched in Figure 1, is defined by the level of the buoyancy flux minimum ${\overline{w^{\prime }{\theta }^{\prime }}|}_i$, and the upper limit of the entrainment zone, $z_2$, is the zero-crossing height of the buoyancy flux profile, using linear interpolation between discrete grid points. The height of the mixed zone, $z_1$, takes the corresponding level of zero buoyancy flux. Due to the large-scale turbulent eddies in the entrainment zone, it is difficult to accurately locate the height of zero buoyancy flux in the upper layer, and Fedorovich et al. [7] and Conzemius and Fedorovich [3,4] adopted the height at which the buoyancy flux returns 90% of the way back to zero from its grid minimum as $z_2$. Similarly, in order to eliminate the influence of the inhomogeneity of thermals in the surface layer, the corresponding position of 5~10% of the surface buoyancy flux is taken as the height of the underlying surface, $z_S$.

In a shear-free CBL, Turner and Benton [8] and Vanzanten et al. [9] concluded that the development of the CBL mainly depends on the turbulence process and the potential temperature change in the entrainment zone. Stull [10] and Zeman and Tennekes [11] established a bulk scaling model for the entrainment zone thickness, based on the conservation of kinetic and potential energy of thermals in the zone, as follows:

$$\mathsf{\Delta }z/z_i\propto R{i}_{*}^{-1},$$

(1)

where the convective Richardson number $R{i}_{*}$ is defined as

$$R{i}_{*}=\frac{g}{{\theta }_0}\mathsf{\Delta }\theta \cdot z_i/w_{*}^2,$$

(2)

and the convective velocity scale $w_{*}$ is given by

$$w_{*}={\left(\raisebox{1ex}{$g$}\!\left/ \!\raisebox{-1ex}{${\theta }_0$}\right.\overline{{w^{\prime }{\theta }^{\prime }|}_s}{z}_i\right)}^{1/3},$$

(3)

in which ${\theta }_0$ is the ambient virtual potential temperature and g is the gravitational acceleration.

With their convective water tank experiments, Deardorff et al. [12] obtained the following correlation for the above bulk scaling model:

$$\mathsf{\Delta }z/z_i=1.31R{i}_{*}^{-1}+0.21.$$

(4)

This correlation is denoted as Model 1 in this study. However, there is no consensus on such a bulk scaling model. For example, Boers [13] obtained a (−1/2) power-law correlation, while Beyrich and Gryning [14] concluded that there is no such correlation through their analysis of lidar measured data.

Gryning and Batchvarova [15] adopted the entrainment velocity $w_e$ ($w_e=d{z}_i/dt$) instead of $w_{*}$ as the velocity scale, and constructed the following bulk scaling model in terms of the entrainment Richardson number $R{i}_E$:

$$\mathsf{\Delta }z/z_i=3.3R{i}_E^{-1/3}+0.2,$$

(5)

where $R{i}_E$ is defined as

$$R{i}_E=\frac{g}{{\theta }_0}\mathsf{\Delta }\theta \cdot z_i/w_e^2.$$

(6)

This model is denoted as Model 2 in this study. Model 2 was found to agree well with the measured atmospheric data of Boers and Eloranta [16].

When there is velocity difference at the entrainment zone, a significant shear effect will enhance the entrainment process. For sheared CBL, Kim et al. [17] suggested that the bulk scaling model for $\mathsf{\Delta }z/z_i$ should include the influence of the shear effect, and a modified convective Richardson number defined below, which takes into account the shear effect, should be used instead:

$$R{i}_{*}=\frac{g}{{\theta }_0}\mathsf{\Delta }\theta \cdot z_i/w_m^2,$$

(7)

in which a new velocity scale is given by

$$w_m^2=w_{*}^2+a{u}_{*}^2+b\left(\mathsf{\Delta }{U}^2+\mathsf{\Delta }{V}^2\right),$$

(8)

where a and b are constants, $\mathsf{\Delta }U$ and $\mathsf{\Delta }V$ are the jumps of the x- and y-component velocity across the entrainment zone, and the friction velocity scale is

$$u_{*}=\sqrt{-{\overline{u^{\prime }{w}^{\prime }}|}_s}.$$

(9)

The modified $R{i}_{*}$, using $w_m$ in place of $w_{*}$ as the velocity scale, takes full account of the buoyancy and shear effects on entrainment. With such a modified convective Richardson number, Kim et al. [17] gave the following three bulk scaling models:

$$\mathsf{\Delta }z/z_i=0.075R{i}_{*}^{-1}+0.0925,\mathrm{with}{w}_m^2=w_{*}^2+180.21u_{*}^2+0.15\left(\mathsf{\Delta }{U}^2+\mathsf{\Delta }{V}^2\right);$$

(10)

$$\mathsf{\Delta }z/z_i=2.68R{i}_{*}^{-1},\mathrm{with}{w}_m^2=w_{*}^2+2.99u_{*}^2+0.03\left(\mathsf{\Delta }{U}^2+\mathsf{\Delta }{V}^2\right);$$

(11)

$$\mathsf{\Delta }z/z_i=1.12R{i}_{*}^{-1}+0.08,\mathrm{with}{w}_m^2=w_{*}^2+4u_{*}^2+0.1\left(\mathsf{\Delta }{U}^2+\mathsf{\Delta }{V}^2\right),$$

(12)

which are denoted in this study as Model 3, Model 4, and Model 5, respectively.

Heng et al. [18] validated above models with their large-eddy simulation (LES) numerical results, and showed that the entrainment zone thickness predicted by Model 4 is consistent with the data of the barotropic and baroclinic atmosphere; however, the values of the proportional coefficients lead to deviations between the predicted results by the models and the numerical results.

As the potential temperature gradient in the free atmosphere layer plays a key role in the entrainment process, Hageli et al. [19], Fedorovich et al. [7], and Jianning et al. [20] proposed the following three models (denoted here as Model 6, Model 7, and Model 8, respectively) in terms of the buoyancy Richardson number Ri_N:

$$\mathsf{\Delta }z/z_i=C_1\cdot R{i}_N^{-1/2},$$

(13)

$$\mathsf{\Delta }z/z_i=C_2\cdot R{i}_N^{-2/5},$$

(14)

$$\mathsf{\Delta }z/z_i=C_3\cdot R{i}_N^{-1/3},$$

(15)

where C₁, C₂, and C₃ are constants, and $R{i}_N^{}$ is defined as

$$R{i}_N=\frac{z_i^2{N}^2}{w_{*}^2},$$

(16)

in which N is the Brunt-Väisälä frequency, and $N^2=\left(g/{\theta }_0\right){\mathsf{\Gamma }}_{\theta }$.

Extensive field observations and lab water tank experiments have shown that the fluctuations of the buoyancy flux near zero in the upper entrainment zone cause large deviations in the determination of the upper limit of $z_2$ [7,12,19,21], which makes the determination of $\mathsf{\Delta }\theta$ very difficult, and weakens the linearity of the fitted results of $\mathsf{\Delta }z/z_i$. Hence, Jianning et al. [20] took ${\mathsf{\Gamma }}_{\theta }$ instead of $\mathsf{\Delta }\theta$ as the atmospheric stability parameter, in order to improve the accuracy of the fitted relationship, and obtained another form of the bulk scaling model (denoted here as Model 9) based on the energy relationship of thermals in the entrainment zone:

$$\frac{\mathsf{\Delta }z}{z_i}=C_4{\left(\frac{{\theta }_{*}}{{\mathsf{\Gamma }}_{\theta }\cdot z_i}\right)}^{1/2},$$

(17)

where C₄ is a constant, and the potential temperature scale ${\theta }_{*}$ is defined as

$${\theta }_{*}={\overline{w^{\prime }{\theta }^{\prime }}|}_s/w_{*}.$$

(18)

The authors obtained $C_4=1.44$ with the data of Deardorff et al. [12] and Boers and Eloranta [16], but it was found that Model 9 failed to reflect the effect of the shear on the entrainment process.

In summary, shear stress in the CBL makes the bulk scaling models for the entrainment zone thickness based on shear-free CBL deviate significantly, because the models fail to take into account the compounding influence of characteristic parameters of the sheared CBL. In this paper, numerical cases using LES are carried out in an atmospheric domain, with different combinations of $\mathsf{\Delta }U$, ${\mathsf{\Gamma }}_{\theta }$, and ${\overline{w^{\prime }{\theta }^{\prime }}|}_i$, and the nine models identified above are examined. Finally, a modified bulk scaling model for the entrainment zone thickness for sheared CBL is developed with the LES results obtained.

2. Numerical Methodology

2.1. LES Description

LES decomposes the instantaneous motion of a turbulent eddy in the flow regime by some filtering function into two parts: a large-scale eddy and a small-scale eddy. LES computes the large-scale components directly and models the subgrid scales via a subgrid-scale (SGS) model that represents the effect of the unresolved scales on the resolved ones. As LES can capture the large-scale effects and coherent structures in unsteady processes, adequately reflecting the turbulent characteristics and overcoming the huge computational cost of direct numerical simulation, it has been widely used in the research of the atmospheric boundary layer [7,22,23].

For a horizontally homogeneous turbulent flow, the filtered governing equations are as follows:

$$\frac{\partial \overline{u_i}}{\partial x_i}=0,$$

(19)

$$\frac{\partial \overline{u_i}}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{u_i}\overline{u_j}\right)=\nu \frac{\partial }{\partial t}\left[\right(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i})-\frac{2}{3}\frac{\partial \overline{u_l}}{\partial x_l}{\delta }_{ij}]-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}-\frac{\partial {\tau }_{ij}}{\partial x_j},$$

(20)

$$\frac{\partial \overline{\theta }}{\partial t}+\frac{\partial \overline{u_i}\overline{\theta }}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{\lambda }{\rho C_p}\frac{\partial \overline{\theta }}{\partial x_i}\right)+\frac{1}{\rho C_p}(\frac{\partial \overline{p}}{\partial t}+\overline{u_j}\frac{\partial \overline{p}}{\partial x_i})-\frac{\partial {\tau }_{\theta i}}{\partial x_j},$$

(21)

where u_i and u_j are the velocity vectors with components along the directions x_i and x_j; $\rho$ is the density; $\nu$ is the kinematic viscosity; ${\tau }_{ij}$ is the subgrid-scale (SGS) stress; δ_ij is Kronecker delta; C_p and $\lambda$ are the specific heat capacity at constant pressure and thermal conductivity, respectively; and ${\tau }_{\theta i}$ is the SGS turbulence flux of the virtual potential temperature. This is modeled by using an SGS turbulent Prandtl number:

$${\tau }_{\theta i}=\frac{{\nu }_{SGS}}{{\mathrm{Pr}}_{SGS}}\frac{\partial \overline{\theta }}{\partial x_j},$$

(22)

in which ${\nu }_{SGS}$ is SGS viscosity and ${\mathrm{Pr}}_{SGS}$ is the SGS Prandtl number. The SGS turbulent stresses are calculated by

$${\tau }_{ij}=\frac{1}{3}{\tau }_{kk}{\delta }_{ij}-2{\nu }_{SGS}\overline{S_{ij}}.$$

(23)

As for incompressible flows in the CBL, the isotropic part of the SGS stress ${\tau }_{kk}$ can be neglected, as the turbulent velocity scale is too small. The rate-of-strain tensor for the resolved scale ${\overline{S}}_i{}_j$ is defined by

$$\overline{S_{ij}}=\frac{1}{2}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right).$$

(24)

In this study, the SGS model adopts eddy viscosity (ν_SGS) as the following [24]:

$${\nu }_{SGS}={[\min (kd,C_{SGS}\mathsf{\Delta })]}^2\sqrt{2\overline{S_i{}_j}\overline{S_i{}_j}},$$

(25)

where $k$ is the von Kármán constant; d is the distance to the closest wall; $C_{SGS}$ is the model coefficient, which is estimated based on the dynamic procedure; and $\mathsf{\Delta }$ represents the local grid scale based on the cell volume $\mathsf{\Delta }=V^{1/3}$.

Recent studies [25,26,27] have addressed that LES has some difficulties in getting entrainment layer dynamics right. These difficulties are mainly attributed to overestimation of eddy viscosity due to an excessive velocity gradient in the entrainment process. That is to say, $C_{SGS}$ cannot cover all cases with a universal constant. In an attempt to address these deficiencies, we therefore perform the dynamic Smagorinsky–Lilly SGS Model [28], in which $C_{SGS}$ is dynamically computed based on the information provided by the resolved scales of motion, which changes with the time of calculation and the spatial position of discrete points in the domain. In addition, other efforts have been made, such as reducing the grid spacing and solving the governing equations by using the central difference to eliminate the oscillation of the solution.

2.2. Large-Eddy Simulation Setup

Fitzjarrald [29] and Hadfield et al. [30] estimated the aspect ratio of the tank to be at least 4.5 to eliminate the sidewall influence on the central eddies in convective water tank experiments. Although the empirical value is based on the water tank, it is still valuable for the determination of computational domain size, and even an indication that even larger values are desirable was provided by some recent large-eddy simulations [31]. This paper carries out the LES runs in a horizontally homogeneous atmospheric domain of 6000 × 6000 × 2000 m³ (X × Y × Z) with an aspect ratio of 6 (the initial entrainment zone height is 1000 m).

The bulk of the computational domain is meshed by the grid size of 50 m × 50 m × 25 m, which is similar to that of Sullivan et al. [32], Fedorovich et al. [7], and Gentine et al. [33]. Because the vertical turbulent kinetic energy (TKE) distribution is not uniform, with buoyancy- and shear-generated TKE mainly concentrated in the underlying surface sublayer and the entrainment zone, the grid size near the surface varies steadily from 3 to 25 m and 15 m of the entrainment zone. Some tests are conducted to examine the numerical solutions’ sensitivity to grid resolution, and the result shows that the means and second-order statistics from the simulations on the chosen grid and a bigger grid (100 m × 100 m × 20 m) differ by no more than 5%, indicating that the current grid resolution provides accurate solutions.

The upper boundary is specified as a symmetrical rigid lid, where the vertical fluxes of physical parameters, such as mass, momentum, heat, and subgrid kinetic energy fluxes, are zero. As statistical characteristics of steady CBL turbulence are almost horizontally uniform, all runs are performed using periodic boundary conditions at the streamwise and side walls, as well as the no-slip boundary condition with constant buoyancy flux at the underlying surface.

2.3. Initial Conditions and Data Post-Processing

The simulated cases are uniformly initialized with an entrainment zone height of 1000 m, below which the potential temperature profile is at a constant 300 K, and above which it increases with a gradient ${\mathsf{\Gamma }}_{\theta }$. The horizontal velocities of the mixed zone and the free atmosphere are $U_1$ and $U_2$, respectively, and the velocity difference at the interface between the upper and lower zones is $\mathsf{\Delta }U$ ($\mathsf{\Delta }U=U_2-U_1$), in order to generate wind shear.

In this study, we consider different sheared CBLs by prescribing a series of combinations of $\mathsf{\Delta }U$, ${\mathsf{\Gamma }}_{\theta }$, and ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$, as listed in

Table 1. Details of the large-eddy simulation (LES) cases with different combinations of $\mathsf{\Delta }U$, ${\mathsf{\Gamma }}_{\theta }$ and ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$. The first letter in the case name represents the surface buoyancy flux (A: 0.05 K·m/s, B: 0.10 K·m/s, and C: 0.15 K·m/s,), the first number after the letter denotes the extent of the potential temperature gradient (1: 0.006 K/m, 2: 0.015 K/m, and 3: 0.003 K/m), and the second number after the letter represents the velocity difference across the entrainment zone (0: 0 m/s, 1: 1.5 m/s, 2: 3.0 m/s, and 3: 4.5 m/s).

${\overline{{\mathit{w}}^{\prime }{\mathit{\theta }}^{\prime }}|}_{\mathit{s}}\mathbf{K}·\mathbf{m}/\mathbf{s}$	${\mathsf{\Gamma }}_{\mathit{\theta }}=0.006\mathbf{K}/\mathbf{m}$				${\mathsf{\Gamma }}_{\mathit{\theta }}=0.015\mathbf{K}/\mathbf{m}$			${\mathsf{\Gamma }}_{\mathit{\theta }}=0.003\mathbf{K}/\mathbf{m}$
	△U = 0 m/s	△U = 1.5 m/s	△U = 3 m/s	△U = 4.5 m/s	△U = 1.5 m/s	△U = 3 m/s	△U = 4.5 m/s	△U = 1.5 m/s	△U = 3 m/s	△U = 4.5 m/s
0.05		A11	A12		A21	A22
0.1	B10	B11	B12	B13	B21	B22	B23	B31	B32	B33
0.15		C11	C12		C21	C22

, of which the main characteristic data are consistent with the actual atmospheric range, as presented in

Table 2. Typical field-observed values of the key parameters in CBL, obtained by Hibberd and Sawford [36] and their counterparts in the LES cases of the current study. Note: Pr, Re and Ra_f are the Prandtl number, Reynolds number, and the flux Rayleigh number, respectively. $\mathrm{Pr}=\nu /\alpha$, $\mathrm{Re}=U{z}_i/\nu$, and ${\mathrm{Ra}}_{\mathrm{f}}=B_s{z}_i^4/\nu D_s^2$, where $D_s$ is the coefficient of diffusion.

Convective System	U m/s	${\mathsf{\Gamma }}_{\mathit{\theta }}\mathbf{K}/\mathbf{m}$	${\overline{{\mathit{w}}^{\prime }{\mathit{\theta }}^{\prime }}|}_{\mathit{s}}\mathbf{K}\mathbf{m}/\mathbf{s}$	${\mathit{z}}_{\mathit{i}}\mathbf{m}$	${\mathit{w}}_{*}\mathbf{m}/\mathbf{s}$	Pr	Re × 108	Raf × 1022
Field observation	0–10	0.0003–0.03	0.03–0.3	500–2500	0.8–3.0	~0.71	0.3–5	1–7000
LES cases of this study	3–7.5	0.003–0.015	0.05–0.15	1000–1300	1.2–1.7	~0.71	0.21–5.1	25–80

; thus, the simulation results are comparable with the actual atmospheric data.

The momentum equation is discretized by the finite volume method, and the pressure–velocity coupling algorithm is based on the PISO (pressure implicit with splitting of operators) scheme [34]. The unsteady calculation time step of all simulations was 0.2 s, and the quasi-steady calculation time satisfies the experiential requirement of $t>6t_{*}$ [35], where the time scale $t_{*}$ is $t_{*}=z_i/w_{*}$, ensuring the TKE does not change significantly with time, and the flows are in a stable development stage at the end of the simulation.

Three-dimensional data volumes were archived from the final stage of the simulations to gather statistics and present flow visualization. The ensemble mean and turbulence statistics shown here are the result of spatial averaging—that is, at a fixed z location averaged over a x–y plane.

3. Evaluation of the Existing Bulk Models

3.1. Thickness of the Entrainment Zone

Since horizontal averages may be taken as substitutes for the ensemble averages in the horizontally homogeneous CBL, the horizontally averaged flow statistics are assumed to converge to corresponding ensemble means. From this point forward, the buoyancy flux of all grid nodes in the computational domain is averaged horizontally, in order to obtain the vertical buoyancy flux profiles of all simulated cases, as illustrated in Figure 2.

Figure 2 presents the evolution of the vertical profiles of buoyancy flux, which is made dimensionless by the surface buoyancy flux ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$, obtained from the LES results. The surface buoyancy flux is a major turbulence production mechanism in atmospheric boundary layer, and a stronger surface buoyancy flux enlarges the temperature difference between the thermals from the surface sublayer and the ambience, which determines the depth of the thermals penetrating into the free atmosphere and decisively impacts the entrainment development. From Figure 2a–c, it is seen that the increase of ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$ from 0.05 K·m/s to 0.15 K·m/s expands $\mathsf{\Delta }z$ by 10 to 70 m, indicating that the promotion of entrainment and its influence is particularly strong in cases with a weak stratified environment (the cases with smaller ${\mathsf{\Gamma }}_{\theta }$, i.e., A11, B11, C11), which is consistent with the analysis for shear-free CBLs [37].

As the behavior of thermals in the CBL is dominantly affected by negative buoyancy, the inhibition of entrainment by the potential temperature gradient of the free atmosphere is very obvious by comparing Figure 2d–f, and $\mathsf{\Delta }z$ significantly shrinks with the strengthening of ${\mathsf{\Gamma }}_{\theta }$. From the figure, it can be seen that the inhibition effect of ${\mathsf{\Gamma }}_{\theta }$ is relatively obvious to $z_2$, especially in the cases of weak surface buoyancy flux (Figure 2d), as ${\mathsf{\Gamma }}_{\theta }$ increases the resistance of the thermals to uprush. This is consistent with the simulated results of Boers and Eloranta [16] and the lidar observations of Hageli et al. [19].

Kim et al. [17] and Pino et al. [38] elucidated that the vertical shear stress of the horizontal wind destroyed the stability of the Kelvin–Helmholtz wave at the interface and enhanced entrainment. The thickness of the entrainment zone is directly related to the strength of the velocity difference under the same ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$ and ${\mathsf{\Gamma }}_{\theta }$ conditions, as shown in Figure 2g–i, especially in the case of weak inversion strength (smaller ${\mathsf{\Gamma }}_{\theta }$), where $\mathsf{\Delta }z$ increases more remarkably. For example, when $\mathsf{\Delta }U$ is increased from 1.5 m/s to 3 m/s, the $\mathsf{\Delta }z$ of the B31 case increased 65 m over that of the B32 case, while the B22 case increased only 28 m over the B21 case. This reveals that the turbulence intensity in the entrainment zone with weak inversion strength is smaller, and the shear-generated TKE consumed by the negative buoyancy in the entrainment zone increases accordingly, which makes the thickness increase more.

From the above observations, it is found that the characteristic factors of the boundary layer, i.e., $\mathsf{\Delta }U$, ${\mathsf{\Gamma }}_{\theta }$, and ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$ have significant influence on the development of the entrainment zone, with both $\mathsf{\Delta }U$ and ${\overline{w^{\prime }{\theta }^{\prime }}|}_s$ promoting entrainment and ${\mathsf{\Gamma }}_{\theta }$ inhibiting the entrainment process. Therefore, these parameters must be taken into account collectively, in order to obtain the bulk scaling models that are able to provide accurate predictions.

3.2. Evaluation of the Existing Bulk Models

The numerical results of the LES cases obtained in this study can be used to evaluate the accuracy of the bulk scaling models presented above, i.e., Models 1–9. The results of the evaluation are presented in Figure 3. The results show that the use of the shear-free velocity scale $w_{*}$ in $R{i}_{*}$ of Model 1, which does not includes the effect of shear in the entrainment zone, making the predictions by the model deviate substantially from the numerical results. Model 2 uses the entrainment velocity $w_e$ as the velocity scale, which implicitly reflects the influence of the CBL characteristic factors on entrainment. However, the predicted results are not in agreement with the numerical results.

Models 3–5 use the fraction velocity scale of the underlying surface and the velocity difference of the entrainment zone in the velocity scale, involving the shear effects occurring in both discrete sublayers of the sheared CBL. Deardoff and Willis [39] and Moeng and Sulliwan [40] concluded that shear-generated TKE is quickly dissipated by molecular viscous at a short depth before it can be utilized to drive the entrainment process through the entire mixed zone. Conzemius and Fedorovich [3,4] argued that the surface friction affects the vertical evolution of momentum and it has almost no impact on the entrainment process, through the simulation of typical daytime atmospheric CBLs. The overestimation of the contribution of the surface friction velocity in the velocity scale of Model 3 results scatters the numerical data more, and the correlations between $\mathsf{\Delta }z/z_i$ and $R{i}_{*}^{-1}$ cannot be clearly seen from the figure. The predictions from Models 4 and 5 agree reasonably well with the numerical results, although they systematically underestimate the results.

There is a large statistical deviation in the determination of the upper limit of the entrainment zone, which reduces the reliability of the derived models. However, this is supposed to be improved by replacing $\mathsf{\Delta }\theta$ with ${\mathsf{\Gamma }}_{\theta }$ in the bulk scaling model. Models 6–8 proposed the scaling relationship between the dimensionless thickness and the buoyancy Richardson number $R{i}_N$ with such a replacement. From the figure, it is clear that $\mathsf{\Delta }z/z_i$ and $R{i}_N$ are in good agreement with the −1/2 to −1/3 power-law correlations, which is in agreement with the conclusions from relevant studies, as mentioned above. However, because $w_{*}$ is still used as the velocity scale in these models, it excludes the effect of the shear; whether the accuracy of the models can be improved requires further study.

The temperature scale ${\theta }_{*}$, which implicitly reflects the influence of the surface buoyancy flux, is taken as the entrainment index in Model 9. The LES results presented in Figure 3 show that Model 9 improves the predictions. It is noteworthy that Model 9 cannot describe a clear entrainment process in sheared CBL, due to the lack of consideration of the shear effect on entrainment.

In conclusion, compared with the shear-free $R{i}_{*}^{-1}$ model (Model 1), the velocity scales $w_e$ and $w_m$ in Models 2–5 include the influence of the shear effect on entrainment, which are comparatively more in line with the numerical results (except for Model 3, due to the extensive estimation of the surface shear effect). In Models 6–8, $R{i}_N$ uses ${\mathsf{\Gamma }}_{\theta }$ instead of $\mathsf{\Delta }\theta$ as the inversion strength parameter in the prediction of the models. Although the prediction is much improved, Model 9 does not take into account of the effect of $\mathsf{\Delta }U$ on the entrainment process, resulting in some limitations of its application.

4. Improved Bulk Scaling Models in a Sheared Convective Boundary Layer

4.1. Development of the Improved Bulk Scaling Model

The shear-generated TKE in the entrainment zone is partly dissipated by the molecular viscosity in the turbulent flow, and the remainder is balanced by the negative buoyancy. The vertical motion of thermals in the mixed zone can effectively relax the turbulence of the mean wind in the middle of CBL, so that the shear-generated TKE at the surface can be dissipated within a short distance, thus reducing the upward transport of TKE. The shear at the surface has little effect on entrainment compared with the shear in the entrainment zone, as illustrated in water tank experiments and numerical simulations [7,16,39]. Therefore, the varied coefficient of $u_{*}$ in Equation (8) (i.e., a) is not sensitive enough to change the value of the velocity scale remarkably; therefore, it is feasible to ignore $u_{*}$ in the velocity scale without a significant impact on the model.

Based on the ZOM framework, Conzemius and Fedorovich [3,4] defined C_p as the ratio of the integral of the buoyancy flux difference with shear-free cases to the integral shear production between $z_s$ and $z_2$, and found that C_p was in the range of 0.2~0.5, with an average value of 0.35 in the sheared CBL atmospheric cases:

$$C_p=\frac{\raisebox{1ex}{$g$}\!\left/ \!\raisebox{-1ex}{${\theta }_0$}\right.{\displaystyle {\int }_0^{z_2}\overline{w^{\prime }{\theta }^{\prime }}dz-0.4}\overline{{w^{\prime }{\theta }^{\prime }|}_s}}{{\displaystyle {\int }_{z_s}^{z_2}\left({\tau }_x\frac{\partial u}{\partial z}+{\tau }_y\frac{\partial v}{\partial z}\right)dz}},$$

(26)

where the shear stress components are ${\tau }_x=-\overline{w^{\prime }{u}^{\prime }}$ and ${\tau }_y=-\overline{w^{\prime }{v}^{\prime }}$. We then suggest the following modified $w_m$:

$$w_m^2=w_{*}^2+C_p\left(\mathsf{\Delta }{U}^2+\mathsf{\Delta }{V}^2\right),$$

(27)

where C_p is the fraction of the shear-generated TKE available for entrainment. At the end of the calculation, the horizontal velocity components of the mixed layer are basically constant, due to the vertical turbulent mixing of the rising motion of the thermals, and the horizontal velocity of the free atmosphere is basically unaffected by the entrainment of the lower layer.

According to the thermal parcel theory, the negative buoyancy $\left(\rho -{\rho }_0\right)g$ in the entrainment zone ($z_i$~$z_2$) decelerates the upward movement of the thermals, from $w_m$ at the top of the mixed zone to zero when it reaches the top of the entrainment zone. The vertical energy conservation equation of thermals’ overshooting in the zone to resist negative buoyancy can be written as

$$0-\frac{1}{2}{\rho }_0{w}_m^2={\displaystyle {\int }_{z_i}^{z_2}\left(\rho -{\rho }_0\right)}gdz,$$

(28)

where ${\rho }_0$ is ambient atmospheric density. The left-hand side of the equation is the change of the kinetic energy of thermals, and the right-hand is the integrated buoyancy-generated work budget in the entrainment zone.

As the CBL is finite in height relative to the earth’s atmosphere, the effect of pressure on its density is negligible, but gravity is sufficiently strong. Therefore, the local density mainly depends on the potential temperature profile, and could be written as following with the Boussinesq approximation:

$$\rho ={\rho }_0\left[1-\frac{1}{{\theta }_0}\left(\theta -{\theta }_0\right)\right].$$

(29)

The entrainment zone thickness is mainly extruded by the potential temperature gradient of the free atmosphere in the CBL. Sun et al. [20] found an approximation between the potential temperature jump and its gradient in the entrainment zone, using the experimental data of Deardorff et al. [12] and the numerical data of Fedorovich et al. [7]:

$$\mathsf{\Delta }\theta /\mathsf{\Delta }z=\beta {\mathsf{\Gamma }}_{\theta }.$$

(30)

With Equations (29) and (30), the right-hand side of Equation (28) can be written as

$${\displaystyle {\int }_{z_i}^{z_2}\left(\rho -{\rho }_0\right)}gdz=-\frac{{\rho }_{0}g\beta {\mathsf{\Gamma }}_{\theta }{(\mathsf{\Delta }z)}^2}{2{\theta }_0};$$

(31)

Thus, Equation (28) can be expressed as

$$-\frac{1}{2}{\rho }_0{w}_m^2=-\frac{{\rho }_{0}g\beta {\mathsf{\Gamma }}_{\theta }{(\mathsf{\Delta }z)}^2}{2{\theta }_0}.$$

(32)

That is,

$$\frac{\mathsf{\Delta }z}{z_i}={\left(\frac{g/{\theta }_0\beta {\mathsf{\Gamma }}_{\theta }{z}_i^2}{w_m^2}\right)}^{-1/2}={\beta }^{-1/2}R{i}_{NC}{}^{-1/2},$$

(33)

where $R{i}_{NC}$ is the corrected Richardson number by the shear effect, defined as $R{i}_{NC}=N^2{z}_i^2/w_m^2$, and $w_m$ is in the modified form expressed by Equation (26).

4.2. Comparison of the Improved Model Against the Available Data

Otte and Wyngaard [41] simulated both shear-free and sheared atmospheric CBL in the range of 2~4.2 km, and obtained the vertical distributions of the potential temperature and velocity. The development of remote sensing technology has enabled us to simultaneously measure the parameters of the whole CBL. Boers and Eloranta [16] used lidar and a pilot balloon to measure the dry atmospheric data in the center of Illinois for six uninterrupted days in June 1979, and captured the vertical profiles of temperature and wind velocity in the CBL, as well as the surface heat flux, friction velocity, etc. Lothon et al. [42] used National Oceanic and Atmospheric Administration (NOAA) high-resolution doppler lidar (HRDL) to measure CBL parameters to determine the position of the entrainment zone.

With the aid of the data from these lidar measurements and numerical simulations, we fit the relationship between $\mathsf{\Delta }z/z_i$ and $R{i}_{NC}$ with the developed, improved model, and present the results in Figure 4. The following quantified scaling model is obtained from the presented data:

$$\frac{\mathsf{\Delta }z}{z_i}=1.26R{i}_{NC}^{-1/2}.$$

(34)

In this paper, C_p takes the statistical average of the turbulence simulation, but it is in fact related to the development process and entrainment factors of CBL. To date, the specific mechanism and rules of CBL are not very clear, and it is necessary to fully understand the entrainment process in future studies.

We took the average distance from the discrete points to the fitting lines as the index to describe the quality of models in Figure 3 and the new model (Equation (34), denoted as Model 10) in Figure 4. It can be seen from Figure 5 that the model developed in the current study provides an improved prediction.

5. Conclusions

In this study, LES was carried out to simulate the sheared CBL cases in a horizontally homogeneous atmospheric domain. The behaviors of the nine current bulk scaling models between the dimensionless entrainment zone thicknesses and different Richardson numbers ($R{\mathrm{i}}_{*}$, $R{\mathrm{i}}_E$, $R{\mathrm{i}}_N$) and ${\theta }_{*}/\left({\mathsf{\Gamma }}_{\theta }\cdot z_i\right)$ were evaluated against the LES results. The following conclusion can be drawn:

• The traditional $R{\mathrm{i}}_{*}$ model does not include the shear effect on entrainment, and thus does not perform well in the sheared CBL. Using $w_m$ as the velocity scale to include the shear effect on entrainment, the revised $R{\mathrm{i}}_{*}$ model provides some improvement in the predictions, but still with large deviations due to the large deviations in determining $\mathsf{\Delta }\theta$ in the entrainment zone.
• For the correlation with $\mathsf{\Delta }\theta$, ${\mathsf{\Gamma }}_{\theta }$ is adopted as the influence factor of inversion stratification on entrainment process in the $R{\mathrm{i}}_N$ model, which reduces the deviations of the predictions. However, this model is mainly aimed at a shear-free CBL, and also fails to take into account of the shear effect.
• In the $R{\mathrm{i}}_E$ model, the entrainment velocity $w_e$ at the interface is used as the velocity scale, which implicitly involves the comprehensive influence of buoyancy and shear on the entrainment process. However, each influence factor on the process cannot be quantified clearly, and the predictions are not very accurate.
• Similarly, due to the exclusion of the shear effect on entrainment, the “${\theta }_{*}/\left({\mathsf{\Gamma }}_{\theta }\cdot z_i\right)$” model also has some limitations in the predictions for sheared CBL.

A characteristic velocity scale $w_m$ based on the fraction of TKE used for entrainment, which includes the shear effect, was constructed in this study, and an improved bulk scaling model with this revised velocity scale has been derived from the modification of the previous methods and theoretical analysis of energy conservation in the entrainment process. The model uses $w_m$ and ${\mathsf{\Gamma }}_{\theta }$ as characteristic parameters, which fully reflects the influence of the boundary layer characteristic factors on the entrainment process. The correlation predicted by this model has been shown to provide significantly improved predictions.