Dynamical Perturbation of the Stratosphere by a Pyrocumulonimbus Injection of Carbonaceous Aerosols

.


Introduction
Wildfires release large quantities of water vapour and aerosols in the atmosphere, as well as carbon oxides and other gases (Wiedinmyer et al., 2011). Under favourable atmospheric conditions, the rising of plumes from the wildfire can reach the stratosphere, as shown by Peterson et al. (2018) and Fromm et al. (2000Fromm et al. ( , 2010. This lofting to the stratosphere is driven by different mechanisms. First, as soon as the plume reaches its lifting condensation level, the release of latent heat from the 20 condensation of water vapor results in its further buoyant rise. This process is referred to as pyroconvection and produces pyroCumulus (pyroCu) clouds. Then, if the state of the troposphere is conducive to the formation of dry thunderstorms (Peterson et al., 2017), pyroCu clouds can further develop vertically, reaching the tropopause and evolving into pyroCumulonimbus (pyroCb) clouds, i.e. wildfire-driven thunderstorms loaded with smoke from the wildfire (Fromm et al., 2010).
A distinctive feature of pyroCb clouds is the transport of wildfire by-products to the Upper Troposphere Lower Stratosphere 25 (UTLS): these thunderstorms act as chimneys, promoting the accumulation of aerosols and other gases in the UTLS. Once in the stratosphere, the main depletion mechanisms for aerosols (i.e. wet and dry deposition) are greatly reduced compared to the troposphere, and as a consequence aerosols have a significantly longer lifetime than in the troposphere (Yu et al., 2019). In the UTLS, the optically thick carbonaceous aerosols from wildfires absorb radiation, causing the temperature of the smoke plume to rise. This warming results in a gradual diabatic lofting of the plume itself (Torres et al., 2020;Das et al., 2021). This process 30 is more pronounced in the days following the aerosol injection into the UTLS, when the concentration of aerosol is higher (de Laat et al., 2012;Ditas et al., 2018).
The second largest pyroCb aerosol injection ever recorded is the 2017 Pacific Northwest Event (PNE). The PNE was characterized by an outbreak of seven pyrocumulonimbus clouds in British Columbia, Canada and the state of Washington, USA during the evening of the 12 th of August 2017 and the first hours of the following day (Peterson et al., 2018;Fromm et al., 35 2021). Overall, the pyroCb outbreak injected into the stratosphere an unprecedented smoke plume, whose mass was estimated to be between 0.1 Tg and 0.35 Tg (Peterson et al., 2018;Torres et al., 2020). This record was surpassed by the injection from the December 2019/January 2020 Australian New Year event (ANY), which resulted in the stratospheric injection of up to 1.1 Tg of aerosol .
Several studies (Das et al., 2021;Bourassa et al., 2019;Kloss et al., 2019;Torres et al., 2020;Yu et al., 2019;Christian et al., 40 2019; Baars et al., 2019) have characterized the aerosol plume from the PNE event and its interaction with large scale stratospheric features such as the Asiatic Summer Monsoon Anticyclone. In the days following the PNE, the smoke gradually rose into the stratosphere due to diabatic heating, with peak ascent rates of 2-3 km/day (Das et al., 2021;Torres et al., 2020). The plume eventually reached 22 km in height around 20 days after the injection, as detected in OMPS LP (Ozone Mapping and Profiling Suite, Limb Profiler) and SAGE III (Stratospheric Aerosol and Gas Experiment) observations (Torres et al., 2020; namical features of the PNE plume by inspecting ERA5 reanalysis fields and CALIOP (Cloud-Aerosol Lidar with Orthogonal 50 Polarization) measurements, finding that it evolved into stratospheric anticyclones that persisted for almost two months.

and have been named SWIRLs (Smoke With Induced Rotation and
Lofting) by Allen et al. (2020). These studies underline the stability of the SWIRLs and their resilience against the large scale shear characterizing the background wind field. The main signatures of SWIRLs are a deep potential vorticity anomaly with 55 respect to the zonal mean (negative in the northern hemisphere, positive in the southern) and an associated anticyclonic motion, enhanced optical thicknesses and a vertical temperature anomaly dipole. Also, the SWIRL encases air from the UTLS with low ozone content and this is reflected in a characteristic negative ozone concentration anomaly as the SWIRL moves upwards. As pointed out by Khaykin et al. (2020), the SWIRL effectively traps the carbonaceous aerosol, with the result of efficiently transporting it to higher altitudes. The cited studies underline how the presence of carbonaceous aerosol is the necessary condition 60 in the formation and maintenance of SWIRLs and how the coupling between aerosol, radiation, and dynamics is necessary to reproduce a SWIRL in model simulations (Khaykin et al., 2020;Allen et al., 2020).
In this work we use a chemistry climate model to simulate the impact on the stratosphere of a pyroCb plume from an event of the magnitude of the 2017 Pacific Northwest Event (PNE), as well as the effect of the aerosol radiative interaction in the development of the plume itself. As in Das et al. (2021), we use the Goddard Earth Observing System Chemistry-Climate 65 Model (GEOS CCM), focusing in this study on the GEOS Atmospheric General Circulation Model configuration with prognostic aerosols from the GOCART (Goddard Chemistry, Aerosol, Radiation, and Transport) module to simulate the long term transport of the plume and its impact on the radiative budget of the atmosphere. While Das et al. (2021) used GEOS CCM in a relaxed replay mode, i.e. driven by observations in the troposphere, we present here free running simulations (i.e. weather forecasts). In previous work by Khaykin et al. (2020), Allen et al. (2020) and Lestrelin et al. (2021), the structure of the SWIRL 70 was reproduced and maintained thanks to the assimilation of the observations and not by the local radiative heating caused by the carbonaceous aerosols. Here, we show how the GEOS CCM is capable of forming and maintaining a SWIRL following a stratospheric aerosol injection such as the one from the PNE event, thanks to its representation of the radiative impact of the aerosol on the dynamics.
The choice of the free running configuration is driven by the need to unambiguously resolve the dynamical and thermodynam-75 ical impact of the radiative heating by the aerosol. Indeed, the absence in the simulations of replay (or nudging) to a reanalysis ensures that any perturbation connected to the presence of the aerosol comes from the aerosol itself and not from underlying reanalysis fields or measurements. However, simulations in the free running configuration will not necessarily reproduce the real world smoke transport, because the simulated meteorological fields will differ from the observed ones. This impedes a thorough comparison between observed SWIRLs and the simulated one, but we show in Section 3 how a free running simula-80 tion still reproduce SWIRLs with realistic characteristics.
The article is structured as follows. In section 2 we briefly describe the GEOS model as well as the setup used for the simulation.
In section 3 we present the results from the simulations, analyzing the dynamical signatures and thermodynamical characteristics of the reproduced SWIRL and comparing them to previous works; more specifically, we will analyze the anomalies of potential vorticity, wind, temperature, and geopotential that describe the SWIRL as well as their evolution during the SWIRL 85 life; moreover, we will show how the diabatic heating provided by the presence of the aerosol is correlated to the dynamics of the SWIRL itself, and how the SWIRL undergoes a daily expansion/compression cycle following this diabatic heating.
Also, the geometrical characteristics of the SWIRL will be pointed out, and we will compare the resulting dimensions with the SWIRLs observed in previous works. In sections 4 and 5 we will present some considerations about the reproduced SWIRLs and draw some conclusions about the novelty and limitation of this work in comparison to other studies.

Model Description and simulation setup
The Goddard Earth Observing System (GEOS) can simulate atmospheric circulation and chemistry, oceanic circulation and biogeochemistry, land surface processes and data assimilation at horizontal resolutions as small as 12 km (Molod et al., 2015;Rienecker et al., 2008). In this study we use GEOS as an atmospheric general circulation model, resolving the atmospheric 95 circulation and composition but using prescribed sea surface temperatures and sea ice fractions. GEOS can be used in replay and free running configurations. In the replay mode, simulations are constrained by the MERRA-2 reanalysis fields, while in the free running mode the evolution of the atmosphere is not constrained to the observations. The results shown in this study were produced in free running simulations with the Icarus 3.3 version of GEOS. The meteorological fields for the model initialization on the 13 th of August (00Z) were obtained from version 2 of the Modern-Era Retrospective analysis for Research 100 and Applications (Gelaro et al., 2017).
Simulations were run on a cubed-sphere horizontal grid with horizontal resolution of ∼ 50 km with 72 hybrid vertical sigma levels extending from the surface to 0.01 hPa (Rienecker et al., 2008). Aerosol concentrations and aerosol processes are simulated with the Goddard Chemistry, Aerosol, Radiation, and Transport module (GOCART), a bulk aerosol module simulating the evolution of black carbon (BC), brown carbon (BrC), organic carbon (OC), nitrates (NO 3 ),sulfates (SO 4 ), dust, and sea salt 105 (Chin et al., 2002(Chin et al., , 2009Colarco et al., 2010Colarco et al., , 2017. In our simulations, aerosols are radiatively coupled to the dynamics and have a direct impact on the meteorological forecast. The optical properties of all the species, except dust and BrC are determined by the OPAC data set (Hess et al., 1998), while dust optics are treated as described in Colarco et al. (2014). As in Das et al. (2021), wildfire emissions are assigned to BC and BrC, which is more absorbing than OC in the near UV. The simulations also include the stratospheric chemistry module StratChem (Considine et al., 2000;Douglass and Kawa, 1999), which allows 110 for the inclusion of background stratospheric sulfate from the oxidation of carbonyl sulfide. Stratospheric chemical reaction rates are impacted by the pyroCb smoke through temperature changes, but not via changes in heterogeneous chemistry.
The PNE plume is represented as a 0.3 Tg injection of carbonaceous aerosol, which is consistent with the observational findings from Peterson et al. (2018) and Torres et al. (2020). Of this, 97.5 % is BrC and the remaining 2.5% is BC. The BC/BrC ratio was set by Das et al. (2021) to reproduce the observed diabatic lofting of the plume, which is sensible to the overall loading of 115 the strongly absorbing BC. The 0.3 Tg are injected in three different locations in British Columbia, in correspondence to the most pronounced three pyroCb clouds of the PNE. The aerosol have been initialized as smeared in regions centered on these locations, spanning horizontally 2°longitude by 2.5°latitude and vertically from 10 to 12 km. The aerosol injections occur during the first hours of the 13 th of August, with the first two puffs occurring from 0 to 3 UTC (0.19 Tg, locations 1A and 1B, Fig. 1) and the third from 4 UTC and 6 UTC (0.11 Tg, location 2, Fig. 1). These time intervals have been obtained using Cloud 120 Top brightness temperature measurements by Das et al. (2021). The modal radius of the BrC aerosol in the model was set to 0.035 µm to achieve an Angstrom Exponent of 1.3, following SAGE III measurements of the event (Das et al. (2021)). In this work we present a simulation with the PNE injection, starting on 13 th August 2017 and ending on 30 th September 2017, that shows the formation of a SWIRL.  125 We define the boundaries of the SWIRL based on the local anomaly of Ertel's potential vorticity (PV) (Ertel, 1942) with respect to the zonal mean at the same altitude and the brown carbon concentration. We identify the SWIRL as the region where the percent PV anomaly with respect to the PV zonal mean at the same altitude is smaller than -10% and the mass mixing ratio of brown carbon is higher than 1.25 µg kg −1 . These criteria constrain PV anomalies to those associated with the presence of high carbonaceous aerosol concentrations. In order to obtain smooth SWIRL contours, both the carbonaceous aerosol concentration 130 fields and the PV anomaly fields have been filtered horizontally using a gaussian filter with a standard deviation of 1 degree in the latitudinal direction and 1.5 in the longitudinal direction. We concentrate on model levels higher or equal to 100 hPa, so as to ensure that the detected SWIRL is above the tropopause. The analysis proposed here to track the SWIRL follows the approach of Lestrelin et al. (2021), which uses the Lait Potential Vorticity anomaly (Lait, 1994), along with the ozone concentration anomaly, to locate the SWIRL from ERA5 reanalysis fields. A similar approach was also used by Khaykin et al. (2020) 135 to track the SWIRL from the ANY event.

Definition of SWIRL boundaries and data analysis
To investigate the characteristics of the SWIRL, we analyzed the following diagnostic quantities: Ertel's potential vorticity (PV), temperature (T), potential temperature (θ), density (ρ), zonal wind speed (U), meridional wind speed (V), vertical wind speed (W) and carbonaceous aerosol (BC and BrC) concentrations. All these quantities are saved on model levels every 6 hours, and will be presented as absolute fields or as anomalies with respect to their zonal means at the same model level and 140 latitude unless specified otherwise. The geopotential of the model levels that is present in our discussion is approximated by the mid-layer geometrical height multiplied by the gravitational acceleration g. In our analysis we have also used the temperature tendencies due to the diabatic processes resolved by the model (radiation, moist physics, friction, turbulence, and gravity wave drag) and to the dynamics. Additionally, the model calculates the contribution of the aerosol on the temperature tendency due to radiation through two calls to the radiation model, with and without the inclusion of aerosols. In order to better capture the 145 relationship between the SWIRL and its closest surroundings, we restricted our analysis to a region centered on the SWIRL spanning 43.75 • longitude by 35 • latitude, and calculate the diagnostic variables as horizontal averages over the portions of this region inside and outside the SWIRL. The SWIRL geometry was described considering its volume, which is calculated by summing the volumes of the grid cells inside the SWIRL contours over the height of the SWIRL, and its average area, which is calculated by vertically averaging the areas of the SWIRL computed for every layer by summing the areas of the cells lower 150 boundaries inside the SWIRL selection. The thickness was obtained by considering the height difference between the highest and lowest levels at which the SWIRL was detected.
To visualize and describe the evolution of the bulk SWIRL properties, we considered averages and extremes over the SWIRL volume of our meteorological fields and temperature tendencies.

SWIRL Tracking
In the days following the PNE aerosol injection in the UTLS, the smoke plume drifted eastward and the diabatic self-lofting promoted its vertical movement. A brief description of the smoke transport is available in appendix A, where we also include a comparison of the plume transport in the free running and relaxed replay configurations of the GEOS model. This comparison 160 shows how the relaxed replay mode is a good proxy for the observed smoke transport, but has the drawback of damping the effects on the dynamics of the radiative heating of the aerosol. On the other hand, the free running simulation does not faithfully reproduce the observed smoke transport because the ambient meteorology in the free running simulation differs from the observed transport. In this work we use the free running configuration to resolve directly (i.e. with no contribution from reanalysis fields or measurements) the dynamical and thermodynamical impact of the aerosol radiative heating, which is crucial 165 in the formation of the SWIRL.
In our free running simulations the plume reached 100 hPa on the 16 th of August, consistent with OMPS LP observations (Das et al., 2021), marking the start of the SWIRL detection using the algorithm described in Sec. 2.2. By the end of the 18 th of August the detection algorithm revealed that the bulk of the SWIRL was above the 100 hPa limit, spanning between the levels 100 hPa and 61 hPa. For our analysis, we mark the 18th of August as the first day of the stratospheric SWIRL. Before this date 170 some signs of SWIRL formation were visible at levels lower than 100 hPa and above the simulated tropopause height (about 150-200 hPa in the days following the injection in the regions corresponding to the plume). We have chosen not to include this initial part in the analysis, since we are interested in the life of the SWIRL when it is well into the stratosphere.
The SWIRL is first detected by the algorithm over the Atlantic on the 18 th of August (Fig.2a) and traveled southeastward until it reached the northern Atlantic. Here it stalled for 7 days (days 6-12 after the PNE, 19 th -25 th August) before starting its 175 westward movement, which eventually brought it to the Arabic peninsula (day 28 after the PNE, 10 th of September). The next day (11 th of September, day 29), the SWIRL started losing compactness and that marked the end of its stratospheric lifetime.
After this date, the detection algorithm was still able to detect its fragmented remains until at least the 16 th of September, when they reached the Caribbean Sea. The contours of the fragmented SWIRL have not been included in Fig. 2a.
The simulated development of the SWIRL differs from the observed one in several aspects. First, the transport is different, Moreover, on the 1 st of September, vortex O split into two vortices (B1 and B2), whose trajectories are not reported here. This 185 splitting is not reproduced by the simulations. Lastly, the simulated SWIRL lasted overall 25 days, which is shorter than in observations Lestrelin et al. (2021), where vortex A remained visible until the 15 th of October. The differences between the simulated and observed development of the SWIRL is due to the underlying meteorology. The generation of offspring vortices and the development of stalling periods are caused by the interaction of the vortex with surrounding meteorological features, such as jets and troughs. Since the underlying meteorology is different in the simulation and observations, its influence on the 190 smoke plume and on the anticyclonic vortices will be different and we do not expect the evolution of the plume to be the same.
Additionally, differences in the amount and properties of the carbonaceous aerosols in the simulated vortex with respect to the observed ones cause differences in its radiative interaction and, therefore, in its evolution.

Analysis of the SWIRL on the 23 rd of August
In this section we will describe the structure of the simulated SWIRL on August 23 rd . On this day the vortex is well into the 195 stratosphere (detected in the model levels between 72 and 43 hPa), with tangential wind speed anomalies of about 10 ms −1 , and well organized, with a clear anticyclonic circulation closed around its center. Inside the borders of the SWIRL lies the deep negative PV anomaly, which indicates the presence of an anticyclone (Figs. 3a-d and 4a-d ). Indeed, the horizontal wind anomaly structure of the SWIRL shows an anticyclonic circulation centered over the PV anomaly minimum, with tangential wind magnitudes of about 10 ms −1 (Figs. 3q-t and 4q-t). The magnitude of the tangential wind speed in the SWIRL here 200 simulated is comparable to the ones presented in literature for the SWIRL developed after the Australian New Year event (Kablick III et al., 2020;Allen et al., 2020). The analysis of the geopotential shows a positive anomaly in correspondence of the SWIRL (Figs. 3q-t and 4q-t): this indicates that the effect of the aerosol is to trigger an outward pressure gradient force 8 Figure 3. Horizontal sections at different heights of the SWIRL on the 23 rd of August, 18z, when the maximum heating occurs. The first to the fourth columns show 43, 52, 61, and 72 hPa altitudes, respectively. The first row (a-d) shows the potential vorticity anomaly [%], the second (e-h) the vertical velocity, the third (i-l) the potential temperature anomaly, the fourth (m-p) the temperature tendency due to radiative heating/cooling of the aerosol, and the fifth row (q-t) the geopotential anomaly (color shades) and the wind field anomaly (white arrows). All anomalies are calculated with respect to the zonal mean. The geopotential anomaly was made positive in the considered region by summing the absolute value of the minimum of the geopotential anomaly to the entire field. This was done for representation purposes. The black contours represent the boundary of the detected SWIRL.
9 Figure 4. Same as Fig. 3, but for the 23 rd August 00 z, when the aerosol is cooling due to the lack of incoming solar radiation. that, together with Coriolis acceleration produces the anticyclonic motion of the SWIRL.
The pyroCb aerosol causes a strong internal heating of the SWIRL during daytime (Fig. 3m-p). During nighttime the aerosol 205 does not absorb SW radiation but emits LW radiation, inducing a net cooling of the plume (Fig 4m-p). Vertical velocities within the SWIRL are comparable to those outside the SWIRL (Fig 4e-h). From this analysis there is no evidence of an organized structure in the vertical wind field in correspondence of the SWIRL and there is no sign of its upward movement. This however is expected considering that the vertical displacement of the SWIRL, which can be quantified in hundreds of meters per day, should be driven by a vertical velocity of the order of fractions of cm s −1 , which cannot be visible in the noisy vertical velocity 210 background. The lofting of the SWIRL will be made evident later on, considering the average of the vertical velocity in the SWIRL (Fig. 7).
While the higher levels of the SWIRL are approximately elliptic, the lowermost levels present a peculiar shape, consisting of an ellipsoid to which a protrusion is attached in the downwind direction. This protrusion can be seen as the tail of the SWIRL, which gradually loses material in time, leaving a trail of air with negative PV anomaly and high aerosol concentration that is 215 classified as SWIRL by our detection method. Such tails have been also observed by Lestrelin et al. (2021).
Another characteristic feature of SWIRLs is the vertical potential temperature anomaly dipole. This is visible considering the maps of potential temperature anomaly (Figs. 3i-l and 4i-l), which indicate a negative anomaly above the center of the SWIRL and a positive anomaly below. As stated by Allen et al. (2020), this signature in the potential temperature is consistent with the vertical expansion of the plume due to the diabatic heating of the aerosol.

220
In the vertical sections of the SWIRL (Figs. 5-6e) a strong anticyclonic motion is visible as well as the vertical structure of the diabatic heating provided by the aerosol (Fig. 5-6c), which is particularly intense in the middle of the SWIRL. The analysis of the meridional wind anomaly curtain (Fig. 6e) reveals also the tilting of the SWIRL, whose axis lies at an angle with respect to the vertical. As described by Allen et al. (2020), the tilting of the vortex is given by the wind shear in the vertical direction.

225
The characteristic vertical temperature and potential temperature anomaly dipoles that characterize SWIRLs are also visible in Figs. 5b-f and 6b-f, with negative anomalies in the upper part of the SWIRL and positive anomalies below. The magnitude of the temperature dipole (10 K) is comparable to that reported for the ANY case (Allen et al., 2020), which reached up to 15 K.
The temperature and potential temperature anomalies extend well below and above the potential vorticity anomaly, as is also reported in Allen et al. (2020); Lestrelin et al. (2021) and Kablick III et al. (2020). As in the horizontal sections, the vertical 230 velocity field does not exhibit any detectable difference between the interior of the SWIRL and its surroundings. The temperature anomaly dipole showed in Figs. 5b and 6b is reflected in a reduced lapse rate inside the SWIRL compared to the surroundings (Fig. 7a), with the internal temperature profile intersecting the external one roughly in the middle of the SWIRL. The same can be observed in the potential temperature anomaly (Fig. 7d). Fig. 7h shows the diurnal heating and 235 nocturnal cooling of the SWIRL due to aerosols and and Fig. 7b its resulting vertical motion. The difference in the average vertical velocity between the interior and exterior of the SWIRL is crucial, since it explains how the SWIRL actually moves upwards despite being immersed in a noisy vertical velocity field.
The profile of the horizontal wind speed anomaly (Fig. 7e) shows stronger wind speeds in the middle part of the SWIRL; this was visible also in Fig 3. The reason for this is that the concentrated heating that takes place inside the SWIRL produces a 240 strong rotation which does not interact with the background flow. Indeed, the similarly strong heating taking place at the bottom of the SWIRL (Fig. 7h) does not result in a similarly strong rotation due to the SWIRL interaction with the background flow.
This fact is also supported by the shape of the SWIRL, which presents tails in its lower levels caused by the dispersive action of the background flow.
The radiative heating/cooling due to the aerosol is the main factor among the physical diabatic processes (Fig. 7g-i). During 245 the day, the SWIRL is subject to strong diabatic heating, that promotes the vertical motion from buoyancy. The vertical motion is accompanied by the expansion of the SWIRL, which partially compensates for this heating; this is evident in the temperature tendency due to the dynamics ( Fig. 7b and i, red and dark red curves). During the night, the situation is the opposite. The diabatic cooling drives a downward vertical motion of the SWIRL and its subsequent heating due to expansion ( Fig. 7b and i, yellow and orange curves). The nocturnal cooling is less intense in absolute value than the diurnal heating. This points toward 250 a net ascent of the SWIRL, given the enhanced intensity of the diurnal lofting with respect to the nocturnal descent (Fig. 7h). (c) Ertel's potential vorticity, (d) potential temperature anomaly with respect to the zonal mean, (e) horizontal speed anomaly with respect to the zonal mean, (f) density anomaly, (g) temperature tendency due to diabatic processes, (h) temperature tendency due to the aerosol, and (i) temperature tendency due to the dynamics. 14 3.3 Evolution in time of the SWIRL properties aerosol loading and (f) absolute temperature. In (d) the red line shows the result of the linear fit between day 5 and 11 of the potential temperature, while the blue line between day 11 and day 30. The steepness of the red line is 10.6 K day −1 , while for the blue line it is 6.2 K day −1 .The green dots and the blue triangles represent the potential temperature of vortices O and A respectively, obtained using ERA5 meteorological fields in Lestrelin et al. (2021). The data were kindly provided by Dr. Bernard Legras.
The volume of the SWIRL increases until day 10-11 (22 nd -23 rd August ) and then decreases steadily in time (Fig. 8a). The daily modulation of the volume, corresponding to the expansion/compression cycle caused by the diurnal radiative cycle, is 255 visible especially in the first 11-12 days. Accordingly, the calculated radius (Fig. 8b) shows the same time evolution, peaking around day 11 at 700 km. The peak radius of the simulated SWIRL is comparable to the maximum radius of the SWIRL observed following ANY (peaking at 750 km, Allen et al. (2020)); instead, the simulated SWIRL is consistently larger in comparison to the main SWIRL observed following PNE, whose horizontal dimensions (diameter) where estimated at 700 km from composite analysisLestrelin et al. (2021). After day 28 (11 th September) the SWIRL volume increases significantly, following its dispersion and fragmentation. The detection algorithm still categorizes the remains of the plume as SWIRL even thought they lack the characteristic compactness of SWIRLs. This process is also visible in the calculated average radius, with a steep increase after day 28. Moreover, the SWIRL thickness (Fig. 8c) indicates that the SWIRL reaches a maximum thickness of 5 km between day 11 (24 th August) and 16 (29 th August); before and after this period, the thickness value fluctuates between 3 and 4 km; the magnitude of the SWIRL thickness is comparable to the one observed in the ANY SWIRL, that peaked at 6 km 265 in the moment of its maximum expansion (Allen et al., 2020), and it is also close to the thickness of 3 km of the PNE SWIRL obtained from the composite analysis of Lestrelin et al. (2021). The step-like fluctuations of the SWIRL thickness is due to the vertical resolution of the model at the considered levels, which is around 1 km. The ascent of the simulated SWIRL can be divided in two phases: the steep initial ascent, during which the volume of the SWIRL increases, and a second slower ascent phase, characterized by a decrease in volume of the SWIRL. The potential temperature θ increases approximately linearly (Fig.   270 8d) during both ascent phases, with a slope of about 10 K day −1 during the first phase and of 6 K day −1 during the second.
The ascent speeds are consistent with the observed ones during the ANY SWIRL: in Allen et al. (2020) the ascent speed during the first 27 days of the SWIRL life is estimated at 7.8 K day −1 and 5.5 K day −1 afterwards. The potential temperature of the simulated vortex is about 50 K higher than in Lestrelin et al. (2021) (Fig. 8 (d)), possibly because of the faster and more consistent simulated ascent of the plume in the first days following the PNE (see also appendix A). Despite this offset, between 275 days 5-29 from the PNE the ascent rates of the simulated SWIRL are compatible to the one of vortex O observed by Lestrelin et al. (2021). In the first 11 days of its life, vortex O has an ascent rate of about 9 K day −1 , while following its first splitting (formation of vortex A) it rises in the stratosphere at about 4.7 K day −1 . These estimates have been obtained with linear fits of the potential temperature curve of vortex O. Before its dissipation 30 days after the PNE, the simulated SWIRL reaches an average potential temperature of 590 K, comparable to the final potential temperature reached by Vortex A (570 K) before its 280 loss 60 days after the PNE Lestrelin et al. (2021). Lastly, the average absolute temperature T within the SWIRL (Fig. 8f) shows little variation in time, with values between 206 and 209 K.
The magnitude of the temperature anomaly dipole in the SWIRL is approximately stable until day 10-11 (22 nd -23 rd August), and then gradually decreases as the SWIRL ages (Fig. 9a). The simulated temperature dipole (peaking at 7 K) is significantly more pronounced than the one calculated for the main PNE SWIRL via composite analysis by Lestrelin et al. (2021) (peaking 285 at 2 K). The maximum horizontal wind speed anomaly in the SWIRL (Fig. 9b) first decreases for three days, in the same period in which the radius of the SWIRL increases significantly; then it stabilizes around 12-13 m s −1 before decreasing again after day 15 (28 th August) and reaching figures between 6-8 m s −1 . The same behaviour is observed in the average horizontal velocity anomaly that decreases in the first days from 10 m s −1 to 5 m s −1 , then remains constant at 5 m s −1 until day 15 and subsequently decreases gradually reaching 3 m s −1 at the end of the SWIRL lifetime.

290
The minimum of the PV anomaly (Fig. 9c) abruptly increases at day 6-7 from about -140% to about -90% and remains roughly stable until day 11, when it starts increasing at a slower rate. The average vertical velocity of the SWIRL (Fig.   9d) is positive (indicating upwelling) at the end of the daily heating period and negative during the night, with maximum values peaking at 0.7 cm s −1 . The average vertical velocity of the SWIRL reaches more pronounced peaks than those of the surroundings, indicated by a green line in Fig. 9d.

295
The diabatic temperature tendency (Fig. 9e) is initially dominated by the radiative heating of the aerosol. The mean aerosol heating rate within the SWIRL is significantly larger than the average background diabatic temperature tendency until day 22-23 (4 th -5 th September). After this day, the diabatic temperature tendency inside the SWIRL becomes progressively closer to the diabatic temperature tendency outside. The quick loss of coherence of the SWIRL that takes place after day 28 is driven by the loss of this differential heating between its inner region and the surroundings. The dynamical temperature tendency is 300 anticorrelated with the diabatic temperature tendency (Fig. 9f), which is dominated by the aerosol heating, and is directly linked to the vertical velocity of the SWIRL: when the SWIRL is heated by radiation, it rises and expands, meaning that its average vertical velocity is positive and the dynamical temperature tendency is negative. During the night this tendency is reversed along with an average negative vertical velocity. This means that while the SWIRL is cooling, it moves downward while being subject to compression.

4 Discussion and Limitations
Overall, GEOS simulates a SWIRL with intensity and characteristics similar to what was observed after the ANY and PNE events (Lestrelin et al., 2021;Allen et al., 2020). The GEOS-simulated SWIRL has a shorter lifetime than in observations and produces a different trajectory. A key point in maintaining the SWIRL appears to be the difference between the diabatic heating within the SWIRL and the surroundings. This depends on the SWIRL's aerosol loading and its radiative properties, as well as 310 the local meteorology and the location of the vortex. Therefore, a role in determining the development of the SWIRL must be played by the specific meteorological conditions after the event, which are not expected to be reproduced in a free-running simulation. A whole ensemble of free-running simulations should be used to analyze the importance of specific background meteorological conditions on the formation and dissipation of the SWIRL, as well as the interaction of aerosol heating and background conditions in the initial and final stage of the SWIRL lifetime.

315
Despite being able to reproduce many characteristics observed in SWIRLs, the configuration of the GEOS model might impact the results. First, the vertical resolutions at the stratospheric altitudes here analyzed is about 1 km. A finer resolution could impact the vertical ascent, as clear in the step-wise behaviour of the SWIRL thickness visible in Fig. 8. Additionally, the GOCART aerosol module includes a parameterization of the trasformation of carbonaceous particles from hydrophobic to hydrophylic and the subsequent hygroscopic growth of hygrophilic particles, but does not include changes in optical properties 320 due to coagulation of aerosols, condensation of gaseous material, or the transfer from external to internal mixture. Lastly, GOCART assumes that black and brown carbon aerosols are spherical. Changes in depolarization ratio (both in space and time) as shown in Christian et al. (2019), or fractal structures such as in Yu et al. (2019) are not simulated. Since (Yu et al., 2019) showed that fractal structures are more absorbing at visible wavelengths, the radiative impact that we simulate might be underestimated. This might explain the shorter SWIRL lifetime in our simulations than in observations.

Conclusions
In this work, we present how the GEOS CCM is able to reproduce a SWIRL, i.e. a stratospheric aerosol-induced anticyclone.
Previous studies used reanalyses and observations to characterize the SWIRLs generated by the PNE (Lestrelin et al., 2021) and ANY (Allen et al., 2020;Kablick III et al., 2020;Khaykin et al., 2020) events. The simulations presented here reproduce several of their findings, and in particular the structures of the anomalies characterizing the SWIRLs. The use of a free-running 330 model with aerosol-radiation coupling allows us for the first time to resolve the role that the aerosol radiative interaction plays in the development of the SWIRL, as well as to characterize the diurnal cycle of the SWIRL and associated transport effects.
The presence of the heating of the aerosol is crucial, and Khaykin et al. (2020) showed that, without this contribution, SWIRLs dissipate in around 6-7 days.
In this work, the simulated SWIRL has been first recognized by using a simple detection algorithm, based on the Ertel's PV 335 anomaly and on the brown carbon concentration. The analysis of the SWIRL has been carried out by considering its structure on the 23 rd of August and then by evaluating how its bulk properties evolved during its lifetime. This analysis shows that the simulated anticyclonic disturbance reproduces the magnitudes of the anomalies of the observed ones; for instance, the magnitude of the flow of the anticyclone is about 10 m s −1 , which is close to the magnitude the flow in the SWIRL observed following ANY (Allen et al., 2020). The same is true for the magnitude of the temperature anomaly dipole, which is observed 340 together with a steeper lapse rate inside the SWIRL than in the surroundings. We also showed that the anticyclonic circulation is maintained by a pressure gradient force, triggered by a positive anomaly of the geopotential of the isobaric levels. This enhanced pressure in the center of the SWIRL is given by the diabatic heating of the aerosol, that also leads the expansion/compression cycle that the SWIRL undergoes and its consequent net diabatic lofting.
The diurnal cycle of the SWIRL is modulated by the daily radiation cycle. During daytime, the aerosol heats up the atmosphere 345 through the absorption of solar radiation. This triggers a radial expansion of the plume and a consequent vertical motion due to buoyancy. The expansion of the plume is also accompained by its dynamical cooling, which effectively counteracts the diabatic heating provided by the aerosol; as a result, the overall process is almost isothermal.
During the night, when the plume cools down radiatively, the process is inverted, and the SWIRL moves downward. The upward daytime motion is more pronounced than the nighttime one, leading to a net ascension of the plume. The analysis of the 350 evolution of the radius/volume of the SWIRL during its lifetime reveals this daily expansion-compression cycle of the SWIRL.
In conclusion, this work shows that free running simulations with a chemistry climate model such as GEOS are suitable to study the formation and development of stratospheric vortices following large injections of carbonaceous aerosols in the UTLS by pyrocumulonibus clouds. Indeed, despite the limitations of the configuration and of the models itself, we were able to 355 exploit them to reproduce such an event and to study several dynamical and thermodynamical features of SWIRLs. Several future developments of this work are possible; for example, model simulations with higher vertical resolution could help in determining more accurately the vertical structure of the vortex. Also, an ensemble of simulations such as the ones proposed here could help in better understanding the mechanisms behind the generation, maintenance and collapse of the vortex.
Data availability. The dataset used to obtain the results presented in this article are available for download in the online repository "GEOS CCM free-running simulation data of the Pacific-Northwest pyrocumulonimbus Event-like aerosol injection, SWIRL selection", url: https: //doi.org/10.5281/zenodo.6366106, (Doglioni, 2022).

20
Appendix: Aerosol Transport 365 Figure A1. Comparison of the zonal mean stratospheric aerosol optical thickness (SAOD) at 550 nm from the free running and replay simulations for four different days in the month following the PNE.
We compare here the transport of the aerosol in the free running and replay configurations of the GEOS model. The replay simulations here presented are thoroughly described in Das et al. (2021), and were performed using a relaxed replay configuration.
Following the injection from the PNE, the aerosol released in UTLS started its journey in the stratosphere. The diabatic self lofting promoted the gradual ascension of the plume while it was advected and dispersed by the large scale circulation.

370
The horizontal transport of aerosol is visualized considering maps of SAOD (Stratospheric Aerosol Optical Depth) at 550 nm (Fig. A1). The SAOD is calculated summing the simulated extinction coefficients fields due to the presence of Black and Brown carbon, and vertically integrating them from the tropopause up to the model's top of the atmosphere. The computation of the SAOD has been carried out both for the free running and replay simulations. As shown in Das et al. (2021), the latter well capture the aerosol transport as observed by OMPS LP; thus, the replay simulation is used here as a reference of the actual 375 aerosol transport.
In the first 2-3 days following the injection, the aerosol is advected and gradually dispersed by the wind field in the UTLS, that moves the plume towards the Atlantic both in the free running and replay simulations ( Fig. A1 (a)); at this stage the difference is limited between the different simulations.
As time passes, the free running simulation begins to significantly differ from the one in replay configuration( fig. A1, (d)). Of 380 particular interest during this phase is the splitting of a portion of the plume that can be observed in the free running simulation around day 11 (24th August 2017) above the northern Atlantic ocean ( fig. A1, compare (b) with (c)). After this event, a section of the plume reaches lower latitudes and starts moving towards the Gulf of Mexico.
This southern part of the aerosol plume exhibits interesting properties such as the confinement of the plume itself and its resilience against the disruption operated by the background flow. This portion of the plume generated the SWIRL, which is 385 presented in this article. The splitting of the plume does not occur in the replay configuration, in which the bulk of the aerosol reaches Europe. In order to compare the vertical displacement in time of the plume, we consider the time evolution of the horizontal average over the northern hemisphere of the extinction coefficient due to carbonaceous aerosols at 550 nm (Fig. A2), which reveals a fast 22 ascent of the aerosol plume in the first days following the injection. Comparing the free running with the replay configurations it 390 is evident that in the former the vertical transport is more pronounced: in the free running simulation the aerosol that reaches 100 hPa in 5 days is consistently more than in the replay simulation. Only from days 9-10 we find comparable extinction coefficients above 100 hPa in the replay simulation. This difference in the ascent rate in the replay and free running configurations could be given by the different background meteorology. Also, in the free running configuration the model's physics is free from the observational constraints, so that the aerosol diabatic heating can have a direct impact on it; instead in the replay configuration 395 the action of the aerosol in the dynamics might be dampened by the replay procedure itself. However, caution is needed when exploring this argument, since we have a single free running simulation; the study of an ensemble of free running simulations might be the correct instrument to better quantify and explain the differences in the diabatic self-lofting between the replay and free running simulations.
Author contributions. PRC and SD designed the modeling approach and performed the simulations. GD performed the data analysis and 400 analyzed the results under the supervision of VA and DZ. GD wrote the manuscript with VA and DZ. All the authors reviewed the manuscript before the submission.
Competing interests. The authors declare that they have no conflict of interest.