## Abstract

Two zero-order bulk models (ZOMs) are developed for the velocity, buoyancy, and moisture of a cloud-free barotropic convective boundary layer (CBL) that grows into a linearly stratified atmosphere. The models differ in the entrainment closure assumption: in the first one, termed the “energetics-based model,” the negative and positive areas of the buoyancy flux are assumed to match between the model and the actual CBL; in the second one, termed the “geometric-based model,” the modeled CBL depth is assumed to match different definitions of the actual CBL depth. Parameterizations for these properties derived from direct numerical simulation (DNS) are employed as entrainment closure equations. These parameterizations, and hence the resulting models, are free from the potential singularity at finite wind strength that has been a major limitation in previous bulk models. The proposed ZOMs are verified using the DNS data. Model results show that the CBL depths obtained from the energetics-based model and previous ZOMs correspond to the height that marks the transition from the lower to the upper entrainment-zone sublayer; this reference height is few hundred meters above the height of the minimum buoyancy flux. It is also argued that ZOMs, despite their simplicity compared to higher-order models, can accurately represent CBL bulk properties when the relevant features of the actual entrainment zone are considered in the entrainment closures. The vertical structure of the actual entrainment zone, if required, can be constructed a posteriori using the available relationships between the predicted zero-order CBL depth and various definitions of the actual CBL depth.

## 1. Introduction

Bulk, or integral, models of a convective boundary layer (CBL) have been developed over the last few decades to parameterize bulk properties such as the CBL depth, the inversion strength, and the entrainment fluxes in atmospheric models whose grid spacings are much larger than the dynamically relevant scales of CBLs (Haltiner and Williams 1980; Suarez et al. 1983; Ayotte et al. 1996). Equally important, bulk models have broadly been used to investigate the sensitivity of the evolution of CBLs to changes in environmental conditions (Pelly and Belcher 2001; De Roode et al. 2014), and even to study process interaction (Naumann et al. 2017). Nonetheless, uncertainties still remain in some key aspects associated with the surface and entrainment closures. The work presented here focuses on the entrainment closure and is motivated by challenges identified in previous work, namely, the lack of agreement on the minimum complexity of the bulk model that is necessary to accurately represent sheared CBLs (Pino et al. 2006; Liu et al. 2016) and a singularity in the entrainment closure that can appear at finite wind strength (Driedonks 1982; Conzemius and Fedorovich 2004; Liu et al. 2018). In this work, we address these two issues.

Bulk models are classified based on their degree of complexity in the representation of the transition layer between the mixed layer and the free atmosphere. The simplest is the zero-order model (ZOM) (Lilly 1968) in which the transition layer is considered as an infinitesimally thin layer with discontinuous variations of velocity, buoyancy, and moisture. Alternatively, the first-order model (FOM) (Betts 1974) and higher-order models (Deardorff 1979) have been proposed, arguing that the transition layer between the mixed layer and the free atmosphere plays a key role in the dynamics of the CBL, and therefore, a better representation of the transition layer is required to accurately reproduce CBL bulk properties (Mahrt and Lenschow 1976; Sullivan et al. 1998; vanZanten et al. 1999; Kim et al. 2003). These models consider the transition layer as a layer of finite thickness with, respectively, linear and high-order polynomial variations of velocity, buoyancy, and moisture. Although the dependence of entrainment-zone properties on the entrainment-zone Richardson number, already shown by Mahrt and Lenschow (1976), might suggest that one needs at least the FOM to adequately represent the effect of entrainment in sheared CBLs, recent work has found no substantial differences between the overall ability of the ZOM and FOM to predict sheared CBL bulk properties (Pino et al. 2006; Conzemius and Fedorovich 2007). Nonetheless, Conzemius and Fedorovich (2007) found that the FOM largely mitigates—though not completely removes—the singularity of the ZOM at finite wind strength. Because of this advantage, they argued that the FOM is superior to the ZOM. Following this line of argumentation and also the necessity of predicting the finite thickness of the transition layer for some applications, most recent work made the effort to further develop a FOM (Sun and Xu 2009; Huang et al. 2011; Gentine et al. 2015; Liu et al. 2016). In this paper, we show that the infinitesimal transition-layer representation of the ZOM is sufficient to precisely reproduce bulk properties in the cloud-free sheared CBL, as long as the entrainment closure appropriately represents the local effects of wind shear on entrainment. If required, the actual finite thickness of the transition layer can be constructed a posteriori at the top of the predicted zero-order CBL depth using either the relationships between the zero-order CBL depth and various actual CBL depths provided in Haghshenas and Mellado (2019) or the transition-layer parameterizations proposed in previous work (Pino et al. 2006; Kim et al. 2006; Liu et al. 2016). (By the term “actual CBLs,” we explicitly mean atmospheric CBLs or three-dimensional simulations of them.)

The second aspect that we address here is the singularity in the entrainment closure that can appear at finite wind strength in previous bulk models. A parameterization for the entrainment-flux ratio, defined as the negative of the ratio of a buoyancy flux at the CBL top to the surface buoyancy flux, is commonly used as the entrainment closure in the integral equations. This parameterization, also referred to as the entrainment parameterization or the entrainment equation, is developed by either a local analysis of the turbulence kinetic energy (TKE) budget (Zeman and Tennekes 1977; Tennekes and Driedonks 1981; Driedonks 1982) or an integral analysis of the TKE budget (Boers et al. 1984; Batchvarova and Gryning 1994; Conzemius and Fedorovich 2006b). Previous entrainment parameterizations suffer from a potential singularity at finite wind strength, which is a major long-standing limitation of previous zero-order and first-order models (Driedonks 1982; Conzemius and Fedorovich 2004; Liu et al. 2018). This singularity arises when the entrainment parameterization is derived in the idealized framework of the bulk models and the CBL depth is used in the scaling of the shear production at the CBL top in the local TKE approach (see, e.g., Tennekes and Driedonks 1981), or is used in the scaling of the integral of the negative buoyancy flux in the integral TKE approach (see, e.g., Boers et al. 1984). This is physically inconsistent with the observation that, under strong-wind-shear conditions, the entrainment zone, defined as the region of negative buoyancy flux at the boundary layer top, is characterized by a local length scale that is different from the CBL depth (Zeman and Tennekes 1977; Kim et al. 2003; Pino and Vilà-Guerau De Arellano 2008; Haghshenas and Mellado 2019). Applying this local length scale in the integral analysis of the TKE budget, Haghshenas and Mellado (2019) derived nonsingular parameterizations for different CBL properties. In the present work, we exploit these parameterizations to develop nonsingular zero-order bulk models for the velocity, buoyancy, and moisture.

We structure the paper as follows. After describing the formulation in section 2, we summarize the derivation of the set of equations in the zero-order bulk model and briefly discuss the closures used in previous work in section 3. In section 4, we introduce two new entrainment closure equations and develop two ZOMs based on them. Evaluation of the proposed models is done in section 5 by comparing their predictions with data from direct numerical simulation (DNS). In section 6, we compare the predictions of the proposed ZOMs with those of previous work, and numerically and analytically investigate the potential singularity observed in the previous entrainment parameterization in section 7. One of the developed models is then used in section 8 to address the dependence of the sheared CBL on environmental conditions. We finally summarize these results and draw conclusions in section 9.

