Stratospheric aerosol layer perturbation caused by the 2019 Raikoke and Ulawun eruptions and climate impact

In June 2019 a stratospheric moderate eruption occurred at Raikoke (48◦N, 153◦E). Satellite observations show the injection of ash and SO2 into the lower stratosphere and an early entrainment of the plume into a cyclone. Following the Raikoke eruption stratospheric Aerosol Optical Depth (sAOD) values increased in the whole northern hemisphere and tropics and remained enhanced for more than one year, with peak values at 0.040 (shorter-wavelength visible, higher northern latitudes) to 0.025 (shorter-wavelength visible, average northern hemisphere). Discrepancies between observations and models 5 indicate that ash has played a role on evolution and sAOD values. Top of the atmosphere radiative forcings are estimated at values between -0.3 and -0.4 W/m (clear-sky), and of -0.1 to -0.2 W/m (all-sky), comparable to what was estimated for the Sarychev eruption in 2009. Almost simultaneously two significantly smaller stratospheric eruptions occurred at Ulawun (5◦S, 151◦E) in June and August. Aerosol enhancements from the Ulawun eruptions had mainly an impact on the tropics and southern hemisphere. The Ulawun plume circled the Earth within one month in the tropics. Peak shorter-wavelength sAOD 10 values at 0.01 are found in the tropics following the Ulawun eruptions, and a radiative forcing not exceeding -0.15 (clear-sky) and -0.05 (all-sky). Compared to the Canadian Fires (2017), Ambae eruption (2018), Ulawun (2019) and the Australian fires (2019/2020) highest sAOD values and RF are found for the Raikoke eruption.

radiation flux simulations in the spectral range from 300 to 3000 nm, with a 0.1 nm spectral resolution. The radiative transfer equation is parameterized and solved as follows: -the solar flux spectra used to drive the simulations are taken from Kurucz 150 (Kurucz, 2005). -vertical profiles of temperature, pressure, humidity and gas concentration come from the climatological standards of the Air Force Geophysics Laboratory (AFGL). Mid latitude standard profiles are used for simulations of the Raikoke plume, while tropical standard profiles are used for Ulawun. -The molecular absorption is parameterized with the LOWTRAN band model (H. Pierluissi and S. Peng, 1985) (as adopted from the SBDART code). -The radiative transfer equation is then solved with the SDISORT method (the pseudo-spherical approximation of the discrete ordinate method (DISORT)). The vol-155 canically perturbed simulations are carried out by adding average SAGE III/ISS profile observations of the volcanic aerosol extinction coefficient (details on the spatio-temporal identification of the volcanic perturbations are described in Sect. 4). As baseline, SAGE III/ISS aerosol extinction profiles are taken for background conditions, i.e. without volcanic aerosols (details on the spatio-temporal identification of the background are described in Sect. 4). For both setups (background and volcanically perturbed) we carry out multiple runs with varying solar zenith angles (SZA). Finally, the daily-average shortwave TOA ra-160 diative forcing is calculated by integrating the SZA-averaged upward diffuse irradiance for the background scenario over the whole shortwave spectral range. The shortwave surface radiative forcing is calculated with the SZA-averaged downward global irradiance with aerosols minus the background scenario, integrated over the whole spectral range.

