Moderate climate sensitivity due to opposing mixed-phase cloud feedbacks

Introduction

How much will Earth ultimately warm in response to a doubling of atmospheric CO₂ concentrations? The answer to this fundamental question about the climate sensitivity of Earth has plagued climate science for decades. Early estimates of Earth’s climate sensitivity spanned 3 °C, ranging from 1.5 °C to 4.5 °C¹. Global climate models (GCMs) participating in the Fifth Phase of the Coupled Model Intercomparison Project (CMIP5) have since estimated a narrower range spanning from 2.0 °C to 4.6 °C². However, the range in climate sensitivity estimated by GCMs participating in the most recent Sixth Phase of the Coupled Model Intercomparison Project (CMIP6) has increased by 1.2 °C relative to CMIP5 and now spans from 1.8 °C to 5.6 °C³. Moreover, of the GCMs participating in CMIP6, 9 out of 27 of them estimate a climate sensitivity value that exceeds the upper range of 4.6 °C estimated by the GCMs participating in CMIP5. This “hot model problem”⁴ has since been reconciled with expert judgment based on data from the historical period and paleoclimate records resulting in a very likely (5–95%) range of 2.3 °C to 4.7 °C⁵. This range is comparable to but 0.6 °C narrower than the very likely range of 2 °C to 5 °C derived in the IPCC’s most recent Sixth Assessment Report⁶.

The representation of the surface temperature-mediated response of clouds in GCMs is a longstanding and dominant source of uncertainty in climate sensitivity^6,7. Many studies have since focused on identifying sources of errors in the representation of clouds in GCMs that impact climate sensitivity in various latitudinal regions^{8,9,10,11,12,13}. In the extratropical regions, transitions in the thermodynamic phase of hydrometeors within mixed-phase clouds in the ~−38 °C to 0 °C temperature range are a robust feature of global warming experiments performed by GCMs¹⁴. This cloud-phase feedback involves the replacement of ice with supercooled liquid within mixed-phase clouds due to the deepening of the troposphere in response to global warming¹⁵. Since liquid droplets are typically more abundant and smaller in size compared to ice crystals¹⁶, the transition to more liquid-dominated mixed-phase clouds in response to global warming reflects more sunlight back out to space¹⁷. The net result is a negative feedback that counteracts global warming, with pre-warmed climate states exhibiting lower supercooled liquid fractions within mixed-phase clouds leading to a more negative feedback due to the larger replacement of ice with liquid due to the cloud-phase feedback. Many of the GCMs participating in CMIP5 were found to underestimate the supercooled liquid fraction within mixed-phase clouds relative to satellite observations on a global scale^18,19. When the supercooled liquid fractions of these GCMs were composited in temperature bins over the Southern Ocean, the temperature at which half of the clouds glaciated—termed T50—exhibited a spread of nearly 40 °C, with the majority of GCMs producing excessively warm T50 values that indicate an underestimate in the supercooled liquid fraction within mixed-phase clouds²⁰.

Results

Evolution of cloud phase partitioning from CMIP5 to CMIP6

Several GCMs participating in CMIP6 have since improved their representation of mixed-phase clouds^21,22,23. To examine how T50 evolved from CMIP5 and CMIP6, we systematically calculate T50 for the full ensemble of 61 GCMs with available data from experiments representative of preindustrial conditions (piControl)^24,25. These calculations use cloud mixing ratios and local temperature in model grid cells. We focused on the extratropics (latitudes poleward of 60° in both hemispheres) as mixed-phase clouds are particularly ubiquitous in these regions^{26,27,28,29,30,31}. We also calculate T50 produced by experiments holding sea surface temperature and sea ice constant to 1978 − 2006 conditions (AMIP)³². We additionally calculated T50 using occurrences of cloud liquid and ice produced by a lidar simulator designed to be consistent with an active satellite product, GCM-Oriented CALIPSO Cloud Product (GOCCP)³³ derived from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument³⁴ within the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) framework³⁵.

