Online treatment of eruption dynamics improves the volcanic ash and SO2 dispersion forecast: case of the Raikoke 2019 eruption

In June 2019, the Raikoke volcano, Kuril Islands, emitted 0.4-1.8 x 10 kg of very fine ash and 1-2 x 10 kg of SO2 up to 14 km into the atmosphere. The eruption was characterized by several phases or puffs of different duration and eruption heights. Resolving such complex eruption dynamics is required for precise volcanic plume dispersion forecasts. To address this issue, we coupled the atmospheric model system ICON-ART (ICOsahedral Nonhydrostatic – Aerosols and Reactive Trace gases) with the 1-D plume model FPlume to calculate the eruption source parameters (ESPs) online. The main inputs are the 5 plume heights for the different eruption phases that are geometrically derived from satellite data. An empirical relationship is used to derive the amount of very fine ash (particles <32μm), which is relevant for long range transport in the atmosphere. On the first day after the onset of the eruption, the modeled ash loading agrees very well with the ash loading estimated from AHI (Advanced Himawari Imager) observations due to the resolution of the eruption phases and the online treatment of the ESPs. In later hours, aerosol dynamical processes (nucleation, condensation, coagulation) explain the loss of ash in the atmosphere 10 in agreement with the observations. However, a direct comparison is partly hampered by water and ice clouds overlapping the ash cloud in the observations. We compared 6-hourly means of model and AHI data with respect to the structure, amplitude, and location (SAL-method) to further validate the simulated dispersion of SO2 and ash. In the beginning, the structure and amplitude values differed largely because the dense ash cloud leads to an underestimation of the SO2 amount in the satellite data. On the second and third day, the SAL values are close to zero for all parameters indicating a very good agreement of model 15 and observations. Furthermore, we found a separation of the ash and SO2 plume after one day due to particle sedimentation, chemistry, and aerosol-radiation interaction. The results confirm that coupling the atmospheric model system and plume model enables detailed treatment of the plume dynamics (phases and ESPs) and leads to significant improvement of the ash and SO2 dispersion forecast. This approach can benefit the operational forecast of ash and SO2 especially in case of complex and non-continuous volcanic eruptions like the 20

with respect to mass loading and structure, amplitude, and location of the plume. In addition, we discuss the separation of the ash and SO 2 plumes due to aerosol-radiation interaction. Finally, Sect. 4 concludes the paper.
2 Methods 2.1 ICON-ART modeling system 95 In this study, we performed simulations with the global weather and climate model ICON (ICOsahedral Nonhydrostatic model) together with the module for Aerosol and Reactive Trace gases (ART). ICON solves the full three-dimensional non-hydrostatic and compressible Navier-Stokes equations on an icosahedral grid and allows seamless predictions from local to global scales ; Heinze et al. (2017); Giorgetta et al. (2018)).
ART, being part of ICON, supplements the model by including emissions, transport, gas phase chemistry, and aerosol dy-100 namics in the troposphere and stratosphere (Rieger et al., 2015;Schröter et al., 2018;Weimer et al., 2017). Muser et al. (2020) reported and demonstrated latest improvements in ICON-ART with respect to the AEROsol DYNamics module (AERODYN), which is also used in the present paper. In AERODYN, aerosols are organized in 7 log normal distributions considering Aitken (as soluble), accumulation (as soluble, insoluble and mixed), coarse (as insoluble and mixed) and a giant mode (as insoluble).
For each mode, the prognostic equations for number density and mass concentration are solved keeping the standard deviations 105 of the modes constant. For the Aitken mode, nucleation, condensation, and coagulation are considered, while accumulation mode and coarse mode are affected by condensation and coagulation only. Shifting of particles into another mode occurs either when a threshold diameter is exceeded (shift into larger mode) or when a mass threshold of soluble coating on insoluble particles is exceeded (shift from insoluble to mixed mode) (Muser et al., 2020). In AERODYN, water and sulfate (also ammonium and nitrate) can condense on ash particles and therefore change the physical properties of ash, e.g., size, density, optical 110 properties. Changes of particle optical properties can further feedback on the radiation and atmospheric state. However, effects of aerosol dynamics on atmospheric humidity and clouds is not considered yet.

