The 2019/20 Australian wildfires generated a persistent smoke-charged vortex rising up to 35 km altitude

The Australian bushfires around the turn of the year 2020 generated an unprecedented perturbation of stratospheric composition, dynamical circulation and radiative balance. Here we show from satellite observations that the resulting planetary-scale blocking of solar radiation by the smoke is larger than any previously documented wildfires and of the same order as the radiative forcing produced by moderate volcanic eruptions. A striking effect of the solar heating of an intense smoke patch was the generation of a self-maintained anticyclonic vortex measuring 1000 km in diameter and featuring its own ozone hole. The highly stable vortex persisted in the stratosphere for over 13 weeks, travelled 66,000 km and lifted a confined bubble of smoke and moisture to 35 km altitude. Its evolution was tracked by several satellite-based sensors and was successfully resolved by the European Centre for Medium-Range Weather Forecasts operational system, primarily based on satellite data. Because wildfires are expected to increase in frequency and strength in a changing climate, we suggest that extraordinary events of this type may contribute significantly to the global stratospheric composition in the coming decades. The 2019/2020 Australian wildfires generated a smoke cloud that organized itself into a persistent vortex structure and ascended to 35 km altitude through solar heating, according to satellite tracking.

T he impact of wildfire-driven thunderstorms on the global stratosphere has been deemed small until the North American wildfires in August 2017. Pyro-cumulonimbus (pyroCb) clouds from that event caused stratospheric perturbations an order of magnitude larger than the previous benchmarks of extreme pyroCb activity and approached the effect of moderate volcanic eruption 1,2 . Volcanic eruptions inject ash and sulphur which is oxidized and condenses to form submicron-sized aerosol droplets in the stratosphere. With the PyroCb, intense fire-driven convection lifts combustion products in gaseous form as well as particulate matter including organic and black carbon, smoke aerosols and condensed water. The solar heating of the highly absorptive black carbon propels the smoke-laden air parcels upward 1 , which, combined with horizontal transport 3,4 , leads to a more efficient meridional dispersion of these aerosols and prolongs their stratospheric residence time 5 .
The Australian bushfires that raged in December 2019-January 2020 have put a new benchmark on the magnitude of stratospheric perturbations. In this study, we use various satellite observations to quantify the magnitude of hemispheric-scale perturbation of stratospheric gaseous compounds and aerosol loading caused by these wildfires. The radiative forcing of the stratospheric smoke is estimated using a radiative transfer model supplied by satellite observations of aerosol optical properties. Finally, using the operational forecasting system of the European Centre for Medium-Range Weather Forecasts (ECMWF) 6 we show that the solar heating of an intense smoke patch has led to generation of a quasi-ellipsoidal anticyclonic vortex which lofted a confined bubble of carbonaceous aerosols and water vapour up to 35 km altitude in about 3 months.

