Introduction

The intertropical convergence zone (ITCZ) is a band of intense deep convection in the Earth’s tropics, forming the ascending branch of the zonally symmetric Hadley circulation. In the longitudinal direction, tropical deep convection is organized on the planetary scale with secondary maxima over near-equatorial South America and Africa, and a maximum in the regions of the so-called Maritime Continent (giving rise to the zonal Walker circulation over the tropical Pacific) and South-East Asia. Changes in the location and structure of the deep convection in the tropics under global warming have profound consequences for tropical climate and the world’s population, and correspondingly the problem has been extensively investigated in models and observations1,2,3,4,5,6,7,8,9,10,11,12. Tropical deep convection is one of the leading sources of uncertainty in future projections of the Earth’s temperature, with major uncertainties in the radiative response of convective clouds, which can have negative or positive radiative effects depending on their optical depth, and the response of convective velocities with warming13,14,15,16,17. The latter is particularly difficult to quantify owing to the small scales involved and is the focus of this paper. Uncertainties linked to convective velocities are pervasive in many aspects of climate projections, from the response of large-scale circulations, where convective velocities modulate the propensity to export energy, to the response of precipitation, which may be enhanced or suppressed by a dynamical feedback on moisture convergence.

Centers of tropical deep convection are characterized by the ubiquitous presence of high ice clouds18 that are formed by detrainment from precipitating deep convection19. At high levels, this detrainment leads to so-called anvil clouds which can persist for several hours after the decay of active convection. Remote sensing techniques that measure ice in the atmosphere thus provide a means to study deep convection globally even though the vertical wind field is not measured directly. Based on such measurements, it has been argued that the tropics respond to surface warming by a narrowing of tropical convergence zones and a decrease of tropical high-cloud fraction11,20,21,22,23,24,25. Generally, previous studies analyzed changes in the occurrence frequency of optically thick ice clouds, which is dominated by the anvil clouds, while the ascending motion in active deep convection is concentrated in a small fraction of the cloudy region, where ice loads exceed about 1 kg m−2 26. Cloud dynamical and microphysical processes, however, render the relation between cloud amount and active deep convection complex27,28. Here, we address the question of how the structure of active deep convection may change under global warming. We explore how, given sufficient horizontal resolution, changes of active convection may be inferred from the ice signal, with the potential to inform on changes of convective velocities with warming. We use measurements and models that resolve km-scale deep convective updrafts, namely active sensor measurements of ice clouds and year-long global warming simulations using a global storm-resolving model (GSRM).

GSRMs are an emerging class of global models with horizontal grid spacings of 2–5 km capable of resolving individual storms29,30. GSRMs can explicitly simulate the range of scales over which tropical deep convection organizes, from individual precipitating deep convective cells to large-scale organization of convection and its interaction with midlatitude dynamics30,31. While GSRMs still rely on parameterizations for turbulence and shallow convection, they are a leap forward in our capability to simulate the response of tropical clouds to global warming32,33,34,35,36,37,38,39. The high horizontal resolution of GSRMs allows direct comparison with active sensor observations, in particular from spaceborne radar such as CloudSat40,41,42. The existence of a continuous record of global CloudSat observations extending over multiple years and with 1.1 km horizontal resolution means that responses to interannual temperature variations can be identified in the observations and compared with results from GSRMs at comparable resolution.

We analyze year-long GSRM simulations of global warming performed at 3 km resolution, made possible by rapid advances in computational resources (see Methods). These simulations use the eXperimental System for High-resolution prediction on Earth-to-Local Domains (X-SHiELD), developed at the Geophysical Fluid Dynamics Laboratory (GFDL)43. Global warming experiments were performed by comparing a year-long control run, forced by sea surface temperatures (SSTs) observed during 2020, with a pair of year-long warming simulations39. These explored the effect of uniformly increasing SST by 4 K without changing CO2 (henceforth, +4K experiment), and then additionally quadrupling CO2 concentration (henceforth, global warming experiment and regarded as a simple analogue of a CO2-warmed climate). The +4K experiment yields a smaller increase of convective available potential energy (CAPE) and peak updraft strength than the global warming experiment, as will be discussed later on, and can therefore be seen as a surrogate of global warming where the response of vertical velocities would be muted. That is, the combination of +4K and global warming experiments allows us to separate the effects of vertical velocities from those of surface warming in the response of convective clouds. Note that the global warming simulation is not in perfect radiative equilibrium (the model has a slightly different climate sensitivity than the implied 4 K for a quadrupling of CO2) and is not coupled to a dynamical ocean.

