A linear stability analysis of two-layer moist convection with a saturation interface

Abstract The linear convective instability of a mixture of dry air, water vapour and liquid water, with a stable unsaturated layer residing on an unstable saturated layer, is studied. It may serve as a prototype model for understanding the instability that causes mixing at the top of stratocumulus cloud or fog. Such a cloud-clear air interface is modelled as an infinitely thin saturation interface where radiative and evaporative cooling take place. The interface position is determined by the Clausius–Clapeyron equation, and can undulate with the evolution of moisture and temperature. In the small-amplitude regime two physical mechanisms are revealed. First, the interface undulation leads to the undulation of the cooling source, which destabilizes the system by superposing a vertical dipole heating anomaly on the convective cell. Second, the evolution of the moisture field induces non-uniform evaporation at the interface, which stabilizes the system by introducing a stronger evaporative cooling in the ascending region and vice versa in the descending region. These two mechanisms are competing, and their relative contribution to the instability is quantified by theoretically estimating their relative contribution to buoyancy flux tendency. When there is only evaporative cooling, the two mechanisms break even, and the marginal stability curve remains the same as the classic two-layer Rayleigh–Bénard convection with a fixed cooling source.

H. Fu 2. The physical model 2.1. The governing equations A two-layer set-up with idealized radiation and moist processes is built. It is modified from the fully nonlinear model of Siems et al. (1990) and Mellado et al. (2009) by adding a Dirac-delta function radiative cooling, as well as changing the basic state buoyancy profile from a local staircase reversal to a global piecewise linear one. The governing equations are expressed with explicit temperature and water species, rather than the mixing fraction formulation. The following four processes are omitted.
(i) The background vertical shear which enhances the interfacial mixing (Schulz & Mellado 2018). (ii) The droplet settling which reduces the interfacial evaporation by removing droplets near the cloud top (Bretherton, Blossey & Uchida 2007;Schulz & Mellado 2019). (iii) The lightening effect of water vapour and the loading of liquid water on buoyancy (Austin 1995). (iv) The adiabatic expansion and compression of a parcel in vertical motion which is only significant for a sufficiently deep domain, such as the whole stratocumulus layer (Siems et al. 1990).
The symbol of all the dimensional variables (but not dimensional parameters) will include a '*'. The basic state saturation interface is at z * = 0, the lower boundary is a rigid lid at z * = −H and the upper layer is infinitely deep. Figure 1 shows a schematic diagram of the set-up. The flow is assumed to be Boussinesq. The continuity equation, momentum equation, Clausius-Clapeyron equation, water vapour diagnostic equation, total water equation and thermodynamic equation are ∇ * · u * = 0, (2.1) ∂u * ∂t * + u * · ∇ * u * = − 1 ρ 0 ∇ * p * + g T 0 T * k + ν∇ * 2 u * , (2.2) q * vs (T * ) = q vs0 +λ(T * − T 0 ), where u * = u * i + v * j + w * k is the three-dimensional velocity vector, i, j and k are the three unit vectors of Cartesian coordinate, x * = x * i + y * j + z * k is the position vector, ∇ * = i∂/∂x * + j∂/∂y * + k∂/∂z * is the gradient operator, t * is time and min{} is the minimum operator. The variable p * is perturbation pressure, T * is temperature, q * vs is the saturation vapour content (the mass of saturated vapour in 1 kg of air), q * t is the total water content (the total water mass in 1 kg of air) which is the sum of vapour content q * v and liquid water content q * l . The constants ρ 0 is the reference density, ν is the kinematic viscosity, κ is the diffusivity shared by T * and all water species, L v is the evaporation latent heat and c p is the isobaric specific heat. For applications to stratocumulus, ν and κ should be viewed as eddy viscosity and diffusivity, respectively. The saturated water vapour content q * vs is linearized to T * (Randall 1980;Spyksma, Bartello & Yau 2006), withλ as the slope, T 0 and q vs0 as the reference temperature and saturated water vapour content. They are the radiative-diffusive equilibrium value at z * = 0. The constant Q * rad is the longwave radiative cooling strength (strictly speaking, the radiative flux density expressed in temperature) that depends on the Stefan-Boltzmann law. It is the small penetration depth limit of the simple radiation transfer model of de Lozar & Mellado (2013). The validity of this radiative representation is discussed in Appendix A. The vector n is the unit vector perpendicular to the saturation interface and points from the saturated region to the unsaturated region. The superscript '−' in z * − s means a tiny distance below z * s , and vice versa for '+' which will appear later.
The adjustment to thermodynamic equilibrium is assumed to be instantaneous, so the supersaturation (Chandrakar et al. 2020) is not allowed, as shown in (2.4). The assumption of identical diffusivity for T * , q * v , q * l and q * t was proposed by Bretherton (1987), and has been used in many cloud-top mixing simulations (Siems et al. 1990;Mellado et al. 2009;Mellado González 2010;de Lozar & Mellado 2015;Schulz & Mellado 2018). It leads to a drastic simplification: evaporative cooling is determined by the liquid water diffusive flux at the interface. An equivalent argument is shown in (8) of de Lozar & Mellado (2015) in their mixing fraction formulation. If the diffusivity of temperature and water species are different, diffusion will cause a phase change inside the saturated region. In the real atmosphere the dispersion of droplets is much more complicated than Fickian diffusion (Bois & Kubicki 2003).