Overall, the ensemble mean T50 of the GCMs participating in CMIP6 is −21 °C and is 4.2 °C colder than the ensemble mean T50 of —15.8 °C of the GCMs participating in CMIP5 (Fig. 1). Observations of T50 derived from GOCCP yield a T50 estimate of −19.7 °C in the extratropics. This observed T50 value is comparable to the multimodel ensemble average of −21 °C in CMIP6, however, arises from averaging 18 GCMs that overestimate the fraction of supercooled liquid water within mixed-phase clouds and 17 GCMs that underestimate this fraction. This contrasts with the 21 out of 26 GCMs that underestimate the supercooled liquid fraction of mixed-phase clouds in participating in CMIP5. The spread in T50 in both generations of GCMs is equally large, spanning ≈30 °C.

Moderate climate sensitivity due to opposing mixed-phase cloud feedbacks — **Fig. 1: Histogram quantifying the evolution of T50 from CMIP5 to CMIP6.**

The generally enhanced proportion of supercooled liquid water within mixed-phase clouds simulated by the most recent generation of GCMs participating in CMIP6 (Fig. 1) is consistent with previous literature^21,22,36,37; there are more GCMs with cold T50 values in CMIP6 than CMIP5, and several GCMs simulate T50 values close to the homogeneous freezing temperature of ≈−38 °C. This result is consistent across both piControl and AMIP experiments as most GCMs produce similar T50 values for both experiments (Fig. S1a). However, the COSP lidar simulator produces a dissimilar distribution of T50 values, with values that are ≥9 °C colder than those calculated using the raw cloud mass variables in half of the GCMs with data available (Fig. S1b). Indeed, a recent reexamination of the thresholds of layer-integrated depolarization ratio and layer-integrated attenuated backscatter of fully attenuating clouds that are incorporated in the CALIOP phase discrimination algorithm suggests a greater proportion of ice hydrometeors within mixed-phase clouds than originally reported³⁸. For example, equatorward of 50°S in the Southern Ocean where mixed-phase clouds are ubiquitous, a T50 value of ≈ −26 °C was derived from the vantage point of ship-based lidar observations while the T50 value based on spaceborne lidar observations shows substantially more liquid and is ≈−31 °C³⁹. Therefore, the lidar simulator potentially overestimates the fraction of supercooled liquid within mixed-phase clouds and we use the raw GCM cloud mass variables in our main analysis to maximize the number of GCMs in our analysis.

Cloud phase and cloud feedback relationships in GCMs

A negative linear relationship between T50 in the pre-warmed climate state and the shortwave (SW) cloud optical depth feedback (λ_τ) in response to global warming was demonstrated in an ensemble of perturbed cloud microphysics experiments in two GCMs^40,41. The linear relationship was particularly prominent in the extratropics and was linked to the cloud-phase feedback. Given the greater proportion of mixed-phase cloud liquid mass in the piControl simulations in CMIP6 (Fig. 1), we examine whether the relationship between T50 and the SW λ_τ previously simulated by the ensemble of perturbed microphysics experiments extends to the full ensemble of GCMs participating in CMIP5 and CMIP6.

A negative linear correlation between T50 and the SW λ_τ in the extratropics emerges in 61 GCMs participating in CMIP5 and CMIP6 (Fig. 2a). This negative correlation is consistent with the cloud-phase feedback: GCMs with warmer T50 values contain relatively higher proportions of ice within mixed-phase clouds that become replaced with optically thicker and thus more reflective liquid droplets in response to global warming that in turn induce a more negative λ_τ. The negative relationship is stronger in CMIP5 than in CMIP6 (Fig. S2a, S2b). However, at the same time, a competing positive linear relationship between T50 and changes in cloud amount in response to a quadrupling of CO₂, i.e. the SW cloud amount feedback (λ_amt) also emerges (Fig. 2b). This positive linear relationship is also stronger among the GCMs participating in CMIP5 (Fig. S2c) compared to the GCMs participating in CMIP6 (Fig. S2d) and has recently been attributed to a “supercooled cloud feedback”³⁷. In this feedback, the supercooled liquid within clouds shift to warm liquid water clouds at temperatures above freezing in response to global warming. Since supercooled liquid water-containing clouds precipitate more efficiently than warm liquid water-containing clouds^42,43,44, transitions from the former to the latter in response to global warming would result in increases in the coverage of reflective clouds and therefore a negative SW λ_amt. As such, T50 is positively correlated with λ_amt. This supercooled cloud feedback associated with λ_amt dominates the “cloud lifetime feedback”^45,46 that would otherwise result in a negative relationship between T50 and λ_amt. This latter cloud lifetime feedback is associated with the extended longevity of mixed-phase clouds that are composed of larger amounts of cloud liquid due to the lower precipitation efficiency of liquid droplets compared to ice crystals.

