Stratocumulus clouds cover 20% of the low-latitude oceans and are especially prevalent in the subtropics. They cool the Earth by shading large portions of its surface from sunlight. However, as their dynamical scales are too small to be resolvable in global climate models, predictions of their response to greenhouse warming have remained uncertain. Here we report how stratocumulus decks respond to greenhouse warming in large-eddy simulations that explicitly resolve cloud dynamics in a representative subtropical region. In the simulations, stratocumulus decks become unstable and break up into scattered clouds when CO₂ levels rise above 1,200 ppm. In addition to the warming from rising CO₂ levels, this instability triggers a surface warming of about 8 K globally and 10 K in the subtropics. Once the stratocumulus decks have broken up, they only re-form once CO₂ concentrations drop substantially below the level at which the instability first occurred. Climate transitions that arise from this instability may have contributed importantly to hothouse climates and abrupt climate changes in the geological past. Such transitions to a much warmer climate may also occur in the future if CO₂ levels continue to rise.

The stratocumulus clouds that top turbulent boundary layers over large swaths of subtropical oceans are unusual among Earth’s clouds: the turbulence that sustains them is driven by longwave radiative cooling of the cloud tops, rather than by heating at the surface¹. The reason is that cloud droplets absorb and emit longwave radiation so efficiently that the clouds are essentially opaque to longwave radiation. By contrast, the free troposphere above them is dry and has few clouds, making it relatively transparent to longwave radiation. Therefore, the longwave flux downwelling just above the cloud tops originates in higher and cooler atmospheric layers and is 50–90 W m⁻² weaker than the upwelling longwave flux emitted by the cloud tops (Fig. 1)¹. The resulting longwave cooling of the cloud tops drives convection, which penetrates the boundary layer and supplies the clouds with moisture from the surface. The convective moistening and radiative cooling of the cloud layer are balanced by turbulent entrainment of dry and warm free-tropospheric air across the trade inversion—the usually sharp interface between the boundary layer and the above-lying free troposphere (Fig. 1). Stratocumulus decks break up when the longwave cooling of the cloud tops becomes too weak to propel air parcels to the surface or when the turbulent entrainment of dry and warm free-tropospheric air across the inversion becomes too strong².

In stratocumulus clouds, longwave radiative cooling of the cloud tops propels air parcels downward, which convectively connects the clouds to their moisture supply at the surface. Turbulence entrains warm and dry air across the inversion, which counteracts the radiative cooling and convective moistening of the cloud layer. When the atmospheric concentration of greenhouse gases (for example, CO₂ and H₂O) increases (1,200 ppm), the longwave cooling of the cloud tops weakens, because the downwelling longwave radiation that reaches the cloud tops from above emanates at lower levels with higher temperatures relative to the cloud-top temperatures. Eventually, at sufficiently high greenhouse gas concentrations (1,300 ppm in our simulation without subsidence changes), stratocumulus decks break up into cumulus clouds, which leads to dramatic surface warming. Evaporation then strengthens, and the average longwave cooling at the level of the cloud tops drops to less than 10% of its value in the presence of stratocumulus decks.

Full size image

Key processes in stratocumulus-topped boundary layers, such as entrainment across the inversion, occur at scales of tens of metres and smaller^3,4,5,6. This is much too small to be resolvable in global climate models (GCMs), which currently have horizontal resolutions of tens of kilometres⁷. GCMs therefore resort to parameterizations, which relate stratocumulus occurrence to resolved large-scale variables, such as temperature or humidity. However, the parameterizations are notoriously inaccurate: GCMs severely underestimate the prevalence of stratocumulus decks^8,9, and confidence in the simulations of their climate change response is low¹⁰. As stratocumulus clouds cover 20% of the tropical oceans¹¹ and critically affect the Earth’s energy balance (they reflect 30–60% of the shortwave radiation incident on them back to space¹), problems simulating their climate change response percolate into the global climate response. For example, uncertainties in the response of stratocumulus and other low clouds lead to large uncertainties in the predictions of global temperatures and climate impacts^12,13,14.

Here we exploit advances in high-performance computing and large-eddy simulation (LES) of clouds, which enable us to faithfully simulate statistically steady states of stratocumulus-topped boundary layers in restricted regions^{15,16,17,18,19,20,21}. Whereas GCMs resolve the most energetic large-scale dynamics and parameterize cloud-scale dynamics, our LES resolves the most energetic cloud-scale dynamics and parameterizes large-scale dynamics (Methods). The LES explicitly simulates cloud-scale dynamics over an ocean patch, with large-scale conditions representative of summertime in subtropical regions that currently have persistent stratocumulus decks, such as off the coasts of California or Peru (Fig. 1)²². Whereas previous LES studies prescribed the sea surface temperature (SST) and thereby precluded feedbacks between SST and cloud cover^{4,15,16,17,18}, the SST in our subtropical LES domain is controlled by the surface energy balance; hence, it varies and interacts with cloud cover changes^23,24. It equilibrates to a realistic 290 K in a baseline simulation with 400 ppm atmospheric CO₂—approximately today’s level.

