# Well-Mixed Boundary-Layer Top Entrainment Instability: Hydrodynamic Analysis

Jun-Ichi Yano

<figure><img src="image_1.png" /></figure>

► To cite this version:

Jun-Ichi Yano. Well-Mixed Boundary-Layer Top Entrainment Instability: Hydrodynamic Analysis. Journal of the Atmospheric Sciences, 2022, 10.1175/JAS-D-21-0246.1. [https://hal.science/hal-03716668v1](https://hal.science/hal-03716668v1)

HAL Id: hal-03716668

https://hal.science/hal-03716668v1

Submitted on 11 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Well-Mixed Boundary Layer-Top Entrainment Instability: Hydrodynamic Analysis

JUN-ICHI YANOª

ª CNRM, UMR 3589 (CNRS), Météo-France, Toulouse, France

(Manuscript received 17 September 2021, in final form 24 March 2022)

**ABSTRACT:** The present study shows by a linear hydrodynamic stability analysis that an unstable mixed-layer deep circulation can be generated in the dry convective well-mixed layer by the entrainment from the top. The newly identified instability arises under the two competing processes induced by the top entrainment: the destabilization by generating thermal perturbations and the damping by mechanical mixing. The former and the latter, respectively, dominate over the other in the limits of large and small scales. As a result, the instability is realized at the horizontal scales larger than the order of the mixed-layer depth (ca. 1 km), and the time scale for the growth is about 1 day. This study has been motivated from a question of whether the cloud-top entrainment instability (CTEI) can induce a transition of the stratocumulus-topped well-mixed boundary layer into trade cumulus. The present study intends to extend the previous studies based on the local parcel analyses to a full analysis based on the hydrodynamics. Unfortunately, being based on a dry formulation, the present result does not apply directly to the CTEI problem. Especially, the evaporative cooling is totally neglected. Nevertheless, the present result can still be applied to moist systems, to some extent, by redefining certain terms in the formulation.

**KEYWORDS:** Atmosphere; Buoyancy; Convective clouds; Entrainment; Boundary layer; Stratiform clouds

# 1. Introduction

The cloud-top entrainment instability (CTEI; Deardorff 1980) is considered a major potential mechanism for the transition of the stratocumulus to the trade cumulus over the marine subtropics [cf. Stevens (2005) as an overview]. The basic mechanism of CTEI resides on a possibility that an environmental air entrained into the cloud from the top can be dry enough so that its mixing with the cloudy air leads to evaporation of the cloud water, and induces a sufficient negative buoyancy, leading to further entrainments of the environmental air from the cloud top. The process is expected to finally lead to a transition of stratocumulus into cumuli. A critical review of this process is provided by Mellado (2017), with the review itself even refuting CTEI as further discussed in the end in section 5. Bretherton and Wyant (1997) and Lewellen and Lewellen (2002) propose the decoupling as an alternative theoretical possibility.

However, the existing literature examines CTEI, mostly, in terms of a local condition, such as a buoyancy anomaly at the cloud top (inversion height). Such a parcel-based analysis leads to a criterion for instability in terms of a sign of buoyancy (e.g., Deardorff 1980; Randall 1980; MacVean and Mason 1990; Duynkerke 1993). This type of approaches does not provide a full dynamical picture of the instability, including a quantitative estimate of a growth rate as a function of a horizontal scale (or a wavenumber), and a spatial structure of a preferred instability mode.

The qualitative nature of the existing criteria for CTEI makes it also difficult to test these criteria observationally (cf. Albrecht et al. 1985; Albrecht 1991; Kuo and Schubert 1988; Stevens et al. 2003; Mathieu and Lahellec 2005; Gerber et al. 2005, 2013, 2016). Most fundamentally, a finite time would be required for CTEI to realize. Unfortunately, bulk of existing theories does not tell how long we have to wait to observe CTEI.

A fundamental limitation of existing CTEI studies arises from a fact that these analyses concern only with a sign of a local buoyancy (or vertical eddy buoyancy flux), without properly putting it into a framework of the hydrodynamic instability (cf. Drazin and Reid 1981). Such a dynamically consistent theoretical analysis of the instability couples a given local instability with a full hydrodynamics. It is a standard approach in the midlatitude large-scale dynamics to interpret the synoptic cyclones in this manner in terms of the baroclinic instabilities (cf. Hoskins and James 2014). In the author's knowledge, a hydrodynamic stability analysis is still to be performed for CTEI, probably an exception of Mellado et al. (2009; cf. section 2d). That is the basic approach of the present study.

The hydrodynamics stability-analysis method adopted here treats the evolution of the height of the inversion at the top of the mixed layer explicitly with time so that, in principle, its evolution until an ultimate transform into a cumulus regime can be evaluated. For preparing a way for such full analyses, the present study introduces a linear analysis method by taking the dry atmosphere as a demonstrative example. Thus, an important purpose of the study is to show how dynamically consistent instability analyses can be performed in problems of atmospheric boundary layers. The author expects that more studies will follow along this line for better elucidating the dynamics of the cloud-topped boundary layers. Importantly, the study is going to show that even in absence of

<footnote>ª Denotes content that is immediately available upon publication as open access.</footnote>

<footnote>Corresponding author: Jun-Ichi Yano, jun-ichi.yano@cnrs.fr</footnote>
evaporative cooling, the mixed layer can be destabilized by the entrainment from the top.

It may be considered questionable to perform a linear stability analysis on a fully turbulent system such as the boundary layers. To circumvent this difficulty, the present study assumes that the main role of fully developed convective turbulence is to maintain a vertically well-mixed state of the boundary layer, and that a boundary layer–deep explicit perturbation can be considered as a linear superposition on the mean state maintained by these turbulent flows, however, without explicitly taking into account of the latter.

For this reason, especially, for describing the buoyancy in the well-mixed layer, only an equation averaged over the mixed layer is considered. Note that as a consequence of active vertical mixing, the well-mixed layer is neutrally stratified, in average; thus, no linear buoyancy anomaly can be generated by linear perturbation flows. In this manner, the role of convection remains completely implicit in the present study. However, it is important to keep in mind that the entrainment, that drives the instability, is also driven by convection. Thus, the instability considered herein is ultimately driven by convection.

More specifically, the present study examines a perturbation growth of a mixed-layer deep circulation. This approach is contrasted with some studies, dealing CTEI primarily as a process of generating kinetic energy for smaller-scale eddies, that directly contribute to vertical eddy transport at the top of the well-mixed layer associated with entrainment (e.g., Lock and MacVean 1999; Katzwinkel et al. 2012). An overall approach of the present study may be compared with that for the mesoscale entrainment instability by Fiedler (1984; see also Fiedler 1985; Rand and Bretherton 1993). As a major difference, the entrainment induces negative buoyancy by a downward displacement of inversion in the present study, whereas Fiedler considered an enhancement of cloudy-air positive buoyancy by entrainment of stable upper-level air. At a more technical level, the present study considers a change of the buoyancy jump crossing the inversion with time, but fixing the entrainment rate. In Fiedler (1984), in contrast, the main role of the inversion jump is to constrain the entrainment rate.

More general words may be required for some readers who are not familiar with the basics of the hydrodynamic stability analysis. To perform a hydrodynamic stability analysis in a general manner, certain simplifications are always necessary. In this respect, the hydrostatic stability does not pursue any “realism” in the same sense as with both operational and research models widely available today. However, our experience says that those simplified theoretical studies provide useful, and often quantitative information on the process in concern (cf. Pedlosky 1987; Hoskins and James 2014).

The formulation, that couples a standard mixed-layer description with a full hydrodynamics, is introduced in the next section. A perturbation problem is developed in section 3, and some simple solutions are presented in section 4. The

paper concludes with the discussion in the last section. An alternative formulation is considered separately in appendix A.

## 2. Formulation

A dry well-mixed boundary layer is considered. Nevertheless, as we remark from time to time, to some extent, the formulation may also be, at least, conceptually applied to the stratiform-topped mixed layer.

### a. Motivations

Essence of CTEI is that a mixing of the free-troposphere air from the above with a cloudy air within stratocumulus leads to evaporation of cloud water due to a dry and relatively high temperature of the entrained free-atmospheric air, but the evaporative cooling, in turn, makes the entrained air colder than the surrounding stratocumulus-cloud air, leading to a convective instability that drives the evaporated mixed air farther downward (Deardorff 1980; Randall 1980). Though less frequently considered, a possible reverse process is an intrusion of the cloudy air from the stratocumulus cloud into the free troposphere (e.g., MacVean and Mason 1990; Duynkerke 1993). In this case, when the detrained air is moist enough, it can be more buoyant than the environment due to the virtual effect. Buoyancy induces a further ascent, the ascent leads to adiabatic cooling, the cooling may lead to further condensation of water vapor, and resulting condensative heating can drive the cloudy air farther upward.

Being motivated by investigating this type of instability fully dynamically, first of all, the present study explicitly describes the deformation of the inversion height with time, associated with the entrainment of warm and drier air from the free atmosphere above. The resulting deformation may ultimately lead to transform into a cumulus regime. We will consider the associated processes under a drastically simplified dry mixed-layer formulation. In spite of these drastic simplifications, we somehow recover some basic features of the CTEI just described. The drastic simplification facilitates the analysis of the coupling of these processes with a full dynamics in a form of hydrodynamic stability analysis.

A simple dry mixed-layer formulation is introduced in the next two subsections. It is coupled with a full hydrodynamics introduced in sections 2d and 2e.

### b. A mixed-layer formulation for the buoyancy

We consider a well-mixed boundary layer with a depth (inversion height) *h*. The basic model configuration is shown in Fig. 1. As the most drastic simplification here, we adopt the standard formulation for the dry boundary layer, in which only the buoyancy *b* vertically averaged over the full mixed layer is considered. This approach is well justified for the dry boundary layer, because the buoyancy phenomenologically is known to be vertically well mixed, as also suggested in Fig. 1.

However, this assumption clearly breaks down for the cloud-topped well-mixed boundary layer. Under standard
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic configuration of the model.</figcaption></figure>

formulations (e.g., Deardorff 1980; Schubert et al. 1979), the buoyancy anomaly is expressed by a linear relationship with the two conservative quantities, say, the equivalent potential temperature and the total water, which are expected to be vertically well mixed. However, the buoyancy is not expected to be vertically well mixed, because the coefficients for this linear relationship are height dependent [cf. Eq. (3.15) of Schubert et al. (1979), Eqs. (15) and (22) of Deardorff (1976)].

More seriously, although the buoyancy in the dry convective layer may vertically be well mixed in undisturbed state, once a perturbation is applied, a nonvanishing buoyancy disturbance is generated, which we will consider explicitly in the following. Under the latter situation, the buoyancy homogeneity assumption no longer applies. Thus, in more general situations with presence of clouds as well as disturbances, the buoyancy is no longer vertically homogeneous distributed. Nevertheless, as going to be shown in the following, a self-consistent formulation of the problem is still possible in terms of the vertically averaged buoyancy $\langle b \rangle$ even under an explicit presence of the buoyancy perturbations.

Thus, we describe the buoyancy, $b$, under its vertical average,

$$h\left(\frac{\partial}{\partial t} + \langle u \rangle \frac{\partial}{\partial x}\right)\langle b \rangle = \overline{w'b'}_0 - \overline{w'b'}_{-} - hQ_R, \quad (2.1)$$

by following a standard formulation for the well-mixed boundary layer [e.g., Eqs. (3.1) and (3.3) of Schubert et al. 1979; Eq. (21) of Stevens 2006]. Here, the bracket, $\langle \cdot \rangle$, designates a vertical average over the well-mixed layer. As already emphasized above, a deviation from a vertical average may actually exist, but we simply neglect these contributions in the formulation. A two-dimensional configuration has been assumed for a sake of simplicity. A full three-dimensional analysis would be substantially more involved without any practical benefits.

Here, we have introduced the variables as follows: $t$ the time, $x$ a single horizontal coordinate considered, $u$ the horizontal wind velocity, $\overline{w'b'}$ the vertical buoyancy flux with the subscripts 0 and - designating the values at the surface and at the level just below the inversion (i.e., $h_-$), respectively; $Q_R$ is the loss of buoyancy due to the radiative cooling over the well-mixed layer. Note that the buoyancy flux is discontinuous over the inversion associated with a discontinuity of the buoyancy (cf. Fig. 1).

Under a standard formulation [cf. Eqs. (1) and (2) of Deardorff 1980], the vertical eddy flux just below the inversion level may be expressed in terms of the entrainment rate, $w_e$ ($> 0$), and a jump, $\Delta b = b_+ - \langle b \rangle$, of the buoyancy over the inversion (with $b_+$ the free troposphere value at $z = h_+$) as

$$\overline{w'b'}_{-} = -w_e \Delta b. \quad (2.2)$$

Note that the radiative cooling over the inversion layer can also be included as a part of the buoyancy flux on the left-hand side above.

Here, the standard CTEI criteria (Deardorff 1980; Randall 1980) require $\overline{w'b'}_{-} > 0$ or $\Delta b < 0$. When this condition is satisfied, the induced negative buoyancy is expected to induce further cloud-top entrainment, which induces further negative buoyancy; that is an essence of CTEI as described in the last subsection. Extensive CTEI literature focuses on defining this condition carefully due to a subtle difference between the inversion buoyancy jump and an actual buoyancy anomaly generated by a cloud-top mixing (cf. Duynkerke 1993). The dry formulation herein does not take into account those aspects in any direct manner. However, a similar tendency somehow arises for perturbation variables as seen immediately below.

### c. Perturbation formulation and instability mechanism

In the following, we only consider the perturbations by setting

$$\begin{align*} h &= \bar{h} + \eta, \\ \langle b \rangle &= \langle \bar{b} \rangle + \langle b \rangle', \\ b_+ &= \bar{b}_+ + b'_+, \end{align*}$$

where a bar and a prime designate equilibrium and perturbation values, respectively. An exception to this rule is the perturbation inversion height designated as $\eta$.

By linearizing Eq. (2.1), the perturbation equation for the buoyancy is given by

$$\bar{h}\left(\frac{\partial}{\partial t} + U_1 \frac{\partial}{\partial x}\right)\langle b \rangle' = -\overline{w'b'}_{-} - \eta Q_R. \quad (2.3)$$

Here, for simplicity, we have assumed that there is no perturbation surface flux, i.e., $\overline{w'b'}_0 = 0$. It is also assumed that $Q_R$ do not change by perturbations. The basic instability mechanism of this system is inferred by multiplying $\langle b \rangle'$ on both sides of this equation; thus, we obtain an equation for the potential-energy conservation law:
$$ \bar{h} \left( \frac{\partial}{\partial t} + U_1 \frac{\partial}{\partial x} \right) \frac{\langle b \rangle'^2}{2} = - \langle b \rangle' (\overline{w' b'_{-}}) - Q_R \langle b \rangle' \eta. \quad (2.4) $$

It is seen that the potential energy can be generated when

$$ - \langle b \rangle' (\overline{w' b'_{-}}) > 0, \quad (2.5) $$

i.e., the downward entrainment flux $-(\overline{w' b'_{-}})$ is positively correlated with the buoyancy perturbation, $\langle b \rangle'$. This is the basic mechanism of the entrainment-induced instability to be investigated in the following. Thus, an essential ingredient to induce an instability is to maintain the perturbation flux $-(\overline{w' b'_{-}})$ to be nonvanishing, regardless of additional assumptions to be introduced.

Here, there is no simple closed expression for the entrainment rate $w_e$; thus, a simple assumption to satisfy the instability condition is to set the perturbation entrainment as $w'_e = 0$ so that the perturbation flux at the mixed-layer top reduces to

$$ \overline{w' b'_{-}} = - w_e \Delta b' \quad (2.6) $$

from Eq. (2.2). Adding a term with $w'_e$ to the above only makes the following analysis more involved without adding anything substantial. An oppositive extreme of possible options is to assume the perturbation buoyancy flux at the mixed-layer top to vanish, i.e., $\overline{w' b'_{-}} = 0$, and to take into account of the resulting perturbation entrainment $w'_e$. This case is considered separately in appendix A, by following a suggestion of one of anonymous reviewers. See the next subsection for the discussions on the basic state, $\bar{h}$, and $\langle \bar{b} \rangle$.

A perturbation on the buoyancy jump may be given by

$$ \Delta b' = b'_{+} - \langle b \rangle'. \quad (2.7a) $$

Here, we assume that the perturbation to the free-atmosphere buoyancy is solely induced by a displacement $\eta$ of the interface (cf. section 2e). Thus, with the help of the Taylor expansion,

$$ b'_{+} = b|_{z=\bar{h}_{+}+\eta} - b|_{z=\bar{h}_{+}} = \left(\frac{d\bar{b}}{dz}\right)\eta. \quad (2.7b) $$

It shows that a positive displacement, $\eta > 0$, of the inversion induces a positive buoyancy perturbation, $\Delta b' > 0$, and vice versa. Downward extrapolation of the Taylor expansion could be problematic, because the free-atmosphere profile may not simply follow downward as the inversion height is distorted downwards. It may be more likely that the inversion layer simply thickens. However, the proposed formula (2.7b) is consistent even under this situation by following an explicit consideration of the finite-depth effect by Betts [1974, see especially his Eq. (5)].

By substituting Eq. (2.7b) into the definition of the buoyancy-jump perturbation, (2.7a), we obtain

$$ \Delta b' = \left(\frac{d\bar{b}}{dz}\right)\eta - \langle b \rangle' \quad (2.8) $$

[cf. Eq. (2) of Tennekes and Driedonks 1981]. The first term on the right-hand side of Eq. (2.8) shows that a displacement of the inversion induces a buoyancy anomaly, which further enhances the displacement by accelerating the vertical velocity in the displacement direction. This positive feedback chain leads to an instability. To see this process more explicitly, the buoyancy equation must be coupled with a hydrodynamic system, as going to be introduced in next two subsections. The second term on the right-hand side, in turn, simply states how a buoyancy perturbation $\langle b \rangle'$ of the mixed layer modifies the buoyancy jump $\Delta b'$ at the inversion. It leads to a damping tendency as seen immediately below.

Substitution of Eq. (2.8) into Eq. (2.6) reduces Eq. (2.3) into

$$ \left[ \bar{h} \left( \frac{\partial}{\partial t} + \langle u \rangle \frac{\partial}{\partial x} \right) + w_e \right] \langle b \rangle' = \alpha \eta, \quad (2.9a) $$

where

$$ \alpha = w_e \left( \frac{d\bar{b}}{dz} \right) - Q_R. \quad (2.9b) $$

Measures a feedback of the inversion height anomaly $\eta$ on the buoyancy anomaly $\langle b \rangle'$. Here, we expect $\alpha > 0$. As already discussed above, the first term on the right-hand side in Eq. (2.9b) shows that displacements of the inversion tend to enhance the buoyancy perturbation. The second term on the right-hand side is a negative radiative feedback, arising from the fact the total radiative cooling rate of the mixed layer changes by the inversion-height displacement. Negative feedback of radiation on CTEI has been pointed out by, e.g., Moeng and Schumann (1991) and Moeng et al. (1995).

Equation (2.9a) contains the two competitive processes arising from the mixed-layer-top entrainment: the first is a mechanical mixing as its direct consequence, that leads to a damping, as indicated by the last term on the left-hand side. The second is a consequence of the inversion-height displacement, as seen on the right-hand side, which may induce instability. The first effect is independent of scales, whereas the second depends on scales, as further discussed with Eq. (3.8a) below. The scale dependence of the latter leads to a scale dependence of the instability growth as will be shown in section 4.

### d. Basic state

To introduce a hydrodynamics, we adopt a two-layer system with constant densities (cf. Fig. 1), closely following a standard formulation for the analysis of the Kelvin–Helmholtz instability as presented, e.g., in chapter 4 of Drazin and Reid (1981). The first layer with a density $\rho_1$ represents the well-mixed layer below, and the second with a density $\rho_2$ the free troposphere above. To some extent, this formulation can be considered a local description of the dynamics around the top of the well-mixed layer (the inversion height), $z = h$, although the bottom (surface: $z = 0$) and the top ($z \to +\infty$) boundary conditions are considered explicitly in the following. A height dependence of the density can be introduced to this system, and so long as the density-gradient scale is much larger than a vertical scale of the
interest, the given system is still considered a good approximation. Under this generalization, for the most parts in the following, the density values $ρ_1$ and $ρ_2$ refer to those at the inversion height, $z = h$. We also assume that the horizontal winds, given by $U_1$ and $U_2$, are constant with height in each layer. Thus, we may reset $U_1 = \langle u \rangle$ in the formulation of the last subsection.

Here, an assumed sharp interface is a necessary simplification for treating the essential features of the CTEI in lucid manner, although both recent observational (Lenschow et al. 2000; Katzwinkel et al. 2012) and modeling (Moeng et al. 2005) studies show that the inversion actually constitutes a finite-depth layer with rich morphologies. Mellado et al. (2009) consider a Rayleigh-Taylor instability problem by inserting a positive density anomaly over this thin inversion layer. Their study may be considered an extension to three layers of the present formulation. However, in contrast to the present study, the fluid density is assumed a passive scalar and no possibility of its change associated with the entrainment.

We assume that the basic state is under a hydrostatic balance; thus, the pressure field is given by

$$p = \begin{cases} p_h - \rho_1 g (z - h), & 0 \le z \le h, \\ p_h - \rho_2 g (z - h), & z > h, \end{cases} \quad (2.10)$$

where $p_h$ is a constant pressure value at the inversion height.

The inversion height $h$ is described by [cf. Eq. (4) of Stevens (2002), Eq. (31) of Stevens (2006)]

$$\left(\frac{\partial}{\partial t} + u_j \frac{\partial}{\partial x}\right)h = w + w_e \quad (2.11)$$

for both layers with $j = 1, 2$. Its steady basic state $\bar{h}$ is defined by the balance:

$$\bar{w} + w_e = 0. \quad (2.12)$$

Here, $\bar{w}$ is a height-dependent background vertical velocity defined below. When $\bar{w} < 0$, we identify an equilibrium state at a certain height. Especially, when $\bar{w}$ is a monotonous function of the height, the equilibrium inversion height is unique. On the other hand, when $\bar{w} > 0$, there is no equilibrium height for the inversion; thus, we may generalize above as

$$\dot{\bar{h}} = \bar{w} + w_e$$

with the rate $\dot{\bar{h}}$ of change of the basic inversion height. In the latter case, the perturbation is applied against an unsteady state with $\dot{\bar{h}} \neq 0$. In the following, we further assume a constant background divergence $D$; thus,

$$\bar{w} = -Dz.$$

Finally, the basic state $\langle \bar{h} \rangle$ for the mixed-layer buoyancy is defined from Eq. (2.1) assuming a steady and homogeneous state. It transpires that the basic state is obtained from a balance between three terms on the right-hand side.

Unfortunately, deriving the basic-state explicitly for $\langle \bar{h} \rangle$ is rather involved with a need of specifying the dependence of $\overline{w'b'_{0}}$ and $Q_R$ on $\langle b \rangle$ (i.e., specifications of physical processes). Here, we do not discuss this procedure, because this problem is, for the present purpose, circumvented by simply prescribing a mean state $\langle \bar{h} \rangle$. As it turns out, the value of $\langle \bar{h} \rangle$ does not play any direct role in the instability problem.

### e. Perturbation dynamics

For developing a perturbation problem, we assume that the perturbations satisfy the following boundary conditions (with the prime suggesting perturbation variables):

$$(i) \quad u' \to 0 \text{ as } z \to +\infty, \quad (2.13a)$$

$$(ii) \quad w' = 0 \text{ at the bottom surface } z = 0, \quad (2.13b)$$

$$(iii) \quad the pressure is continuous by crossing the inversion, $z = h$; thus,

$$p'_{1} - \rho_{1}g\eta = p'_{2} - \rho_{2}g\eta \quad (2.13c)$$

at $z = \bar{h}$ after linearization. Furthermore, we may note that the perturbation equation for the inversion height is given by

$$\left(\frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x}\right)\eta = -D\eta + w'\bigg|_{z=\bar{h}} \quad (2.13d)$$

for $j = 1$ and $2$.

The perturbation equations for the dynamics are given by

$$\left(\frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x}\right)w'_{j} = -\frac{1}{\rho_j} \frac{\partial p'_{j}}{\partial z} + b'_{j}, \quad (2.14a)$$

