Introduction

Aerosol-induced changes in cloud structure have clear implications for cloud radiative effects and precipitation1,2,3,4,5 via aerosol–cloud interactions.6,7 Until now most attention has focused on low clouds (i.e., non-precipitating or warm-rain), where reduction of cloud droplet size is well established8,9,10,11,12 but the net effect on the albedo and lifetime of the low clouds remains controversial.13,14,15,16 For deep-convective clouds, which contain frozen and liquid hydrometeors, these effects are extremely complex and poorly understood due to the complicated dynamics, thermodynamics and microphysics of these clouds.

It has been hypothesized that aerosols can affect deep convective clouds via a process, here called “Aerosol–cloud invigoration,” where an increase in aerosol loading deepens the convective clouds due to the strong coupling between the cloud microphysics and cloud system dynamics.17,18 Numerous observation-based studies have reported higher cloud tops and cloud cover with greater aerosol loading;19,20,21,22 however, difficulty in discriminating aerosol–cloud interaction from meteorological covariations23 and uncertainties in satellite retrievals24,25,26,27,28 are major limitations for quantifying these effects.29 Additionally, there has been no observational estimate of radiative forcing by aerosol–cloud interactions in deep convective clouds.

These major limitations in observing aerosol–cloud effects call for the use of numerical modelling. A number of studies have examined aerosol effects on deep convection in cloud-resolving, global and regional models.18,30,31,32,33 However, representation of cloud microphysics and aerosol–cloud interactions in numerical models remains a major challenge, particularly in global models, while the lack of feedbacks between clouds and their large-scale environment may compromise findings from cloud-resolving models.34,35,36 Several studies have employed regional model downscaling over continents.37,38,39 Such simulations represent meteorology realistically, but a limitation is that some state variables are typically relaxed to a global field, which could damp perturbations involving changes to such variables. Some model studies suggest aerosol invigoration is minimal when averaged over time and space, due to larger-scale adjustments, even if localized instantaneous effects do occur.40,41

In this study, we use a novel approach of simulating aerosol effects in a closed, idealised, meteorologically heterogeneous domain with a large-scale overturning circulation and an organized zone of deep convection. To our knowledge this is the first model study to consider aerosol impacts on organised deep convection in a closed domain with large-scale circulation. This experimental setup allows the representation of the environmental feedbacks related to changes in large-scale circulations. Furthermore, in order to address uncertainties about aerosol process representation, we employ multiple approaches for representing aerosol effects in the model. Finally, we show that model-predicted effects are supported by observations and suggest and test an explanatory mechanism.

Our numerical experiments use the Weather Research and Forecasting model version 3.5.1.42 To represent a large-scale tropical overturning circulation, we adopted a setup used in previous studies43,44 with domain boundaries at the Equator and 30° latitude (see Methods). We represented aerosol effects in the model using two different approaches. In the first or “proxy heating” approach, we assumed40 that additional latent heat is released in the presence of larger number concentration of cloud concentration nuclei (CCN) due to the freezing of more cloud droplets lifted above the freezing level17 and used this to artificially represent the aerosol effect. These proxy heating and cooling perturbations are applied in updrafts and downdrafts respectively, such that the net heat input is zero. In the second approach, we used a microphysical scheme with specified cloud droplet and CCN concentrations. For polluted cases, we increased the CCN or cloud droplet concentration (depending on the microphysical scheme used) by a factor of eight (see Methods), without changing the size spectra.

Experiments were performed at two different resolutions: 2 km (convection-permitting) and 20 km (using convective parameterization). We moreover used three different combinations of microphysical (Double-Moment Six-Class (WDM6)45 or Morrison Double Moment Four-class (Morr2)46) and shortwave radiation (Dudhia47 or Rapid Radiative Transfer Model (RRTM)48) physics parameterizations, denoted a/b/c (respectively: Morr2+RRTM, WDM6+RRTM, and WDM6+Dudhia). Experiments with proxy heating perturbations and high CCN/cloud droplet concentrations are denoted as PERT and POLL respectively, followed by number 2 or 20 denoting the spatial resolution. Detailed explanations of the experiments and the idealized model setup are given in Table 1. Results shown in this paper are based on the last three months of each simulation.

Table 1 Table of experiments.
Full size table

