Improved estimation of volcanic SO2 injections from satellite observations and Lagrangian transport simulations: the 2019 Raikoke eruption

Monitoring and modeling of volcanic plumes is important for understanding the impact of volcanic activity on climate and for practical concerns, such as aviation safety or public health. Here, we applied the Lagrangian transport model Massive-Parallel Trajectory Calculations (MPTRAC) to estimate the SO2 injections into the upper troposphere and lower stratosphere by the eruption of the Raikoke volcano (48.29◦N, 153.25◦E) in June 2019 and its subsequent long-range transport and dispersion. First, we used SO2 observations from the AIRS (Atmospheric Infrared Sounder) and TROPOMI (TROPOspheric 5 Monitoring Instrument) satellite instruments together with a backward trajectory approach to estimate the altitude-resolved SO2 injection time series. Second, we applied a scaling factor to the initial estimate of the SO2 mass and added an exponential decay to simulate the time evolution of the total SO2 mass. By comparing the estimated SO2 mass and the observed mass from TROPOMI, we show that the volcano injected 2.1±0.2 Tg SO2 and the e-folding lifetime of the SO2 was about 13 to 17 days. The reconstructed injection time series are consistent between the AIRS nighttime and the TROPOMI daytime measure10 ments. Further, we compared forward transport simulations that were initialized by AIRS and TROPOMI satellite observations with a constant SO2 injection rate. The results show that the modeled SO2 change, driven by chemical reactions, captures the SO2 mass variations observed by TROPOMI. In addition, the forward simulations reproduce the SO2 distributions in the first ∼10 days after the eruption. However, diffusion in the forward simulations is too strong to capture the internal structure of the SO2 clouds, which is further quantified in the simulation of the compact SO2 cloud from late July to early August. Our 15 study demonstrates the potential of using combined nadir satellite observations and Lagrangian transport simulations to further improve SO2 timeand height-resolved injection estimates of volcanic eruptions.