WACCM model
Model simulations were performed using the global CESM1 (Community Earth System Model 1) using its Whole Atmosphere 165 Community Climate Model (WACCM) module linked to the CARMA (Community Aerosol and Radiation Model for Atmospheres) module, involving the sulfur cycle with a sectional aerosol scheme (English et al., 2011). Land, sea ice, and rivers were active modules, whereas oceans were prescribed. The spatial resolution was a longitude/latitude grid of 144 points by 96, respectively (i.e. approximately 2 • resolution), and over 88 levels of altitude ranging from the ground to approximately 150 km altitude with approximately 20 levels in the troposphere. Specified dynamics were used, with a nudging towards the 170 Modern-Era Retrospective analysis for Research and Applications 2 (MERRA2) meteorological data (Randles et al., 2017) at every time step (30 min) with a weight factor of 0.1 towards the analysis, for temperature and wind fields. Anthropogenic surface emissions were prescribed for SO 2 using the MACCity data set (e.g., Diehl et al., 2012). Carbonyl sulfide (OCS) was prescribed using data from (Kettle et al., 2002). The simulation presented in this study deals with a multi-annual model experiment starting on January 1 st 2013 using the CESM1 initial atmosphere state file at that date. The Raikoke and Ulawun 175 eruptions have been simulated by injecting a volcanic SO 2 mass burden into model grid boxes corresponding to the location of the volcanoes (Raikoke: 48 • N and 153 • E, Ulawun 5 • S and 151 • E), over 6 hours, spread evenly between a certain altitude range for each eruption (see Table 1 for a summary of the model setup) following the method of (Mills et al., 2016). The chosen SO 2 burden of 1.5 Tg is in fairly good agreement with Muser et al. (2020), who calculate 1.37 ±0.07 x 10 9 kg with TROPOMI and estimate 1-2 x 10 9 kg with HIMAWARI data. The model's 2.5 • longitude x 1.875 • latitude grid resolution means that 180 the volcanic plumes are initially too dilute in the model compared to reality. This is nevertheless a typical methodology used in the literature (e.g., Lurton et al., 2018). The timing and altitude injection of the SO 2 emissions is based on information  (2015)), and MLS/Aura (Krotkov et al., 2008). However, an OMPS aerosol extinction profile 185 shortly after the Raikoke eruption, supports the chosen altitude range of 9-16 km (see supplementary material, Figure A1).  (Carn, 2019a). Airplanes flying over the North Pacific had to be redirected (Crafford and Venzke, 2019).

Ulawun
The Ulawun volcano in Papua New Guinea (5 • S, 151 • E) was identified as one of the 16 'decade volcanoes' by the International Association of Volcanology and Chemistry of the Earth's Interior (IAVCEI) and is therefore known as one of the most 205 potentially destructive volcanoes on Earth (Cas, 2019). Two eruptions occurred during summer 2019, on June 26 th and August For the second and larger eruption, IASI/Metop-B data indicate SO 2 altitudes of around 14-17 km for August 3 rd and 4 th (Aeris, 2018). For the first eruption Sentinel5P/TROPOMI data suggest a SO 2 load of ∼0.14 Tg of the plume, while the second one was a bit larger and data suggest ∼0.2Tg of SO 2 (Carn, 2019b). With its tropical location, the eruptions at Ulawun have 210 the potential to have an impact on the lower stratosphere of both hemispheres within the BDC, once injected into the UTLS (Butchart, 2014). Ulawun remained in an active phase with e.g. observed ash plumes in October 2019 up to 3 km altitude (Bennis and Venzke, 2019). By February 2020 only water vapor plumes were observed and the Alert Level remained at Stage 1 (Sennert, 2020).

Injection and early dispersion of the Raikoke plume
Using a similar method as in Kloss et al. (2020), we attempt an estimation of the injection height using Himawari infrared brightness temperature information, at the moment of the main eruption, and comparing with coincident temperature profiles from ERA5 reanalyses. The brightness temperature of the plume core (not shown) exhibits a plateau at about 225 K within a few hours after the eruption. However, the exact injection altitude could not be identified due to the fact that the temperature profile for the second eruption. The initial evolution of the Raikoke plume is shown with the Himawari Dust RGB images starting from the 21/06 at 19:00 (Fig. 1). The Dust RGB product is used, instead of the Ash RGB product, because it is more sensitive for large satellite viewing angles, which is the case for the region of interest for Raikoke. This plume is initially composed of ash (reddish colors, in Fig. 1), with also some evidence of SO 2 (yellowish colors, in Fig. 1). The remaining greenish and pinkish colors indicate the presence of water clouds around the volcanic plume. Over June 22 nd , the plume disperses eastward  Fig. 1). In the following days, the ash plume is rapidly diluted or sediments, and cannot be further followed. The SO 2 plume instead persists and, from June 23 rd , stops mowing eastward to wrap upon itself and get trapped for several days within the cyclonic circulation of the Aleutian low which was unusually strong for this summer period. As a consequence, the confined plume remains compact and exhibits a number of dense patches and filaments that are well defined in the Himawari images, reaching locations as far as Alaska and central Russia, as visible from IASI D SO2 observations (Fig.   A3). CALIOP sections of these patches on June 25 th and 26 th (not shown) exhibit aerosol plumes up to 15.5 km. We found no confirmation of the rise to 22 km within a few days reported in the modelling study of (Muser et al., 2020). After June 25 th , the SO 2 plume gets more diluted and is possibly converted to sulfate aerosols. The presence of a compact SO 2 plume, after ash removal, is supported by the strong detection of SO 2 , i.e. D SO2 significantly larger than 1.0, that is obtained with the 245 high-spectral-resolution observations of IASI, starting from 23/06/2019, at about 9:00 am (morning overpass, Fig. 2b). The intensity of the D SO2 detection parameter decreases in the following days (Fig. 2c,d), as the plume dilutes and possibly a part of the SO 2 converts to sulfate aerosols.