The unperturbed Control20c simulation was already evaluated43 and found to qualitatively resemble the observed meridional overturning and cloud structure in the eastern Pacific. The simulated cloud fraction in the convergence zones however overestimated observed (CALIPSO-GOCCP) cloud fraction by 30%, perhaps because the strong and linear SST gradient produces a very narrow ITCZ.43 The other Control experiments produce qualitatively similar results, except that some (Control20a and Control20b) underestimate shallow cloud cover (not considered here).

Stronger updraft velocities (Supplementary Fig. 2) occur when aerosol effects are represented in the model, irrespective of whether the aerosols are represented as the proxy perturbations to latent heating (PERT20 and PERT2) or as high CCN concentration (POLL2). On average across physics versions, the updrafts are stronger by 0.31, 0.14, and 0.19 m/s in the POLL2, PERT2, and PERT20 experiments, respectively (Supplementary Table 4). For the PERT experiments this invigoration can be attributed to the heating perturbations directly, which make updrafts more buoyant allowing the cloud parcels to rise faster and reach higher in the atmosphere. On the other hand, in the POLL experiments, the nucleation of smaller cloud droplets delays warm precipitation so that the clouds do not precipitate before reaching the super-cooled levels (~4 km altitude), whereupon freezing releases latent heat of fusion.1,6 This additional latent heating (Supplementary Table 4) in the upper troposphere makes the updrafts buoyant, invigorating the clouds.

Previous work has suggested that such intensification produces higher and more extensive cirrus, implying a net warming effect on climate.6,49 We find however that once the atmospheric temperature and dynamics adjust to the aerosol, there is a net top-of-atmosphere cooling effect (discussed later in the paper). Intensification of updrafts in the perturbed experiments is instead accompanied by the lateral spreading and an overall increase in the cloud fraction of the mid and high-level clouds (Fig. 1 and 2a, b) and latitudinal shrinking and concentration of the near-equatorial rain band (Fig. 2c, d). Clouds initially concentrated near the Equator now spread farther away from it, accompanied by a slight reduction at the Equator itself (Fig. 1). This reduction can be attributed to the concentration and enhancement of rain which desiccates the cloud. The results are robust for different resolutions, radiation and cloud microphysical schemes (Supplementary Figs. 4 and 5).

Fig. 1
figure 1

Impact of aerosols on simulated cloud amount. Zonal and temporal averages of difference in cloud fraction between a POLL2b and Control2b, b PERT2b and Control2b, and c PERT20b and Control20b, all using the WDM6-RRTM physics version. (see Supplementary Fig. 3 for other physics versions).

Full size image
Fig. 2
figure 2

Impact of aerosols on simulated zonal-mean precipitation and total cloud. Zonally, temporally and vertically (from 6–11 km) averaged spatial profiles of a, b cloud fraction and c, d precipitation, for a, c 2-km and b,d 20-km resolution experiments. Each profile represents the average of the results using the three physics versions (see Supplementary Figs. 4–5 for results separated by version).

Full size image

Simulations using the Morr2 microphysics scheme (case a) produced twice the cirrus cloud cover as those using the other scheme (cases b,c). The corresponding perturbed experiments (PERT20a, PERT2a and POLL2a) show strong changes in cirrus (Supplementary Fig. 3a, d, g). Another (PERT20c, Supplementary Fig. 3c) shows a reduction in the shallow clouds, due to the behavior of the cumulus parameterization scheme40 (Supplementary material). However, neither of these changes are robust across the suite of model experiments and they will be ignored henceforth.

Past studies have suggested that factors like mass flux largely control the convective cloud amount.50,51 However, all the experiments except PERT2a showed a decrease in mass flux (weaker by 0.09, 0.05, and 0.02 kg/m2 s on average in the POLL2, PERT2 and PERT20 experiments, respectively), see Supplementary Fig. 6 and Supplementary Table 4 for more details). Hence the increase in convective cloud amount (between 6 and 12 km altitude) found here is not due to any increase in updraft mass flux.

