Nonlinearity of the cloud response postpones climate penalty of mitigating air pollution in polluted regions

Aerosol–cloud interactions contribute substantially to uncertainties in anthropogenic forcing, in which the sensitivity of cloud droplet number concentration (Nd) to aerosol plays a central role. Here we use satellite observations to show that the aerosol–Nd relation (in log–log space) is not linear as commonly assumed. Instead, the Nd sensitivity decreases at large aerosol concentrations due to the transition from aerosol-limited to updraft-limited regime, making the widely used linear method problematic. A sigmoidal transition is shown to adequately fit the data. When using this revised relationship, the additional warming that arises from air pollution mitigation is delayed by two to three decades in heavily polluted locations, compared to the linear relationship. This cloud-mediated climate penalty will manifest markedly starting around 2025 in China and 2050 in India after applying the strongest air quality policy, underlining the urgency of mitigating greenhouse gas emissions. Cloud droplet number concentrations are often assumed to depend linearly on atmospheric aerosols (in log–log space). Here the authors show that this relationship is instead sigmoid, which delays additional warming due to air pollution mitigation by 20–30 years in heavily polluted regions.

Aerosol particles affect Earth's energy balance through exerting effective radiative forcings (ERF aer ) via both aerosol-radiation interactions and aerosol-cloud interactions 1 . Specifically, aerosols can alter radiative balance directly by scattering and absorbing sunlight (ERF ari ) and also indirectly by acting as cloud condensation nuclei (CCN), in turn modifying cloud properties and precipitation (ERF aci ; [2][3][4] ). ERF aci has been known as the most uncertain component in anthropogenic forcing, in which the sensitivity of cloud droplet number concentration (N d ) to aerosol plays a fundamental role 5,6 .
The latest report by the Intergovernmental Panel on Climate Change provided a best estimate (5 to 95% confidence interval) of the total ERF aer at −1.1 (−1.7 to -0.4) W m −2 (ref. 1). The resulting cooling due to aerosols in turn counteracts about half of the ERF due to CO 2 from 1750 to 2019 7 . To mitigate risks of air pollution, most industrial regions have implemented air quality regulations, for example, emissions reductions in Europe (EU) and North America (NA) since the mid-1980s 8,9 and in China since 2006-2010 10,11 . This could bring additional 'unmasked' warming due to mitigated cooling from aerosols, leading to what has been called a climate 'penalty' 12 . Modelling studies showed that environmental policies induced warming of 0.25 ± 0.12 °C globally from 2015 to 2055 13 16 . Interestingly, more recent evidences from multi-satellite observations together with models also show a clear increase in negative ERF aer by 0.1 to 0.3 W m −2 from 2001 to 2019, confirming that the climate penalty has now clearly started to happen 17 .
To confidently predict the future climate penalty, a better observational constraint on ERF aci is vital because ERF aci contributes most to the uncertainty of ERF aer 18,19 . In this regard, long-term trends in aerosol loading and cloud variables observed by satellites may be particularly helpful. There have been some attempts to link long-term trends in cloud variables with emissions change. Along with the emissions reduction in EU and the United States, N d declines were observed correspondingly [20][21][22][23] , while trends in liquid water path and cloud fraction (CF) were not always consistent with emissions change 22,23 due to their large natural variability compared to the variability driven by N d . Consistent trends in N d and aerosols were also reported over India (IN) 21,23 and the East China (EC) Sea 22,24,25 but were not evident over China 21,23 . Those findings support theoretical hypotheses 2,3 but only in terms of qualitative consistency in signs across regions 17 .