For the analysis of the Raikoke eruption we used the TROPOMI Level 2 offline (OFFL) V01.01.07 SO 2 data product for the 125 time period between 20 June 2019 and 16 August 2019. The TROPOMI SO 2 data product provides four total vertical columns of SO 2 in mol m −2 , one for the total atmospheric column between the surface and the top of the atmosphere and three columns assuming an SO 2 layer at 1, 7, and 15 km altitude in the retrieval. The details of the retrieval are given in Theys et al. (2017Theys et al. ( , 2021. In studies investigating volcanic plumes, it is common to use the vertical column densities retrieved for distinct plume heights (e. g., Theys et al., 2019). In this study, we used the total vertical SO 2 column of the 15 km retrieval, as we considered it 130 to provide the best approximation for the Raikoke eruption as in other studies (e. g., Muser et al., 2020;de Leeuw et al., 2021).
Compared with AIRS, the lower detection limit for TROPOMI is 0.3 DU (Theys et al., 2021) and data below this threshold are discarded in this study.

The MPTRAC model
Massive-Parallel Trajectory Calculations (MPTRAC) is a Lagrangian particle dispersion model for the analysis of atmospheric 135 transport processes in the troposphere and stratosphere . It calculates particle trajectories by solving the kinematic equation of motion using given wind fields from reanalysis or forecast meteorological data. The MPTRAC model currently uses the midpoint method to solve the equation of motion, which gives the optimized balance between accuracy and computational efficiency (Rößler et al., 2018). Besides vertical motion driven by the vertical velocity (i. e., kinematic trajectories), the MPTRAC model provides options to constrain the pressure of the air parcels to either constant pressure 140 (isobaric surface), constant density (isopycnic surface), potential temperature (isentropic surface), or pressure time series from Figure 2. Spatial distribution of SO2 vertical column density (DU) during the 24 hour period between 26 June 2019, 12:00 UTC and 27 June 2019, 12:00 UTC from AIRS nighttime (a), AIRS daytime (b), and TROPOMI daytime (c) observations. Note that AIRS observations of SO2 vertical column density less than 5 DU are not shown here as those data are not actually used in the analysis because they are affected by background noise. The AIRS nighttime observations have a ∼12 hour time shift compared with the AIRS daytime and TROPOMI observations. balloon measurements (Hoffmann et al., 2017). In addition, the model also includes turbulent diffusion and subgrid scale wind fluctuations to simulate the diffusion. The turbulent diffusion is described by fixed diffusivity coefficients. Following the FLEXPART model (Stohl et al., 2005), the MPTRAC model uses a constant horizontal diffusivity of 50 m 2 s −1 for the troposphere and a vertical diffusivity of 0.1 m 2 s −1 for the stratosphere as default values. Subgrid scale wind fluctuations are 145 simulated using the Langevin equation to add time-correlated stochastic perturbations to the trajectories. The subgrid scale wind standard deviations are downscaled from the grid scale standard deviations by using a default scaling factor of 40 % (Stohl et al., 2005). To investigate the effect of parameterizations of turbulent diffusion and subgrid scale wind fluctuations on the simulated SO 2 dispersion for the Raikoke case, we varied the diffusivity and the scaling factor of the subgrid variance separately. As the actual diffusivity can vary by several orders of magnitude (e. g., Ishikawa, 1995;Desiato et al., 1998;Legras 150 6 https://doi.org/10.5194/acp-2021-874 Preprint. Discussion started: 9 November 2021 c Author(s) 2021. CC BY 4.0 License. et al., 2005;Pisso et al., 2009), we tested the turbulent diffusion by varying the diffusivities from 10 −2 to 10 3 times the default values. For subgrid scale wind fluctuations, we varied the scaling factor from 0 to 100 %.
Additional modules are implemented to simulate convection, sedimentation, radioactive decay, hydroxyl chemistry, dry deposition, and wet deposition. In this study, we used the hydroxyl chemistry module to simulate the loss of SO 2 by its reaction with the hydroxyl radical (OH). The MPTRAC model also provides variable output methods. In this study, we implemented a new module for "sample output", which allows us to sample the model data at the exact time and location of the satellite overpasses/footprints. In addition to the trajectory and gridded outputs, the model also provides ways to directly evaluate the performance of the simulations, such as calculating the critical success index (CSI) (Wilks, 2011). Basically, the CSI is based on the counts of detection by the observation and simulation on a regular grid basis. If the vertical column density in a grid cell passed a user specified threshold, it will be counted as "yes", otherwise it will be counted as "no". The CSI is the ratio between 160 the count of hits and the total number of hits, false alarms, and misses. Along with CSI, the probability of detection (POD) and the false alarm rate (FAR) are also calculated.
In this study, the main meteorological data used to drive the MPTRAC simulations is taken from the ERA5 reanalysis.
The ERA5 is the ECMWF's (European Centre for Medium-Range Weather Forecasts) fifth generation reanalysis (Hersbach et al., 2020), which is meant to replace its predecessor ERA-Interim (Dee et al., 2011). ERA5 provides hourly outputs of a 165 comprehensive set of variables at 31 km horizontal resolution and 137 levels spanning from the surface up to 0.01 hPa. In this study, the ERA5 data are interpolated to 0.3 • × 0.3 • horizontal resolution. In comparison, the ERA-Interim data have a horizontal resolution of 80 km, 60 model levels, and output every 6 hours, i. e., at 00:00, 06:00, 12:00, and 18:00 UTC. The differences between ERA5 and ERA-Interim in driving Lagrangian transport simulations have been assessed by Hoffmann et al. (2019), finding that the choice of data has a considerable impact on the simulations, in particular due to better spatial and 170 temporal resolutions of the ERA5 data. We also considered both, ERA5 and ERA-Interim data in this study, with a major focus on results derived from the ERA5 reanalysis.

Estimation of volcanic SO 2 injections
To reconstruct the altitude-and time-resolved injection parameters of the Raikoke eruption, being represented by the altitude, time, and SO 2 mass of each air parcel over the volcano, we used a method based on backward trajectories released from the 175 columns of the AIRS and TROPOMI SO 2 measurements Wu et al., 2017Wu et al., , 2018. The analysis was done separately for AIRS and TROPOMI data, covering time periods from a few days up to weeks after the eruption, and the results were compared against each other. As both, AIRS and TROPOMI provide information on the horizontal location and time of the SO 2 observations, but lack information on the vertical distribution of the SO 2 , we released multiple air parcels between 0 and 25 km altitude at each individual satellite footprint with volcanic SO 2 detections. In contrast to our earlier work, 180 the vertical profile of the number of air parcels has been made to follow the mean kernel function of the satellite measurements ( Fig. 1) in order to take into account their different vertical sensitivity. The total number of air parcels at each location was linearly proportional to the total column density of the satellite observations. At the same time, a Gaussian scatter of the air parcels with 15 and 5 km full width at half maximum (FWHM) was introduced to represent the horizontal footprint size for AIRS and TROPOMI, respectively.

185
In total, 5 million air parcels were released to calculate backward trajectories. If a backward trajectory passed the Raikoke volcano within a search radius of 15 km, the location and time of the air parcel was saved to reconstruct the injection parameters.
We note that, based on our sensitivity tests, the results are not very sensitive to FWHM (i. e., between 1 and 50 km) and the search radius around the volcano (i. e., between 1 and 100 km). All backward trajectories that met the selection condition were re-sampled to a total number of 5 million particles and an initial total mass of 1.5 Tg was assigned to them. After the initial 190 relative SO 2 distribution has been estimated from the backward trajectories, we conducted forward simulations and applied a scaling factor to the SO 2 total mass for further calibration. To calibrate the total mass, we assumed that the SO 2 starts to decay exponentially with a fixed e-folding lifetime immediately after the injection and compared the change of the SO 2 total mass from the simulations with the change of the SO 2 total mass derived from TROPOMI observations.  Figure 3a shows the final reconstruction of the Raikoke SO 2 injections based on TROPOMI observations and Lagrangian transport simulations using ERA5 winds. The mass in the reconstruction has been turned to achieve a total injection of 2.1 Tg.
The altitude-and time-resolved injection, and the integrated vertical profile are also shown in Fig. 3. A major SO 2 emission was 200 reconstructed during the first two days of the time series, i. e., between 21-22 June 2019. After this major eruption, significantly smaller amounts of SO 2 were continuously injected by the volcano until the end of June with a prominent second and third plume during 24-25 June and 27-28 June, respectively. The first plume crossed the tropopause and injected SO 2 between 5 and 15 km of altitude, with ∼45 percent of the SO 2 mass reaching the stratosphere ( Fig. 3a and b). The second and the third plumes mainly injected material into the troposphere. As the Raikoke eruption is dominated by the first plume, the overall injected 205 SO 2 (Fig 3c) distributes around the tropopause with peak injections at an altitude of 11 km.
As an intercomparison as well as a validation, the vertical profiles integrated over the entire eruption period (21 to 30 June 2019) of our different injection estimations and the profile derived by the VolRes team (de Leeuw et al., 2021) are shown together in Fig. 4. The profile derived by the VolRes team is mainly based on IASI observations during the first two days of the Raikoke eruption (de Leeuw et al., 2021). Compared with the VolRes profile, the altitude of peak injections is about 210 1 km above the VolRes profile, no matter which satellite data, i. e., AIRS nighttime or TROPOMI daytime, and reanalysis data, i. e., ERA5 or ERA-Interim, have been used. Our reconstructions also show enhanced injections between 12 and 14 km, being consistent with de Leeuw et al. (2021) that the injections reached higher into the stratosphere than indicated by the VolRes estimation. When excluding the second and third plume, the vertical profile for the first major eruption (figure not shown) is similar with the overall injection profile ( Fig. 4) with slightly reduced emission rate in the stratosphere and larger reduction in the troposphere part. The VolRes profile also indicates a small peak at low altitudes around 2 km. However, this part is not found in our reconstruction. The most likely reason is that both, AIRS and TROPOMI, have a limited sensitivity in the lower troposphere, i. e., below 5 km (

Calibration of the total mass of the SO 2 injections
For the initial reconstruction, we estimated the SO 2 injection rates under the assumption of a SO 2 total mass of 1.5 Tg, as 220 found by the VolRes team. Figure 5 shows the time series of the vertically integrated SO 2 injections. The comparison of the reconstructions based on different satellite observations ( Figure 5) show that the results derived from AIRS nighttime and TROPOMI observations agree well. The reconstruction derived from the AIRS daytime observations shows weaker injections in the first plume, but the second and third plume are stronger compared to the reconstructions based on AIRS nighttime and TROPOMI measurements. As pointed out in Sect. 2.1, we will focus our analyses on the AIRS nighttime and TROPOMI 225 results in the following parts.
To estimate the final total injected SO 2 mass, we derived the daily SO 2 mass from the TROPOMI observations (Fig. 6). The TROPOMI observations shows that a total SO 2 mass of ∼1.4 Tg peaked during 24-26 June, while the cumulative SO 2 injection from the initial reconstruction at 26 June is only 1.2 Tg. When the cumulative SO 2 injection in the initial reconstruction reached 1.5 Tg, the total SO 2 mass from TROPOMI decreased to 1.2 Tg due to the removal of SO 2 . To better represent the evolution of the total SO 2 mass in the simulations, we scaled our initial mass reconstruction to the TROPOMI data and applied an exponential decay to account for the removal of SO 2 (Fig. 6). In this experiment, we found that a total injection of 1.9 to 2.3 Tg SO 2 and an e-folding lifetime of 13 to 17 days best represents the temporal evolution of total SO 2 mass in the atmosphere.
Therefore, we re-scaled the initial reconstruction to a total mass of 2.1 Tg. Note that although the e-folding lifetime of 13-17 days well represents the overall removal of SO 2 injections for the Raikoke case, SO 2 removal rates in general are very sensitive 235 to the altitude of the SO 2 injections and the atmospheric background conditions. Figure 6. Temporal change of total SO2 mass from TROPOMI measurements (blue), and calculated total mass with an e-folding lifetime of 15 days (orange) and a mass scaling factor of 2.1/1.5. Orange shadings show the combination of the scaling factor ranging between 1.9/1.5 and 2.3/1.5 and e-folding lifetime ranging between 13 and 17 days. The black curve shows the accumulated SO2 injection with a total injection of 1.5 Tg (the initial reconstruction).

Sensitivities of estimated SO 2 injections on data and model parameters
Investigating the sensitivity of the reconstructed injection time series on the underlying input data, we ran the reconstruction using AIRS SO 2 measurements together with ERA5 winds as well as TROPOMI measurements together with ERA-Interim winds. In comparison to TROPOMI and ERA5, the overall patterns of injection estimations based on AIRS nighttime obser- weaker (stronger) injections at the beginning (late) stage of the first plume, respectively. Differences do exist during the second and the third plumes, but they are relatively small. 245 We conducted more than two hundred simulations to test the sensitivity on the injection reconstructions. Among the tested parameters, we found that the coverage of the satellite observations has the largest impact on the injection reconstruction. More specifically, it matters how many days of satellite observations are used for the reconstruction and how close to the location of the volcano satellite observations are discarded. Here, we only describe the sensitivity tests on the temporal and spatial coverage of the satellite observations, whereas the sensitivity tests on other parameters are not shown.   Figure 8 shows the SO 2 mass change in forward simulations initialized by using TROPOMI observations covering different numbers of days since the beginning of the eruption. The total SO 2 mass in all the simulations was assigned to 2.1 Tg. As shown in Fig. 8, when using just a few days of observations, the simulation produces a too strong peak at the beginning of the volcanic eruption. Increasing the time period of the satellite data, a gradual decrease of the first peak and redistribution of SO 2 to a later stage of the eruption is observed. Therefore, using short term observations will lead to a more pronounced first plume, 255 and on the contrary, longer term observations will increase the amplitude of the second and third plume. The sensitivity test shows that using 12 days of observations gives an optimal representation of the SO 2 mass. As the backward trajectory method heavily relies on the quality of the trajectories, satellite observations too close to the volcano can not provide enough information to separate between different altitudes. Therefore, we defined a circle with a certain distance to the Raikoke volcano, were the satellite observations falling inside the circle are discarded. When using 260 TROPOMI observations and the distance being set to a very small value, such as a few kilometers, most of the reconstructed SO 2 is injected at the beginning of the eruption. When the distance is being set to several hundred kilometers, similar to increasing the temporal duration of the trajectories, the injection of SO 2 at the beginning of the eruption is weakened and more SO 2 is injected within a few to several days after the beginning of the eruption. When using the AIRS measurements, this effect becomes less pronounced. Overall, the AIRS and TROPOMI reconstructions agreed better with each other when setting 265 a larger distance. In the final reconstruction, we used a distance of 750 km, which gave the most consistent results. Unless noted differently, all forward simulations that were initialized by a constant injection rate in the following sections have the same setup as described here. In the following subsections, we will present and compare the forward simulations of SO 2 in terms of total mass burden and spatial distributions. In addition, we also performed different forward simulations driven by ERA5 and ERA-Interim data. However, the overall patterns of simulated SO 2 were generally similar between ERA5 and 275 ERA-Interim. Therefore, if not specified otherwise, the forward simulations driven by the ERA5 data are shown.
In the most recent version of the MPTRAC model, a hydroxyl (OH) chemistry module has been implemented to simulate the removal of SO 2 due to chemical reaction with OH. This module enables the direct comparison of total SO 2 mass change in model simulations with the simple exponential decay experiments and the satellite observations by TROPOMI (Fig. 9).
In the forward simulations, we have used different injection parameters with total injected SO 2 mass ranging from 1.9 to 280 2.3 Tg. As shown in Fig. 9, the SO 2 mass in the forward simulation initialized with a 2.1 Tg total injection, either using injection parameters derived from TROPOMI daytime (Fig. 9a) or AIRS nighttime (Fig. 9b) observations, agrees well with the exponential decay experiment of a 14-day e-folding lifetime. In addition, all the experiments are consistent with the total SO 2 mass derived from TROPOMI observations. Figure 9 also shows the SO 2 mass change in the forward simulation initialized by a constant injection time series. In contrast to the forward simulation initialized by our reconstructed injection time series, the 285 simulation initialized with a constant injection rate produced an SO 2 mass peak, which is comparable with the maximum SO 2 mass in TROPOMI observations, at the beginning of the eruption. Then, the gradual removal of SO 2 leads to lower mass in the model simulation than observed by TROPOMI (Fig. 9). From these comparisons, we conclude that the June 2019 Raikoke eruption produced a total injection of 2.1 Tg SO 2 , which has an overall e-folding lifetime of 14 days in the UTLS region during the first two weeks after the eruption.

290
The comparison of the temporal changes of the SO 2 total mass among the different forward simulation settings and TROPOMI observations ( Fig. 9) suggests that our estimation of 2.1 Tg SO 2 injection is reasonable. The initial estimation of 1.5 Tg mainly reflects SO 2 injections of the major eruption during the first two days. Consistent with this estimate, the total mass in our estimation for the first plume is about 1.5 Tg (Fig. 3b). However, additional injections after the first plume are required to reproduce the observed SO 2 mass in the model simulations (Fig. 9).

295
Although an e-folding lifetime of 14 days well captures the overall mass reduction of injected SO 2 in the atmosphere, the real removal rates of SO 2 at different altitudes are different. Figure 10 shows the remaining mass of SO 2 injected to 1 km thick layers during the first 12 hours of the eruption (21 June 2019, 18:00 UTC to 22 June 2019, 06:00 UTC). In general, the removal rate decreases with altitude, mostly because of lower OH concentrations in the lower stratosphere compared to the troposphere.
In the troposphere, the SO 2 mass is reduced to less than 50 % within several days to a week, while the stratospheric injections

Simulation of SO 2 transport during the first ∼10 days
In general, the forward simulations initialized by both the TROPOMI daytime and AIRS nighttime observations well reproduce the spatial distribution of SO 2 during the first ∼10 days of simulation time, especially in terms of spatial location and extent.
We note that the forward simulations initialized by the TROPOMI and AIRS nighttime observations are highly consistent with each other, as the injection parameters estimated from these two datasets do not differ fundamentally (Fig. 7). Therefore, 310 the results from the forward simulation initialized by the AIRS nighttime observations are not shown here. To illustrate the performance after major SO 2 injections of Raikoke, we selected three satellite overpasses on 23 June, 25 June, and 28 June to show the SO 2 distribution in observations and model simulations (Fig. 11).
The TROPOMI observations and the MPTRAC simulations show that the SO 2 injections separated into two major clouds ( Fig. 11). We added mean trajectories for injections between 7-8 km and 11-15 km in Fig. 11 to indicate the major movements 315 of the two clouds. Both of the clouds are moving cyclonically. A smaller SO 2 cloud, which is represented by the mean trajectory for injections between 7-8 km, moves faster and the SO 2 column density also decreased very fast to less than 10 DU on 28 June.
The SO 2 column density in the other major cloud, which mainly reflects the stratospheric injections, decreased much slower compared with the SO 2 cloud that reflects tropospheric injections. This observation by TROPOMI is consistent with the faster removal of SO 2 in the troposphere (Fig. 10).

320
After the first injection, TROPOMI observations show that the SO 2 clouds are located to the east of the Raikoke volcano, but split into two branches with one branch being in the north and the other branch being in the south (Fig. 11a). The forward simulations initialized by the TROPOMI observation and a constant injection rate reproduce the main northern branch, which locates just to the east of the Raikoke (Fig. 11b and 11c). However, both simulations only reproduce a part of the southern branch and the part reproduced in the simulation initialized by a constant injection rate is too strong. Comparing with the major 325 northern branch SO 2 cloud, note that the southern branch is very weak with SO 2 column densities mostly less than 10 DU. A sensitivity test suggested that the southern branch is mainly associated with transport of SO 2 in the lower troposphere (between 0 and 5 km), which was not represented in both initializations.
After the second plume by 25 June (Fig. 11d -f), most of the SO 2 injections were moved to the northwest direction over the Asian continent and to the northeast direction over the northwest Pacific ocean (Fig. 11d). In addition to these two major SO 2 330 clouds, there is a weaker SO 2 cloud to the east of the Raikoke volcano, which is probably related to the injections between 23-25 June (Fig. 3). The forward simulation initialized by TROPOMI observations reproduced the general pattern of the three clusters. However, the forward simulation initialized by a constant injection rate only reproduces the two major SO 2 clouds in the northwest and the northeast directions.
After the third plume by 28 June (Fig. 11g), the SO 2 cloud that was over the Asian continent in the northwest direction 335 now moved back to the east of Raikoke along a cyclonic circulation. In contrast, the SO 2 cloud that was over the northwest Pacific ocean showed a slower movement and is now located over the Asian continent (Fig. 11g). Both forward simulations reproduced the two major SO 2 clouds. The forward simulation initialized by the TROPOMI observations reproduced a stronger SO 2 cloud to the east of the Raikoke volcano (Fig. 11h) due to injections during the third plume. Similarly, the simulated SO 2 cloud over the volcano after the second plume is also stronger than in the observations (Fig. 11e). This result indicates that our reconstructed injection parameters potentially overestimate the second and the third plume. However, totally removing the second and/or the third plume would severely reduce the ability to correctly simulate the total SO 2 mass burden as shown in Fig. 9.
Although the forward simulations can reproduce the observed SO 2 distributions with relatively high performance during the first week (Fig. 11), the simulation starts gradually losing the ability to capture the structures of the SO 2 cloud thereafter.
345 Fig. 12 shows the observed and simulated spatial distribution of the SO 2 at the beginning of 1 July. Overall, the SO 2 distributes like a strip pattern with major peaks over the Sea of Okhotsk and the west coast of the Bering Sea (Fig. 12a). The model simulation captured this overall strip like pattern and even some fine details over northern high latitudes (Fig. 12b). However, the simulated SO 2 distribution does not correctly reproduce the peaks over the Sea of Okhotsk and the west coast of the Bering Sea. Several days later, the observed SO 2 over the west coast of the Bering Sea gradually spreads out and its vertical column 350 density gradually attenuates (not shown). In contrast, the two SO 2 peaks over the Sea of Okhotsk retained their compact structures and relatively high vertical column density. From mid to late July, the two peaks over the Sea of Okhotsk eventually comprise the main parts of the remaining SO 2 of the Raikoke eruption and developed into two compact SO 2 clouds. As the simulated SO 2 did not reproduce the two peaks over the Sea of Okhotsk, however, the forward simulation lost its ability to capture the observed SO 2 distribution during the first week of July.

Assessment of forward simulations by means of the Critical Success Index
We performed analyses of the Critical Success Index (CSI) to evaluate the forward simulations at five different detection thresholds ranging from 0.3 to 50.0 DU (Figs. 13 and 14). Figure 13 shows the CSI, POD, and FAR time series for the forward simulation initialized with the TROPOMI observations. The reference for calculating the CSI, POD, and FAR are also the TROPOMI observations to keep consistency of data between observations and simulations. The smallest threshold considered 360 here represents the lower detection limit of the TROPOMI observations, which means that a threshold of 0.3 DU includes all available TROPOMI observations of the Raikoke event in the Northern Hemisphere. For the detection threshold of 0.3 DU the forward simulation produced a very high POD, being around 80% or larger. Such a high POD suggests that the overall spatial extent of the SO 2 distribution is well reproduced by the forward simulation. However, the FAR shows a significantly increasing trend towards the end of the simulation, which suggests that the forward simulation transported some SO 2 outside 365 of the observed SO 2 clouds. Due to the increasing trend of the FAR, the CSI values peak at the beginning of the simulation with a maximum value of 77% and gradually decrease to ∼20% after 10 days. Compared with CSI analyses in previous studies on other volcanic eruption events Hoffmann et al., 2016), in which the CSI decreased to below 10% after 10 days of forward simulations, our study for the Raikoke eruption shows improved performance of the model, meteorological input data, and satellite observations.

370
When increasing the detection threshold, however, the model performance decreases. For instance, the POD shows a clear decreasing trend when the detection threshold is increased from 0.3 to 50 DU (Fig. 13a). The differences of the FAR between the different detection thresholds are smaller (Fig. 13b). This result suggests that although the forward simulation well reproduced the overall spatial extent of the SO 2 clouds, it has less ability to reproduce the internal structure and the location of the maxima of the SO 2 clouds. For example, the POD at all SO 2 thresholds except for the 0.3 DU level decreased notably during the first 375 week of July (Fig. 13a) agreeing with our earlier findings that the forward simulation has lost the ability to capture structures of the SO 2 clouds at this time.
For comparison, we also assessed the performance of the forward simulation initialized by the AIRS nighttime observations with reference to the AIRS nighttime observations (Fig. 14). As the AIRS observations have a higher background level (about 5 DU), the lowest detection threshold of 0.3 DU to assess the CSI is not very meaningful in this case. The trend and magnitude 380 of the POD from AIRS observations are very similar to the simulation with TROPOMI observations, but the FAR has lower values of 20 -40 %. During the first 10 days of the simulation, the CSI values for a detection threshold of 5.0 DU are between 40-80%, which is about 1.5 times higher than for the simulation with TROPOMI observations.
To make a comparison between forward simulations with different initializations, we used the TROPOMI observations as a common reference and the detection threshold was set to 5.0 DU. Figure 15 shows the POD, FAR, and CSI time series of 385 forward simulations initialized with TROPOMI observations, AIRS nighttime observations, and a constant injection rate. In comparison to the TROPOMI observations, the POD, FAR, and CSI are very consistent between forward simulations initialized with TROPOMI and AIRS nighttime observations. Compared with the simulation with a constant injection rate, the POD is constantly higher for simulations initialized by observations. However, the simulations initialized with the satellite observations  suffer from higher FAR between 29 June and 4 July. In summary, the overall trends of POD, FAR, and CSI are generally 390 similar between forward simulations with different initialization settings. This result indicates that the quality of the forward simulations is less affected by the injection parameters as estimated by the backward trajectory method, probably due to the fact that the major SO 2 injection occurred during a small time window at the beginning of the eruption in this case.  all the forward simulations reproduce the trajectory of the observed AC SO 2 cloud (Fig. 18a). The MAD of the AC SO 2 in the TROPOMI observations is ∼1 degree in the longitude dimension and ∼0.6 degree in the latitude dimension. A common problem in the simulations is that the simulated dispersion is much stronger than the observed dispersion, even when the dispersion parameters were set to zero and just the initial horizontal spread of the cloud was taken into account. In the longitude dimension, the MAD in simulations is roughly one order of magnitude larger than the MAD in observations after 5 days of We also tested the role of diffusivity ranging from 10 −2 to 10 3 times of the default diffusivity values (see Sect. 2.3 for reference). Similar to the results in Fig. 18, the simulated dispersion is too strong, but the simulation is more sensitive to diffusivity than subgrid scale variance (Fig. 19). When the diffusivity is on the order of ≥ 10 times of the default diffusivity values, the simulation also does not capture the trajectory of the AC SO 2 cloud anymore.
We performed forward trajectory simulations with vertical motion driven by constant potential temperature to avoid poten-450 tial vertical velocity fluctuations due to data assimilation. The degree of dispersion shows less sensitivity to the dispersion parameters and simulated MAD for longitude and latitude are very close to the results from the default setting (Fig. 18). As the 29 https://doi.org/10.5194/acp-2021-874 Preprint. Discussion started: 9 November 2021 c Author(s) 2021. CC BY 4.0 License.
vertical position is adjusted at every time step to retain a constant potential temperature, the dispersion in the vertical direction is suppressed in the isentropic mode. Further, as the AC SO 2 cloud is located in the stratosphere, the turbulent diffusion is driven by vertical diffusivity, which is also switched off in the isentropic mode. Therefore, dispersion in the isentropic trajec-455 tory simulation is less sensitive to the choice of dispersion parameters and is weaker than the dispersion in kinematic mode.
Taken together, our tests on the dispersion parameters suggest that the reason why the simulated SO 2 cloud generally is too dispersive is not only due to too strong diffusion in the model, but the stretching effect associated with horizontal wind shear also seems to play an important role.