Boundary condition
The boundary condition is identical to the simulation of Siems et al. (1990) and Mellado et al. (2009), except for an infinitely deep upper layer (rather than an upper lid). The lower boundary is a free-slip rigid lid with fixed temperature T * and total water content q * t , ∂u * ∂z * where T and q t are the temperature and water content drop across the saturated layer, respectively, and are both positive. In the real atmosphere the downward penetration depth of the mixture produced at the cloud top is limited by its moisture content. The adiabatic compression can make it dry out and then quickly gain buoyancy in the dry adiabatic process, even before reaching the cloud bottom (Siems et al. 1990). This effect, which is not explicitly considered in this model, can be implicitly carried by the lower lid. At z * → ∞, the velocity vanishes, and T * and q * t asymptotically approach a reference linear profile, where Γ T∞ (negative) and Γ q∞ = q t /H (positive) are the lapse rate of T * and q * t at z * → ∞. Strictly speaking, the upper-layer thickness should be smaller than q vs0 /Γ q∞ to make q * t positive. In the real atmosphere Γ q∞ is changing with height, and is small outside of the cloud-top mixing layer (Mellado 2017). In this model, switching the water lapse rate to a small value above the convective penetration depth H p should not influence the result. Thus, we only require q * t to be still positive at z * = H p . Considering that H p is the neutrally buoyant level of a parcel which ascends adiabatically from the lower boundary, we get H p ∼ − T/Γ T∞ . This poses a physical constraint on our model: H p q vs0 /Γ q∞ or − T/Γ T∞ q vs0 /Γ q∞ . Furthermore, as q vs0 only influences the basic state and does not enter the linear stability analysis, we can always choose a q vs0 to satisfy this physical constraint.

The non-dimensional group
The basic state is in radiative-diffusive equilibrium, without any motion. The basic state variables are denoted with an overline. The Dirac-delta function cooling at the saturation interface makes the basic state temperature profile T * piecewise linear, which decreases with height in the lower saturated layer and increases with height in the upper unsaturated layer. The basic state water content q * t decreases linearly with height. These basic state profiles will be used to non-dimensionalize the equations.
The governing equation of q * t is obtained from (2.5). Further using the boundary condition in (2.7) and (2.8), we get The governing equation of T * and its solution are obtained with (2.3), (2.4), (2.6) and (2.9), where Γ Ts = T/H (positive) is the T * lapse rate in the saturated layer. The parameters Γ Ts and Γ T∞ are linked with the radiative and evaporative cooling rate at the interface: (2.11) We choose T and q t as the temperature and water content scale, the saturated layer depth H as the length scale and H 2 /ν as the time scale. The dimensionless variables (without '*') obey The problem is controlled by five non-dimensional parameters: Rayleigh number Ra measures the relative strength of the destabilizing effect of unstable stratification and the stabilizing effect of viscosity and temperature diffusivity. Prandtl number Pr is the ratio of viscosity to temperature and water diffusivity. The Q rad measures the contribution of radiative cooling. Now, we consider λ and M that involve the moisture effect. It is worth noting that their ranges are constrained. The parameter λ represents the Clausius-Clapeyron equation which influences the interface height z s , as well as the partition between q v and q l in the saturated layer which is purely diagnostic under the thermodynamic equilibrium assumption. There must be λ < 1 to guarantee the lower layer is saturated, and λ > 0 to guarantee the monotonicity of saturated vapour pressure on T. The parameter M is the ratio of latent to dry enthalpy change across the saturated layer. It can be represented as M = (L v /c p )λλ −1 . The parameter (L v /c p )λ ranges from 1 to 2 in the 280-290 K temperature range at 900 hPa level which represents the stratocumulus (Randall 1980). This, together with 0 < λ < 1, indicates M 1.
Before linearizing the problem, the basic state non-dimensional profilesT, q t , q vs and q l are derived from their dimensional expression. The cooling strength at the interface can be further depicted with a parameter γ T which measures the temperature gradient ratio of the upper and lower layers: (2.14) Here () u and () s denote the unsaturated (upper) and saturated (lower) layer property, respectively. The parameter γ T is negative in our case where temperature decreases with height in the lower layer and increases with height in the upper layer. Using (2.3), (2.9), 928 A13-7

H. Fu
(2.10), (2.12) and (2.14), we get The parameter γ T can be linked to the radiative cooling rate Q rad and the basic state evaporative cooling rate Q evap by substituting (2.11) into (2.14), and using the non-dimensional treatment in (2.12) and (2.13): (2.16) 2.4. The linearized governing equation The disturbance non-dimensional temperature and water content is obtained by subtracting their basic state values from their total values, i.e.
(2.17a-e) The non-dimensional gradient operator is defined as ∇ = i∂/∂x + j∂/∂y + k∂/∂z. With (2.12), (2.13), (2.15) and (2.17a-e), the dimensional equations (2.1)-(2.6) are non-dimensionalized and linearized, (2.20) where H(z) is the Heaviside function. The variables q t , q vs , q l , T , u and z s all have small amplitude. The non-dimensional boundary conditions are obtained from (2.7) and (2.8), Three remarks on the governing equations are made below. First, the dδ(z)/dz term on the right-hand side of (2.23) indicates that there is a discontinuity of the small-amplitude T at z = 0, which is strictly proved later in (3.6). Figure 2. A schematic diagram of (a) the effect of the undulating interface and (b) the effect of non-uniform evaporation at the interface. The solid blue line denotes the perturbed saturation interface. The dashed black line denotes the basic state saturation interface. The vertical blue arrows denote the vertical motion at the interface. The blue patch denotes the cooling anomaly and the red patch denotes the warming anomaly. The circulating red arrows in (a,b) denote the baroclinic torque produced by the interface undulation and non-uniform evaporation, respectively. The solid red lines denote the basic state temperature profiles (T), and the dotted red lines denote the temperature profiles changed by the interface undulation or the non-uniform evaporation alone.
Physically, it is because the shift of the cooling source produces a small amplitude yet systematic deviation of the temperature profile from the basic state, as illustrated in figure 2(a). The finite-amplitude T is continuous everywhere, with a sharp transition within z ∈ (−z s , z s ). The thickness of the T transition region reduces to the thickness of the saturation interface at the small-amplitude limit, which is infinitely thin. Second, (2.21) is derived by first transforming (2.4) to q l = [1 − H(z − z s )](q t − q vs ), and then decomposing it into the basic state and perturbation variables where we have used q vs = λT in (2.15) and q vs = λT , as well as the Taylor expansion of the Heaviside function H(z − z s ) ≈ H(z) − z s δ(z). Because q t − λT = 0 at z = 0, the third term on the right-hand side of (2.25) vanishes. Because q l = [1 − H(z)](q t − λT) as is indicated by (2.15), the perturbation part of (2.25) becomes (2.21).
Third, (2.23) is derived by first linearizing ∇q l | z=z − s · n to ∂q l /∂z| z=z − s , and then using Taylor expansion of the Dirac-delta function δ(z − z s ) ≈ δ(z) − z s dδ(z)/dz to linearize the interfacial cooling term around z = 0. The derivative and Taylor expansion of a Dirac-delta function is defined by considering it as a generalized function (a distribution). See appendix C.4 of the book of Pope (2000) for a reference. The Taylor expansion of the Heaviside function also exists, because it is the integral of the Taylor expansion of a Dirac-delta function. The detail of the expansion of the Dirac-delta function term in the non-dimensional version of (2.6) is as follows: Here we have used (2.16) and (2.21). The quantities ∂q l /∂z| z=0 − , ∂q t /∂z| z=0 , ∂T /∂z| z=0 − and z s all have small amplitude, so the product between them vanish in the linear analysis. On the right-hand side of (2.26), the first term denotes the basic state part, the second term denotes the interface undulation which produces a vertical dipole heating pattern, and the third term denotes the horizontally non-uniform evaporation which produces a monopole heating pattern. The basic state part leads to a piecewise linearT and, therefore, a piecewise constant dT/dz, which corresponds to the −w term in (2.23).

