# Exploring subtropical stratocumulus multiple equilibria using a mixed-layer
model

###### Abstract

Stratocumulus clouds cover about a fifth of Earth’s surface, and due to their
albedo and low-latitude location, they have a strong effect on Earth’s
radiation budget. Previous studies using Large Eddy Simulations have shown
that multiple equilibria (both cloud-covered and cloud-free states) exist as a
function of fixed-SST, with relevance to equatorward advected air masses.
Multiple equilibria have also been found as a function of atmospheric CO2,
with a subtropical SST nearly 10 K higher in the cloud-free state, and with
suggested relevance to warm-climate dynamics.
In this study, we use a mixed-layer model with an added surface energy balance
and the ability to simulate both the cloud-covered and cloud-free states to
study both types of multiple equilibria and the corresponding hysteresis. The
model’s simplicity and computational efficiency allow us to explore the
mechanisms critical to the stratocumulus cloud instability and hysteresis as
well as isolate key processes that allow for multiple equilibria via mechanism
denial experiments not possible with a full-complexity model. For the
hysteresis in fixed-SST, we find that decoupling can occur due to either
enhanced entrainment warming or a reduction in cloud-top longwave cooling. The
critical SST at which decoupling occurs is highly sensitive to precipitation
and entrainment parameterizations. In the CO2 hysteresis, decoupling occurs in
the simple model used even without the inclusion of an active SST and SST-
cloud cover feedbacks, and the width of the hysteresis displays the same
sensitivities as the fixed-SST case. Overall, the simple model analysis and
results motivate further studies using higher complexity models.

##  1 Introduction

Stratocumulus clouds are low-altitude clouds that form beneath sharp
temperature inversions, where convection between the surface and the cloud
layer is driven by cloud-top longwave cooling (Turton and Nicholls,, 1987;
Bretherton and Wyant,, 1997; Wood,, 2012). On an annual average, stratocumulus
clouds cover about 20% of Earth’s surface and are particularly concentrated
over subtropical eastern ocean margins (Warren et al.,, 1988). Their location,
combined with their prevalence and high albedo, means that they have a
considerable effect on Earth’s radiation budget, providing a net local forcing
of up to $-100$ W/m2 in subtropical regions (Klein and Hartmann,, 1993).

Cloud-top longwave cooling causes turbulent convective motions that homogenize
the boundary layer and couple surface moisture fluxes to the cloud layer.
These fluxes sustain the cloud layer against strong entrainment warming and
drying through mixing with the inversion layer above (Lilly,, 1968; Nicholls,,
1984; Bretherton and Wyant,, 1997). Turbulence is also enhanced by latent heat
fluxes from the surface and latent heat release within the cloud layer
(Bretherton and Wyant,, 1997). Convection between the surface and cloud layer
is interrupted when cloud-top longwave cooling weakens (Lilly,, 1968) or when
entrainment warming becomes too strong (Bretherton and Wyant,, 1997). Under
these conditions, the cloud layer is effectively cut off (decoupled) from its
surface moisture supply and the stratocumulus clouds give way to much less
reflective scattered cumulus clouds. The decrease in albedo associated with
this cloud transition can cause considerable surface warming (Schneider et
al.,, 2019).

Bellon and Geoffroy, (2016) found a hysteresis in stratocumulus cloud cover as
a function of fixed-SST using a Large Eddy Simulation (LES), where multiple
equilibria (that is, both cloudy and non-cloudy states) are possible between
SST values of 288 K and 293 K. The presence of multiple equilibria was
attributed to the cloud radiative effect, in which a very dry cloud would have
a low emissivity and weak longwave cooling while a moist cloud would have a
high emissivity and strong longwave cooling. The demise of the stratocumulus
equilibrium with increased SST was attributed to an enhancement of entrainment
drying from the inversion layer above the cloud. The enhanced entrainment is
driven at high prescribed SST by turbulence produced by increased latent heat
fluxes from the surface. This enhanced entrainment drying leads to a lower
liquid water path in the cloud layer and therefore weaker longwave cooling at
the cloud-top. This fixed-SST scenario ignores the cloud cover-SST feedback.

More recent work has investigated how stratocumulus cloud cover in the
presence of an interactive SST can change due to CO2-driven climate warming,
and experience bistability and hysteresis. Schneider et al., (2019) analyzed a
hysteresis of stratocumulus clouds as a function of atmospheric CO2. For CO2
values between 300 ppm and 1200 ppm, both cloud-covered (coupled) and cloud-
free (decoupled) states were possible. They attributed the breakup of the
stratocumulus equilibrium at high CO2 to the increased opacity of the free
troposphere at higher CO2, which increases the downwelling longwave radiation
and therefore weakens the net longwave cooling at the cloud-top, ceasing
convection. As the CO2 is then decreased, it needs to be lowered to CO2 values
smaller than the original critical CO2 value in order to regain the coupled
equilibrium. The abrupt breakup of the cloud layer due to rising CO2 levels
may be relevant to future climate change as well as to past periods of
hothouse climates such as the Eocene (56–48 Myr). This adds to previous work
on cloud feedbacks in warm climates, including a possible role for polar
stratospheric clouds (Sloan and Pollard,, 1998; Kirk-davidoff et al.,, 2002),
a convective Arctic cloud feedback that can lead to multiple equilibria and
keep the arctic ocean ice-free (Abbot and Tziperman, 2008b, ; Abbot and
Tziperman,, 2009), and low-clouds over land due to the advection of moist air
from the ocean, which can prevent the surface from reaching sub-freezing
temperatures during winter time (Cronin and Tziperman,, 2015).

Stratocumulus topped boundary layers (STBLs) have been studied with multiple
models, ranging from simple mixed-layer models (MLM) (Lilly,, 1968; Turton and
Nicholls,, 1987; Bretherton and Wyant,, 1997; Pelly and Belcher,, 2001) to
Large Eddy Simulations (LES, Uchida et al.,, 2010; Bretherton and Wyant,,
1997; Bellon and Geoffroy,, 2016; Schneider et al.,, 2019). Simple models have
not been used to investigate the hysteresis behavior found in LES as a
function of SST or CO2 (Schneider et al.,, 2019; Bellon and Geoffroy,, 2016),
because the assumption of a well-mixed boundary layer breaks down when the
layer decouples from the surface. However, the diurnal cycle of stratocumulus
cloud was studied by Turton and Nicholls, (1987), who used a “stacked” mixed-
layer approach to model decoupled solutions.

Our objective is to explore the processes behind the multiple equilibria and
hysteresis behavior of stratocumulus cloud layers as a function of both CO2
and SST. For this purpose, we adapt a mixed-layer model to include a surface
energy balance and the ability to model decoupled solutions. With this model,
we reproduce the instability that causes a break-up of stratocumulus clouds
with increasing SST following Bellon and Geoffroy, (2016) and in CO2 based on
Schneider et al., (2019). We then explore (i) what nonlinearities specifically
allow for multiple equilibria, and (ii) what factors influence the critical
SST or CO2 at which the cloud layer onset or breakup occurs. Using a simple
model allows us to both isolate key processes and efficiently explore a wide
parameter space.

In Section 2 we describe the modified MLM that includes a surface energy
balance and allows decoupled solutions. This model is then used to explore the
dynamics of multiple equilibria and hysteresis in CO2 and fixed-SST in Section
3. We conclude in Section 4.

##  2 Model Description

In this section, we describe the mixed-layer model including the interactive
ocean mixed-layer component. We assume that moist static energy,
$\text{MSE}=c_{p}T+Lq_{v}+gz$, and total (vapor plus liquid) moisture content,
$q_{t}=q_{v}+q_{l}$, are well mixed throughout the layer in the coupled state
and are therefore constant in height (Lilly,, 1968). In Section 22.1, we
describe the coupled mixed-layer model equations that follow Bretherton and
Wyant, (1997), with a surface energy balance added. In Section 22.3, we follow
Turton and Nicholls, (1987) and model decoupled solutions by stacking a
“surface mixed layer” beneath the “decoupled mixed layer”. The difference in
vertical structure of conserved variables is summarized in Figure 1.

Figure 1: Vertical profiles of total (vapor plus liquid) moisture, $q_{t}$,
and moist static energy, MSE, in both the coupled and decoupled states. The
blue patch represents the cloud layer within the mixed-layer which is
infinitely thin when the model is decoupled. Figure 2: A schematic of the
fluxes in the mixed-layer model where red arrows correspond to terms in the
energy budgets for the surface and mixed-layer(s). Purple arrows correspond to
terms in the moisture budget. When the model is decoupled, the cloud layer is
infinitely thin and there is no precipitation.