The subtropical LES domain is coupled to a tropical atmospheric column in radiative–convective equilibrium (RCE), as in classical conceptual models of the tropical atmosphere²⁵. The SST in the tropical column is controlled by the imbalance R_t between net incoming and outgoing radiation at the top of the atmosphere (TOA). For the baseline simulation, a tropical TOA imbalance R_t within the range of observed values yields a realistic tropical SST of 300 K. In the subtropical LES domain, however, the TOA imbalance R_s varies with atmospheric composition and cloud cover. It increases by ΔR_s relative to the baseline simulation if the outgoing longwave radiation weakens, or if more incoming shortwave radiation is absorbed because reflective cloud cover decreases. Such an increased R_s in nature needs to be balanced by an enhanced energy export out of the subtropics. In our simulations, we assume that the implied energy export is distributed homogeneously across the globe, so that our tropical TOA imbalance R_t decreases when R_s increases, by an amount that depends on ΔR_s and the area fraction γ of the globe represented by the subtropical LES domain. As stratocumulus clouds cover 18.5% of the oceans between 5° and 35° latitude in both hemispheres, and this subtropical ocean area makes up 35% of Earth’s surface area¹¹, we estimate γ = 0.185 × 0.35 = 6.5%. Two-way coupling between the subtropical LES domain and the tropical column results by driving the free-tropospheric temperature in the LES domain towards the moist adiabatic temperature profile of the tropical column, as is common in two-column models of the tropical atmosphere^24,25,26.

The baseline simulation with 400 ppm CO₂ produces stratocumulus decks like those observed (Fig. 2). When CO₂ levels are increased, SST first increases as in current GCMs²⁷, for example, in the tropics by 3.6 K for the CO₂ doubling from 400 to 800 ppm (Fig. 3d). (CO₂ levels here should be understood as the equivalent CO₂ levels that correspond to a change in the concentration of all well-mixed greenhouse gases.) Cloud cover remains dense (Fig. 3a), but the amount of liquid water in the clouds decreases slightly (Fig. 3b), as seen in previous LES studies with prescribed SSTs^17,19. However, when a CO₂ threshold is crossed at around 1,200 ppm, the stratocumulus decks abruptly become unstable and break up into scattered cumulus clouds (Figs. 2 and 3a,b, Supplementary Fig. 1 and Supplementary Movie 1). When CO₂ levels are lowered again after the stratocumulus breakup, the stratocumulus decks only reform once the CO₂ levels drop below 300 ppm (Fig. 3a,b). That is, there is bistability as a function of CO₂ levels, and this results in hysteresis.

Stratocumulus decks prevail in the baseline simulation at 400 ppm CO₂, which has a subtropical SST of 290 K. At 1,600 ppm CO₂, the stratocumulus decks have been replaced by scattered cumulus clouds, which leads to strong warming and a subtropical SST of 308 K because cloud shading of the surface is reduced.

Full size image

a, Subtropical cloud fraction. b, Cloud liquid water path (LWP). c, Subtropical SST. d, Tropical SST. Red upward arrows indicate the simulations that started from the baseline simulation with 400 ppm CO₂; blue downward arrows indicate the simulations that started from 1,600 ppm. The CO₂ axis is logarithmic (ticks every 200 ppm) because the radiative forcing of CO₂ is logarithmic in concentration. Departures of the SST changes from straight lines indicate the degree to which the climate sensitivity is state dependent. Higher climate sensitivities at high CO₂, after stratocumulus breakup, occur because the tropical column is approaching a local runaway greenhouse state²⁵.

Full size image

The subtropical SST jumps by 10 K and the tropical SST by 8 K across the stratocumulus instability (Fig. 3c,d). The tropical warming is a plausible estimate of the global-mean warming triggered by the instability. Subtropical marine stratocumulus clouds cover about 6.5% of the Earth’s surface and, where they occur, reduce the solar radiative energy flux absorbed in the climate system by ~110 W m⁻², compared to about a 10 W m⁻² reduction by scattered cumulus^22,28. If we assume a climate sensitivity parameter of 1.2 K (W m⁻²)⁻¹ (as for the more sensitive among current GCMs²⁷), this implies (110 − 10) W m⁻² × 6.5% × 1.2 K (W m⁻²)⁻¹ ≈ 8 K global-mean surface warming when subtropical marine stratocumulus break up.

