Articles | Volume 22, issue 5
Research article
16 Mar 2022
Research article |  | 16 Mar 2022

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

Julia Bruckert, Gholam Ali Hoshyaripour, Ákos Horváth, Lukas O. Muser, Fred J. Prata, Corinna Hoose, and Bernhard Vogel

In June 2019, the Raikoke volcano, Kuril Islands, emitted 0.4–1.8×109 kg of very fine ash and 1–2×109 kg of SO2 up to 14 km into the atmosphere. The eruption was characterized by several eruption phases of different duration and height summing up to a total eruption length of about 5.5 h. 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 with the Aerosols and Reactive Trace gases module) with the 1D plume model FPlume to calculate the eruption source parameters (ESPs) online. The main inputs are the 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, and coagulation) explain the loss of ash in the atmosphere 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 for SO2 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 (except for the structure value of ash), indicating a very good agreement of the model and observations. Furthermore, we found a separation of the ash and SO2 plume after 1 d 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 the case of complex and noncontinuous volcanic eruptions like that of Raikoke in 2019.

1 Introduction

Explosive volcanic eruptions inject particulate matter and gases into the atmosphere, which are then dispersed by atmospheric transport processes. Volcanic ash can remain airborne for up to a few months (e.g., Robock2000; Niemeier et al.2009; Jensen et al.2018) and drift away several thousand kilometers from the eruption point when emitted into the lower stratosphere or higher. Moreover, oxidation of volcanogenic SO2 leads to the formation of sulfate particles as secondary volcanic aerosols. These aerosols enhance the scattering of sunlight in the atmosphere and, thus, reduce the incoming shortwave radiation reaching the ground (e.g., Robock2000). Therefore, while the dispersion of ash particles mainly matters for aviation safety (Casadevall1994), regional public health (Horwell and Baxter2006), and local environment and infrastructure (e.g., Spence et al.2005; Stewart et al.2006; Wardman et al.2012), the fate of SO2 is crucial for predicting the impacts of volcanism on weather and climate (Robock2000; Mather2008; Timmreck2012; von Savigny et al.2020).

Forecasting the dispersion of volcanic aerosols in the atmosphere relies on the representation of both the source and sink parameters and processes. It has been shown that eruption source parameters (ESPs) such as the mass eruption rate (MER), the initial plume height, the emission profile, and the duration of the eruption can strongly influence the quality of the forecast of the spatial distribution of the volcanogenic gases and particles (e.g., Scollo et al.2008; Harvey et al.2018). The plume height can be estimated instantaneously by visual, radar- and lidar-based, or satellite observations. Until a few hours after the onset of a volcanic eruption when such plume height observations become available, the MER usually remains uncertain. Estimates of the MER include empirical parametrizations based on plume height (e.g., Mastin et al.2009) partially corrected by wind effects (e.g., Degruyter and Bonadonna2012; Woodhouse et al.2013) or are derived with 1D plume models (e.g., Folch et al.2016). Further uncertainties arise from the choice of the eruption profile (e.g., De Leeuw et al.2021), i.e., the vertical distribution of mass. Different approaches exist to parameterize the emission profile, e.g., idealized profiles (Stuefer et al.2013), plume-theory-based profiles (Marti et al.2017), Gaussian-shaped profiles derived from backward-trajectory modeling (Rieger et al.2015), constant profiles (e.g., Beckett et al.2020; Muser et al.2020), or more complex ones derived from the observations (e.g., De Leeuw et al.2021). However, in all parametrizations of the ESPs the volcanic plume dispersion remains decoupled from unresolved volcanic eruption dynamics including also the influence of the atmosphere on the emission height. This accounts for large uncertainties in modeling studies at regional to global scales (Textor et al.2005; Timmreck2012; von Savigny et al.2020). Marti et al. (2017) overcame this issue by coupling the NMMB-MONARCH-ASH transport model (Nonhydrostatic Multiscale Model on the B grid–Multiscale Online Nonhydrostatic AtmospheRe CHemistry model–ash) with the 1D plume model FPlume, which calculates the MER and the mass distribution in the column online. Another example is the study by Collini et al. (2013), who combined the WRF–ARW (Weather Research and Forecasting–Advanced Research WRF) forecast system with FALL3D and highlighted a good agreement in ash transport simulations with satellite observations for the 2011 Cordón Caulle eruption. Plu et al. (2021) simulated the 2010 Eyjafjallajökull eruption with the MOCAGE model (Modèle de Chimie Atmosphérique de Grande Echelle) and hourly changing MER from FPlume. They highlighted more concentrated ash concentrations in the horizontal and vertical scale, which more realistically represents the horizontal dispersion compared to parameterized MERs.

Only ash particles smaller than 32 µm (hereafter referred to as very fine ash) are relevant for long-range transport in the atmosphere (Rose and Durant2009). However, the amount of very fine ash emitted by a volcanic eruption is uncertain and depends on different parameters such as the strength and height of an eruption (Gouhier et al.2019), the composition of magma (Rose and Durant2009), and the availability of water (e.g., van Eaton et al.2012; Prata et al.2017). Gouhier et al. (2019) analyzed data of past volcanic eruptions with respect to the fraction of very fine ash in the whole mass erupted. They found that strong volcanic eruptions are less efficient in emitting very fine ash into the atmosphere possibly due to higher sedimentation within the plume. Most forecast models assume a fixed value for the fraction of very fine ash between 1 % (e.g., Muser et al.2020) and 5 % (e.g., Webster et al.2012; Beckett et al.2020) regardless of the strength of the eruption and lava composition. Throughout this paper, we use the term “plume” or “volcanic plume” to describe the part of the volcanic material dispersed in the atmosphere as commonly used in the meteorology instead of “volcanic cloud” to maintain a clear distinction from the meteorological clouds.

Once emitted into the atmosphere, aerosol dynamics (including aggregation) lead to a faster growth of particles and, thus, a quicker removal from the atmosphere (e.g., Brown et al.2012,  and references therein). Muser et al. (2020) investigated the impacts of aerosol dynamics and radiation interactions on the ash dispersion after the Raikoke eruption in June 2019. They showed that aerosol dynamical processes such as nucleation, condensation, and coagulation enhance the removal of the ash particles from the atmosphere. On the other hand, the absorption of incoming shortwave radiation by internally mixed aerosols leads to the lofting of the aerosol plume. The simulated lofting effect for the Raikoke eruption resulted in a 6 km rise of the plume top after the first 4 d (Muser et al.2020). Zhu et al. (2020) confirmed that coagulation of mixed particles with ash and sulfate is required to produce the evolution of the size distribution of mixed particles following the 2014 Kelud eruption. However, they further found that the initial SO2 lifetime is determined by direct SO2 uptake on ash, rather than its oxidation by OH.

