The evolution and dynamics of the Hunga Tonga-Hunga Ha’apai sulphate aerosol plume in the stratosphere

. We use a combination of space-borne instruments to study the unprecedented stratospheric plume after the Tonga eruption of 15 January 2022. The aerosol plume was initially formed of two clouds at 30 and 28 km mostly composed of sub-micronic sulphate particles, without ashes washed-out within the first day following the eruption. The large amount of injected water vapour led to a fast conversion of SO 2 to sulphate aerosols and induced a descent of the plume to 24-26 km over the first three weeks by radiative cooling. Whereas SO 2 has returned to background levels by the end of January, volcanic sulphates 5 and water still persisted after 6 months, mainly confined between 35 ◦ S and 20 ◦ N until June due to the zonal symmetry of the summer stratospheric circulation at 22-26 km . Sulphate particles, undergoing hygroscopic growth and coagulation, sediment and gradually separate from the moisture anomaly entrained in the ascending branch Brewer-Dobson circulation. Sulphate aerosol optical depths derived from the IASI infrared sounder show that during the first two months the aerosol plume was not simply diluted and dispersed passively but rather organized in concentrated patches. Space-borne lidar winds suggest that 10 those structures, generated by shear-induced instabilities, are associated with vorticity anomalies that may have enhanced the duration and impact of the plume.