Two key mechanisms and their interaction with the surface shape the stratocumulus instability and hysteresis. First, as the atmosphere above the clouds becomes more opaque to longwave radiation when CO₂ levels rise, the downwelling longwave flux at the cloud tops increasingly originates in lower and warmer atmospheric layers. Hence, the difference between upwelling and downwelling radiative energy fluxes at the cloud tops decreases, and the cloud-top long-wave cooling weakens (Fig. 1 and Supplementary Fig. 2a). This weakens the convection between the cloud layer and the surface and promotes the breakup of stratocumulus decks through decoupling from their surface moisture supply^2,17,29. Second, evaporation at the surface strengthens under warming (Fig. 1 and Supplementary Fig. 2b), which, other things being equal, enhances the generation of turbulence through latent heat release in the cloud layer. This strengthens the turbulent entrainment of dry and warm air across the inversion and likewise promotes stratocumulus breakup²—unless it is, as is the case in our and other simulations¹⁹, compensated by other processes, such as thinning of the cloud layer over which the turbulence generation can be realized, or weakening of the entrainment by a reduction in the cloud-top longwave cooling. A minimal conceptual model² that combines these interacting processes suggests that stratocumulus decks break up when the instability parameter S = (LHF/ΔL) × (h_c/h) exceeds a critical value around S_c ≈ 0.6, where LHF is the latent heat flux at the surface, ΔL is the longwave cooling of the cloud tops, h_c is the thickness of the cloud layer and h is the cloud-top height (Fig. 1). In our simulations, the instability parameter S increases from 0.4 at 400 ppm to 0.7 at 1,200 ppm. The increase in S arises because ΔL decreases by 17% from 400 ppm to 1,200 ppm, LHF increases by 13% and h_c/h increases by 24% (Supplementary Fig. 2). (However, the cloud thickness h_c itself decreases by 5%.)

Amplifying cloud cover–SST feedbacks are crucial for the abrupt stratocumulus breakup: as the stratocumulus cloud cover decreases, the SST increases, which strengthens surface evaporation and enhances the atmosphere’s longwave opacity through water vapour feedback. Indeed, if the water vapour concentrations that the radiative transfer scheme sees are fixed at their values in the baseline simulation, the reduction in cloud-top radiative cooling at 1,200 ppm relative to the baseline is 55% weaker than with water vapour feedback; that is, water vapour feedback accounts for about half of the reduction in cloud-top cooling. Both the effect of water vapour feedback on cloud-top radiative cooling and the strengthened surface evaporation lead to a sharp increase in S across the instability (for example, LHF alone jumps 73%, see Fig. 1 and Supplementary Fig. 2). The non-linear changes in the thermodynamic state of the atmosphere and surface cannot immediately reverse when CO₂ levels drop again (Supplementary Fig. 2). This leads to the bistability and hysteresis. Since previous LES studies have prescribed the SST and thus suppressed the surface feedbacks, the abrupt stratocumulus instability, as well as the bistability and hysteresis as a function of CO₂ levels, remained undiscovered, although the governing cloud-layer mechanisms were known^2,17,18,19. (Multiple equilibria of stratocumulus clouds have been demonstrated in a mixed-layer model as a function of large-scale subsidence³⁰ and in LES as a function of the initial condition³¹. However, these involve mechanisms and phenomenology that differ from those here, where the focus is on bistability as a function of CO₂ levels.)

The CO₂ level at which the instability occurs depends on how large-scale dynamics change with climate, which is heuristically parameterized in our simulations and hence is uncertain. In particular, the large-scale subsidence in the troposphere weakens under warming³², which lifts the cloud tops and counteracts the instability^15,19,24. Indeed, when we weaken the parameterized large-scale subsidence by 1 or 3% per Kelvin of tropical SST increase (within the range of GCM responses to warming³³), the stratocumulus instability occurs at higher CO₂ levels: around 1,400 ppm with 1% K^–1 subsidence weakening, and around 2,200 ppm with 3% K^–1 (Fig. 4). The hysteresis when the CO₂ levels drop thereafter remains, but it narrows: stratocumulus decks reform once the CO₂ levels drop below 500 ppm for a 1% K^–1 subsidence weakening, and once they drop below 1,900 ppm for one of 3% K^–1.