Here, we aim to link complex ESPs during the first hours of volcanic eruptions to the fate of volcanogenic gases and aerosols. As a case study, we investigate the ash and SO2 dispersion of the Raikoke eruption (48.29 N, 153.24 E) on 21 and 22 June 2019 during the first 3 d after the eruption onset. The eruption was characterized by 10 eruption phases of 5 to 14 km height lasting between 5 min and 3 h (Horváth et al.2021b). Such complexity leads to further difficulties in deriving reasonable ESPs for plume dispersion forecasts. Throughout the paper, we define “eruption phase” as one distinct time period in which the volcano was erupting. The constant emission rate and profile used by Muser et al. (2020) caused an approximately 6 h time lag between the time series of modeled and observed ash mass loading. This gap might be filled by improving the representation of the ESPs and a varying very fine ash fraction according the relationship by Gouhier et al. (2019). Moreover, the impacts of source and sink processes on the fate of SO2 erupted from Raikoke remain unexplored in the study by Muser et al. (2020). Kloss et al. (2021) investigated SO2 transport following the 2019 Raikoke eruption with observations and models. They found enhanced stratospheric aerosol optical depths in the whole Northern Hemisphere for more than 1 year following the Raikoke eruption when using an SO2 setup which realistically represents the transport of volcanic compounds during the first hours after the Raikoke eruption. De Leeuw et al. (2021) found that simulating the correct burden of SO2 is sensitive to the fraction emitted into the lower stratosphere and therefore depends strongly on the emission profile chosen. In this work, we want to answer the following research questions. How large is the influence of resolving the eruption phases on the predicted ash mass loading after the Raikoke eruption? Can an online treatment of volcanic ESPs improve the predicted mass loading and dispersion of ash and SO2 plumes? And what is the impact of aerosol–radiation interaction on the dispersion of the SO2 plume? The paper is structured as follows: in Sect. 2, the methodology, including the model setup, the inputs, observations, and the validation method used, is described. In Sect. 3, we evaluate our experiments with respect to the mass loading and structure, amplitude, and location of the plume. In addition, we discuss the separation of the ash and SO2 plumes due to aerosol–radiation interaction. Finally, Sect. 4 concludes the paper.

2 Methods

2.1 ICON-ART modeling system

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 3D nonhydrostatic and compressible Navier–Stokes equations on an icosahedral grid and allows for seamless predictions from local to global scales (Zängl et al.2015; 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 dynamics in the troposphere and stratosphere (Rieger et al.2015; Weimer et al.2017; Schröter et al.2018). Muser et al. (2020) reported and demonstrated the 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 into seven lognormal 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 of the modes constant. For the Aitken mode, nucleation, condensation, and coagulation are considered, while the accumulation mode and coarse mode are affected by condensation and coagulation only. The shifting of particles into another mode occurs either when a threshold diameter is exceeded (shift into a larger mode) or when a mass threshold of soluble coating on insoluble particles is exceeded (shift from an insoluble to a 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, and optical properties. Changes of particle optical properties can further feed back on the radiation and atmospheric state. However, the effects of aerosol dynamics on atmospheric humidity and clouds are not considered yet.

For detailed descriptions of ICON, ART, and AERODYN, we here refer to the works by Zängl et al. (2015), Rieger et al. (2015) and Schröter et al. (2018), and Muser et al. (2020), respectively.

2.2 Coupling ICON-ART with FPlume

For a better estimation of the ESPs, we coupled ICON-ART online with the 1D 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 the 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 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., the total MER in the case of a given plume height (as here) or plume height in the case of a given MER. Thirdly, the fraction of very fine ash is determined based on plume height and the 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 from the normalized Suzuki distribution, which is the same as the one used by Marti et al. (2017):

(1) S ( z ) = E ( ( 1 - z H p ) exp ( 4.0 ( z H p - 1 ) ) ) 5.0 MAX ( S ( z ) ) ,

where S(z) is the numerator of the equation and describes the vertical emission profile, E is the emission rate of very fine ash, Hp is the plume top height, and z refers to the height in the plume. Equation (1) explains the shape of the emission profile used here, which is also plotted in Fig. A1 in comparison with other profiles. To ensure the correct total ash mass emission and units when the particles are released into ICON-ART at discrete point sources in each model layer between the bottom and top height of the plume, we further normalized Eq. (1) by the vertical integral of S(z) (Rieger et al.2015). 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 Durant2009). We only used the MER from FPlume and calculated the very fine ash fraction and emission profile independently instead of using the vertical distribution of mass from FPlume due to two main reasons. (1) Based on offline analysis we figured out that the mass profiles for the predefined bin sizes strongly depends on the assumption of the initial total grain size distribution (TGSD). As information on the TGSD is often lacking, using FPlume mass profiles leads to a less generic approach and large uncertainties. (2) The definition of ash modes in ICON-ART is only relevant for long-range transport in the atmosphere and differs from the TGSD at the vent. Thus, we would have to convert the FPlume size bins into ICON-ART modes which requires several assumptions and increases the uncertainty of the emissions.

More details on the initialization of the ash particles is given in the next section. Besides ash, we also emitted SO2. Different from the ash emission, we prescribed the MER of SO2 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 temporal SO2 emission is available. Yet, during volcanic eruptions in general it is possible that the ash and SO2 are emitted at different phases of the eruption (e.g., Thomas and Prata2011).

Figure 1Schematic of the setup with the coupling between the global ICON-ART model and the 1D plume rise model FPlume.


Besides meteorological data, FPlume needs estimates of the exit temperature, exit velocity, exit volatile fraction, and plume height to solve for the total MER. Our setup enables the definition of these parameters for multiple eruption phases meaning that the MER used for ICON-ART depends on the exit conditions for each plume phase and the meteorological conditions. When solving the plume dynamics for the MER knowing the plume height, FPlume performs the calculations iteratively for a range of possible MERs to reach the plume height wanted (Folch et al.2016).

2.2.1 Separation of the volcanic plume from the background

To study the evolution of the ash and SO2 plume in the atmosphere following the Raikoke eruption, we had to separate the ash and SO2 volcanic plume from the background mixing ratios due to technical reasons. In ICON-ART, we initialized the ash modes with 100 particles per kilogram air to avoid a division by zero in the diameter calculation routines. Thus, to especially study the behavior of the plume top, we used the following mixing ratio thresholds above which a grid cell is considered inside the plume and which are based on Muser et al. (2020): 0.01, 1, and 100 µg kg−1 for the accumulation, coarse, and giant modes, respectively. The corresponding threshold of SO2 is 10 ppm.

2.3 Eruption source parameters

2.3.1 Vent conditions

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 (Mastin2007); the exit velocity of the individual phase veph was a linear function of the plume height above the vent Hph between 14 000 and 4000 m, where the exit velocity was set to 150 and 90 m s−1, respectively. Thus, the following equation calculates veph based on Hph in meters:

(2) ve ph = 6 s - 1 H ph 1000 + 66 m s - 1 .

The resulting MERs are insensitive to the input vent conditions (temperature, velocity, and volatile fraction) in the range of 10 %.

The equation for the very fine ash fraction by Gouhier et al. (2019) depends on whether the SiO2 content is high or low and whether the conduit was opened or closed. As no information on the conduit has been available so far, we averaged the very fine ash fraction for low SiO2–closed conduit and low SiO2–opened conduit.

2.3.2 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 Earth ellipsoid of the individual eruption phases was estimated by a recently developed geometric technique (Horváth et al.2021a), which exploits the near-limb views provided by Geostationary Operational Environmental Satellite 17 (GOES-17). Such oblique observations offer close to orthogonal side views of vertical columns protruding from the Earth ellipsoid and thereby facilitate a simple height-by-angle method to derive point estimates of eruption column height in the vicinity of the vent. The GOES-17 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).

