Strong control of effective radiative forcing by the spatial pattern of absorbing aerosol

Over the coming decades, it is expected that the spatial pattern of anthropogenic aerosol will change dramatically and the global aerosol composition will become relatively more absorbing. Yet, the climatic impact of this evolving spatial pattern of absorbing aerosol has received relatively little attention, in particular its impact on global-mean effective radiative forcing. Here, using model experiments, we show that the effective radiative forcing from absorbing aerosol varies strongly depending on their location, driven by rapid adjustments of clouds and circulation. Our experiments generate positive effective radiative forcing in response to aerosol absorption throughout the midlatitudes and most of the tropical regions, and a strong ‘hot spot’ of negative effective radiative forcing in response to aerosol absorption over the tropical Western Pacific. Further, these diverse responses can be robustly attributed to changes in atmospheric dynamics and highlight the importance of this ‘aerosol pattern effect’ for transient forcing from regional biomass-burning aerosol. Changes in the spatial pattern of aerosol could influence climate through effects on radiative forcing. Model experiments show that while aerosol absorption in the midlatitudes and regions of tropical descent can warm the planet, aerosol absorption in regions of tropical ascent can cool the planet.

A nthropogenic aerosol plays a substantial role in determining the influence of human activities on climate 1-3 . Due to their short lifetimes, the spatial pattern of anthropogenic aerosol is highly heterogeneous, and this has been shown to impact the climate through changes in atmospheric circulation 4,5 , tropical precipitation [5][6][7][8] and sea-surface temperature (SST) patterns 6,9,10 , with implications for the transient climate response 11 .
Over the coming decades, emissions of anthropogenic aerosol are expected to change dramatically both in their spatial pattern and composition 3 , particularly over South-East Asia and Central Africa 2,3 . These trends are further complicated by the fact that policy measures often target sources of scattering and absorbing aerosol differently, such that even with no net change in aerosol optical depth (AOD), there could be climatic impacts from changes in aerosol composition 12 . Such a trend is already being seen globally; while the global anthropogenic AOD has been roughly constant since the 1980s, the contribution from aerosol absorption has almost doubled over this period 3,13 . In the future, particularly under weak air-quality policies, these trends towards a more absorbing aerosol distribution will probably continue 2 .
Given these observations, a key question arises: how does the evolving spatial pattern of anthropogenic aerosol (particularly absorbing aerosol, hereafter ' AA') impact climate? Crucially, how does the global-mean effective radiative forcing (ERF)-defined as the net change in radiative fluxes at the top-of-atmosphere in a fixed-SST experiment following a change in forcing agent, allowing for rapid adjustments 14 -respond to this changing pattern of AA? A robust answer to this question would have substantial ramifications because the impact of the spatial pattern of AA on global-mean ERF is not considered either in the simple energy-balance models 15,16 or in integrated assessment models 17 which form an essential part of the IPCC process 18 .
Although the magnitude of, and mechanisms for, a dependence of ERF on aerosol location are not yet clear, there are reasons why we might expect location-dependence. For instance, the dependence of IRF on the effective albedo of the below-aerosol layer means that for AA there is a 'critical surface albedo' above which their IRF switches from negative to positive 19,20 . Additionally, localized AA (and scattering aerosol) can influence atmospheric circulation 4,9,[21][22][23][24][25][26][27] , with subsequent impacts on regional precipitation 5,22,23 and cloudiness 28,29 . The component of ERF arising from changes in clouds and circulation is referred to as the 'rapid adjustments' (Adj) and is known to be a larger fraction of the ERF for AA than scattering aerosol 14 due to the strong influence of AA-induced diabatic heating on moist processes 30 . However, understanding the link between AA, circulation and ERF remains challenging 28 . Aerosol-cloud interactions are also known to depend on both the cloud regime and meteorological conditions [31][32][33] . Thus, through the co-location of aerosols with different cloud types, we may also expect a location-dependence through microphysical changes (however, this pathway is not examined here).
Recent work using emissions perturbations in key geopolitical regions [34][35][36] has also shown that the global-mean ERF is a function of the location of aerosol emissions. However, due to differences in aerosol lifetimes and transport between regions, emissions-driven simulations can exhibit order-of-magnitude differences in regional aerosol burden, despite identical emissions [37][38][39] . Hence, it is not clear whether this 'location-dependence' is because of differences in aerosol transport and lifetimes between regions or whether there is a more fundamental physical mechanism linking aerosol absorption in different regions to different global-mean ERF. Furthermore, perturbing emissions of different aerosol species simultaneously 35,36 makes it difficult to attribute the changes to AA alone. Some of these issues can be circumvented using models driven by aerosol concentrations rather than emissions 40,41 . However, such approaches can still generate regional variations in aerosol radiative properties despite identical concentrations due to changes in the mixing state 38,39,[42][43][44] . Furthermore, as these studies have often been targeted at regions of historical emissions 41 , or made use of climatological distributions 9,29 , it is still difficult to determine the distinct physical mechanisms by which AA in different locations can impact global-mean ERF.
To tackle this, here we present a large ensemble of simulations with a state-of-the-art, atmosphere-only global climate model forced by idealized AOD perturbations (with fixed magnitude and radiative properties) placed systematically across the globe. This approach allows us to explore the sensitivity of ERF to localized aerosol absorption, without being biased to specific regions and without the added uncertainties of aerosol transport and deposition. Our results demonstrate a strong unambiguous dependence of global-mean ERF on the spatial pattern of AA, with positive ERF in response to midlatitude AA and a 'hot spot' of negative ERF for AA in the tropical Western Pacific. We show that these varying ERF responses can be understood by considering the dynamical response of atmospheric circulation and clouds to aerosol heating in different locations and propose a constraint on the ERF of midlatitude AA due to vorticity conservation. Finally, we highlight the importance of our results for large transient ERF from Indonesian peatland wildfires.