###  2.1 Modelling the coupled state

As in Bretherton and Wyant, (1997), the prognostic equations for moist static
energy in the mixed layer, MSEM, and total moisture in the mixed layer,
$q_{tM}$, are

$\displaystyle\frac{d\text{MSE}_{M}}{dt}$ $\displaystyle=-\frac{\partial
E}{\partial z}=-\frac{\partial}{\partial z}\left(\left\langle
w^{\prime}\text{MSE}^{\prime}\right\rangle+F_{R}/\rho_{a}\right)$ (1)
$\displaystyle\frac{dq_{tM}}{dt}$ $\displaystyle=-\frac{\partial W}{\partial
z}=-\frac{\partial}{\partial z}\left(\left\langle
w^{\prime}q_{tM}^{\prime}\right\rangle-F_{P}\right).$ (2)

where $E$ and $W$ are the total energy and moisture vertical fluxes. Angle
brackets indicate horizontally averaged fluxes and $F_{R}$ (W/m2) and $F_{P}$
(kg/kg times m/s) are the net upward radiative flux and the precipitation
flux, respectively. The conservation equations for the full mixed layer are
found by integrating Equations (1) and (2) from the surface to the top of the
layer, $z_{i}$.

Variable | Value | Description
---|---|---
$\rho_{a}$ | 1.225 kg/m3 | Air density
$C_{T}$ | $0.001(1+0.07V)$ | Transport parameter
V | 7.1 m/s | 10-m wind
$T_{io}$ | 290 K | Initial inversion temperature
$\Gamma_{i}$ | 6.5 K/km | Lapse rate in inversion layer
RHi | 0.4 | Relative humidity of inversion layer
C | 0.5 | Entrainment efficiency
$w^{inv}_{ed}$ | 0.0056 m/s | Decoupled cloud-top entrainment, CO2-hysteresis
$w^{inv}_{ed}$ | 0.002 m/s | Decoupled cloud-top entrainment, SST-hysteresis
$\tau_{\mathrm{precip}}$ | 2 hrs | Precipitation
$\Gamma_{p}$ |  $0.86\times 10^{-5}$ kg/kg 1/km | Evaporation rate beneath cloud base
$\tau_{\mathrm{export}}$ | 3.96 years | Moisture export
$\epsilon_{c0}$ | 0.45 | Minimum cloud layer emissivity
D | 1.66 | Diffusivity factor
$\kappa$ | 90.36 m2/kg | Absorption coefficient of water droplets
$A$ | 0.05 | Climate sensitivity to CO2
$B$ | 0.1 | Climate sensitivity to H2O
$\epsilon_{i0}$ | 0.5 | Minimum inversion layer emissivity
Div |  $6\times 10^{-6}$ 1/s | Large scale subsidence
$p_{o}$ | 1000 hPa | Surface pressure
S | 471 W/m2 | Diurnally averaged insolation
$\mathrm{OHU}$ | 10 W/m2 | Ocean heat uptake
$\alpha_{w}$ | 0.06 | Albedo of water
$\alpha_{c}$ | 0.6 | Maximum albedo of cloud
$\mu$ | 0.93 | Thermodynamic factor (Equation (27))
$w^{SML}_{e}$ | 0.013 m/s | Entrainment in SML, CO2-hysteresis
$w^{SML}_{e}$ | 0.003 m/s | Entrainment in SML, SST-hysteresis
$q^{DML}_{t0}$ | 8 g/kg | Initial DML total moisture, CO2-hysteresis
$q^{DML}_{t0}$ | 3.5 g/kg | Initial DML total moisture, SST-hysteresis
$z_{SML}$ | 200 m | Top of SML
$z_{id}$ | 1000 m | Decoupled inversion height
Table 1: Table of constants in order of appearance in the text.

In order for the MSE and total moisture in the mixed layer to remain
vertically uniform, the energy and moisture flux convergences appearing in the
above equation must be independent of height and therefore the corresponding
fluxes must be linear functions of height,

$\displaystyle E(z)$
$\displaystyle=\left(1-\frac{z}{z_{i}}\right)E(0)+\frac{z}{z_{i}}E(z_{i})$ (3)
$\displaystyle W(z)$
$\displaystyle=\left(1-\frac{z}{z_{i}}\right)W(0)+\frac{z}{z_{i}}W(z_{i}).$
(4)

Surface fluxes, defined at $z=0$, are calculated using bulk aerodynamic
formulas,

$\displaystyle E(0)$
$\displaystyle=C_{T}V(\text{MSE}_{s}-\text{MSE}_{M})+F_{R}(0)/\rho_{a}$ (5)
$\displaystyle W(0)$ $\displaystyle=C_{T}V(q_{s}-q_{tM})-F_{P}(0),$ (6)

where MSEs and $q_{s}$ are the saturation surface moist static energy and
specific humidity. At the surface, $F_{R}(0)$ is the net upwards radiation at
the surface and $F_{P}(0)$ is the precipitation that reaches the surface. At
the inversion above the mixed layer, entrainment mixing lead to energy and
moisture fluxes of the form,

$\displaystyle E(z_{i})$
$\displaystyle=-w^{inv}_{e}(\text{MSE}_{i}-\text{MSE}_{M})+F_{R}(z_{i})/\rho_{a}$
(7) $\displaystyle W(z_{i})$ $\displaystyle=-w^{inv}_{e}(q_{i}-q_{tM}),$ (8)

where MSEi and $q_{i}$ are the moist static energy and specific humidity in
the inversion layer. We assume the warm air in the inversion layer originates
from updrafts in the tropics and will warm with increasing CO2 levels
(Schneider et al.,, 2019) assuming a climate sensitivity of 3 °C. We therefore
prescribe the temperature of the inversion layer as a function of CO2 and
inversion height,

$T_{i}=T_{io}+3\log_{2}\left(\frac{\text{CO}_{2}}{280}\right)-\Gamma_{i}(z_{i}-z_{io}),$
(9)

where $T_{io}$ and $z_{io}$ are the temperature and altitude of the inversion
layer at 280 ppm. The lapse rate in the inversion layer, $\Gamma_{i}$, is
assumed to follow a moist adiabat. We assume a constant relative humidity of
the inversion layer, RHi, so that the specific humidity in the inversion layer
can be calculated from the Clausius-Clapeyron relationship.

Following Pelly and Belcher, (2001), we define the entrainment rate,
$w^{inv}_{e}$, of warm, dry air from above the inversion into the cloud layer
as

$\rho_{a}c_{p}w^{inv}_{e}=\begin{cases}C\frac{\Delta F_{R}-\mu L\Delta
F_{p}}{\Delta\theta^{i-c}_{v}}&\text{Coupled}\\\
\rho_{a}c_{p}w^{inv}_{ed}&\text{Decoupled},\end{cases}$ (10)

where $C$ is a prescribed constant entrainment efficiency, $\Delta F_{R}$ is
the net longwave cooling by the cloud layer, $\Delta F_{P}$ is the net
precipitation flux across the mixed-layer, and $\Delta\theta^{i-c}_{v}$ is the
jump in virtual potential temperature across the inversion. Entrainment is
driven by convective turbulence in the mixed-layer and is inhibited by a
stable stratification across the inversion layer. Therefore, sources and sinks
of turbulence (net longwave cooling and evaporative cooling of drizzle beneath
the cloud-base) appear in the numerator of Equation (10) and the inversion
strength is in the denominator. This parameterization is very similar to the
energy balance entrainment rate described in Bretherton and Wyant, (1997) and
the “minimum entrainment” closure by Lilly, (1968). When the cloud layer is
decoupled from the surface, we set the entrainment to a background value,
$w^{inv}_{ed}$.

The precipitation flux, $F_{P}$ is assumed constant throughout the cloud and
decreases linearly with height beneath the cloud base at a rate of
$\Gamma_{P}$, as in Turton and Nicholls, (1987). By defining a constant
timescale of precipitation, we write the precipitation flux as a function of
height, $z$, and cloud liquid water path, LWP,

$F_{p}(z)=\begin{cases}\frac{\text{LWP}}{\rho_{a}\tau_{\mathrm{precip}}}&z_{i}>z\geq
z_{b}\\\
\frac{\text{LWP}}{\rho_{a}\tau_{\mathrm{precip}}}-\Gamma_{p}(z_{b}-z)&z_{b}<z,\end{cases}$
(11)

where $z_{b}$ is the cloud-base. Combining Equations (1)–(8), we write the
time evolution of MSEM and $q_{tM}$ as,