The perturbed experiments show a warmer atmosphere (by 0.5 °C averaged from 6–12 km, the upper troposphere) and weaker lapse rate in comparison to the control experiments (Fig. 3). Warmer updrafts enhance upward heat transport, with warmer convective cores rapidly warming the environment via gravity waves.52 There is also however some warming of even the lower troposphere, which (Supplementary Fig. 7) appears to be associated with the strengthening of shallow circulations. The 0.5 °C upper-troposphere warming at fixed surface temperature represents a change in stability that is significant compared to typical time variations.53 We propose that this stabilization is responsible for enhancing the cloud fraction throughout most of the troposphere. A more stable environment should lead to more extensive and long-lived clouds, because instability tends to break up clouds via convection which leads to water removal by precipitation and entrainment of dry air into the cloud. A strong positive correlation between static stability and cloud cover in the lower troposphere is well known.54,55 Our results show a similar-magnitude relationship for deeper clouds: for a ~0.04 K/km change in tropospheric stability (between surface and 14 km altitude), the POLL2, PERT2 and PERT20 experiments show a 0.95%, 1.08%, and 0.60% increase in the convective cloud fraction (between 6 and 12 km), respectively (Supplementary Table 4). This represents roughly 2% per 0.1 K/km lapse rate change, in the same direction but much greater than the reported 0.28% increase per 0.1 K/km between surface and 3 km;54 however the two are of closer magnitude if compared on the basis of bulk temperature change across the respective cloud layer.

Fig. 3
figure 3

Impact of aerosol invigoration on atmospheric temperatures. Difference in temperature between perturbed and control simulations in the convective region (0°–10° latitude) for the last 3 months of the simulation.

Full size image

To test the lapse-rate hypothesis, a suite of stratification experiments is performed at both resolutions, to investigate how the lifetime and amount of clouds respond to change in stability. In two perturbed runs the model is restarted with a dry boundary layer to prevent development of further convection; in one of these (DBL) this is the only change, while in the other (DBLW), the lapse rate was stabilized by 0.04 K/km between surface and 14 km altitude, without altering other prognostic variables. These experiments may be compared to each other and to the unperturbed control run (CTRL).

Cloud water is continuously removed by precipitation formation and replenished by convective water transport. Both of these effects are driven by moist turbulence. The drying of the boundary layer in both experiments inhibits this turbulence, suppressing precipitation (Supplementary Fig. 10) but with little immediate effect on relative humidity above the boundary layer. Due to this, both experiments show immediate increases in clouds above the boundary layer (during the first few hours) (Supplementary Fig. 8), but these are stronger and longer-lived in the DBLW experiment than in DBL.

The stabilised experiment (DBLW) shows greater cloud fraction (relative to unstabilised DBL) by 0.54% and 0.82% for the 20 and 2 km resolution experiments respectively (see Supplementary Fig. 8b, d, f, h). This is a change of 1.3–2% per 0.1 K/km, approximately comparable to those noted above for the aerosol experiments, supporting the hypothesis that the lapse rate changes are causing the cloud changes. A role for large-scale static stability in promoting free-tropospheric cloudiness has been reported previously,56 where absorbing aerosols increased GCM-simulated cloud cover when placed high in the troposphere and reduced it when placed low. Thus, our idealised stratification experiments, and GCM studies,56 offer support for our hypothesis that stabilisation of the atmosphere is the means by which aerosol affects cloud cover.

A more stable troposphere can also explain the contraction of the rain zone seen in our experiments. A weaker lapse rate means more boundary-layer moisture is needed for convective instability. To accumulate additional moisture, equatorward-moving air would have to move closer to the equator, an effect that has previously been dubbed the “upped ante” mechanism.57 This leads to enhanced precipitation very close to the Equator but a reduction in precipitation at the margins (Fig. 2c, d), i.e a narrowing of the region of deep convection rain-band. All the perturbed experiments except PERT20a and PERT20b show a narrower rain band and a reduction in mean precipitation (Supplementary Table 4). Note that the stabilisation does not reduce average convection or precipitation, which are ultimately controlled by the atmospheric energy budget; in effect the extra stability acts as a buffer that limits net invigoration, consistent with previous results.40