The Raikoke plume heights cannot be unambiguously determined by the traditional infrared brightness temperature method. For most eruption phases, the minimum 11 µm brightness temperature (BT11) falls within the narrow temperature range of the quasi-isothermal 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 eruption phases, on the other hand, the BT11 has a warm bias due to contributions from the warmer lower-level marine stratocumulus cloud 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. As FPlume requires the plume height above the vent, we converted the GOES-17 above-ellipsoid heights by subtracting a vent height of 550 m.

2.4 Model configuration

We performed global simulations with ICON-ART using a horizontal grid size of roughly 13.2 km (R3B07 grid) and 90 vertical levels up to 75 km. The global icosahedral grid of ICON ensures a uniform resolution across the globe. For each experiment, we simulated 72 h starting from 21 June at 12:00 UTC with initialized analysis data provided by the German Weather Service (DWD). During active eruption periods, the ESPs of Raikoke are calculated online with FPlume.

The 2019 Raikoke eruption was characterized by nine shorter eruption phases between 18:00 UTC on 21 June and 07:00 UTC on 22 June and one more or less continuous eruption phase between 22:40 UTC and 01:55 UTC. We performed three experiments. (1) In the reference experiment, FPlume calculates the ESPs with a varying very fine ash fraction and aerosol–radiation interaction activated in ICON-ART (FPlume-rad, Table 1). (2) The second experiment calculates the ESPs in the same way as above but neglects the interaction of aerosols and radiation (FPlume-norad). The comparison of FPlume-rad and FPlume-norad allows for quantifying the lifting of the volcanic plume due to radiation. (3) The third experiment derives the ESPs with the empirical relationship by Mastin et al. (2009), and it emits volcanic compounds with a prescribed very 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). The experiments FPlume-rad and FPlume-norad calculate the ESPs online within the simulation, whereas in Mastin-rad the ESPs are derived offline independent of the atmosphere and vent conditions. Table 1 summarizes the prescribed input parameters for the FPlume-rad experiment associated with the different eruption phases, which are fixed for the individual phases. The time limits for the phases and plume heights above the vent are based on satellite images from GOES-17, which are described in Sect. 2.3.2. Due to the 10 min temporal resolution of the GOES-17 data, the uncertainty in the start and end time of each individual eruption phase is smaller than ±5 min.

Table 1Model setup and input parameters for the individual eruption phases (FPlume-rad). The definition of the phases and plume heights above sea level (a.s.l.) are based on GOES-17 satellite observation as described in Sect. 2.3.2. The exit conditions are based on typical values of basaltic eruptions as described in Sect. 2.3.1. The SO2 mass emission rate is based on an observational estimate of the total SO2 mass following the 2019 Raikoke eruption from Muser et al. (2020), which was distributed over the individual phases with Eq. (3). This table only shows the values that are predefined and fixed for the individual phases. The temporally varying MER of the very fine ash, which is derived with FPlume and the relationship by Gouhier et al. (2019) and which is released into ICON-ART, is shown in Fig. 2.

Download Print Version | Download XLSX

Figure 2 shows the 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 on 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 relationship by Mastin et al. (2009) neglects atmospheric conditions and the very fine ash fraction is fixed within one phase, the MERs of the very fine ash are constant within each phase. The vertical profiles of the meteorological variables in Fig. A2a indicate increasing temperatures in most levels below 10 km during the long eruption phase between 9 and 14 h after simulation start. Additionally, the specific humidity increases by up to 1 g kg−1 in the lower 2 km (Fig. A2c). When warmer and moist air is entrained into the plume, the plume density reduces faster due to the lower ambient air density and the release of latent heat. This effect results in a higher positive buoyancy and a lower MER to reach a fixed height. In addition, the wind speed decreases in the lower 4 km between 9 and 14 h after simulation start, which reduces the plume bending, and subsequently the MER needed to reach a fixed height.

According to the MER values in Fig. 2, the total mass of very fine ash emitted in the model for all eruption phases together is about 1.21×109 kg using FPlume and 1.48×109 kg using Mastin-derived MERs.

Figure 2Mass eruption rate E for very fine ash (t s−1) calculated with the FPlume MER times the very fine ash fraction from Gouhier et al. (2019) (red) and calculated with the Mastin MER times the very fine ash fraction derived in the FPlume experiment (blue). The very fine ash fraction is the same in both experiments to allow for a direct comparison of the FPlume- and Mastin-derived MER on the dispersion in the atmosphere. Active eruption phases are indicated by the gray shading. Please note that the date format in this and following figures is month day, year.


The total mass of very fine ash is evenly distributed as insoluble tracers over the accumulation, coarse, and giant modes. The three insoluble modes are emitted as lognormal distributions with median diameters of 0.8, 2.98, and 11.35 µm, respectively. The standard deviation is 1.4 for each mode.

Following previous studies, we emitted a total of 1.5×109 kg SO2 (Muser et al.2020; Kloss et al.2021; De Leeuw et al.2021). However, the SO2 release is linearly adjusted to the eruption heights and length of each phase (Table 1) as follows:

(3) E ph = E SO 2 H ph H T ,

where Eph is the phase-dependent MER of SO2, ESO2 is the mean MER based on the observed amount of SO2 and the sum of the duration of all phases, Hph is the phase plume height (above the vent), and HT=11 571.2 m is the phase duration-weighted (t) mean plume height derived as

(4) H T = i = 1 10 H ph , i t i i = 1 10 t i .

Finally, SO2 release is vertically distributed according to the Suzuki profiles (comparison to previously used profiles in Fig. A1). Here, ICON-ART treats SO2 as a chemical tracer that can be oxidized by a simplified OH-chemistry scheme as presented in Weimer et al. (2017).

2.5 Himawari-8 ash and SO2 retrievals

To validate the model results, we used column SO2 and ash mass loadings estimated from the 16-band visible and infrared Advanced Himawari Imager (AHI) onboard the Himawari-8 geostationary satellite at every full hour. Himawari-8 is operated by the Japan Aerospace Exploration Agency (JAXA) and the Japan 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, SO2 is retrieved by the AHI band centered near 7.3 µm, where the absorption of SO2 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 11.2 and 12.4 µm are considered. The lower detection threshold is <0.2g m−2 for volcanic ash. The ash retrievals were corrected by a mask that accounts for pixels that contain meteorological clouds but which were classified as completely cloud covered. Hereby, only pixels inside a 0.1 g m−2 contour line are considered, and a 9×9 median filter smooths out “spikes”.

2.6 SAL method