115
For a better estimation of the ESPs, we coupled ICON-ART online with the one-dimensional volcanic plume rise model FPlume (Folch et al., 2016;Macedonio et al., 2016). FPlume solves the equations of the buoyant plume theory (Morton et al., 1956) along the vertical plume axis. It includes processes like ambient air entrainment, plume bending due to wind, particle wet aggregation, energy supply due to water phase changes, particle fallout, and re-entrainment of particles (Folch et al., 2016). Figure 1 summarizes the procedures performed at every time step in which FPlume is active: First, vertical profiles for 120 wind, temperature, pressure, and humidity simulated with ICON serve as meteorological inputs for FPlume. In the second step, FPlume calculates the plume properties, i.e, total MER in case of a given plume height (as here) or plume height in case of a given MER. Thirdly, the fraction of very fine ash is determined based on plume height and total MER by using the relationship of Gouhier et al. (2019).
In the last step, ash is emitted into ICON-ART by multiplying the MER of very fine ash with the vertical profile derived 125 from the normalized Suzuki distribution, which is the same as the one used by (Marti et al., 2017): Here, S(z) describes the vertical emission profile, E is the emission rate of very fine ash, H p is the plume top height, and z refers to the height in the plume. We completely disregarded the mass of particles larger than 32 µm as this fraction has been shown to be irrelevant for long-range transport Rose and Durant (2009).
More details on the initialization of the ash particles is given in the next section. Besides ash, we also emitted SO 2 . Different from the ash emission, we prescribed the MER of SO 2 based on satellite estimates, but we released it into ICON-ART with the same profile and phases as the ash. This simplification was necessary as no further information on temporally SO 2 emission is available. Yet, during volcanic eruptions in general it is possible that the ash and SO 2 are emitted at different phases of the eruption.  Raikoke emits primary basaltic lava, and, therefore, we assumed the following exit conditions for FPlume: The exit temperature of 1273 K and exit water mass fraction of 3% were the same for all eruption phases (Mastin, 2007); the exit velocity was linearly changed depending on the plume top height relative to 14 km, where the exit velocity was set to 120 m/s. Our resulting MERs are insensitive to the input values in the range of 10%.
The equation for the fine ash fraction by Gouhier et al. (2019) depends on whether the SiO 2 content is high or low and 150 whether the conduit was opened or closed. As no information on the conduit has been available so far, we averaged the fine ash fraction for low SiO 2 -closed conduit and low SiO 2 -opened conduit.

Geometric plume heights
In addition to temperature and exit velocity, FPlume requires the plume height as input to calculate the MER. The height above the ellipsoid of the individual puffs was estimated by a recently developed geometric technique (Horváth et al., 2021a), which side view heights were in good agreement with independent geometric estimates derived from plume shadows and GOES-17-Himawari-8 stereo observations (Horváth et al., 2021b).

160
The Raikoke plume heights cannot be unambiguously determined by the traditional infrared brightness temperature method.
For most puffs, the minimum 11 µm brightness temperature (BT 11 ) falls within the narrow temperature range of the quasiisothermal layer above the tropopause and leads to multiple height solutions within a wide altitude range of 10-24 km. At certain times (e.g., 23:50 UTC on 21 June or 01:20 UTC on 22 June), the massive eruption plume is undercooled, even precluding the application of the temperature method. For the smaller plumes produced by less energetic puffs, on the other 165 hand, the BT 11 has a warm bias due to contributions from the warmer lower-level marine Sc layer around the volcano, resulting in underestimated heights. A detailed analysis of the Raikoke plumes, including a comparison of the various height estimates, is given in Horváth et al. (2021b). The uncertainty of the plume heights lays within a range of ±500 m.

Model configuration
We performed global simulations with ICON-ART using a horizontal grid size of roughly 13.2 km (R3B07 grid) and 90 vertical 2) with a varying fine ash fraction and aerosol-radiation interaction activated in ICON-ART (FPlume-rad); 2) the second experiment calculates the ESPs same as above, but neglects the interaction of aerosols and radiation (FPlume-norad).
The comparison of FPlume-rad and FPlume-norad allows to quantify the lifting of the volcanic plume due to radiation; 3) the third experiment derives ESPs with the empirical relationship by Mastin et al. (2009), and it emits volcanic compounds with a 180 prescribed fine ash fraction from the reference case (mean value for each phase) along a Suzuki profile (i.e., Eq. 1). It further assumes aerosol-radiation interaction (Mastin-rad).  Figure 2 shows the plume height and MER of very fine ash that is released into ICON-ART by FPlume (red dots) and the Mastin relationship (blue dots). The fraction of very fine ash relative to the total MER predicted by FPlume is in the order of 1.5 -3 % (not shown). In most phases, the MER calculated with FPlume is lower than the MER calculated with the Mastin equation, and the difference tends to be higher for larger plume heights. Since the exit parameters are fixed during each phase in the reference case, variation of the MER derived by FPlume must be due to changes in the atmospheric conditions. As the 190 relationship by Mastin et al. (2009) neglects atmospheric conditions and the fine ash fraction is fixed within one phase, the MERs of the very fine ash are constant within each phase. According to these MER values, the total mass of very fine ash emitted in the model for all puffs together is about 1.37 · 10 9 kg using FPlume and 1.75 · 10 9 kg using Mastin-derived MERs.
The total amount of very fine ash is evenly distributed as insoluble tracer over the accumulation, coarse, and giant modes. Muser et al. (2020) explained the detailed size distribution parameters of these modes. Following previous studies, we emitted a total of 1.5 · 10 9 kg SO 2 (Muser et al., 2020;Kloss et al., 2021;De Leeuw et al., 2020). However, the SO 2 release is linearly adjusted to the eruption heights and length of each phase (Table 1) as follows:

Himawari-8 Ash and SO 2
To validate the model results, we used column SO 2 and ash mass loadings estimated from the 16-band visible and infrared 205 Advanced Himawari Imager (AHI) onboard the Himawari-8 geostationary satellite at every full hour. Himawari-8 is operated by the Japanese Space Agency (JAXA) and the Japanese Meteorological Agency (JMA). A detailed description of the data product and methods used here is already given in Muser et al. (2020) and references therein. In short, SO 2 is retrieved by the AHI band centered near 7.3 µm, where the absorption of SO 2 is high. A further retrieval scheme, as described in Prata et al. (2004), was applied to minimize the interference with vapor and clouds. For volcanic ash retrievals, the AHI bands near 210 11.2 µm and 12.4 µm are considered. The ash retrievals were corrected by a mask that accounts for the effect of pixel containing meteorological clouds but were detected as completely cloud covered. Hereby, only pixels inside a 0.1 g m −2 contour line are considered and a 9x9 median filter smooths out 'spikes'.

SAL Method
The SAL method is an object-based quality measure originally developed to verify precipitation forecasts (Wernli et al., 2008, 215 2009). However, it has also been successfully applied for transport forecasts of volcanic compounds (e.g., Muser et al., 2020;De Leeuw et al., 2020) and was used in this study for volcanic ash and SO 2 as well. The method assesses predefined objects For SAL comparison of Himawari-8 and ICON-ART data, we derived 6-hour averages and interpolated ash and SO 2 values onto a regular grid between 120 • W and 80 • E and 20 • N and 85 • N with a resolution of 0.1 • . However, before interpolation, 230 we applied a 5x5 pixel mean averaging to fill gaps in the satellite data whereby only values different from zero are considered.
Otherwise, the linear interpolation would have led to a loss of information when mapping on a coarser grid, because the regular grid is about 4 to 5 times coarser than the retrieval grid.
3 Results and discussions

Validation of mass loading 235
The Raikoke eruption 2019 injected ash and SO 2 up to 14 km into the atmosphere. Figure 3 shows mean ash (left) and SO 2 (right) column on June 22 (top row) and 23 (bottom row) in our reference simulation FPlume-rad. The volcanic plume first spreads with westerly winds and is then dragged into a low pressure system over the Northern Pacific Ocean. In the mass loadings of both compounds, no clear horizontal separation of the ash and SO 2 plume is visible (compare Fig. 3 left and right side). However, we will further investigate the separation of ash and SO 2 due to radiation in Sect. 3.3 after we validated our   (Muser et al., 2020). We have closed these gaps as follows.
The maximum of total ash derived with ICON-ART coupled with FPlume in both experiments with and without radiation-250 aerosol interaction coincides very well with the Himawari-8 data (compare red and yellow curve with black curve). The total ash derived with Mastin (different MER but the same fine ash fractions and emission profile as in the FPlume experiments) overestimates the amount of ash during the first 12 h after the onset of the eruption (blue curve). Thus, neglecting meteorological effects and other plume-related processes in the case of the Raikoke eruption results in higher MER in the long continuous phase of the eruption and subsequently to an overestimation of the ash emissions into ICON-ART (Fig. 2).