The fact that this mechanism causes opposing changes in rain and cloud area provides a unique opportunity to detect it in observations by directly examining the ratio of these two areas, which normalizes out some of the local meteorological noise that has hobbled past efforts to detect subtle aerosol impacts. We analyzed 10 years of de-seasonalised, monthly mean aerosol, cloud, and rain observations. For aerosol we used three alternative measures: 550-nm aerosol optical depth (AOD) and aerosol index (AICAMS) from the Copernicus Atmospheric Monitoring Service (CAMS) global aerosol reanalysis,58,59,60 and the aerosol index (AITOMS) from the Total Ozone Mapping Spectrometer (TOMS) satellite instrument.61 No available aerosol measure is an ideal proxy for CCN, but many recent studies indicate that AICAMS is better than AOD in this regard.62,63 Note that the above two AI quantities (AICAMS and AITOMS) are not defined equivalently (see Methods) and AITOMS is more sensitive to absorbing aerosols. Precipitation and cloud-area data came from the University of Utah Tropical Rainfall Measuring Mission (TRMM) satellite feature database,64 which reports raining and cloud areas based on the TRMM Microwave Image (TMI) and Visible and Infrared Sensor (VIRS), respectively, where a 273K IR brightness temperature threshold (to detect clouds above the boundary layer) is used to define cloud features. In addition to this, additional meteorological variables used include the number of precipitation features from the TRMM database, and relative humidity, zonal wind and meridional wind obtained from the CAMS.

For our statistical analysis we divided the region between 40o N and 40° S into 5° latitude by 5° longitude boxes (see Supplementary Fig. 1), as performing analysis over smaller regions would help in restricting the influence of regional meteorological conditions on aerosol–cloud-precipitation relationship in a particular box. The boxes between 30° N and 30° S are considered as Tropical and the others mid-latitudinal. We averaged each variable within each box, took the ratios of the box-mean rain and cloud areas for each month, and calculated a slope between log2(AOD) or log2(AICAMS)) or log2(AITOMS)) and area ratio, considering that aerosol effects are expected to be logarithmic in CCN concentration.65 Finally, we took the median (in order to avoid the influence of outliers) slope over all (or ocean-only) boxes as our estimate of the aerosol effect.

Although our model simulations represent an idealised Tropical mean circulation, our mechanism as proposed would not depend on atmospheric properties (e.g. Coriolis parameter) specific to the Tropics. It would presumably require that rain be produced by mixed-phase, organized deep convective cloud systems; these dominate in the Tropics at all times of year, but also dominate in mid-latitudinal regions during warm seasons. Hence, we might expect a similar impact of aerosol on area ratios in the warm-season mid-latitudes and Tropics, but not in mid-latitudes during the cold season.

Results (Fig. 4a–c) show large regional variations in the slopes, indicating that substantial regional meteorological contamination/noise effects remain even when looking at area ratios. Nonetheless, the medians across regions are robustly negative across land, ocean, tropical, and mid-latitude domains (Fig. 4d), ranging from −0.2 to −0.17, roughly consistent with the simulations, regardless of the aerosol data used.

Fig. 4
figure 4

Relationship between aerosol and area ratios from the observations. Map of correlation between rain-to-cloud area ratio and a CAMS-AOD, b AITOMS, and c AICAMS; d overall change in rain-to-cloud area ratio per doubling of CCN. In d the light blue, pink, yellow and green boxes represent model, observed area ratios vs. CAMS-AOD, observed ratios vs. AITOMS and observed ratios vs. AICAMS, respectively. The black lines inside the boxes denote the median over all locations, and the colored boxes show the 5–95% range from bootstrap sampling. Values shown are from the basic model without meteorological predictors (other results summarized in Supplementary Tables 1–3).

Full size image

We perform two statistical significance tests based on the null hypothesis of no systematic aerosol effect (i.e., the predominance of negative r being due to chance). First, we define q = P (r < 0) as the probability that the correlation between aerosol and area ratio within any box is negative. Under the null hypothesis, with only random influences on r, q = 0.5. For n boxes, if X is the number that have r < 0, then its 95% confidence interval is calculated as \(X = nq \pm \sqrt {nq\left( {1 - q} \right)}\). The calculated X for all the three datasets (AOD, AICAMS, AITOMS) were outside this 95% confidence interval, so the null hypothesis is rejected i.e. q cannot be 0.5, and it is concluded that there is a systematic relationship between aerosol and the area ratio.

Second, we computed 5th–95th percentile confidence intervals for the median slope via 2000 bootstrap samples, created from the original sample with replacement and of the same size as the original sample. The resulting statistical uncertainties (Fig. 4d) are small compared to the discrepancies between aerosol products in midlatitudes; in the Tropics there is better agreement. These results imply decreases in the observed rain-to-cloud area ratio per aerosol doubling of 5.3 [6.8, 3.0]% in the Tropics (5 [6.6, 2.4]%, ocean-only); and 6.3 [12.8, 1.5]% in mid-latitudes (11 [17.6, 2.6]%, ocean-only Table 2), where the uncertainty ranges account for both the spread between aerosol products and the statistical uncertainty.