Introduction
The phreato-magmatic eruption of the Hunga Tonga-Hunga Ha'apai (hereafter HTHH) on 15 January 2022 was exceptional in several respects. Its explosive intensity is close to that of the eruption of Mount Pinatubo in 1991, with a Volcanic Explosivity 15 Index of ∼6 (Poli and Shapiro, 2022). The induced atmospheric Lamb wave circled the globe at least 4 times with an amplitude comparable to that of the 1883 Krakatau eruption (Matoza et al., 2022;Vergoz et al., 2022;Wright et al., 2022). Within a few hours, several successive events injected material up to the mesosphere, with the bulk of the plume being detrained 1 between 26 and 34 km (Carr et al., 2022;Khaykin et al., 2022;Podglajen et al., 2022;Proud et al., 2022;Taha et al., 2022).
A further remarkable fact is that the plume carried an unprecedented amount of water vapour into the stratosphere, increasing 20 instantaneously its overall water vapour content by ∼10% (Millán et al., 2022;Khaykin et al., 2022). Quite surprisingly, the satellite data gathered after the event reported a stratospheric SO 2 injection of only 0.5 Tg, on par with much smaller and less explosive eruptions (Millán et al., 2022;Carn et al., 2022). This led to an early estimate of negligible climatic impact (Witze, 2022;Zhang et al., 2022). Here, we report on the evolution of the stratospheric plume during the first months after the eruption and we advocate that, due to the amount of water vapour and of the sulphate aerosols which have resulted from 25 a fast conversion, its climatic effect is very significant. We focus on the circumnavigation of the plume and proceed from the large-scale to the local patterns.
2 Six-month evolution of the zonal mean  Fig. 1a). The diabatic heating 30 rate is positive everywhere except at a narrow region near 27 km over the equator (Fig. 1b). These conditions are stable during the whole January-March period . In April-June, the angular speed weakens and changes sign (Fig. 1d&f) while the warming turns to cooling (Fig. 1e&g) as a combined effects of the Quasi Biennial Oscillation (QBO) and the seasonal cycle. Figure 2 shows that after an initial rapid meridional dispersion in the first days after the eruption , the aerosol and water plume stay mostly confined within the latitude band 35 • S-20 • N until June when wave activity increases, 35 and evolve slowly in the zonal mean. By mid-February, aerosols and water vapour have already spread through all longitudes (Khaykin et al., 2022, , and Fig. 6 below). The OMPS-LP extinction ratio increases in time at the core of the cloud and reaches a maximum in mid-April. The simultaneous increase of the CALIOP scattering ratio suggests particle growth. The MLS water vapour distribution initially coincides with that of the aerosols but they progressively move apart vertically (see also Schoeberl et al., 2022) due to the sedimentation of the aerosols. By early June, aerosols and moisture appear fully separated, respectively 40 below and above 25 km.
The plume vertical motion (Appendix A3) is analysed (Fig. 3) for two latitude bands and the apparent aerosol radius is estimated by interpreting the aerosol plume motion as a fall speed of the scattering particles. The descent of aerosols separates in two subsequent phases. For the first phase, lasting until about 20 February, Figure 3a-d show a fast descent in the two latitude bands which would imply unrealistically large aerosol sizes (Fig. 3e). During this phase, the water vapour follows the aerosol 45 downward motion. Sellitto et al. (2022) (S2022 hereafter) explain this behaviour by the cooling effect of water vapour infrared emission which is strong as long as the water vapour is concentrated and is located well above its neutral radiative level.
Sedimentation is then a secondary effect.
In the second phase, after 20 February, the diluted water vapour is ascending but is still producing a cooling detected by the departure of its ascent rate from the ERA5 ascent rate (Appendix A2), especially in the 25 • S-15 • S band, which persists 50 until June (see also Schoeberl et al., 2022;Coy et al., 2022) where it eventually vanishes as water vapour gets too diluted to produce significant cooling. The estimated sedimentation rate of the scattering aerosols and the corresponding particle radius (Appendix A3) suggest (Fig. 3e) that the aerosol particle size grows up to about 1.4 µm in April-May and starts shrinking in May-June in the 15 • S-5 • S band whereas it stays near 1 µm until June before shrinking in the 25 • S-15 • S band.
The extinction-to-backscatter ratio, obtained by combining OMPS-LP and CALIOP data (Fig. 3g) exhibits a growth followed 55 by a decay which are qualitatively consistent with the aerosol size evolution (Fig. 3g) and the expected behavior of the ratio ( Fig. 3f) in the 1-2, µm size range. Considering the saturation and decay of the extinction (Fig. 2), and the progressive vertical separation of aerosols and moisture, we suggest that the initial growth of the particles was by hygroscopic growth until April, where the extinction culminates, and was followed by coagulation over April-May and then by a decay due to evaporation as the ambient air gets drier and the aerosol plume is diluted. Coagulation and evaporation are obviously not exclusive and their 60 competition depends on the ambient conditions that vary over space and time (Hamill et al., 1977). It is also apparent from Fig. 2 that the moist layer is less confined than the aerosol layer and extends in latitude beyond the limits of the figure. The extinction-to-backscatter ratio is also smaller on the periphery of the aerosol plume (Fig. 3g). Therefore, we expect evaporation of the transported sulphate aerosols to occur at such latitudes.
3 Inferred composition of the plume 65 We now consider the history of the aerosol composition of the plume. The sequence in Fig. 4a-d shows, in agreement with Carr et al. (2022), S2022 and Khaykin et al. (2022), that the ash and ice cloud (brown and deep blue) is rapidly removed within the first day following the eruption likely via sedimentation of large ice particles which condensed water in excess of the saturated profile up to 35 km on 15 January . Taha et al. (2022) mention that ashes are missing in the plume on 17 January from UV satellite observations. What emerges on the west side are two greenish clouds (C1 and C2 70 on Fig. 4b-d) without any hint of ash (ashes would appear as yellow/reddish). The CALIOP cross section through these clouds ( Fig. 4e-f) shows high-scattering-ratio without depolarization, hence indicative of predominantly small spherical particles. The two clouds C1 and C2 are well separated in altitude. A few days later, the Light Optical Aerosol Counter flight and ground lidar observations, both from La Réunion, confirm this by showing sub-micron size, mainly non-absorbing, particles Baron et al., 2022).