**Fig. 2: Relationship between extratropical T50 and SW extratropical cloud feedbacks in GCMs participating in CMIP5 and CMIP6.**

Adjustment of climate sensitivity via mixed-phase clouds

Given the linear T50-λ_τ and T50-λ_amt relationships (Fig. 2), we now exploit these relationships to adjust the climate sensitivity of each individual GCM based on its bias in T50 relative to GOCCP observations. Our method is based on the forcing-feedback framework⁴⁷, where the change in the net global mean radiative flux at the top of the atmosphere in response to radiative forcing, ΔN, can be expressed as

$$Delta N=F+lambda Delta {T}_{s},$$

(1)

where F is the global mean radiative forcing perturbation, λ is the climate radiative feedback parameter in response to F, and ΔT_s is the change in global mean surface air temperature in response to F. λ can be decomposed into the summation of surface and air temperature feedbacks, the albedo feedback, the water vapor feedback, and the cloud feedback, each of which can be further decomposed into the sum of their SW and longwave contributions⁴⁸. The cloud feedback can also be decomposed into a sum of contributions from cloud property changes including the SW λ_τ and λ_amt feedbacks that both exhibit a linear relationship with T50 (Fig. 2). However, T50 is a biased estimator of the observed T50 (Fig. 1). The bias of the extratropical T50 for every GCM can be calculated as

$$Delta T5{0}_{i}=T5{0}_{i}-T5{0}_{o},$$

(2)

where T50_o is the observed extratropical T50 derived from GOCCP and the subscript i denotes an individual GCM. For CMIP5 and CMIP6, we analyzed 26 and 35 GCMs, respectively. Through the linear T50-λ_τ and T50-λ_amt relationships established in the GCMs participating in CMIP5 and CMIP6, the biases in λ_τ and λ_amt can correspondingly be calculated as

$$Delta {lambda }_{c,i}=mcdot Delta T5{0}_{i},$$

(3)

where the subscript c represents either τ or amt, and m represents the slope of the least squares regression fit in Fig. 2. Substituting Eq. (3) back into Eq. (1), and solving for ΔT_s while assuming a vanishing value of ΔN at equilibrium yields the climate sensitivity of each individual GCM after adjusting for biases in T50,

$$EC{S}_{adj,i}=-frac{{F}_{i}}{{lambda }_{i}-wleft(Delta {lambda }_{tau ,i}+Delta {lambda }_{amt,i}right)}$$

(4)

where w ≈ 0.15 represents the areal fraction of the extratropics to that of the entire globe. F and λ were calculated as half the value of the y-intercept and the slope, respectively, of the least squares linear regression of the global annual mean change in the net radiative flux at the top of the atmosphere against the change in global annual mean surface temperature for 150 years of the piControl and abrupt4xCO2 experiments⁴⁷. Our method assumes that all other feedbacks in the climate system are held constant in the adjustment process. This contrasts with other methods that have observationally constrained T50 by way of modifying the representations of cloud physical processes within a given model that in turn interacted with other model processes^40,41,49.

The average climate sensitivity after adjusting for cloud feedback biases via their relationship with T50 based on 61 GCMs participating in both CMIP5 and CMIP6 is 3.63 K (Fig. 3). This value results from a downward adjustment of climate sensitivity from 3.28 K to 3.24 K based on the GCMs participating in CMIP5 (Fig. S3a) and essentially no change in the climate sensitivity of 3.86 K based on the GCMs participating in CMIP6 (Fig. S3b). In CMIP5, the downward adjustment in climate sensitivity results from upward adjustments via λ_τ (Fig. S3c) being overwhelmed by downward adjustments via λ_amt (Fig. S3e) while negligible adjustments in climate sensitivity occur due to λ_τ (Fig. S3e) and λ_amt in CMIP6 (Fig. S3f). Since the majority of the GCMs participating in CMIP5 underestimate the supercooled liquid fraction within mixed-phase clouds relative to the observations (Fig. 1) such that ΔT50 is positive for these GCMs, ECS_adj,i in Eq. (4) for these GCMs is therefore adjusted to larger values from the adjustment of λ_i to more positive values via λ_τ,i. The upward adjustment of the climate sensitivity of the GCMs participating in CMIP5 due to the λ_τ alone, i.e. assuming a vanishing value of λ_amt, is consistent with previous studies of the influence of T50 on the cloud-phase feedback^15,40. However, due to the steeper slope in the T50-λ_amt relationship of the opposite sign, ECS_adj,i of these same GCMs will be adjusted to small values.