Results
Large-scale perturbation of the stratosphere. The Australian wildfire season 2019/2020 was marked by an unprecedented burn area of 5.8 million hectares (21% of Australia's temperate forests) 7 and exceptionally strong PyroCb activity in the southeast of the continent 8 . The strongest PyroCb outbreak occurred on the New Year's Eve (Fig. 1a) and on the 1st of January the instantaneous horizontal extent of the stratospheric cloud amounted to 2.5 million km 2 as inferred from nadir-viewing TROPOMI 9 satellite measurement (Fig. 1b). On that day, an opaque cloud of smoke was detected in the stratosphere by the CALIOP space-based laser radar (lidar) 10 at altitudes reaching 17.6 km (Fig. 2). Another PyroCb outbreak with stratospheric impact, although less vigorous, took place on 4 January 2020 and on 7 January, the horizontal extent of the stratospheric smoke cloud peaked at 6.1 million km 2 (Fig. 1b) extending over much of the Southern midlatitudes ( Fig. 2 and Supplementary Fig. 4c).
The high-altitude injections of smoke rapidly tripled the stratospheric aerosol optical depth (SAOD) in the southern extratropics. The SAOD perturbation has by far exceeded the effect on stratospheric aerosol load produced by the North American wildfires in 2017, putting the Australian event on par with the strongest volcanic eruptions in the last 25 years (Fig. 3), i.e. since the leveling off of stratospheric aerosol load after a major eruption of Mount Pinatubo in 1991 11 . Three months after the PyroCb event, the SAOD perturbation has remained at the volcanic levels, gradually decreasing with a rate similar to the decay of stratospheric aerosol produced by moderate volcanic eruptions.
Using aerosol extinction profiles retrieved from the limbviewing NASA OMPS-LP instrument 12 we find the total aerosol particle mass lofted into the so-called stratospheric "overworld" 13 (above 380 K isentropic level corresponding to~12-17 km altitude) is 0.4 ± 0.2 Tg (Fig. 4), which is nearly three times larger than the estimates for the previous record-high North American wildfires 2 . The increase in the stratospheric abundance of the gaseous combustion products, derived from the NASA Microwave Limb Sounder (MLS) satellite observations 14 , is as remarkable as the aerosol increase. Figure 4 puts in evidence that the stratospheric masses of carbon monoxide (CO) and acetonitrile (CH 3 CN) bounded within the southern extratropics Latitude-altitude evolution of the smoke plumes in the stratosphere. The pixels, colour coded by date, indicate doubling of aerosol extinction with respect to December 2019 levels for data where aerosol to molecular extinction ratio is 1 or higher. The black circles with date-colour filling indicate the locations of high amounts of water vapour and/or carbon monoxide detected by MLS (see "Methods"). The black contour encircles the locations of aerosol bubble detections by CALIOP lidar (see Fig. 5a and Supplementary Fig. 3). The cross marks the latitude-altitude extent of the stratospheric cloud detected by CALIOP on the 1st January (see Fig. 5a). The grey solid and black dashed curves indicate respectively the zonal-mean 380 K isentrope and the lapse-rate tropopause for the January-March 2020 period. Fig. 4 Time evolution of the daily total mass of CO, CH 3 CN, H 2 O and aerosols above the 380 K potential temperature, between 20°S and 82°S. The dotted and solid lines correspond to daily data and 1-week smoothed data, respectively. Envelopes represent two standard deviations over the 1week window (see "Methods"). As shown in this figure, the levels of CO, CH 3 CN, H 2 O, and aerosols started to increase simultaneously and kept increasing during~2-3 weeks, a duration corresponding probably to the time taken by products injected in the lowermost stratosphere to ascend above 380 K. The stratospheric masses of carbon monoxide (CO) and acetonitrile (CH 3 CN) bounded within the southern extratropics increase abruptly by 1.5 ± 0.9 Tg and 3.7 ± 2.0 Gg, respectively during the first week of 2020. This gives a CO/CH 3 CN mass ratio of 0.0025, consistently with previous estimates for temperate Australian wildfires 42 . The injected mass of water was estimated at 27 ± 10 Tg that is about 3% of the total mass of stratospheric overworld water vapour in the southern extratropics. The shading shows that the amplitude of fluctuations increases sharply during the sharp rise of species masses, reflecting the fact that sampling of the bubble by MLS is more random than on a more homogeneous field. The lagging increase of the aerosol mass is due to the fact that the OMPS-LP extinction retrieval saturates at extinction values above 0.01 km −1 . Profiles are, therefore, truncated below any altitude exceeding this value, which can lead to an underestimation of the early aerosol plume when it is at its thickest. This artifact, which explains the slower increase of aerosol mass than gases, persists until mid-February when the plume is sufficiently dispersed so that OMPS-LP extinction measurements no longer saturate. increase abruptly by 1.5 ± 0.9 Tg (~20% of the pre-event levels) and 3.7 ± 2.0 Gg (~5%), respectively, during the first week of 2020. The injected mass of water was estimated at 27 ± 10 Tg that is about 3% of the total mass of stratospheric overworld water vapour in the southern extratropics (see "Methods").
The gases and particles injected by the PyroCbs were advected by the prevailing westerly winds in the lower stratosphere. The patches of smoke dispersed across all of the Southern hemisphere extratropics in less than two weeks with the fastest patches returning back over Australia by 13 January 2020, whereas the carbon-rich core remained bounded within midlatitudes as shown in Fig. 2. During the following months, most of the particulate material dwelled in the lower stratosphere, the larger and heavier particles sedimented to lower altitudes while the carbon-rich fraction ascended from 15 to 35 km due to solar heating of black carbon ( Supplementary Fig. 1).
Radiative forcing. The large amount of aerosols produced a significant radiative forcing (RF), which we quantified using explicit radiative transfer modelling based on the measured aerosol optical properties (Supplementary notes 1). In the latitude band between 25 and 60°S, an average cloud-free reference monthly radiative forcing as large as about −1.0 W/m 2 at the top of the atmosphere (TOA) and −3.0 W/m 2 at the surface is found in February 2020 ( Supplementary Fig. 2). This can be attributed to perturbation to the stratospheric aerosol layer by the Australian fires plumes. The area-weighted global-equivalent cloudfree RF is estimated (Supplementary Table 1) to values as large as −0.31 ± 0.09 W/m 2 (TOA) and −0.98 ± 0.17 W/m 2 (at the surface). It is important to notice that these estimations don't take the presence of clouds into account and are to be taken as purely reference values. For typical average cloud cover in the area affected by the plume 15 , the surface all-sky RF can be reduced tõ 50% and the TOA all-sky RF to~30-50% of the clear-sky RF estimations 16 (see Supplementary notes 1 for details). From the perspective of the stratospheric aerosol layer perturbation, the global TOA RF produced by the Australian fires 2019/2020 is larger than the RF produced by all documented wildfire events and of the same order of magnitude of moderate volcanic eruptions during the last three decades (that have an integrated effect estimated at 17 −0.19 ± 0.09 W/m 2 , or smaller 18 ). In contrast to the non-absorbing volcanic sulphates, the carbonaceous wildfire aerosols absorb the incoming solar radiation, leading to yet more substantial radiative forcing at the surface, due to the additional large amount of energy absorbed in the plume. This can be linked to the ascent of a smoke cloud in the stratosphere, which is dicussed in the next section.
Rising bubble of smoke. The primary patch of smoke originating from the New Year's Eve PyroCb event followed an extraordinary dynamical evolution. By the 4 January 2020, en route across the Southern Pacific, the core plume started to encapsulate into a compact bubble-like structure, which was identified using CALIOP observations on 7 January 2020 as an isolated 4-km tall and 1000 km wide structure (Fig. 5a). Over the next 3 months, this smoke bubble crossed the Pacific and hovered above the tip of South America for a week. It then followed a 10-week westbound round-the-world journey that could be tracked until the beginning of April 2020 (Supplementary notes 2 and Supplementary Fig. 3), travelling over 66,000 km.
The large amount of sunlight-absorbing black carbon contained in the smoke cloud provided a localized heating that forced Vertical evolution of the smoke bubble, its chemical composition and thermal structure. a Selection of CALIOP attenuated scattering ratio profiles for clear intersections of the bubble by the orbit except the first panel of 1st of January that shows the dense and compact plume on its first day in the stratosphere. The attenuated scattering ratio is calculated by dividing the attenuated backscattering coefficient by the calculated molecular backscattering. The data are further filtered horizontally by an 81-pixel moving median filter to remove the noise. The crosses show the projected interpolated location of the vortex vorticity centroid from the ECMWF operational analysis onto the orbit plane at the same time and the white contour shows the projected contour of the half maximum vorticity value in the pane passing by the vortex centroid and parallel to the orbit plane (see Fig. 6). b Evolution of the water vapour mixing ratio within the rising bubble based on MLS bubble detections (see "Methods"). The dashed contours show the equivalent mixing ratio of ice water derived from MLS ice water content vertical profiles collocated with the bubble. The thick dashed curve marks the top altitude of the aerosol bubble determined as the level where OMPS-LP extinction triples that of the nearest upper altitude level. c Evolution of carbon monoxide (MLS) within the bubble. The centroid and the vertical boundaries of the aerosol bubble determined using CALIOP data are overplotted as circles and bars respectively. d Composited temperature perturbation within the smoke bubble from Metop GNSS radio occultation (RO) temperature profiles collocated with the smoke bubble (see "Methods"). The black line shows the centroid of vortex detected from ECMWF data (see Supplementary notes 4.1 and Supplementary Fig. 7b).
the air mass to rise through the stratosphere. With an initial ascent rate of about 0.45 km/day, the bubble of aerosol continuously ascended during the three months with an average rate of 0.2 km/day. While remaining compact, the bubble was leaking material from its bottom part, leaving an aerosol trail that was progressively dispersed and diluted, filling the whole midlatitude austral stratosphere up to 30 km ( Supplementary Fig. 1, see also Fig. 10b). The rise ceased in late March 2020 when the top of the bubble reached 36 km altitude ( Fig. 5b and Supplementary Fig. 3). This is substantially higher than any coherent volcanic aerosol or smoke plume observed since the major eruption of Pinatubo in 1991.
Along with the carbonaceous aerosols, the bubble entrained tropospheric moisture in the form of ice aggregates injected by the overshooting PyroCbs. In the warmer stratosphere, the ice (detected by MLS sensor as high as 22 km, cf. Figure 5b) eventually evaporated, enriching the air mass with water vapour. This led to extraordinary high water vapour mixing ratios emerging across the stratosphere within the rising smoke bubble (Fig. 5b). The decay of CO within the bubble (Fig. 5c) was faster than that of water vapour, reflecting the fact that, unlike water vapour, the carbon monoxide was subject to photochemical oxidation whose efficiency increases sharply with altitude 19 .
Temperature profiles from GNSS radio occultation sensors exhibit a clear dipolar anomaly within the bubble with a warm pole at its bottom and a cold pole at its top (Fig. 5d). Although counterintuitive from the pure radiative transfer perspective, the observed temperature dipole within the heated cloud represents an expected thermal signature of a synoptic-scale vortex.
The vortex. The compact shape of the smoke bubble could only be maintained through an efficient confinement process. The meteorological analysis of the real-time operational ECMWF integrated forecasting system (IFS) 6 reveals that a localized anticyclonic vortex was associated with the smoke bubble during all its travel, moving and rising with it (Figs. 5, 6a, b and Supplementary Fig. 5). With a peak vorticity of 10 −4 s −1 (Figs. 6c, 7b) and a maximum anomalous wind speed of 13 m s −1 (Fig. 6d) during most of its lifetime, the vortex had a turnover time of about 36 h. It has therefore survived about 60 turnover times demonstrating a remarkable stability and resilience against perturbations. The ascent was surprisingly linear in potential temperature at a rate of 5.94 ± 0.07 K day −1 (Fig. 7a). This corresponds to a heating rate dT/dt which varies from about 3 K day −1 at the beginning of January to 1.5 K day −1 at the end of March. The altitude rise is from 16 to 33 km for the vortex centroid and from 17 to 36 km for the top of the bubble according to CALIOP. The upper envelope of the OMPS detection of the bubble is seen as the cyan curve on Fig. 7a. This envelope is always above the top detected from CALIOP. Such a bias is expected as OMPS-LP is a limb instrument that scans a much wider area than the narrow CALIOP track. Both CALIOP and OMPS-LP detect that the top of the bubble rises initially faster than the vortex core, by about 10 K day −1 . This period corresponds to the initial travel of the bubble to the tip of South America. In terms of altitude ascent, the rates 10 and 5.94 K day −1 translate approximately as 0.45 and 0.2 km day −1 .
The confining properties of the vortex are confirmed by the colocated anomalies in tracers and aerosol from the TROPOMI instrument. The satellite observations (Supplementary notes 3.3 and Supplementary Fig. 4a, b) reveals an isolated enhancement of the aerosol absorbing index and of the CO columnar content, as well as the presence of a deep mini ozone hole, depleted by up to 100 DU. All three features were captured in the same position by the ECMWF analyses ( Fig. 6e and Supplementary Figs. 4 and 7).
Companion vortices. It is worth noticing that the vortex was not a single event. It had several companions, albeit of smaller magnitude and duration, also caused by localized smoke clouds. The most noticeable lasted one month and travelled the hemisphere. Another one found a path across Antarctica where it was subject to the strong aerosol heating of permanent daylight and rose up to 27 km.
The second vortex is borne from the smoke cloud that found its way to the stratosphere during the PyroCb event of 4-5 January 2020. This cloud initially travelled north east passing north of New Zealand before taking a south easterly direction crossing the path of the cloud emitted on 31 December 2019 and the main vortex. Fig. 8a, b shows that a vortex-like structure can be spotted as early as 7 January, coinciding with the location of a compact bubble according to CALIOP. Subsequently, the bubble crossed the path of the first vortex on 16 January while rising and intensifying (see Fig. 8c, d)  The third bubble has been first detected by CALIOP on 7 January at 69°S/160°W. It then moved over Antarctica (Fig. 9a) until the end of January where it spent a week over the Antarctic Peninsula before moving to the tip of South America, shortly after this region was visited by the main vortex. It eventually moved to the Atlantic where it dissipated by 25 February (Fig. 9b). The bubble was accompanied by a vortex during its whole life cycle as seen in Fig. 9c. Although the magnitude of this vortex was modest compared to the main one and even the second one in terms of maximum vorticity (Fig. 9e), it performed a very significant ascent from 18 to 26 km (Fig. 9d). We attribute this effect to the very effective aerosol heating received during the essentially permanent daylight of the first period over Antarctica. The simultaneous rise of the main vortex and the third vortex is very clear from the OMPS-LP latitude-altitude sections in Fig. 10c, d.
The combined trajectories of the main vortex and its two companions are shown in Fig. 10e. Figure 10a shows how the trajectories of the vortices form the skeleton of the dispersion path of the smoke plume in the stratosphere. Figure 10b shows that OMPS-LP follows closely the evolution of the top altitude of all the vortices under the shape of welldefined branches in a longitude-time Hovmöller diagram. It is worth noticing that both the main and the third vortex spent some time wandering in the vicinity of the Drake passage at the beginning of February. Such a stagnation situation is prone to sensitivity. The IFS forecasts during the end of January predicted that the main vortex would cross to the Atlantic, while instead it did not and began moving westward over the Pacific as it reached a higher altitude where easterly winds prevail. Several secondary branches that seem to separate from the main one followed by the main vortex are also visible in this diagram. A detailed inspection reveals that they are indeed associated with patches left behind by the main bubble as it moved upward. It is apparent from several of the panels of Fig. 7a that the top part of the bubble remained always compact while the bottom part was constantly leaking material. Figure 10c, d shows latitude-altitude cross sections of OMPS aerosols on 16 January and 1 February. The fast rise of the main vortex during that period is visible near 50S while the second vortex is located at lower altitude and the third vortex corresponds to the towering structure by 75-80S.  . In a, b, the black lines show the ECMWF 10-day forecast evolution, plotted every four days. The forecast evolution is only shown for the period where it maintains the vortex. The slight discrepancy between the analysis and the initial point of the forecast is because the 10-day forecast is produced from a slightly inferior 6-h analysis, due to real time constraints.

Discussion
Long-lived anticyclones have also been observed as very rare events in the summer Arctic stratosphere [20][21][22] , but they are much larger structures of about 2000 km, very near the pole, and do not display any of the specific characters of the self-generated smokecharged vortex. According to geophysical fluid dynamics theory 23 , a local heating in the austral stratosphere is expected to produce positive potential vorticity aloft and destroy it beneath. The positive potential vorticity is partially realized as anticyclonic rotation and partially as temperature stratification. The negative potential vorticity is apparently dispersed away with the tail of the bubble. The ECMWF analysed thermal structure ( Fig. 6f and Supplementary Fig. 7c) shows the same dipole as in the GNSS-RO satellite profiles with the same amplitude. The observed vortex is quite similar to known isolated ellipsoidal solutions of the quasigeostrophic equations 24 which explain some of the long-lived vortices in the ocean 25 , but such a structure is described here for the first time in the atmosphere.
The ECMWF analysis uses climatological aerosol fields, so it does not take the aerosol emissions by the Australian wildfires into account, nor the satellite measurements of aerosol extinction. Thus, the replication of the vortex by the analysis was due to assimilation of temperature, wind and ozone measurements from operational satellites. We found that the radio-occultation temperature profiling satellite constellation played the key role in the successful reconstruction of this unusual stratospheric phenomenon by ECMWF analyses (Supplementary notes 4.2). The analyses were also influenced by a few radiosonde profiles and wind profiles in the lower stratosphere by the ESA Aeolus satellite's Doppler wind lidar 26 , which provided an observational evidence of the anticyclonic vortex.
The ECMWF IFS produce analyses using a 4-dimensional variational data assimilation method 27 , combined with a highresolution time-evolving forecast model that predicts the atmospheric dynamics. The assimilation system updates the state of the atmosphere using satellite data and in situ observations. As aerosols are not assimilated in the IFS, the ECMWF forecast was not expected to maintain the vortex as observed. Indeed, as shown in Fig. 7a, b, the forecast consistently predicts a decay of the vortex amplitude and fails to predict its rise except on 3rd March where the vorticity centroid underwent a jump following the stretching and breaking of the vortex under the effect of vertical shear a few days earlier. As the nature of this event was purely dynamical, it was correctly predicted. This trend would lead to a rapid loss of the vortex and is corrected by the assimilation of new observations, providing an additional forcing that ensures maintenance and rise of the structure.
The ozone hole is forced by the assimilation of satellite observations of ozone ( Supplementary Fig. 6b, f). The physical mechanism leading to the ozone hole must be a combination of the uplift of ozone-poor tropospheric air and ozone-depleting chemistry in the smoke cloud.
The generation of the smoke-charged vortex is reported in another study 28 , of which we were made aware during the review process. They used MLS, OMPS and CALIOP satellite instruments as well as the US Navy Global Environmental Model (NAVGEM) analysis to characterize the evolution of chemical composition and structure of the main vortex until 10th March 2020, that is three weeks before it has reached its apogee and collapsed. Our analysis not only covers the entire lifetime of all the smoke-generated vortices, but provides a more comprehensive description of their dynamical evolution, quantifies the largescale perturbation of the stratospheric composition and radiative balance, and describes the impact of satellite data assimilation. While the ref. 28 reveals a number of remarkable similarities with respect to our approach and analysis, they report an average diabatic ascent rate of 8 K day −1 whereas we find 5.9 K day −1 . An important point regarding the vortex dynamics and maintenance, which is not discussed in ref. 28 , is the need of a mechanism to suppress the negative potential vorticity produced by the isolated heating as discussed hereinbefore.

Conclusions
The observed and modelled planetary-scale repercussions of the Australian PyroCb outbreak around the turn of 2020 revolutionize the current understanding and recognition of the climatealtering potential of the wildfires. A single stratospheric overshoot of combustion products produced by the New Year's eve PyroCb event has led to an unprecedented hemispheric-scale perturbation of the stratospheric gaseous and aerosol composition, radiative balance and dynamical circulation with a prolonged effect. Whilst rivalling the volcanic eruptions in terms of stratospheric aerosol load perturbation, this exceptionally strong wildfire event had a substantial impact on a number of other climate-driving stratospheric variables such as water vapour, carbon monoxide and Fig. 10 Three-dimensional evolution of the three vortices from OMPS-LP and ECMWF IFS. a Time-latitude section of zonal-mean stratospheric aerosol optical depth (above the tropopause) from OMPS-LP measurements. The markers show locations of smoke-charged vortices identified using ECMWF vorticity fields. b Longitude-temporal evolution (Hovmöller diagram) of the maximum altitude of smoke plume inferred from OMPS-LP extinction data within 30°S-60°above 15 km where aerosol to molecular extinction ratio exceeds 5. The markers indicate the locations of the smoke-charged vortices identified using ECMWF vorticity fields. c Latitude-altitude section of zonal-mean aerosol to molecular extinction ratio above the local tropopause from OMPS-LP measurements for 16 January 2020. The thick and the thin curves indicate, respectfully, the zonal-mean lapse-rate tropopause and the 380 K potential temperature level. The markers show the positions of the main vortex and of the two companion vortices. d Same as (c) for 1 February 2020. e Trajectories of the three vortices with colour-coded date in the longitude-latitude plane.
ozone. As the frequency and intensity of the Australian wildfires is expected to increase in the changing climate 29 , it is possible that this type of extraordinary event will occur again in the future eventually becoming a significant contributor to the global stratospheric composition.
This work reports the self-organization of an absorbing smoke cloud as a persistent coherent bubble coupled with a vortex that produces confinement preserving the compactness of the cloud. This structure is maintained and rises in the calm summer stratosphere due to its internal heating by solar absorption. The intensity, the duration and the extended vertical and horizontal path of this event certainly ranks it as extraordinary. More detailed studies will be necessary to understand fully the accompanying dynamical processes, in particular how a single sign vorticity structure emerges as a response to heating. Whether stratospheric smoke vortices have already occurred during previous large forest fires is to be explored.  30 . Aerosol extinction coefficient and ozone number density profiles are retrieved from the limb radiance using a two-dimensional tomographic inversion 12 and a forward model that accounts for multiple scattering developed at the University of Saskatchewan 31 . This OMPS-LP USask aerosol product is retrieved at 746 nm and has a vertical resolution of 1-2 km throughout the stratosphere. The aerosol extinction profiles are exploited to analyze the spatiotemporal evolution of the smoke plumes and for computing the mass of particulate matter lofted above the 380 K potential temperature level corresponding to~12-17 km altitude in the extratropics, see Fig. 2). The aerosol mass is derived from the aerosol extinction data assuming a particle mass extinction coefficient of 4.5 m 2 g −1 (ref. 2 ). SAGE III. The Stratospheric Aerosol and Gas Experiment (SAGE) III provides stratospheric aerosol extinction coefficient profiles using solar occultation observations from the International Space Station (ISS) 32 . These measurements, available since February 2017, are provided for nine wavelength bands from 385 to 1550 nm and have a vertical resolution of approximately 0.7 km. The SAGE III/ISS instrument and the data products have characteristics nearly identical to those from the SAGE III Meteor mission 33 . Here we use the 754 nm wavelength band for quantifying the error of OMPS-LP aerosol extinction retrieval (using 16-84 percentiles) and the wavelength pair 1019/521 nm for deriving the Angstrom exponent, which is used for the radiative forcing calculations.
CALIOP. The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is a two-wavelength polarization lidar on board the CALIPSO mission 10 that performs global profiling of aerosols and clouds in the troposphere and lower stratosphere. We use the total attenuated 532 nm backscatter level 1 product V3. 40 14 instrument on the NASA Aura satellite has been measuring the thermal microwave emission from Earth's atmospheric limb since July 2004. With~15 orbits per day, MLS provides day and night near-global (82°S-82°N) measurement of vertical profiles of various atmospheric gaseous compounds (including H 2 O, CO and CH 3 CN) cloud ice, geopotential height, and temperature of the atmosphere. The measurements yield around to 3500 profiles per day for each species with a vertical resolution of~3-5 km. For tracking the smoke bubble, we selected profiles bearing CO enhancements in the stratosphere exceeding 400 ppbv and/or H 2 O enhancement exceeding 12 ppmv with respect to the pre-event conditions (late December 2020). Stratospheric mass loads of H 2 O, CO and CH 3 CN are derived from MLS volume mixing ratio measurements of species in log pressure space, molecular mass of the compound and the air number density derived from MLS temperature profile on pressure levels above the 380 K isentropic level and between 20°S and 82°S. The error bars on the mass of injection are estimated by combining accuracies on the measurements and the mean standard deviations over 20-day periods before and after the sharp increase. The error bar on the CH 3 CN mass above 380 K is only calculated on the standard deviation because accuracies on CH 3 CN measurements are extremely large. The error bar on the aerosol mass takes also into account the uncertainty on the particle mass extinction coefficient (1.5 m 2 g −1 ).
TROPOMI. The TROPOspheric Monitoring Instrument (TROPOMI) is a nadirviewing shortwave spectrometer jointly developed by the Netherlands Space Office and the European Space Agency on board the Sentinel 5 Precursor mission. The instrument has spectral bands in the ultraviolet (270-500 nm), the near-infrared (710-770 nm) and the shortwave infrared (2314-2382 nm) providing therefore observation of key atmospheric constituents, among which we use O 3 , CO and aerosol index at high spatial resolution (7 × 3.5 km 2 at nadir for the UV, visible and near-infrared bands, 7 × 7 km 2 at nadir for shortwave infrared bands) 34 . Here we exploit the Aerosol Index, and the CO and O 3 columnar values from the offline Level 2 data. The aerosol absorbing index (AI) is a quantity based on the spectral contrast between a given pair of UV wavelength (in our case the 340-380 nm wavelengths couple). The retrieval is based on the ratio of the measured top of the atmosphere reflectance (for the shortest wavelength) and a pre-calculated theoretical reflectance for a Rayleigh scattering-only atmosphere (assumed equal for both wavelengths). When the residual value between the observed and modelled values is positive, it indicates the presence of UV-absorbing aerosols, like dust and smoke. Negative residual values may indicate the presence of non-absorbing aerosols, while values close to zero are found in the presence of clouds. AI is dependent upon aerosol layer characteristics such as the aerosol optical thickness, the aerosol single scattering albedo, the aerosol layer height and the underlying surface albedo. Providing a daily global coverage and a very high spatial resolution (7 × 3.5 km 2 at the nadir), the AI from TROPOMI is ideal to follow the evolution of smoke, dust, volcanic ash, or aerosol plumes. GNSS-RO. We use Global Navigation Satellite System (GNSS) radio occultation (RO) dry temperature profiles acquired onboard Metop A/B/C satellites and processed in near real time mode at EUMETSAT RO Meteorology Satellite Application Facility (ROM SAF) 35 . For computing the composited temperature perturbation within the smoke bubble we use temperature profiles collocated with the vortex centroid as identified using IFS analyses (8 h, 400 km collocation criteria). The perturbation is computed as the departure from a mean temperature profile within the corresponding spatiotemporal bin (3-day, 3°latitude, 40°l ongitude).
ECMWF IFS. ECMWF IFS is the operational configuration of the ECMWF global Numerical Weather Prediction system (46R1, https://www.ecmwf.int/en/ publications/ifs-documentation). It consists of an atmosphere-land-wave-ocean forecast model and an analysis system that provides an accurate estimate of the initial state. The forecast model has a 9 km horizontal resolution grid and 137 vertical levels, with a top around 80 km altitude. The analysis is based on a 4dimensional variational method, run twice daily using more than 25 million observations per cycle, primarily from satellites. The IFS produces high-resolution operational 10-day forecasts twice daily.
Radiative transfer calculations. The equinox-equivalent daily-average shortwave (integrated between 300 and 3000 nm) surface and top of the atmosphere (TOA) direct radiative forcing (RF) are estimated using the UVSPEC (UltraViolet SPECtrum) radiative transfer model in the libRadtran (library for Radiative transfer) implementation 36 and a similar methodology as in ref. 37 and ref. 4 Baseline and fire-perturbed simulations are carried out with different aerosol layers: the average OMPS-LP aerosols extinction coefficient profiles, for January and February 2019 (baseline simulation) and January and February 2020 (fire-perturbed simulation). The spectral variability of the aerosol extinction is modeled using the measured Ångström exponent from SAGE III for January 2020 (fireperturbed simulation) and typical background values inferred from SAGE III (baseline simulation). Different hypotheses have been considered for the nonmeasured optical parameters of fire aerosols: single scattering albedo from 0.85 to 0.95 (typical of wildfire aerosols, see e.g., ref. 38 ) and a Heyney-Greenstein phase function with an asymmetry parameter of 0.70. More information about the UVSPEC runs can be found in the Supplementary Material. The daily-average shortwave TOA radiative forcing for the fire-perturbed aerosol layer is calculated as the SZA-averaged upward diffuse irradiance for a baseline simulation without the investigated aerosols minus that with aerosols, integrated over the whole shortwave spectral range. The shortwave surface radiative forcing is calculated as the SZAaverage downward global (direct plus diffuse) irradiance with aerosols minus the baseline, integrated over the whole spectral range.
Mass estimation. The stratospheric masses of CO, CH 3 CN, and H 2 O are derived from the MLS species mixing ratios vertical profiles combined with MLS vertical profiles of pressure and temperature on the standard 37 pressure levels. MLS data are filtered according to the recommendations from the MLS team (https://mls.jpl. nasa.gov/data/v4-2_data_quality_document.pdf). CH 3 CN measurements are not recommended below 46 hPa however have already been used successfully in a study on combustion products from Australian bushfires in the stratosphere 39 . The mass calculation is performed as follows. First, we calculate for each profile the partial column of species (CO, CH 3 CN, H 2 O) and air for each measurement layer (vertical resolution of ∼3 km); the air partial column is derived from the difference in pressure between the top and bottom of the layer. Then, adding the partial columns in a profile, we derive the total column of species and air above the 380 K potential temperature (∼12. 17 km altitude), the so-called stratospheric "overworld" 13 . The ratio of the species and air total columns gives us the mean volume mixing ratios (VMR) of the species over the stratospheric profile. Finally, we calculate the mean VMR of all the profiles between 20°S and 82°S (which is the southernmost latitude of the MLS sampling) and multiply it by the molecular mass of the species and the total number of air molecules above 380 K and between 20°S and 82°S following ref. 40 to obtain the total mass burden of the species plotted in Fig. 4. The aerosol mass is derived from the OMPS satellite aerosol extinction data assuming a particle mass extinction coefficient of 4.5 m 2 g −1 following ref. 2 Standard deviations of the injected mass are estimated by combining accuracies on the measurements and the mean standard deviations over 20-day periods before and after the sharp increase. The standard deviation of the CH 3 CN mass above 380 K is only calculated from the standard deviation because accuracies on CH 3 CN measurements are extremely large. The standard deviation of the aerosol mass takes also into account the uncertainty on the particle mass extinction coefficient (error = 1.5 m 2 g −1 ).