Results

Response of ice water path frequencies in active convection

Changes of active convection with warming have a signature in cloud ice which can be characterized by the change in the area frequencies of ice water path (IWP; the vertically integrated ice mass per unit area). For context, the IWP of a cloud in the ITCZ is well correlated with the cloud type. Precipitating deep convection carries the highest ice loads, and clouds with IWP > 1 kg m−2 are typically deep convective clouds, with higher IWP indicating deeper convection44,45. IWP in excess of 10 kg m−2 may be found in the most intense convective cores. Convective anvils typically carry lower ice loads (IWP < 1 kg m−2), decreasing as they evolve into thin remnant cirrus.

Figure 1a shows the frequencies of IWP in the control and global warming simulations of X-SHiELD. The frequencies are shown for the tropics (20°S to 20°N) and are defined as

$$f\left({{{\rm{IWP}}}}\right)=\frac{1}{{S}_{{{{\rm{tropics}}}}}}\frac{\partial S}{\partial {\log }_{10}{{{\rm{IWP}}}}},$$
(1)

where S is the surface area of columns with ice water path below IWP (i.e., the cumulative distribution function of IWP) and Stropics is the surface of the tropics from 20°S to 20°N. IWP statistics have been computed at 3-hourly resolution and results have been averaged over the full seasonal cycle in each simulation. In order to deal with the considerable volume of data generated by year-long simulations, the statistics have been coarse-grained (horizontally averaged) from the native 3 km resolution to a (25 km)2 grid. We present limited results at native resolution in the next section to estimate the impact of such averaging. In Fig. 1a, active convection (IWP > 1 kg m−2) occupies only 2.7 % of the tropical area (in the control simulation). Most of the ice is found in recently detrained anvils, with a mean IWP integrated over the distribution of 0.1 kg m−2 (vertical red line in Fig. 1a). IWP frequencies increase for smaller IWP because anvil clouds spread horizontally after detraining from convection45. The simulated annual-mean changes in IWP with increasing surface temperature and CO2 are reductions for much of the tropics with some regions of increase (Supplementary Fig. 1).

Fig. 1: Change in the frequencies of IWP in the global warming GSRM experiment and with interannual temperature variations in the CloudSat observations.
figure 1

a Simulations in present-day climate (control) and in the warmer future climate (global warming). The ranges of IWP corresponding to active convection and convective anvils are indicated in the panel. The vertical red line with the cross underneath marks the mean IWP integrated over the control distribution. The plotted range of IWP frequencies, i.e. above 8 × 10−2 kg m−2, represents about 14 % of the tropical area. b Interannual anomalies of tropical-mean surface temperature TS in the ERA5 reanalysis. c Interannual anomalies of IWP frequency from July 2006 to December 2010 in CloudSat observations (2B-CWC-RO). d Fractional change of IWP frequency with TS in the simulation and with interannual temperature variations in CloudSat observations. The red line is for observations from 2006 to 2010 (2B-CWC-RO) and the yellow line is for observations from 2006 to 2017 (DARDAR). See text and Supplementary Fig. 6 for discussion of the differences between 2B-CWC-RO and DARDAR. The error bars indicate the standard errors of regression of the IWP interannual frequency anomalies on the TS interannual anomalies (see Methods; the confidence interval on the correlations and associated p-values are shown in Supplementary Fig. 4).

Full size image