a, Subtropical cloud fraction. b, Cloud liquid water path. c, Subtropical SST. d, Tropical SST. The left column shows simulations in which subsidence weakens at a rate of 1% K^–1 with tropical SST increase. The right column shows simulations in which subsidence weakens at a rate of 3% K^–1. Red upward arrows indicate simulations that started from 400 ppm CO₂; blue downward arrows indicate simulations that started from 1,600 ppm (left column) and 2,400 ppm (right column). As in Fig. 3, the CO₂ axis is logarithmic.

Full size image

We expect the width of the hysteresis loop in nature to be reduced further by sources of noise neglected in our simulations, such as seasonal or synoptic variations. Nonetheless, although the precise extent to which stratocumulus decks can be bistable remains to be investigated, the abrupt instability itself appears to be robust for the physical reasons outlined above and as seen in the simulations here. The instability will probably occur first in regions and seasons in which the stratocumulus decks are close to the stability threshold, for example, at the margins of current stratocumulus regions. One may expect large transient fluctuations in cloud cover between the states with and without stratocumulus decks near the stability threshold—the flickering phenomenon common near critical transitions in complex dynamical systems³⁴. If the stratocumulus decks in different subtropical regions differ in their proximity to the stability threshold, the global effects of the instability as a function of CO₂ levels may also be smoothed out.

Current GCMs run with CO₂ levels up to 9,000 ppm appear not to exhibit the stratocumulus instability^27,35,36, probably because interactions of turbulence, cloud processes, and radiation are inadequately parameterized in them¹⁹. Inadequacies in cloud parameterizations may account for the difficulties GCMs have in simulating warm climates of the past³⁷. For example, the Arctic was frost-free during the early Eocene, around 50 million years ago. However, current GCMs only reproduce a frost-free Arctic at CO₂ levels above 4,000 ppm^35,38—much higher than the levels below 2,000 ppm reconstructed for the Eocene³⁹. The stratocumulus instability may explain how such hothouse climates can exist without implausibly high CO₂ levels. It may also have played a role in past climate transitions, such as the Eocene–Oligocene transition to a cool climate 34 million years ago⁴⁰.

For the future, our results suggest that stratocumulus decks may break up if CO₂ levels continue to rise. Equivalent CO₂ concentrations around 1,300 ppm—the lowest level at which the stratocumulus instability occurred in our simulations—can be reached within a century under high-emission scenarios⁴¹. However, it remains uncertain at which CO₂ level the stratocumulus instability occurs because we had to parameterize rather than resolve the large-scale dynamics that interact with cloud cover. To be able to quantify more precisely at which CO₂ level the stratocumulus instability occurs, how it interacts with large-scale dynamics and what its global effects are, it is imperative to improve the parameterizations of clouds and turbulence in climate models.

Experimental design

We performed the numerical experiments using the Python Cloud Large Eddy Simulation code, which solves the anelastic equations of motion with total water specific humidity q_t and specific entropy s as the prognostic thermodynamic variables²⁰. The lower boundary of the subtropical LES domain is a slab ocean whose surface temperature T_s,0 evolves according to the surface energy balance^23,24. The subtropical LES domain is coupled to a tropical column in RCE, with weak horizontal temperature gradients in the free troposphere between the tropical column and the subtropical LES domain^25,26. We incorporated simple but physically plausible formulations for the energy and moisture exports from the subtropical LES domain into other regions. Our simulation set-up builds on Tan et al.^23,24 and is described in what follows.

Simulation numerics and resolution

The simulations employ an implicit LES approach, with nominally fifth-order weighted essentially non-oscillatory (WENO) advection schemes⁴² for all the prognostic variables. Subgrid-scale turbulent fluxes are modelled implicitly by the numerical dissipation inherent to WENO schemes, except in the surface layer, where a Smagorinsky–Lilly subgrid-scale closure allows surface fluxes to penetrate into the flow domain. We previously found this approach to lead to high-fidelity stratocumulus simulations because it minimizes spurious numerical mixing near the inversion²¹. Indeed, in benchmark tests⁴³ our LES approach at a relatively coarse resolution (5–10 m in the vertical) reproduces observations as well as LES and direct numerical simulations at a higher resolution^44,45 and greater computational expense (Supplementary Fig. 3).