$\displaystyle z_{i}\frac{d\text{MSE}_{M}}{dt}$
$\displaystyle=C_{T}V(\text{MSE}_{s}-\text{MSE}_{M})+w^{inv}_{e}(\text{MSE}_{i}-\text{MSE}_{M})-\Delta
F_{R}/\rho_{a}$ (12) $\displaystyle z_{i}\frac{dq_{tM}}{dt}$
$\displaystyle=C_{T}V(q_{s}-q_{tM})-w^{inv}_{e}(q_{tM}-q_{i})-\Delta
F_{P}-z_{i}\frac{q_{s}}{q_{so}}\frac{1}{\tau_{\mathrm{export}}}.$ (13)

The last term contains an added moisture export out of the subtropics and into
other regions of the form used in Schneider et al., (2019) that scales with
the saturation specific humidity at the surface, $q_{s}$ and has a constant
timescale of export, $\tau_{\mathrm{export}}$

The net longwave cooling across the layer, $\Delta F_{R}$, is the difference
in the net upwards longwave flux between the cloud-top and the surface. Since
there is no precipitation flux entering from above at the cloud-top, the net
precipitation flux across the layer is simply the precipitation flux at the
surface, $F_{P}(0)$. We use a simple radiation parameterization, illustrated
in Figure 2, rather than the two-stream approximation used by Bretherton and
Wyant, (1997). At the cloud-top, the net upward longwave radiation is,

$F_{R}(z_{i})=\epsilon_{c}\sigma T^{4}_{top}+(1-\epsilon_{c})\sigma
T_{s}^{4}-\epsilon_{i}\sigma T_{i}^{4},$ (14)

where $T_{top}$ is the temperature of the cloud top, found by making vertical
profiles of temperature using the MSEM and $q_{tM}$. The cloud base, $z_{b}$,
is defined as the altitude where the saturation point is reached. The
emissivity of the cloud, $\epsilon_{c}$, is found as a function of LWP by
assuming the cloud acts as a grey body,

$\epsilon_{c}=1-(1-\epsilon_{co})\exp(-D\kappa\text{LWP}),$ (15)

where $D$ is a nondimensional diffusivity factor and $\kappa$ is the
absorption coefficient of liquid water (McCaa,, 2003). In the absence of
liquid water (no cloud), the emissivity is set to a minimum value,
$\epsilon_{co}$, due to the CO2 and water vapor in the mixed layer.

The emissivity of the inversion layer is written as a function of CO2 and
$q_{i}$,

$\epsilon_{i}=\epsilon_{io}+A\log_{2}\left(\frac{\text{CO}_{2}}{280}\right)+B\log_{2}\left(\frac{q_{i}}{q_{io}}\right),$
(16)

where $A$ is tuned by running the full MLM model to give a surface warming of
1.5 °C for a doubling of CO2 when $B=0$. The value of $B$ is tuned to give a
surface warming of 3 °C for a doubling of CO2 when including the water vapor
feedback.

At the surface, the net upwards radiation is,

$F_{R}(0)=\sigma T_{s}^{4}-\epsilon_{c}\sigma
T_{b}^{4}-(1-\epsilon_{c})\epsilon_{i}\sigma T_{i}^{4},$ (17)

where $T_{b}$ is the temperature at the base of the cloud. Combining Equations
(14) and (17), the net longwave cooling across the mixed layer is,

$\Delta
F_{R}=\epsilon_{c}\sigma\left(T_{top}^{4}+T_{b}^{4}\right)-\epsilon_{c}\sigma
T_{s}^{4}-\epsilon_{c}\epsilon_{i}\sigma T_{i}^{4}.$ (18)

Finally, we write the time evolution of the inversion height, $z_{i}$, and the
surface temperature, $T_{s}$.

$\displaystyle\frac{dz_{i}}{dt}$ $\displaystyle=w^{inv}_{e}-(\text{Div})z_{i}$
(19) $\displaystyle C_{s}\frac{dT_{s}}{dt}$
$\displaystyle=S(1-\alpha)+\epsilon_{c}\sigma
T_{b}^{4}+(1-\epsilon_{c})\epsilon_{i}\sigma T_{i}^{4}-\sigma
T_{s}^{4}-\rho_{a}C_{T}V(\text{MSE}_{s}-\text{MSE}_{M})-\mathrm{OHU},$ (20)

where $S$ is the diurnally-averaged insolation in June in the subtropics and
$\alpha$ is the albedo of the cloud layer,

$\alpha=\alpha_{c}-(\alpha_{c}-\alpha_{w})\exp{(-D\kappa\text{LWP)}},$ (21)

a form qualitatively consistent with Gettelman, (2015). We set a maximum cloud
albedo, $\alpha_{c}$, and a minimum albedo representing the reflectivity of
the sea surface, $\alpha_{w}$. The second to last term in the surface energy
budget Equation (20) is a combination of sensible and latent heat fluxes and
the last term represents ocean heat uptake.

The model requires the calculation of the cloud base altitude, $z_{b}$, as
well as the temperature of the cloud base, $T_{b}$. These may be determined by
solving for $z_{b}$ and $T_{b}$ from the moist static energy,
$\mathrm{MSE}_{M}=c_{p}T_{b}+Lq_{tM}+gz_{b}$, and from the condition that
saturation is reached at the cloud base $z_{b}$, $q_{tM}=q^{*}(z_{b},T_{b})$,
where $q^{*}$ is the saturation specific humidity. If $z_{b}$ is greater than
$z_{i}$ (as in the decoupled, dry case), we set the cloud base equal to the
cloud top to represent an infinitely thin cloud layer. At the cloud top, the
specific humidity may be smaller than the total moisture, $q_{tM}$, when
liquid water exists, and we similarly solve for the temperature at the cloud
top, $T_{top}$, from the moist static energy,
$MSE_{M}=c_{p}T_{top}+L\min\left(q_{tM},q^{*}(z_{i},T_{top})\right)+gz_{i}$,
where $z_{i}$ is calculated using Equation (19).

###  2.2 Diagnosing Decoupling

To begin, we solve Equations (12), (13), (19), and (20) to obtain steady-state
values of MSEM, $q_{tM}$, $z_{i}$, and $T_{s}$. When a stratocumulus cloud
layer is coupled to the surface, kinetic energy is produced within the cloud
layer by longwave cooling and entrainment drying which cools parcels such that
they are negatively buoyant and sink through the mixed layer and to the
surface. Decoupling occurs when negative buoyancy fluxes (downward motion of
light parcels) convert turbulent kinetic energy to potential energy. The
weaker turbulence weakens the convective motions and inhibits the parcel from
reaching the surface (Bretherton,, 1997). This means that the surface is
decoupled from the mixed layer.

Following Bretherton and Wyant, (1997), we write the vertical structure of the
buoyancy flux as,

$\langle w^{\prime}b^{\prime}\rangle(z)=(g/s_{vo})\langle
w^{\prime}s_{v}^{\prime}\rangle(z),$ (22)

where $s_{v}$ is the virtual static energy (equivalent to virtual potential
temperature, a buoyancy variable),

$s_{v}=c_{p}(T+T_{ref}(0.61q_{v}-q_{l}))+gz.$

To obtain a profile of the buoyancy flux with height, we must write the
virtual static energy in terms of known variables ($\text{MSE}_{M}$ and
$q_{tM}$).

$s_{v}=\text{MSE}_{M}-q_{v}L\left(1-\frac{0.61c_{p}T_{ref}}{L}\right)-c_{p}T_{ref}q_{l}$
(23)

Decoupling occurs when there are negative buoyancy fluxes beneath the cloud
base. There $q_{l}=0$ and therefore $q_{tM}=q_{v}$. Defining
$\mu=\left(1-0.61c_{p}T_{ref}/{L}\right)$,

$\displaystyle s_{v}(z<z_{b})$ $\displaystyle=\text{MSE}_{M}-\mu Lq_{tM}$ (24)
$\displaystyle\langle w^{\prime}s_{v}^{\prime}\rangle(z<z_{b})$
$\displaystyle=\langle w^{\prime}\text{MSE}_{M}^{\prime}\rangle-\mu L\langle
w^{\prime}q_{tM}^{\prime}\rangle.$ (25)

Equation (25) gives the sub-cloud buoyancy flux in terms of fluxes of the
known variables in the model. Rearranging Equations (3–6), the sub-cloud
buoyancy flux can be written as,