The change with warming is characterized by a broad decrease of IWP frequency below IWP 4 kg m−2 and an increase of IWP frequency above this value. This signature of warming indicates a decrease of anvil cloud amount and a decrease of the frequency of all but the most intense convective cores. (See also Fig. 1d).

Impact of coarse graining

Given the volume of data generated by year-long simulations at kilometric resolution (1600 TB per year per simulation in the case of X-SHiELD), horizontal averaging of three-dimensional output fields such as ice mass is necessary. The strategy deployed is to perform an online coarse-graining from the native resolution to 25 km. Note that the 25 km output of a 3 km model is different from the output of a 25 km model. In particular, convective accelerations are explicitly resolved by a 3 km model, so that its simulated sensitivity to environmental factors, such as convective available potential energy (CAPE), is preserved in the 25 km output. To estimate the impact of coarse graining, we perform limited 10-day simulations of the control and global warming runs and compare the response of IWP frequencies at 3 km and 25 km (Supplementary Fig. 2). Those runs are not fully equilibrated and cannot be compared to those of Fig. 1, but a similar structure of response is observed at high IWP, albeit the switchover from lower to higher occurrence frequency shifts to a higher IWP in the 3 km data. A similar shift is observed when comparing output of the year-long simulations coarse-grained at 250 km, 100 km and 25 km (Supplementary Fig. 3). Based on this comparison, we would speculate that coarse graining to a (25 km)2 grid preserves the basic structure of IWP response to global warming, but translates the response to lower IWP values. That is, we expect that the transition from decrease to increase of IWP frequency occurs above IWP 4 kg m−2 at 3 km resolution.

Evidence of a similar response in the observations

There is evidence that a similar signature is observed in the response of convective cloud ice with interannual temperature variations in the tropics. The observations use IWP estimates from CloudSat cloud radar in the tropics, namely from the Level 2B Radar-only Cloud Water Content retrievals (2B-CWC-RO)46 in the same latitude range as in the model. Only for the last 15 years have accurate observational estimates of IWP from such sensors been possible, and the period of nominal CloudSat radar observations is only 4.5 years, which is too short to directly assess the sensitivity of IWP frequencies to greenhouse warming. The observations do cover, however, periods of different mean temperature due to internal variability from El Niño/Southern Oscillation (ENSO), as shown in Fig. 1b. As an intriguing proxy, we look at the response of IWP frequencies with interannual variability in tropical-mean surface temperature TS, derived from ERA5 reanalysis, shown in Fig. 1c. The proxy is imperfect because ENSO is associated with large changes in the geographic distribution of deep convection and underlying sea surface temperatures (see Supplementary Fig. 1), but we again find evidence of a rectified dependence of IWP frequencies on TS (which tends to warm during El Niño and cool during La Niña) comparable to the GSRM sensitivity to global warming.

We infer the sensitivity with warming, shown in Fig. 1d, from the slope of a linear regression of IWP frequency anomaly on TS anomaly (see Methods). IWP frequencies have been accumulated monthly in the tropics (from 20°S to 20°N) from radar profiles with a resolution of 1.1 km and 0.16 seconds along track. No attempt has been made to coarse grain the observations because horizontal averaging along track and on a grid have different properties. The frequencies follow the definition in Eq. (1), ensuring intercomparability with the model. Because CloudSat cannot retrieve profiles with IWP less than 10−4 kg m−2 47, we keep track of tropical columns with no detectable IWP to ensure the correct normalization of IWP frequencies. Interannual anomalies are then computed by detrending, deseasonalizing and applying a 3-month moving mean to the monthly frequencies. Error bars of the observational slopes are also plotted in each IWP bin in Fig. 1d. Because the observations span less than two ENSO cycles, the IWP observations have very few effective degrees of freedom, so the error bars are large and the slopes are not statistically significantly different from zero at 95% confidence (Supplementary Fig. 4). However, the functional form of the model response and the observations are remarkably similar for IWP > 1 kg m−2, i.e., in the range assumed to be active convection.