The SAL (structure, amplitude, and location) method is an object-based quality measure originally developed to verify precipitation forecasts (Wernli et al.2008, 2009). However, it has also been successfully applied for transport forecasts of volcanic compounds (e.g., Muser et al.2020; De Leeuw et al.2021) and was used in this study for volcanic ash and SO2 as well. The SAL method evaluates modeled and observed data according to their structure (S), amplitude (A), and location (L). While it assesses predefined objects based on a threshold value for the S and L component, the A component is a normalized domain-averaged quantity. The structure component S compares model and observations with respect to the volume of the defined objects. The value ranges between 2 and 2. Positive values indicate objects that are too large and/or too flat, whereas negative values indicate objects that are too small and/or too peaked. A value of zero refers to a perfect forecast with respect to the structure. The amplitude component A evaluates the domain-averaged relative deviation of the forecasts from observations, and it is positive when the model overestimates the predicted quantity and vice versa (it also ranges between 2 and 2). For a perfect forecast of the amplitude, A is zero. The location component L consists of two parts: L1 describes the agreement between the forecast and observation in terms of the normalized difference between the centers of mass, whereas L2 refers to the average distance between the center of mass of all objects and the individual objects. Both L1 and L2 can reach values between 0 and 1 so that L in total can have values between 0 and 2 with a perfect forecast with respect to the location at L=0 (Wernli et al.2008, 2009).

For the SAL comparison of Himawari-8 and ICON-ART data, we derived 6 h averages from both datasets at every full hour. Furthermore, we interpolated ash and SO2 values onto a regular grid at 120 W–80 E and 20–85 N with a resolution of 0.1. However, before interpolation, we applied a 5×5 pixel mean averaging to fill gaps in the mapped satellite data, considering only values different from zero. 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. These gaps in the satellite data arise during mapping from the native format onto a regular latitude–longitude grid as needed for the SAL analysis and are due to the increasing pixel sizes towards the edges of the retrieval domain.

To define objects in the SAL analysis, we used a threshold of 0.2 g m−2 for modeled and observed ash because this is the detection threshold for the Himawari-8 ash retrievals. For SO2, a threshold of 2.5 g m−2 for the model and observations is used to remove background SO2 concentrations in Himawari-8 data. This was necessary because we did not initialize the model with realistic background conditions and, therefore, can only compare the observed and modeled SO2 plume from the eruption.

3 Results and discussions

3.1 Validation of mass loading