Implications for N d sensitivity
It is important to note that S AI is calculated based on decadal variations of regional means, in turn limiting aerosol variability, which ties S AI to a certain aerosol condition. There is a strong negative dependence of S AI on regional mean AI (Fig. 2a). In very polluted regions, which is generally updraft limited, S AI is close to 0, while for relatively clean oceans, which is more likely aerosol limited, it can be almost 1, that is, the maximal value we would expect. This reflects a transition from aerosol-limited to updraft-limited regime with increasing aerosol 28 .
However, what we have been doing is to compute the N d sensitivity as a fixed value with aerosol varying for a specific regime. Here the regimes are mainly used to constrain the co-variability of aerosol type and meteorological background 29 . The problem is, no matter how accurately those regimes can be constrained, as long as S is kept as a constant with the aerosol regime, the aerosol-updraft transition is naturally fixed. But according to the finding here, it is apparently not true. By continuously increasing aerosol, we will eventually go to the updraft-limited regime, for example, the insensitive N d seen in EC and IN. A reduced S with aerosol increasing or a saturated N d at high aerosol loading even in log-log scale is thus theoretically expected.
By considering the linear aerosol-N d relation in logarithmic space, the OLS method can capture the tendency towards saturation as reflected by a power-law relation in regular space, but it is insufficient to capture the absolute saturated N d under high aerosol loading unless the sensitivity is set to zero. Although this nonlinear ln aerosol-ln N d relation has been observed in previous studies that used different CCN proxies [30][31][32] , the OLS method was still used to fit the relation. As shown in Fig. 2b (the upper panel), the AI-N d joint histogram constructed from non-precipitating clouds over global oceans is well characterized by an S shape. Note that the flat N d under clean conditions is not physically meaningful and instead is caused by the problematic aerosol retrieval 33 . The conventional OLS method is unable to capture the 'shallow' relation on both low and high ends (Fig. 2b). Besides this, the OLS method is sensitive to the probability distribution of aerosol as it puts more weight to commonly occurring aerosol values. This explains why the low sensitivity was often observed over the pristine Southern Ocean 31,34 , where a strong N d response is expected instead.
To overcome the issues raised by the OLS method, we propose using a sigmoidal function to fit the ln AI-ln N d relation. Here the median N d in each AI bin is taken as the input of regression, ensuring that the equal weight is put to each AI. It is exciting to see that the sigmoid method perfectly captures the N d variation. The peak S AI of sigmoid method is 1 at AI = 0.16, which is a demarcation between the aerosol-limited and transition regime, assuming first-order unchanged updraft with aerosol. In the aerosol-limited regime, the reduced S AI was reported as an artefact of the aerosol-retrieval issue 33 , thus the plausible value should be 1, given the nature of the aerosol-limited regime 28 . Then, the plausible N d prediction is inferred accordingly (dotted cyan/ grey line in Fig. 2b). This corrected sigmoid curve will be used in the rest of the paper. Surprisingly, in the aerosol-updraft transition regime (AI > 0.16), the S AI from daily variations are overall well consistent with ones derived from decadal trends-decreasing sensitivity as aerosol increases, except for EU, where the AI variation cannot reflect the CCN decadal trend. The two lines of evidence confirm the nonlinearity of ln AI-ln N d relation.
Capturing this nonlinearity is crucially important when predicting N d evolution and/or estimating radiative forcing due to aerosol-cloud interactions (RF aci ). According to the observations, N d starts to be largely saturated around AI = 0.35, but the OLS method gives a continuously increasing N d (Fig. 2c). This means that the OLS method tends to overestimate N d change from pre-industrial (PI) to present day (PD) Here our study provides quantitative estimates of the N d sensitivity from decadal variations via filtering out high-frequency signals. Together with the analyses based on daily variations of aerosol and N d , this study comes up with a robust depiction of the N d sensitivity that the aerosol-N d relation is nonlinear with a saturated N d under polluted conditions even in log-log scale. This, in turn, suggests that the commonly used ordinary least-squares (OLS) linear fit on ln aerosol-ln N d (OLS method) might be problematic to accurately predict the N d change on the basis of aerosol changes. To tackle this, a sigmoidal fit is proposed as a better alternative and then used to describe the past and predict the timing of cloud-mediated climate penalty with future emissions reduction. Furthermore, we show that current climate models have a limited ability to reproduce the past N d variation, mainly owing to the insufficient representation of the saturation effect. Our findings highlight the importance of improving the estimation of aerosol-cloud interactions from both observational and modelling sides, especially in terms of the nonlinearity of N d sensitivity.