$$\left(\frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x}\right)u'_{j} = -\frac{1}{\rho_j} \frac{\partial p'_{j}}{\partial x}, \quad (2.14b)$$

for $j = 1$ and $2$. Here, the buoyancy perturbation equation for the lower layer ($j = 1$) is given by setting $b'_1 = b'$ in Eq. (2.9a). In the upper layer ($j = 2$), we simply set $b'_2 = 0$. Nonvanishing buoyancy perturbation in the upper layer (free troposphere) would contribute to the gravity wave dynamics (cf. Fiedler 1984). We simply neglect this contribution.

We further introduce the perturbation vorticity $\zeta'$ and streamfunction $\psi'$ so that

$$\zeta' = \frac{\partial u'}{\partial z} - \frac{\partial w'}{\partial x} = \nabla^2 \psi', \quad (2.15a)$$

$$w' = -\frac{\partial \psi'}{\partial x}, \quad (2.15b)$$

$$u' = \frac{\partial \psi'}{\partial z}, \quad (2.15c)$$

and for a later purpose, it is useful to note from Eqs. (2.15a) and (2.15b),

$$\frac{\partial \zeta'}{\partial x} = - \nabla^2 w'. \quad (2.15d)$$
The perturbation equations for the vorticity in both layers are obtained from Eqs. (2.14a) and (2.14b):

$$ \left(\frac{\partial}{\partial t} + U_1 \frac{\partial}{\partial x}\right) \zeta_1' = - \frac{\partial b_1'}{\partial x}, \qquad (2.16a) $$