Dependence of erF on the location of absorbing aerosol
The global-mean ERF in response to the 144 separate aerosol plume experiments is illustrated in Fig. 1a. The AOD perturbations are prescribed using two-dimensional Gaussian distributions in the horizontal, with an amplitude of 1 at their centre (Methods). The global-mean AOD perturbation is the same in each experiment and only the location is varied (Extended Data Fig. 1). The single-scattering albedo at 550 nm is set to 0.8 in each experiment, corresponding to strongly absorbing aerosol 43 . The variations in ERF in Fig. 1a illustrate that the ERF from a certain global-mean AA perturbation depends strongly on how the AA is spatially distributed. For example, when the AA is concentrated in the extratropics, the ERF is generally positive, whereas when the AA is located over the tropical Western Pacific, there is a strong negative ERF in the global-mean. The two insets below Fig. 1a show that the spatial pattern of the ERF also varies markedly depending on the plume location, with an extratropical plume generating a positive local increase in ERF, whereas the plume over the tropical Western Pacific generates a global-scale ERF perturbation, with large regions of positive and negative response. The spatial pattern of ERF for six more representative plumes is shown in Supplementary Fig. 2.
Next, we decompose the ERF into contributions from IRF and rapid adjustments (Adj) and analyse these in clear and cloudy sky. This allows us to disentangle the direct, radiative effect of the AA (the IRF) from its influence on circulation and cloudiness (the adjustments). For most of the plumes, the contribution from IRF is a balance between the changes in clear and cloudy sky, resulting in an overall small contribution from IRF in most cases (Fig. 1d,e). The cloudy-sky IRF is generally positive because low clouds serve as a high-albedo surface, whereas the clear-sky IRF is negative for most perturbation sites except for over high-albedo surfaces such as the Sahara or the high latitudes.
Variations in global-mean ERF in response to AA are mainly caused by variations in the global-mean rapid adjustments, as shown in Fig. 1b,c. The rapid adjustments in clear sky are well-explained by changes in clear-sky outgoing longwave radiation arising from global-mean temperature changes ( Supplementary Fig. 3), as land temperatures are free to evolve in our simulations and only the SSTs are fixed. The clear-sky adjustments contribute negatively to the total ERF for most locations except the Sahara and the tropical Western Pacific, where they amplify the cloudy-sky contribution. Overall, the cloudy-sky rapid adjustments (Fig. 1b) account for over 63% of the total variation in ERF in response to AA in different locations and are the key determinant of the 'hot spot' of ERF in the tropical Western Pacific. Hence, we mainly focus our attention on the cloudy-sky adjustments.