In the anvils (IWP < 1 kg m−2), a robust anticorrelation is observed between interannual TS and IWP frequency, i.e. a broad decrease of anvil clouds during periods of anomalously high tropical average surface temperature (i.e. July 2006–July 2007 and July 2009–July 2010), and the opposite for average and colder conditions. This response is qualitatively consistent with the decrease of IWP frequency with warming simulated by the GSRM experiments, but of larger magnitude. Possible reasons for this discrepancy include unconstrained microphysical degrees of freedom, since anvil cloud dynamics is strongly dependent on cloud microphysical processes which remain poorly resolved in GSRMs36, and the differing spatial patterns associated with uniform warming in the model and interannual warming in the observations (Supplementary Fig. 1).

In active convection (IWP > 1 kg m−2), on the other hand, the fractional changes of IWP frequency in the CloudSat observations show the signature identified in the GSRM experiments, with quantitative consistency (although a lot of uncertainty). Similar results are obtained with other active remote sensing datasets, including the Level 2B Radar-only Cloud Water Content assuming Ice Only retrievals (2B-CWC-RO IO) and the raDAR/liDAR retrievals (DARDAR) (see Supplementary Fig. 5)46,48.

In the above discussion, we have restricted the observations to the period of nominal CloudSat operations (when the radar operated both in daytime and nighttime) prior to a battery anomaly in 2011 (Fig. 1d, red line). An extension of this analysis to 2012 – 2017 using the DARDAR data48 (Fig. 1d, yellow line), after the battery anomaly limited the radar to daytime observations only, shows very similar changes of IWP frequencies in active convection, with smaller error bars and correlations significant at the 80 % level (see Supplementary Fig. 6) which increases the confidence in the analysis. We are largely focused on the basic structure of response of active convection (transition from a decrease to an increase in IWP frequency, with a change of sign of the correlation) rather than in the quantitative sensitivities. The differing spatial patterns of warming complicate comparisons between the model and the observations. Differences in vertical wind shear environment between El Niño and La Niña, or in mechanisms of freezing mediated by aerosol interaction, can affect the sensitivity to warming of IWP frequencies. Note also that CloudSat observations use fixed overpass of the satellite at 01:30 and 13:30 local time, missing the late afternoon maximum of land convection. The identification of a similar structure of response in model and observations suggests the existence of a mode of change with interannual variability that is independent of the spatial reorganization of convection. In the following, we focus on the changes of IWP in active convection (IWP > 1 kg m−2) where model and observations show agreement.

Interpretation of repartitioning of active convection

We further analyze and interpret how IWP in and near active convection (IWP > 1 kg m−2) changes under warming in the X-SHIELD simulations. An easily visualized approach, following a technique used in self-aggregation studies for water vapor path49, is to sort the tropical 25 km coarse-grained IWP into percentiles, each accounting for 1 % of the fractional area of the tropics (we use a bin size of 0.02 %). Thus we transform the data from IWP space to IWP percentile space, which ensures that comparisons are performed over regions of equal surface area, with the surface always ordered by decreasing ice load.

Figure 2a shows the IWP in IWP percentile space (sorted in ascending order) in the control and global warming simulations of X-SHiELD. About 2.7 % of the tropics have IWP > 1 kg m−2, and we therefore truncate the data below the 95th percentile. Only about 0.02 % of the tropics have IWP > 10 kg m−2. The response of IWP at fixed IWP percentile under warming (see Fig. 2a, comparing the control and global warming simulations; and Fig. 2b showing the fractional change per unit warming) is characterized by an increase of IWP above the 99th IWP percentile and a decrease of IWP below the 99th IWP percentile. That is, the ice loading in the very active deep convective locations is increasing, but is decreasing under less active conditions. We note that these changes are on the scale of convective systems in the tropics and that precipitation shows a similar response in the IWP coordinate (see Supplementary Fig. 7).

Fig. 2: Response of IWP at fixed IWP percentile in the GSRM experiments.
figure 2