$\displaystyle\langle w^{\prime}s_{v}^{\prime}\rangle(z=z_{b})$
$\displaystyle=\overbrace{(1-\frac{z_{b}}{z_{i}})C_{T}V(\text{MSE}_{s}-\text{MSE}_{M})}^{\text{Convective
heat flux}}+\overbrace{\frac{z_{b}}{z_{i}}\Delta F_{R}/\rho}^{\text{Net
longwave cooling}}$
$\displaystyle-\overbrace{\frac{z_{b}}{z_{i}}w^{inv}_{e}(\text{MSE}_{i}-\text{MSE}_{M})+\frac{z_{b}}{z_{i}}\mu
Lw^{inv}_{e}(q_{i}-q_{tM})}^{\text{Entrainment warming}}$ $\displaystyle-$
$\displaystyle\overbrace{(1-\frac{z_{b}}{z_{i}})C_{T}VL(q_{s}-q_{tM})}^{\text{Latent
heat flux}}-\overbrace{\mu L(F_{p}(z_{b})-(1-\frac{z_{b}}{z_{i}})\Delta
F_{p})}^{\text{Sub-cloud drizzle evaporative cooling}}.$ (26)

We invoke a strict decoupling criterion, similar to the minimal decoupling
criterion in Bretherton and Wyant, (1997), such that if there are negative
buoyancy fluxes beneath the cloud, the mixed layer decouples,

$\langle w^{\prime}s_{v}^{\prime}\rangle(z=z_{b})\ \
\begin{cases}>0&\text{Coupled}\\\ \leq 0&\text{Decoupled.}\end{cases}$ (27)

Looking at individual terms in the vertical flux profile, we can learn about
the physical processes that can lead to decoupling. For example, net longwave
cooling is a source of positive buoyancy flux and therefore a producer of
total kinetic energy (TKE). Similarly, evaporative cooling due to cloud-top
entrainment is a source of TKE. The convective heat flux at the cloud base
also provides TKE by encouraging mixing. Evaporation of drizzle beneath the
cloud destroys TKE by working to stabilize the stratification there. Cloud-top
entrainment warming combats cooling and also destroys TKE. The competition of
all of these terms determines the total sub-cloud buoyancy flux.

We check the decoupling criterion at every time step. If it is positive, we
say the mixed layer is coupled, and we accept the solution from Equations
(12), (13), (19), and (20). If the decoupling criterion is negative, the mixed
layer is decoupled and the well-mixed assumption is no longer valid, so we
transition to the stacked mixed-layer model described below.

###  2.3 Decoupled Model

In Figure 1, we show the model configuration for a decoupled state. Instead of
a single well-mixed layer, the model is split into two stacked mixed layers, a
surface mixed layer (SML) and a decoupled mixed layer (DML), following Turton
and Nicholls, (1987). The DML is isolated from direct surface fluxes and
receives surface moisture and energy only through mixing with the SML. Each
mixed-layer has its own energy and moisture budgets, similar to Equations (12)
and (13).

The SML is coupled to the surface and experiences mixing at the top with the
DML,

$\displaystyle z_{SML}\frac{d\text{MSE}_{SML}}{dt}$
$\displaystyle=C_{T}V(\text{MSE}_{s}-\text{MSE}_{SML})+w^{SML}_{e}(\text{MSE}_{DML}-\text{MSE}_{SML})$
(28) $\displaystyle z_{SML}\frac{dq_{SML}}{dt}$
$\displaystyle=C_{T}V(q_{s}-q_{SML})-w^{SML}_{e}(q_{SML}-q_{DML})-z_{SML}\frac{q_{s}}{q_{so}}\frac{1}{\tau_{\mathrm{export}}},$
(29)

where $w^{SML}_{e}$ is the entrainment rate between the SML and the DML, set
as a constant in this model. We assume a constant thickness of the SML,
$z_{SML}$. Note that we also include the moisture export term out of the SML
as in the coupled state.

The energy and moisture budgets for the DML are,

$\displaystyle(z_{id}-z_{SML})\frac{d\text{MSE}_{DML}}{dt}$
$\displaystyle=w^{inv}_{e}(\text{MSE}_{i}-\text{MSE}_{DML})$
$\displaystyle-w^{SML}_{e}(\text{MSE}_{DML}-\text{MSE}_{SML})-\Delta
F_{R}/\rho_{a}$ (30) $\displaystyle(z_{id}-z_{SML})\frac{dq_{DML}}{dt}$
$\displaystyle=w^{SML}_{e}(q_{SML}-q_{DML})-w^{inv}_{e}(q_{DML}-q_{i})$
$\displaystyle-(z_{i}-z_{SML})\frac{q_{DML}}{q_{DMLo}}\frac{1}{\tau_{\mathrm{export}}}-\Delta
F_{P},$ (31)

where $z_{id}$ is the top of the DML, set as a constant in a decoupled case,
and $z_{id}-z_{SML}$ is the thickness of the DML. The precipitation flux,
$\Delta F_{P}$, is calculated as in Equation (11), although we note that in
all the results shown in this paper there is no liquid water in the DML, so
the precipitation flux is calculated to be zero.

When the mixed layer is decoupled from the surface, there is stable
stratification between the SML and DML which inhibits the moisture and heat
fluxes from the surface from reaching the DML. In other words, the
stratification,

$\Delta\theta_{v}=\theta^{SML}_{v}-\theta^{DML}_{v}>0.$ (32)

At every time step, we check the sign of $\Delta\theta_{v}$, and if it is
positive we say that the mixed layer is decoupled. When $\Delta\theta_{v}\leq
0$, the mixed layer recouples and Equations (12), (13) are used.

The surface energy balance is nearly identical to Equation (20), except now
sensible and latent heat fluxes are calculated between the surface and the
SML,

$C_{s}\frac{dT_{s}}{dt}=S(1-\alpha)+\epsilon_{c}\sigma
T_{b}^{4}+(1-\epsilon_{c})\epsilon_{i}\sigma T_{i}^{4}-\sigma
T_{s}^{4}-\rho_{a}C_{T}V(\text{MSE}_{s}-\text{MSE}_{SML})-\mathrm{OHU}.$ (33)

##  3 Results

In this section, we present and analyze the results of the mixed-layer model
for bi-stability as a function of CO2 (Section 33.1) and as a function of
fixed-SST (Section 33.2). In both cases, we explore the factors that determine
where the critical CO2 or SST that leads to coupling and decoupling occurs,
and what is the minimum physics that allows for multiple equilibria. We do
this using mechanism-denial experiments, turning off various physics and
examining the results. Throughout, we emphasize new lessons learned thanks to
the use of the simple model. We discuss in the conclusions the value and
expected applicability of these new lessons given the simplicity of the model.

###  3.1 Hysteresis in CO2

We solve Equations (12), (13) by setting the time derivative to zero and using
a root solver to find steady states of the prognostic variables as a function
of CO2. If the decoupling criteria (Equation (27)) is negative, we assume the
mixed-layer is decoupled, and we solve Equations (28)–(31) to obtain the
steady state of $q_{DML}$, $\text{MSE}_{DML}$, and $T_{s}$ and we check that
$\Delta\theta_{v}$ is positive.

Figure 3 shows several variables in the CO2 hysteresis. We solve the system
for CO2 values ranging from 280 ppm to 2500 ppm. Coupled solutions (plotted in
red) are characterized by lower surface temperatures due to the albedo effect
of the cloud layer, as in Schneider et al., (2019). Additionally, cloud-top
longwave cooling is higher in the coupled case, driving convection and
providing the cloud layer with a moisture source from the surface (Lilly,,
1968; Bretherton and Wyant,, 1997).

As CO2 increases, the emissivity of the inversion layer ($\epsilon_{i}$)
increases, both directly from its dependence on CO2 as well as from the water
vapor feedback (Equation (16)). As in Schneider et al., (2019), this causes a
decrease in net longwave cooling across the mixed layer, $\Delta F_{R}$
(Figure 3B), until at around 1750 ppm, the longwave cooling is no longer
strong enough to drive convection and the model predicts only decoupled
solutions at higher CO2.

Figure 3: Hysteresis in CO2. Red dots indicate coupled solutions while blue
dots indicate decoupled solutions. A. Sea surface temperature (SST, $T_{s}$).
B. Net longwave cooling across the mixed layer. C. Surface latent heat flux.
D. Decoupling/coupling criteria (Equation (27) for coupled solutions and the
negative of Equation (32) for decoupled solutions). E. Total albedo, including
cloud and surface, and cloud emissivity. F. Cloud liquid water path. G. Height
of mixed-layer top (inversion height, $z_{i}$). H. Cloud base height, where in
the decoupled case when there is no cloud, the cloud base is equal to the
cloud top. I. Cloud-top entrainment rate.