Circulation changes control the location-dependence of erF
Cloudy-sky adjustments are strongly anti-correlated with global-mean liquid-water path (LWP) changes across all our experiments ( Supplementary Fig. 4) because higher LWP is associated with optically thicker clouds which reflect more incoming solar radiation. But why does AA in some locations cause such an increase in LWP while others cause a decrease? To see why, consider that shortwave absorption from AA introduces a positive temperature tendency to the local atmospheric column, and in the absence of terms that balance this tendency, the temperature in this region would increase without bound. However, in the extratropics (Fig. 2a), temperature anomalies tend to generate a downstream low-pressure system due to vorticity conservation 45 which draws down cold air from higher latitudes ( Supplementary Fig. 5). This cold-air advection in response to AA in the extratropics is a negative feedback that balances the temperature tendency from the AA, generating a finite localized temperature increase. Other local processes such as radiative cooling or precipitation feedbacks also contribute to balancing the temperature tendency from the AA (Supplementary Fig. 7). This temperature increase causes a decrease in relative humidity, which inhibits cloud formation and decreases LWP in the region of the AA perturbation (Fig. 2c). This mechanism explains the broad pattern of positive cloudy-sky adjustments (Fig. 1b) and negative LWP changes ( Supplementary  Fig. 4) throughout the extratropics.
A similar localized decrease in LWP occurs when the AA is in a region of tropical descent ( Fig. 2d-f); however, the mechanism is quite distinct. In the tropics, as the Coriolis force is weak, the diabatic heating from AA is balanced mostly by vertical motion, generating a thermally-direct overturning circulation 46,47 . However, when the AA perturbation is in a region of large-scale descentin this case, the descending branch of the Pacific Walker cell (Fig. 2j)-the anomalous overturning from the diabatic heating interacts destructively with the existing large-scale descent and becomes confined to local scales, causing a decrease in low-level liquid-water content (Fig. 2k).
Conversely, when the AA perturbation is placed in a region of large-scale tropical ascent (Fig. 2g), the anomalous vertical motion generated by the diabatic heating projects onto this existing mode of variability, amplifying the Walker circulation and increasing the upper-level divergence. This pattern of low-level convergence and upper-level divergence excites equatorial waves which generate large positive liquid-water content anomalies throughout the tropics (Fig. 2i,l), and planetary-scale Rossby waves which communicate the perturbation to the midlatitudes 45 (Fig. 2i). A schematic representation of our results is given in Fig. 3.
To distill the key physical mechanisms, we conducted another set of AA perturbation experiments in a simplified aquaplanet model forced by an SST dipole along the equator (a 'mock-Walker' circulation; Methods). This generates an overturning circulation with similar characteristics to the Walker circulation in our control experiment ( Supplementary Fig. 6). Again, we find that when the AA is located over the warm pool region, there is a strengthening of the circulation and a large increase in LWP which drives a shortwave cloudy-sky adjustment, although this is somewhat offset by increases in upper-level ice, which contributes negatively to the cloudy-sky adjustments. Similar increases in upper-level cloud ice are found in our main experiments, albeit to a lesser extent ( Supplementary Fig. 7).
Over land, there is greater convective inhibition due to lower boundary-layer relative humidity 48 , which inhibits the upward motion associated with diabatic heating in the tropics. This helps explain why similarly large ERFs are not seen in response to AA over South America and the African monsoon regions (which are also associated with large-scale ascent). Additionally, AA causes surface cooling in our experiments when located over land 27 . This drives reductions in the surface sensible heat flux and in the diabatic heating perturbation for AA over land ( Supplementary  Fig. 8c), explaining the land-ocean contrast in Fig. 1b. Adjustments of clouds and precipitation can also impact the diabatic heating perturbation through changes in latent heat release or by reductions in outgoing longwave radiation ( Supplementary Fig. 8b,f).
To assess how well this explanation captures the simulated variations in cloudy-sky adjustments (and thus, ERF) across the tropics, we quantify it with a simple scaling, where δQtot is the total diabatic heating from the plume (including latent heating) and (∇ · u) control 850−400 is the difference in the horizontal divergence between 850 hPa and 400 hPa in the control run. Each filled circle denotes the global-mean results from the experiment with the idealized AA plume centred on that location, and the size of the plumes is illustrated in Supplementary Fig. 1. Enclosing black circles are present where the timeseries of global-and annual-mean forcing is statistically significant from zero with P < 0.05 according to a one-sided t-test. Two insets at the bottom show the spatial pattern of ERF changes for two illustrative plume experiments. In each of these subplots, the centre of the AA perturbation in that experiment is marked with an x and its size illustrated by an ellipse which contains 95% of the total AOD perturbation in that experiment. Ellipses are used in the insets to represent the distortion introduced by the map projection. b-e, Decomposition of the total ERF into contributions from Adj (b,c) and IRF (d,e) in clear and cloudy sky.   Both quantities are averaged within a fixed radius of the centre of the AOD perturbation in that experiment. Similar results are obtained if we replace (∇ · u) control 850−400 with the local vertical velocity at 850 hPa ( Supplementary Fig. 9).
Our scaling captures the variation in simulated cloudy-sky adjustments in response to AA perturbations across the tropics (Fig. 4) and demonstrates that this can be understood quantitatively in terms of simple tropical dynamics, as illustrated in Fig. 3. Notably, the scaling also successfully predicts the 'hot spot' of negative ERF over the tropical Western Pacific.