a IWP in present-day climate (control) and in the warmer future climate (in the global warming simulation). b Fractional increase of IWP with tropical-mean TS at constant percentile in the global warming and +4K experiments.

Full size image

In order to better understand this signal, we also analyzed an X-SHIELD simulation in a +4K configuration where the surface temperature is raised by 4 K, but atmospheric CO2 is held fixed at the level of the control run. Comparison of the fractional changes with warming between the global warming and the +4K experiment (Fig. 2b; blue and yellow lines) shows that the basic structure of the response is similar, but noticeably weaker in the +4K experiment. Hence, the signal results both from the warming and the effect of increasing CO2.

Thermodynamic and convective mechanisms for ice water path changes

Additional insights into the mechanism may be obtained by considering the vertical structure of the ice mixing ratios as a function of IWP percentile and temperature, shown in Fig. 3. The use of temperature coordinates is motivated by the fact that, while tropical ice clouds rise to lower pressures with warming, their vertical structure is more invariant in these coordinates. This is explained by the temperature-invariance of several microphysical processes, such as the onset of freezing occurring at 0°C and the temperature constraint on the existence of supercooled liquid water droplets, combined with the radiative cooling constraint on the outflow level of deep convection known as the Fixed Anvil Temperature (FAT) hypothesis50.

Fig. 3: Response of cloud ice in IWP percentile and temperature space in the GSRM experiments.
figure 3

Fractional increase per TS warming of ice mixing ratio at fixed temperature and fixed IWP percentile in (a) the +4K experiment and (b) the global warming experiment. The white dashed line shows the 7 % K−1 isoline where the increase in ice mixing ratio approximately corresponds to the increase in boundary layer specific humidity. The IWP response from Fig. 2 is reproduced in the upper panel of a, b.

Full size image

The IWP in temperature coordinates is given by

$${{{\rm{IWP}}}}=\int-{q}_{{{{\rm{i}}}}}\left(\frac{\partial p}{\partial T}\right)\frac{dT}{g},$$
(2)

where qi is the mixing ratio of cloud ice, T is temperature, p is pressure and g is the acceleration of gravity. It follows that a change in IWP may arise from changes in cloud ice mixing ratio qi and/or from changes in the pressure-temperature relation ∂p/∂T in the atmospheric column under consideration. When the atmosphere warms (roughly) as a moist adiabat, temperature increases faster at low pressures than at high pressures, thus IWP is subject to a negative tendency due to the term ∂p/∂T, while qi changes on isotherms must have a non-local origin such as a change in the vertical advection of ice. Figure 3 shows that, for both the global warming and +4K experiments, the ice mixing ratio increases over the full vertical profile at the highest IWP percentiles at a rate of several percent per Kelvin.

The increase of cloud ice mixing ratio qi at an isotherm under warming may arise from two main effects. To understand this, one may define an ice conversion efficiency ϵc relating qi to its sources from condensation (generating liquid water which subsequently freezes) and deposition:

$${\epsilon }_{{{{\rm{c}}}}}:= \frac{{q}_{{{{\rm{i}}}}}}{{q}_{{{{\rm{sat}}}}}({T}_{{{{\rm{S}}}}},{p}_{{{{\rm{S}}}}})-{q}_{{{{\rm{sat}}}}}(T,p)},$$
(3)

where the difference of saturation mixing ratio in the denominator is the ice source term on a path connecting the isotherm to surface temperature TS. For ϵc = 1, qi is equal to the amount of water repartitioned from the vapor to the ice phase in a closed system, whereas the limit of zero conversion efficiency means that all ice has precipitated out. The value of ϵc reflects effects of vertical advection, sedimentation, microphysics and other small-scale processes (such as entrainment of dry air into convective clouds). The source term in the denominator of Eq. (3) reduces to qsat(TS, pS) at low temperatures, which is sensitive to surface temperature.