$$ \left(\frac{\partial}{\partial t} + U_2 \frac{\partial}{\partial x}\right) \zeta_2' = 0. \qquad (2.16b) $$

### 3. Stability analysis

The perturbation problem is solved for the dynamics and the buoyancy separately in the following two subsections. Each leads to an eigenvalue problem.

#### a. Dynamics problem

The solution for the upper layer is obtained in a relatively straightforward manner. From Eq. (2.16b), we find an only solution satisfying the condition of the vanishing perturbation flow toward $z \to +\infty$ (2.13a) is $\zeta_2' = 0$; thus,

$$ \nabla^2 \psi_2' = 0, $$

whose solution consistent with the boundary condition (2.13a) is

$$ \psi_2' = \hat{\psi}_2 e^{ikx-k(z-\bar{h})+\sigma t}. $$

Here, both the horizontal and the vertical scales are characterized by a single parameter $k$, which is assumed to be positive; $\sigma$ is a growth rate. It immediately follows that we may set

$$ w_2' = \hat{w}_2 e^{ikx-k(z-\bar{h})+\sigma t}, \qquad (3.1a) $$

$$ p_2' = \hat{p}_2 e^{ikx-k(z-\bar{h})+\sigma t}, \qquad (3.1b) $$

where $\hat{\psi}_2$, $\hat{w}_2$, and $\hat{p}_2$ are the constants to be determined. The same conventions for the notation are also applied to the lower-layer solutions below.

The treatment of the lower layer is slightly more involved, because the vorticity is forced by the buoyancy. Nevertheless, by taking into account of the bottom boundary condition (2.13b), we may set

$$ \zeta_1' = \hat{\zeta}_1 \sin mz e^{ikx+\sigma t}, \qquad (3.2a) $$

$$ w_1' = \hat{w}_1 \sin mz e^{ikx+\sigma t}, \qquad (3.2b) $$

$$ p_1' = \hat{p}_1 \cos mz e^{ikx+\sigma t}, \qquad (3.2c) $$

$$ b_1' = \hat{b}_1 \sin mz e^{ikx+\sigma t}. \qquad (3.2d) $$

Here, in the lower layer, the horizontal and the vertical scales are characterized by different wavenumbers $k$ and $m$. Note that at this stage, a possibility that the vertical wavenumber $m$ is purely imaginary as in the upper layer is not excluded, but it is only excluded a posteriori.

From Eq. (2.16a), we find

$$ \hat{\zeta}_1 = - \frac{ik\hat{b}_1}{\sigma + ikU_1}. $$

It immediately follows from Eq. (2.15d) that

$$ \hat{w}_1 = \frac{k^2}{(k^2 + m^2)(\sigma + ikU_1)} \hat{b}_1 \qquad (3.3a) $$

or

$$ \hat{b}_1 = \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \hat{w}_1. \qquad (3.3b) $$

Note that Eq. (3.3a) corresponds to Eq. (2.53) of Fiedler (1984). Substitution of Eq. (3.3b) into Eq. (2.14a) further finds

$$ \hat{p}_1 = - \frac{\rho_1 m}{k^2} (\sigma + ikU_1) \hat{w}_1. \qquad (3.4a) $$

A similar procedure applied to the upper layer leads to

$$ \hat{p}_2 = \frac{\rho_2}{k} (\sigma + ikU_2 + k\bar{h}) \hat{w}_2. \qquad (3.4b) $$

Application of the height perturbation equation, Eq. (2.13d), to both layers leads to

$$ \hat{w}_1 = \frac{\sigma + ikU_1 + D}{\sin m\bar{h}} \hat{\eta}, \qquad (3.5a) $$

$$ \hat{w}_2 = (\sigma + ikU_2 + D)\hat{\eta}, \qquad (3.5b) $$

and further substitution of Eqs. (3.5a) and (3.5b), respectively, into Eqs. (3.4a) and (3.4b) results in

$$ \hat{p}_1 = - \frac{\rho_1 m}{k^2} (\sigma + ikU_1)(\sigma + ikU_1 + D) \frac{\hat{\eta}}{\sin m\bar{h}}, \qquad (3.6a) $$

$$ \hat{p}_2 = \frac{\rho_2}{k} (\sigma + ikU_2 + D)(\sigma + ikU_2 + k\bar{h})\hat{\eta}. \qquad (3.6b) $$

Finally, substitution of Eqs. (3.6a) and (3.6b) into the pressure boundary condition (2.13c) leads to an eigenvalue problem to be solved:

$$ -\rho_1 \frac{m}{k^2} (\sigma + ikU_1)(\sigma + ikU_1 + D) \cot m\bar{h} - \frac{\rho_2}{k} (\sigma + ikU_2 + D) \\ \times (\sigma + ikU_2 + k\bar{h}) - (\rho_1 - \rho_2)g = 0. \qquad (3.7) $$

#### b. Buoyancy problem

Another eigenvalue problem is obtained from the buoyancy equation, Eq. (2.9a). By substitution of the general solutions, we obtain

$$ [\bar{h}(\sigma + ikU_1) + w_e] (\sin mz) \hat{b}_1 = \alpha \hat{\eta}. $$
Here, the vertical average, $\langle \sin mz \rangle$, is evaluated by

$$ \langle \sin mz \rangle = \frac{1}{\bar{h}} \int_{0}^{\bar{h}} \sin mz \, dz = - \frac{1}{m\bar{h}} \cos mz \Big|_{0}^{\bar{h}} = \frac{1 - \cos m\bar{h}}{m\bar{h}}. $$

Thus,

$$ \hat{\eta} = \frac{1}{\alpha m \bar{h}} [\bar{h}(\sigma + ikU_1) + w_e](1 - \cos m\bar{h})\hat{b}_1. \quad (3.8a) $$

On the other hand, by combining Eqs. (3.3b) and (3.5a), we obtain

$$ \hat{b}_1 = \frac{(k^2 + m^2)(\sigma + ikU_1)(\sigma + ikU_1 + D)}{k^2 \sin m\bar{h}} \hat{\eta}. \quad (3.8b) $$

By substituting Eq. (3.8b) into Eq. (3.8a), we obtain the second eigenvalue problem

$$ (k^2 + m^2)(\sigma + ikU_1)(\sigma + ikU_1 + D)[\bar{h}(\sigma + ikU_1) + w_e] \times (1 - \cos m\bar{h}) - \alpha m k^2 \bar{h} \sin m\bar{h} = 0. \quad (3.9) $$

As it turns out from the result of section 4, a main balance in Eq. (3.8b) that controls the system is

$$ (k^2 + m^2)\hat{\eta} \sim \hat{b}_1; \quad (3.10) $$

thus, the interface is displaced by the buoyancy more efficiently for larger horizontal scales (i.e., the smaller $k^2$). A larger interface displacement $\hat{\eta}$, in turn, leads to a further buoyancy perturbation through Eq. (2.9a); thus, the system becomes more unstable for the larger scales as will be found in section 4.

### c. Eigenvalue problems

As the analysis of the last two subsections show, the stability problem reduces to that of solving the two eigenvalue problems given by Eqs. (3.7) and (3.9). Here, the problem consists of defining two eigenvalues: the growth rate $\sigma$ and the vertical wavenumber $m$ of the mixed layer for a given horizontal wavenumber $k$. Thus, two eigen equations must be solved for these two eigenvalues.

In the following, we first nondimensionalize these two eigen equations, then after general discussions, derive a general solution for the growth rate obtained from a nondimensionalized version of Eq. (3.7). This solution has a general validity. It also constitutes a self-contained solution when a coupling of the dynamical system considered in sections 2c and 3a with the buoyancy is turned off by setting $\alpha = 0$ in Eq. (2.9a).

We note in Eq. (3.7) that a key free parameter of the problem is

$$ \mu = \frac{m}{k} \cot m\bar{h}. \quad (3.11a) $$

A key parameter in Eq. (3.9) is $\alpha$, which is nondimensionalized into

$$ \tilde{\alpha} = (kg^3)^{-1/2} \alpha. \quad (3.11b) $$

Nondimensional versions of Eqs. (3.7) and (3.9) are given by

$$ \mu(\tilde{\sigma} + i\tilde{U}_1)(\tilde{\sigma} + i\tilde{U}_1 + \tilde{D}) + \tilde{\rho}(\tilde{\sigma} + i\tilde{U}_2 + \tilde{D})(\tilde{\sigma} + i\tilde{U}_2 + \tilde{h}) + (1 - \tilde{\rho}) = 0, \quad (3.12a) $$

$$ (1 + \tilde{m}^2)(\tilde{\sigma} + i\tilde{U}_1)(\tilde{\sigma} + i\tilde{U}_1 + \tilde{D})[\tilde{h}(\tilde{\sigma} + i\tilde{U}_1) + \tilde{w}_e] \times (1 - \cos \tilde{m}\tilde{h}) - \tilde{\alpha}\tilde{m}\tilde{h} \sin \tilde{m}\tilde{h} = 0, \quad (3.12b) $$

where the nondimensional parameters and variables are introduced by

$$ \tilde{\sigma} = (kg)^{-1/2} \sigma, \quad (3.13a) $$

$$ \tilde{U}_j = (k/g)^{1/2} U_j, \quad (3.13b) $$

$$ \tilde{D} = (kg)^{-1/2} D, \quad (3.13c) $$

$$ \tilde{\rho} = \rho_2 / \rho_1, \quad (3.13d) $$

$$ \dot{\tilde{h}} = (k/g)^{1/2} \tilde{h}, \quad (3.13e) $$

$$ \tilde{w}_e = (k/g)^{1/2} w_e, \quad (3.13f) $$

$$ \tilde{m} = m/k, \quad (3.13g) $$