The model equations are integrated forward in time using an explicit second-order strong stability-preserving Runge–Kutta scheme⁴⁶, with a timestep dynamically adjusted to maintain a Courant–Friedrichs–Lewy (CFL) number of approximately 0.7. The model is integrated until a statistically steady state is reached, which takes 40–350 simulated days (Supplementary Fig. 1). To reduce the large computational expense of such long simulations, a mean-state acceleration is applied to the horizontal velocities u and v, scalars q_t and s, and T_s,0: the horizontally averaged tendencies of these variables are boosted by a factor 8 to approach a statistically steady state more quickly⁴⁷. (As the surface temperature T_s,0 is the same throughout the subtropical LES domain, for T_s,0 this amounts to reducing the heat capacity of the slab ocean.) In tests of a few cases, we found the same statistically steady states without the mean-field acceleration, so we expect our results to be insensitive to the mean-field boosting. The LES domain extends 4.8 km in each horizontal direction and 2.25 km in the vertical. The horizontal and vertical grid spacings are 50 m and 10 m, respectively, for a total of 2 million grid points. We conducted additional simulations at a coarser resolution (75 m × 15 m), with essentially unchanged results (Supplementary Fig. 4). Therefore, although our LES resolution is not sufficient to have reached numerical convergence, we are confident in the numerical robustness of the results.

Radiative transfer

Radiative energy fluxes are calculated with the Rapid Radiative Transfer Model for GCMs (RRTMG)⁴⁸, with radiative temperature tendencies updated every 20 s of simulated time. The concentrations of all the atmospheric trace gases except CO₂ and water vapour are specified as the model’s default profiles (for example, the ozone concentration increases from 30 ppb near the surface to about 9 ppm above 10 hPa; the methane concentration decreases from 1.7 ppm near the surface to 0.15 ppm at the TOA). CO₂ is assumed to be well mixed. As the subtropical LES domain is limited to the lower troposphere, the profiles of temperature and specific humidity required for the radiative transfer calculations are extended above the LES domain following the tropical moist adiabat and assuming a fixed relative humidity of 30%. The atmosphere is assumed to be cloud free above the LES domain.

Insolation

TOA insolation is specified as the diurnally averaged insolation in July. For a total solar irradiance of 1,365 W m⁻², this gives an effective solar zenith angle of 69.9° and an incoming shortwave flux of \(S_{\rm{s}}^ \downarrow = 471\,{\mathrm{W}}\,{\mathrm{m}}^{ - {\mathrm{2}}}\) into the subtropical LES domain (assumed to be located at 30 °N). It gives an effective solar zenith angle of 73.4° and an incoming shortwave flux of \(S_{\rm{t}}^ \downarrow = 390\,{\mathrm{W}}\,{\mathrm{m}}^{ - {\mathrm{2}}}\) into the tropical column (assumed to be located at the Equator). Since shortwave and longwave cloud radiative effects nearly cancel in the deep tropics, the RCE calculation in the tropical column does not explicitly consider cloud radiative effects²⁵. However, we reduced the incoming shortwave flux at the TOA by a factor of (1 − α_t), where α_t = 0.09 is an effective albedo owing to tropical clouds, to obtain a realistic tropical SST. This results in an effective incoming shortwave flux of 356 W m⁻² into the tropical column.

Surface energy balance

The surface temperature T_s,0 in the subtropical LES domain evolves according to the surface energy balance^23,24:

where ρ_w = 10³ kg m⁻³ is the density of liquid water, C_w = 4.19 × 10³ J kg⁻¹ K⁻¹ is the specific heat capacity of water and H_w = 1 m is the thickness of the slab ocean. (This slab ocean is thinner than typical ocean mixed layers in the subtropics. The thinness of the slab ocean affects the adjustment timescale of ocean temperatures T_s,0 and thus the timescale over which the stratocumulus decks break up in the simulations; however, we do not expect it to affect the statistically steady states on which we focus.) The quantities on the right-hand side are horizontal averages over the LES domain, denoted by the overbar: \({\mathrm{RAD}} = \left( {1 - \alpha _0} \right)S_{{\rm{s}},0}^ \downarrow + L_{{\rm{s}},0}^ \downarrow - L_{{\rm{s}},0}^ \uparrow\) is the net downward radiative energy flux into the surface, which consists of the downwelling shortwave flux \(S_{{\rm{s}},0}^ \downarrow\) and the down- and upwelling longwave fluxes \(L_{{\rm{s}},0}^ \downarrow\) and \(L_{{\rm{s}},0}^ \uparrow\), with surface albedo α₀ = 0.1; SHF and LHF are the sensible and latent heat fluxes into the atmosphere, which are determined at each point according to Monin–Obukhov similarity theory^20,49, assuming a surface roughness height of 10⁻³ m; OHU is the ocean heat uptake in the subtropics. We view OHU as an adjustable parameter to obtain a realistic SST in the baseline simulation; with \(\overline {{\mathrm{OHU}}} = - 3\,{\mathrm{W}}\,{\mathrm{m}}^{ - {\mathrm{2}}}\), we obtain an SST of 290 K—close to those observed in subtropical stratocumulus regions. We then kept OHU fixed in all the other simulations, which means that we assume the corresponding ocean energy flux convergence and/or seasonal heat storage to be fixed—an assumption that may need to be revisited in a more realistic set-up.