Implications for Indonesian biomass-burning events
The strong 'hot spot' of ERF over the tropical Western Pacific region raises the question: could this have real world implications? One potential candidate is Indonesian wildfire events. Wildfires have been recorded in this region since the late Pleistocene 49 ; however, the worst known fires have occurred since the 1980s due to changes in land use and population density 50 . Large peat deposits make Indonesian wildfires uniquely impactful since peat is highly combustible and difficult to extinguish, and peat fires have repeatedly led to high AOD 'haze' events 51 . Time series of ground-based observations of mid-visible AOD are plotted in Fig. 5a using data from five Aerosol Robotic Network (AERONET) sites located in Indonesia. There is substantial biomass-burning in the dry seasons in most years, particularly in 2015 (associated with a strong El Nino) and 2019 (without an El Nino, but with a positive Indian Ocean dipole).
To investigate whether these events could generate large ERF as suggested by our idealized experiments, we conduct experiments using a realistic AOD climatology of Indonesian biomass-burning (Methods), with a seasonal cycle that peaks in the dry season (September-October-November, SON) ( Supplementary Fig. 10). We perform five experiments in which we successively double the climatological Indonesian plume AOD until the peak dry-season AOD is comparable to the 2015/2016 observations 51 , corresponding to 16 times the climatological AOD (Methods). Despite the use of a more realistic AOD distribution and a seasonal cycle, Fig. 5b shows that the spatial pattern and magnitude of the resultant ERF are remarkably similar to our idealized simulations over the 'hot spot' (for example, Fig. 2h).
The magnitude of the ERF response varies seasonally in response to seasonal shifts in tropical circulation, as it also does in the idealized experiments ( Supplementary Fig. 11). However, the lack of AOD outside of the dry season (SON) in the realistic experiment means that there is little ERF response in March-April-May (MAM), in contrast to the idealized experiment. The magnitude and spatial pattern of the annual-mean ERF response are also robust to realistic variations in the vertical distribution of the total column AOD ( Supplementary Fig. 12).
The global-mean ERF scales linearly with the AOD burden in each of the experiments (Supplementary Fig. 13). This provides evidence that the reason similarly large ERF has not been noted in previous studies which use climatological AOD distributions 29 or emissions from a single year 30 is not because we have misunderstood the physics, but rather because these large biomass-burning events are under-represented in climatological distributions and have strong interannual variability.
Tropical SST variability such as the El Nino Southern Oscillation (ENSO) also plays a role by modulating the pattern of tropical overturning and generating anomalous subsidence over this region, which would inhibit the remote influence on LWP and ERF we see in our experiments (Fig. 3). This helps to explain why previous studies on the ERF of the 1997 Indonesian wildfire 52 (also associated with a strong El Nino) have found a positive ERF, localized to small scales. Our proposed mechanism predicts such a response (Fig. 3).

