Particle Aging and Aerosol–Radiation Interaction Affect Volcanic Plume Dispersion: Evidence from Raikoke Eruption 2019

A correct and reliable forecast of volcanic plume dispersion is vital for aviation safety. This can only be achieved by representing all responsible physical and chemical processes (sources, sinks, and interactions) in the forecast models. The representation of the sources has been enhanced over the last decade, while the sinks and interactions have received less attention. In particular, aerosol dynamic processes and aerosol–radiation interaction are neglected so far. Here we address this gap by further developing the ICON-ART (ICOsahedral Nonhydrostatic – Aerosols and Reactive Trace gases) global modelling 5 system to account for these processes. We use this extended model for the simulation of volcanic aerosol dispersion after the Raikoke eruption in June 2019. Additionally, we validate the simulation results with measurements from AHI (Advanced Himawari Imager), CALIOP (CloudAerosol Lidar with Orthogonal Polarization), and OMPS-LP (Ozone Mapping and Profiling Suite – Limb Profiler). Our results show that around 50 % of very fine volcanic ash mass (particles with diameter d < 30 μm) is removed due to particle growth 10 and aging. Furthermore, the maximum volcanic cloud top height rises more than 6 km over the course of 4 days after the eruption due to aerosol–radiation interaction. This is the first direct evidence that shows how cumulative effects of aerosol dynamics and aerosol–radiation interaction lead to a more precise forecast of very fine ash lifetime in volcanic clouds.

Based on the AHI measurements the total emitted very fine ash mass (d < 32 µm) ranges between 0.4-1.8 × 10 9 kg, the SO 2 mass between 1-2 × 10 9 kg. The latter agrees well with the TROPOMI measurement in Sect. 2.1.1. and 03:10 UTC (Aqua and NOAA-20) when a brownish-colored and still localized plume was largely distinguishable from white/gray meteorological clouds in visible channel images. MODIS cloud top height, available at 1 km horizontal resolution, is obtained by matching the retrieved cloud top pressure to a Numerical Weather Prediction (NWP) geopotential height profile (Menzel et al., 2015). For the Raikoke plume, classified essentially as ice phase with a few liquid phase pixels, cloud top pressure was mostly determined by the CO 2 -slicing technique from channels near 13 µm and to a lesser degree by the infrared scheme (Gasch et al., 2017). This approach ensures full coupling and feedback between aerosol processes, radiation and the atmospheric state (Shao et al., 2011). Besides, a forward operator is implemented in the model to diagnose the attenuated 190 backscatter at the wavelengths 532 and 1064 nm (Hoshyaripour et al., 2019). To account for secondary aerosol formation and internally mixed aerosols, a new aerosol dynamics module is currently developed and implemented in ICON-ART. Details of this module are described in the following section.

Aerosol dynamics
The aerosol dynamics module (AERODYN) includes 10 log-normal modes that consider Aitken, accumulation and coarse 200 Figure 1. Chemical composition of the soluble (first row) and insoluble (second row) modes, mixing state of the modes (third row) and particle size distribution (giant mode is not shown). The dotted line represents a particle size distribution of soluble particles, the dashed line of mixed particles, and the solid line of insoluble particles, respectively. POM: primary organic matter, SOA: secondary organic aerosols, BC: black carbon. DU: desert dust, VA: volcanic ash. Upper panel adopted from Kaiser et al. (2014). In the current work, insoluble mode contains volcanic ash only while soluble mode contains only SO 2− 4 and H2O.
where M 0,i and M 3,i describe the zeroth (number density) and third (mass concentration) moment of mode i, respectively.

205
The terms Ca, Co and N u refer to coagulation, condensation and nucleation, respectively. The terms Ca m,ii and Ca m,ij are intra and inter-modal coagulation in the moment m, respectively. Nucleation is considered for the Aitken mode only.
Condensation and coagulation affect all modes except the giant mode. The nucleation, condensation and coagulation terms are calculated following Riemer et al. (2003) and Vogel et al. (2009). Furthermore, ISORROPIA II model is used to calculate the gas-aerosol partitioning according to thermodynamic equilibrium (Fountoukis and Nenes, 2007).