$$ \tilde{h} = k\tilde{h}, \quad (3.13h) $$

for $j = 1, 2$. Note that a tilde $\tilde{\cdot}$ is added for designating the nondimensional variables.

A convenient general strategy for solving this set of eigen equations would be to first solve Eq. (3.12a) for $\tilde{\sigma}$, and by substituting this result, solve Eq. (3.12b) for $\tilde{m}$. Note that Eq. (3.12a) is only the second order in respect to $\tilde{\sigma}$; thus, an analytical solution for the latter is readily obtained. On the other hand, the resulting equation by substituting this result into Eq. (3.12b) is transcendental in respect to $\tilde{m}$. Thus, the solution for $\tilde{m}$ must be sought numerically in general cases.

The general solution for the growth rate $\tilde{\sigma}$ obtained from Eq. (3.12a) is

$$ \tilde{\sigma} = -i\tilde{U}_1 \frac{\mu + \tilde{\rho}\hat{U}}{\mu + \tilde{\rho}} - \frac{(\mu + \tilde{\rho})\hat{D} + \tilde{\rho}\hat{h}}{2(\mu + \tilde{\rho})} \tilde{U}_1 \pm \frac{(\mu\tilde{\rho})^{1/2}\tilde{U}_1}{\mu + \tilde{\rho}} \left[ (1 - \hat{U})^2 (1 - \tilde{R}\tilde{i}) + \frac{\tilde{\rho}}{4\mu} \left( \hat{h} - \frac{\mu + \tilde{\rho}}{\tilde{\rho}} \hat{D} \right)^2 + i(1 - \hat{U})\hat{h} \right]^{1/2}. \quad (3.14) $$

Here, for simplifying the final expression, some nondimensional parameters have been normalized by $\tilde{U}_1$:

$$ \hat{U} = \tilde{U}_2 / \tilde{U}_1, \quad (3.15a) $$

$$ \dot{\hat{h}} = \dot{\tilde{h}} / \tilde{U}_1, \quad (3.15b) $$

$$ \hat{D} = \tilde{D} / \tilde{U}_1. \quad (3.15c) $$
Furthermore, a Richardson number $\tilde{Ri}$, is introduced by

$$ \tilde{Ri} = \frac{(\mu + \bar{\rho})(1 - \bar{\rho})}{\mu \bar{\rho} \tilde{U}_1^2 (1 - \tilde{U})^2} = \left(\frac{g}{k}\right) \frac{(\mu \rho_1 + \rho_2)(\rho_1 - \rho_2)}{\mu \rho_1 \rho_2 (U_1 - U_2)^2}. \quad (3.16) $$

Note especially that the system is unstable when $\tilde{Ri} < 1$ and the shear is strong enough. However, both the deepening of the mixed layer, $\tilde{h} > 0$, and the divergence, $\tilde{D} > 0$, tend to suppress the destabilization tendency.

## 4. Simple solutions

### a. Simplest case

The general solution (3.14) is clearly a rich source of instabilities, including a contribution of the shear with $\tilde{Ri}$, that is clearly worthwhile for further investigations (cf. Brost et al. 1982; Kurowski et al. 2009; Mellado et al. 2009; Katzwinkel et al. 2012; Malinowski et al. 2013). However, for focusing on the entrainment-driven instability, we turn off here the background winds $\tilde{U}_1 = \tilde{U}_2 = 0$. In this subsection, we consider the simplest case by further setting $\tilde{h} = \tilde{D} = 0$. As a result, the growth rate obtained from Eq. (3.12a) reduces to

$$ \tilde{\sigma}^2 = -\frac{1 - \bar{\rho}}{\mu + \bar{\rho}}. \quad (4.1a) $$

It suggests that when the system is unstable [i.e., $R(\tilde{\sigma}) > 0$], the mode is purely growing with no imaginary component. These simplifications also make the structure of the solution much simpler: we find immediately from Eq. (3.3a) that the mixed-layer vertical velocity $w'_1$ is in phase with the buoyancy perturbation $b'_1$ with the same sign, i.e., $w'_1 \sim b'_1$. Same wise, we find $w'_1 \sim w'_2 \sim \eta$ from Eqs. (3.5a) and (3.5b), and $-p'_1 \sim p'_2 \sim \eta$ from Eqs. (3.6a) and (3.6b).

The remainder of this subsection provides a self-contained mathematical description of how a closed analytic solution is derived. Readers who wish only to see the final results may proceed directly to the last two paragraphs of this subsection.

Equation (3.12b) reduces to

$$ (1 + \tilde{m}^2)\tilde{\sigma}^2(\tilde{h}\tilde{\sigma} + \tilde{w}_e)(1 - \cos \tilde{m}\tilde{h}) - \tilde{\alpha}\tilde{m}\tilde{h}\sin \tilde{m}\tilde{h} = 0. \quad (4.1b) $$

We immediately notice that by substituting an explicit expression (4.1a) for $\tilde{\sigma}^2$ into Eq. (4.1b), the latter further reduces to

$$ -(1 + \tilde{m}^2) \frac{1 - \bar{\rho}}{\mu + \bar{\rho}} (\tilde{h}\tilde{\sigma} + \tilde{w}_e)(1 - \cos \tilde{m}\tilde{h}) - \tilde{\alpha}\tilde{m}\tilde{h}\sin \tilde{m}\tilde{h} = 0. \quad (4.1c) $$

Here, a term with $\tilde{\sigma}$ is left unsubstituted for an ease of obtaining a final result later.

When the dynamics is not coupled with the buoyancy anomaly with $\tilde{\alpha} = 0$, there are three possible manners for satisfying Eq. (4.1c): setting $\tilde{m}^2 = -1$, $\tilde{\sigma} = -\tilde{w}_e/\tilde{h}$, or $\cos \tilde{m}\tilde{h} = 1$. The first possibility leads to

$$ \mu = \coth \tilde{h}. $$

In this case, $\mu$ is always positive so long as $\tilde{h} > 0$. Thus, the system is always stable so long as it is stably stratified with $\tilde{\rho} < 1$ according to Eq. (4.1a). The second gives a damping mode with the value of $\mu$ to be defined from Eq. (4.1a) by substituting this expression for $\tilde{\sigma}$. The last possibility leads to $\mu \to +\infty$; thus, the system becomes neutrally stable.

On the other hand, when the dynamics is coupled with the buoyancy anomaly with $\tilde{\alpha} \neq 0$, the parameter $\mu$ may turn negative; thus, the solution (4.1a) may become unstable. Here, recall the definition (3.11a) of this parameter, in which $\cot \tilde{m}\tilde{h}$ is a monotonously decreasing function of $\tilde{m}\tilde{h}$, and it changes from $+\infty$ to $-\infty$ as $\tilde{m}\tilde{h}$ changes from 0 to $\pi$, passing $\cot \tilde{m}\tilde{h} = 0$ at $\tilde{m}\tilde{h} = \pi/2$. For focusing on the state with $\cot \tilde{m}\tilde{h}$ negative enough, we take the limit toward $\tilde{m}\tilde{h} \to \pi$, and set

$$ \tilde{m}\tilde{h} = \pi - \Delta\tilde{m}\tilde{h}. \quad (4.2) $$

We expect that $(0 <) \Delta\tilde{m}\tilde{h} \ll 1$.

Note that $\tilde{m}\tilde{h} = \pi$ corresponds to a solution that the perturbation vertical velocity vanishes exactly at the inversion height, $z = \tilde{h}$, and as a result, the disturbance is strictly confined to the mixed layer without disturbing the inversion interface. In this case, no buoyancy anomaly is induced. Equation (4.2) with $\tilde{m}\tilde{h} < \pi$ suggests that the perturbation vertical velocity slightly intrudes into the free atmosphere.

Under the approximation (4.2), we obtain

$$ \sin \tilde{m}\tilde{h} \approx \Delta\tilde{m}\tilde{h}, \quad (4.3a) $$

$$ \cos \tilde{m}\tilde{h} \approx -1, \quad (4.3b) $$

as well as

$$ \mu \approx -\tilde{m}(\Delta\tilde{m}\tilde{h})^{-1}, \quad (4.4) $$

where

$$ \tilde{m} \approx \pi/\tilde{h} = \pi/k\tilde{h} \quad (4.5) $$

from the leading-order expression in Eq. (4.2). Note that from Eq. (4.4) and an assumption of $|\Delta\tilde{m}\tilde{h}| \ll 1$, we also expect $|\mu| \gg 1$. As a result, in the growth rate (4.1a), $\mu$ becomes dominant in denominator, and it reduces to

$$ \tilde{\sigma}^2 \approx -\frac{1 - \bar{\rho}}{\mu} \approx \frac{1 - \bar{\rho}}{\tilde{m}} \Delta\tilde{m}\tilde{h}. \quad (4.6) $$

By substituting all the approximations introduced so far into Eq. (4.1c):

$$ 2(1 + \tilde{m}^2) \frac{1 - \bar{\rho}}{\tilde{m}} (\Delta\tilde{m}\tilde{h})(\tilde{h}\tilde{\sigma} + \tilde{w}_e) - \tilde{\alpha}\tilde{m}\tilde{h}\Delta\tilde{m}\tilde{h} \approx 0. $$

Two major terms share a common factor $\Delta\tilde{m}\tilde{h}$ that can simply be dropped off, and a slight rearrangement gives

$$ \tilde{\sigma} + \frac{\tilde{w}_e}{\tilde{h}} \approx \frac{\tilde{\alpha}}{2(1 - \bar{\rho})} \frac{\tilde{m}^2}{1 + \tilde{m}^2}. $$

It leads to a final expression:

$$ \tilde{\sigma} = -\mathcal{D} + \mathcal{A}, \quad (4.7) $$
where

$$ \mathcal{D} = \frac{\tilde{w}_e}{\bar{h}} = k^{-1/2} \mathcal{D}_0, \qquad (4.8a) $$

$$ \mathcal{A} = \frac{\tilde{\alpha}}{2(1 - \bar{\rho})} \left( \frac{\tilde{m}^2}{1 + \tilde{m}^2} \right) = k^{-1/2} \tilde{\omega}(k) \mathcal{A}_0, \qquad (4.8b) $$

with the coefficients, $\mathcal{D}_0$ and $\mathcal{A}_0$, and a function $\tilde{\omega}(k)$ defined by

$$ \mathcal{D}_0 = \frac{w_e}{g^{1/2} \bar{h}} \sim 10^{-4} \text{ km}^{-1/2}, \qquad (4.9a) $$

$$ \mathcal{A}_0 = \frac{\alpha}{2(1 - \bar{\rho})g^{3/2}} \sim 10^{-4} \text{ km}^{-1/2}, \qquad (4.9b) $$

$$ \tilde{\omega}(k) = [1 + (k\bar{h}/\pi)^2]^{-1}. \qquad (4.9c) $$

Here, the order of magnitude estimates above are based on the values listed in appendix B. By further substituting the expressions (4.8a) and (4.8b) into Eq. (4.7):

$$ \tilde{\sigma} \approx [-\mathcal{D}_0 + \tilde{\omega}(k)\mathcal{A}_0]k^{-1/2}. \qquad (4.10) $$

Finally, the growth rate of the instability is given by

$$ \sigma = g^{1/2}[-\mathcal{D}_0 + \tilde{\omega}(k)\mathcal{A}_0] \qquad (4.11) $$