TOA energy balance and tropical–subtropical coupling

The subtropical LES domain is coupled to the tropical column through the free-tropospheric temperature profile and the TOA energy balance. In the subtropical LES domain, the TOA imbalance \(R_{\rm{s}} = S_{\rm{s}}^ \downarrow - S_{\rm{s}}^ \uparrow - L_{\rm{s}}^ \uparrow\) between net incoming shortwave radiation \(S_{\rm{s}}^ \downarrow - S_{\rm{s}}^ \uparrow\) and outgoing longwave radiation \(L_{{{\rm{s}}}}^ \uparrow\) is determined by the surface energy balance and the explicitly simulated dynamics in the domain, which includes cloud radiative effects. In the baseline simulation, this results in R_s(400 ppm) = −10 W m⁻². In the tropical column, the radiative imbalance R_t and the composition of the atmosphere determine the SST. With the relative humidity in the RCE calculation fixed as described below, a TOA imbalance of R_t(400 ppm) = 31 W m⁻² gives a realistic SST of 300 K in the baseline simulation. This imbalance is similar to that observed in the equatorial region (about 30–100 W m⁻² in the zonal mean, depending on the season and the precise latitude considered²⁸). The climate sensitivity increases as the TOA imbalance R_t decreases and the tropical column moves closer to a local runaway greenhouse state²⁵.

When the TOA imbalance in the subtropical LES domain increases relative to the baseline simulation by ΔR_s = R_s − R_s(400 ppm), we change the TOA imbalance in the tropical column proportionately by ΔR_t = R_t − R_t(400 ppm) = −γ(1 − γ)⁻¹ΔR_s. What motivates this choice is the assumption that the anomalous subtropical imbalance is balanced by the divergence of an anomalous atmosphere–ocean energy flux F, so that ΔR_s = ∇ · F. The divergence of this anomalous energy flux F must integrate to zero over the globe. Our formulation for how ΔR_t changes results from assuming that the subtropical LES domain is representative of a fraction γ of the globe, and that everywhere outside this area (that is, in the complementary fraction 1 − γ of the globe) there is a uniform convergence of the anomalous atmosphere–ocean energy flux.

RCE calculation

The RCE temperature profile T_t(z) in the tropical column is computed using a time marching approach⁵⁰. Convective adjustment is performed with a direct algorithm⁵¹, using the moist-pseudoadiabatic temperature profile of a parcel lifted from the surface with an 80% relative humidity. The water vapour specific humidity profile q_v,t(z) is set to maintain a 60% relative humidity below the height of the temperature minimum H in the upper troposphere/lower stratosphere; the specific humidity is held constant above this cold point, so that q_v,t(z > H) = q_v,t(H). This procedure is iterated until all the temperature tendencies in the column fall below 0.001 K d⁻¹ and the computed TOA radiative flux imbalance R_t is within 0.01 W m⁻² of its current desired value. The TOA imbalance itself is linearly relaxed to its target value R_t(400 ppm) + ΔR_t on a timescale of ten simulated days to avoid rapid fluctuations that otherwise may arise because ΔR_t evolves with the dynamics of the subtropical LES domain. The RCE calculation is updated whenever the difference between the current desired value and the last computed value of TOA radiative flux imbalance exceeds 0.2 W m⁻².

The equilibrium climate sensitivity in RCE with these parameters and a fixed TOA imbalance of R_t = R_t(400 ppm) = 31 W m⁻² is 3.9 K. This is based on a quadrupling of CO₂ levels from 400 to 1,600 ppm, which yields a 7.9 K warming. Doubling the CO₂ levels from 400 to 800 ppm yields a weaker warming of 3.5 K, whereas the next doubling to 1,600 ppm yields a stronger warming of 4.4 K. (The climate sensitivity increases with temperature at fixed relative humidity because of strengthening water vapour feedback⁵².) The reduction in R_t when the shortwave reflection by subtropical low clouds diminishes leads to additional climate sensitivity increases with warming.

Free troposphere in the subtropical LES domain

In the free troposphere of the subtropical LES domain, horizontal-mean profiles of variables ϕ(z) ∈ {u_s(z), v_s(z), q_v,s(z) and T_s(z)} are relaxed towards the reference profiles \(\bar \phi _{{\mathrm{ref}}}(z)\) (the overbar indicates the horizontal mean)^15,16,23,24. Relaxation is applied above a height of 1.2h_r, where h_r is determined as the lowest model level at which the specific humidity \(\bar q_{{\rm{v}}, {\rm{s}}}\) falls below 1.2q_v,ref at the same level^15,23. The relaxation is linear, with horizontally uniform tendencies:

and a relaxation coefficient:

Here, τ = 6 h is a relaxation timescale, and η is a height-dependent factor given by:

For z ≥ 1.5h_r, this gives a constant relaxation coefficient ξ ≈ (26.1 h)⁻¹. (A coding error led to the omission of a factor π in the cosine in equation (3), which would have made the transition in relaxation coefficient at z = 1.5h_r smoother¹⁵.)

The reference profiles are specified as follows:

1. (1)

   Horizontal velocities: the reference velocities are:

   $$u_{{\mathrm{ref}}}(z) = {\mathrm{min}}\left[ { - 4,2 - 1.2p_0(z){\mathrm{/}}(10^4\,{\mathrm{Pa}})} \right]{\mathrm{m}}\,{\mathrm{s}}^{ - {\mathrm{1}}}$$

   (5)
   $$v_{{\mathrm{ref}}}(z) = 0$$

   (6)

   where p₀ is the anelastic reference pressure.

2. (2)

   Temperature: the reference temperature is the moist-adiabatic tropical temperature, T_ref(z) = T_t(z).

3. (3)

   Specific humidity: the reference specific humidity q_v,ref(z) corresponds to a relative humidity of 30% given the reference temperature T_ref(z). As in the tropical column, we set q_v,ref(z > H) = q_v,ref(H) above the cold point at height H.

Geostrophic winds

Horizontal winds in the subtropical LES domain are forced by Coriolis accelerations at 30° N. The geostrophic velocities are the reference velocities given by equations (5) and (6).

Energy and moisture export out of subtropical boundary layer

Large-scale dynamics (synoptic eddies and the Hadley circulation) transport energy and moisture out of the subtropical boundary layer into other regions. As a simple but physically plausible representation of these large-scale fluxes^23,24, we specify temperature and (total water) specific humidity tendencies at pressures above 900 hPa:

Here, \(q_{{{\rm{v}}}}^ \ast \left( {T_{{{{\rm{s}}}},0}} \right)\) is the saturation specific humidity, evaluated at the subtropical surface temperature T_s,0. These advective tendencies are constant at pressures above 900 hPa and go linearly in pressure to zero between 900 and 800 hPa. The temperature tendency corresponds to a fixed advective cooling of atmospheric columns of 21 W m⁻². The moisture tendency corresponds to a latent energy flux of 27 W m⁻² at 290 K, the subtropical temperature in the baseline simulation. The moisture tendency increases with the saturation specific humidity at the surface because large-scale moisture fluxes scale approximately with specific humidity, and relative humidity changes are comparatively weak under global warming^32,53,54,55. Additional energy and moisture fluxes in the subtropical free troposphere are implicit in the relaxation toward the reference profiles.

Subsidence

Tendencies due to large-scale subsidence are applied to the specific entropy s and specific humidity q_t. The large-scale subsidence velocity is specified as^23,24:

where p_* is the surface pressure, ρ₀(z) is the anelastic reference density, g is the gravitational acceleration and D is a divergence rate near the surface. For the standard simulations (Fig. 3), we chose D = 6 × 10⁻⁶ s⁻¹. The tendencies owing to large-scale subsidence are then given by:

We do not specify a large-scale subsidence that balances radiative cooling, as is sometimes done^26,56. The reason is that in, addition to radiative cooling, dynamic cooling associated with the divergence of eddy energy fluxes plays an essential role in the energy budget of the subtropical free troposphere⁵⁷. Assuming subsidence to be balanced by radiative cooling gives a too weak lower-tropospheric subsidence velocity that, moreover, weakens too rapidly as the climate warms (it would weaken by 20–40% per CO₂ doubling in our simulations).

Microphysics

As the clouds simulated in this study only occur at temperatures above freezing, a one-moment warm-rain bulk microphysics scheme suffices to represent the formation and sedimentation of precipitation in the clouds. We use a slightly modified version of the scheme in Kaul et al.⁵⁸, omitting ice-phase interactions and using more accurate, temperature-dependent formulae for vapour diffusivity and thermal conductivity of the liquid phase⁵⁹. The treatment of microphysics does not include the effects of cloud droplet sedimentation. The accounting of entropy and moisture source terms due to microphysical processes follows Pressel et al.²⁰.

Initial conditions