The entrainment closure (10) predicts that, all other things being equal,
entrainment should weaken with longwave cooling. However, in Figure 3I, we see
a slight strengthening of entrainment with CO2. This is because a warmer
mixed-layer leads to a weaker inversion strength, $\Delta\theta^{i-c}_{v}$,
which allows for stronger mixing between the cloud-top and the inversion
layer.

Figure 3D plots the decoupling criterion in red and the coupling criterion in
blue. At low CO2, the coupled solution predicts positive buoyancy fluxes
beneath the cloud from Equation (27), but as CO2 increases, the buoyancy
fluxes beneath the cloud approach zero. Since we have enacted a strict
decoupling criterion, negative buoyancy fluxes beneath the cloud lead to
decoupling, so once the buoyancy flux reaches zero only decoupled solutions
exist.

In the decoupled state, without a moisture source from the surface, the LWP of
the DML is zero and the cloud base is equal to the cloud top (Figures 3G & H).
This causes a low albedo ($\alpha$) and cloud emissivity ($\epsilon_{c}$),
leading to higher surface temperatures and weaker longwave cooling (Figure
3E). Once the mixed layer is decoupled, decreasing CO2 does not lead to re-
coupling at 1750 ppm. Instead, CO2 must be decreased down to 500 ppm in order
for the stratification between the SML and DML to become unstable and allow
for convection between the surface and DML. In other words, for CO2 between
500 and 1750 ppm, there are multiple equilibria in which the layer could be
either coupled or decoupled depending on its initial condition.

In general, our mixed-layer model yields similar results to those of Schneider
et al., (2019), who used an LES. In both this model and the LES, increasing
CO2 causes i) an increase in latent heat release from the surface due to
enhanced evaporation, ii) a decrease in longwave cooling due to a more opaque
free troposphere, iii) a decrease in LWP, iv) a decoupling of the cloud layer
at high CO2, and v) an abrupt increase in SST following decoupling. In their
hysteresis in CO2, Schneider et al., (2019) cited the importance of amplifying
cloud cover-SST feedbacks as well as the water vapor feedback for the abrupt
stratocumulus cloud breakup. To test this assertion as well as explore the
importance of other processes such as entrainment, we perform a suite of
mechanism denial experiments in Figure 4.

In the experiments shown in the first row of Figure 4, we create a non-
interactive SST by making $T_{s}$ a function of CO2 only to remove cloud
cover-SST feedbacks,

$T_{s}=T_{so}+3\log_{2}\left(\frac{\text{CO}_{2}}{280}\right).$

The solution is very similar to the baseline solution in Figure 3, plotted as
faint lines in Figure 4 for easy comparison. The weakening of longwave cooling
and strengthening of cloud-top entrainment are still seen with increasing CO2,
and decoupling occurs at high CO2, with multiple equilibria. This indicates
that amplifying cloud cover-SST feedbacks is not essential for decoupling in
this model.

In the experiments shown in the second row of Figure 4, we restore the
interactive SST but remove the contribution of water vapor to the inversion
emissivity by setting $B=0$ in Equation (16). The most notable difference
between the no water vapor feedback case and the baseline case is that it
takes over twice as much CO2 for decoupling to occur. As in Schneider et al.,
(2019), at 1000 ppm, the reduction in longwave cooling is only half as much in
the no water vapor feedback case than in the baseline case. However, abrupt
decoupling does still occur, just at a higher CO2 than the baseline case.
Therefore, in this model the water vapor feedback affects the timing of
decoupling but is not necessary for decoupling to occur.

Figure 4: Summary of mechanism denial experiments for hysteresis in CO2. 1.
Non-interactive SST, only a function of CO2. 2. No water vapor feedback, where
inversion emissivity is only a function of CO2 (B = 0 in Equation (16)). 3.
Both 1 and 2. 4. Constant inversion emissivity, $\epsilon_{i}$. 5. Constant
cloud-top entrainment, $w^{inv}_{e}$. 6. Constant cloud emissivity,
$\epsilon_{c}=1$. 7. Constant cloud emissivity, $\epsilon_{c}=0.45$. Columns
plot A, SST. B, net longwave cooling. C, cloud-top entrainment rate. D,
decoupling criteria for each experiment. Transparent lines show the baseline
simulation from Figure 3.

In the experiment shown in the third row of Figure 4, we perform the above two
tests simultaneously, turning off both the interactive SST and the water vapor
feedback. We conclude that even with both cloud cover-SST feedbacks and the
water vapor feedback disabled, abrupt stratocumulus breakup still occurs in
this model. Interestingly, Schneider et al., (2019) cited cloud cover-SST
feedbacks and the water vapor feedback as crucial components of the abrupt
stratocumulus cloud break-up in their model, but in our model decoupling
occurs without either of these feedbacks. This suggests that these feedbacks
may not be critical for the abrupt breakup of stratocumulus clouds with
increasing CO2. However, we recognize that the simplicity of our mixed-layer
model may mean that these processes may still be important in the parameter
regime relevant to the presently observed climate. Further studies of the
importance of these feedbacks should be done with more realistic models.

In Schneider et al., (2019), decoupling at higher CO2 was attributed to the
reduction in longwave cooling caused by the increased emissivity of the
inversion layer. To test this assertion, in the fourth row of Figure 4, we set
the inversion emissivity to a constant ($A=B=0$ in Equation (16)). In panel B
in the fourth row of Figure 4, we see no reduction in longwave cooling which
translates to no decrease in the decoupling criteria in panel D. We ran the
model out to 10,000 ppm, and still, we found no critical CO2 where the model
decouples and the stratocumulus clouds dissipate. This is consistent with the
results from Schneider et al., (2019) that stratocumulus breakup is caused by
the increase in emissivity of the free troposphere with CO2.

In the baseline solution, along with a weakening of longwave cooling with
increasing CO2, there is also a slight strengthening of the cloud-top
entrainment which could also lead to decoupling (Bretherton and Wyant,, 1997).
In the fifth row of Figure 4, we set the cloud-top entrainment to a constant
and find very little difference between this test and the baseline solution,
indicating that enhanced entrainment is not an important part of the dynamics
for a hysteresis in CO2 in our model. In all of the above-outlined sensitivity
tests, the existence of multiple equilibria remains robust. In their
hysteresis in SST, Bellon and Geoffroy, (2016) attributed the existence of
multiple equilibria to the cloud radiative effect. That is, clouds with high
LWP would have large longwave emissivity, leading to strong longwave cooling
which maintains the cloud layer. Clouds with a low LWP, on the other hand,
would have a small emissivity and weak longwave cooling, again maintaining the
low LWP. To test this, we set the cloud emissivity in our two final mechanism
denial experiments to a constant, independent of LWP. In the sixth row of
Figure 4, $\epsilon_{c}=1$ and in the seventh row of Figure 4,
$\epsilon_{c}=\epsilon_{co}$. For the high (low) emissivity case, only coupled
(decoupled) solutions exist because the longwave cooling is never small (high)
enough to support a decoupled (coupled) solution. Therefore, in our model the
dependence of cloud emissivity on liquid water path is necessary for multiple
equilibria to exist, consistent with Bellon and Geoffroy, (2016).

In the above experiments, we see that even processes not essential for
decoupling to occur can still affect the timing of decoupling. In Schneider et
al., (2019), the critical CO2 where decoupling occurs as well as the width of
the hysteresis was found to be very sensitive to model parameters such as
large-scale subsidence, but the computational cost of the model precluded
further sensitivity tests.

In Figure 5, we explore how the width of the hysteresis (that is, the width of
the SST or CO2 range for which multiple equilibria exist) varies with
uncertain model parameters. Of specific interest are the precipitation
timescale, $\tau_{\mathrm{precip}}$, and the entrainment efficiency, $C$, both
of which are ill-constrained because of the complexity of cloud microphysics
and aerosol-cloud interactions (Gerber,, 1996; Jiang et al.,, 2002; Bretherton
et al.,, 2007; Wang and Feingold,, 2009; Ackerman et al.,, 2009; Uchida et
al.,, 2010).

Figure 5: Critical CO2 where transition occurs when increasing CO2
(CO${}^{inc}_{2}$) and decreasing CO2 (CO${}^{dec}_{2}$) as a function of A,
precipitation timescale ($\tau_{\text{precip}}$). B, cloud-top entrainment
efficiency ($C$). C, entrainment between SML and DML ($w^{SML}_{e}$). D,
entrainment between DML and inversion layer ($w^{inv}_{ed}$). The difference
between the red and blue lines represents the width of the hysteresis.