Linear stability analysis
The standard practice for treating a two-layer convection problem like this is to manipulate all the equations into a single high-order equation of w for each layer respectively, and then find the proper matching condition (Ogura & Kondo 1970). It is going to be shown that the moisture and radiative effects, which are the novel contributions of this work, are interfacial processes that only require a special configuration of the matching condition. We call the case with the full non-uniform evaporation and interface undulation the 'free TRBC' (TRBC stands for two-layer Rayleigh-Bénard convection). For comparison, the classic TRBC where the cooling source is horizontally uniform and not undulating (Ogura & Kondo 1970) is called the 'fixed TRBC'. Whenever we compare these two cases, the basic state temperature profile and, therefore, γ T is identical. The physics of the fixed TRBC reported by Whitehead & Chen (1970) is briefly introduced here. The main body of the convective cell is constrained in the lower layer where there is a warm (cold) anomaly in the ascending (descending) region, so the vertical buoyancy flux is positive and kinetic energy is produced. In the upper layer the weak ascending (descending) flow generates a cold (warm) anomaly, so the vertical buoyancy flux is negative and the kinetic energy of the cell produced in the lower layer is partly diminished. The newly added radiation and evaporation should influence the instability with the buoyancy anomaly they induce. Whether such additional buoyancy anomalies favour the original cell or not should determine whether the new factors are stabilizing or destabilizing.

The eigenvalue problem
As the non-uniform evaporation and interface undulation take the form of a Dirac-delta function and its derivative in (2.23), their influence is concentrated near the saturation interface. Thus, compared with the fixed TRBC, the moisture effect only influences the matching condition. The procedure in this subsection is a standard practice that follows Ogura & Kondo (1970), except the matching conditions ofT and DT which are novel.
The normal mode solution is where k x is the x-direction wavenumber, k y is the y-direction wavenumber, σ ≡ σ r + iσ i is the complex growth factor, where σ r denotes the growth rate and σ i denotes the oscillation frequency. The z-dependent variables with a 'hat' denote the eigenfunctions. By substituting (3.1) into (2.18), (2.19) and (2.23), a sixth-order ordinary differential equation (ODE) ofŵ is obtained: The derivative of the eigenfunctions will use notation D rather than d/dz. At z = −1, the non-penetrative and free-slip velocity boundary condition, as well as the Dirichlet temperature boundary condition T | z=−1 = 0 that are adapted from (2.7) to yieldŵ As (3.2) have constant coefficients, obey the homogeneous lower boundary condition shown in (3.3) and the natural boundary condition at z → ∞, we havê Here a m , m = 1, 2, 3, 4, 5, 6 are the amplitude coefficients, p m are the eigenvalues. Substituting (3.4) into (3.2), and viewing σ , Ra and K 2 as coefficients, a complex cubic equation of p 2 m is obtained and analytically solved using the method of Lebedev (1991). As the sign of p m is not constrained by the cubic equation, we choose the p m with a positive real part.
The six coefficients a m need to be determined by six matching conditions at z = 0. The first four are the continuity of vertical velocity w, horizontal velocity u and v, the non-dimensional tangential stress components (∂u/∂z + ∂w/∂x) and (∂v/∂z + ∂w/∂y), and pressure p. Using the normal mode form of (2.18) and (2.19), a few lines of algebra yield (3.5) The remaining two conditions are about the disturbance temperature T and heat flux Pr −1 ∂T /∂z. Unlike the fixed TRBC whereT and DT are continuous at z = 0, § § 3.2 and 3.3 will show that it is not the case for the free TRBC due to the interface undulation and the non-uniform evaporation there. The full solution procedure of the eigenvalue problem is shown in Appendix B.
3.2. The effect of undulating interface The first term on the right-hand side of (2.23) denotes the interface undulation. First, we study how the small-amplitude interface displacementẑ s produces a jump ofT at z = 0. Let be a small positive constant. The following double integral is performed on (2.23), within a small slot encompassing the interface. All terms vanish except the vertical 928 A13-11 H. Fu diffusion term which has the highest order, and the dδ(z)/dz term which changes most abruptly: Another more straightforward way to derive (3.6) is to use the interfacial temperature continuity T| z=z + s = T| z=z − s and perform Taylor expansion of T near z = 0, as has been employed by Hsieh (1972) in studying liquid-gas flow. The parameter Pr appears on both terms of the second line of (3.6) and can be eliminated.
Then, we consider the Clausius-Clapeyron equation that determines z s . At the saturation interface, there is q vs = q t . The normal mode form of (2.20) yieldŝ Fielder (1984) has derived an expression similar to (3.7) to find the cloud bottom undulation magnitude, where he considered the matching between a well-mixed sub-cloud layer and a well-mixed cloud layer. The normal mode form of (2.18) and (2.19) are manipulated to representT withŵ: (3.8) Substituting (3.7) and (3.8) into (3.6), we obtain the matching condition ofT: where η T is a negative constant. Theŵ| z=0 + andŵ| z=0 − on the second line of the right-hand side denote the derivative value ofŵ at z = 0 + and z = 0 − . The quantityq t | z=0 is obtained by solving q t as an advective-diffusive tracer driven by w. The normal mode form of (2.22) is a second-order ODE ofq t : It is solved with the method of variation of parameter in Appendix C. For σ = 0, (3.10) can be viewed as the normal mode form of a Poisson equation with −ŵPr as the source term. Thus, whenŵ is a real function (no phase tilt with height),q t is real and has the same sign asŵ. In other words, there is always positive q t | z=0 at the updraft region for a neutral mode due to vertical advection. The physical influence of the undulating interface is briefly analysed here, and illustrated in figure 2(a). The updraft at the saturation interface pushes the air at the interface upward. These parcels are the coldest air in the vertical air column. The temperature field adjusts to the rising dome through diffusion, so the result is a cold anomaly at the upper layer and a warm anomaly at the lower layer, and vice versa in a valley produced by a downdraft. This produces a vertical dipole buoyancy anomaly, which corresponds to two undulation-induced torques which are in the opposite sign. A qualitatively similar yet more singular pattern can be produced by the interface undulation in a three-layer Rayleigh-Taylor instability model (Mellado et al. 2009). We note that the torque below the interface is in the same direction as the lower-layer convective cell (also a torque), which is generated by the unstable stratification. Because the main body of the convective cell is below the interface, the cell receives a more positive influence from the interfacial torque right below the interface than the negative influence from the torque right above the interface. Thus, the interface undulation is a destabilizing factor. A stricter explanation based on energetics is given in § 4.5.