75
A further source of information is from the Infrared/Microwave Sounder (IMS) retrieval (Appendix A1.2) of SO 2 column and Sulphate Aerosols (SA) optical depth. Figure 4g-h shows that the conversion to sulphates started immediately after the eruption with an SA optical depth reaching 0.1 one day after the eruption suggesting that the two clouds seen by CALIOP are composed of almost pure sulphate droplets. The fast conversion of SO 2 to sulphate aerosols is also discussed by S2022 and Zhu et al. (2022), using observations and chemical/transport modelling, respectively. The presence of a significant amount of 80 gas sulphates is not expected under the ambient conditions of the plume (Hamill et al., 1977;Tsagkogeorgas et al., 2017).
Four days later (Fig. 5b), the two clouds are still separated but have elongated under the zonal shear forming a pair of long strips. Comparing Figs. 5a and b, makes apparent that the conversion to sulphates is almost complete in the western strip generated from C1 while it is incomplete in the eastern strip generated from C2. S2022 show that the western cloud C1 is much moister than the eastern cloud C2, offering a likely reason for faster conversion, as also discussed by Zhu et al. (2022). A cloud 85 C3 of almost pure SO 2 is located between Australia and Indonesia (Fig. 5a), at lower altitudes than the other two clouds, as inferred from its low traveling angular speed. Comparing the IMS products to RGB-Ash ( Fig. 5c) demonstrates that RGB-Ash shows sulphates rather than SO 2 as usually assumed since both C1 and C2 are present but C3 is absent. The sensitivity of geostationary broad-band products, like RGB-Ash, to sulphates is shown by Sellitto and Legras (2016).
The conversion of remaining SO 2 to sulphates proceeds until SO 2 returns to background conditions by late January (Fig. 5d).

90
The sulphates persist for at least six months (Fig. 2) and the comparison of Fig. 5e and f shows that zonal averages of IMS and CALIOP products exhibit very similar patterns. The CALIOP depolarization ratio never exceeds its initial value (Fig. 4f) until July. Figure 6 shows the circumnavigation of the sulphate plume from a series of SA optical depth maps over one month and half (an 95 extended view until 30 April is provided by the supplement movie (Appendix B)). Due to the differential rotation, the fastest patches near 5 • S caught the slowest by 30 • S by mid-February and the plume filled the whole latitude circle. As time proceeds the components of the plume kept elongating and mixed together towards a zonal uniformity (see movie).

Circumnavigation and instabilities
However, Fig. 6 shows a number of localized concentrated patches which persist and keep forming in the plume one month after the eruption. Figure 7 investigates the structure of some of them and compares the SA optical depth to the observations  A remarkable feature in the SA optical depth maps is the train of compact elliptical structures linked together by filaments which is visible all along February and early March in Fig. 6 and the supplemental movie. This peculiar shape is reminiscent of shear-induced instabilities (Juckes, 1995), leading to the formation of a chain of vortices. The suspicion is reinforced by the pattern of a wrapping-up tripolar structure seen near 180 • E on 11 February (Fig. 7g) and perfectly captured by CALIOP as a core surrounded by two arms at the same level (Fig. 7h). This comparison also reveals the ability of the IMS product to retrieve 115 small-scale details.
Barotropic shear instability requires a reversal of the meridional gradient of absolute vorticity. The mean flow in ERA5 does not satisfy this criterion. A generalized baroclinic instability requires a reversal of the potential vorticity gradient but the mean flow again hardly satisfies this criterion at the required altitude of 25 km (Fig. 8). The very fact that the instability produces aerosol patches suggests that they are related to the generation of vorticity. The detection of an anomalous anticyclonic shear 120 across the concentrated patches of the plume by ALADIN (Fig. 7b-c) supports this hypothesis. However, sulphates are poor absorbers and neither these vortical structures nor their thermal signature have been detected by our present investigation of the ERA5. Therefore, this observation still requires an explanation that we leave for future studies.