Unsurprisingly, we find that the CO2 value where decoupling occurs when
increasing (CO${}^{inc}_{2}$) is very sensitive to both of these processes. In
Figure 5A, decreasing $\tau_{\mathrm{precip}}$ (increasing the precipitation
rate), leads to decoupling at a lower CO2. An order of magnitude change in
precipitation timescale leads to nearly a 2000 ppm increase in
CO${}^{inc}_{2}$. A decreased precipitation timescale leads to enhanced
evaporative cooling beneath the cloud base which stabilizes the layer and
inhibits convection (Stevens et al.,, 1998). Precipitation has no effect on
CO${}^{dec}_{2}$ because, in our model, LWP and therefore precipitation is
zero in a decoupled case. The timescale of precipitation in reality depends on
the density of cloud condensation nuclei (CCN) as well as the rate of cloud
droplet merging due for example to turbulent processes (Falkovich et al.,,
2002), both of which affect the average size of the water droplet (Ackerman et
al.,, 2009). Future work could investigate how these processes affect the
width of the hysteresis in a more realistic model.

The cloud-top entrainment efficiency, $C$, is another uncertain parameter
known to play an important role in simulations of stratocumulus cloud dynamics
using mixed layer models, therefore in Figure 5B, we test how $C$ affects the
width of the CO2-hysteresis. The model demonstrates strong sensitivity to the
entrainment efficiency, showing nearly a 2000 ppm decrease of CO${}^{inc}_{2}$
for a doubling of $C$. Once again, CO${}^{dec}_{2}$ shows no sensitivity,
since in decoupled solutions we set a constant background entrainment rate,
$w^{inv}_{ed}$. For models run under prescribed SST, cloud-top entrainment
from the inversion layer plays a key role in decoupling with increasing SST by
stabilizing the air column through entrainment warming (Bretherton and Wyant,,
1997). LES have been shown to overestimate entrainment rates compared to
observational best guesses, which has complicated efforts to quantify how
aerosols impact entrainment (Bretherton et al.,, 2007; Uchida et al.,, 2010).
When sedimentation feedbacks are included, high CCN concentrations lead to
more efficient entrainment that decreases cloud LWP, which can cancel up to
50% of the change in albedo due to changes in cloud droplet size (Twomey
effect, Uchida et al.,, 2010).

While our model does not include cloud microphysics, we can learn about its
effects on the hysteresis and multiple equilibria indirectly from our model
results, using knowledge of how CCN concentration affects precipitation and
entrainment. In a parameter regime with abundant CCNs, decreasing CCN
concentration should lead to an increase in precipitation efficiency
(corresponding to lower values of $\tau_{\mathrm{precip}}$ in Figure 5A), by
increasing the average size of cloud droplets (Twomey,, 1977; Ackerman et
al.,, 2009). In our model, this corresponds to a lower CO${}^{inc}_{2}$ as
seen in Fig. 5A. Additionally, decreasing CCN concentration could lead to a
decrease in entrainment efficiency (corresponding to lower values of $C$ in
Figure 5B) by creating larger cloud droplets that sediment out of the
entrainment zone rather than evaporating and causing evaporative cooling, and
lead to less efficient entrainment (Bretherton et al.,, 2007; Uchida et al.,,
2010). Based on the results of our model, this could lead to an increase in
CO${}^{inc}_{2}$. The question of which process would dominate for varying CCN
concentration is outside the scope of this study but would be interesting to
study using more complex models.

Our treatment of the uncoupled case with two stacked mixed-layers (surface and
decoupled mixed layers, denoted SML and DML, correspondingly), following
Turton and Nicholls, (1987), is highly idealized, especially with our added
simplifications that entrainment between the SML and DML as well as cloud-top
entrainment is constant. In Figure 5C and D, we explore how sensitive the
lower CO2 boundary of the multiple-equilibria region, CO${}^{dec}_{2}$, is to
these assumptions. Our model shows high sensitivity of CO${}^{dec}_{2}$ to the
entrainment rate between the SML and DML ($w^{SML}_{e}$) and comparatively
little sensitivity to the cloud-top entrainment ($w^{inv}_{ed}$). For
decreasing CO2, the critical value at which coupling occurs, CO${}^{dec}_{2}$,
occurs when the stratification between the SML and DML is no longer stable
($\theta_{DML}=\theta_{SML}$ at $z_{SML}$). This can be written as,

$T_{DML}\left(\frac{P_{SML}}{P_{o}}\right)^{R/c_{p}}(1+0.61q_{DML})=T_{SML}\left(\frac{P_{SML}}{P_{o}}\right)^{R/c_{p}}(1+0.61q_{SML}),$
(34)

where $T_{DML}$ is the temperature in the DML at $z_{SML}$ and $T_{SML}$ is
the temperature in the SML at $z_{SML}$, calculated from MSE and $q_{t}$ in
each layer. Increasing $w^{SML}_{e}$ moistens the DML, such that CO2 must be
decreased to very low levels in order for longwave cooling to become strong
enough to decrease $T_{DML}$ and have $\theta_{DML}\leq\theta_{SML}$. In an
uncoupled state, the DML is very dry and therefore cloud-top entrainment does
not lead to as much evaporative cooling and has a much smaller effect on the
dynamics of the DML compared to mixing with the SML. This leads to a much
smaller sensitivity to $w^{inv}_{ed}$ (compare Figs 5C,D). Again,
CO${}^{inc}_{2}$ does not depend on either entrainment rate since that
parameterization is only included in the uncoupled state.

###  3.2 Hysteresis in SST

Previous work by Bellon and Geoffroy, (2016) used an LES to demonstrate how
the cloud radiative effect allowed for multiple equilibria in fixed-SST. As a
test for our model, we replaced the surface energy balance with,

$T_{s}=T_{s}^{fixed},$ (35)

with $T_{s}^{fixed}$ representing the prescribed SST with which the model is
forced. CO2 is given a constant value of 280 ppm throughout this series of
model experiments. In previous sections, we demonstrated that increased
inversion emissivity was vital for decoupling to occur. In the model
configuration used above, the inversion temperature was a function of CO2 to
represent warming in the tropics due to CO2 changes. Therefore, for the fixed-
SST case, we rewrite the parameterization of inversion temperature as
$T_{i}=T_{s}-\Gamma_{i}(z_{i}-z_{io})$, an approximation consistent with
current observed relation between the inversion and surface temperatures
(Nicholls,, 1984). Once again assuming constant relative humidity,
$\epsilon_{i}$ now increases with SST through the water vapor feedback alone.

Figure 6 shows the results of our model in the fixed-SST case. Though the SST
is the same for both a coupled and decoupled case, there are multiple
equilibria between surface temperatures of 287 K and 305 K. In the LES study
of Bellon and Geoffroy, (2016), the SST range considered did not reveal the
entire hysteresis because the low SST at which coupling occurs was not
reached. Their SST range extended down to 288 K and in our model coupling
occurs at 286 K, suggesting that simply expanding their SST range may have
revealed the left side of the hysteresis. In addition, the critical SST where
decoupling occurred was around 294 K while ours is at 305 K. In order to
understand this 10 K difference, first we examine the mechanisms behind
decoupling in our mixed-layer model under prescribed increasing SST.

Starting from 280 K in a coupled state, increasing SST both weakens longwave
cooling (Figure 6B) due to the water vapor feedback and strengthens cloud-top
entrainment (Figure 6I), both of which lead to a decreasing decoupling
criterion (Figure 6D). As in the LES study of Bellon and Geoffroy, (2016), LWP
has a maximum at about 292 K and then rapidly decreases due to the enhanced
entrainment. The combination of enhanced entrainment and weakened longwave
cooling decreases the buoyancy fluxes at the cloud base (decoupling criterion,
Equation (27)), but at 305 K decoupling happens very abruptly. This is because
of the decrease in cloud emissivity (Figure 6E) when cloud LWP dips past a
threshold value set by the liquid water LW absorption coefficient, $\kappa$
appearing in Equation (15). This decrease in cloud emissivity results in a
large decrease in net longwave cooling and results in negative buoyancy fluxes
beneath the cloud base and therefore in decoupling. This effect was also noted
by Bellon and Geoffroy, (2016), who saw decoupling at 294 K. The difference
from Bellon and Geoffroy, (2016) in critical SST at which decoupling occurs
can be partially explained by the fact that their $\kappa$ was nearly half of
ours, leading to a reduction in cloud emissivity at a higher LWP. Other
differences likely exist in parameterizations of entrainment and
precipitation, as well as the difference in model complexity.