after dimensionalizing the result (4.10) by following Eq. (3.13a). Here, $\tilde{\omega}(k)$ is a decreasing function of $k$, and asymptotically $\tilde{\omega}(k) \to 1$ and $0$, respectively, toward $k \to 0$ and $+\infty$. Thus, the growth rate is asymptotically $\sigma \to g^{1/2}(-\mathcal{D}_0 + \mathcal{A}_0)$ and $\sigma \to -g^{1/2}\mathcal{D}_0$, respectively, as $k \to 0$ and $+\infty$. It is seen that the sign of the growth rate with $k \to 0$ is defined by relative magnitudes of the entrainment-induced mechanical damping $\mathcal{D}_0$ and buoyancy feedback $\mathcal{A}_0$. When the latter dominates the system is unstable in the large-scale limit, whereas when the former dominates it is damping. As the horizontal scale decreases (toward $k \to +\infty$), contribution of the entrainment-induced buoyancy feedback gradually decreases, and the system becomes simply stable due to the mechanical damping. These points are visually demonstrated in Fig. 2 by plotting the growth rates for selected values of $\mathcal{A}_0/\mathcal{D}_0$. Here, the order of magnitude of the growth rate is estimated as $\sigma \sim g^{1/2}\mathcal{D}_0 \sim g^{1/2}\mathcal{A}_0 \sim 10^{-5} \text{ s}^{-1}$.

Recall that this solution is derived under an approximation of Eq. (4.2). Under this approximation, we seek a solution with convective plumes in the mixed layer slightly intruding into the free troposphere (cf. Fig. 3), as inferred by examining the assumed solution forms (3.2a)–(3.2d). By combining this fact with the phase relations between the variables already identified Eqs. (3.3a), (3.3b), (3.4a), (3.4b), (3.5a), (3.5b), (3.6a), and (3.6b), we can easily add spatial distributions of the other variables to Fig. 3, as already outlined after Eq. (4.1a) in section 4a.

b. Large-scale divergence effect

The simplest case considered in the last subsection illustrates how an instability can be induced by the mixed-layer-

<figure><img src="image_4.png" /><figcaption>FIG. 2. Nondimensional growth rate $\sigma/g^{1/2}\mathcal{D}_0$ [Eq. (4.11)] as a function of the horizontal wavenumber $k$ (km$^{-1}$). The curves are with the fractional contribution of the cloud-top buoyancy feedback of $\mathcal{A}_0/\mathcal{D}_0 = 0.5$ (solid), $\mathcal{A}_0/\mathcal{D}_0 = 1$ (long dashed), and $\mathcal{A}_0/\mathcal{D}_0 = 2$ (short dashed). Note that the dimensional order of the growth rate is $g^{1/2}\mathcal{D}_0 \sim 1 \text{ day}^{-1}$.</figcaption></figure>

top entrainment. However, the setting is rather unrealistic by neglecting a contribution of the large-scale divergence rate $\tilde{D}$ to the problem. An existence of a positive finite divergence rate $\tilde{D}$ defines the equilibrium height $\bar{h}$ of the inversion under

<figure><img src="image_5.png" /><figcaption>FIG. 3. Schematic structure of the perturbation solution: the streamfunction $\psi$ (contours) and the inversion-height deformation (thick solid curve).</figcaption></figure>
its balance with the entrainment is a crucial part of the well-mixed boundary layer problem. Thus, in this subsection, we consider the modification of the problem by including a contribution of nonvanishing $\tilde{D}$.

Equation (3.12a) for the growth rate is modified to

$$ \tilde{\sigma}(\tilde{\sigma} + \tilde{D}) = -\frac{1 - \tilde{\rho}}{\mu + \tilde{\rho}}, \quad (4.12a) $$

and its solution is

$$ \tilde{\sigma} = -\frac{\tilde{D}}{2} \pm \left[ \left( \frac{\tilde{D}}{2} \right)^2 - \frac{1 - \tilde{\rho}}{\mu + \tilde{\rho}} \right]^{1/2} . \quad (4.12b) $$

Note that as suggested by the first term of the growth-rate expression, (4.12b), a primarily role of the environmental descent is to damp the inversion-interface instability. However, as seen below, the full role of the environmental descent is subtler than just seen here.

The second eigenvalue equation, Eq. (3.12b), reduces to

$$ (1 + \tilde{m}^2)\tilde{\sigma}(\tilde{\sigma} + \tilde{D})(\tilde{h}\tilde{\sigma} + \tilde{w}_e)(1 - \cos \tilde{m}\tilde{h}) - \tilde{\alpha}\tilde{m}\tilde{h}\sin \tilde{m}\tilde{h} = 0. \quad (4.12c) $$

Note that the first two appearance of $\tilde{\sigma}$ in Eq. (4.12c) exactly constitutes the expression of the left-hand side of Eq. (4.12a). The direct substitution of this expression into the corresponding place in Eq. (4.12c) reduces the latter into Eq. (4.1c), as obtained for the case without the background divergence $\tilde{D}$. In other words, the effect of the environmental descent cancel out under the inversion-interface buoyancy condition. It immediately follows that we obtain the identical growth rate as the case without background divergence.

### c. Under steady deepening by entrainment

Alternative consistent treatment is to turn off the environmental descent, i.e., $\tilde{D} = 0$, but instead, to assume that the well-mixed layer deepens steadily by entrainment; thus, $\tilde{h} \neq 0$ (and we will set $\tilde{h} = \tilde{w}_e$ at the last stage). In this case, Eq. (3.12b) still reduces to Eq. (4.1b) as in section 4a. On the other hand, Eq. (3.12a) leads to

$$ \tilde{\sigma}^2 = -\frac{1}{\mu + \tilde{\rho}} \left[ \tilde{\rho} \dot{\tilde{h}} \tilde{\sigma} + (1 - \tilde{\rho}) \right]. \quad (4.13) $$

Substituting this expression for $\tilde{\sigma}^2$ into Eq. (4.1b), and only where $\tilde{\sigma}^2$ itself is found, leads to

$$ -\frac{\tilde{\rho}\tilde{h}\dot{\tilde{h}}}{\mu + \tilde{\rho}}(1 + \tilde{m}^2)\left[\tilde{\sigma}^2 + \left(\frac{1 - \tilde{\rho}}{\tilde{\rho}\tilde{h}} + \frac{\tilde{w}_e}{\tilde{h}}\right)\tilde{\sigma} + \frac{1 - \tilde{\rho}}{\tilde{\rho}\tilde{h}\tilde{h}}\tilde{w}_e\right] \\ \times (1 - \cos \tilde{m}\tilde{h}) - \tilde{\alpha}\tilde{m}\tilde{h}\sin \tilde{m}\tilde{h} = 0. $$

Finally, as before, we introduce approximations (4.3a), (4.3b), and (4.4) obtained under $\Delta\tilde{m}\tilde{h} \ll 1$. We retain only the terms with $O(\Delta\tilde{m}\tilde{h})$. Thus, the term with $\tilde{\sigma}^2$ drops off in the above,

because it is expected to be $O(\Delta\tilde{m}\tilde{h})$ by itself. After further reductions, we obtain

$$ \tilde{\sigma} = \left( 1 + \frac{\tilde{\rho}}{1 - \tilde{\rho}} \frac{\tilde{w}_e \dot{\tilde{h}}}{\tilde{h}} \right)^{-1} (\tilde{A} - \tilde{D}). \quad (4.14) $$

The result is the same as before apart from a prefactor containing $\tilde{h} \neq 0$ to the front. The growth rate diminishes by this prefactor. The order of this correction is

$$ \frac{\tilde{\rho}}{1 - \tilde{\rho}} \frac{\tilde{w}_e \dot{\tilde{h}}}{\tilde{h}} = \frac{\tilde{\rho}}{1 - \tilde{\rho}} \frac{w_e^2}{g\tilde{h}} \sim 10^{-6}; $$

Thus, the contribution of the prefactor is negligible, and the same conclusion as before holds.

## 5. Discussion

A hydrodynamic stability analysis has been applied to the dry convective well-mixed boundary layer with an ultimate application of the methodology to the CTEI in mind. The key difference of this approach from the more conventional parcell-based analysis is that it can derive the growth rate of an instability as a function of the horizontal wavenumber as well as its spatial structure.

The analysis has identified a new type of instability associated with the mixed-layer-top entrainment. This instability arises under a competition between the destabilization tendency due to the entrainment-induced buoyancy perturbation and the stabilization tendency due to the mechanical damping associated with entrainment. Importantly, these two entrainment effects can be separated into these the two terms in the buoyancy equation, Eq. (2.9a). Damping tendency of the buoyancy perturbation is directly proportional to the entrainment velocity $w_e$, whereas the destabilization tendency by buoyancy perturbation is generated by the displacement of the interface, which is more directly controlled by another parameter $\alpha$ [Eq. (2.9b)]. This instability is, more fundamentally, driven by the entrainment, which is in turn, driven by convection in the well-mixed layer.

Obtained growth-rate tendencies with changing horizontal scales are consistent with qualitative arguments in section 3 associated with Eq. (3.10). In the small-scale limit, the mechanical damping effect dominates over the destabilization tendency by buoyancy perturbation, and as a result, the perturbation is always damping. In the large-scale limit, instability may arise when the magnitude of the destabilization tendency by buoyancy feedback is stronger than the mechanical damping as measured by a ratio between the two parameters $A_0$ and $D_0$ defined by Eqs. (4.9a) and (4.9b). A transition from the small-scale damping regime to the large-scale unstable regime is defined by the scale $k\tilde{h}/\pi \sim 1$, where the horizontal scale, $\pi/k$, of the disturbance is comparable to the mixed-layer depth $\tilde{h}$ (~1 km) with an exact transition scale depending on the ratio $A_0/D_0$. It can easily be shown that this ratio is essentially proportional to the vertical gradient of the buoyancy in the free troposphere, and a contribution of the entrainment rate is completely removed when a radiative feedback is set $Q_R = 0$ in Eq. (2.9b).
Thus, rather ironically, the identified entrainment-induced instability does not strongly depend on the entrainment rate.

The identified instability is inherently of a large-scale nature, and a reasonably large domain is required to numerically realize it, as suggested by Fig. 2. If this instability had any implication for the CTEI, it could explain why the evidence for the CTEI by LES studies so far is rather inconclusive (e.g., Kuo and Schubert 1988; Siems et al. 1990; MacVean 1993; Yamaguchi and Randall 2008). In these simulations, relatively small domain sizes (5 km square or less) are taken, that may prevent us from observing a full growth of the CTEI. The obtained growth time scale is also very slow, about an order of a day. With typically short simulation times with LESs (about few hours), that could be another reason for a difficulty for realizing a CTEI with these simulations. Direct numerical simulations (DNSs) by Mellado (2010), in spite of an advantage of resolving everything explicitly, are even in less favorable position for simulating a full CTEI due to an even smaller modeling domain. Unfortunately, dismissal of a possibility of CTEI by Mellado (2017) in his review is mostly based on this DNS result.

In contrast to these more recent studies, it may be worthwhile to note that an earlier study by Moeng and Arakawa (1980) identifies a reasonably clear evidence for CTEI over a high sea surface temperature (SST) region of their two-dimensional non-hydrostatic experiment with a 1000-km horizontal domain, assuming a linear SST distribution. A preferred scale identified by their experiment is 30–50 km, qualitatively consistent with the present linear stability analysis, although it is also close to the minimum resolved scale in their experiment due to a crude resolution. A time scale estimated from the present study is also consistent with a finding by Moeng and Arakawa (1980) that their CTEI-like structure develops taking over 24 h. However, due to limitations of their simulations with parameterizations of eddy effects, a full LES is still required to verify their result. From an observational point of view, an assumption of horizontal homogeneity of the stratocumulus over such a great distance may simply be considered unrealistic in respect of extensive spatial inhomogeneity associated with the stratocumulus as realized in LESs (e.g., Chung et al. 2012; Zhou and Bretherton 2019).

In this respect, it may be interesting to note that a recent observational study by Zhou et al. (2015) suggests a possibility of a certain cloud-top instability, if not CTEI, leading to a decoupling, which ultimately induces a transition to trade cumulus regime. We should realize that a slow time scale suggested for CTEI by the present study may be another reason for difficulties of identifying it observationally. Previous observational diagnoses on CTEI criterions have been based on instantaneous comparisons (e.g., Albrecht et al. 1985; Albrecht 1991; Kuo and Schubert 1988; Stevens et al. 2003; Mathieu and Lahellec 2005; Gerber et al. 2005, 2013, 2016). A finite time lag could be a key missing element for a successful observational identification of CTEI. If that is the case, data analyses from a point of view of the dynamical system as advocated by Yano and Plant (2012) as well as Yano et al. (2020) becomes a vital alternative approach.