Table 2 Slopes (Percentage per doubling) of observed rain-to-cloud area vs. aerosol changes, showing sensitivity to several data selection criteria.
Full size table

As another test we repeated the analysis for only the (a) coldest three months of the year, (b) hottest three months, and (c) months with precipitation <5 mm/day. Consistent with expectations discussed previously, we obtained no consistent signal in the midlatitudes during winter, and signals were weakened somewhat when stronger rain cases were excluded (Table 2). Thus, our signal comes predominantly from situations expected to be dominated by mixed-phase convective precipitation.

Finally, it is likely that meteorological factors (humidity, wind etc.) explain some of the weak or inconsistent r values; there is also a danger of confounding aerosol and non-aerosol causes. To quantify this, we tested seven regression models (Supplementary Tables 2, 3, and 4) that include up to four meteorological predictors in addition to one aerosol measure. The meteorological predictors considered were mean relative humidity (RH) in the mid-atmosphere (between 700 and 500 hPa), mean zonal (UA), and meridional (VA) wind shear in the lower troposphere (between 950–500 hPa) and the number of precipitation features (P_Area) (Supplementary Tables 2, 3, and 4).

We found that adding meteorological predictors (especially upper-level relative humidity) improves the goodness of fit of the model (r2), confirming that these predictors do affect the cloud/rain area ratio. But too many predictors in a model will lead to overfitting. To check for this we performed cross validation, randomly dividing each set of data into two equal training and test parts. Models with four predictors were not robust, as the r2 over the test dataset declined significantly, whereas the models with fewer predictors were robust. The robust models all support an aerosol effect of comparable magnitude to that shown in the aerosol-only model. This demonstrates that meteorological variables affect the area ratios independently from aerosols. While we cannot rule out that our aerosol signal could be due to some undiscovered, systematic confounding meteorological factor, none of the factors we tested appears to account for it even though they do affect the area ratio. In summary, after performing resampling tests, exclusion tests, and meteorological confounder tests, results remain qualitatively consistent with our hypothesis.

The observed signal is comparable to that from the WRF model. In the POLL2 experiments, the increase from 100 to 800 per cc represents approximately three doublings. To compare simulated results quantitatively with observations, we assume AOD and CCN change proportionally, and hence divide the simulated changes by three (to get a value per doubling of CCN). The perturbed model experiments then show an average 4 ± 3.8% decrease in the cloud-area ratio per doubling, which lies within the spread of the observations (Fig. 4c).

The effect found here has implications for the top-of-atmosphere (TOA) cloud radiative effect (CRE) and hence global climate. With an increase in mid to high-level cloud fraction, there is also observed a TOA net cooling effect in all the perturbed experiments (Supplementary Table 4). The average POLL2, PERT2, and PERT20 experiments respectively show TOA coolings of −1.59 ± 0.46, −2.16 ± 0.83, and −2.66 ± 0.72 W/m2 (between the Equator and 15o latitude), associated with increases of 1.40 ± 0.43, 0.96 ± 0.35, and 0.75 ± 0.34% in cloud fraction (from 0 to 15° latitude). As the proxy heating and CCN experiments produce similar changes in cloud and precipitation, we regard them as equivalent, and compute mean TOA CRE by averaging all the perturbed experiments and dividing by three. This yields a cooling of −0.71 ± 0.25 W/m2 per CCN doubling, associated with an increase of 0.34 ± 0.13% in cloud fraction. Note the quoted means and uncertainties are based only on the model ensemble and do not invoke the observations.

What is the global radiative impact of this mechanism? Here we roughly estimate this, assuming that the albedo response occurs at the magnitude given by WRF throughout the tropics and warm-season midlatitudes, in proportion to local precipitation rate, and that the forcing impact scales with albedo × insolation. Accordingly, we multiply the WRF (tropical) forcing value by 0.63 (the ratio of global to tropical mean rainfall) and again by 0.8 (ratio of global to tropical mean insolation).66,67 This will be an overestimate since we did not exclude winter midlatitudes nor polar regions, but those receive little insolation, and in the opposite direction we ignored the diurnal cycle which would increase the radiative effect since deep convection is somewhat greater during daylight hours. This crude estimate yields a global mean annual forcing of −0.36 W/m2 per doubling (i.e., roughly −0.2 to −0.5 W/m2 given our assumptions). These assumptions are clearly provisional so at best this range represents the rough magnitude expected for this mechanism.