Discussion and conclusion
The very intense and unusual HTHH eruption generated an intense and unusual stratospheric plume with a huge amount of 125 injected water vapour that remains well above normal 6 months after the eruption. After a fast initial removal of ice and ashes, the bulk of the remaining plume consisted of two main clouds between 26 and 32 km traveling westward due to the prevailing phase of the QBO. The ensuing zonal transport dispersed the plume through all longitudes in less than a month (see also Khaykin et al., 2022). The initial SO 2 is fully converted into sulphates in less than two weeks under the influence of water vapour. Notice that the absence of detection does not entirely rule out the possible presence of very thin ash as nucleous within 130 sulphate liquid droplet without optical signature.
The fast initial descent of the upper part of the plume induced by the radiative water vapour cooling has concentrated the aerosols within a fairly narrow layer, about 2 km thick as seen from CALIOP measurements (Figs. 2&5). Within the limit of MLS resolution, the water vapour distribution then coincides with the aerosols. The aerosols later continued subsiding at a slower rate under the effect of gravitational sedimentation, whereas the moist layer entrained by the Brewer-Dobson circulation 135 was simultaneously ascending, so that the two layers progressively separated (as also seen by Schoeberl et al., 2022). The spurious warming in the ERA5 that overlaps the moist layer suggests that radiative cooling by water vapour persists until May (see also Schoeberl et al., 2022;Coy et al., 2022). Although a precise sequencing is difficult without quantitative modelling, it is likely that the sulphate aerosols first grew by hygroscopic growth, then by coagulation and ended by dwindling under evaporation. Our estimation of fall speed and extinction-to-backscatter ratio trends is consistent with a growth up to about 140 1.4 µm and then a decrease in mean size.
The fast conversion of SO 2 suggests that the initial sulphur injection might have been underestimated. Consistently, S2022 showed that the HTHH eruption produced the largest stratospheric aerosol perturbation since the Pinatubo eruption in 1991, and suggested a large potential for climatic impacts (see also Khaykin et al., 2022). By June, the hemispheric stratospheric aerosol optical depth perturbation of the HTHH plume is twice as large as the peak perturbation of the 2019 Raikoke eruption, 145 and the tropical impact is at least three times as large as any volcanic perturbation since Pinatubo 1991 (S2022, Khaykin et al., 2022). As the SO 2 emissions for the Raikoke eruption have been estimated at 1.5 Tg (de Leeuw et al., 2021), we assume this value as the lower limit for the HTHH eruption, three times larger than early estimates (Witze, 2022). The young aerosols seem mostly made of sub-micronic liquid sulphate particles then growing to 1.4 µm due to hygroscopic growth and coagulation. The dispersion of the plume questions the magnitude and the duration of the impact. An early estimate of the resulting radiative 150 forcing by S2022 shows that stratospheric aerosol and water vapour perturbations from the eruption may significantly impact the climate system. Given the large greenhouse potential of stratospheric water vapour (e.g. Solomon et al., 2010), it was proposed that the dispersed plume has a net warming effect (S2022), in contrast with the cooling expected from stratospheric aerosols.
Finally, we have shown that the dynamics repetitively generates compact aerosol structures in a process that bears similarities 155 with shear instability and that some structures carry anomalous anticyclonic vorticity. That points to the possible role of such processes in extending the life time and the impact of the plume.
Appendix A: Data and methods

A1 Observations
We use data from the following instruments and products.

A1.1 CALIOP
The Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP) is a spaceborne lidar onboard the CALIPSO satellite (Vaughan et al., 2004;Winker et al., 2010). We use the L1 532 nm attenuated backscatter which is filtered in the horizontal with a median filter of width 102 km. In particular, this filter removes the noise associated with the South Atlantic Anomaly (SAA) in the Earth's magnetic field which perturbs CALIOP data between 30 • W and 80 • W (Noel et al., 2014). In practice, a 165 limited amount of data are usable in this region and only at night. After filtering, the data are further averaged at a resolution of 34 km for compactness. The other channels are processed in the same way.
Due to solar activity, CALIOP was not operating on 18 January and between 20 and 26 January. Hence, our CALIOP series start on 27 January. We use only night data in this work. The molecular backscatter is calculated following Hostetler et al. (2006). For each day, the backscatter ratio is zonally averaged over all available orbits of that day (14 to 15 for a nominal day).

170
The native vertical resolution in the 20-30 km range is 180 m.

A1.2 IMS
The RAL (Rutherford Appleton Laboratory) Infrared/Microwave Sounder (IMS) retrieval core scheme (Siddans, 2019) uses an optimal estimation spectral fitting procedure to retrieve atmospheric and surface parameters jointly from co-located measurements by IASI (Infrared Atmospheric Sounding Interferometer), AMSU (Advanced Microwave Sounding Unit) and MHS (Mi-175 crowave Humidity Sounder) on MetOp-B spacecraft, using RTTOV 12 (Radiative Transfer for TOVS) (Saunders et al., 2017) as the forward radiative transfer model. The use of RTTOV 12 enables the quantitative retrieval of volcanic-specific aerosols (sulphate aerosol) and trace gases (SO 2 ). The present paper uses IMS SO 2 and sulphate aerosols observations from its near-real time implementation. The IMS scheme retrieves the SO 2 in the sensitive region around 1100-1200 cm −1 , in ppbv assuming a uniform vertical mixing ratio. It retrieves sulphate-specific optical depth at 1170 cm −1 (i.e. the peak of the mid-infrared 180 extinction cross section (Sellitto and Legras, 2016)), assuming a Gaussian extinction coefficient profile shape peaking at 20 km altitude, with 2 km full-width half-maximum. The bulk of the spectroscopic information on SO 2 and sulphate aerosols, in the IMS scheme, thus comes from the Infrared Atmospheric Sounding Interferometer (IASI) (Clerbaux et al., 2009). We refer to the two retrieved product as IMS SO 2 and SA optical depth (SA OD in the figure titles) in this work. The data are provided daily on a regular grid with 0.25 • resolution in latitude and longitude with one image collecting the day swaths and another 185 collecting the night swaths.