This combined sensitivity to changes in surface temperature and ice conversion efficiency explains the differing responses in the global warming and +4K experiments. Two modes of change are operating. First, the warming corresponds to a deepening of the atmospheric temperature range and the saturation mixing ratio at the cloud base increases (or, more intuitively, the subcloud layer is warmer and can hold more water vapor), giving rise to a thermodynamic mode. Because the warming of the subcloud layer is similar in the global warming and +4K experiments, the thermodynamic mode is similar. Clausius-Clapeyron scaling (at constant relative humidity) of subcloud water vapor mixing ratio would imply a fractional change of cloud ice mixing ratio at a rate of about 7% K−1 in the global warming and +4K experiments. The different responses in the global warming and +4K experiments must therefore be explained by differences in ice conversion efficiency. Conversion efficiency may change because the updraft velocity in the convective cores may increase8,10, generating and suspending more ice. In the +4K experiment, the increase of ice mixing ratio in the most intense convective percentiles is less than 7 % K−1 (about 4 % K−1), indicating that vertical advection cannot transport all the additionally available subcloud moisture. In the global warming experiment, there is an entire domain in temperature and IWP (at the top at around − 70°C, see white dashed line in Fig. 3b) where the fractional increase is above 7 % K−1 (up to 10 % K−1). The larger fractional increase indicates that in this experiment vertical updrafts are more vigorous than in the +4K experiment. Indeed, analysis of the convective available potential energy (CAPE) and column-maximum updraft velocity (\({w}_{\max }\)) (see Table 1) shows that CAPE and \({w}_{\max }\) increase more in the global warming than in the +4K experiment, due to slightly enhanced boundary layer humidity across the tropics.

Table 1 Amplification of CAPE and extreme column-maximum updraft velocity at 3 km resolution in the GSRM experiments from 20°S to 20°N.
Full size table

At less extreme IWP percentiles, there is a weaker response of cloud ice to surface warming, with even a decrease between 0°C and −20°C, in the temperature range where mid-level cumulus clouds are detraining. This is consistent with less intense convection being diminished in a warmer climate, even as the most intense convection is strengthening. At temperatures below −40°C in the upper troposphere, the increase of cloud ice with surface warming at less extreme IWP percentiles is likely due to detrainment from nearby active convection.

Discussion

The GSRM results (Figs. 2 and 3) suggest that the increase of IWP at the highest IWP percentiles in response to warming in the model simulations and observational data is strongly dependent on changes in the vertical wind at the convective scale. The global warming simulation (+4K SST and increased CO2) and the +4K simulation have identical SSTs and as such similar boundary layer-specific humidity, but the former yields a much stronger increase in ice mixing ratios over the control simulation, which exceeds at upper levels the approximate Clausius-Clapeyron scaling of boundary layer moisture. Figure 4 provides additional context to interpret the changes in response to warming.

Fig. 4: Circulation in convective clouds computed at 3 km resolution in the tropics (20°S to 20°N) from limited model output (20 days) of the control simulation in August.
figure 4

a Column-maximum updraft velocity \({w}_{\max }\) (black dashed line) and column ice conversion efficiency (gray line with diamond markers) at a given IWP. b Mass streamfunctions in IWP and temperature space (with corresponding IWP percentile indicated below). The black dot marks the position of maximum streamfunction. The streamfunctions have units of 1011 kg s−1 and the contours are equal to 0.01, 0.05 and then equally spaced between 0.25 and 5 with 0.25 spacing. The column-integrated ice conversion efficiency is sensitive to the definition of the upper bound of the ice cloud (see Methods). The upper bound is defined as the highest level where qi exceeds 2.5 g kg−1 (gray line with diamond markers).

Full size image