To estimate the consequent present-day anthropogenic indirect forcing requires quantifying the CCN from anthropogenic emissions. Studies have reported an effective increase in CCN due to anthropogenic sources between a factor of 1.2 to 1.7 since preindustrial, or roughly half a doubling, with huge temporal and spatial variation.68,69,70,71 Therefore, to estimate the anthropogenic forcing we multiply the above estimate by 0.5, which gives a best-guess TOA cooling of roughly 0.2 W/m2. This highly provisional estimate is smaller than the direct aerosol radiative effect but is not negligible. An accurate quantification of this effect is needed, which would require more realistic and elaborate modeling of the climate system and the full aerosol distribution.

While there is clear evidence that aerosols affect deep-convective cloud droplet numbers,72 the possible climate consequences of this have been unclear. Our results indicate a new type of indirect aerosol forcing mechanism which acts by altering the dynamics of middle and deep convective clouds. This new forcing mechanism would act separately from previously known direct and indirect aerosol effects. Its non-local behavior also highlights the complexity of aerosol effects on climate, and the difficulty of narrowing uncertainties in total anthropogenic aerosol radiative forcing and hence climate sensitivities to forcings as deduced from the historical climate record.

Methods

WRF setup

The idealized Hadley-cell setup of refs. 43,44 we adopt can equally represent either hemisphere. We chose an idealised model with fixed SST because this is more efficient and useful in answering questions about atmospheric processes. For example, uncoupled atmospheric models are able to reproduce anthropogenic climate change processes simulated by coupled models.73,74,75,76 Past studies have commonly used such models to simulate aerosol indirect effects.77,78

Zonally periodic domains of 3200 (latitude) × 4000 (longitude) km2 and 3200 (latitude) × 100 (longitude) km2 are used for the 20 km and 2 km resolution experiments, respectively. The model top is at 20 km, with 40 or 50 vertical levels used at 20 km or 2 km resolution, respectively. The time step is 60 or 6 sec at 20 or 2 km, respectively.