## 2. Formulation

We consider a cloud-free CBL forced by a constant and homogeneous buoyancy flux at the surface *B*_{0} growing into a linearly stratified dry free atmosphere with a Brunt–Väisälä frequency *N*_{0} (Fig. 1). The background profile of buoyancy is

where *z* is the vertical distance from the surface. Henceforth, the subscript “bg” denotes background and the symbol ≡ indicates a definition. The buoyancy is approximated as *b* ≃ *g*(*θ*_{υ} − *θ*_{υ,0})/*θ*_{υ,0}, where *θ*_{υ} is the virtual potential temperature and *θ*_{υ,0} is its constant reference value obtained by extrapolating the linear variation of *θ*_{υ} in the free atmosphere toward the surface. The CBL develops over an aerodynamically rough surface with constant surface roughness *z*_{0}. In addition, we consider a barotropic case, which implies that the wind strength in the free atmosphere *U*_{0} is constant with height, and we consider the limit of zero Coriolis parameter.

The background profile of specific humidity is

where *γ*_{q} ≥ 0 is the lapse rate, and *q*_{bg,0} is the background specific humidity at the surface obtained by extrapolating the linear variation of specific humidity in the free atmosphere toward the surface. In addition, we assume that the surface kinematic flux of specific humidity, *F*_{q,0} ≥ 0, is constant and homogeneous.

### a. Governing equations

The set of governing equations comprises the conservation equations for mass, momentum, energy, and moisture in the Boussinesq approximation. Under the assumptions of statistical homogeneity in the horizontal directions, no subsidence, and neglecting condensation and radiation, and in the limit of zero Coriolis parameter, the horizontally averaged equations for streamwise kinematic momentum *u*, buoyancy *b*, and specific humidity *q* read

We have chosen the streamwise coordinate *x* aligned with the wind in the free atmosphere, so that the mean wind in the spanwise direction is zero. The variables *τ*_{x}, *B*, and *F*_{q} are, respectively, the mean vertical fluxes of streamwise kinematic momentum, buoyancy, and specific humidity. Angle brackets denote averaging along horizontal planes. Formulating the system in terms of buoyancy and moisture instead of temperature and moisture facilitates the study of the sensitivity of moisture properties to changes in environmental conditions. The reason is that this formulation along with the linearization of the equation of state mathematically cause the moisture to become a passive scalar, meaning that changing moisture without changing buoyancy does not alter the CBL dynamics (Mellado et al. 2017). The energy variable (e.g., potential temperature or static energy) can be recovered from the buoyancy, the specific humidity, and the linearized equation of state.

### b. Dimensional analysis

In the limit of a high Reynolds number and once the initial conditions have been sufficiently forgotten, the dynamics of the sheared CBL is completely governed by the control parameters {*B*_{0}, *N*_{0}, *U*_{0}, *z*_{0}} and the independent variables {*z*, *t*}, where *t* represents the time. We focus on the quasi-steady (equilibrium) entrainment regime under which CBL properties evolve on time scales much larger than the large-eddy turnover time, and the profiles of various properties, when appropriately normalized, behave approximately self-similarly (Fedorovich et al. 2004). Hereafter, we will use the term “quasi-steady regime” for simplicity.

The system in the quasi-steady regime, hence, depends on two nondimensional parameters: a reference Froude number,

and a normalized surface roughness, *z*_{0}/*L*_{0}. Here *L*_{0} is the reference Ozmidov length defined as

which provides a relevant measure for the thickness of the upper region of the entrainment zone in the shear-free and sheared CBL (Garcia and Mellado 2014; Haghshenas and Mellado 2019).

Statistical properties of moisture in the quasi-steady regime depend on three parameters {*F*_{q,0}, *γ*_{q}, *q*_{bg,0}} in addition to the aforementioned nondimensional parameters. Mellado et al. (2017) have shown that moisture statistics can be conveniently analyzed by the nondimensional parameter

The parameter *F*_{q,1} is a reference scale for the entrainment flux of the specific humidity and is defined as

which can be interpreted as the product of a moisture variation *γ*_{q}*L*_{0} and a velocity scale *N*_{0}*L*_{0} in the upper region of the entrainment zone. The quantity *φ* is a flux-ratio parameter that varies, by definition, between 0, which corresponds to the pure-drying regime, and 2, which corresponds to the pure-moistening regime.

We express the dependence of statistical properties on time in terms of the nondimensional variable *z*_{enc}/*L*_{0}. The variable *z*_{enc} is the encroachment length scale (Lilly 1968; Carson and Smith 1975) defined as

where *z*_{∞} is the height sufficiently far above the CBL top so that the integral is approximately independent of *z*_{∞}. The integral analysis of the buoyancy equation in the limit of a high Reynolds number yields

where *t*_{0} is a constant of integration, which quantifies the dependence on the initial buoyancy profile.

The logic behind using *z*_{enc}/*L*_{0} instead of *tN*_{0} to represent the state of the CBL development is that it facilitates the comparison between atmospheric measurements and results from numerical simulations conducted with different initial conditions. Notice that the encroachment length scale provides a measure for the depth of the mixed layer in shear-free and sheared CBLs growing into a linearly stratified atmosphere (van Heerwaarden and Mellado 2016; Mellado et al. 2016; Haghshenas and Mellado 2019), and it can be easily calculated from the mean buoyancy profile, obtained from atmospheric measurements or numerical simulations, according to Eq. (8).

For typical midday conditions of the sheared CBL over land, one finds *N*_{0} ≃ 0.006–0.018 s^{−1}, *B*_{0} ≃ 0.001–0.01 m^{2} s^{−3}, *U*_{0} ≃ 0–20 m s^{−1}, *z*_{0} ≃ 0.01–0.1 m, *γ*_{q} ≃ 0–0.002 g kg^{−1} m^{−1}, *F*_{q,0} ≃ 0.03–0.1 g kg^{−1} m s^{−1}, and *z*_{enc} ≃ 500–2000 m (Conzemius and Fedorovich 2006a; Garcia and Mellado 2014; Mellado et al. 2017), which yields the parameter space Fr_{0} ≃ 0–85, *z*_{0}/*L*_{0} ≃ (0.05–5) × 10^{−3} and *z*_{enc}/*L*_{0} ≃ 5–50. The parameter *φ* can change between its theoretical limits, 0 and 2.

## 3. Zero-order bulk model

In this section, we summarize the derivation of ZOM equations for a barotropic CBL without the Coriolis force, and we discuss the basic form of previous surface and entrainment closures and their corresponding limitation and uncertainty as needed for the discussion in the following sections. Further details of the derivation of the equation set and closures can be found, for example, in Fedorovich (1995) and Conzemius and Fedorovich (2006b).

### a. Derivation of zero-order model equations

The set of zero-order model equations is derived by approximating the actual properties with the ZOM properties, by vertically integrating Eq. (3) from the surface *z* = 0 up to a height that is slightly above the CBL depth *z* = *h*^{(0)} + *ε* and taking the limit *ε* → 0 after the integration, and by evoking basic assumptions of the zero-order representation of the CBL vertical structure. This analysis yields

The superscript “(0)” indicates the zero-order bulk model, and we use a prefix “zero-order” to distinguish quantities in the model from those in the actual CBL. The zero-order fluxes of buoyancy, kinematic momentum, and specific humidity at the CBL top are related to the zero-order increments of these properties at the CBL top and the growth rate of the CBL depth as

These equations are derived by vertically integrating Eqs. (3) over the height from *z* = *h*^{(0)} − *ε* up to *z* = *h*^{(0)} + *ε* and taking the limit *ε* → 0 after the integration. In addition to the three fluxes introduced above, the unknown variables are *h*^{(0)}, $u*\u2061(0)$, Δ*u*^{(0)} (or alternatively the mixed-layer velocity $uml\u2061(0)=U0\u2212\Delta u\u2061(0)$), Δ*b*^{(0)} (or alternatively the mixed-layer buoyancy $bml\u2061(0)=N02h\u2061(0)\u2212\Delta b\u2061(0)$), and Δ*q*^{(0)} (or alternatively the mixed-layer specific humidity $qml\u2061(0)=qbg,0\u2212\gamma q\u2009h\u2061(0)+\Delta q\u2061(0)$). Therefore, two more equations are required to close the system of equations.

### b. Surface closure equation

Previous work has often considered the surface-drag relationship

as the surface closure equation (see, e.g., Boers et al. 1984; Garratt 1992). The parameter *C*_{D} is the surface-drag coefficient, which is derived from the Monin–Obukhov similarity theory as a function of the surface roughness, Obukhov length, and surface-layer depth [see Garratt et al. (1982) for a review]. This functional relationship for the CBL that grows over an aerodynamically rough or smooth surface is explained in detail in appendix A for completeness. A constant value for the surface-drag coefficient is, however, usually taken for simplification (Flamant et al. 1999; Kim et al. 2006; Conzemius and Fedorovich 2007), in which case the dependence of CBL properties on the normalized surface roughness *z*_{0}/*L*_{0} (cf. section 2b) translates into a dependence on the surface-drag coefficient *C*_{D}.

### c. Entrainment closure equation

Previous work has often developed a parameterization for the entrainment-flux ratio as the entrainment closure equation. The basic form of this parameterization in previous work reads

which is derived by either a local analysis of the TKE budget (Zeman and Tennekes 1977; Tennekes and Driedonks 1981; Driedonks 1982) or an integral analysis of the TKE budget (Boers et al. 1984; Batchvarova and Gryning 1994; Conzemius and Fedorovich 2006b), assuming that, respectively, the local or bulk energetics in the model and in the actual CBL match. The variable

is the convective velocity scale (Deardorff 1970), and the variables

and

are the bulk Richardson numbers associated, respectively, with the accumulation term and the entrainment-zone wind shear.

The parameters *C*_{1}, *C*_{T}, *C*_{P}, and *A* are empirical constants. The constant *C*_{1} corresponds to the zero-order entrainment-flux ratio in the shear-free limit and its most commonly used value is 0.2 (see Table 1). The constant *C*_{T} corresponds to the contribution of the accumulation term, which is negligibly small with respect to the other terms in the TKE budget equation once the quasi-steady regime is reached. The main differences among entrainment parameterizations in previous work are the constants associated to wind shear effects on entrainment, namely, *C*_{P} and *A*, which correspond to the contributions from the entrainment-zone and surface wind shear, respectively, to the entrainment flux.

We note that Sun and Xu (2009) and Liu et al. (2016) have derived the entrainment parameterization in the FOM framework, but we obtain their corresponding parameterization in the ZOM framework by setting to zero the variable that represents the thickness of the transition layer between the mixed layer and the free atmosphere. This approach has also been applied in Pino et al. (2006) and Conzemius and Fedorovich (2007). As explained in Conzemius and Fedorovich (2007), the same assumption regarding the TKE source and sink terms (e.g., the dissipation terms are scaled in the same manner in both model frameworks) has been made in previous work to derive the entrainment parameterization in ZOM and FOM. This indicates that the empirical constants *C*_{1}, *C*_{T}, *C*_{P}, and *A* have the same physical meaning in both model frameworks.

Aside from the large uncertainty in the empirical constants *C*_{P} and *A*, the main limitation of the previous entrainment parameterization is the potential singularity at finite wind strength. The contribution of the entrainment-zone shear to the entrainment flux represented by a negative sign term in the denominator of Eq. (13) could cause the denominator to become zero and the entrainment-flux ratio to become unbounded. Such a singularity occurs not only under strong-shear conditions characterized by large Froude numbers, but also under moderate-shear conditions with initial conditions that are far from the quasi-steady regime (Driedonks 1982; Conzemius and Fedorovich 2004, 2007; Liu et al. 2018). As explained in the introduction, this singularity arises when the CBL depth, instead of a local scale in the entrainment zone, is used to estimate the TKE budget equation in the entrainment zone. We will further discuss this singularity in section 7.

## 4. Nonsingular entrainment closure equations

The work presented here focuses on the entrainment closure equation and, following previous work, uses the surface-drag relation as the surface closure equation [cf. Eq. (12)]. As entrainment closure equations, we introduce two new nonsingular equations by making two different closure assumptions and by employing the nonsingular parameterizations for different CBL properties derived in Haghshenas and Mellado (2019). In contrast to previous work, these authors have considered the actual CBL structure instead of the bulk-model structure, and have used the local length scale to characterize the entrainment zone, which led to nonsingular parameterizations.

### a. Energetics-based closure

As the first option for the entrainment closure, we assume that the negative and positive areas of the buoyancy flux in the model equal the ones in the actual CBL. This assumption is similar to that used in previous work, where the bulk or local energetics between the model and the actual CBL were assumed to be equal (cf. section 3c).

Mathematically, the zero-order entrainment-flux ratio can be written in terms of the energetics as (Conzemius and Fedorovich 2006a)

evoking basic assumptions of the zero-order representation of the CBL vertical structure (cf. Fig. 1). Here $AN\u2061(0)$ and $AP\u2061(0)$ are, respectively, the negative and positive areas of the buoyancy flux in the zero-order model framework.

A parameterization for the ratio of the negative and positive buoyancy flux in the actual CBL is obtained in appendix B using the results of Haghshenas and Mellado (2019). This parameterization, Eq. (B4), along with Eq. (17) and the closure assumption that the negative and positive areas of the buoyancy flux in the actual CBL equal the ones in the model yields

We will refer to this entrainment closure as energetics-based closure and to the model that uses this closure as energetics-based model.

The energetics-based closure indicates that the shear-free entrainment flux of buoyancy is solely characterized by the surface buoyancy flux as the only source of turbulence in this case, and that the entrainment flux of buoyancy in the sheared CBL increases due to extra turbulence generated by the entrainment-zone shear. Explicitly representing the extra shear-generated turbulence, which is a relevant feature of the dynamics of actual sheared CBLs (see review in Fedorovich and Conzemius 2008), raises the expectation that the energetics-based ZOM, despite the simplification in the CBL structure, should be able to faithfully represent the CBL bulk properties.

We note that we obtain $\u2212Bh\u2061(0)/B0\u22430.21$ in the shear-free limit, which agrees with Liu et al. (2016) and is only slightly larger than the value 0.2 that is commonly used in previous work (cf. Table 1). This result indicates that we start to determine this empirical constant with an error of 5% or less.

### b. Geometric-based closure

The CBL depth predicted by the energetics-based model and by models in previous work cannot be a priori associated to any actual CBL depth, such as the height of the minimum buoyancy flux or the height of the maximum buoyancy gradient. These heights differ by approximately 100 m for typical midday atmospheric conditions over land in the shear-free limit, and might increase to 200 m under strong-shear conditions in the barotropic CBL. This uncertainty about the CBL depth might be relevant for the parameterization of other processes in the CBL, for instance, cloud formation. To reduce this uncertainty, we propose a new model in which the zero-order CBL depth is assumed to match different definitions of the actual CBL depth, and then we compare the results of this model with those of the energetics-based model and models in previous work.

Several definitions of the actual CBL depth might be associated with the zero-order CBL depth. Two common definitions are the height of the minimum buoyancy flux (Fedorovich et al. 2004; Pino et al. 2003, 2006) and the height of the maximum buoyancy gradient (Sullivan et al. 1998). These heights fall within the lower and upper entrainment-zone sublayers, respectively (Garcia and Mellado 2014; Haghshenas and Mellado 2019). A third option is the height that marks the transition between these two sublayers. Using Eqs. (4.12) and (5.11) in Haghshenas and Mellado (2019) as parameterizations for these heights, we obtain

where *α* is a free parameter and is discussed below. We will refer to this entrainment closure as geometric-based closure and to the model that uses this closure as geometric-based model. Comparing the CBL depth for the shear-free limit obtained from Eq. (19) with those reported in Garcia and Mellado (2014) and Mellado et al. (2017), where the Reynolds number in simulations was approximately 4 times larger than in Haghshenas and Mellado (2019), shows Reynolds-number dependence of less than 2% in Eq. (19).

Taking *α* as a free parameter enables us to consider two different definitions of the CBL depth as the zero-order CBL depth: *α* ≃ 0.8 corresponds to the height of the minimum buoyancy flux, and *α* ≃ 1.0 corresponds to the height that marks the transition from the lower to the upper entrainment-zone sublayer. The second term in the right-hand side of Eq. (19) represents the local length scale that characterizes the lower entrainment-zone sublayer. Given that the height of the maximum buoyancy gradient is located inside the upper entrainment-zone sublayer and that the characteristic length scale of this sublayer differs from the one of the lower entrainment-zone sublayer, the parameterization for the height of the maximum buoyancy gradient has an additional contribution to Eq. (19). This parameterization is discussed in appendix C but not in the main text, for conciseness. We note that, however, the height of the maximum buoyancy gradient as well as the other actual upper depth of the CBL can be constructed a posteriori at the CBL top obtained from both energetics- and geometric-based models using the relationships between the predicted zero-order CBL depth and these actual CBL depths as explained in appendix C.

The geometric-based closure indicates that the zero-order CBL can be interpreted as a two-layer entity, namely, a buoyancy-driven layer that represents the actual mixed layer, and a buoyancy- and shear-driven layer that represents part of the actual entrainment zone. This is important because, as mentioned before, the finite thickness of the entrainment zone is a relevant feature of actual sheared CBLs (Mahrt and Lenschow 1976; vanZanten et al. 1999; Kim et al. 2003), and explicitly representing this feature as in Eq. (19) raises the expectation that this ZOM should faithfully represent the CBL bulk properties.

### c. Closed set of zero-order model equations

The closed set of ZOM equations in nondimensional form for the buoyancy and velocity are derived from Eqs. (10a) and (10b) as follows:

plus either the energetics-based closure, Eq. (18), which can be written as

or the geometric-based closure, Eq. (19). Three steps have been taken to derive the equations above. First, we have substituted Eqs. (12) and (11a) in Eqs. (20a) and (21), respectively. Second, we have integrated Eq. (10b) with respect to time. Third, we have changed the variable from *t* to *z*_{enc} in all three equations above, where $d/dt=\u2061(N0L02zenc\u22121)\u2061(d/dzenc)$ [cf. Eq. (9)]. Recall that Fr_{0} and *C*_{D} are the nondimensional control parameters.

The ZOM equation for the moisture reads

which is obtained by integrating Eq. (10c). Here *φ* is the flux-ratio parameter that characterizes the moisture [Eq. (6)], and *q*_{ref} is the reference moisture scale

which is defined as a linear combination of *F*_{q,0} and *F*_{q,1} normalized by a velocity scale *N*_{0}*L*_{0} (Mellado et al. 2017). Normalization of Δ*q*^{(0)} with *q*_{ref} instead of *F*_{q,0}/(*N*_{0}*L*_{0}) or *F*_{q,1}/(*N*_{0}*L*_{0}) allows us to study the whole theoretical range of the flux-ratio parameter, since *q*_{ref} remains nonzero for both limits of *φ* = 0, which corresponds to *F*_{q,0} = 0 (pure-drying regime), and *φ* = 2, which corresponds to *F*_{q,1} = 0 (pure-moistening regime).

A key property to characterize the moisture is the critical zero-order flux-ratio parameter $\phi cr\u2061(0)$ that marks the boundary between drying and moistening regimes. A functional relationship for $\phi cr\u2061(0)$ can be readily determined from the condition $dqml\u2061(0)/dt=0$, which, by definition [obtained from Eq. (10c)], corresponds to (Mahrt 1991)

Substituting $Fq,h\u2061(0)$ from Eq. (11c) in the equation above, and using Eq. (11a) to rewrite *dh*^{(0)}/*dt* in the resulting equation in terms of Δ*b*^{(0)} and $Bh\u2061(0)$ yields

as the condition that marks the boundary between drying and moistening regimes. If the left-hand side (lhs) of Eq. (25) is larger than the right-hand side (rhs), the entrainment flux of drying air is dominant and the CBL is in the drying regime. When the lhs is smaller than the rhs, the surface flux of moisture dominates and the CBL is in the moistening regime. When the lhs equals the rhs, the mean moisture $qml\u2061(0)$ remains constant in time and the water vapor introduced at the surface is used to moisten the entrained dry air toward the mixed-layer value.

The critical zero-order flux-ratio parameter $\phi cr\u2061(0)$ can be written in terms of the CBL depth and the entrainment rate as

by substituting $Bh\u2061(0)$, Δ*b*^{(0)}, and Δ*q*^{(0)} from Eqs. (11a), (20b), and (22) in Eq. (25) and by solving the resulting equation for *φ*. The idea behind providing this functional relationship is that it allows us to determine whether the CBL is in the drying or in the moistening regime. The CBL is in the moistening regime when the flux-ratio parameter, calculated from Eq. (6), is larger than the critical value, determined from Eq. (26), and the CBL is in the drying regime when the flux-ratio parameter is smaller than the critical value.

The critical flux-ratio parameter is constant in time for the shear-free limit, when *h*^{(0)} ∝ *z*_{enc}, but varies in time for the sheared CBL, according to Eq. (26). The reason is that the entrainment enhancement due to the wind shear—causing the critical flux-ratio parameter to increase with respect to the shear-free limit—diminishes as the CBL grows. Therefore, the critical flux-ratio parameter decreases with time and asymptotes toward the corresponding shear-free value. We further discuss this behavior in section 8.

## 5. Validation of ZOM predictions with DNS data

We validate the ZOMs proposed in this work with the DNS data presented in Haghshenas and Mellado (2019) for a reference Reynolds number: $Re0\u2261B0/\u2061(\nu N02)=25$, where *ν* is the kinematic viscosity. Although this DNS data has been used in the derivation of the entrainment closures and therefore a good agreement between ZOM predictions and DNS data is expected, this validation already provides insight into how the CBL depth predicted from the energetics-based model relates to various definitions of the CBL depth obtained from the DNS data and used to construct the geometric-based models. We consider the shear-free case Fr_{0} = 0 and the shear case Fr_{0} = 25, which corresponds to mean wind velocities in the free atmosphere of 15 m s^{−1} in typical midday atmospheric conditions. We initialize the ZOMs using the DNS data at *z*_{enc}/*L*_{0} = 15, which corresponds to conditions inside the quasi-steady regime. Given that the simulated CBLs develop over an aerodynamically smooth surface, we consider the corresponding functional relationship for the surface-drag coefficient [Eq. (A1) with $z0=0.13\nu /u*$] in the surface closure equation, Eq. (12).

Figure 2 shows that the ZOM predictions agree with the DNS data, as anticipated. In particular, Fig. 2a illustrates that the CBL depths obtained from the geometric-based model with *α* = 0.8 and *α* = 1, respectively, collapse on top of the height of the minimum buoyancy flux *z*_{i,f} and the height that marks the transition from the lower to the upper entrainment-zone sublayer *z*_{i,s}. In addition, Fig. 2b shows that the zero-order entrainment-flux ratio obtained from the energetics-based model falls on top of the square root of the ratio between the negative and positive areas of the buoyancy flux acquired from the DNS data. All curves in Fig. 2c approximately fall on top of each other, because the corresponding surface closure is the same among the models and the deviation in the CBL depth is small (approximately 10% or less).

Interestingly, Fig. 2a also shows that the CBL depth obtained from the energetics-based model better corresponds to the height that marks the transition from the lower to the upper entrainment-zone sublayer, rather than the height of the minimum buoyancy flux. Concomitantly, the energetics obtained from the geometric-based model with *α* = 1.0 match with those of the energetics-based model and the DNS data (see Fig. 2b).

These relationships between the CBL depths represented in different ZOMs become clearer in the shear-free limit, where the set of ZOM equations has an analytical solution of the following form (see, e.g., Fedorovich et al. 2004; Liu et al. 2018):

the model coefficients satisfying the following relationships:

One of the three model coefficients {*C*_{1}, *C*_{2}, *C*_{3}} remains free and has to be prescribed to close the system. In the energetics-based model, we prescribe *C*_{1} = 0.21 according to Eq. (18) and obtain *C*_{2} ≃ 1.19 and *C*_{3} ≃ 0.18. In the geometric-based model with *α* = 1.0, we prescribe *C*_{2} ≃ 1.19 according to Eq. (19) and obtain *C*_{1} ≃ 0.21 and *C*_{3} ≃ 0.18. Hence, the predictions using the energetics-based model coincides with the predictions using the geometric-based model based on the CBL height that separates the lower and upper entrainment-zone sublayers.

In contrast, the geometric-based model with *α* = 0.8, where the zero-order CBL depth matches the height of the minimum buoyancy flux, prescribes *C*_{2} ≃ 1.14 according to Eq. (19) and yields *C*_{1} ≃ 0.15 and *C*_{3} ≃ 0.13. These coefficients are different than those of the energetics-based model. This difference helps explain the controversy in *C*_{1} in some previous work. Fedorovich et al. (2004) and Mellado et al. (2017) estimated *C*_{1} ≃ 0.17 and *C*_{1} ≃ 0.16 ± 0.01 assuming that the zero-order CBL depth matches the height of the minimum buoyancy flux in large-eddy simulations and direct numerical simulation, respectively. These values are smaller than 0.2 simply as a result of the aforementioned assumption, and not because of statistical uncertainty.

Similarly, the ZOM moisture properties in the shear-free limit can be obtained explicitly as

where the model coefficients satisfy the following relationships:

We obtain *C*_{4} ≃ 1.19 − 0.17*φ* and *C*_{5} ≃ 1.4 − 0.2*φ* for the energetics-based model and the geometric-based model with *α* = 1, and *C*_{4} ≃ 1.14 − 0.13*φ* and *C*_{5} ≃ 1.3 − 0.15*φ* for the geometric-based model with *α* = 0.8. The transition from drying to moistening regimes occurs for $\phi cr\u2061(0)\u22431.17$ for the energetics-based and the geometric-based models with *α* = 1, as obtained from Eqs. (26) and (27). One also obtains $\phi cr\u2061(0)\u22431.13$ for the geometric-based model with *α* = 0.8. From the definition of the flux-ratio parameter by Eq. (6), the condition *φ* = *φ*_{cr} ≃ 1.13–1.17 implies that *F*_{q,0} ≃ (1.3–1.4)*F*_{q,1} at the boundary between drying and moistening regimes in a shear-free CBL. Hence, the reference entrainment flux *F*_{q,1} provides a good estimate for the actual entrainment flux of moisture at the boundary between drying and moistening regimes, and its definition by Eq. (7) allows us to estimate a priori if a shear-free CBL is in the drying or moistening regime from estimates of the environmental conditions, before doing a simulation or a bulk model analysis (Mellado et al. 2017).

## 6. Comparison with previous work

In this section, we compare the proposed ZOMs with previous work to better understand the importance of the assumptions and idealizations used in the different approaches. From the previous work presented in Table 1, we consider only those that take into account the entrainment-zone shear (Conzemius and Fedorovich 2006b; Pino et al. 2006; Sun and Xu 2009; Liu et al. 2016), since a consensus has been reached in the literature that the shear-generated turbulence in the entrainment zone, and not the one generated in the surface layer, accounts for the entrainment enhancement in sheared CBLs (Fedorovich and Conzemius 2008; Pino and Vilà-Guerau De Arellano 2008).

As reference case for the comparison among the different models, we choose the strongest-shear case investigated in Pino et al. (2006) and initialize the ZOMs using the values listed in Table 1 of that work (see Table 2). The reference Froude number and the reference Ozmidov length for the selected case are, respectively, Fr_{0} ≃ 41 and *L*_{0} ≃ 34 m, according to Eqs. (4) and (5). We determine the encroachment length using Eq. (8) and integrating the vertical profile of the mean buoyancy presented in Fig. 2 of Pino et al. (2006), obtaining *z*_{enc} ≃ 510 m. We also determine *t*_{0} ≃ −35 s from Eq. (9), which is consistent with the fact that the corresponding numerical simulation has been started from a very shallow initial CBL. The surface is aerodynamically rough and we use a constant surface-drag coefficient in the surface closure equation.

Figure 3 illustrates that, qualitatively, all considered ZOMs appropriately represent relevant features of sheared CBLs that have been documented in previous observational and numerical studies (Mahrt and Lenschow 1976; Pino and Vilà-Guerau De Arellano 2008; Haghshenas and Mellado 2019). First, the CBL depth and the entrainment-flux ratio increase with respect to the shear-free limit. Second, as time goes by and the CBL depth grows, wind shear effects diminish, and the CBL properties asymptote toward the shear-free values. Last but not least, wind shear effects on buoyancy-related properties remain constrained to the CBL top. For instance, the entrainment-flux ratio grows substantially (by up to 100%, as shown in Fig. 3b), while the CBL depth increases only slightly (by up to 10%, as shown in Fig. 3a). Note, however, that this slight increase in the CBL depth implies an order-of-one change of the actual entrainment-zone thickness.

Quantitatively, the variation in the entrainment-flux ratio predicted by previous ZOMs is ≈30%–40% (see Fig. 3b). This variation is explained by the different values of *C*_{P} used in the models (cf. Table 1), since the contribution of the surface wind shear to the entrainment flux, represented by the empirical constant *A* in Eq. (13) in Pino et al. (2006), Sun and Xu (2009), and Liu et al. (2016), is around 15% at *z*_{enc}/*L*_{0} ≃ 15 and decreases to less than 5% at *z*_{enc}/*L*_{0} ≃ 25 for the case studied here. This small effect of the surface contribution supports the conclusion in recent work that the surface wind shear affects entrainment mainly indirectly by changing the mean velocity in the mixed layer and thus the velocity jump at the CBL top (Fedorovich and Conzemius 2008). A very good agreement is observed between the recent results of Liu et al. (2016) and results of both the energetics-based model and the geometric-based model with *α* = 1.0, with only ≃5% difference in the entrainment-flux ratio. This agreement indicates that the error of ZOM coefficients derived from numerical simulations starts to be on the order of 5% or less.

The coincidence of the results of the geometric-based model with *α* = 1.0 with those of the energetics-based model once again indicates that the CBL depth predicted by the energetics-based model can be associated to the actual CBL depth that marks the transition from the lower to the upper entrainment-zone sublayer. The CBL depth predicted from the geometric-based model with *α* = 0.8, which corresponds to the height of the minimum buoyancy flux, is ≃5% smaller than the one obtained from the energetics-based model and the one obtained in previous ZOMs (see Fig. 3a). This finding helps explain the reported ≃5% deviation of the zero-order CBL depth from the height of the minimum buoyancy flux in Conzemius and Fedorovich (2007) (see Table 2 in that reference).

We have shown that previous works that appropriately estimate the contribution of the entrainment-zone wind shear to the entrainment flux predict the CBL bulk properties similarly to the energetics-based model and the geometric-based model with *α* = 1. This agreement is remarkable because Liu et al. (2016) have simulated a variety of CBLs in middle latitudes including shear-free, barotropic sheared, and equivalent-barotropic sheared CBLs over an aerodynamically rough surface, while Haghshenas and Mellado (2019) have simulated only shear-free and barotropic sheared CBLs without the Coriolis force over an aerodynamically smooth surface. The observed agreement is, hence, promising in two aspects: first, it confirms that the parameterizations derived in Haghshenas and Mellado (2019) are independent of the surface properties, as they are expressed in terms of the velocity increment at the CBL top. Second, it suggests that they would most likely apply to sheared CBLs with Coriolis force and also to equivalent-barotropic CBLs in which the velocity varies linearly with height in the free atmosphere, although a proof of concept is necessary to draw a definitive conclusion.

The observed agreement between the prediction of the present and previous models might also sound surprising because of the differences in entrainment closures, in particular, differences in the length scale used to estimate the various terms of the TKE budget equation in the entrainment zone. The reason for such an agreement is that under weak- and moderate-shear conditions, the CBL depth (applied in previous work) and the local length scale of the entrainment zone (applied in the present work) are approximately proportional to each other (Haghshenas and Mellado 2019), which results in equally good predictions of the CBL bulk properties from the present and previous models for the moderate-shear conditions considered in this section. Under strong-shear conditions, however, these two length scales are not proportional to each other, and a constant fraction of the CBL depth is not an appropriate proxy of the local length scale of the entrainment zone. This different scaling eventually leads to the emergence of the singularity in models developed in previous work for strong-shear conditions. The models proposed in the present work do not suffer from this limitation.

## 7. Singularity at finite wind strength in previous work

A singularity in previous work takes place when the denominator of Eq. (13) equals zero, that is, when

wherein we have already taken *C*_{T} = 0 (cf. Table 1) and substituted Eq. (20b). Under strong-shear conditions, one can write

where the first approximation follows from Eq. (19) with *α* = 1, and the second approximation holds for Δ*u*^{(0)}/(*N*_{0}*z*_{enc}) ≳ 1.0 with less than 10% error. Substituting the expression for Δ*u*^{(0)} in terms of *h*^{(0)} and *z*_{enc} from Eq. (32) in Eq. (31) and solving the resulting equation for *h*^{(0)} gives

as the singularity condition. For *C*_{P} = 0.43, which is the constant used in Liu et al. (2016), this condition yields *h*^{(0)} ≃ 2*z*_{enc} as the critical value of the CBL depth at which the singularity takes place. This shear condition corresponds to Δ*u*^{(0)}/(*N*_{0}*z*_{enc}) ≃ 1.8, according to Eq. (31). Such an extreme shear condition might happen during the morning in a windy day within the quasi-steady regime, when the well-mixed CBL is still shallow, and the velocity increment at the CBL top is strong.

Another relevant aspect is that, even though the entrainment closures in previous and present work have been derived for the quasi-steady regime, bulk models are likely initialized in atmospheric models with conditions that are far from the quasi-steady regime. This might be problematic in models proposed in previous work because they can also develop a singularity or become unrealistic when the initial conditions are far from the quasi-steady regime even in moderate-shear conditions. This occurs when

according to Eq. (13) with *C*_{T} = 0.

To address this issue, we evaluate our ZOMs and the one by Liu et al. (2016) (as a representative of previous models) with different initial velocity increments at the CBL top, ranging from 5 m s^{−1}, which corresponds to the quasi-steady regime retrieved from Pino et al. (2006), to 8 m s^{−1} in intervals of 1 m s^{−1}. Results are shown in Fig. 4. We observe that all models smoothly relax toward the quasi-steady solutions, with deviations on the order of 10% in the CBL depth or less, except for Liu et al.’s (2016) model with Δ*u*^{(0)} = 8 m s^{−1}, which predicts unrealistically small CBL depths (out of the shown scale). This result illustrates that the singularity or unrealistic results might occur in previous models even for moderate-shear conditions.

## 8. Dependence on environmental conditions of sheared CBL properties

As discussed in the introduction, one main aim of developing bulk models is to investigate the sensitivity of the evolution of CBLs to changes in environmental conditions. Having developed and verified the ZOMs, we explicitly discuss here the dependence of sheared CBL bulk properties on environmental conditions using one of the ZOMs developed in this work, for the purpose of illustration. We employ the energetics-based model and scan the whole parameter space corresponding to midday atmospheric conditions over land (cf. section 2b). As indicated before, dynamical properties of the described CBL, once the initial conditions are sufficiently forgotten and the quasi-steady regime is reached, depend on two nondimensional parameters, namely Fr_{0} and *C*_{D}, and the nondimensional independent variable *z*_{enc}/*L*_{0}. In addition, the nondimensional parameter *φ* characterizes moisture properties.

We consider the case {Fr_{0} = 41, *C*_{D} = 0.002} at *z*_{enc}/*L*_{0} = 40 as reference state and vary one nondimensional parameter at a time. The model is always initialized with the initial condition provided in section 6. Because this initial condition is not the exact one corresponding to the quasi-steady regime for all the parameter space, we illustrate the CBL properties at *z*_{enc}/*L*_{0} = 40, which is sufficiently beyond the initial state of the CBL development, *z*_{enc}/*L*_{0} = 15 (cf. Fig. 4). This assures that the initial conditions have been sufficiently forgotten, and the discussed data are in the quasi-steady regime.

Although the validation of the proposed ZOMs against DNS data has only been done up to Fr_{0} = 25, we consider up to Fr_{0} = 60 in this parametric study. We have already seen in section 6 that the proposed ZOMs compare favorably with previous work where Fr_{0} is substantially larger than 25. To the best of our knowledge, the largest Froude number considered in previous work is 63 in Conzemius and Fedorovich (2006a). Even under this strong-shear condition, these authors report a well-mixed profile of buoyancy and velocity in the mixed layer, and accordingly, the independence between entrainment properties and surface properties. This observation suggests that the proposed closures should remain valid even for such a strong-shear condition, although further analysis should be done to draw a definitive conclusion.

### a. Buoyancy

The dependence of the normalized buoyancy and buoyancy flux on Fr_{0}, *C*_{D}, and *z*_{enc}/*L*_{0} are provided graphically in Fig. 5. This dependence is very small for the buoyancy and buoyancy flux within the mixed layer. The dependence is, however, on the order of one for the entrainment-flux ratio and the buoyancy increment at the CBL top. This behavior illustrates one of the main features of the barotropic CBL, namely, wind shear effects on the CBL structure and buoyancy-related properties remain localized at the CBL top.

We also observe that, with increasing *z*_{enc}/*L*_{0}, the entrainment enhancement due to wind shear diminishes, and hence, wind shear effects on CBL properties reduce and we approach the shear-free limit (cf. section 5). The CBL depth and the mixed-layer buoyancy reduce, respectively, by ≃7% and ≃1%, and the entrainment-flux ratio decreases by ≃25% over the interval of *z*_{enc}/*L*_{0} shown in Figs. 5a and 5b.

Wind shear effects on CBL properties, as expected, grow when the reference Froude number or the surface-drag coefficient increase, because both nondimensional parameters directly and indirectly lead to the larger velocity increment at the CBL top. The CBL depth and the mixed-layer buoyancy increase, respectively, by ≃20% and ≃4% and the entrainment-flux ratio grows by ≃125% for the interval of Fr_{0} shown in Figs. 5c and 5d. For the case Fr_{0} = 60, *C*_{D} = 0.002 and *z*_{enc}/*L*_{0} = 40, the independent variable Δ*u*^{(0)}/(*N*_{0}*z*_{enc}) is approximately 0.8, which is much below the critical condition to observe the singularity in previous ZOMs. Further analysis (not shown) indicated that the critical condition, Δ*u*^{(0)}/(*N*_{0}*z*_{enc}) ≃ 1.8, takes place, for instance, for Fr_{0} = 80 and *C*_{D} = 0.005 in the early state of the CBL development.

The effect of Fr_{0} and *C*_{D} on the CBL evolution differs in that the CBL depth grows with *C*_{D} asymptotically toward a finite value (see Figs. 5e and 5f), whereas the growth of the CBL with Fr_{0} is unbounded (see Figs. 5c and 5d). The reason is that growing *C*_{D} with fixed Fr_{0} causes wind shear effects to emerge earlier (at a smaller *z*_{enc}/*L*_{0}) because the velocity increment at the CBL top increases fast. Wind-shear effects are, however, limited since the reference Froude number, or equivalently the velocity in the free atmosphere, is fixed.

### b. Moisture

To address the dependence of moisture properties of the sheared CBL on environmental conditions, we first consider *φ* = 0, which corresponds to the pure-drying regime. Figure 6 shows graphically the dependence of the normalized moisture and moisture flux on Fr_{0}, *C*_{D}, and *z*_{enc}/*L*_{0}. There are two features worth mentioning in this figure. First, given that the free atmosphere is dry and the wind shear enhances entrainment, the CBL dries more when the Froude number or the surface-drag coefficient increases (see Figs. 6c and 6e). The entrainment enhancement due to the wind shear, however, diminishes as the CBL grows, and hence, the sheared CBL dries less as *z*_{enc}/*L*_{0} increases (see Fig. 6a). Second, wind shear effects on the mixed-layer specific humidity, $qml\u2061(0)$, are much larger than their effects on mixed-layer buoyancy, $bml\u2061(0)$ (cf. Figs. 5 and 6). The reason is that *φ* = 0 corresponds to *F*_{q,0} = 0, that is, there is no surface flux of moisture but only the entrainment flux, so the entrainment enhancement due to the wind shear is more relevant in the moisture field than in the buoyancy field.

Figures 7a and 7b illustrate graphically the dependence of moisture bulk properties on the flux-ratio parameter for Fr_{0} = 41, *C*_{D} = 0.002, and *z*_{enc}/*L*_{0} = 40. For *φ* ≲ 1.2, the entrainment flux of drying air is larger than the surface flux of moisture, causing the CBL to dry. The surface flux becomes equal to the entrainment flux for *φ* ≃ 1.2 (Fig. 7b); this condition corresponds to the boundary between drying to moistening regime. For *φ* ≳ 1.2, the surface flux of moisture is dominant over the entrainment flux of drying air that causes the CBL to moisten.

The crossover value, $\phi cr\u2061(0)\u22431.2$, for the considered sheared CBL is slightly larger than $\phi cr\u2061(0)\u22431.17$ for the shear-free limit, since entrainment of dry air increases with the wind shear. Figure 7c illustrates this behavior more clearly, as the critical flux-ratio parameter is enhanced with increase of the Froude number. This enhancement, however, diminishes as the CBL grows.

## 9. Summary and conclusions

Two zero-order bulk models (ZOMs) with different entrainment closures have been developed for a cloud-free barotropic convective boundary layer (CBL) that grows into a linearly stratified atmosphere. In the first one, the negative and positive areas of the buoyancy flux were assumed to match between the model and the actual CBL. In the second one, the CBL depth was the variable chosen to match between the model and the actual CBL. Nonsingular parameterizations for these properties derived from direct numerical simulation (DNS) in Haghshenas and Mellado (2019) have been employed as the entrainment closure equation in each model. We referred to these models as energetics- and geometric-based models, respectively. The proposed ZOMs have been verified by comparing their predictions with the DNS data.

Under moderate-shear conditions, predictions from previous ZOMs that appropriately estimated the contribution of the entrainment-zone shear in the entrainment closure equation agree well with predictions of the energetics-based model. Under strong-shear conditions, however, previous ZOMs can develop a singularity at finite wind strength that has often been reported as a major limitation for their use. Increasing the complexity of the bulk model to the first-order or higher-order models does not remove this singularity as long as the CBL depth is used to estimate the various terms of the TKE equation in the entrainment zone. The proposed ZOMs in the present work are free from this potential singularity because they consider the local length scale of the entrainment zone, instead of a fraction of the CBL depth, when estimating the various terms of the TKE budget equation. We discussed the potential singularity in the entrainment closure equation of previous work analytically and numerically. Using the parameterizations for different CBL properties, we have shown that the singularity takes place under strong-shear conditions when the sheared CBL depth becomes nearly 2 times larger than the encroachment length scale, or equivalently when Δ*u*^{(0)}/(*N*_{0}*z*_{enc}), as the independent variable that characterizes wind shear effects, equals 1.8. In addition, we have shown that considering initial conditions far away from the quasi-steady regime also leads to the singularity or unrealistic results in the models proposed in previous work even for moderate-shear conditions.

We have also explained how the different zero-order CBL depths relate to actual CBL depths, which might become important when the bulk model is intended to include more complexity like cloud formation. We developed the geometric-based model to precisely address this issue. We considered three options for the CBL depth in the geometric-based model, namely, the height of the minimum buoyancy flux, the height that marks the transition from the lower to the upper entrainment-zone sublayer, and the height of the maximum buoyancy gradient. These heights differ by few hundred meters under typical midday atmospheric conditions over land. Predictions of the geometric-based model suggested that the CBL depths in the energetics-based model and also models in previous work correspond better to the height that marks the transition from the lower to the upper entrainment-zone sublayer, rather than the height of the minimum buoyancy flux. This finding helps explain the approximately 5% deviation of the zero-order CBL depth from the height of the minimum buoyancy flux reported in Conzemius and Fedorovich (2007).

An important conclusion of this study is that the zero-order bulk model, despite its simplicity, can appropriately represent bulk properties of sheared CBLs, meaning that a finite transition layer between the mixed layer and the free atmosphere is not explicitly required. This is because the relevant shear-induced features of the actual entrainment zone are considered in the entrainment closure. If needed, the vertical structure of the actual entrainment zone of the sheared CBL can be constructed a posteriori using the zero-order CBL depth predicted from any of the ZOMs and using the relationships between the zero-order CBL depth and various actual CBL depths provided in Haghshenas and Mellado (2019).

## Acknowledgments

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). Funding was provided by the Max Planck Society through its Max Planck Research Groups program. Primary data and scripts used in the analysis and other supporting information that may be useful in reproducing the author’s work are archived by the Max Planck Institute for Meteorology and can be obtained by contacting publications@mpimet.mpg.de.

### APPENDIX A

#### The Surface-Drag Coefficient

The surface-drag coefficient is a complicated function of the aerodynamic roughness length of the surface *z*_{0}, the surface-layer depth *h*_{s1}, and the Obukhov length *L*_{Ob}. Such expression is derived from the Monin–Obukhov similarity theory as (see, e.g., Garratt et al. 1982)

where *κ* is the von Kármán constant, and *h*_{s1} is usually considered as 10% of the CBL depth (Tennekes 1973b). Here we consider the height that separates the lower from the upper entrainment-zone sublayer as the CBL depth. The function *ψ*_{m}(*h*_{s1}/*L*_{Ob}) is the Businger–Dyer (Kansas) formulation defined as (see, e.g., Grachev et al. 2000)

where *x* = [1 − 16*h*_{s1}/*L*_{Ob}]^{1/4}. For the aerodynamically smooth surface in which the viscous sublayer is deeper than surface roughness protuberances, the surface roughness is proportional to $\nu /u*$ (Hinze 1975). Our analysis provides evidence that the surface roughness in the weakly to strongly unstable CBL with the Reynolds number considered in Haghshenas and Mellado (2019) is well obtained by $z0\u22430.13\nu /u*$ (see Fig. A1a). The obtained proportionality constant is in a good agreement with the one corresponding to the canonical channel flow (Pope 2000).

### APPENDIX B

#### Parameterization for the Area Ratio of the Negative and Positive Buoyancy Flux in the Actual CBL

To derive the parameterization for the area ratio of the negative and positive buoyancy flux, we consider Eq. (5.4) of Haghshenas and Mellado (2019):

which implies that the shear-generated turbulence within the entrainment zone accounts for the entrainment enhancement in sheared CBLs. Here *z*_{i,0} is the reference height at which the buoyancy flux becomes zero. Subscript *c* denotes the shear-free limit, primes indicate turbulent-fluctuation fields, and *w* is vertical velocity. In agreement with previous work, the negative area of the buoyancy flux in the shear-free limit is scaled by the convective scales as $\u2212\u222bzi,0z\u221e\u2329b\u2032w\u2032\u232ac\u2009dz\u22430.02B0zenc$, where the coefficient of proportionality is obtained from the DNS data. Thus, the area of the negative buoyancy flux in sheared CBLs can be approximated as

using Eq. (B1) and the well-known scaling argument for the integral of the shear production term (see, e.g., Boers et al. 1984). Here *w*_{e} ≡ *dz*_{i,f}/*dt* is the growth rate of the CBL depth, where *z*_{i,f} is the height of the minimum buoyancy flux. To avoid the poor statistical convergence associated with determining *w*_{e} by taking the time derivative of *z*_{i,f}, we use the approximation provided in Haghshenas and Mellado (2019) as *w*_{e} ≃ (*z*_{i,f}/*z*_{enc})*dz*_{enc}/*dt*. The DNS data support the ansatz in Eq. (B2) and show *c*_{1} ≃ 0.09.

The fact that wind shear only modifies the vertical structure of the entrainment-zone indicates that the positive area of the buoyancy flux in sheared CBLs, consistent with shear-free CBLs, is scaled as

where the coefficient of proportionality is obtained from the DNS data. The parameterization for the area ratio then reads as

We obtain a value ≃0.044 for the area ratio in the shear-free limit that corresponds to the zero-order entrainment-flux ratio ≃0.21 in the shear-free limit, according to Eq. (17) (see Fig. A1b). Figure A1c supports the derived parameterization for the area ratio in the sheared CBLs with DNS data, showing that the dependence of this parameterization on the Reynolds number and also the dependence of the proportionality coefficient *c*_{1} on the choice of the reference definition of the CBL depth in *w*_{e} is smaller than the achieved statistical convergence. This sensitivity analysis allows us to consider Eq. (B4) as the entrainment closure in the energetics-based model, regardless of knowing a priori to which definition of the actual CBL depth, the modeled CBL depth can be associated.

### APPENDIX C

#### Geometric-Based Model Corresponding to the Height of the Maximum Buoyancy Gradient

As explained in the main text, one can also develop a geometric-based model with the modeled CBL depth chosen to be equal to the actual height of the maximum buoyancy gradient. The corresponding closure equation reads (Haghshenas and Mellado 2019)

where

Term I in Eq. (C1) corresponds to the depth of the actual mixed layer, term II corresponds to the thickness of the lower entrainment-zone sublayer, and term III corresponds to half of the thickness of the upper entrainment-zone sublayer. Even though employing the geometric-based model with the aforementioned closure is straightforward, it is more convenient to use those two geometric-based models explained in the main text and construct the height of the maximum buoyancy gradient from the results of those models and Eq. (C1). It is also noted that the actual upper height of the CBL, where the buoyancy flux is ≃15% of the minimum, can also be constructed a posteriori using Eq. (C1) by multiplying the third term by 2.

## REFERENCES

## Footnotes

^{a}

Current affiliation: Division of Aerospace Engineering, Department of Physics, Universitat Politècnica de Catalunya, Barcelona, Spain.

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