**Fig. 3: Histogram of the climate sensitivity of GCMs and their climate sensitivity values after adjusting for biases in the partitioning of cloud thermodynamic phase in mixed-phase clouds in GCMs.**

In contrast, the lack of change in climate sensitivity in the GCMs participating in CMIP6 results from adjustments that are approximately canceled due to the equal number of GCMs that overestimate and underestimate the supercooled liquid fraction within mixed-phase clouds relative to observations. In contrast to the 26 GCMs participating in CMIP5, of the 35 GCMs considered (Fig. S4a), λ_τ was adjusted downward for approximately half of the GCMs while λ_τ was adjusted upward for the other half of the GCMs (Fig. S4b). The approximately equal number of GCMs that overestimate and underestimate T50 relative to observations results in a near-zero adjustment in the mean climate sensitivity from their original values via biases in λ_τ. Conversely, λ_amt was correspondingly adjusted upward for the GCMs that produced λ_τ values that were adjusted downward while λ_amt was correspondingly adjusted downward for the GCMs that produced λ_τ values that were adjusted upward. This similarly also resulted in a near-zero adjustment in the mean climate sensitivity from their original values via biases in λ_amt. We also note that despite the fact that we chose to analyze raw model cloud mass instead of COSP lidar simulator data to maximize the number of GCMs in our analysis, the fact that λ_τ and λ_amt have opposing effects on climate sensitivity ultimately leads to negligible material difference in our adjusted climate sensitivity estimate.

Discussion

The average adjusted climate sensitivity of 3.63 °C arising from corrections to biases in the thermodynamic phase partitioning within the mixed-phase clouds in GCMs participating in CMIP5 and CMIP6 suggests a moderate climate sensitivity that is within the very likely range based on two different assessments^5,6. Our estimate is also in line with recent estimates of 3.56 ±0.72 K and 3.47 ±0.33 K derived from observational constraints on low-level clouds equatorward of 60° in both hemispheres¹¹ and those confined to the subtropics¹², respectively. Our results thus provide additional support for a moderate climate sensitivity based on complementary observational constraints on extratropical mixed-phase clouds in a total of 61 GCMs. While seemingly contradictory to previous literature showing a rise in climate sensitivity in GCMs when mixed-phase clouds are constrained by observations, the rise in climate sensitivity due to the cloud-phase feedback is consistent; rather, it is the role of λ_amt that unexpectedly emerges in GCMs causing a decrease in the average climate sensitivity in CMIP5. Namely, the positive relationship between T50 and λ_amt hypothesized to be due to the supercooled cloud feedback³⁷ strongly opposes and dominates the negative relationship between T50 and λ_τ in CMIP5. This positive relationship between T50 and λ_amt albeit with a weaker magnitude was shown to exist in a single GCM (Fig. S2⁴⁰); however, is opposite in sign to the response attributed to entrainment of dry air into the boundary layer as a result of increasing the mixed-phase supercooled liquid fraction achieved by perturbing the GCM’s shallow convective detrainment scheme⁴⁹. This suggests that the response of λ_amt to perturbations in T50 is sensitive to physical parameterizations in GCMs.

Our findings highlight the need to consider changes in the horizontal coverage of clouds in addition to their changes in opacity in response to global warming as they can have a strong opposing effect on climate sensitivity in some GCMs depending on their treatment of cloud physical processes.

Methods

Calculation of T50

The mass fraction of cloud water (‘clw’) and the mass fraction of cloud ice (‘cli’) were first vertically interpolated to the 17 and 19 standard pressure levels of CMIP5 and CMIP6, respectively. The ratio of ‘clw’ to the sum of ‘clw’ and ‘cli’ was then calculated and composited in 1 °C air temperature bins for grid cells at latitudes poleward of 60°. T50 was calculated as the median air temperature at which 0.48 ‘clcalipsoliq’ and ‘clcalipsoice’ in AMIP experiments in place of ‘clw’ and ‘cli’, respectively, however, these variables were only available for a limited number of GCMs participating in CMIP6. Since ‘clcalipsoliq’ and ‘clcalipsoice’ are provided in height coordinates, air temperature was instead vertically interpolated to height coordinates before calculating the ratio of ‘clcalipsoliq’ to the sum of ‘clcalipsoliq’ and ‘clcalipsoice’ and compositing this quantity in 1 °C air temperature bins for all grid cells at latitudes poleward of 60°.