The WRF double-moment six-class (WDM645) and Morrison Double-Moment four-class (Morr246) microphysical schemes are used for microphysics. Second-order horizontal diffusion is used for turbulent diffusion, and Yonsei University79 and Betts Miller (BM80) schemes are used for the planetary boundary layer and cumulus parameterizations, respectively. For shortwave radiation, the Dudhia47 and Rapid Radiative Transfer Model (RRTM48) parameterization schemes are used, whereas for the longwave, the Rapid Radiative Transfer Model (RRTM48 scheme is used. The simulations include the effects of Earth’s rotation. We note that averages over individual months, while all of the same sign, varied significantly in magnitude such that a three-month average was necessary for stable and accurate estimates.

Aerosol representation in WRF

In this study we employ two different approaches in order to represent aerosol effects in the model. In the first approach, aerosol effects are represented via proxy heating perturbations whereas in the second approach aerosols effects are represented via explicit CCN increase.

A key pathway for convective invigoration by aerosols is thought to be that additional latent heat is released due to the freezing of more cloud droplets lifted above the freezing level.17 In order to test this hypothesis, Morrison and Grabowski40 performed experiments where they increased the local latent heating within updrafts by 20% within an altitude range where freezing rates might be most affected (6–8 km). In order to balance the added latent heating, cooling was applied to downdrafts in such a manner that the net heat input was zero. Here, in our first approach we adopt the same methodology and apply heating and cooling perturbations between the same altitude range, but only in the convective part of the domain i.e. from the Equator to 5° latitude. Perturbations are applied in such a manner that there is net instantaneous balance between the horizontally averaged heating and cooling in each model layer.

The local latent heating QH and cooling perturbations QC, are given by:

$$\begin{array}{*{20}{c}} {{\mathrm{Q}}_{\mathrm{H}} = \left( {{\mathrm{c}} - 1} \right)\,{\mathrm{Q}}_{\mathrm{L}}} & {{\mathrm{if}}\;\;{\mathrm{w}} \,>\, 0\;\;{\mathrm{and}}\;\;{\mathrm{Q}}_{\mathrm{L}} \,>\, 0} \end{array}$$
(1)
$$\begin{array}{*{20}{c}} {{\mathrm{Q}}_{{\mathrm{c}}}\;{=}\;\langle {\mathrm{Q}}_{\mathrm{H}}\rangle \left({{\mathrm{A}}_{\mathrm{H}}/{\mathrm{A}}_{\mathrm{C}}} \right)} & {{\mathrm{if}}\;\;{\mathrm{w}} \,<\, 0\;\;{\mathrm{and}}\;\;{\mathrm{A}}_{\mathrm{C}} \,>\, 0,{\mathrm{A}}_{\mathrm{H}} \,>\, 0} \end{array}$$
(2)

where QL is the local latent heating rate, 〈 〉 is the horizontal average of the heating rate in K/day on a model level weighted by the local pressure thickness of the model level. Updraft velocity is represented by w, and the coefficient c is defined to be 1.2. AH and AC are defined as the horizontal extent of the heating and cooling areas respectively at a given level and time step. These perturbations are applied within the 6–8 km altitude range where most freezing would take place. In the 20-km experiments, the heating rate associated with condensation includes both a resolved component (associated with resolved condensation) and a sub-grid component (associated with the parameterized convection scheme). We applied perturbations to the sum of parameterized and resolved latent heating (PERT20). In the 2-km experiments, there is no convective scheme, so the heating perturbation is applied only to the resolved latent heating (PERT2).

In the second approach, we exploit the complete microphysical schemes in which cloud condensation nuclei (CCN) concentrations (WDM6) or cloud droplet concentrations (Morr2) are explicitly considered. We assume that all the CCN are of the same type and with equal activation potential as considering the effect of aerosol type on convection was beyond the scope of this study. The prognostic water substance variables included in the WDM6 scheme are water vapor, cloud, rain, ice, snow and graupel. In addition to the cloud water mixing ratio, the activated CCN, cloud, and rain number concentrations are prognostic variables. The activated CCN number concentration is predicted from total CCN and diagnosed supersaturation.45 Complete evaporation of cloud droplets is assumed to return the corresponding CCN particles to the total CCN count. In this scheme, an initial fixed value of CCN number concentration 100 cm−3, representing clean conditions, was used in the control and proxy-heating experiments. In the polluted case experiment, the CCN number concentration is set to 800 cm−3 in the convective part of the domain, i.e., between Equator and 5° latitude (the same region where proxy heating perturbations were applied). These CCN concentrations were used using the references from past studies.2,81

The Morr2 microphysical scheme contains kinetic equations for the mixing ratio and number concentration of each hydrometeor species in the model: cloud droplets, cloud ice, snow, and rain.46 The scheme considers prognostic supersaturation, droplet activation, heterogeneous and homogeneous ice nucleation, and the spectral index (width) of the droplet size distribution. Similar to the WDM6 scheme, we in this scheme set an initial value for cloud droplet concentration as 100 cm−3 for control and proxy heating experiments and for the polluted experiment increase this value to 800 cm−3 in the convective part of the domain. A second 2-km experiment (POLL2) was done by increasing either the CCN or the cloud droplet concentration (based on the microphysical scheme used) to 800 cm−3. Name and details of all the simulations are listed and defined in Table 1.

Observations

We also analyse 10 years of monthly mean precipitation, cloud and aerosol data in the Tropics and mid-latitudes. As accurate aerosol retrievals from space are currently not possible underneath clouds, and may be affected by clouds even where retrievals are attempted, we examine monthly mean observed aerosol optical depth (AOD) data at 550 nm from the Monitoring Atmospheric Composition and Climate (MACC) II dataset which is now a part of Copernicus Atmospheric Monitoring Service (CAMS) global atmospherics composition data.58,59,60 CAMS use a four-dimensional variational data assimilation technique to combine satellite observations with chemistry-aerosol modelling to obtain a gridded continuous representation of the mass mixing ratios of atmospheric gases and aerosol. The global model and data assimilation system of CAMS is based on the ECMWF’s integrated forecast system, and the atmospheric chemical system is represented by the Model of Ozone and Related Chemical Tracers (MOZART)82 chemical transport model.

Recent studies indicate that Aerosol Index (AI) is a better proxy for CCN than is AOD, and that the variability in AOD only explains a small fraction of the CCN variance.62,63 For these reasons, AICAMS is used as a second product for the aerosol observational analysis. AICAMS is calculated using AOD data at 550 and 865 nm from the CAMS dataset. The following formula is used to calculate AICAMS:

$${\mathrm{AI}}_{{\mathrm{CAMS}}} = {\mathrm{AOD}}_{550} \ast \left( { - \frac{{\ln \left( {\frac{{{\mathrm{AOD}}_{550}}}{{{\mathrm{AOD}}_{865}}}} \right)}}{{\ln \left( {\frac{{{\mathrm{\lambda }}_{550}}}{{{\mathrm{\lambda }}_{865}}}} \right)}}} \right)$$
(3)

We also use aerosol index (AI) from Total Ozone Mapping Spectrometer (TOMS) instrument61 as a third aerosol product. AITOMS is a measure of how much the wavelength dependence of backscattered UV radiation from an atmosphere containing aerosols (Mie scattering, Rayleigh scattering, and absorption) differs from that of a pure molecular atmosphere (pure Rayleigh scattering). Quantitatively, AITOMS is calculated from the ratio of measured to calculated 360 nm TOMS radiances so it is not equivalent to AICAMS.

The reason for using TOMS AI in this study is because TOMS AI uses UV wavelengths to determine aerosol properties (the other two aerosol products CAMS AOD and AI are retrieved in the infrared channel). Measurement in the UV range from space is a useful method for detecting absorbing aerosols (smoke and dust), which cannot be effectively discriminated in the visible range. Measurement in the UV range has another advantage of being able to detect aerosols above backgrounds that are bright in the visible range, such as low clouds (and snow although not relevant to our study). Several studies in past have shown that AITOMS is proportional to AOD measured using different techniques over various parts of the world.83,84

Monthly mean precipitation and cloud-area data are taken from the University of Utah Tropical Rainfall Measuring Mission (TRMM) satellite feature database,64 which reports raining and cloud areas based on the TRMM Microwave Image (TMI) and Visible and Infrared Sensor (VIRS), respectively. This dataset is created by temporally and spatially collocating the measurements from these instruments and then by grouping contiguous pixels using surface rain, cold infrared or microwave brightness temperature criteria. We use data from the TMI Precipitation Feature (TPF) and Cloud features at 273 K (C273F). TPF’s are defined as the TMI pixels with surface rainfall greater than 1 mm, whereas the C273F features are defined as the VIRS pixels with brightness temperature less than 273 K. We also re-did the latter analysis with 235 K(C235F) and 210 K (C210F) and found similar results. We also recalculated the observed result on a variety of grid sizes and by restricting the analysis to a particular season (summer) and did not find any significant dependence on any of the two parameters. However, we obtained no consistent signal in the midlatitudes during winter, and signals were weakened somewhat when stronger rain cases were excluded (results not shown).

In order to investigate aerosol effect on cloud and precipitation in observations, we divide the area between 40° S and 40° N into 1152 square boxes of equal size (Supplementary Fig. 1). Each box is 5° in latitude and 5° in longitude. Tests with larger box sizes (30° in latitude and 30° in longitude) gave similar results (not shown). Each box containing only the ocean pixels are considered as purely ocean box.

For quality control, we capped monthly mean AOD values at 0.4 since values above this could be spurious, and in any case, may excessively influence linear regressions. Because seasonal cycles of cloud and rain area may be strongly influenced by non-aerosol meteorological factors, we removed the seasonal cycle and trend from each quantity (AOD, AICAMS, AITOMS, rain area, cloud area), by first removing the long-term linear time trend for each calendar month and then subtracting the difference between the mean of that month and the annual mean.

Over each box the correlation and slope parameter are calculated between the log2 of aerosol optical depth (AOD) or log2 of aerosol index (AICAMS) or log2 of aerosol index (AITOMS) and area ratio between rain and cloud. In order to compare the observed area ratios to the model, we calculate the area ratio between the precipitation area fraction and cloud area fraction for all the experiments. Here, precipitation area fraction is calculated by averaging the grid points with precipitation greater than 1 mm, and an average cloud fraction is calculated for the altitudes excluding the boundary layer (similar to the observation) by simply averaging the model cloud-fraction variable.