255
All simulation experiments in Fig. 4a include aerosol dynamics and have correctly reproduced the fallout of particles as indicated by the decrease of ash after two days. From Fig. 4b, where the temporal development of the different modes is shown, we can conclude that the decrease of ash after two days is mainly due to coarse and giant mode particles. The total ash from the simulations with FPlume display the best agreement with the Himawari-8 data in this analysis (Fig. 4a). However, the other curves remain mostly within the error range of Himawari-8 data as well (gray shading). Thus, we conclude that the 260 online treatment of plume development improves the ash loading prediction during the first hours and days after the eruption.
Afterwards, the aerosol dynamical processes become more important, and the differences between the experiments decrease.

Validation of dispersion using SAL
For the quantitative validation of the forecast quality, we performed a SAL analysis based on the values of 6h averaged SO 2 and ash mass loadings. We compare the results of FPlume-rad experiment and Himawari-8 satellite data. Figure 5 shows the 265 values for the structure on the abscissa, the amplitude on the ordinate, and the location in colors.
The location of the SO 2 plume agrees very well between model and observations throughout the whole simulation period. This is shown by the location values which are close to zero. The structure and amplitude values are close to zero between 24 to 72 hours after the beginning of the simulation on June 21 at 12 UTC. Thus, there is a high agreement between model and observations during this period. However, during the first 24 hours, the model prediction shows higher amplitude values and 270 a low structure value, indicating a larger mass loading in the model and a less diffuse SO 2 cloud in the model compared to the satellite estimates. We argue that the discrepancy in the amplitude between model and observation during the first hours of the Raikoke eruption stems from the possible underestimation of SO 2 by the satellite retrievals due to the dense ash cloud covering the region around the volcano. Moreover, several studies found that including parametrizations for the gravitational spreading of volcanic plumes in dispersion models is important to reproduce main features of the volcanic eruption cloud (e.g.,

275
Costa et al., 2013;Webster et al., 2020). We cannot rule out that the gravitational spreading is inadequately in the beginning of the simulation due to our relatively coarse horizontal resolution of 13 km. This might partly explain the large difference in the structure between ICON-ART and Himawari SO 2 during the first hours. In addition, our simulation is also affected by the uncertainties of input parameters (e.g., start and end time of individual puffs, plume heights, exit conditions). The model also predicts the location of the ash plume very well (Fig. 5). However, the positive amplitude during the first 280 day and positive structure values throughout the whole simulation indicate that the modeled ash loading is too high at the beginning and more diffusive over the domain. Figure A2 (first and second column) compares all 6 h-mean ash loadings. The large spread of the modeled ash plume across large parts of the northern Pacific Ocean is not seen in the observations. We argue that Himawari-8 measurements of ash at this time are might be hampered by water and ice clouds overlapping and obscuring the ash plume. As for SO 2 , the S values for ash between 12 and 18 h are slightly negative and, thus, different in sign from the 285 other time steps.

Vertical separation of SO 2 and ash plume
In this section, we discuss the evolution of the ash and SO 2 plume top heights and focus on the radiative effects on the plume dynamics. Figure 6 shows the top height of the ash and SO 2 plume for the FPlume-rad and FPlume-norad experiment (a), the resulting 290 vertical temperature difference on June 23, 12 UTC (b), and the vertical distribution of SO 2 and ash mixing ratios on the same date (c). We picked this particular day, because it allows a direct comparison to Fig. 8 in Muser et al. (2020), which only shows the ash plume top height. Our definition of the plume top height is the same as in Muser et al. (2020): we separated the ash and SO 2 plume from the background concentrations by applying a threshold value above which a model grid box is considered as volcanic plume. The threshold values are 0.01 µg kg −1 , 1 µg kg −1 , and 10 ppm for the accumulation ash mode, coarse ash 295 modes, and SO 2 , respectively. The lines for the plume top height and mass averaged height are smoothed by a Savitzky-Golay filter to remove 'steps' due to the low vertical resolution in upper atmospheric model levels. During the first hours after the beginning of the eruption, the plume top height for ash and SO 2 mainly rises due to the higher eruption heights of later puffs. The gray bars, which indicate the eruption height, coincide well with the top height (Fig. 6a).
Shortly after the end of the long eruption phase, we clearly see a separation of the ash plume top height between the FPlume-300 rad and the FPlume-norad experiments. The effect of the ash lofting due to radiation was already investigated in detail by Muser et al. (2020) with the same model system. They found that the absorption of shortwave and longwave radiation by the coated ash particles leads to the warming and rising of the ash plume. We compare the vertical profile of the temperature difference between the FPlume-rad and FPlume-norad case here with the vertical temperature differences in Muser et al. (2020) on June 23, 12 UTC. A single large positive anomaly of approximately 0.3 K near 11 km occurs in our simulation (Fig. 6b).