Observed decadal trends
Unlike qualitatively comparing the consistency between trends in aerosol and cloud variables 17 , here we compute the N d sensitivity from decadal variations in a quantitative sense. To isolate the signal that is solely relevant to decadal emissions change, a de-seasonalization and a locally weighted scatterplot smoothing were used (Methods).
Decadal variations of normalized anthropogenic aerosol emissions, aerosol index (AI) and N d from 2001 to 2020 over major industrial regions ( Supplementary Fig. 1) are shown in Fig. 1. Continuous declines in sulfur (SO 2 ) emissions are seen with changes, obtained as linear trends, of -5% per year over EU and -9% per year over NA. In contrast, the SO 2 emissions over IN increase throughout the whole period by 4% per year; in EC it has been increasing until around 2008 (11% per year) and decreasing thereafter (-9% per year). Black carbon (BC) and organic carbon (OC) emissions have similar trends as SO 2 but are less pronounced.
Due to the short lifetime of aerosols, the variations of AI-a proxy of fine-mode aerosol concentration-are tightly following the emissions changes ( Fig. 1e-l). However, N d changes are not always consistent with AI, which is mostly aerosol-environment dependent. For relatively clean regions (AI < 0.25), the N d variations in general well fit aerosol changes; over the rather clean ocean of NAO, a 1:1 N d -to-AI sensitivity , unitless) is even observed; NA and EUO also have high S AI of 0.49 and 0.57, respectively. An exception is EU, where an increasing AI is seen until 2010, presumably due to the increasing anthropogenic OC emissions (Fig. 1a) and continuous wildfire events 26 . In this case, the AI variation is linked more to increasing non-CCN aerosols (BC and OC) rather than the expected declining CCN (sulfate), hence hampering the S AI estimate. In contrast, for polluted regions, N d is less responsive to aerosol changes; over IN, a 53% increase in AI from 2001 to 2020 results only in 13% enhancement in N d . This is particularly severe over EC, where N d does not change at all with such a strong switch in aerosol trends (135% increase before 2008 and then −95% decline afterwards). Similar behaviours also appear in their adjacent oceans ( Fig. 1k-l). The corresponding S AI for those polluted regions are lower than 0.15, much smaller than the lower bound (0.3) given by ref. 6.
In addition to aerosols, clouds also respond to global warming and internal climate variability on interannual and decadal timescales 1 , which were reported to largely affect observed trends in CF and liquid water content 27 . To confirm if the N d changes are solely a manifestation of aerosol impacts so that the decadal variations can be confidently used to calculate S AI , we examine Coupled Model Intercomparison Project Phase 6 (CMIP6) multi-model simulations with external radiative forcing only from greenhouse gases (CMIP6-GHG) and only from natural solar variations and volcanic aerosol (CMIP6-NAT; Supplementary Table 2 includes model information). The negligible Article https://doi.org/10.1038/s41558-023-01775-5 over very polluted regions, while underestimating it over moderately polluted regions, which would bias the spatial patterns of N d change 30 and hence RF aci [34][35][36] .