Figure 6: Hysteresis in SST. Red dots indicate coupled solutions while blue
dots indicate decoupled solutions. A. Sea surface temperature (SST, $T_{s}$).
B. Net longwave cooling across mixed layer. C. Latent heat flux from surface.
D. Decoupling/coupling criteria (Equation ((27)) for coupled solutions and the
negative of Equation ((32)) for decoupled solutions). E. Albedo and cloud
emissivity. F. Cloud liquid water path. G. Mixed-layer depth. H. Cloud base.
I. Cloud-top entrainment rate.

For example, Bellon and Geoffroy, (2016) prescribed the radiative effect of
clouds as a function of liquid water path alone, assuming constant fluxes from
the surface and free troposphere. This neglects the effect of the water vapor
feedback on the emissivity of the free troposphere, which may also contribute
to differences in our results.

As in the analysis of the CO2-hysteresis, we perform a series of mechanism
denial experiments, shown in Figure 7. As enhanced entrainment and weakening
longwave cooling can both decouple the mixed layer, we begin by shutting off
each of these mechanisms in turn. In the first row of Figure 7, we set the
entrainment to a constant and find that the weakening of longwave cooling due
to the water vapor feedback does not prevent decoupling and only shifts the
decoupling SST from 305 K to 307 K. Alternatively, in the second row of Figure
7, setting the emissivity of the inversion layer to a constant does not
inhibit decoupling either, as enhanced entrainment eventually decreases the
LWP such that longwave cooling dramatically decreases. It is only when both of
these mechanisms are shut off that the mixed layer stays coupled, even up to
320 K. This indicates that in the fixed-SST hysteresis, both enhanced
entrainment and weakening longwave cooling can cause decoupling on their own.

Figure 7: Summary of mechanism denial experiments for hysteresis in SST. 1.
Constant cloud-top entrainment, $w^{inv}_{e}$. 2. Constant inversion
emissivity, $\epsilon_{i}$. 3. Both 1 and 2. 4. Constant cloud emissivity,
$\epsilon_{c}=1$. 5. Constant cloud emissivity, $\epsilon_{c}=0.45$. Columns
plot A, SST. B, net longwave cooling. C, cloud-top entrainment rate. D,
decoupling criteria for each experiment. Transparent lines show the baseline
simulation from Figure 6.

In both cases, decoupling happens very abruptly due to the decrease of LWP
caused by both mechanisms (cloud-top entrainment and inversion layer
emissivity). Even when entrainment is constant, the weakening of longwave
cooling in the first row of Figure 7 leads to a warmer mixed layer, which
increases the saturation point and therefore decreases the liquid water ratio
in the mixed layer. Enhanced entrainment decreases LWP both by enhanced drying
but also by warming the mixed layer and allowing more water to exist in vapor
form. Abrupt decoupling occurs when the LWP passes a threshold at which the
emissivity of the cloud layer is no longer near-unity and longwave cooling
begins to weaken.

Finally, we test the assertion by Bellon and Geoffroy, (2016) that multiple
equilibria are caused by the cloud radiative effect by setting the emissivity
of the cloud to unity (experiments shown on fourth row of Figure 7) and to the
minimum value of $\epsilon_{c0}=0.45$ (fifth row of Figure 7), regardless of
LWP. As in the CO2-hysteresis, only coupled (decoupled) solutions exist when
$\epsilon_{c}=1$ ($\epsilon_{co}$). Once again, our model requires the
dependence of cloud emissivity on LWP to show multiple equilibria, consistent
with Bellon and Geoffroy, (2016).

Figure 8: Critical $T_{s}$ where transition occurs when increasing $T_{s}$
($T^{inc}_{s}$) and decreasing $T_{s}$ ($T^{dec}_{s}$) as a function of A,
precipitation timescale ($\tau_{\text{precip}}$). B, cloud-top entrainment
efficiency ($C$). C, entrainment between SML and DML ($w^{SML}_{e}$). D,
entrainment between DML and inversion layer ($w^{inv}_{ed}$). The difference
between the red and blue lines represents the width of the hysteresis.

We also use the computational efficiency of the model to study the sensitivity
of the SST hysteresis to various model parameters in Figure 8. In Figure 8A,
decoupling occurs at lower SST values when the precipitation timescale is very
small (precipitation rate is very fast), because of the enhanced evaporative
cooling beneath the cloud base. At higher $\tau_{\mathrm{precip}}$, the SST at
which decoupling occurs, $T^{inc}_{s}$, asymptotes to around 305 K because at
these timescales all of the precipitation evaporates before it reaches the
surface, and the net precipitation flux out of the mixed layer vanishes,
$\Delta F_{P}=0$. In Figure 8B, we vary the entrainment efficiency, $C$, and
find similarly to the CO2-hysteresis case that increasing $C$ leads to
decoupling at lower $T_{s}$ due to enhanced entrainment warming at the cloud-
top. In both Figures 8A, B $T^{dec}_{s}$ shows no sensitivity to either
parameter since there is no LWP and a constant cloud-top entrainment in the
uncoupled state.

In Figure 8C, we test the sensitivity of $T^{dec}_{s}$ to entrainment mixing
between the SML and the DML. Increasing $w^{SML}_{e}$ decreases $T^{dec}_{s}$
by moistening the DML and requiring very low surface temperatures to decrease
$\theta^{DML}_{v}$ beyond the threshold in (32), similar to the CO2-hystersis.
In Figure 8D, we test the sensitivity of the model to $w^{inv}_{ed}$. Unlike
in the CO2-hysteresis, there is a strong sensitivity to the cloud-top
entrainment rate due to its larger size relative to $w^{SML}_{e}$. The
hysteresis in CO2 and SST demonstrate similar sensitivities to these four
model parameters. Though the prescribed SST case ignores cloud cover-SST
feedbacks, it does provide an interesting case study for advected columns over
warm ocean water. In both cases, precipitation and cloud-top entrainment
parameterization play a large role in determining where decoupling occurs and
should be studied with higher complexity models.

##  4 Conclusions

Due to their wide coverage in the subtropics and strong albedo effect,
stratocumulus clouds are an important component of Earth’s energy budget.
Previous work by Bellon and Geoffroy, (2016) found multiple equilibria in
stratocumulus cloud cover as a function of fixed-SST, with applications to air
masses advected over meridional SST gradients. Similarly, Schneider et al.,
(2019) demonstrated how elevated CO2 concentrations could cause stratocumulus
clouds to dissipate with multiple equilibria in cloud cover. The dissipation
of the clouds corresponded to a local surface warming of up to 10 K. Both
studies used a Large Eddy Simulation (LES) approach, and the model complexity
precluded an in-depth study of important physical mechanisms and the
sensitivity to model parameters of the critical SST or CO2 that define the
range of multiple equilibria.

In this work, we developed a mixed-layer model adapted to include a surface
energy balance and the ability to model decoupled solutions using “stacked”
mixed-layers following Turton and Nicholls, (1987). Using this model, we were
able to reproduce the hysteresis as a function of SST (Bellon and Geoffroy,,
2016) and CO2 (Schneider et al.,, 2019). The simplicity of the model allowed
us to perform mechanism denial experiments by easily switching off various
physics. The computational efficiency of the model allowed us to perform a
suite of sensitivity tests to find which parameters control the hysteresis
behavior.

Schneider et al., (2019) deduced from their experiments that a dynamic SST,
cloud-SST feedbacks, and the water vapor feedback were crucial for the abrupt
stratocumulus cloud breakup to occur for higher CO2 values. Interestingly, we
find that neither cloud-cover SST feedbacks nor the water vapor feedback are
required for decoupling to occur in our model. Instead, only an increase in
inversion layer emissivity with CO2 is required. We found in extensive
sensitivity experiments that the main nonlinearity needed to support the
existence of multiple equilibria either as a function of fixed SST or of CO2
is provided by the dependence of cloud emissivity on liquid water path, as
suggested by Bellon and Geoffroy, (2016) for they SST-driven hysteresis. When
this nonlinearity is turned off the multiple equilibria cease to exist
regardless of the many other nonlinearities that are still present. The width
of the hysteresis (that is, the range of CO2 or SST in which multiple
equilibria exist) in our model also demonstrates strong sensitivity to
processes which depend in reality on cloud microphysics (precipitation and
entrainment parameterizations in our model). In their study of the hysteresis
in fixed-SST, Bellon and Geoffroy, (2016) attributed the decoupling and
stratocumulus dissipation at higher SST to the strong enhancement of
entrainment with increasing SST, but did not examine the role of the water
vapor feedback. We find that decoupling with increasing SST can occur either
via the enhanced entrainment or via increased inversion emissivity due to the
water vapor feedback in our model.