The global dispersion of the Raikoke and Ulawun plumes with OMPS observations and WACCM simulations
After the first atmospheric processing following the injection in the UTLS, including the entrainment into the storm discussed   sphere and can result from a horizontal tropopause crossing of the aerosol plume towards the south (Fig. 4a). This hypothesis is confirmed by the model simulation in Fig. 4b and c, where only volcanic sources of stratospheric aerosols are considered.
origin tracer with CLaMS (Fig. 5). The simulation cannot be taken for quantitative estimations for the following reasons. First, the chosen initialization is given in a box shape, whereas the real injection does not appear in the shape of a box. Therefore, many trajectories in this simulation 365 do not necessarily correspond to an actual plume air parcel during injection. Second, in this simulation we use a passive tracer, with no chemical/microphysical processes being taken into account. Finally, the injected burden and related quantitative factors are not accounted for in the CLaMS simulations, as the Raikoke and Ulawun air mass tracers are equally represented. However, as CLaMS transport is driven by the newest reanalysis (ERA5) the simulation provides a reliable diagnostic of the air mass transport from the volcano region (initialization box).

370
Once initialized after the Raikoke eruption, the air mass tracer is transported towards the East, which is consistent with OMPS observation (see Fig. 3). By mid-July (roughly within 3 weeks after the eruption), the plume tracer has circled the Earth on latitudes mostly north of the Raikoke location. At the beginning of July the main bulk of the air mass tracer remains west of the Atlantic Ocean. Therefore, the sAOD enhancement above Europe observed by OMPS in Fig. 3b does not originate from Raikoke, but rather from forest fires in Alberta, Canada. The plume air mass transport is largely consistent with OMPS 375 observations, as by the end of July (Fig. 3) enhanced AOD values are apparent throughout all longitudes, mostly north of the Raikoke position. For the CLaMS simulation a clear signal of the tracer is visible around the area of the AMA from end-July until mid-September, which is also consistent with OMPS data (Fig. 3c-e). By mid-August a small percentage of the initialized Raikoke tracer has reached the tropics in the CLaMS simulations. Such a transport can also be seen from OMPS and WACCM data in Fig. 4a and b in July/August 2019 (with sAOD values below 0.01 for OMPS). As seen for OMPS data, the plume tracer Even though CLaMS simulations neither take any chemical/microphysical processes into account nor possible lifting due to   Figure 6 shows the vertical distribution of aerosol extinction values, and its evolution, around the location of the volcano. The initial injection phase after the Raikoke eruption is more evident for the WACCM simulation than for OMPS observations ( Fig. 6a and b). In the model, the aerosol plume rises from around 15-20 km altitude during the month following the eruption, (second eruption), directly after the respective eruptions ( Fig. 6c and d). A subsequent transport to ∼21 km in the area around the volcano is also shown in observations and reproduced in the model. One month after the eruption, the signal of the dispersed plume is at higher altitudes in the observations than in the model. This can potentially reflect an underestimation of the amount of SO 2 initially injected in the model. As seen in Fig. 4, OMPS reveals increased aerosol extinction values even 10 months after the second Ulawun eruption, while WACCM values seem almost back at background conditions within 5 months. The large 405 differences between OMPS observations and the WACCM simulation seen in the troposphere can be explained by clouds and other tropospheric sources of aerosols, which are not included in the model. We focus on the transport in the lower stratosphere, rather than the troposphere, therefore, those differences are of no interest in this study.