210
Shifting between modes is performed using two mechanisms. The first mechanism is activated when a threshold diameter is exceeded. Then, a shift to a corresponding mode with larger median diameter is performed. The second mechanism shifts mass and number concentration from insoluble modes to mixed modes if a mass threshold of soluble coating on insoluble particles (currently 5 %) is exceeded (Weingartner et al., 1997).

215
The RRTM requires the mass extinction coefficient k e , single scattering albedo ω, and asymmetry parameter g in 30 wavelength bands to account for the radiative effect of aerosols (Gasch et al., 2017). In this connection, k e can be interpreted as the extinction cross-section per aerosol mass in the units m 2 kg −1 .The wavelength bands range between 0.2 and 100 µm. The calculation of the optical properties is based on the wavelength-dependent refractive indices of volcanic ash (Walter, 2019), water, and sulfuric acid (Gordon et al., 2017).

220
No study so far has treated volcanic ash as a core in an internal mixture. It is suggested, but not proven, that most volcanic ash particles are coated to some degree (Bagnato et al., 2013;Hoshyaripour et al., 2015). Therefore, the core-shell treatment is physically more realistic than the external-mixture treatment even though the reality lies between the externally mixed and core treatments (Jacobson, 2000;Riemer et al., 2019). Hence, this study deploys both externally mixed (in the soluble and insoluble modes) and internally mixed (in the mixed mode) treatments. For the mixed mode, we use the core-shell model in which the 225 core and shell consist of well-mixed volcanic ash and H 2 O-H 2 SO 4 solution, respectively. To calculate the optical properties, the Mie code for coated spheres is used which has been developed by Mätzler (2002) and Bond et al. (2006) based on Bohren and Huffman (1983). Based on the ICON-ART simulations the shell fraction (increased diameter due to coating) is assumed to be 0.2 with 50 % H 2 O-H 2 SO 4 solution. The volume-average mixing rule is used to compute the complex refractive index of each layer, which then serves as input for the core-shell calculation. The results of the Mie calculations for the ash-containing modes are shown in Fig. 2. It can be seen that the mixed modes (coated ash) have higher k e and ω in the visible range than the insoluble modes (uncoated ash). This is caused by the H 2 O-H 2 SO 4 coating which is a strong scatterer. Particles with a strongly absorbing core coated by a weak absorber generally absorb more sunlight than an external mixture of the same components, which is caused by the increase of the core cross section due to coating (Jacobson, 2000). This is not the case for volcanic ash as it is not a strong absorber compared to soot particles. This 235 can be seen in the imaginary part of refractive indices, i.e., absorbing part, at 500 nm that are 0.00092 and 0.74 for volcanic ash and soot, respectively. The Mie theory assumes that the particles have spherical shapes. In reality, volcanic ash particles are exclusively nonspherical particles (Bagheri and Bonadonna, 2016). Therefore, their optical properties may be better represented by spheroids, ellipsoids or even more complex structures (Gasteiger et al., 2011;Vogel et al., 2017). However, the liquid coating can lead 240 to spherical particle surfaces, which justifies the assumption of the particle sphericity in the mixed mode. For consistency reasons, the sphericity assumption is also applied to the insoluble mode that contains uncoated ash particles. Implementing The plume height estimate is based on the MODIS and VIIRS data shown in Fig. 3. The dedicated ash algorithm (lower panel) is much more restrictive than the standard cloud-top height algorithm (upper panel), but produces similar heights where it is applied. In general, both of these brightness temperature-based products indicate maximum plume heights in the 12-255 12.6 km range for the time period 7-9 h after the eruption. The estimated height uncertainty is ∼ 1.5 km. Based on this plume height estimate and also other studies (Sennet, 2019), the Raikoke eruption emits ash and SO 2 in our simulations at a constant eruption rate between 8 and 14 km above sea level.
The eruption rate of SO 2 is derived from measurements of the total emitted SO 2 mass. According to the TROPOMI (Sect. 2.1.1) and AHI data (Sect. 2.1.2), in our simulation 1.5 × 10 9 kg of SO 2 is emitted over the eruption period. To estimate the 260 total mass eruption rate of volcanic ash, several 1D plume simulations using Plumeria (Mastin, 2007) and FPlume (Folch et al., 2016) are conducted assuming the following parameter ranges: plume height 12-14 km, vent diameter 90-110 m, exit velocity 100-120 m s −1 , exit temperature 900-1100 • C, and exit gas mass fraction 3 %. For this purpose, atmospheric profiles are obtained from ERA-Interim (Dee et al., 2011) and introduced in the 1D models as wind and no-wind atmospheres. By this method, the key sources of uncertainty are considered in the estimation of mass eruption rate. The results are in the range of 265 1.45-9.95 × 10 6 kg s −1 . Taking the mean value 5.7 × 10 6 kg s −1 suggests that about 190 × 10 9 kg tephra is emitted within 9 hours. Assuming that 1 % of the erupted mass is very fine ash with d < 30 µm (relevant for long range transport) Gouhier et al., 2019), we estimate that 1.9 × 10 9 kg very fine ash is injected into the atmosphere during the eruption. The estimates by the 1D models are in agreement with AHI data (Sect. 2.1.2).
The estimated 1.9 × 10 9 kg of very fine ash are used in the ICON-ART simulations and distributed equally between accu-270 mulation, coarse, and giant modes. The number concentration of the log-normal distribution is calculated based on the median diameter d e and standard deviation σ e of the emitted particle distribution. Table 1 lists details about these emitted particle size distributions.
We study the effect of aerosol dynamic processes and the radiative effect of internally mixed particles on the volcanic plume dispersion with the help of three different simulation scenarios summarized in Table 2  aerosol formation and particle aging are switched off. However, the volcanic ash still interacts with solar and thermal radiation.
The third scenario (AERODYN-no_rad) considers the effects of aerosol aging without any radiative feedback of these particles.
The two scenarios with AERODYN treat SO 2 as a chemical substance which can be oxidized. The chemical reaction scheme 280 is a simplified OH-chemistry scheme that has been implemented into ICON-ART by Weimer et al. (2017). The no_AERODYN scenario treats SO 2 as a passive tracer without any gas phase chemistry.   Fig. 4 (c) and (e). For the first day after the eruption, the aerosol-radiation interaction has no significant influence on the volcanic ash mass loading. On 23 June the averaged AHI measurements show a more fragmentary ash distribution in Fig. 4 (b). This might be a result of volcanic cloud dilution in 295 combination with deficiencies in the volcanic ash measurement of opaque regions. Most of the ash is measured between 50-55 • N and around 180 • E. The simulation results in Fig. 4 (d) and (e) support the assumption of the diluted volcanic cloud, as the mass loading only shows values smaller than 2 g m 2 . For both simulated scenarios, the overall structure of the volcanic cloud is similar. However, differences prevail in location and absolute values of maximum mass loading. These differences are due to radiative effects which are addressed in more detail in Sect. 3.3. Compared to these two simulations, the averaged AHI 300 measurements ( Fig. 4 (b)) shows slightly higher values for the maximum ash mass loading. This could be an artifact of the averaging approach, as it favors single pixels with high values for patchy retrievals. Each of these three graphs is a composite of several satellite orbits, chosen from a batch of 14 consecutive orbits (approximately 24 h coverage). Those orbits that directly detect the volcanic cloud in Fig. 5 (a) intersected with the area of interest (see Sect. of the SO 2 mass loading agrees well between model results and observations. This is especially true for the two earlier dates 310 when the modelled atmospheric state can be assumed to be closer to reality than for later dates. But also the model result 3.5 days after its initialization in Fig. 5 (f) shows very good agreement with the TROPOMI measurement in (c). A main difference between satellite retrieval and model result is the location of the maximum SO 2 mass loading. Although the magnitude of the maximum SO 2 mass loading is in good agreement, in the model results its location appears further downstream compared to the satellite measurement. One reason could be the different time of measurement and model result. However, a greater 315 influence can be expected by uncertainties of the emission profile parametrization and of the simulated wind velocities. In case more SO 2 is emitted in altitudes with higher wind speeds in the model, it will be transported faster. The same applies for the case that in some altitudes wind speeds in the model are slightly higher than they are in reality. Furthermore, the TROPOMI measurements can also be erroneous. The TROPOMI sensor might not capture all of the SO 2 due to deficiencies of the measurement technique in opaque regions. Assumptions about a vertical SO 2 profile made for the retrieval can also result 320 in incorrect SO 2 mass loadings. a similar cloud top height which will be addressed in more detail in Sect. 3.3. Also the height of the volcanic cloud measured by CALIOP on 23 June 2019, agrees well with the model result. This will be addressed in more detail in the following section.

Effect of aerosol dynamics
So far we only discussed the ICON-ART model result of the AERODYN-rad scenario. In this section, we compare it with the no_AERODYN-rad scenario to study the influence of secondary aerosol formation and particle aging on volcanic aerosol For example, the peak in the simulated attenuated backscatter (Fig. 6 (c)) at around 44 • N up to 3 km is also present in the 345 CALIOP signal at a comparable order of magnitude. This suggests that the elevated CALIOP signal in this region is due to volcanic aerosols. Other features in the modeled attenuated backscatter, north of 51 • N, also collocate with structures in the CALIOP signal. This suggests that part of the elevated CALIOP signal in these regions is due to the volcanic aerosol cloud. It nicely shows the advantage of considering model results for the interpretation of satellite retrievals.
Comparing AERODYN-rad in Fig. 6 (c) with no_AERODYN-rad in Fig. 6 (d) shows the distinct effect of aerosol dynamics 350 on vertical distribution of the volcanic cloud. No_AERODYN-rad catches the main feature between 49 • N and 51 • N at a height up to 17 km. However, the volcanic aerosol layer extends significantly further north, up to 54 • N. This is in contrast to the CALIOP signal in Fig. 6 (b). Also the smaller patterns in lower altitudes and higher latitudes are missing in the no_AERODYNrad scenario. The same applies for the feature at around 44 • N and 3 km height. Without aerosol dynamics, most of the aerosol stays at one height level, whereas with aerosol dynamics, the particles get also mixed down to lower altitudes. Coagulation 355 of particles and condensation of sulfate and water onto existing particles increases the aerosol mass. Hence, these particles sediment faster and therefore, are removed from the atmosphere more efficiently.
To further investigate the effect of aerosol dynamics on the residence time of very fine ash, we examine the temporal variation of ash concentration in the atmosphere. This is illustrated in Fig. 7. The graph shows how the normalized total ash mass m ash two displayed simulation scenarios (panels (c) to (f)). However, the differences in the mass loading patterns can be explained by radiative effects. In order to investigate the influence of aerosol-radiation interaction on volcanic plume dispersion in more detail, we look at the maximum height that the volcanic cloud reaches over the course of time. A volcanic cloud that is lifted up in the atmosphere has a longer lifetime. Hence, it can be transported over longer distances, remains a hazard for aircraft over a longer period of 385 time, and has longer lasting climatic effects. Additionally, the height of the volcanic cloud in the atmosphere also influences its transport, as wind speed and direction can differ between height levels. Figure 8  Prata, A., O Brien, D., Rose, W., and Self, S.: Global, long-term sulphur dioxide measurements from TOVS data: A new tool for studying explosive volcanism and climate, GEOPHYSICAL MONOGRAPH-AMERICAN GEOPHYSICAL UNION, 139, 75-92, 2003.