Figure 4 shows the column-maximum vertical velocity \({w}_{\max }\) (panel a) and the time-mean mass streamfunction (panel b) in IWP percentile space of the base simulation over a limited period of 20 days where data was stored at the full horizontal resolution of about 3 km. Figure 4a further shows the column-integrated ice conversion efficiency, as defined in Eq. (3), calculated as the ratio of the model total column ice mass divided by the approximate ice mass if it came from perfect conversion of subcloud-specific humidity (see Methods). The figure shows that the time-mean mass streamfunction switches from ascending to descending (black solid circle) around IWP ~ 1 kg m−2, qualitatively consistent with the frequently used classification of IWP in excess of 1 kg m−2 to be indicative of deep convection. Figure 4a further shows—consistent with expectation—the strong connection between updraft speed and ice conversion efficiency, with the ice conversion efficiency approaching unity in the highest IWP percentiles of about 50 kg m−2 (consistent with the expectation based on a closed system, see Methods), where the mean column-maximum vertical velocity in X-SHiELD reaches about 20 m s−1.

Returning to Figs. 2 and 3, the base state shown in Fig. 4 supports the notion that the higher response of IWP at the most convective percentiles in the global warming vs the +4K experiment is explained by a greater increase of column ice conversion efficiency following a greater increase of \({w}_{\max }\). More generally, the results show that column ice conversion efficiencies will change with warming as a result of the changing distribution of vertical velocities and that the changes will be registered in the ice field.

The broader implications of these results are twofold. First, there is an implication that the optical depth of convective clouds, a quantity closely related to IWP, and its changes with climate change, are partially controlled by convective dynamics and may not simply scale with subcloud specific humidity. This sensitivity may be difficult to capture with parametrizations of deep convection used in current General Circulation Models, and requires further investigation due to its impact on the cloud radiative feedback. Second, the model results are important for the interpretation of the information provided by cloud ice observations. The GSRM simulation suggests that convective dynamics have a strong bearing on cloud ice, and consequently observations of cloud ice may be used to extract observational constraints on convective dynamics and the changes thereof with warming.

The computational costs of long-term integrations of GSRMs presently remain very high even for large high-performance computer systems. The year-long X-SHiELD model simulations used here with a resolution of 3 km are a leap forward in terms of their ability to resolve deep convection self-consistently embedded in the general circulation. However, the simulations may be subject to biases due to, for example, the simplicity of the model’s microphysics scheme. Our analysis shows the potential of these model simulations to interpret high-resolution observations such as provided by the CloudSat radar measurements. These observations offer an important lens to analyze kilometer-scale simulations, which are a key avenue to improve quantitative understanding of the processes controlling deep convective clouds and their changes under global warming. The results presented here indicate that deep convective ice clouds will experience changes at the kilometer-scale and mesoscale, and that these changes are manifest in the presently available high-resolution CloudSat radar measurements.

Methods

Model

X-SHiELD is the experimental global storm-resolving model configuration of the unified forecast system SHiELD, developed at NOAA/GFDL43. It is powered by the Finite-Volume Cubed-Sphere (FV3) Dynamical Core51,52. The horizontal resolution is 3.25 km globally and the model uses a “vertically Lagrangian" discretization with 79 levels in the vertical. The physical parameterizations include the in-line GFDL microphysics scheme53, the turbulent kinetic energy-based moist eddy-diffusivity mass-flux vertical turbulence mixing scheme54, a simplified Arakawa-Schubert scheme for shallow convection55, the Noah-MP land surface model56 and a mixed layer ocean model57. SSTs are nudged to European Centre for Medium-Range Weather Forecasts (ECMWF) analyses with a timescale of 15 days, following the DYAMOND protocol30. Three sets of simulations are compared: a control simulation (control), a warming simulation where SSTs are uniformly increased by 4 K without variation of CO2 (+4K) and a warming simulation where CO2 is increased to 1270 ppmv (global warming). Both simulations are integrated for 15 months from October 2019 to January 2021. The present study uses model output that has been coarse-grained at 25 km resolution in the tropics from 20°S to 20°N in the period from Jan 2020 to Jan 2021. IWP and ice mixing ratio are computed by summing the values for cloud ice, snow and graupel.

Satellite observations