Vertical distribution
The panel series in Fig. 6e shows sampling, reaching high latitudes from OMPS gives confidence in the representation of the overall AOD evolution (Fig. 7a).
While we present 3-day averages for the OMPS data set, we calculate 30-day averages for SAGE III/ISS, to account for the much sparser sampling of SAGE III/ISS.
The timing and total value of sAOD enhancements for OMPS and SAGE III/ISS ( Fig. 7a and b) Fig 7a and b). Particular sAOD enhancements from the  Figure 7a and b. Furthermore, this is consistent with Figure 3, which also shows enhanced sAOD values above France in August 2019. For LOAC, only partial AOD (in terms of particles size) have been derived for LOAC in situ data, i.e. in the range from 0.2-0.7 µm, to avoid spurious aerosol extinction enhancements resulting from the presence of low-concentrated micrometer-sized particles, for instance coming from the balloon flight chain above the instrument or corresponding to the "background" meteoritic population (e.g., Renard et al., 2008;Murphy et al., 2014). As a result, the LOAC AOD values cannot be directly compared with OMPS. The in situ AODs reveal a significant enhancement over the 2019 summer-autumn period above France. Following the Raikoke eruption, the in situ data present an oscillating behaviour with some low values in late 2019 (especially the October measurement in Fig. 7c). This could reflect the sparse and very local sampling of in situ observations and could also be explained by a still inhomogeneous volcanic plume at this period. The slight 450 increase in the observed AOD in April 2019 can be related to remnants of the midlatitude signature of the Ambae eruption (Kloss et al., 2020) and could reflect that background aerosol conditions were not reached in the stratosphere for the period before the Raikoke eruption, which is consistent with OMPS and SAGE III/ISS observations in Fig. 7a and b.

455
The multispectral SAGEIII/ISS observations are used to further characterize the optical properties of the Raikoke and Ulawun plumes and to estimate their radiative forcing (RF). Despite their sparser spatiotemporal sampling, with respect to OMPS, the solar occultation geometry of SAGE III/ISS observations is associated with a better signal-to-noise ratio. Figure 8a September 2019 the stratosphere is expected to be still somewhat perturbed by late Raikoke plume, the selected period is chosen to be representative for both peak and declining volcanic perturbation (see Fig. 4a). To get a more detailed characterization of the plume and its impact, we subdivided the overall latitude range chosen for Raikoke into two sub-intervals: 40-55 • N and 55-70 • N. It is important to mention that latitudes higher than 70 • N are very sparsely sampled with the SAGEIII/ISS  Fig. 8c,d. Regional RF values as large as -2 to -3 W/m 2 are found for Raikoke, at both TOA and surface, in the 40-55 • N and 55-70 • N, respectively, for the assumption of very small (g=0.5) and very reflective (SSA=1.0) particles, which is linked to a significant cooling of the regional climate system and a very limited, if any, energy absorbed by the plume. The TOA RF at the highest northern latitudes (70-90 • N) is found to values as large as -5 W/m 2 but this estimation has to be taken with caution (as discussed above, it is based on an extrapo-510 lation). For smaller SSA, the TOA and surface RF start to deviate significantly (larger surface than TOA RF), thus indicating a significant absorption of radiative energy into the plume. This energy imbalance and the possible resulting radiative heating of the plume can be a possible reason for the observed lifting, shown in Fig. 6e; this hypothesis requires further investigation.
The assumption on the asymmetry parameter g dominates the uncertainty of the RF estimations (compare the error bars of Fig.   8c and d). It is important to mention that all these RF estimations are based on the assumption of clear-sky, so these are just a reference and have to be scaled down to take into account the impact of clouds in reducing the effective RF.
Based on the above mentioned regional clear-sky RF estimations in the shortwave (Table 2), the equinox-equivalent daily average shortwave global TOA radiative forcing of Raikoke and Ulawun plumes, based on their stratospheric aerosol layer perturbations, can be estimated. We calculated this as a latitude-weighted mean of the regional RF, extended over the whole globe, by considering a zero-impact outside the regions defined in this section. Because we know that the Raikoke plume had 520 an influence on the tropics (which is here considered as a 'zero impact region'), the calculated global clear-sky RF values are