On the other hand, although the present analysis has been performed by assuming a dry atmosphere, it is less likely than

in the marine stratiform-topped boundary layers that this instability is to be seen in dry well-mixed convective boundary layers. The latter typically go through very pronounced diurnal cycles with the boundary layer itself becomes stably stratified during nights; thus, a good stationarity of the system required to observe such a slow growth of instability is hardly satisfied.

The present study focuses on the instability induced by the boundary layer-top entrainment. Nevertheless, a basic formulation is presented in fully general manner. Thus, its simple extension can consider rich possibilities of the mixed-layer inversion-interface instabilities under a coupling with the buoyancy anomaly. Especially, the present formulation allows us to explicitly examine a possibility of the Kelvin–Helmholtz instability over the mixed-layer observationally suggested by Brost et al. (1982), Kurowski et al. (2009), Katzwinkel et al. (2012), and Malinowski et al. (2013).

Furthermore, the present analysis of the dry convective well-mixed layer constitutes a first step to fully examine the CTEI as a hydrodynamic-instability problem, most importantly, by explicitly introducing the evaporative cooling effect (cf. de Lozar and Mellado 2015). Other types of possible instabilities in the cloud-topped boundary layers, such as decoupling (Bretherton and Wyant 1997; Lewellen and Lewellen 2002), can equally be addressed by the present framework.

Extensive physics can also be incorporated. In this respect, LESs by Yamaguchi and Randall (2008) can be instructive: although their idealized version of LES leads to a positive feedback suggesting CTEI, the tendency is overcompensated by longwave radiation and surface heat flux in simulations with full physics. LES studies also show that the cloud-top entrainment rate is sensitively modified under aerosol–cloud (Xue et al. 2008; Hill et al. 2008, 2009) and cloud–radiation interactions (Zhou and Bretherton 2019). The present formulation provides a basis for elucidating those various feedbacks between physics under the framework of the linear-stability analysis.

A crucial aspect of the present formulation is to treat a deformation process of the inversion interface explicitly, that could ultimately transform the well-mixed layer into a cumulus regime. The main original contribution of the present study is, under a crude representation of the well-mixed layer, to present its linear growth rate as a function of the horizontal scale. More elaborated studies would certainly be anticipated, and the present study suggests that they are actually feasible. A more elaborated entrainment formulation (cf. Stevens 2002) is just one example. The most challenging step is to proceed to a fully nonlinear formulation, probably, by taking an analogy with the contour dynamics for the vortex dynamics (cf. Dritschel 1989; Dritschel and Ambaum 1997), but by considering a full nonlinear evolution of the inversion height as a contour. Such an extension would be able to simulate a transformation of stratocumulus into trade cumulus in terms of a finite amplitude deformation of the inversion height. Both modeling and observational studies are further expected to follow.

*Acknowledgments.* Chris Bretherton led my attention to Fiedler (1984), Bjorn Stevens to Mellado (2017), and Szymon
Malinowski to Zhou et al. (2015). I would also like to thank the editor in charge, David Mechem, for his enduring efforts leading to publication of the present manuscript. Many reviewers have also contributed to this process, but I would like to specifically thank an anonymous reviewer who followed the series of revisions of the manuscript to the end. I specifically refer to this reviewer at several places of the text for this reason.

*Data availability statement.* The present study does not use any data either generated numerically nor by observation nor by laboratory experiments. Programs used for generating graphics are available by request.

## APPENDIX A

### Alternative Formulation with Entrainment Perturbation

An alternative formulation of the perturbation to Eq. (2.2) with $w'_e \neq 0$ is obtained by assuming that no perturbation to the buoyancy flux $\overline{w'b'}_-$ is generated at the top of the mixed layer, i.e.,

$$ \overline{w'b'}_- = 0. \qquad (A.1a) $$

This assumption can partially justified, because it is phenomenologically suggested that the mixed-layer-top flux tends to be constrained by the surface flux (cf. Betts 1973; Carson 1973; Tennekes 1973). However, it is also emphasized that the known tendency is established only in average sense, and there is no existing observation suggesting that the same also applies to the perturbations. Nevertheless, one of the anonymous reviewers strongly prefers the assumption (A.1a) over the assumption of $w'_e = 0$ in the main text. This is the main motivation of the present appendix.

Here, note that the perturbation buoyancy flux, $\overline{w'b'}_-$, is controlled by the two perturbation variables as

$$ \overline{w'b'}_- = -w_d \Delta b' - w'_e \Delta \bar{b}. \qquad (A.1b) $$

Substitution of Eq. (A.1a) into Eq. (A.1b) leads to

$$ w'_e = -\beta \Delta b', \qquad (A.2a) $$

where

$$ \beta = w_e / \Delta \bar{b} \sim 10^{-1} \, \text{s}^{-1}. \qquad (A.2b) $$

Due to the assumption of (A.1a), the perturbation buoyancy is no longer affected by entrainment, and Eq. (2.9a) is replaced by

$$ \bar{h} \left( \frac{\partial}{\partial t} + U_1 \frac{\partial}{\partial x} \right) \langle b \rangle' = \alpha \eta \qquad (A.3a) $$

with the definition (2.9b) of the parameter $\alpha$ replaced by

$$ \alpha = -Q_R. \qquad (A.3b) $$

Thus, the perturbation buoyancy is no longer damped by the entrainment, but instead the downward displacement of the

inversion warms rather than cools the mixed layer by a radiative feedback. Consequently, the entrainment-driven instability considered in the main text is no longer available in this case.

In turn, due to the modification of the formulation, the entrainment perturbation modifies the evolution of $\eta$, with Eq. (2.13d) modified into

$$ \left( \frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x} \right) \eta = -D\eta + w' \Big|_{z=\bar{h}} + w'_e. \qquad (A.4) $$

Note that, in this case, the basic mechanism of the instability is an increase of the potential energy due to the inversion-height deformation by a entrainment forcing. It is seen by multiplying $\eta$ on Eq. (A.4) so that the conservation law for the potential energy in concern is obtained:

$$ \left( \frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x} \right) \frac{\eta^2}{2} = -D\eta^2 + w'\eta \Big|_{z=\bar{h}} + w'_e \eta. $$

Thus, the condition for the instability is given by

$$ w'_e \eta > 0. $$

A more explicit instability process is seen by substituting Eq. (A.2a) into Eq. (A.4):

$$ \left( \frac{\partial}{\partial t} + U_j \frac{\partial}{\partial x} + D^* \right) \eta = w' \Big|_{z=\bar{h}} + \beta \langle b \rangle' \qquad (A.5a) $$

with the divergence, $D$, replaced by the effective value:

$$ D^* = D + \beta \frac{d\bar{b}}{dz}. \qquad (A.5b) $$

Note that with the two terms on the right-hand side of Eq. (A.5b) are of the same order $(10^{-5} \, \text{s}^{-1})$, the modified formulation finds a larger effective divergence $D^*$ than the actual. Also note an additional term proportional to $\beta$, arising from entrainment-rate feedback, found on the right-hand side of Eq. (A.5a).

Here, a positive perturbation buoyancy, $\langle b \rangle' > 0$, leads to a decrease of the inversion jump, i.e., $\Delta b' < 0$, by Eq. (2.7a); $\Delta b' < 0$, in turn, leads to an enhanced entrainment, $w'_e > 0$, by Eq. (A.2a). As a result, a positive buoyancy perturbation in the mixed layer tends to move the inversion height farther upward, and vice versa. This tendency provides a source of new instabilities.

As a result, Eqs. (3.5a) and (3.5b) are modified into

$$ (\sigma + ikU_1 + D^*)\hat{\eta} = \hat{w}_1 \sin m\bar{h} + \beta\hat{b}_1 \langle \sin mz \rangle, \qquad (A.6a) $$

$$ (\sigma + ikU_2 + D^*)\hat{\eta} = \hat{w}_2 + \beta\hat{b}_1 \langle \sin mz \rangle. \qquad (A.6b) $$

To obtain closed expressions for $\hat{w}_j$ ($j = 1, 2$) in terms of $\hat{\eta}$, we need a help of Eq. (3.3b). Its substitution into Eq. (A.6a) leads to

$$ \hat{w}_1 = \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right]^{-1} (\sigma + ikU_1 + D^*) \frac{\hat{\eta}}{\sin \hat{m}}. \qquad (A.7a) $$
It also follows that

$$ \hat{w}_2 = \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right]^{-1} \left[ (\sigma + ikU_1 + D^*) + \frac{i(k^2 + m^2)(\sigma + ikU_1)}{k} \beta \gamma (U_2 - U_1) \right] \hat{\eta}, \quad (A.7b) $$

$$ \hat{b}_1 = \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right]^{-1} \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \times (\sigma + ikU_1 + D^*) \frac{\hat{\eta}}{\sin \hat{m}}, \quad (A.7c) $$

where

$$ \gamma = \frac{\langle \sin m\bar{z} \rangle}{\sin m\bar{h}} = \frac{1 - \cos \hat{m}}{\hat{m} \sin \hat{m}}, \quad (A.8) $$

and $\hat{m} = m\bar{h}$. Note especially that Eq. (A.7c) replaces Eq. (3.8b). Equation (A.7b) may also be worthwhile to pay particular attention to, because it suggests that the shear effect here potentially leads to a further possibility of instabilities.

With the help of those modified solutions, Eqs. (3.4a) and (3.4b) are replaced by

$$ \hat{p}_1 = - \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right]^{-1} \frac{(\sigma + ikU_1)}{k^2} \times \frac{\rho_1 m \hat{\eta}}{\sin \hat{m}}, \quad (A.9a) $$

$$ \hat{p}_2 = \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right]^{-1} \frac{(\sigma + ikU_2 + k\bar{h})}{k^2} \times \left[ (\sigma + ikU_2 + D^*) + \frac{i(k^2 + m^2)(\sigma + ikU_1)}{k} \beta \gamma (U_2 - U_1) \right] \rho_2 \hat{\eta}. \quad (A.9b) $$

By substituting Eqs. (A.9a) and (A.9b) into Eq. (2.13c), we obtain the first revised dispersion relation in place of Eq. (3.7):

$$ \begin{aligned} & -\frac{\rho_1 m}{k^2} (\sigma + ikU_1)(\sigma + ikU_1 + D^*) \cot \hat{m} \\ & -\frac{\rho_2}{k} (\sigma + ikU_2 + k\bar{h}) \left[ (\sigma + ikU_2 + D^*) + \frac{i(k^2 + m^2)(\sigma + ikU_1)}{k} \beta \gamma (U_2 - U_1) \right] \\ & = \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right] (\rho_1 - \rho_2) g. \end{aligned} \quad (A.10) $$

In the buoyancy problem, we obtain from Eq. (A.3a),

$$ \hat{\eta} = \frac{1}{\alpha m} (\sigma + ikU_1)(1 - \cos \hat{m}) \hat{b}_1. \quad (A.11) $$

Consistency between Eqs. (A.7c) and (A.11) leads to the second dispersion relation, in place of Eq. (3.9):

$$ \begin{aligned} & (k^2 + m^2)(\sigma + ikU_1)^2(\sigma + ikU_1 + D^*)(1 - \cos \hat{m}) \\ & - \alpha m k^2 \sin \hat{m} \left[ 1 + \beta \gamma \frac{(k^2 + m^2)(\sigma + ikU_1)}{k^2} \right] = 0. \end{aligned} \quad (A.12) $$

Nondimensionalization of the dispersions (A.10) and (A.12) is performed in the identical manner as in the main text except for $D^*$ is nondimensionalized into $\tilde{D}$ in place of $D$. Additionally, we set

$$ \tilde{\beta} = (gk)^{1/2} \beta. \quad (A.13) $$