The model used here is highly idealized and is not meant for quantitative
predictions such as the precise SST or CO2 at which stratocumulus coupling or
decoupling occurs. Instead, it is designed as a tool to identify zeroth-order
physical mechanisms and dependencies and to explore a wide parameter space. It
is important to understand the applicability of the above-summarized lessons
learned using the mechanism denial experiments based on an idealized model
such as used here. A finding that a certain process is or is not essential to
explain a given observation (e.g., that the water vapor feedback can cause
decoupling at higher SST, or that dynamic SST is not required for the
existence of multiple equilibria as a function of CO2) adds to the fundamental
qualitative understanding of the processes involved. A feedback that is found
not to be essential for decoupling to occur, for example, may still contribute
quantitatively by affecting the critical SST or CO2 at which decoupling
occurs, as shown above. The findings of the idealized model motivate further
experiments with more realistic models that can then address these issues more
quantitatively.

As discussed in Schneider et al., (2019), the stratocumulus cloud instability
with CO2 presents an interesting cloud mechanism for past hothouse climates
with elevated CO2 levels. This joins other cloud feedbacks such as polar
stratospheric clouds (Sloan and Pollard,, 1998; Kirk-davidoff et al.,, 2002),
convective winter-time Arctic clouds in a warmer climate that can keep the
arctic from freezing (Abbot and Tziperman, 2008a, ), or low clouds that can
suppress polar air formation and keep winter-time continental surface
temperatures above freezing Cronin and Tziperman, (2015) that have been
proposed as possible mechanisms for the Eocene equable climate period. A study
of how these cloud mechanisms interact under rising CO2 levels could provide a
mechanism for the equable climate problem. Such a study could benefit from
simple models like the one developed and employed here.

We thank Zhiming Kuang for very helpful comments and suggestions. ET thanks
the Weizmann Institute for its hospitality during parts of this work. The work
has been supported by the NSF Climate Dynamics program (joint NSF/NERC) grant
AGS-1924538.

## References

  * (1) Abbot, D. S. and Tziperman, E. (2008a).  A high latitude convective cloud feedback and equable climates.  Quarterly Journal of the Royal Meteorological Society, 134:165–185.
  * (2) Abbot, D. S. and Tziperman, E. (2008b).  Sea ice, high latitude convection, and equable climates.  Geophysical Research Letters, 35:L03702.
  * Abbot and Tziperman, (2009) Abbot, D. S. and Tziperman, E. (2009).  Controls on the activation and strength of a high-latitude convective cloud feedback.  Journal of the atmospheric sciences, 66(2):519–529.
  * Ackerman et al., (2009) Ackerman, A. S., VanZanten, M. C., Stevens, B., Savic-Jovcic, V., Bretherton, C. S., Chlond, A., Golaz, J.-C., Jiang, H., Khairoutdinov, M., Krueger, S. K., et al. (2009).  Large-eddy simulations of a drizzling, stratocumulus-topped marine boundary layer.  Monthly Weather Review, 137(3):1083–1110.
  * Bellon and Geoffroy, (2016) Bellon, G. and Geoffroy, O. (2016).  Stratocumulus radiative effect, multiple equilibria of the well-mixed boundary layer and transition to shallow convection.  Quarterly Journal of the Royal Meteorological Society, 142(697):1685–1696.
  * Bretherton, (1997) Bretherton, C. (1997).  Convection in stratocumulus-topped atmospheric boundary layers.  In The physics and parameterization of moist atmospheric convection, pages 127–142. Springer.
  * Bretherton et al., (2007) Bretherton, C., Blossey, P. N., and Uchida, J. (2007).  Cloud droplet sedimentation, entrainment efficiency, and subtropical stratocumulus albedo.  Geophysical research letters, 34(3).
  * Bretherton and Wyant, (1997) Bretherton, C. S. and Wyant, M. C. (1997).  Moisture transport, lower-tropospheric stability, and decoupling of cloud-topped boundary layers.  Journal of the atmospheric sciences, 54(1):148–167.
  * Cronin and Tziperman, (2015) Cronin, T. W. and Tziperman, E. (2015).  Low clouds suppress Arctic air formation and amplify high-latitude continental winter warming.  Proceedings of the National Academy of Sciences, 112(37):11490–11495, 10.1073/pnas.1510937112.
  * Falkovich et al., (2002) Falkovich, G., Fouxon, A., and Stepanov, M. (2002).  Acceleration of rain initiation by cloud turbulence.  Nature, 419(6903):151–154.
  * Gerber, (1996) Gerber, H. (1996).  Microphysics of marine stratocumulus clouds with two drizzle modes.  Journal of Atmospheric Sciences, 53(12):1649–1662.
  * Gettelman, (2015) Gettelman, A. (2015).  Putting the clouds back in aerosol–cloud interactions.  Atmospheric Chemistry and Physics, 15(21):12397–12411.
  * Jiang et al., (2002) Jiang, H., Feingold, G., and Cotton, W. R. (2002).  Simulations of aerosol-cloud-dynamical feedbacks resulting from entrainment of aerosol into the marine boundary layer during the atlantic stratocumulus transition experiment.  Journal of Geophysical Research: Atmospheres, 107(D24):AAC–20.
  * Kirk-davidoff et al., (2002) Kirk-davidoff, D. B., Schrag, D. P., and Anderson, J. G. (2002).  On the feedback of stratospheric clouds on polar climate.  Geophysical Research Letters, 29(11).
  * Klein and Hartmann, (1993) Klein, S. A. and Hartmann, D. L. (1993).  The seasonal cycle of low stratiform clouds.  Journal of Climate, 6(8):1587–1606.
  * Lilly, (1968) Lilly, D. K. (1968).  Models of cloud-topped mixed layers under a strong inversion.  Quarterly Journal of the Royal Meteorological Society, 94(401):292–309.
  * McCaa, (2003) McCaa, J. (2003).  Cloud emissivity.
  * Nicholls, (1984) Nicholls, S. (1984).  The dynamics of stratocumulus: Aircraft observations and comparisons with a mixed layer model.  Quarterly Journal of the Royal Meteorological Society, 110(466):783–820.
  * Pelly and Belcher, (2001) Pelly, J. and Belcher, S. (2001).  A mixed-layer model of the well-mixed stratocumulus-topped boundary layer.  Boundary-layer meteorology, 100(1):171–187.
  * Schneider et al., (2019) Schneider, T., Kaul, C. M., and Pressel, K. G. (2019).  Possible climate transitions from breakup of stratocumulus decks under greenhouse warming.  Nature Geoscience, 12(3):163–167.
  * Sloan and Pollard, (1998) Sloan, L. C. and Pollard, D. (1998).  Polar stratospheric clouds: A high latitude warming mechanism in an ancient greenhouse world.  Geophysical Research Letters, 25(18):3517–3520.
  * Stevens et al., (1998) Stevens, B., Cotton, W. R., Feingold, G., and Moeng, C.-H. (1998).  Large-eddy simulations of strongly precipitating, shallow, stratocumulus-topped boundary layers.  Journal of the atmospheric sciences, 55(24):3616–3638.
  * Turton and Nicholls, (1987) Turton, J. and Nicholls, S. (1987).  A study of the diurnal variation of stratocumulus using a multiple mixed layer model.  Quarterly Journal of the Royal Meteorological Society, 113(477):969–1009.
  * Twomey, (1977) Twomey, S. (1977).  The influence of pollution on the shortwave albedo of clouds.  Journal of the atmospheric sciences, 34(7):1149–1152.
  * Uchida et al., (2010) Uchida, J., Bretherton, C., and Blossey, P. (2010).  The sensitivity of stratocumulus-capped mixed layers to cloud droplet concentration: Do les and mixed-layer models agree?  Atmospheric Chemistry and Physics, 10(9):4097–4109.
  * Wang and Feingold, (2009) Wang, H. and Feingold, G. (2009).  Modeling mesoscale cellular structures and drizzle in marine stratocumulus. part i: Impact of drizzle on the formation and evolution of open cells.  Journal of the Atmospheric Sciences, 66(11):3237–3256.
  * Warren et al., (1988) Warren, S., Hahn, C., London, J., Chervin, R., and Jenne, R. (1988).  Global distribution of total cloud cover and cloud types over ocean. ncar tech.  Technical report, Note NCAR/TN-317+ STR, 42 pp.+ 170 maps.
  * Wood, (2012) Wood, R. (2012).  Stratocumulus clouds.  Monthly Weather Review, 140(8):2373–2423.