We use the Level 2B Radar-only Cloud Water Content product (2B-CWC-RO) Version P1_R05 to compute monthly-aggregated IWP frequencies in the tropics from 20°S to 20°N46. We use the record of observations from July 2006 to December 2010, corresponding to nominal operations prior to a battery anomaly of CloudSat in 2011. We also use the Level 2B Radar-only Cloud Water Content product assuming ice-only retrieval (2B-CWC-RO IO) Version P1_R05 and the raDAR/liDAR ice cloud property retrieval product (DARDAR) Version 3.1048 for comparison and to extend the analysis up to December 2017, after the battery anomaly limited the radar to daytime observations only (Supplementary Fig. 6). Interannual anomalies are computed by detrending, deseasonalizing and applying a 3 month moving mean to the monthly time series. The sensitivity to tropical-mean surface temperature TS is computed by regressing against the interannual anomalies of 2-meter temperature derived from ERA5 reanalyses in the tropics58. The p values of the correlations between TS and IWP frequency anomalies (probability of observing a correlation as extreme if the true correlation were zero) are computed using a two-sided Student’s t-test (Supplementary Figs. 4 and 6). The confidence intervals on the correlations are computed at the 95 % level using a Fisher’s z-transformation (Supplementary Figs. 4 and 6). In computing the p values, confidence intervals and the standard errors of regression, we estimate the effective number of degrees of freedom due to autocorrelation in the timeseries by using the methodology described in Bretherton et al.59. We estimate the reduction in degrees of freedom to be a factor 3.

Mass streamfunctions in IWP and temperature space

The streamfunction Ψ is computed in IWP and pressure space as

$$\Psi \left(p,{{{\rm{IWP}}}}\right)={S}_{{{{\rm{tropics}}}}}\int\nolimits_{{{{\rm{IWP}}}}}^{+\infty }f\left({{{\rm{IWP}}}}\right)\left\langle \frac{-\omega }{g}\right\rangle \left(p,{{{{\rm{IWP}}}}}^{{\prime} }\right)d\,{{{\rm{IWP}}}}^{\prime} ,$$
(4)

where ω is pressure velocity, g is gravitational acceleration, \(f\left({{{\rm{IWP}}}}\right)\) is the IWP area frequency defined in Eq. (1) and Stropics is the surface of the tropics from 20°S to 20°N. The angle brackets 〈〉 denote a conditional average in IWP and pressure space. \(\Psi \left(p,{{{\rm{IWP}}}}\right)\) is thus the net vertical mass flux from 20°S to 20°N at level p of all air parcels with ice water path greater than IWP. Ψ is then interpolated in temperature to yield results in IWP and temperature space.

Column ice conversion efficiency calculation

The theoretical IWP from repartitioning of subcloud moisture is computed using a simple adiabatic model of updraft, assumed to be a closed system. The model assumes the conservation of frozen moist static energy:

$${c}_{{{{\rm{p}}}}}T+gz+{L}_{{{{\rm{v}}}}}{q}_{{{{\rm{v}}}}}-{L}_{{{{\rm{f}}}}}{q}_{{{{\rm{i}}}}},$$
(5)

where T is temperature, z is geopotential height, qv and qi are vapor and ice mixing ratios, cp is the specific heat at constant pressure for air and Lv and Lf are the latent heats of vaporization and freezing. qv assumes a saturation adjustment over liquid water for T > 0°C, following Bolton 198060, and over ice for T < 0°C, following Murphy and Koop 200561. Geopotential height is computed using the hydrostatic approximation. Liquid water and ice mixing ratios are computed by repartitioning the subcloud water vapor for a closed system, assuming instantaneous freezing at 0°C. The model is integrated from the surface (assumed at 1000 hPa and z = 0) with an initial temperature of 300 K and relative humidity of 80%. The integration is stopped at a level assumed to be the top of convective ice, diagnosed from the ice field in the model as the highest level where qi exceeds 2.5 g kg−1 in IWP and pressure space. Column ice conversion efficiency is defined as the ratio of the IWP of the atmospheric column to the theoretical IWP computed using the adiabatic model.