The baseline simulation with an atmospheric CO₂ concentration of 400 ppm is initialized from a well-mixed layer with a specific humidity of q_t = 8.449 × 10⁻³ and liquid potential temperature of θ_l = 290.95 K at pressures greater than 920 hPa. Liquid potential temperature fluctuations drawn from a uniform distribution between ±0.1 K are superimposed within the well-mixed layer to break the symmetry. Above the well-mixed layer, the initial profiles are the reference profiles with R_t = 31 W m⁻² and T_t,0 = 300 K. The SST in the subtropical LES domain is initialized to T_s,0 = 290 K and is held fixed for the first three simulated days, after which it evolves freely according to the surface energy balance. All other simulations with varying CO₂ concentrations (both in the tropical column and in the subtropical LES domain) are initialized from nearly statistically steady states of either the baseline simulation or a warm simulation after stratocumulus breakup (Supplementary Fig. 1). The warm simulation is the 1,600 ppm simulation for cases with a fixed subsidence and with a 1% K^–1 weakened subsidence, and it is the 2,400 ppm simulation for the cases with 3% K^–1 weakened subsidence.

Experiments with weakened subsidence

To examine the sensitivity of the results to the weakening of the large-scale subsidence expected under global warming^32,54, we carried out two series of simulations: one with the divergence D reduced by 1% K^–1 of the tropical SST increase relative to the baseline simulation, and one with D reduced by 3% K^–1. (D is correspondingly increased when the tropical SST falls below 300 K.) Such subsidence weakening is broadly representative of that seen in GCMs³³ and can be accounted for by theoretical arguments^32,54, keeping in mind that the changes in the gross vertical mass flux in the tropics are generally larger than the changes in net subsidence in the subsiding branches of tropical overturning circulations⁵⁴.

Sensitivity to experimental design

We experimented with variants of this set-up, which yielded qualitatively similar results, with quantitative differences in the CO₂ level of the stratocumulus instability. A broader exploration of the configuration space would be desirable but it is complicated by the large computational expense (the simulations reported here required about 2 million core hours of computing). Variants of the experimental setup with which we experimented include:

• Setups in which the SST difference between the tropical column and the subtropical LES domain was either fixed or gradually reduced as the shortwave cloud radiative effect in the subtropics weakened. In this setup, the stratocumulus instability tended to occur at lower CO₂ levels (as low as 800 ppm). However, the TOA energy balance in the tropical column (free to adjust in this setup) implied unrealistic lateral energy fluxes.

• Setups in which OHU increased as the subtropical clouds thinned and eventually broke up. This yielded qualitatively similar results to those reported here.

• A few simulations with a diurnal cycle (without boosting of the mean-field tendencies, which is difficult to justify in the presence of a diurnal cycle). Again, the stratocumulus instability occurred. However, because the overall surface energy balance with a diurnal cycle is not the same as with diurnally averaged insolation (principally because stratocumulus cover is reduced during the day), OHU or other parameters need to be adjusted in such simulations to obtain a realistic SST. This makes one-on-one comparisons with the results reported here difficult.

The CO₂ level at which the stratocumulus instability occurs depends on how large-scale dynamics are parameterized, as is evident in Fig. 3. But the instability itself appears physically robust.

The source code for the simulations is available at climate-dynamics.org/software/#pycles.

The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Information files. The raw data in the figures are available from the corresponding author upon request.

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This research was supported by C. Trimble. The computations were performed on ETH Zurich’s Euler cluster and on Caltech’s High Performance Cluster, which is partially supported by a grant from the Gordon and Betty Moore Foundation. We thank M. Hell for assistance with the figures. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

Affiliations

1. California Institute of Technology, Pasadena, CA, USA

   • Tapio Schneider
   • , Colleen M. Kaul
   •  & Kyle G. Pressel

2. Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA

   • Tapio Schneider

Contributions

T.S. designed the study, analysed results, and wrote the paper. C.M.K. and K.G.P. implemented the study design numerically, conducted the LES, analysed and visualized the results, and contributed to the writing.

Competing interests

The authors declare no competing interests.

Corresponding author

Correspondence to Tapio Schneider.

1. Supplementary Information

   Supplementary Figures

2. Supplementary Movie

   Time evolution of the stratocumulus instability. The shading shows the liquid water path in the subtropical LES domain in the 1300 ppm simulation between days 255 and 275, which is during the time of the stratocumulus breakup (Fig. 1). Blue for zero LWP, and the colorscale saturates white at a LWP of 300 μm. The time series at the bottom show the cloud fraction (blue) and SST (red) in the LES domain. The breakup of the stratocumulus clouds is more rapid than it would be in nature because of the unrealistically small thermal inertia of the underlying slab ocean.

To obtain permission to re-use content from this article visit RightsLink.

Publication history

DOI

https://doi.org/10.1038/s41561-019-0310-1