Past and future trends in N d and RF aci
The OLS method has been often employed to estimate the PD-PI N d change (ΔN d ) and RF aci , however, it is difficult to evaluate its accuracy due to the impossibility of observing the PI N d . In this regard, satellite-observed N d trends over past decades provide a great opportunity to evaluate the predictability of OLS and sigmoid methods. Moreover, the future trends are also of particular interest. As recently documented by ref. 17, on a global scale, the reduction in aerosol emissions since 2000 brought additional warming radiative effects through both direct and indirect effects. Nevertheless, at a regional scale over China, we do not actually see the climate penalty mediated by the indirect effects, as N d nearly did not change because of the saturation effect at high aerosol. It is expected that with the future emissions reduction, the climate penalty will eventually manifest at increasing rates. Thus, what is interesting now is when this will happen.
The agreement between the AI predicted by the historical SO 2 emissions and observed by MODerate Resolution Imaging Spectroradiometer (MODIS) lends credibility to future AI predictions (Fig. 3). From these, in turn, the transient changes in N d can be inferred. The OLS-diagnosed N d strongly scales with aerosol changes (Fig. 3b,e) and thereby is unable to reproduce the saturated (slightly increasing) N d over EC (IN) in the past. Also, for the SSP2-4.5 and SSP3-7.0 pathways in IN, where the AI remains higher than the critical AI ( Fig. 3d and Extended Data Fig. 1), the updraft-limited behaviours are expected, however, the OLS-diagnosed N d is still highly sensitive to aerosol.
It is interesting to see that the sigmoid method performs very well in reproducing the past trends over each region-not only the saturated N d in EC and IN (Fig. 3b,e) but also strong declines in NA and EU (Extended Data Fig. 2 (Fig. 3a,d).
The N d changes are consequently reflected in RF aci evolutions. The strong changes in aerosol emissions in EC over the past two decades did not exert any appreciable RF aci (Fig. 3c). Similarly in IN, the strong sulfur emissions resulted in a regional mean RF aci of merely -0.45 W m −2 by 2020 relative to 2001 (Fig. 3f). However, for the regions that already transited out of the updraft-limited regime due to air quality policy implementions (including NA and EU in Extended Data Fig. 2), aerosol reductions efficiently led to stronger warming with RF aci of 1.37 W m −2 and 0.91 W m −2 by 2020, respectively.
According to the sigmoid-based predictions, if the strongest air quality policy (SSP1-2.6) will be applied, the rapid warming via RF aci are estimated to occur around 2025 in EC and 2050 in IN (the same time as N d falling), with ~17 and ~30 years delayed compared to the conventional OLS predictions. Note that we considered aerosol-limited S AI = 1 in the sigmoidal prediction. Sensitivity experiments further show that the diverse treatments of aerosol-limited S AI lead to spreads in the RF aci of ~1-2 W m −2 at 2100 ( Supplementary Fig. 2), but this does not affect the prediction of timing when the climate penalty will occur. Similarly, the warming in EC will be postponed by ~22 years under SSP2-4.5; but under SSP3-7.0, no additional forcing will be imposed to EC and IN owing to the saturated N d . Beyond this delayed warming in polluted regions, the sigmoidal predictions reveal an accelerated warming in clean regions (NA and EU), with nearly double increase rates in RF aci relative to the OLS results (Extended Data Fig. 2). As cloud adjustments approximately scale with RF aci 6,37 , the predicted evolutions are relevant for ERF aci , too.

Performance of CMIP6 models in reproducing N d trends
Modelling studies demonstrated enhanced future warming and precipitation from aerosol reductions 14,16 . However, the lack of a proper representation of the observed nonlinearity of the N d response to aerosols (in log-log space) in climate models may imply they overestimate the warming from aerosol reductions in polluted regions that are still in the updraft-limited regime. For example, the climate models suggested a regional ERF aer of ~3 W m −2 induced by emissions reductions during 2006-2017 over EC, along with strong enhancements of surface air temperature (~0.3 °C) and precipitation (~70 mm per year) (ref. 16). Such strong climate effects are unlikely to occur in the absence of aerosol indirect forcing, as it contributes the largest part of ERF aer 1,6 . Therefore, it is essential to evaluate the performance of current climate models in reproducing observed N d trends. We see that none of the CMIP6 models can reproduce the observed unchanged N d trends in EC (Fig. 4). Though some models perform better for the regions where N d is not completely saturated (IN, NA and EU), the majority of models still tend to systematically overestimate the relative changes compared to observations and sigmoid-based predictions. Under the future polluted scenario (SSP3-7.0), the saturated N d is expected over EC and IN (Supplementary Fig. 3), according to our analysis, but the CMIP6 models instead show fairly strong N d changes.
When focusing on models that have the same dynamical core, turbulent diffusion and cloud microphysics but only differ in activation treatments (green shaded areas), it is clear that the detailed activation scheme consistently performs better than the OLS empirical parameterization. The detailed activation scheme accounts for the depletion of the maximum supersaturation under high aerosol loading 38 , thus allowing the simulation of the reduced N d sensitivity with more aerosol abundance 39,40 to some extent. This, however, is still insufficient to fully capture the nonlinearity as observed by satellites.

Discussion
Our results based on both decadal and daily variations demonstrate that the aerosol-N d relation even in log-log space is not linear as commonly assumed. Specifically, the N d sensitivity is shown to be negatively correlated with aerosol loading in a region, indicating a transition from aerosol-limited to updraft-limited regime with aerosol increasing. This finding suggests that the widely used OLS method and corresponding RF aci estimate 30,34-36 are problematic because they assume a constant N d sensitivity despite the mean aerosol varying.
As a step forward, a sigmoidal function is proposed and shown to fit the S-shaped ln AI-ln N d relation seen by satellite statistics. The optimized sigmoid method overcomes two major issues: (1) the nonlinearity of the ln aerosol-ln N d relation and (2) the poor retrieval capability under clean conditions, which drives part of the uncertainty in aerosol-cloud radiative forcing 39,41 . As a result, the sigmoid method performs very well in reproducing the past N d trends-not only the saturated N d in EC and IN but also strong declines in NA and EU, lending credibility to the near-term N d and RF aci predictions.
Although the cloud-mediated climate penalty of air quality improvements was not observed in EC over past decades because of the saturation effect, it is expected to manifest eventually with future mitigation of air pollution (scenario SSP1-2.6). According to the sigmoid-based predictions, the rapid warming via RF aci is estimated to start occurring from around 2025 in EC and from 2050 in IN. In relatively clean regions (NA and EU), the sigmoid-based predictions show much larger increasing rates in RF aci compared to the conventional OLS method. This highlights the urgency of mitigating greenhouse gas emissions to avoid strong temperature increase rates.
ERF aci can be faithfully projected only if the historical N d evolution is faithfully simulated. However, we show that none of the CMIP6 models can reproduce the saturated N d over EC satisfactorily. Generally, the majority of models tend to overestimate the relative changes of N d . Although the detailed activation scheme is better than the OLS empirical parameterization, it is still insufficient. Failure to simulate the  Fig. 2b, below which the updraft-limited regime is transiting towards the aerosol-limited regime, that is, the aerosol effect starts to appear. Here the critical AI is defined as the AI value when S AI equals to 0.2, in which case the aerosol effect on N d is assumed to be negligible.
Article https://doi.org/10.1038/s41558-023-01775-5 saturation effect might explain why climate models suggest stronger cooling from aerosols in the 1970s 42 . The findings here emphasize that further improvements to the current activation scheme are needed to achieve better ERF aci projections. A few potential caveats to the sigmoid predictions should be noted. Considering the large uncertainties in retrievals over land, the sigmoidal curve was inferred from well-filtered single-layer liquid clouds over the global ocean. The cost of retrieval quality is to discard a portion of clouds, which is acceptable compared to the large confounding effects caused by retrieval errors 32,43 . Additionally, the differing backgrounds of meteorological conditions and aerosol types among regions might lead to a difference in the threshold of AI for the updraft-limited regime that cannot be reflected in the current sigmoidal fit. However, this effect appears small (Supplementary Fig. 4). The reproduction of historical N d trends over land largely supports the utility of this sigmoidal relationship over land. However, to thoroughly fix the above issues, using improved CCN and cloud retrievals over land 31 and constraining updraft regimes 32 would be useful ways forward. Additionally, to eliminate the confounding effect of precipitation on the aerosol-N d relation, only non-precipitating clouds were analysed, so the sigmoid-based N d and RF aci projections did not account for the impact of changing precipitation under global warming. Although those caveats might add uncertainties or alter the predicted year when the climate penalty will occur, it does not hamper the robustness of the optimized sigmoidal fit as a promising alternative for the conventional OLS method that has been shown to be problematic.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-023-01775-5. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.

Aerosol and N d observations
Level 3 1° × 1° products of the MODIS Collection 6.1 (MOD08 and MYD08) provide aerosol optical depth (AOD) and Ångström exponent (AE) retrieved at several wavelengths globally 47,48 . In this study, the AOD is from Dark Target and Deep Blue combined product. Compared to AOD, the AI is considered a better proxy for CCN because it is more weighted towards smaller particles 49,50 . The AI is derived from AOD and Ångström exponent (AI = AOD × AE), in which AE is calculated from AOD at wavelengths of 460 and 660 nm. To minimize the influence of heterogeneous aerosol distribution, we discarded aerosol retrievals where the relative standard deviation is larger than unity 51 . Cloud droplet effective radius (CER) and cloud optical depth from the MODIS Level 2 products (MOD06 and MYD06; 52 ) are used to calculate N d based on the adiabatic approximation 53 . The retrievals at 3.7 μm are used as expected to produce more accurate CER retrievals in inhomogeneous conditions 54 . The following sampling strategy is applied to obtain confident retrievals, that is, only including single-layer liquid clouds with (1) cloud-top temperature higher than 268 K, (2) CER > 4 μm and cloud optical depth > 4 to reduce the uncertainty in weak optical signal conditions 53 , (3) CF at 5 km resolution > 0.9 and sub-pixel inhomogeneity index (cloud_mask_SPI) < 30 to reduce the retrieval errors in broken cloud inhomogeneous conditions 54 and (4) solar zenith angle < 65° and sensor zenith angle < 41.4° to minimize the uncertainties raised by cloud 3D effects and multiple scattering 55 .
The 2B-CLDCLASS product 45 from CloudSat radar for the year 2008 is employed to identify precipitation. The precipitation data at a 1.4 × 2.5 km 2 resolution are matched to the nearest MYD06 5 × 5 km 2 pixels to obtain co-located precipitation flag and N d . To exclude the possible confounding effect of aerosol-precipitation interactions on N d sensitivity estimate 32 , only non-precipitating clouds (with the flag of 'no precipitation') are analysed to generate the AI-N d joint histogram, which is then used for linear and sigmoidal fits between AI and N d .
For the analysis on long-term trends, aerosol and cloud observations in the morning (MOD08 and MOD06 from Terra for 2001-2020) and afternoon (MYD08 and MYD06 from Aqua for 2002-2020) are averaged to represent daily-mean conditions. However, for the AI-N d joint histogram (Fig. 2b) that is used for the predictions, the simultaneous precipitation observation from CloudSat is necessary for excluding precipitating clouds, thus only observations from the A-Train constellation of satellites in the afternoon (2B-CLDCLASS, MYD08 and MYD06 from CloudSat and Aqua for 2008) are used. We show that one-year data are sufficient to capture the sigmoidal relationship with only a slight change in the statistics compared to the analysis on a longer period ( Supplementary Fig. 5).

Emissions data
For the historical period (2001-2019), anthropogenic emissions of sulfur dioxide (SO 2 ), OC and BC are obtained from the newest version of the Community Emissions Data System (CEDSv_2021_04_21; 44 ). An important update compared to the old version of CEDS used in CMIP6 56 is a correction on aerosol emissions trends in China. Comparing to somewhat unchanged SO 2 in the old CEDS from 2006 onwards, the newest version is able to capture strong declines in reality 44 .
For the future period (2020-2100), anthropogenic aerosol emissions from different shared socio-economic pathway (SSP) scenarios are used 57 , including the low-emission SSP1-2.6 ('sustainability'), intermediate SSP2-4.5 ('middle of the road') and high-emissions SSP3-7.0 ('regional rivalry') scenarios. Notably, the officially released SSP emissions are harmonized to old CEDS emissions (CEDS-v2017-05- 18) in the year 2015 58 , which causes a discontinuity in the emissions between the newest CEDS data and the future scenarios. To perform a consistent analysis, we thus utilize the software Aneris 59 to harmonize regional-integrated SSP emissions trajectories to the newest CEDS; the methodology of Aneris are described by ref. 60.

Climate model data
The N d simulated by climate models are taken from CMIP6 'historical', 'hist-GHG' and 'hist-NAT' experiments for the period 2001-2014 61,62 and then extended by the Scenario Model Intercomparison Project (Sce-narioMIP) 'sspxxx', 'ssp245-GHG' and 'ssp245-NAT' experiments until 2100 63 . Here 'sspxxx' means different SSP scenarios (SSP1-2.6, SSP2-4.5 and SSP3-7.0). The 'historical' and 'sspxxx' experiments are driven by all forcings, whereas 'hist-GHG' and 'ssp245-GHG' experiments only by greenhouse gas forcing and 'hist-NAT' and 'ssp245-NAT' experiments only by natural forcing. Due to the difficulty in determining the top height of liquid clouds from CMIP6 outputs, alternatively, we make use of the maximal N d in a vertical atmospheric column to compare with the cloud-top N d from satellite 64 . To investigate the impact of aerosol activation treatments on the simulated trend in N d , the models are further divided into two subsets that using: (1) empirical approach that directly links N d to aerosol concentration and (2) detailed activation schemes, respectively. The details of CMIP6 models analysed in this study are given in Supplementary Table 2.

Isolating decadal signal from high-frequency noise
The systematic change in anthropogenic aerosol emissions occurred over a quasi-decadal timescale. To identify the aerosol-cloud signals solely related to the decadal emissions change, possible noises (that is, high-frequency variability) need to be excluded, specifically (1) the seasonality caused by co-varied meteorological fields and emissions and (2) the occasional fluctuation due to natural emissions, for example, volcano eruptions, wild fires and dust storms. To eliminate (1), we first aggregate the filtered daily AI and N d into monthly means at each grid point and then remove the monthly cycle from the monthly series by subtracting the difference value between monthly climatology and full-period mean. After regional averaging, we obtain a de-seasonalized monthly series with sample size of 240 (12 months × 20 years) for each grid point. Furthermore, the locally weighted scatterplot smoothing (LOWESS) method is applied to the de-seasonalized series to exclude (2). At each data point, LOWESS conducts a weighted regression within a prescribed span width that determines the number of nearby data used in each local fit 65 . The choice of span width depends on the timescales of interest; here a seven-year time window is chosen to smooth out high-frequency noise. It is illustrated in Supplementary Fig. 6 that the derived N d sensitivities only change slightly if a different time window is chosen (for example, 4 and 14 years).

Regression models for capturing aerosol-N d relation
Two regression models are employed to fit aerosol-N d relations in this study. The first model is the widely used OLS linear regression where β is the sensitivity of N d to CCN proxy, α. The OLS method describes a power-law relationship between aerosol and N d , which has been broadly observed by in situ aircraft measurements, groundand satellite-based remote sensing, however, mostly over clean to moderately polluted regions 66 . The second model is a sigmoidal regression ln N d = a 1 1 + e a 2 ln α+a 3 + a 4 (2) which has a characteristic S-shaped curve, thereby is able to capture the markedly saturated N d at high aerosol as seen in satellite observations [30][31][32]41 .

N d prediction
The key idea of predicting N d from an observational perspective is to bridge two connections: (1) anthropogenic aerosol emissions to aerosol Nature Climate Change Article https://doi.org/10.1038/s41558-023-01775-5 and (2) aerosol to N d , by which future aerosol and N d can be inferred from different SSP scenarios. It has been documented that sulfate (thus CCN) concentrations are less sensitive to SO 2 precursor emissions under polluted conditions because of the limited availability of oxidants in the atmosphere and the intensive condensation sink of nucleating H 2 SO 4 vapour onto existing aerosols 67,68 . Therefore, we use a power-law relationship to describe this nonlinearity. Because non-sulfate aerosols have been found to contribute only minimally to the long-term trends in AOD 69,70 and therefore to AI, only SO 2 emissions are used to predict AI (treated as CCN proxy). In this case, a linear fit on a logarithmic scale is applied: For each individual region, the regression models are trained by annual-mean CEDS emissions and MODIS aerosol (2001-2019) and then applied on SSP emissions to predict AI in different future scenarios (2020-2100). Using the CMIP6 model output as a direct source for obtaining aerosol in the near-term future is also an alternative approach; however, considering the incorrect emissions used in CMIP6 models ('Emissions data' section), we choose not to do so here.
Regarding the aerosol-to-N d relations, both sigmoidal and OLS methods (as detailed earlier) are utilized but with the preference being put on the former. Data used to construct the regressions are one-year (2008) MODIS 5 × 5 km 2 daily N d and AI (projected from 1° × 1° without altering the original values) observations over global oceans between 60° S and 60° N, which are carefully filtered for retrieval errors and precipitation to reduce uncertainty. For the OLS method, the data are grouped into 20 AI bins, with each bin having an equal number of samples. The median values of N d and AI in these bins are then used in OLS regression as in refs. 31 and 32 . As for the sigmoid fit, an AI-N d joint histogram is constructed first with fixed intervals in a logarithmic scale; then the medians of N d in each AI interval together with central AI are considered as inputs for the regression.
Given the AI-N d relationship over land is quite unreliable due to the retrieval issues; for example, even the 'anti-Twomey effect' can be seen from satellites in most continental regions 43,71,72 , here we apply the aforementioned AI-N d relationship over global oceans to different continental regions. Note that using a single global relation would cause differences between observed and predicted N d for individual regions. This is probably attributable to the varying (1) cloud adiabaticity, (2) surface reflectance and (3) meteorological background across regions, in which (1) and (2) appear to impact N d retrieval, while (3) reflects the capability of droplet activation. It is also interesting to note that the difference between observed and predicted N d was found to be close to a certain constant for each region and independent of long-term temporal evolution 30 . Building on this idea, we add a correction constant-that is, defined as the difference in multi-year averaged N d between observations and predictions-to predicted N d for each individual region to overcome this issue.

Radiative forcing
From predicted N d , the change in ln N d relative to 2001 (Δ ln N d ) is computed. On the basis of Δ ln N d , the change in cloud albedo is derived using the Twomey formula 73 ; then radiative forcing due aerosol-cloud interactions (the Twomey effect; RF aci ) relative to 2001 can be thereby computed: where F ↓ and α cld are mean incoming solar radiation flux and the cloud albedo from the Clouds and the Earth's Radiant Energy System and f liq is the CF of liquid clouds from the MODIS Level 3 product.