305
Subsequently, the profile of the vertical mixing ratio, which in the FPlume-norad case still shows the overlap of the different phase dependent profiles, smooths out in the FPlume-rad case. Additionally, the whole ash plume rises to higher altitudes. In contrast, Muser et al. (2020) found two distinct temperature anomaly peaks around 10 km and 14 km in the order of 0.25 K each, which result in the formation of two maxima in the ash mixing ratios near 10 km and 15 km. The resulting uplift during the first 12 h in the ash plume is reduced by about 24% in our simulation compared to Muser et al. (2020).

310
In the first hours during and after the eruption, the absorption-induced warming of the ash cloud also causes the SO 2 plume to rise in FPlume-rad ( Fig. 6a and c). However, as SO 2 itself absorbs neither solar nor terrestrial radiation in our model setup, the ash plume top height clearly separates from the SO 2 plume top height with increasing time (Fig. 6a). The vertical profiles of the SO 2 mixing ratio in the FPlume-rad and FPlume-norad case indicate that radiation interaction smooths and reduces the vertical gradient of the SO 2 mixing ratios in the troposphere. In the stratosphere, a second peak occurs above the maximum 315 emission height (Fig. 6c).
The evolution of the mass-averaged height of the ash and SO 2 plume indicate an opposite behavior than the plume top height. The mass average of the SO 2 plume is generally higher than for the ash plume. In Fig. 7a and c, the vertical distribution of the total ash mass and the total SO 2 confirm that the SO 2 plume is about 4 km higher on average than the ash plume after three days. This is in agreement with several existing studies (e.g., Robock, 2000;Timmreck, 2012), which emphasized a fast 320 removal of ash after volcanic eruptions related to the higher weight of ash particles compared to SO 2 . The step-wise reduction of the ash in the mass averaged height is related to the loss of the giant mode during the first 24 h and the a large fallout of the coarse mode until about 50 h relative to simulation start ( Fig. 4b and Fig. A3).
For the overall median radius R, we consider the five ash modes i = 1, 5 (insoluble and mixed accumulation mode, insoluble and mixed coarse mode, and giant mode) and calculated R at every grid cell. r m is the median radius from the log normal distribution and N is the number of particles per grid box. The vertical distribution of the horizontally averaged median radius in Fig. 7 also shows the loss of the larger particles (coarse and giant mode) during the first 24 h. Afterwards, the values of the 330 mean radius are below 1.5 µm with a maximum around 5 to 6 km. In comparison to the FPlume-norad experiment, the radius is higher on average because aerosol radiation-interaction slows down the removal of larger particles from the atmosphere (Fig.   7d). This effect is especially visible for the removal of the coarse mode which is reduced and delayed of about 2 h in FPlumerad (maximum around 24 h in Fig. 7d and Fig. 4b). The temporal removal of the giant mode only shows small differences between FPlume-rad and FPlume-norad (Fig. 4b). The larger mean radius after approximately 42 h is related to an increasing 335 removal of accumulation mode particles in FPlume-norad compared with FPlume-rad ( Fig. 7d and Fig. 4b). However, we will leave a detailed analysis of the processes changing particle radii to further work.
Both the main SO 2 mass and the main ash mass are restricted to a narrower vertical range after three simulation days compared the end of the eruption (around 13 h after simulation start). The location of the SO 2 is between 8 and 14 km and between 4 and 10 km for ash (Fig. 7). Thus, for initializing long-range and climate simulations, a release of SO 2 and ash at 340 these altitudes is justified if the sedimentation during the first hours is considered in the total emission rate.
Despite of the clear vertical separation of the ash and SO 2 plume, the horizontal separation in model and observations is only small in the first three days after the eruption. Nevertheless, a strong vertical wind shear can result in the horizontal separation of the ash and SO 2 plume for longer time scales as in Kloss et al. (2021).

345
We investigated the Raikoke 2019 eruption, which was characterized by nine puffs and one continuous eruption phase of almost 3-hour-duration. Here, we describe a model setup in which the eruption source parameters were improved by (1) coupling ICON-ART with FPlume to account for the effect of changing volcanic and meteorological conditions and (2) a delineation of eruption phases ('puffs'). We further investigated the effect of radiation-aerosol interaction on the SO 2 plume due to a warming of the ash plume. The main outcomes are: 350 1. We demonstrated a large improvement of the total ash burden forecast in the first 12 h by resolving the individual puffs of the Raikoke eruption, which reduces the ash mass overestimation from 37% to 25%. Additionally, the online calculation