The observed T50 value was also calculated following the same method using the frequency of liquid and ice cloud occurrences and the air temperature at the corresponding locations at 3 °C resolution that are provided by the GOCCP product. Both daytime and nighttime observations were used and the calculations were performed for the time period from 2007-2019, excluding 2016 due to missing data.

Calculation of λ
_τ and λ
_amt

The SW λ_τ and SW λ_amt values were calculated using the approximate partial radiative perturbation (APRP) method⁵⁰ applied to the piControl and abrupt4xCO2 experiments. This method is a simplified version of the partial radiative perturbation method⁵¹. This simplified model assumes that the atmosphere is a single layer that scatters and absorbs radiation, but that incident radiation is only absorbed in its first transit through the layer. SW cloud feedbacks are obtained from the APRP method by first calculating the planetary albedo as a function of seven parameters detailed in ref. ⁵⁰ that are based on the following variables: the cloud area fraction (including stratiform and convective clouds) defined for the whole atmospheric column, as seen from the surface or the top of the atmosphere (CMIP variable “clt”), incident SW radiative flux at the top of the atmosphere (CMIP variable “rsdt”), outgoing SW radiative flux at the top of the atmosphere (CMIP variable “rsut”), outgoing clear-sky SW radiative flux at the top of the atmosphere (CMIP variable “rsutcs”), downwelling SW radiative flux (CMIP variable “rsds”), downwelling clear-sky SW radiative flux (CMIP variable “rsdscs”), upwelling SW radiative flux (CMIP variable “rsus”), and upwelling clear-sky SW radiative flux (CMIP variable “rsuscs”). Akin to the PRP method, the process of calculating a feedback parameter with the APRP method involves running an offline radiative transfer model to calculate the change in radiative fluxes when the variable associated with the feedback parameter of interest from the perturbed climate state is substituted into the control climate state. For example, changes in cloud amount, C, in a control climate (state 1) due to a quadrupling of CO₂ (state 2) would result in a cloud amount radiative feedback that would be calculated as R(C₂, x1, y1, …) − R(C₁, x1, y1, …) normalized by the change in surface air temperature, where x and y are other state variables on which the radiative fluxes, R depend.

To be consistent with previous studies, we take λ_τ to be the summation of the SW scattering and absorption feedbacks³ for clouds at all vertical levels. λ_amt represents the change in cloud amount at all vertical levels due to an abrupt quadrupling of atmospheric carbon dioxide concentrations. Monthly anomalies of the extratropical cloud-induced SW radiative fluxes at the top of the atmosphere were linearly regressed against monthly anomalies of global mean surface air temperature to calculate SW λ_τ and SW λ_amt. The SW scattering feedback is calculated as the sum of the average values of a “forward calculation” and a “backward calculation” linearly regressed onto global mean surface air temperature anomalies. The “forward calculation” is the difference between the planetary albedo calculated with the scattering coefficient from the abrupt4xCO2 experiment and the planetary albedo calculated with all other fields in the piControl state. The “backward calculation” is the difference between the planetary albedo calculated with the scattering coefficient from the piControl experiment and the planetary albedo calculated with all other fields in the abrupt4xCO2 state. The absorption feedback and λ_amt are similarly calculated except with the absorption coefficient and cloud area fraction (“clt”), respectively, in place of the scattering coefficient. The reader is referred to Taylor et al.⁵⁰ for details.

Quantifying uncertainty

The 95% confidence interval of the slope of the linear regression was used to resample the slope 2000 times using a bootstrapping method. This produces a distribution of 2000 additional Δλ_τ values following Eq. (3). The same procedure was repeated for the slope of the T50-λ_amt relationship to obtain a distribution of 2000 additional Δλ_amt values following Eq. (3). These Δλ_τ and Δλ_amt were then inserted into Eq. (4) to calculate a distribution of 2000 additional adjusted climate sensitivity values. The maximum and minimum of these 2000 additional adjusted climate sensitivity values were taken to represent the 5 to 95% confidence interval in Fig. 3.