The effect of non-uniform evaporation at the interface
The second term on the right-hand side of (2.23) denotes the non-uniform liquid water diffusion and, therefore, evaporative cooling at the interface. It induces a jump of perturbation heat flux and, therefore, DT at z = 0. This term is the key to reproduce the deformation-induced anomalous evaporative cooling at a 'flame front' (the interfacial ridge lifted by the updraft) described by Siems et al. (1990) and Mellado et al. (2009). The author is unaware of any previous work that addresses the influence of the deformation-induced anomalous evaporation on the interfacial instability.
Equation (2.23) indicates that the matching condition for heat flux in normal mode form is (3.11) Substituting (3.8) into (3.11), we get the matching condition represented withŵ: As in § 3.2, theŵ| z=0 + andŵ| z=0 − on the left-hand side denote the value of the derivatives ofŵ at z = 0 + and z = 0 − . The expression of Dq t | z=0 is documented in Appendix C.
In § 4 the eigenfunctions of the neutral mode show that in most cases there is always ∂q l /∂z| z=0 − < 0 at the updraft region (w| z=0 > 0), primarily due to the squeezing of the q t contour surfaces near the interface. The detailed physical understanding and constraint are to be discussed in § 4.4. Figure 2(b) illustrates the physical influence of non-uniform evaporation. There is more evaporative cooling in the updraft region than in the downdraft region, so the baroclinic torque produced by the non-uniform evaporation is opposite to the torque produced by the lower-level unstable stratification. Thus, the non-uniform evaporation is stabilizing. It is worth noting that M only influences the instability via (3.11), not directly relevant to the interface undulation. Thus, a larger M directly means a stronger non-uniform evaporation. Sensitivity tests are performed around Ref-ER, with an emphasis on studying the sensitivity to Q rad /Q evap by changing λ and γ T . The Ref-E test can be regarded as a special sensitivity test of the Ref-ER test with a larger M, which directly means a stronger non-uniform interfacial evaporation. Though this parameter set does not correspond to a specific state of real stratocumulus which is turbulent, it is an idealized state where the convective instability in the saturated layer, interfacial instability and upper-layer stratification all play a role. Only the marginal stability curve is studied in this section to understand the physical mechanism. Thus, Ra is to be solved, rather than prescribed. In § 5 the weakly supercritical (linearly unstable) regime is briefly explored with nonlinear numerical simulation. The neutral mode might represent the large eddy component of a turbulent flow, with the small-scale eddies parameterized as the diffusivity (Zhou, Simon & Chow 2014;Thuburn & Efstathiou 2020). We leave the careful investigation of the relevance to the turbulent stratocumulus for future works.
The critical (marginally stable) Ra is defined as Ra c . The problem is horizontally isotropic, so there is no need to split the horizontal total wavenumber K into k x and k y . On a marginal curve, the minimum Ra c is defined as Ra cm and the corresponding K is K cm . For the single layer Rayleigh-Bénard convection at Pr = 1, the existence of minimum Ra c can be heuristically understood as the competition between the destabilizing effect of the inviscid unstable gravity wave component, and the stabilizing effect of the diffusive and viscous component (Thuburn & Efstathiou 2020). The supposed growth rate of the unstable gravity wave increases mildly with K and asymptotically approaches (Ra/Pr) 1/2 , which is the modulus of the non-dimensional imaginary Brunt-Väisälä frequency. The damping rate of diffusion increases as K 2 . Thus, Ra cm must be achieved at a finite K = K cm . It will be shown that in this problem, a smaller Ra cm is mostly accompanied with a smaller K cm , similar to the fixed TRBC problem (Ogura & Kondo 1970). Qualitatively speaking, a smaller Ra cm means a smaller imaginary Brunt-Väisälä frequency, so less diffusion and therefore a lower K cm is needed to balance it.
By calculating the growth rate of various Ra in the parameter space around the reference test, we have not found any growing or neutral oscillatory mode (σ r ≥ 0, σ i / = 0). Thus, we let the solver (introduced in Appendix B) only search across stationary modes (σ i = 0) in all the tests. Whitehead (1968) has shown that an unstable or neutral oscillatory mode does not exist for the fixed TRBC problem (the principle of exchange of stability), but a rigorous proof is hard for our free TRBC problem. For the stationary mode, we arbitrarily set the phase by lettingŵ be a real function with positive value in the lower layer. With this setting, the normal mode form of (2.18)-(2.23) show that all the eigenfunctions in (3.1) are real functions exceptû(z) andv(z). These real-value eigenfunctions are proportional to the physical space variables at an updraft column. In Appendix D we introduce the novel finding that the purely evaporative cooling case yields an identical marginal stability curve (for the neutral mode) and growth rate (for the growing mode) to the corresponding fixed TRBC case (using the same γ T and Pr but forT| z=0 + =T| z=0 − and DT| z=0 + = DT| z=0 − ). In § 4.5 it is explained as the cancellation between the effect of the interface undulation and non-uniform evaporation with an argument of energetics. Whether it is flatter than the q t contour lines is hard to see from this figure due to a 'technical issue' in the plotting: too large an amplitude violates the small-amplitude assumption, and too small an amplitude is hard to observe. This point will be carefully discussed in § 4.3 and the interface will be shown to be no steeper. As there is Dq l | z=0 − < 0 in figure 4 for the Ref-E and Ref-ER tests, the perturbation liquid water gradient is negative at the updraft. This corresponds to a larger q l gradient at the updraft ridge of the interface, which produces more evaporation and explains the stabilizing effect of the non-uniform evaporation.