Discussion
In summary, our knowledge of how the ERF due to AA depends on its spatial pattern is still uncertain because of uncertainty in the translation of AA emissions to regional AA concentrations, and the subsequent interactions of heterogeneous AA heating on atmospheric circulation and clouds. To address this, we have conducted a In b and c, the scaling was evaluated in S.I. units (W m -2 s -1 ) and then, for plotting purposes, each scaling estimate was multiplied by the same dimensional constant such that the results match the magnitude of the cloudy-sky adjustments in a (this does not affect the correlation).
series of experiments with a state-of-the-art climate model in which identical absorbing AOD perturbations are placed across the globe. We have shown that the ERF of AA depends strongly on its spatial pattern, with even the sign of the forcing varying markedly between different perturbation sites. In the extratropics, we find that meridional advection of high-latitude air confines the heating anomaly to local scales, driving a positive ERF through reductions in LWP local to the plume. On the other hand, in the tropics the ERF is determined by interactions between the diabatic heating and the structure of the tropical overturning circulation. In regions of large-scale descent, the diabatic heating anomaly is confined to local scales, generating a similar response to the extratropical case. However, in regions of large-scale ascent, the AA excites substantial wave activity which generates positive LWP anomalies across the globe and a strong negative ERF. Of particular interest is the 'hot spot' of negative ERF we identify in response to AA over the tropical Western Pacific, which suggests that there could be large transient ERF from seasonal biomass-burning over this region. We explored this question using a realistic climatology of Indonesian biomass-burning, constrained by AERONET observations during the 2015/2016 wildfire event, and showed that significant global-mean ERF of ~O(−1 W m −2 ) is possible via our proposed mechanism. Our results demonstrate that Indonesian biomass-burning could represent a large source of negative ERF on interannual timescales, particularly during neutral ENSO or La Niña years 53 .
Although non-local circulation adjustments are included in the definition of ERF 14,54 and are known to play some role in the ERF of well-mixed greenhouse gases 54,55 , their importance for the global-mean aerosol ERF is highly uncertain 30,41 . The inherent uncertainty in linking circulation changes to global-mean ERF is amplified for aerosol due to the additional challenge of understanding how spatially heterogeneous forcers impact circulation 5,[22][23][24]30,36,56 . By taking an idealized approach to this problem, we have shown that non-local circulation adjustments can play a first-order role in determining the ERF of AA, even on global-mean scales. This contrasts with a recent study by Johnson et al. which showed that suppressing circulation changes with nudging did not impact estimates of global-mean ERF 30 . However, this experiment used emissions from a single year (2014) which did not have substantial biomass-burning over the tropical Western Pacific. Overall, this suggests that caution is needed when using nudged experiments to determine the ERF from regional aerosol, particularly for AA. However, quantitative comparisons between our study and that of Johnson et al. 30 are complicated by the fact that the atmospheric response to regional aerosol can be non-linear 57 , and this should be investigated in future work.
One feature of our approach is the use of fixed SSTs (as is required to calculate ERF), which prevents the ocean from responding to the imposed aerosol forcing. Dynamical ocean feedbacks can play a crucial role in determining the overall climate response to an imposed perturbation 5,7,40,56 and sea-surface temperatures can play a role in preconditioning the atmospheric response to aerosol changes 52,58 . For the case of transient aerosol perturbations (such as wildfires on seasonal timescales), one might expect ocean dynamics/thermodynamics to play a smaller role in the real world than is found in previous studies that used long integrations 34,56 . However, this is still uncertain and future work should approach this question using a hierarchy of idealized experiments as illustrated here.
Although the responses we identify are robust across a hierarchy of model configurations and interpretable in terms of robust atmospheric dynamics, we have only used one model and the uncertainties associated with parameterized convection suggest that it would be fruitful to perform similar idealized experiments in other global climate models, or in global storm-resolving models 59 . Our work could also be extended by conducting analogous experiments with scattering aerosol. Both would help to build confidence in the physical mechanisms driving this location-dependence and provide a unique community resource for investigating inter-model differences in aerosol ERF, which are particularly large for AA 29,60 , potentially allowing this location-dependence to be parameterized in simpler models.
Finally, our experiments focus on the ERF from aerosol-radiation interactions and neglect the indirect impacts of aerosol on cloud microphysical processes. However, given the uncertainty surrounding the microphysical interactions between clouds and AA, we consider this to be a reasonable limitation of our study, especially as the uncertainty in ERF from aerosol-radiation interactions contributes considerably to the total uncertainty in aerosol ERF 54 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements and 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-022-01415-4.   Fig. 6). methods Idealized simulations. An ensemble of 145 fixed-SST simulations were conducted using the ICON (ICOsahedral Non-hydrostatic) general circulation model which has a newly developed dynamical core and has been tuned to give a well-balanced top-of-atmosphere radiation budget and reasonable representation of clouds and precipitation compared with observations 61,62 . This configuration of the ICON model uses the ECHAM6 physics package 63 , including a bulk mass-flux convection scheme 64,65 and cloud cover calculated locally in each grid box using a relative humidity scheme 66 . Cloud microphysics is based on the Lohmann and Roeckner scheme 67 ; however, the droplet number is fixed to prescribed profiles over land and ocean 63 and the model does not simulate aerosol-cloud interactions. Radiation is parameterized using the PSrad radiative transfer scheme 68 , which implements the gas radiative transfer and solvers from the widely used Rapid Radiation Transfer Model 69 . We have extended the model so that the total aerosol radiative effect is calculated at every radiation time-step by performing a double-call to the radiation scheme with and without the aerosol perturbation and then subtracting the results. The IRF is then calculated by differencing the total aerosol radiative effect between a perturbed run and the control.
The simulations were run on a triangular grid (R02B04 specification) with a 223 km edge length, corresponding to an average grid spacing of approximately 157.8 km on a Cartesian grid. The model uses a terrain-following vertical sigma-height grid 70 with 47 levels between the surface and the model top at 83 km. Each simulation, including the control run, is integrated for 20 years with timesteps of 15 min and 90 min for the dynamics and radiation, respectively, after which the first year is discarded as spin-up and analysis is conducted over the remaining 19 years. The simulation length was chosen as a trade-off between computational constraints and the requirement that our runs are long enough to give robust results in the face of internal variability.
As a further attempt to reduce variability due to trends in atmospheric composition, each simulation is run with greenhouse gases and ozone fixed at their 1979 levels. Additionally, the simulations are also forced by climatological monthly SSTs and sea-ice concentrations derived from the Atmospheric Model Intercomparison Project (AMIP) 71 boundary conditions over 1979-2016, eliminating SST-driven variability such as the ENSO, which could otherwise bias the results. Supplementary Fig. 1c,d illustrates the annual-mean SST field used to force the model, demonstrating neutral ENSO conditions with a warm pool in the Western Pacific and a cold pool in the Eastern Pacific.
Aside from the control run, the only difference between the remaining 144 simulations is the inclusion of an idealized plume of absorbing aerosol at 144 different geographical locations ( Supplementary Fig. 1a), parameterized using the Max Planck Institute Aerosol Climatology version 2, Simple Plume (MACv2-SP) model 72 . The plume's horizontal structure is specified to be a 2-dimensional Gaussian function with a standard deviation of 10°. The vertical distribution of aerosol is based on the kernel of Euler's β-function (see details in Stevens et al. 72 ), with most of the aerosol concentrated in the lower 3 km of the atmosphere and a peak at ~1 km ( Supplementary Fig. 1b). We note that differences in the vertical plume profile may cause different responses of the climate system and general circulation 73,74 ; however, a vertical distribution similar to ours has been found to be broadly representative of the vertical structure of fine-mode AOD 75 . Furthermore, sensitivity tests for the Indonesian wildfire perturbation show that our results are robust to realistic changes in the vertical aerosol profile ( Supplementary Fig. 11). Both the horizontal and vertical structure of the plume are constant in time and the shape of the plume is not influenced by the local meteorology.
In each of the aerosol perturbation experiments, the plume is prescribed with a total column AOD = 1 and a single-scattering albedo (SSA) of 0.8, both defined at a wavelength of 550 nm. This corresponds to a strongly absorbing plume of aerosol, generating a large diabatic heating perturbation in the lower atmosphere, even while generating cooling at the surface. Furthermore, to isolate the role of the radiative impact of the plume, we do not include the microphysical effects of the aerosol plume, which are also described in Stevens et al. 72 , so the radiative forcing we calculate is the effective radiative forcing from aerosol-radiation interactions (ERF ARI ) in the framework of the IPCC Sixth Assessment Report 1, 54 .
The location of the 144 idealized aerosol plumes was chosen to provide broad coverage over the entire globe while also ensuring that there are several plumes in semi-realistic locations. For example, our setup includes plumes over Europe, North America, South America, Australia, China, India, Central Africa and the Maritime Continent alongside far less realistic plumes of absorbing aerosol over the remote oceans and polar regions. We chose to take this approach for two reasons: first, although aerosol perturbations are highly localized over their source region, it is possible for them to be advected long distances by large-scale circulation. Second, although many of these plumes are less realistic, they allow us to sample many different surface properties and local and large-scale meteorological conditions, which aids us in building a physical picture of what sets the response.