The 2019 Raikoke eruption injected ash and SO2 up to 14 km into the atmosphere. Figure 3 shows mean ash (left) and SO2 (right) column on 22 June at 00:00–23:00 UTC (top row) and 23 June at 00:00–23:00 UTC (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 SO2 plume is visible (compare the left and right sides of Fig. 3). However, we will further investigate the separation of ash and SO2 due to radiation in Sect. 3.3 after we validated our setup.

Figure 3Simulated daily mean column mass loadings for ash (a, c) and SO2 (b, d) on 22 June 2019 at 00:00–23:00 UTC (a, b) and 23 June 2019 at 00:00–23:00 UTC (c, d) (g m−2). The results are based on the FPlume-rad experiment.

Figure 4a shows the temporal evolution of the ash loading in the atmosphere following the Raikoke eruption for different experiments and observations. The temporal resolution of the data is 1 h. The Himawari-8 data reveal a steep increase of ash mass at 22:00 UTC on 21 June until a peak of 1.0 Tg is reached at 05:00 UTC on 22 June, and the curve remains above 1.0 Tg for 5 h. The maximum at 07:00 UTC (22 June) of 1.1 Tg is followed by a descent to 0.3–0.5 Tg.

Muser et al. (2020) (green curve) emphasized that aerosol processes account for the ash removal. Nucleation, condensation, and coagulation increase the size of aerosol particles and, thus, lead to a faster sedimentation. However, Muser et al. (2020) were not able to quantitatively explain the time lag between the model and observations during the first hours of the eruption (18:00 UTC on 21 June until 03:00 UTC on 22 June). Besides, the continuous emission with a constant MER led to a slight overestimation of the ash mass loading (Muser et al.2020). We have closed these gaps as follows.

Figure 4(a) Temporal evolution of the total amount of very fine ash in the atmosphere (Tg) for Himawari-8 observations (black) with an estimated uncertainty range (gray shading). Simulation with one constant emission phase (green, after Muser et al.2020). Simulation with ICON-ART with and without aerosol–radiation interaction, each with phase-dependent FPlume MERs and a very fine ash fraction derived from the relationship by Gouhier et al. (2019) (red and yellow). Simulation with ICON-ART, eruption-phase-dependent Mastin MERs, and the same very fine ash fraction as for the red curve (blue). (b) Temporal evolution of the simulated mass for the different ash modes (dotted for accumulation modes, dash-dotted for coarse modes, and dashed for the giant mode). The colors refer to the model experiment shown in the panel (a).


The maximum of total ash derived with ICON-ART coupled with FPlume (online treatment) in both experiments with and without radiation–aerosol interaction coincides very well with the Himawari-8 data (in Fig. 4, compare the red and yellow curves with the black curve). The total ash derived with Mastin (a 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 (offline treatment), as is often done in volcanic dispersion forecasts, results in a higher MER especially in the long continuous phase of the eruption and subsequently increased ash emissions into ICON-ART (Fig. 2).

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 2 d. From Fig. 4b, where the temporal development of the different modes is shown, we can conclude that the decrease of ash after 2 d 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 online treatment of plume development improves the ash loading prediction during the first hours and days after the eruption. After about 30 h, the aerosol dynamical processes become more important, and the differences between the experiments decrease.

3.2 Validation of dispersion using SAL

For the quantitative validation of the forecast quality, we performed a SAL analysis using 6 h averaged SO2 and ash mass loadings. We compare the results of the FPlume-rad experiment and Himawari-8 satellite data. Figure 5 shows the values for the structure on the abscissa, the amplitude on the ordinate, and the location in colors. We do not discuss the SAL values for the FPlume-norad and the Mastin-rad case here. This is because FPlume-norad only shows very small differences to the FPlume-rad case in all SAL values, and Mastin-rad only changes the amplitude value, as only the MER is higher compared to FPlume-rad. Based on the analysis of hourly to daily mean values, we conclude that 6 h averages provide a reasonable compromise between both reproducing the details and reducing the amount of missing values and noise.

Figure 5SAL values for ash and SO2 in the reference experiment (FPlume-rad) with S (structure) on the abscissa, A (amplitude) on the ordinate, and L (location) in colors. The dotted lines connect consecutive hours to visualize the temporal development (black for SO2 and gray for ash). The labels of the hours are given relative to the start of the simulation on 21 June 2019 at 12:00 UTC.


The location of the SO2 plume agrees very well between the 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 h after the beginning of the simulation on 21 June at 12:00 UTC. Thus, there is a high agreement between the model and observations during this period. However, during the first 24 h, the model prediction shows higher amplitude values and a low structure value, indicating a larger mass loading in the model and a less diffuse SO2 plume in the model compared to the satellite estimates. We argue that the discrepancy in the amplitude between the model and observations during the first hours of the Raikoke eruption stems from the possible underestimation of SO2 by the satellite retrievals due to the dense ash plume covering the region around the volcano. In addition, our simulation is also affected by the uncertainties of input parameters (e.g., start and end time of individual eruption phases, plume heights, and exit conditions).

The model also predicts the location of the ash plume very well (Fig. 5). The positive structure values indicate that the modeled ash loading becomes more diffusive over the domain for most of the time. Figure A3 (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, which is the main reason for the high structure value in Fig. 5. We argue that Himawari-8 measurements of ash at this time might be hampered by water and ice clouds overlapping and obscuring the ash plume. This argument further explains why the temporal evolution of the Himawari-8 measurements in Fig. 4a shows variations between 30 and 60 h, although the emission from Raikoke ceased.

The high amplitude value for ash between 12 and 36 h, despite the almost perfect agreement in the total mass in Fig. 4, also stems from the larger spread of the ash plume in the beginning. The reason is that the background values are considered zero and that the amplitude in the SAL analysis, unlike the object-based structure and location values, is a domain-averaged quantity.

We have shown that the model setup realistically represents the amount of ash following the 2019 Raikoke eruption and that the dispersion of the ash and SO2 agrees well with the observations. In the next step, we analyze the vertical distribution of the ash and SO2 plume.

3.3 Vertical separation of the SO2 and ash plume

In this section, we discuss the evolution of the ash and SO2 plume top heights and focus on the radiative effects on the plume dynamics.

Figure 6(a) Temporal evolution of the SO2 (red) and ash (yellow) plume top height and mass-averaged height for the FPlume-rad (solid) and FPlume-norad (dashed) experiment smoothed by a Savitzky–Golay filter. The gray bars indicate the duration and height of the 10 individual eruption phases. They are based on the analysis of the GOES-17 data, which serve as inputs for FPlume. (b) Vertical profile of the temperature difference between FPlume-rad and FPlume-norad in the ash plume 48 h after the start of the simulation. (c) Vertical SO2 (red) and ash (yellow) profile after 48 h for the FPlume-rad (solid) and FPlume-norad (dashed) experiment.


Figure 6 shows the top height of the ash and SO2 plume for the FPlume-rad and FPlume-norad experiments (a), the resulting vertical temperature difference on 23 June at 12:00 UTC (b), and the vertical distribution of the SO2 and ash mixing ratios on the same date (c). The plume top height in (a) is defined as the maximum height of all grid cells in the plume that were separated from background mixing ratios as explained in Sect. 2.2.1. The average plume height in (a) is the mean height weighted by the mass of all grid cells that are considered inside the plume. The values in (b) and (c) were horizontally averaged over the whole detected plume, again excluding grid cells outside the plume. In (b) and (c), we picked 23 June at 12:00 UTC because it allows for a direct comparison to Fig. 8 in Muser et al. (2020), which only shows the ash plume top height. 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 at upper-atmospheric model levels. The difference in height between FPlume-rad and FPlume-norad remains similar, regardless of the use of this filter. However, the increasing plume height already starting before the beginning of the eruption is a result of the filtering.

During the first hours after the beginning of the eruption, the plume top height for ash and SO2 mainly rises due to the higher eruption heights of later eruption phases. 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-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 cases here with the vertical temperature differences in Muser et al. (2020) on 23 June at 12:00 UTC. A single large positive anomaly of approximately 0.3 K near 11 km occurs in our simulation (Fig. 6b). Subsequently, the whole ash plume rises to higher altitudes. In contrast, Muser et al. (2020) found two distinct temperature anomaly peaks around 10 and 14 km on the order of 0.25 K each, which result in the formation of two maxima in the ash mixing ratios near 10 and 15 km. The resulting uplift during the first 12 h in the ash plume in our simulation is about 33 % of the resulted lifting in Muser et al. (2020).

In the first hours during and after the eruption, the absorption-induced warming of the ash plume also causes the SO2 plume to rise in FPlume-rad (Fig. 6a and c). However, as SO2 itself absorbs neither solar nor terrestrial radiation in our model setup, the ash plume top height clearly separates from the SO2 plume top height with increasing time (Fig. 6a). The vertical profiles of the SO2 mixing ratio in the FPlume-rad and FPlume-norad cases indicate that radiation interaction smooths and reduces the vertical gradient of the SO2 mixing ratios in the troposphere. In the stratosphere, a second peak occurs above the maximum emission height (Fig. 6c).

The evolution of the mass-averaged height of the ash and SO2 plume indicates behavior opposite to that of the plume top height. The mass average of the SO2 plume is generally higher than for the ash plume. In Fig. 7a and c, the vertical distribution of the ash and SO2 mass concentrations (kg m−3) confirm that the SO2 plume is about 5 km higher on average than the ash plume after 3 d. This is in agreement with several existing studies (e.g., Timmreck2012; Robock2000), which emphasized a fast removal of ash after volcanic eruptions related to the higher weight of ash particles compared to SO2. The stepwise 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 large fallout of the coarse mode until about 50 h relative to simulation start (Figs. 4b and A4).

Figure 7(a, c) Temporal evolution of the horizontally averaged vertical distribution of ash mass concentration (a) and SO2 mass concentration (c) in the FPlume-rad experiment. (b, d) Number-weighted characteristic ash radius (Eq. 5) for FPlume-rad (b) and difference between FPlume-rad and FPlume-norad (d).


In the next step, we want to compare the vertical distribution of a characteristic ash particle radius R, which we calculated as follows:

(5) R = i = 1 5 r m , i N i i = 1 5 N i .

For the overall characteristic radius R, we consider the five ash modes i=1,5 (insoluble and mixed accumulation modes, insoluble and mixed coarse modes, and giant mode) and calculated R at every grid cell. rm is the median radius from the lognormal distribution, and N is the number of particles per grid box. The vertical distribution of the horizontally averaged characteristic radius in Fig. 7b also shows the loss of the larger particles (coarse and giant modes) during the first 24 h. Afterwards, the values of the mean characteristic radius are below 1.5 µm with a maximum around 5 to 6 km. In comparison to the FPlume-norad experiment, the characteristic radius is higher on average compared to the FPlume-rad experiments because aerosol–radiation interaction slows down the removal of larger particles from the atmosphere (Fig. 7d). This effect is visible in the removal of the accumulation mode, which is reduced in FPlume-rad (Fig. 4b). Compared to FPlume-norad, the removal of the coarse modes in FPlume-rad is delayed by about 1 h between 18 and 25 h after simulation start (Fig. 4b). After about 30 h the amount of the coarse mode is slightly lower in FPlume-rad than in FPlume-norad because the particle radius is larger at most altitudes (Figs. 4b and 7d). The temporal removal of the giant mode only shows small differences between FPlume-rad and FPlume-norad (Fig. 4b). The larger mean characteristic radius after approximately 48 h is related to an increasing removal of accumulation mode particles in FPlume-norad compared with FPlume-rad (Figs. 7d and 4b). However, we will leave a detailed analysis of the processes changing particle radii to further work.

Both the main SO2 mass and the main ash mass are restricted to a narrower vertical range after 3 simulation days compared to the end of the eruption (around 13 h after simulation start). The location of the SO2 is between 8 and 14 km and between 4 and 7 km for ash (Fig. 7). Thus, for initializing long-range and climate simulations, a release of SO2 and ash at these altitudes is justified if the sedimentation during the first hours is considered in the total emission rate.

Despite the clear vertical separation of the ash and SO2 plume, the horizontal separation in both model and observations remains small in the first 3 d after the eruption. Nevertheless, a strong vertical wind shear can result in the horizontal separation of the ash and SO2 plume on longer timescales as in Kloss et al. (2021).

4 Conclusions

We investigated the 2019 Raikoke eruption, which was characterized by nine shorter eruption phases and one continuous eruption phase of almost 3 h duration. Here, we describe a model setup in which the ESPs 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. We further investigated the effect of radiation–aerosol interaction on the SO2 plume due to a warming of the ash plume. The main findings are the following.

  1. We demonstrated a large improvement of the total ash burden forecast in the first 12 h by resolving the individual eruption phases of the Raikoke eruption, which reduces the ash mass overestimation from 37 % to 18 %. Additionally, the online calculation of ESPs by FPlume further improves temporal evolution of the simulated ash mass, which shows an almost perfect agreement with the observed evolution of ash mass loading.

  2. In addition to the mass loading, the predicted spatial dispersion of the ash and SO2 plume also agrees well with observations from Himawari-8 as our SAL analysis reveals. However, we hypothesize that the validation of the simulated ash and SO2 dispersion is partially hampered by a dense ash plume in the beginning of the eruption and by overlapping water and ice clouds later on.

  3. As already demonstrated in Muser et al. (2020), aerosol–radiation interaction leads to a warming of the volcanic ash plume and, therefore, to a lofting of ash particles during the first hours. Here, we additionally found a lofting of the SO2 plume during the first 12 h after the eruption caused by a warming of the plume. However, with increasing time, the SO2 plume becomes more and more vertically separated from the ash plume, and the lofting slows down. This is related to a faster sedimentation of the ash particles compared to SO2 and the fact that SO2 does not absorb solar radiation in our model.

Appendix A

Figure A1Comparison of the Suzuki profiles used in this work with vertical emission profiles from previous studies (“VolRes” and “StratProfile” from De Leeuw et al.2021, and constant profiles from Muser et al.2020, and Kloss et al.2021). (a) Normalized with respect to the MER (E) but with actual emission heights (z). (b) Normalized with respect to the MER and emission height (z).


Figure A2Meteorological conditions above the vent in the first 20 h after simulation start to explain variations in the FPlume-derived MER. Four different atmospheric variables are shown: (a) temperature (C), (b) air density (kg m−1), (c) specific humidity (g kg−1), and (d) wind speed (m s−1). For each variable the temporal development in the vertical axis above the vent is given over time in contours and the difference between the two time steps indicated by the vertical dashed lines in the contour plot (later step minus earlier time step).


Figure A3Comparison plots for seven ash and SO2 6 h averaged column loadings (row 1 to 7) in order to explain the discrepancy between simulated ICON-ART data and Himawari-8 observed data in the SAL analysis. First and second column: ash column loadings (g m−2) from ICON-ART and Himawari-8. Third and fourth column: SO2 column loadings (g m−2) from ICON-ART and Himawari-8.

Figure A4Horizontally averaged ash particle characteristic radius of the lognormal distribution (µm) (left, Eq. 5) and horizontal sum of the particle number (right) of the different ash modes: insoluble accumulation (a, b), mixed accumulation (c, d), insoluble coarse (e, f), mixed coarse (g, h), and giant (i, j).


Code and data availability

The output from ICON-ART simulations performed in this study can be provided upon request by the corresponding author. The ICON-ART code is license protected and can be accessed upon request to the corresponding author. Himawari-8 AHI datasets that have been analyzed in the scope of this study can be provided upon request by the corresponding author.

Author contributions

JB, GAH, LOM, CH, and BV developed the ICON-ART code and carried out simulations. ÁH provided plume heights and durations based on GOES-17 data. FJP processed the Himawari-8 data and performed the ash and SO2 retrievals. JB and GAH prepared the paper with significant contributions and comments on the original draft from all authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Satellite observations, in situ measurements and model simulations of the 2019 Raikoke eruption (ACP/AMT/GMD inter-journal SI)”. It is not associated with a conference.


This research has been funded by the Deutsche Forschungsgemeinschaft (DFG) as part of the research unit VolImpact (FOR2820, DFG grant no. 398006378). The contributions are within the VolImpact sub-projects VolPlume (Julia Bruckert, Gholam Ali Hoshyaripour, Bernhard Vogel, and Ákos Horváth) and VolCloud (Corinna Hoose). JMA and JAXA are acknowledged for providing Himawari-8 data. We thank Fabio Crameri for the development of scientific color maps (Crameri2021) to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies (Crameri et al.2020). We also thank our reviewers, Arnau Folch, Sara Barsotti, and Leonardo Mingari, for their helpful suggestions and comments on the manuscript.

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (FOR2820, grant no. 398006378).

Review statement

This paper was edited by Peter Haynes and reviewed by Arnau Folch, Sara Barsotti, and Leonardo Mingari.


Beckett, F. M., Witham, C. S., Leadbetter, S. J., Crocker, R., Webster, H. N., Hort, M. C., Jones, A. R., Devenish, B. J., and Thomson, D. J.: Atmospheric dispersion modelling at the London VAAC: A review of developments since the 2010 eyjafjallajökull volcano ash cloud, Atmosphere, 11, 352,, 2020. a, b

Brown, R. J., Bonadonna, C., and Durant, A. J.: A Review of Volcanic Ash Aggregation, Phys. Chem. Earth Pt. A/B/C, 45–46, 65–78,, 2012. a

Casadevall, T. J.: Volcanic ash and aviation safety: Proceedings of the first international symposium on volcanic ash and aviation safety, U.S. Geological Survey, report,, 1994. a

Collini, E., Osores, M. S., Folch, A., Viramonte, J., Villaosa, G., and Salmuni, G.: Volcanic ash forecast during the June 2011 Cordón Caulle eruption, Nat. Hazards, 66, 389–412,, 2013. a

Crameri, F.: Scientific colour maps Version 7.0.0 (February 2021),, 2021. a

Crameri, F., Shephard, G. E., and Heron, P. J.: The misuse of colour in science communication, Nat. Commun., 11, 5444,, 2020. a

de Leeuw, J., Schmidt, A., Witham, C. S., Theys, N., Taylor, I. A., Grainger, R. G., Pope, R. J., Haywood, J., Osborne, M., and Kristiansen, N. I.: The 2019 Raikoke volcanic eruption – Part 1: Dispersion model simulations and satellite retrievals of volcanic sulfur dioxide, Atmos. Chem. Phys., 21, 10851–10879,, 2021. a, b, c, d, e, f

Degruyter, W. and Bonadonna, C.: Improving on mass flow rate estimates of volcanic eruptions, Geophys. Res. Lett., 39, L16308,, 2012. a

Folch, A., Costa, A., and Macedonio, G.: FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation, Geosci. Model Dev., 9, 431–450,, 2016. a, b, c, d

Giorgetta, M. A., Brokopf, R., Crueger, T., Esch, M., Fiedler, S., Helmert, J., Hohenegger, C., Kornblueh, L., Köhler, M., Manzini, E., Mauritsen, T., Nam, C., Raddatz, T., Rast, S., Reinert, D., Sakradzija, M., Schmidt, H., Schneck, R., Schnur, R., Silvers, L., Wan, H., Zängl, G., and Stevens, B.: ICON-A, the Atmosphere Component of the ICON Earth System Model: I. Model Description, J. Adv. Model. Earth Sy., 10, 1613–1637,, 2018. a

Gouhier, M., Eychenne, J., Azzaoui, N., Guillin, A., Deslandes, M., Poret, M., Costa, A., and Husson, P.: Author Correction: Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere, Sci. Rep., 9, 1–12,, 2019. a, b, c, d, e, f, g, h

Harvey, N. J., Huntley, N., Dacre, H. F., Goldstein, M., Thomson, D., and Webster, H.: Multi-level emulation of a volcanic ash transport and dispersion model to quantify sensitivity to uncertain parameters, Nat. Hazards Earth Syst. Sci., 18, 41–63,, 2018. a

Heinze, R., Dipankar, A., Henken, C. C., Moseley, C., Sourdeval, O., Trömel, S., Xie, X., Adamidis, P., Ament, F., Baars, H., Barthlott, C., Behrendt, A., Blahak, U., Bley, S., Brdar, S., Brueck, M., Crewell, S., Deneke, H., Di Girolamo, P., Evaristo, R., Fischer, J., Frank, C., Friederichs, P., Göcke, T., Gorges, K., Hande, L., Hanke, M., Hansen, A., Hege, H.-C., Hoose, C., Jahns, T., Kalthoff, N., Klocke, D., Kneifel, S., Knippertz, P., Kuhn, A., van Laar, T., Macke, A., Maurer, V., Mayer, B., Meyer, C. I., Muppa, S. K., Neggers, R. A. J., Orlandi, E., Pantillon, F., Pospichal, B., Röber, N., Scheck, L., Seifert, A., Seifert, P., Senf, F., Siligam, P., Simmer, C., Steinke, S., Stevens, B., Wapler, K., Weniger, M., Wulfmeyer, V., Zängl, G., Zhang, D., and Quaas, J.: Large-eddy simulations over Germany using ICON: a comprehensive evaluation, Q. J. Roy. Meteor. Soc., 143, 69–100,, 2017. a

Horváth, Á., Carr, J. L., Girina, O. A., Wu, D. L., Bril, A. A., Mazurov, A. A., Melnikov, D. V., Hoshyaripour, G. A., and Buehler, S. A.: Geometric estimation of volcanic eruption column height from GOES-R near-limb imagery – Part 1: Methodology, Atmos. Chem. Phys., 21, 12189–12206,, 2021a. a

Horváth, Á., Girina, O. A., Carr, J. L., Wu, D. L., Bril, A. A., Mazurov, A. A., Melnikov, D. V., Hoshyaripour, G. A., and Buehler, S. A.: Geometric estimation of volcanic eruption column height from GOES-R near-limb imagery – Part 2: Case studies, Atmos. Chem. Phys., 21, 12207–12226,, 2021b. a, b, c

Horwell, C. J. and Baxter, P. J.: The respiratory health hazards of volcanic ash: A review for volcanic risk mitigation, B. Volcanol., 69, 1–24,, 2006. a

Jensen, E. J., Woods, S., Lawson, R. P., Bui, T. P., Pfister, L., Thornberry, T. D., Rollins, A. W., Vernier, J.-P., Pan, L. L., Honomichl, S., and Toon, O. B.: Ash Particles Detected in the Tropical Lower Stratosphere, Geophys. Res. Lett., 45, 11483–11489,, 2018. a

Kloss, C., Berthet, G., Sellitto, P., Ploeger, F., Taha, G., Tidiga, M., Eremenko, M., Bossolasco, A., Jégou, F., Renard, J.-B., and Legras, B.: Stratospheric aerosol layer perturbation caused by the 2019 Raikoke and Ulawun eruptions and their radiative forcing, Atmos. Chem. Phys., 21, 535–560,, 2021. a, b, c, d

Macedonio, G., Costa, A., and Folch, A.: Uncertainties in volcanic plume modeling: A parametric study using FPLUME, J. Volcanol. Geoth. Res., 326, 92–102,, 2016. a

Marti, A., Folch, A., Jorba, O., and Janjic, Z.: Volcanic ash modeling with the online NMMB-MONARCH-ASH v1.0 model: model description, case simulation, and evaluation, Atmos. Chem. Phys., 17, 4005–4030,, 2017. a, b, c

Mastin, L. G.: A user-friendly one-dimensional model for wet volcanic plumes, Geochem. Geophy. Geosy., 8, Q03014,, 2007. a

Mastin, L. G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., Durant, A., Ewert, J. W., Neri, A., Rose, W. I., Schneider, D., Siebert, L., Stunder, B., Swanson, G., Tupper, A., Volentik, A., and Waythomas, C. F.: A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions, J. Volcanol. Geoth. Res., 186, 10–21,, 2009. a, b, c

Mather, T. A.: Volcanism and the atmosphere: the potential role of the atmosphere in unlocking the reactivity of volcanic emissions, Philos. T. R. Soc. A, 366, 4581–4595,, 2008. a

Morton, B. R., Taylor, G., and Turner, J. S.: Turbulent Gravitational Convection from Maintained and Instantaneous Sources, P. R. Soc. Lond. A, 234, 1–23,, 1956. a

Muser, L. O., Hoshyaripour, G. A., Bruckert, J., Horváth, Á., Malinina, E., Wallis, S., Prata, F. J., Rozanov, A., von Savigny, C., Vogel, H., and Vogel, B.: Particle aging and aerosol–radiation interaction affect volcanic plume dispersion: evidence from the Raikoke 2019 eruption, Atmos. Chem. Phys., 20, 15015–15036,, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y

Niemeier, U., Timmreck, C., Graf, H.-F., Kinne, S., Rast, S., and Self, S.: Initial fate of fine ash and sulfur from large volcanic eruptions, Atmos. Chem. Phys., 9, 9043–9057,, 2009. a

Plu, M., Bigeard, G., Sič, B., Emili, E., Bugliaro, L., El Amraoui, L., Guth, J., Josse, B., Mona, L., and Piontek, D.: Modelling the volcanic ash plume from Eyjafjallajökull eruption (May 2010) over Europe: evaluation of the benefit of source term improvements and of the assimilation of aerosol measurements, Nat. Hazards Earth Syst. Sci., 21, 3731–3747,, 2021. a

Prata, A., Rose, W., Self, S., and O'Brien, D.: Global, Long-Term Sulphur Dioxide Measurements from TOVS Data: A New Tool for Studying Explosive Volcanism and Climate, American Geophysical Union (AGU), 75–92,, 2004. a

Prata, F., Woodhouse, M., Huppert, H. E., Prata, A., Thordarson, T., and Carn, S.: Atmospheric processes affecting the separation of volcanic ash and SO2 in volcanic eruptions: inferences from the May 2011 Grímsvötn eruption, Atmos. Chem. Phys., 17, 10709–10732,, 2017. a

Rieger, D., Bangert, M., Bischoff-Gauss, I., Förstner, J., Lundgren, K., Reinert, D., Schröter, J., Vogel, H., Zängl, G., Ruhnke, R., and Vogel, B.: ICON–ART 1.0 – a new online-coupled model system from the global to regional scale, Geosci. Model Dev., 8, 1659–1676,, 2015. a, b, c, d

Robock, A.: Volcanic Eruptions and Climate, Rev. Geophys., 38, 191–219, 2000. a, b, c, d

Rose, W. I. and Durant, A. J.: Fine ash content of explosive eruptions, J. Volcanol. Geoth. Res., 186, 32–39,, 2009. a, b, c

Schröter, J., Rieger, D., Stassen, C., Vogel, H., Weimer, M., Werchner, S., Förstner, J., Prill, F., Reinert, D., Zängl, G., Giorgetta, M., Ruhnke, R., Vogel, B., and Braesicke, P.: ICON-ART 2.1: a flexible tracer framework and its application for composition studies in numerical weather forecasting and climate simulations, Geosci. Model Dev., 11, 4043–4068,, 2018. a, b

Scollo, S., Folch, A., and Costa, A.: A parametric and comparative study of different tephra fallout models, J. Volcanol. Geoth. Res., 176, 199–211,, 2008. a

Spence, R. J. S., Kelman, I., Calogero, E., Toyos, G., Baxter, P. J., and Komorowski, J.-C.: Modelling expected physical impacts and human casualties from explosive volcanic eruptions, Nat. Hazards Earth Syst. Sci., 5, 1003–1015,, 2005. a

Stewart, C., Johnston, D., Leonard, G., Horwell, C., Thordarson, T., and Cronin, S.: Contamination of water supplies by volcanic ashfall: A literature review and simple impact modelling, J. Volcanol. Geoth. Res., 158, 296–306, 2006. a

Stuefer, M., Freitas, S. R., Grell, G., Webley, P., Peckham, S., McKeen, S. A., and Egan, S. D.: Inclusion of ash and SO2 emissions from volcanic eruptions in WRF-Chem: development and some applications, Geosci. Model Dev., 6, 457–468,, 2013. a

Textor, C., Graf, H.-F., Longo, A., and Neri, A.: Numerical simulation of explosive volcanic eruptions from the conduit flow to global atmospheric scales, Ann. Geophys., 48,, 2005. a

Thomas, H. E. and Prata, A. J.: Sulphur dioxide as a volcanic ash proxy during the April–May 2010 eruption of Eyjafjallajökull Volcano, Iceland, Atmos. Chem. Phys., 11, 6871–6880,, 2011. a

Timmreck, C.: Modeling the climatic effects of large explosive volcanic eruptions, WIRes Clim. Change, 3, 545–564,, 2012. a, b, c

van Eaton, A. R., Muirhead, J. D., Wilson, C. J., and Cimarelli, C.: Growth of volcanic ash aggregates in the presence of liquid water and ice: An experimental approach, B. Volcanol., 74, 1963–1984,, 2012. a

von Savigny, C., Timmreck, C., Buehler, S. A., Burrows, J. P., Giorgetta, M., Hegerl, G., Horvath, A., Hoshyaripour, G. A., Hoose, C., Quaas, J., Malinina, E., Rozanov, A., Schmidt, H., Thomason, L., Toohey, M., and Vogel, B.: The research unit volimpact: Revisiting the volcanic impact on atmosphere and climate – preparations for the next big volcanic eruption, Meteorol. Z., 29, 3–18,, 2020. a, b

Wardman, J. B., Wilson, T. M., Bodger, P. S., Cole, J. W., and Stewart, C.: Potential impacts from tephra fall to electric power systems: A review and mitigation strategies, B. Volcanol., 74, 2221–2241,, 2012.  a

Webster, H. N., Thomson, D. J., Johnson, B. T., Heard, I. P. C., Turnbull, K., Marenco, F., Kristiansen, N. I., Dorsey, J., Minikin, A., Weinzierl, B., Schumann, U., Sparks, R. S. J., Loughlin, S. C., Hort, M. C., Leadbetter, S. J., Devenish, B. J., Manning, A. J., Witham, C. S., Haywood, J. M., and Golding, B. W.: Operational prediction of ash concentrations in the distal volcanic cloud from the 2010 Eyjafjallajökull eruption, J. Geophys. Res.-Atmos., 117, D00U08,, 2012. a

Weimer, M., Schröter, J., Eckstein, J., Deetz, K., Neumaier, M., Fischbeck, G., Hu, L., Millet, D. B., Rieger, D., Vogel, H., Vogel, B., Reddmann, T., Kirner, O., Ruhnke, R., and Braesicke, P.: An emission module for ICON-ART 2.0: implementation and simulations of acetone, Geosci. Model Dev., 10, 2471–2494,, 2017. a, b

Wernli, H., Paulat, M., Hagen, M., and Frei, C.: SAL – A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts, Mon. Weather Rev., 136, 4470–4487,, 2008. a, b

Wernli, H., Hofmann, C., and Zimmer, M.: Spatial Forecast Verification Methods Intercomparison Project: Application of the SAL Technique, Weather Forecast., 24, 1472–1484,, 2009. a, b

Woodhouse, M. J., Hogg, A. J., Phillips, J. C., and Sparks, R. S. J.: Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland, J. Geophys. Res.-Sol. Ea., 118, 92–109,, 2013. a

Zängl, G., Reinert, D., Rípodas, P., and Baldauf, M.: The ICON (ICOsahedral Non-hydrostatic) modelling framework of DWD and MPI-M: Description of the non-hydrostatic dynamical core, Q. J. Roy. Meteor. Soc., 141, 563–579,, 2015. a, b

Zhu, Y., Toon, O. B., Jensen, E. J., Bardeen, C. G., Mills, M. J., Tolbert, M. A., Yu, P., and Woods, S.: Persisting volcanic ash particles impact stratospheric SO2 lifetime and aerosol optical properties, Nat. Commun., 11, 1–11,, 2020. a

Short summary
Volcanic emissions endanger aviation and public health and also influence weather and climate. Forecasting the volcanic-plume dispersion is therefore a critical yet sophisticated task. Here, we show that explicit treatment of volcanic-plume dynamics and eruption source parameters significantly improves volcanic-plume dispersion forecasts. We further demonstrate the lofting of the SO2 due to a heating of volcanic particles by sunlight with major implications for volcanic aerosol research.
Final-revised paper