The dynamics of interface undulation
In the absence of radiative cooling (Q rad /Q evap = 0) the interface has constant T = 0 and q t = 0 (ẑ s =q t | z=0 for any λ), a property due to the synchronous evolution of q t and enthalpy h = (c p T This was pointed out in an equivalent thermodynamic framework by Mellado et al. (2009), and is derived with our thermodynamic framework in Appendix D. Asq t | z=0 must be positive at the updraft, as has been discussed in § 3.2, the saturation interface should rise at the updraft due to the moisture increment.
What is the physical mechanism for T to be fixed in the Q rad /Q evap = 0 case? We explain it as a competition: for an updraft, the warming tendency due to the change of interface height is balanced by the cooling tendency due to the extra evaporative cooling at the interface. The warming tendency can be understood with an idealized purely diffusive problem of T, with the interface cooling rate fixed. Suppose z s moves upward at an updraft due to interface undulation. As the upper level is infinitely thick, the vertical thermal diffusion strongly adjusts the upper level temperature profile to the new diffusive-equilibrium profile, which makes T| z=z s rise. The full adjustment is impossible due to the damping by horizontal diffusion.
What if radiative cooling is considered? Equation (2.16) shows that the Q rad /Q evap is larger for a smaller M, larger λ or larger −γ T . The suppression of the instability by a larger M has been discussed at the end of § 3.3, so we focus on λ and γ T . We hypothesize that the addition of radiative cooling leads to a stronger interface undulation relative to the non-uniform evaporation, so the anomalous cooling at the updraft ridge of the interface is weaker compared with the warming tendency. Thus, the interfacial temperature T| z=z s should rise more above the basic state value T = 0, and the Clausius-Clapeyron equation indicates that the saturation interface should fall behind the q t = 0 contour line (ẑ s /q t | z=0 < 1). The hypothesis is confirmed in figure 5(a), which shows thatẑ s /q t | z=0 decreases as λ or −γ T increases. In particular, there isẑ s =q t | z=0 for the λ = 0 case, because the Clausius-Clapeyron equation (3.7) is independent of temperature when λ = 0. Theẑ s =q t | z=0 is a property shared by the Q rad /Q evap = 0 (yet λ / = 0) case. The analysis leads to an inference: even though radiative cooling enhances interface undulation, the decrease ofẑ s /q t | z=0 means that it is self-limiting.
In fact, a quantitative lower bound ofẑ s /q t | z=0 can be found. Note thatẑ s can be expressed with eitherT z=0 − orT z=0 + in (3.7). As long asẑ s is positive (the interface rises at the updraft), both the movement of the interfacial cooling source and the direct vertical advection of the basic state temperature makeT z=0 − 0 andT z=0 + 0. Substituting this argument into (3.7), and combining it with the physical argumentẑ s q t | z=0 , we obtain a range estimation ofẑ s :q The lower bound is plotted as the dashed lines in figure 5(a). Though this constraint is loose, the lower bound does grasp the decrease ofẑ s with λ, as well as a faster decrease for a higher −γ T (with fixed M and Pr).  the updraft intersects the interface, so the q t contour surfaces are squeezed and there is a negative ∂q t /∂z, as has been discussed in § 3.3. The rising of the interface cooling source also induces a warming anomaly right below the interface (e.g. figure 2a), which locates above the peak warming anomaly induced by the upward advection of the unstable basic state. Thus, the negative DT| z=0 − should rise towards zero, as can be seen by comparing figures 4(c) and 4(e).

The liquid water gradient near the interface
This can be partly quantified. We start from the Q rad /Q evap = 0 case introduced in Appendix D, where (D7) showsq t =T for z < 0. This indicates When radiative cooling is incorporated and gradually raised by increasing −γ T , thê T jump at z = 0 should increase according to (3.6). Thus, we expect that as −γ T increases, Dq l | z=0 − will approach Dq t | z=0 . This argument and (4.2) lead to a constraint of Dq l | z=0 − /Dq t | z=0 : For a −γ T as large as 20 (very strong interfacial undulation), DT| z=0 − only marginally exceeds 0 and Dq l | z=0 − /Dq t | z=0 marginally exceeds 1. We have not figured out the mechanism of this asymptotic behaviour. Equation (4.3) is verified in figure 5(b) which shows that Dq l | z=0 − /Dq t | z=0 does increase mildly with increasing Q rad /Q evap (performed by changing γ T and fixing λ), with (1 − λ) as its lower bound. This indicates a coupling of the two effects: a stronger interface undulation tends to enhance the non-uniform evaporation to some extent. Unfortunately, (4.3) is not accurate enough to include γ T . With (4.1) and (4.3) in hand, we try to answer whether λ is a destabilizing factor or not. First, we look at the marginal curve of different λ for fixed γ T . Figure 3(b) shows that for the reference value γ T = −2.5, a larger λ significantly destabilizes the system. For double γ T , however, a larger λ weakly stabilizes the system. This indicates that a larger −γ T (stronger upper-layer stratification) gradually shifts λ from a destabilizing factor to a stabilizing one. Physically, this is due to two competing factors. For given −γ T , a larger λ suppressesẑ s and therefore weakens the interface undulation as shown in (4.1), but it also weakens the non-uniform evaporation as shown in (4.3) due to a smaller magnitude of liquid water gradient. The dependence of the relative strength of the two factors on γ T is not quantitatively constrained at present.

The relative importance of interface undulation and non-uniform evaporation
The relative strength of the destabilizing effect of the undulating interface and the stabilizing effect of non-uniform evaporation is quantified by comparing their relative contribution to the tendency of the vertical buoyancy flux wT . Here the overline on perturbation quantities denotes horizontal average in a wavelength. Multiplying w to (2.23), multiplying T to the w component of (2.19), and then summing them up and performing a horizontal wavelength average, we get where 'IU' denotes the contribution from interface undulation and 'NE' denotes that from non-uniform evaporation. Other terms on the right-hand side of the buoyancy flux tendency equation (4.4), which do not involve the interfacial process, are not shown. The relative importance of the two mechanisms can be roughly estimated as the ratio of the vertical integral of w(∂T/∂t) IU to w(∂T/∂t) NE . As their contributions are opposite in sign, their ratio with a minus sign is defined as a positive quantity Φ, To make the undulating interface render a positive buoyancy flux tendency and favour the convective cell, the interface must locate at the upper part of the convective cell (Dŵ| z=0 < 0), so that theŵ at the heating (lower) part of the interfacial dipole is larger than that at the cooling (upper) part, as shown in (4.6) and illustrated in figure 2(a). Substituting the estimate ofẑ s (4.1) and that of Dq l | z=0 − (4.3) into (4.6), we get .
(i) The destabilizing effect of the undulating interface can be stronger than the stabilizing effect of non-uniform evaporation only when the radiative cooling is present. (ii) The upper bound of the net destabilizing effect is roughly proportional to the ratio of the basic state interfacial radiative cooling rate to the diffusion-induced evaporative cooling rate. It confirms that raising Q rad /Q evap by decreasing M (fixing λ and γ T ) makes the system more unstable. (iii) The energetics argument is unable to depict the change of the instability behaviour on λ for different γ T (e.g.by calculating ∂ 2 Φ/(∂λ∂γ T )), which is discussed in § 4.4. This is due to a lack of the way to include γ T into the inequality (4.3).
Another invariant property is that changing Pr alone does not influence the marginal curve and the eigenfunctions at all. The former is illustrated in figure 3(c). Both properties are evident by checking the sixth-order ODE in (3.2) and the matching conditions in (B1), or the heuristic comparison in (4.7). To explain this, first we note that Pr itself does not influence the neutral mode of the fixed TRBC where ν and κ are equivalent dissipating factors. Second, we note that Pr does not influence the moist process. This originates from the assumption of identical thermal and water diffusivity, which makes the interfacial moist processes detached from Pr. For Ra > Ra c , Pr still has an influence on the growth rate, because the time evolution of theq t (3.10) involves Pr.

Nonlinear numerical simulation
In this section fully nonlinear numerical simulations of the non-dimensional version of (2.1)-(2.6) are used to validate the linear stability analysis, and demonstrate the flow pattern in the weakly nonlinear regime.
The simulation code is built by the author using a standard two-dimensional Fourier spectral method (Durran 2010). The governing equation is in vorticity-stream function formulation, coupled to an equation of q t and an equation of enthalpy h = T + Mq v . The h equation is ∂h ∂t + u · ∇h = κ∇ 2 h + Q rad δ(z − z s ), (5.1) whose only source term is radiative cooling. In the code Q rad δ(z − z s ) is regularized to Q rad /(π 1/2 l rad ) exp[−(z − z s ) 2 /l 2 rad ], where l rad = 0.05 is a fixed regularization length scale. Using (2.4), the temperature T is diagnosed from q t and h with a maximum operator, which is equivalent to the 'dry and moist buoyancy' formulation of Bretherton (1987): The position where the two quantities equal each other is z s . This formulation avoids the direct treatment of evaporative cooling, so it is convenient for numerical simulation. Imposing a material derivative on (5.2) and then expressing the maximum operator with Heaviside function yields the non-dimensional version of (2.6), which has explicit evaporative cooling and is therefore more suitable for theoretical analysis. A note of the derivation is deposited in the supplemental material available at https://doi.org/10.1017/ jfm.2021.784. In the code, the α − softmax function is used as a smoothed maximum operator: The same as (a), but for the Growth-ER tests. (c) The growth rate against wavenumber K. The theoretical calculation for the Growth-E tests is shown as the solid blue line, and that diagnosed from the nonlinear numerical simulation is shown as the blue '+.' Those corresponding to the Growth-ER tests are shown as the solid red line and red '+.' The Growth-E test at K = 1 has a negative growth rate (very close to zero) and does not appear in this σ > 0 plotting. and φ 2 are arbitrary input variables. The variable z s , which must be accurately determined to calculate the radiative cooling term, is set by vertically interpolating φ 1 and φ 2 on a very fine mesh, and then using the unmodified maximum operator to find the zero-crossing point of φ 1 − φ 2 . We choose α = 20, which leads to a saturation interface thickness of about α −1 = 0.05 for O(φ 1 − φ 2 ) ∼ 1. The second-order Adams-Bashforth scheme is used in time stepping. Zero padding is used to accurately calculate the product terms of all resolved waves in the physical space. No hyperdiffusion is used.
The infinite upper boundary in the linear stability analysis is replaced by a rigid lid in the simulation, which is located at a distance l above the basic state interface z = 0. The length l must be much larger than the convective penetration depth H p = −γ −1 T to make the upper lid a good approximation of the infinite boundary. As a modification from (2.8), the upper lid has free-slip velocity, fixed h and q t , The parameter set is identical to the reference test shown in § 4.1, except that we use a supercritical Ra = 600 to study the unstable mode. The domain depth in the z-direction is 1 + l, with l = 3. The domain width in the x-direction is different from test to test. The vertical grid point number is 256, equivalent to a vertical grid interval of z = 0.0156. This guarantees at least three grid points in the radiative cooling slot whose thickness is l rad , and also three points in the saturation interface whose thickness is α −1 . The horizontal grid point number is 4 for the ten small-amplitude growth rate benchmark tests, and 64 for the two flow pattern tests which extend to the finite-amplitude regime. The time step is t = 0.041( z) 2 . The initial condition is the analytical solution of radiative-diffusive equilibrium (running for 0.1 time units with the advection term off), plus a smoothed noise on h.
Among the ten low-horizontal resolution simulations that aim to benchmark the growth rate, the first five runs, named Growth-E tests, are for the purely evaporative cooling case. They have the same parameter as the Ref-E test except for the given Ra. The horizontal domain width takes L x = 2π/K, where K = 1.0, 1.5, 2.0, 2.5 and 3.0. The second five runs, named Growth-ER tests, are for the mixed evaporative-radiative cooling tests which use the same array of K. The parameters are the same as the Ref-ER test except for the given Ra. The time series of the standard deviation of w are used as the norm for calculating the growth rate, as shown in figure 6(a,b). An exponentially growing range of the stationary mode is clear for each run, though there is some initial oscillation which is probably due to the decaying oscillatory mode. The growth rate of all simulations match well with the theoretical prediction of the horizontally longest wave, as shown in figure 6(c). The Growth-ER tests grow faster than the Growth-E tests, which is not surprising given that the Ref-ER test has a lower Ra cm than the Ref-E test. Two high-resolution tests with 64 grid points in the x-direction are performed for the Growth-E and Growth-ER parameters. Both tests have a domain width of L x = 2π/K, with K = 2.0 which is close to the most unstable wavenumber of each test (e.g. figure 6c). Their vertical profiles at t = 1.4 match well with the respective theoretical eigenfunctions of the Ra = 600 and K = 2.0 modes, as shown in figure 7. The result is not sensitive to the choice of t, so long as it is in the small-amplitude regime. Figure 8 shows that in the finite-amplitude regime of both tests, the cold air produced at the ridge of the interface (the 'flame front') falls towards the valley of the interface. This is similar to that reported in the local staircase initial state simulations of Siems et al. (1990) and Mellado et al. (2009), though the structure here is smoother due to the smoother basic state profile. The radiative cooling influences the growth rate in this framework, but it does not change the flow pattern dramatically.
In the nonlinear regime the moist process is no longer the simple balance between the interface undulation and the non-uniform evaporation. As suggested by one of the reviewers, a detailed comparison between the fixed TRBC and the free TRBC in the weakly nonlinear regime can be made. For the turbulent regime, the moist convection has the distinctive feature that the vortex ring associated with the 'flame front' can engulf the upper-layer air and produce further turbulence through evaporative cooling (Siems et al. 1990). However, in the turbulent regime the strong horizontal mixing might homogenize the horizontal cooling difference, as suggested by de Lozar & Mellado (2013). It is unclear whether the two mechanisms here are still important for the large eddy dynamics of the turbulent cloud.

Conclusion
We present a novel two-layer convective instability problem, with the lower layer thermally unstable and saturated to water vapour, and the upper layer stable and unsaturated. This is a prototype model for understanding the cloud or fog top mixing. Two novel mechanisms which compete with each other are found.
First, the undulation of the interface is a destabilizing effect. The interface as a radiative and evaporative cooling source induces a vertical dipole buoyancy anomaly, which is superposed on the upper part of the lower-layer convective cell and therefore favours further undulation.
Second, the non-uniform evaporation at the interface is a stabilizing effect. There is stronger interfacial evaporative cooling at the updraft region due to the larger q l (liquid water) gradient magnitude near the interface, and vice versa at the downdraft. At the updraft, this is induced by a larger q t gradient due to the squeezing of q t contour surfaces by the horizontally divergence flow of the lower-layer primary cell, as well as a more complicated temperature gradient adjustment.
The competition between the destabilizing effect of undulating interface and the stabilizing effect of non-uniform evaporation is quantified by calculating the ratio of their contribution to the tendency of the vertical buoyancy flux. Without radiative cooling, the two effects break even. With radiative cooling, Q rad /Q evap can roughly measure the net destabilizing effect of radiative cooling, which enhances the interface undulation but does not directly influence the non-uniform evaporation.
As for the flow pattern, the baroclinic torque generated near the interface makes the primary convective cell more confined in the lower layer and favours the secondary cell in the upper layer. The squeezed q t contour surfaces at the interfacial ridge can be viewed as a quantification of the 'flame front' structure seen in the cloud-top mixing simulation (Siems et al. 1990;Mellado et al. 2009), though the front is smoother due to the smoother basic state profile in this theory. The subsequent aggregation of cold air in the valley of the interface is grasped in the nonlinear simulation.
This research is far from complete. Future works may include the following.
(i) The proof of principle of exchange of stability for the case with non-zero radiative cooling, if it exists. (ii) Using numerical simulation to investigate the instability with initial state profiles that lie between the local staircase buoyancy profile of Siems et al. (1990) and our global piecewise linear one. It may better link the real cloud-top mixing problem and this simple model. An investigation of the weakly nonlinear regime, as has been done by Ogura & Yagihashi (1971) for the fixed TRBC, is a necessary intermediate step to link our linear theory to the turbulent regime. (iii) An extension to consider the influence of vapour and liquid water on buoyancy, as well as the droplet settling. The first two factors can be readily included into the current framework. Their coupling to droplet settling is a challenging topic. (iv) An extension to study the instability of the cloud bottom (e.g.mammatus cloud formation), which is a saturation interface that is susceptible to evaporative cooling and the heating of the surface-emitted longwave radiation (Garrett et al. 2010;Ravichandran, Meiburg & Govindarajan 2020). The coupling of the cloud bottom to the cloud-top instability is also an important topic. Appendix A. The representation of radiative cooling In this appendix a physical interpretation of the constant radiative cooling rate formulation used in (2.6) is presented. The non-uniform radiative cooling at the interface due to temperature difference is shown to be a very weak stabilizing factor. Only longwave radiation is explicitly considered. The solar radiation causes heating deep inside the saturated layer (Oliver, Lewellen & Williamson 1978). It is qualitatively considered as a part of the fixed temperature boundary condition at z * = −H which is a boundary heat source. The sharp longwave absorption coefficient contrast between the cloudy and clear air makes the radiative flux diverge near the cloud top, causing radiative cooling that is concentrated near the saturation interface z * s (x * , y * , t * ) (Larson, Kotenberg & Wood 2007;de Lozar & Mellado 2013). The typical longwave radiation penetration depth at the stratocumulus top is around 15 m (de Lozar & Mellado 2013). If the longwave radiation penetration depth in the saturated layer is much shorter than its depth, the saturation interface can be regarded as a blackbody. If we further assume the unsaturated layer is transparent, the radiative cooling can be represented as a Dirac-delta function cooling concentrated at the saturation interface that obeys the Stefan-Boltzmann law, Here Q * rad is a constant radiative flux density expressed in temperature (negative, unit: K s −1 m −2 ), σ S ≈ 5.67 × 10 −8 W m −2 K −4 is the Stefan-Boltzmann constant. The motivation of using the Dirac-delta function is to make the problem more analytically tractable. Equation (A1) is straightforward for radiative cooling in the unsaturated layer. The radiative cooling in the saturated layer is slower due to the compensation by condensation heating (de Lozar & Mellado 2015). However, the extra liquid water caused by radiative cooling, which is close to the interface, is swiftly evaporated by diffusion in this laminar problem. Thus, the net cooling is the same.
In the turbulent regime, de Lozar & Mellado (2013) set a horizontally uniform radiative cooling, which neglects both the undulation of the radiative cooling source and the inhomogeneity of the radiative flux density at the interface. This can be regarded as the fixed TRBC test in the turbulent regime. They invoked the Stefan-Boltzmann law to explain that the characteristic horizontal temperature difference at the interface of a real cloud is too small to produce any difference in emission, and the eddy mixing is fast enough to eliminate the tiny cooling difference. This physical judgment is verified by de Lozar & Mellado (2015) who included both effects in their simulation and found no significant difference. Though the scaling argument by de Lozar & Mellado (2013) is The constants K 1m and K 2m are determined by imposing the lower and upper boundary conditionq t | z=−1 = 0 andq t | z→∞ = 0 on each of the m = 1, 2, 3, 4, 5, 6 mode independently, and their expressions are (D11) The last equality of each line of the equations has used the sixth-order ODE (3.2). The modified matrix A with the new E m and F m is defined asÃ, and the new amplitude vector is f . Here we show that the free TRBC in this case has the identical marginal stability curve with the fixed TRBC. The idea is to manipulateÃ to show that for given Ra and K (and therefore p m ), det(A) = 0 and det(Ã) = 0 occur at the same time. The proof procedure is shown below.
(i) First, we introduce a diagonal matrix B = diag{( p 2 m − K 2 − Prσ )}, and rewritẽ Af = 0 asÃBB −1f = 0. (ii) Second, we note that the fifth and sixth row ofÃB is identical to the original first and second line ofÃ. Some row transformations onÃB can transform it to A, soÃB and A share the row space and have the same rank. As B is non-singular,Ã and A always have the same rank, so det(A) = 0 and det(Ã) = 0 must occur at the same time.
This indicatesf = B f . It means the only difference between the fixed TRBC and the free TRBC is the amplitude coefficients, which makes the eigenfunctions (e.g.ŵ) different. Figure 4 shows that, compared with the Ref-F test, the primary lower cell in the Ref-E test is more confined in the lower layer, and the secondary upper cell is stronger. This can be understood withf = B f . For the neutral mode where σ = 0, B can be viewed as a ∇ 2 operator which acts on the eigenfunction of the fixed TRBC. As ∇ 2 means coarsening, the eigenfunction of the free TRBC should be curvier.
As the marginal curve of the free and fixed TRBC always collapse, the candidate σ for any given Ra should be identical. Thus, the principle of exchange of stability which applies for the fixed TRBC (Whitehead 1968) also applies for this purely evaporative case. Note that the relation shown in this appendix also works for a growing mode.