Calculation of ERF and fast adjustments.
We calculate the ERF as the difference between the net top-of-atmosphere radiative fluxes between a perturbation experiment and the control (after discarding the first year as spin-up and taking averages over the subsequent 19 years), following the suggestion of Forster et al. 76,77 . All experiments are conducted with fixed SSTs and sea-ice concentrations.
We then use the IRF, defined as the difference between aerosol radiative effect in the perturbed and control experiments, to calculate the contributions from fast adjustments in clear and cloudy sky as: where ERF clear and ERF cloudy are defined the same way as the ERF but only using top-of-atmosphere fluxes in clear and cloudy sky, respectively.
Clear-sky rapid adjustments. To understand the variations in clear-sky rapid adjustments, we note that, globally, the outgoing longwave radiation (OLR) can be written as: where σT 4 is the surface blackbody emission, and ϵ is an 'effective emissivity' defined such that equality holds at a given surface temperature. In our control experiment, with OLR ≈ 240 W m −2 and T ≈ 280 K, we find ϵ ≈ 0.8. Given a change in global-mean surface temperature, one can then relate the clear-sky adjustments to the first-order Taylor expansion about the control state as: This scaling is found to reproduce most of the variance in clear-sky rapid adjustments in our experiments ( Supplementary Fig. 3).
Aquaplanet simulations. We run our ICON aquaplanet simulations for 10 years, with the first year discarded as spin-up and analysis conducted on time-averages of the final 9 years. The SST is prescribed in each experiment, with the zonal-mean SST taken as the 'control' distribution from Neale and Hoskins 78 , in addition to a wavenumber-2 SST dipole which spans one half-hemisphere and takes the form Here, λ is the longitude, ϕ is the latitude and ϕ0 = π 6 = 30 o is the latitudinal width of the perturbation. The amplitude of the dipole is thus ±3 K. The AA perturbations are prescribed at regular longitude intervals spanning between the peaks of the SST dipole ( Supplementary Fig. 6) and the aerosol properties are identical to those used in the main experiments.