The resulting nondimensional dispersion relations, obtained in place of Eqs. (3.12a) and (3.12b), are

$$ \mu(\tilde{\sigma} + i\tilde{U}_1)(\tilde{\sigma} + i\tilde{U}_1 + \tilde{D}) + \tilde{\rho}(\tilde{\sigma} + i\tilde{U}_2 + \tilde{h})[(\tilde{\sigma} + i\tilde{U}_2 + \tilde{D}) \\ + i\gamma\tilde{\beta}(1 + \tilde{m}^2)(\tilde{\sigma} + i\tilde{U}_1)(\tilde{U}_2 - \tilde{U}_1)], \\ + (1 - \tilde{\rho})[1 + \gamma\tilde{\beta}(1 + \tilde{m}^2)(\tilde{\sigma} + i\tilde{U}_1) = 0, \quad (A.14a) $$

$$ \begin{aligned} & \hat{m} \gamma (1 + \tilde{m}^2)(\tilde{\sigma} + i \tilde{U}_1)^2 (\tilde{\sigma} + i \tilde{U}_1 + \tilde{D}) \\ & - \tilde{\alpha} \tilde{m} [\gamma \tilde{\beta} (1 + \tilde{m}^2)(\tilde{\sigma} + i \tilde{U}_1) + 1] = 0. \end{aligned} \quad (A.14b) $$

Note especially rich possibilities of the shear instabilities found in Eq. (A.14a), which are to be further investigated.

However, for now, we focus on the case without shear by setting, $\tilde{U}_1 = \tilde{U}_2 = 0$. We also set $\tilde{h} = 0$. Thus, Eqs. (A.14a) and (A.14b) reduce to

$$ (\mu + \tilde{\rho})\sigma^2 + [(\mu + \tilde{\rho})\tilde{D} + \gamma\tilde{\beta}(1 - \tilde{\rho})(1 + \tilde{m}^2)]\tilde{\sigma} + (1 - \tilde{\rho}) = 0, \quad (A.15a) $$

$$ \hat{m} \gamma (1 + \tilde{m}^2) (\tilde{\sigma} + \tilde{D}) \tilde{\sigma}^2 - \tilde{\alpha} \tilde{m} [\gamma \tilde{\beta} (1 + \tilde{m}^2) \tilde{\sigma} + 1] = 0. \quad (A.15b) $$

These two relations lead to two expressions for the growth rates:

$$ \tilde{\sigma}/\tilde{D} = -F/2 \pm (-\Lambda + F^2/4)^{1/2} \equiv \sigma_{1\pm}/\tilde{D}, \quad (A.16a) $$

$$ \tilde{\sigma}/\tilde{D} = -G/2F \pm (G^2 - 4FH)^{1/2}/2F \equiv \sigma_{2\pm}/\tilde{D}. \quad (A.16b) $$

Note that the second solution is obtained by substituting an expression for $\tilde{\sigma}^2$ obtained directly from Eq. (A.15a) into Eq. (A.15b). Here,

$$ F = 1 + \frac{\gamma \lambda_1 \tilde{\rho}}{\mu + \tilde{\rho}} \left( \tilde{k} + \frac{\hat{m}^2}{\tilde{k}} \right), \quad (A.17a) $$

$$ \Lambda = \frac{\lambda_2 \tilde{\rho} \tilde{k}}{\mu + \tilde{\rho}}, \quad (A.17b) $$

$$ G = 1 + \Lambda \left[ 1 + \gamma \lambda_3 \left( 1 + \frac{\hat{m}^2}{\tilde{k}^2} \right) \right] + \lambda_4, \quad (A.17c) $$

$$ H = \Lambda + \frac{\lambda_5}{\gamma} \left( 1 + \frac{\hat{m}^2}{\tilde{k}^2} \right)^{-1}, \quad (A.17d) $$
in which some parameters are introduced by

$$ \lambda_1 = \frac{\bar{\beta}_0}{\bar{D}_0} \frac{1 - \bar{\rho}}{\bar{\rho}} \sim 1, \quad \lambda_2 = \frac{1 - \bar{\rho}}{\bar{\rho} \bar{D}_0^2} \sim 10^6, \quad \lambda_3 = \bar{\beta}_0 \bar{D}_0 \sim 10^{-6}, $$

$$ \lambda_4 = \frac{\bar{\alpha}_0 \bar{\beta}_0}{\bar{D}_0^2} \sim -1, \quad \lambda_5 = \frac{\bar{\alpha}_0}{\bar{D}_0^3} \sim -10^6, $$

with the orders of magnitude estimates also provided based on the values listed in appendix B, also assuming the diverging state (i.e., $\bar{D} > 0$). As a minor modification, the sign of $\bar{\alpha}$ is simply reversed based on Eq. (1.3b). Note that the signs of $\lambda_1$, $\lambda_3$, and $\lambda_5$ must be reversed when the converging state is considered. These values are adopted in the graphic presentation below. The subscript, 0, in the above definitions suggest the values of the nondimensionalized parameters with $k = \bar{h}^{-1}$. Finally, the nondimensional horizontal wavenumber is introduced by $\bar{k} = k\bar{h}$.

Asymptotic expansions as in the main text have also been attempted, but without success. Thus, eigensolutions are numerically sought over the range of $10^{-2} \le \bar{k} \le 10^2$ and $0 < \hat{m}/\pi < 2$. More precisely, solutions are sought for a fixed horizontal wavenumber by dividing the full range of $\hat{m}$ in 400 segments of equal distance, and an interval that two curves, $\sigma_{1\pm}$ and $\sigma_{2\pm}$, cross is identified. A solution with $\sigma_{1\pm} = \sigma_{2\pm}$ is determined iteratively by a bisection method over this identified segment.

The result is summarized in Fig. A1. Here, only the most prominent mode (i.e., fastest growing or least damping) is plotted both for diverging and converging background flows. In the divergent-flow regime (i.e., $\bar{D} > 0$: solid), the nondimensional growth rate increases with the increasing wavenumber, and it reaches $10^2$ (approximately in the unit of day$^{-1}$) above $k = 1$ (km$^{-1}$). This is unrealistic, because no growth of disturbances at these scales is observed in dry convective well-mixed boundary layer. On the other hand, one of the anonymous reviewers believe that the result herein with much simplified assumptions does not necessarily mean that this instability mechanism is unrealistic: presumably, unrealistically large growth rates are attributed to the absence of horizontal smoothing due to turbulence in the present analysis.

Note that in the converging-flow regime (i.e., $\bar{D} < 0$), no instability realizes with the least-damping mode plotted instead.

# APPENDIX B

## Typical Physical Values

Typical physical values (in the orders of magnitudes) of the problem are as follows:

* Acceleration of the gravity: $g \sim 10 \text{ m s}^{-2}$

* Entrainment rate: $w_e \sim 10^{-2} \text{ m s}^{-1}$ (cf. Stevens et al. 2003; Gerber et al. 2013)

* Inversion height: $\bar{h} \sim 10^3 \text{ m}$ (cf. Schubert et al. 1979)

Here, the values for $w_e$ and $\bar{h}$ may be considered upper bounds, but they provide convenient rounded-up values.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. A1. (a) Nondimensional growth rate $\tilde{\sigma}/|\bar{D}|$ in the same format as in Fig. 2 but assuming that the perturbation of the buoyancy flux vanishes at the top of the mixed layer: when the background flow is diverging (i.e., $\bar{D} > 0$: solid) and diverging (i.e., $\bar{D} < 0$: long dashed). (b) As in (a), but for the vertical wavelength $\pi/m$ of the perturbation mode is plotted as a function of the wavenumber $k$.</figcaption></figure>

These two values further provide an estimate of a typical divergence rate:

$$ D = w_e / \bar{h} \sim 10^{-5} \text{ s}^{-1} $$

(cf. Schubert et al. 1979).

The feedback rate $\alpha$ of the inversion height anomaly $\eta$ to the buoyancy anomaly $\langle b \rangle'$ is estimated by substituting these typical values into Eq. (2.9b) as

$$ \alpha \sim w_e \frac{d\bar{b}}{dz} \sim 10^{-2} \text{ m s}^{-1} \times 10^{-4} \text{ s}^{-2} \sim 10^{-6} \text{ m s}^{-3}, $$

where

$$ \frac{d\bar{b}}{dz} \sim \frac{g}{\bar{\theta}} \frac{d\bar{\theta}}{dz} \sim 10 \text{ m s}^{-2} \times \frac{3 \times 10^{-3} \text{ K m}^{-1}}{300 \text{ K}} \sim 10^{-4} \text{ s}^{-2}, $$

and $\bar{\theta}$ is the basic state for the potential temperature. It further provides a rate of the change of buoyancy anomaly by

$$ \alpha \frac{\eta}{\bar{h}} \sim \alpha \sim 10^{-6} \text{ m s}^{-3}, $$

which leads to a buoyancy anomaly of the order $\langle b \rangle' \sim 10^{-2} \text{ m s}^{-2}$ over a period of an hour ($\sim 10^4 \text{ s}$). This value may be considered an underestimate compared with those obtained by local
analyses: $0.01 < b < 0.2 \text{ m s}^{-2}$ (Fig. 3 of Stevens 2002), $-2 < b < 1 \text{ K}$ (Fig. 2 of Duynkerke 1993).

We also set $1 - \bar{\rho} \approx 10^{-2}$ assuming a jump of the temperature $\Delta T \approx 3 \text{ K}$ crossing the inversion in estimating the parameter values in the main text.

REFERENCES

Albrecht, B. A., 1991: Fractional cloudiness and cloud-top entrainment instability. *J. Atmos. Sci.*, **48**, 1519–1525, https://doi.org/10.1175/1520-0469(1991)048<1519:FACTE>2.0.CO;2.

—, R. S. Penc, and W. H. Schubert, 1985: An observational study of cloud-topped mixed layers. *J. Atmos. Sci.*, **42**, 800–822, https://doi.org/10.1175/1520-0469(1985)042<0800:AOSOCT>2.0.CO;2.

Betts, A. K., 1973: Nonprecipitating cumulus convection and its parameterization. *Quart. J. Roy. Meteor. Soc.*, **99**, 178–196, https://doi.org/10.1002/qj.49709941915.

—, 1974: Reply to comment on the paper 'Non-precipitating cumulus convection and its parameterization'. *Quart. J. Roy. Meteor. Soc.*, **100**, 469–471, https://doi.org/10.1002/qj.49710042517.

Bretherton, C. S., and M. C. Wyant, 1997: Moisture transport, lower-tropospheric stability, and decoupling of cloud-topped boundary layers. *J. Atmos. Sci.*, **54**, 148–167, https://doi.org/10.1175/1520-0469(1997)054<0148:MTLTSA>2.0.CO;2.

Brost, R. A., J. C. Wyngaard, and D. H. Lenschow, 1982: Marine stratocumulus layers. Part II: Turbulence budget. *J. Atmos. Sci.*, **39**, 818–836, https://doi.org/10.1175/1520-0469(1982)039<0818:MSLPIT>2.0.CO;2.

Carson, D. J., 1973: The development of a dry inversion-capped convectively unstable boundary layer. *Quart. J. Roy. Meteor. Soc.*, **99**, 450–467, https://doi.org/10.1002/qj.49709942105.

Chung, D., G. Matheou, and J. Teixeira, 2012: Steady-state large-eddy simulations to study the stratocumulus to shallow cumulus cloud transition. *J. Atmos. Sci.*, **69**, 3264–3276, https://doi.org/10.1175/JAS-D-11-0256.1.

Deardorff, J. W., 1976: On the entrainment rate of a stratocumulus-topped mixed layer. *Quart. J. Roy. Meteor. Soc.*, **102**, 563–582, https://doi.org/10.1002/qj.49710243306.

—, 1980: Cloud top entrainment instability. *J. Atmos. Sci.*, **37**, 131–147, https://doi.org/10.1175/1520-0469(1980)037<0131:CTEI>2.0.CO;2.

de Lozar, A., and J. P. Mellado, 2015: Mixing driven by radiative and evaporative cooling at the stratocumulus top. *J. Atmos. Sci.*, **72**, 4681–4700, https://doi.org/10.1175/JAS-D-15-0087.1.

Drazin, P. G., and W. H. Reid, 1981: *Hydrodynamic Stability*. Cambridge University Press, 527 pp.

Dritschel, D. G., 1989: Contour dynamics and contour surgery: Numerical algorithms for extended, high-resolution modelling of vortex dynamics in two-dimensional, inviscid, incompressible flows. *Comput. Phys. Rep.*, **10**, 77–146, https://doi.org/10.1016/0167-7977(89)90004-X.

—, and M. H. P. Ambaum, 1997: A contour-adptive semi-Lagrangian numerical algorithm for simulating fine-scale conservative dynamical fields. *Quart. J. Roy. Meteor. Soc.*, **123**, 1097–1130, https://doi.org/10.1002/qj.49712354015.

Duynkerke, P. G., 1993: The stability of cloud top with regard to entrainment: Amendment of the theory of cloud-top entrainment instability. *J. Atmos. Sci.*, **50**, 495–502, https://doi.org/10.1175/1520-0469(1993)050<0495:TSOCTW>2.0.CO;2.

Fiedler, B. H., 1984: The mesoscale stability of entrainment into cloud-topped mixed layers. *J. Atmos. Sci.*, **41**, 92–101, https://doi.org/10.1175/1520-0469(1984)041<0092:TMSOEI>2.0.CO;2.

—, 1985: Mesoscale cellular convection: Is it convection? *Tellus*, **37A**, 163–175, https://doi.org/10.3402/tellusa.v37i2.11663.

Gerber, H., S. P. Malinwoski, J.-L. Brenguier, and F. Brunet, 2005: Holes and entrainment in stratocumulus. *J. Atmos. Sci.*, **62**, 443–459, https://doi.org/10.1175/JAS-3399.1.

—, G. Frick, S. P. Malinowski, H. Jonsson, D. Khelif, and S. K. Krueger, 2013: Entrainment rates and microphysics in POST stratocumulus. *J. Geophys. Res. Atmos.*, **118**, 12094–12109, https://doi.org/10.1002/jgrd.50878.

—, S. P. Malinowski, and H. Jonsson, 2016: Evaporative and radiative cooling in POST stratocumulus. *J. Atmos. Sci.*, **73**, 3877–3884, https://doi.org/10.1175/JAS-D-16-0023.1.

Hill, A. A., S. Dobbie, and Y. Yin, 2008: The impact of aerosols on non-precipitating marine stratocumulus. I: Model description and prediction of the indirect effect. *Quart. J. Roy. Meteor. Soc.*, **134**, 1143–1154, https://doi.org/10.1002/qj.278.

—, G. Feingold, and H. Jiang, 2009: The influence of entrainment and mixing assumption on aerosol–cloud interactions in marine stratocumulus. *J. Atmos. Sci.*, **66**, 1450–1464, https://doi.org/10.1175/2008JAS2909.1.

Hoskins, B. J., and I. N. James, 2014: *Fluid Dynamics of the Mid-Latitude Atmosphere*. Wiley, 432 pp.

Katzwinkel, J., H. Siebert, and R. A. Shaw, 2012: Observation of a self-limiting, shear-induced turbulent inversion layer above marine stratocumulus. *Bound.-Layer Meteor.*, **145**, 131–143, https://doi.org/10.1007/s10546-011-9683-4.

Kuo, H.-C., and W. H. Schubert, 1988: Stability of cloud-topped boundary layers. *Quart. J. Roy. Meteor. Soc.*, **114**, 887–916, https://doi.org/10.1002/qj.49711448204.

Kurowski, M. J., S. Malinowski, and W. W. Grabowski, 2009: A numerical investigation of entrainment and transport within a stratocumulus-topped boundary layer. *Quart. J. Roy. Meteor. Soc.*, **135**, 77–92, https://doi.org/10.1002/qj.354.

Lenschow, D. H., M. Zhou, X. Zeng, L. Chen, and X. Xu, 2000: Measurements of fine-scale structure at the top of marine stratocumulus. *Bound.-Layer Meteor.*, **97**, 331–357, https://doi.org/10.1023/A:1002780019748.

Lewellen, D. C., and W. S. Lewellen, 2002: Entrainment and decoupling relations for cloudy boundary layers. *J. Atmos. Sci.*, **59**, 2966–2986, https://doi.org/10.1175/1520-0469(2002)059<2966:EADRFC>2.0.CO;2.

Lock, A. P., and M. K. MacVean, 1999: The generation of turbulence and entrainment by buoyancy reversal. *Quart. J. Roy. Meteor. Soc.*, **125**, 1017–1038, https://doi.org/10.1002/qj.4971255513.

MacVean, M. K., 1993: A numerical investigation of the criterion for cloud-top entrainment instability. *J. Atmos. Sci.*, **50**, 2481–2495, https://doi.org/10.1175/1520-0469(1993)050<2481:ANIOTC>2.0.CO;2.

—, and P. J. Mason, 1990: Cloud-top entrainment instability through small-scale mixing and its parameterization in numerical models. *J. Atmos. Sci.*, **47**, 1012–1030, https://doi.org/10.1175/1520-0469(1990)047<1012:CTEITS>2.0.CO;2.

Malinowski, and Coauthors, 2013: Physics of Stratocumulus Top (POST): Turbulent mixing across capping inversion. *Atmos. Chem. Phys.*, **13**, 12171–12186, https://doi.org/10.5194/acp-13-12171-2013.

Mathieu, A., and A. Lahellec, 2005: Comments on ‘On entrainment rates in nocturnal marine stratocumulus’ by Bjorn Stevens, Donald H. Lenschow, Ian Faloona, C.-H. Moeng,
D. K. Lilly, B. Blomquist, G. Vali, A. Bandy, T. Campos, H. Gerber, S. Haimov, B. Morley, and D. Thornton (October B, 2003, 129, 3469–3493). Quart. J. Roy. Meteor. Soc., **131**, 1293–1295, https://doi.org/10.1256/qj.04.32.

Mellado, J. P., 2010: The evaporatively driven cloud-top mixing layer. *J. Fluid Mech.*, **660**, 5–36, https://doi.org/10.1017/S0022112010002831.

—, 2017: Cloud-top entrainment in stratocumulus clouds. *Annu. Rev. Fluid Mech.*, **49**, 145–169, https://doi.org/10.1146/annurev-fluid-010816-060231.

—, B. Stevens, H. Schmidt, and N. Peters, 2009: Buoyancy reversal in cloud-top mixing layers. *Quart. J. Roy. Meteor. Soc.*, **135**, 963–978, https://doi.org/10.1002/qj.417.

Moeng, C.-H., and A. Arakawa, 1980: A numerical study of a marine subtropical stratus cloud layer and its stability. *J. Atmos. Sci.*, **37**, 2661–2676, https://doi.org/10.1175/1520-0469(1980)037<2661:ANSOAM>2.0.CO;2.

—, and U. Schumann, 1991: Composite structure of plumes in stratus-topped boundary layer. *J. Atmos. Sci.*, **48**, 2280–2291, https://doi.org/10.1175/1520-0469(1991)048<2280:CSOPIS>2.0.CO;2.

—, D. H. Lenschow, and D. A. Randall, 1995: Numerical investigations of the role of radiative and evaporative feedbacks in stratocumulus entrainment and breakup. *J. Atmos. Sci.*, **52**, 2869–2883, https://doi.org/10.1175/1520-0469(1995)052<2869:NIOTRO>2.0.CO;2.

—, B. Stevens, and P. P. Sullivan, 2005: Where is the interface of the stratocumulus-topped PBL? *J. Atmos. Sci.*, **62**, 2626–2631, https://doi.org/10.1175/JAS3470.1.

Pedlosky, J., 1987: *Geophysical Fluid Dynamics*. 2nd ed. Springer, 710 pp.

Rand, H. A., and C. S. Bretherton, 1993: Relevance of the mesoscale entrainment instability to the marine cloud-topped atmospheric boundary layer. *J. Atmos. Sci.*, **50**, 1152–1158, https://doi.org/10.1175/1520-0469(1993)050<1152:ROTMEI>2.0.CO;2.

Randall, D. A., 1980: Conditional instability of the first kind upside down. *J. Atmos. Sci.*, **37**, 125–130, https://doi.org/10.1175/1520-0469(1980)037<125:CIOTFK>2.0.CO;2.

Schubert, W. H., J. S. Wakefield, E. J. Steiner, and S. K. Cox, 1979: Marine stratocumulus convection. Part I: Governing equations and horizontally homogeneous solutions. *J. Atmos. Sci.*, **36**, 1286–1307, https://doi.org/10.1175/1520-0469(1979)036<1286:MSCPIG>2.0.CO;2.

Siems, S. T., C. S. Bretherton, M. B. Bakar, S. Shy, and R. E. Breidenthal, 1990: Buoyancy reversal and cloud-top entrainment

instability. Quart. J. Roy. Meteor. Soc., **116**, 705–739, https://doi.org/10.1002/qj.49711649309.

Stevens, B., 2002: Entrainment in stratocumulus-topped mixed layers. Quart. J. Roy. Meteor. Soc., **128**, 2663–2690, https://doi.org/10.1256/qj.01.202.

—, 2005: Atmospheric moist convection. *Annu. Rev. Earth Planet. Sci.*, **33**, 605–643, https://doi.org/10.1146/annurev.earth.33.092203.122658.

—, 2006: Bulk boundary-layer concepts for simplified models of tropical dynamics. *Theor. Comput. Fluid Dyn.*, **20**, 279–304, https://doi.org/10.1007/s00162-006-0032-z.

—, and Coauthors, 2003: On entrainment rates in nocturnal marine stratocumulus. Quart. J. Roy. Meteor. Soc., **129**, 3469–3493, https://doi.org/10.1256/qj.02.202.

Tennekes, H., 1973: A model for the dynamics of the inversion above a convective boundary layer. *J. Atmos. Sci.*, **30**, 558–567, https://doi.org/10.1175/1520-0469(1973)030<0558:AMFTDO>2.0.CO;2.

—, and A. G. M. Driedonks, 1981: Basic entrainment equations for the atmospheric boundary layer. *Bound.-Layer Meteor.*, **20**, 515–531, https://doi.org/10.1007/BF00122299.

Xue, H., G. Feingold, and B. Stevens, 2008: Aerosol effects on clouds, precipitation, and the organization of shallow cumulus convection. *J. Atmos. Sci.*, **65**, 392–406, https://doi.org/10.1175/2007JAS2428.1.

Yamaguchi, T., and D. A. Randall, 2008: Large-eddy simulation of evaporatively driven entrainment in cloud-topped mixed layers. *J. Atmos. Sci.*, **65**, 1481–1504, https://doi.org/10.1175/2007JAS2438.1.

Yano, J.-I., and R. S. Plant, 2012: Finite departure from convective quasi-equilibrium: Periodic cycle and discharge-recharge mechanism. Quart. J. Roy. Meteor. Soc., **138**, 626–637, https://doi.org/10.1002/qj.957.

—, M. H. P. Ambaum, H. Dacre, and A. Manzato, 2020: A dynamical-system description of precipitation over the tropics and the midlatitudes. *Tellus*, **72**, 1–17, https://doi.org/10.1080/16000870.2020.1847939.

Zhou, X., and C. S. Bretherton, 2019: Simulation of mesoscale cellular convection in marine stratocumulus: 2. Nondrizzling conditions. *J. Adv. Model. Earth Sys.*, **11**, 3–18, https://doi.org/10.1029/2018MS001448.

—, P. Kollias, and E. R. Lewis, 2015: Clouds, precipitation, and marine boundary layer structure during the MAGIC field campaign. *J. Climate*, **28**, 2420–2442, https://doi.org/10.1175/JCLI-D-14-00320.1.