A1.3 OMPS-LP
The Ozone Mapping and Profiler Suite Limb Profiler (OMPS-LP) onboard the Suomi-NPP satellite provides along track vertical profiles of aerosol extinction in several visible bands (Loughman et al., 2018;Taha and Loughman, 2020). We use version 2.1 and the 745 nm band as recommended by Taha et al. (2021). Swaths with non zero quality flag are discarded. Basically, this 190 filters data polluted by the SAA but filtered and non filtered results differ very little in our processing. The molecular extinction is calculated from the same formulas as the CALIOP molecular backscatter but for a change of wavelength. The extinction is averaged daily over all available orbits of that day and after a horizontal interpolation to a latitude grid of 1.1 • resolution that corresponds to the mean resolution of OMPS-LP in the considered range of latitudes.
OMPS-LP has a vertical resolution of 1.5 km which is lower than the vertical resolution of CALIOP. It is also sensitive to the 195 arch effect (Gorkavyi et al., 2021) for that tends to shift downward by several kilometers the apparent bottom of an extended aerosol layer.

A1.4 MLS
The Microwave Limb Sounder (MLS) onboard NASA's AURA satellite provides along track vertical profiles of water vapour mixing ratio (Lambert et al., 2015;Schwartz et al., 2020) as well as other trace gases, temperature and cloud ice. We use 200 the version 4 without accounting the quality flag as in Millán et al. (2022). The data are projected and zonally averaged daily onto a fixed latitude grid of 1.45 • resolution in the domain of interest. As they are provided on pressure levels with an approximate vertical resolution of 1.5 km, similar to that of OMPS-LP, they are interpolated to altitudes using the geopotential calculated daily on the ERA5 zonal mean. In order to get estimates of the altitudes by the method described in Appendix A3, the interpolation is made to a resolution of 100 m using a non-oscillating Akima interpolator. The native vertical resolution of 205 the product is about 1.5 km in the relevant range of altitudes.

A1.5 ALADIN
The Atmospheric Laser Doppler Instrument (ALADIN) onboard the Aeolus satellite is the first space-borne Doppler wind lidar. It is designed to measure wind along the line of sight from the Doppler shift of the 355 nm light emitted by the laser and scattered back by molecules (Rayleigh wind) or aerosols (Mie wind). Horizontal line-of-sight wind is retrieved neglecting the 210 7 vertical wind component. The anomaly wind is calculated by removing the background wind at same time and location from ERA5. As the line-of-sight is perpendicular to the heliosynchronous orbit, the measured component at low and mid-latitudes is essentially the zonal wind. The ceiling of Aeolus vertical bins can be adjusted and was increased to 30 km in the area of the HTHH plume (30 • S-0 • ) a few days after the eruption. We use the Mie product which is of better quality than the Rayleigh product inside the plume Zuo et al. (2022).