Stream-function calculation.
To calculate the zonal-mass stream function for tropical overturning, we first use the 'windspharm' library 78 to calculate the divergent component of the zonal wind u div using a Helmholtz decomposition 79 , and then calculate the zonal-mass stream function at each pressure p and latitude ϕ using: u div (p′, ϕ)dp′.
Here, R e is the earth's radius. To plot longitude-pressure cross-sections, we then average ψ between ±20 o of the equator using area-weighting.

Indonesian plume experiments.
We conduct experiments with the realistic climatology of Indonesian biomass-burning from Stevens et al. 72 , which has a simple seasonal cycle peaking in the dry season ( Supplementary Fig. 10) and a prescribed single-scattering albedo of 0.87. We then run a series of five 10 year experiments, with the AOD in each experiment scaled by an extra factor of 2 compared with the climatology. The vertical profile of the aerosol is the same in each experiment and only the magnitude is varied. The 'Climatology_x16' experiment has a seasonal cycle of AOD which compares favourably with the 2015/2016 AERONET observations and the results of Kiely et al. 51 . For more details, see Stevens et al. 72 .

Data availability
Model data from this study are available at https://zenodo.org/record/6501508. Questions about the data can be directed to the corresponding author. Each plume corresponds to a separate experiment and a filled contour illustrating the horizontal structure of the plume AOD is shown for the plume centered over the Indian Ocean. b Vertical profile of fractional aerosol optical depth at 533 nm with respect to the column integral. The vertical axis denotes the approximate altitude above the local surface. The vertical profile is constant for all the perturbation experiments we conduct. c, d Annual-mean SST field used to drive the simulations, along with the corresponding SST anomaly relative to the zonal-mean temperature.