A1.6 RGB-Ash
We use a composite RGB product, denoted as RGB-Ash, that benefits from the sensitivity of the 8.5 µm band of the Advanced Himawari Imager (AHI) and Spanning Enhanced Visible and InfraRed Imager (SEVIRI) onboard the geostationary Himawari-8 and Meteosat-8 satellites. The product is based on the EUMETSAT Ash RGB recipe (https://navigator.eumetsat. This product allows to qualitatively distinguish thick ash plumes or ice clouds (brown), thin ice clouds (dark blue) and sulphur-containing plumes (green). Mixed ash/sulphur-containing volcanic species would appear in reddish and yellow shades.

A2 ERA5 reanalysis and meteorological data
We use the European Center for Medium Range Forecasts ERA5 reanalysis (Hersbach et al., 2020) at 1 • × 1 • resolution and all the model levels with 6-hourly sampling at 0, 6, 12 and 18 UTC. Geopotential, potential temperature and potential vorticity are calculated at full resolution for each time. All the fields are then averaged in longitude and over the four daily samples to provide a daily zonal average. At the stratospheric altitudes which are relevant to this study, the model levels are pure pressure 230 and therefore the averages are made over isobars.
The total all sky radiative heating rate is converted into diabatic vertical velocity from the profile of pressure, geopotential and temperature. The motion of the isentrops with respect to the geopotential in the zonal mean is used to define the adiabatic vertical velocity (see Appendix A3).
The ERA5 does not assimilate the anomalous water vapour or the aerosols in the stratosphere and therefore cannot account 235 for their direct radiative effect, either shortwave absorption or longwave absorption and emission. However, it assimilates the induced temperature perturbation if large enough to be detected and then react to damp it by longwave radiative relaxation with a time-scale of the order of 5 days at 25km (Lestrelin et al., 2021).
In the present case, the water vapour radiative cooling creates a negative temperature anomaly overlapping the plume Schoe-  The Lait potential vorticity (LPV) used in Fig.8 is defined from the Ertel potential vorticity (PV) as where θ is the potential temperature in K.

A3 Vertical motion from CALIOP and MLS
The observed vertical motion is obtained from CALIOP and MLS by applying a second-order Savitsky-Golay filter with a 245 31-day window to the daily mean vertical location of CALIOP scattering ratio and MLS water vapour, retaining data above 2 and 6 ppmv offsets, respectively. The offset are defined to isolate the aerosol and the water plumes from the background.
The 31-day window has been adjusted from several trials with 11, 21, 31 and 41 days as the value beyond which the resulting motion curve was rid of short time fluctuations and did not change any more in shape. This was reached with the 31-day window for MLS and the 21-day window for CALIOP but in the sake of consistency we use the 31-day window for both.

250
The diabatic and adiabatic background vertical velocities are calculated from ERA5 zonal means. The diabatic motion results from the total radiative heating rate DT Dt | RAD multiplied by θ T δz δθ where (T, θ, z) are temperature, potential temperature and geopotential altitude. The adiabatic motion, which is always a small correction, is estimated as -∂θ ∂t | p δz δθ + ∂z ∂t | p . The calculations are made by centered finite differences on the model grid which is in pure pressure in the considered altitude range.

255
Two correction methods are applied to the observed descent of aerosols from CALIOP in order to estimate the sedimentation velocity with respect to the air. In the first method, the ascent of the water vapour plume seen by MLS is used as an estimate of air motion and in the second method, the sum of diabatic and adiabatic ERA5 vertical velocities provides this estimate. Then this air motion is added to the aerosol descent rate to estimate the sedimentation velocity. The first method is applicable to the period during which the aerosol and water vapour distributions overlap and the water vapour cooling perturbs the heating 260 rate estimate of ERA5. The second method applies at a later stage when the aerosol and water vapour distributions are well separated and the radiative effect of water vapour has been defeated by dilution. As the boundary between these two regimes cannot be defined accurately we show the results of the two methods.
The scattering aerosol radius is then estimated using the Stokes fall speed formula for small particles (Seinfeld and Pandis, 2016).

A4 Mie calculations
The theoretical extinction-to-backscatter ratio for the plume has been calculated using the Python-based miepython Mie code, available at: https://miepython.readthedocs.io/en/latest/. The extinction and backscatter coefficients have been estimated at 750 and 532 nm, respectively, to simulate OMPS and CALIOP observations. Typical sulphate aerosols refractive indices have been considered, with the assumption of very weakly absorbing particles (based on the results of ). Log-normal 270 size distributions with varying standard deviation are simulated, to study how this ratio changes with radius. Author contributions. BL, PS and AP conceived the study and conducted the analyses. EC and RS elaborated and provided the IMS product.
JUG, SK and FP were involved in discussions of the results and their interpretation. The paper was written by BL, CD, PS and AP. All authors contributed to the final version.

285
Competing interests. The authors declare no competing interests.