Modeling volcanic ash aggregation processes and related impacts on the April–May 2010 eruptions of Eyjafjallajökull volcano with WRF-Chem

Volcanic eruptions eject ash and gases into the atmosphere that can contribute to significant hazards to aviation, public and environment health, and the economy. Several volcanic ash transport and dispersion (VATD) models are in use to simulate volcanic ash transport operationally, but none include a treatment of volcanic ash aggregation processes. Volcanic ash aggregation can greatly reduce the atmospheric budget, dispersion and lifetime of ash particles, and therefore its impacts. To enhance our understanding and modeling capabilities of the ash aggregation process, a volcanic ash aggregation scheme was integrated into the Weather Research Forecasting with online Chemistry (WRF-Chem) model. Aggregation rates and ash mass loss in this modified code are calculated in line with the meteorological conditions, providing a fully coupled treatment of aggregation processes. The updated-model results were compared to field measurements of tephra fallout and in situ airborne measurements of ash particles from the April–May 2010 eruptions of Eyjafjallajökull volcano, Iceland. WRF-Chem, coupled with the newly added aggregation code, modeled ash clouds that agreed spatially and temporally with these in situ and field measurements. A sensitivity study provided insights into the mechanics of the aggregation code by analyzing each aggregation process (collision kernel) independently, as well as by varying the fractal dimension of the newly formed aggregates. In addition, the airborne lifetime (e-folding) of total domain ash mass was analyzed for a range of fractal dimensions, and a maximum reduction of 79.5 % of the airborne ash lifetime was noted.

1 Introduction

Volcanic eruptions inject gases and ash particles of various sizes into the atmosphere, posing hazards to life, infrastructure, and aviation (Miller and Casadevall, 2000). Volcanic emissions can alter the composition of the atmosphere and affect the Earth's radiation budget and climate (Angell, 1993; Cole-Dai, 2010; Thordarson and Self, 2003). The environmental and economic impacts of past and recent eruptions have spurred increased interest in the inclusion of volcanic ash into numerical weather prediction (NWP) models (Folch et al., 2009, 2016; Lin et al., 2012; Stuefer et al., 2013). Today, forecasters and scientists utilize volcanic ash transport and dispersion (VATD) models for ash hazard mitigation; the development, calibration, and validation of remote sensing tools; and the study of ash physics. Numerical models have been developed to better describe the initial plume characteristics of eruptions, such as the plume height, water content, particle size distribution, and plume shape. A current limitation of most VATD models is their ability to capture volcanic ash aggregation.

Volcanic ash aggregation is important for many reasons. Aggregation affects the atmospheric lifetime of ash, the distance ash is transported from the eruption source, the size and type of tephra observed on the ground, and the duration ash poses a threat to aircraft (Brown et al., 2012; Casadevall, 1994; Rose and Durant, 2011). Aggregation has been observed in several well-studied volcanic eruptions such as those of Mount St. Helens (Washington), Mount Redoubt (Alaska), and Eyjafjallajökull (Iceland). Additionally, aggregation occurs in both proximal (<15 km from the plume edge) and distal ash clouds (Bonadonna et al., 2011; Brown et al., 2012; Carey and Sigurdsson, 1982; Rose and Durant, 2009, 2011; Taddeucci et al., 2011; Wallace et al., 2013).

Proximal volcanic ash aggregates form more rapidly than distal aggregates for a number of reasons. For example, ice and liquid water enhance the sticking of particles and thus increase the rate of aggregation (Brown et al., 2012; Rose and Durant, 2011). This process can occur in a hail-like process with a cycle of freezing and thawing leading to enhanced aggregation (Van Eaton et al., 2015). In addition, the higher concentration of ash in the proximal plume increases the number of collisions.

Water-enhanced aggregation in the proximal plume has been observed in a number of eruptions. Field observations of tephra from the 18 May 1980 eruption of Mount St. Helens detail the formation of large volcanic aggregates (up to 1 mm) closely correlated with the presence of rain, snow, and hail (Waitt et al., 1981). Gilbert and Lane (1994) note that aggregation rates were enhanced by high proximal water vapor concentrations during the eruptions of Sakurajima volcano in the 1990s, and the majority of this water-enhanced aggregation occurred proximally, within the first minutes of the eruption. In addition, studies of the 2009 eruption of Mount Redoubt in Alaska show definitive evidence for aggregation-enhanced sedimentation in the proximal plume (Van Eaton et al., 2015; Wallace et al., 2013). Van Eaton et al. (2015) conclude that the effects of aggregation in the Redoubt eruption resulted in over 95 % of fine ash mass deposited to the ground as aggregates.

Distal aggregation usually occurs at a slower rate than proximal aggregation as the plume ages and diffuses (Rose and Durant, 2009, 2011). Despite a slower rate of aggregation the majority of distal fine ash settles to the ground as larger aggregates (Brown et al., 2012; Carey and Sigurdsson, 1982; Rose and Durant, 2011; Wallace et al., 2013). Both coarse and fine ash particles are known to aggregate in distal clouds by forming dry clusters due to electrostatic attraction, or as liquid or frozen water particles (Brown et al., 2012; Rose and Durant, 2011). Distal aggregate formation has been observed from eruptions such as Mount Etna, Italy, in 1971, Mount St. Helens, US, in 1980, and Mount Redoubt, US, in 1990 (Booth and Walker, 1973; Sorem, 1982; Sparks et al., 1997). For many eruptions, electrostatic aggregation of fine ash is expected to be responsible for the bimodal distribution of volcanic ash fallout (Carey and Sigurdsson, 1982; Cornell et al., 1983; James et al., 2003).

Recently, aggregation processes were observed to play an integral role in the dispersion of the plume generated from the April and May 2010 eruptions of Eyjafjallajökull volcano, Iceland. In situ measurements of ash particle fall velocities using high-speed photography observed aggregation-enhanced sedimentation that increased fallout rates by a factor of 10 (Taddeucci et al., 2011). The effect of ash aggregation caused a significant quantity of additional ashfall across Iceland, rather than being transported further. Ash aggregation overall clearly reduced the atmospheric residency time of the Eyjafjallajökull ash plume (Gudmundsson et al., 2012). In addition, aggregation was observed to cause enhanced fallout over parts of mainland Europe and the United Kingdom (Stevenson et al., 2012).

Aggregation processes not only affect the lifetime of volcanic ash, but also the makeup of volcanic ash cloud particle size distributions (PSDs) which may complicate modeling and remote sensing efforts (Brown et al., 2012; Rose and Durant, 2011). For example, volcanic ash remote sensing algorithms require information regarding particle sizes and extinction coefficients (Stohl et al., 2011; Wallace et al., 2013). Remote sensing methods are also used to estimate eruption parameters and PSDs via extinction coefficients using inverse modeling (Kristiansen et al., 2012; Stohl et al., 2011). Additionally, volcanic PSDs are also important for the study of radiative properties of volcanic ash and their effects on the atmosphere (Hirtl et al., 2019; Young et al., 2012).

The effects imposed on volcanic ash clouds by aggregation processes necessitates their parameterization in volcanic ash transport and dispersion (VATD) models. Despite this, only a few of the existing VATD models capture aggregation processes. For example, a volcanic ash aggregation parameterization scheme has been implemented within the FALL3D model (Folch et al., 2009). In an operational setting, FALL3D runs by ingesting offline meteorological fields from gridded atmospheric models, such as the Weather Research Forecasting (WRF) model, and then calculating volcanic ash advection and sedimentation during the parent model output time step. Another method of capturing volcanic ash aggregation is to initialize VATD models with PSDs that account for volcanic aggregation in the eruptive column by using initial plume models. FPLUME, a one-dimensional (1D) plume model based on buoyant plume theory, constructs initial plume characteristics that account for ash aggregation (Folch et al., 2016). In this case, the 1D plume model develops an initial PSD at the source that accounts for aggregation processes and then keeps this PSD invariant during further plume transport.

In an effort to study and predict volcanic ash aggregation effects using a fully coupled modeling system, where the fate of the airborne ash particles is coupled to the atmospheric environment, a volcanic ash aggregation scheme was incorporated into the Weather Research Forecasting with Chemistry (WRF-Chem) model (Grell et al., 2005). This coupled system requires no temporal or spatial interpolations as it calculates interactions between the meteorology and ash at each modeling time step (on the order of seconds). While many dispersion models require less computing power than WRF, a number of them require a mesoscale model, like WRF, to generate regional, gridded meteorological fields for their initialization. As an example, FALL3D is typically initialized with a WRF model run that is executed prior to the dispersion model. Modeling particle dispersion with WRF-Chem is, therefore, as computationally feasible as running these models since in many cases a mesoscale gridded model must be run for their initialization. In addition, volcanic ash aircraft hazard mitigation typically focuses on limiting commercial aircraft to ash concentration thresholds (Casadevall, 1994). WRF-Chem solves the advection equations such that ash concentration is tracked over time. This ability to track volcanic ash mass, rather than particle number as is done in many VATD particle dispersion models, augments current VATD model guidance and offers another tool to constrain atmospheric ash loading.

The following sections of this paper detail the inclusion of a computationally feasible volcanic ash aggregation scheme into the WRF-Chem model and the impacts of these modifications on model output. The following “Aggregation parameterization and implementation” section (Sect. 2) details the background and incorporation of a mathematical scheme that is physically descriptive of aggregation processes into WRF-Chem, as well as the development of a new methodology for selecting aggregation sticking efficiencies that depend on relative humidity. This newly implemented code is then applied to the April and May 2010 eruptions of Eyjafjallajökull, as well as to a controlled sensitivity study using a single eruption. The setup of these two cases is discussed in Sect. 3 “Methods”, with remarks on the model output in Sect. 4 “Results”. Concluding remarks are then provided in the final Sect. 5 “Conclusions”.

2 Aggregation parameterization and implementation

Smoluchowski (1917) developed the original analytical theory of the process of coagulation of colloid particles based upon Richard Zsigmondy's experiments with gold solutions. The Smoluchowski coagulation equation (Eq. 1) is an integrodifferential population balance equation that describes the evolution of particle number density, n_v(v), in time, t, as primary particles of one volume, v, collide and stick together with particles of different volumes, v^′, to form aggregates (Smoluchowski, 1917). It is physically descriptive of the aggregation process.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial n_{\mathrm{v}}\left(v\right)}{\partial t}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{v}{\int }}K(v-v^{\prime },v)n_{\mathrm{v}}(v-v^{\prime })n_{\mathrm{v}}\left(v^{\prime }\right)d{v}^{\prime }-n_{\mathrm{v}}\left(v\right)\\ & \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}K(v,v^{\prime })n_{\mathrm{v}}\left(v^{\prime }\right)d{v}^{\prime }\end{array}\end{array}$$

Equation (1) describes the number of aggregates of volume v formed, n_v, per unit time, t, on the left, and the loss of primary particles between volumes v and v^′ on the right as particles aggregate based on the collision frequency of the particles. Frequency is weighted by the coagulation kernel, K, which is the product of the collision kernel, A, and a sticking efficiency, α; thus, K=Aα.

Volcanic ash may undergo various processes that result in collisions, such as Brownian motion, differential sedimentation, and fluid shear, and as a result there are many formulations of the coagulation kernel, K (Jacobson, 2005). For example, collisions due to Brownian interactions (A_B) occur randomly during diffusion and are temperature dependent. As temperature increases, the diffusion rate increases, thus increasing their chances of interacting with other particles. Particle collisions due to shear (A_S) occur when ash moving in different horizontal directions collides due to changes in laminar flow. This kernel therefore depends on wind speed and direction. Lastly, differential sedimentation (A_DS) captures particle interactions due to the different fall velocities of different sized particles. The rate at which particles settle is dependent on their size and therefore the differential sedimentation kernel depends on the difference in size between particles. As larger particles fall, they have a greater chance of encountering smaller, slower-moving particles on their descent. In summary, the collision kernels A_B, A_S, and A_DS represent the rate at which ash particles collide based on Brownian motion, fluid shear, and differential sedimentation, respectively. Each kernel depends directly on the number concentration and size distribution of ash particles, and each depends highly on its own set of parameters.

Table 1 Derived coagulation kernel equations used in the calculation of ΔN.

While physically descriptive of the aggregation process, the Smoluchowski equation itself, in addition to the equations governing the coagulation kernel, K, is prohibitively computationally expensive to solve explicitly, even with simple boundary conditions. Advances in simplifying the equation for use in computational volcanic ash modeling resulted in large part from work by Dekkers and Friedlander (2002) and Costa et al. (2010) by assuming a time-independent aggregate size distribution and fractal geometry of volcanic ash aggregates, respectively. Assuming a fractal aggregate geometry greatly simplifies the equations describing the coagulation kernels (A_B, A_S, and A_DS) by establishing a particle size–volume fractal relationship, described by a fractal dimension factor, ξ. In addition, an assumption of fractal geometry allows n_v in Eq. (1) to be described in terms of the total number of particles in a computational space, n_tot, forming aggregates of a certain fractal dimension, D_f, based on a generally accepted fractal relationship (Jullien and Botet, 1987; Lee and Kramer, 2004). The simplified Smoluchowski equation described by Costa et al. (2010) results in a calculation of $\frac{\partial n_{\mathrm{v}}\left(v\right)}{\partial t}$, from Eq. (1), that is much more computationally feasible (Eq. 2):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{n}_{\mathrm{tot}}=\mathit{\alpha }\left(A_{\mathrm{B}}{n}_{\mathrm{tot}}^{\mathrm{2}}+A_{\mathrm{s}}{\mathrm{\Phi }}^{\frac{\mathrm{3}}{D_{\mathrm{f}}}}{n}_{\mathrm{tot}}^{\mathrm{2}-\frac{\mathrm{3}}{D_{\mathrm{f}}}}+A_{\mathrm{DS}}{\mathrm{\Phi }}^{\frac{\mathrm{4}}{D_{\mathrm{f}}}}{n}_{\mathrm{tot}}^{\mathrm{2}-\frac{\mathrm{4}}{D_{\mathrm{f}}}}\right)\mathrm{\Delta }t.\end{array}$$

Here, Δn_tot represents the total number of particles per unit volume lost to aggregation. The equation relies on the solid volume fraction of the aggregates, Φ (Folch et al., 2016), the number densities of the bins, n_tot, and the fractal dimension of the fine ash particles, D_f (Costa et al., 2010). Equations describing the collision kernel, A, were also simplified using a fractal representation of ash geometry and were reduced to Eq. (3) through Eq. (5), shown in Table 1.

New code capable of calculating Eqs. (2)–(5) was developed in this study and integrated into the Fortran 90 module “module_vash_settling.F” file, located in the “chem” subdirectory of the WRF main directory, which is available to download from the WRF home page: https://www2.mmm.ucar.edu/wrf/users/downloads.html (last access: August 2019). Modified code is available upon request. See the “Code availability” section for details.

Most of the source variables necessary to solve Eqs. (2) to (5) are available in WRF-Chem by selecting the appropriate aerosol and chemistry packages. For example, chemistry option (chem_opt) 402 (WRF-Chem User Guide 3.9; Peckham et al., 2018) includes chemistry and humidity variables provided by the Regional Deposition Acid Model Version 2 (RADM2) (Stockwell et al., 1990) and the Goddard Chemistry Aerosol Radiation and Transport (GOCART) models (Chin et al., 2000), as well as the inclusion of volcanic sulfur dioxide (SO₂) and 10 volcanic ash particle size bins (Stuefer et al., 2013). Three variables required by Eqs. (2) to (5), the sticking efficiency, α, fractal dimension, D_f, and fractal dimension factor, ξ, are not, however, included in WRF-Chem and therefore must be calculated or assumed.

The fractal dimension, D_f, relates the number of primary particles N in an aggregate to the size of the aggregate, R, such that N scales proportionally as $N\propto R^{D_{\mathrm{f}}}$. For example, as D_f approaches 3, primary particles in the aggregate use up more and more space such that D_f=3.0 would indicate a solid, filled aggregate. A lack of experimental data adds a degree of uncertainty when selecting the fractal dimension; however previous analysis studies of aggregates selected after the eruptive events from Mount St. Helens and Mount Spurr suggested a dimension D_f=2.99. This favorable fractal dimension resulted from a regression analysis between model output and observed deposits (Folch et al., 2010). The fidelity of confidence in the choice of the fractal dimension is hindered by the fact that it does not necessarily, by its definition, remain constant within a plume.

The fractal dimension factor, ξ, used to simplify the coagulation kernel equations relates the fractal dimension, D_f, to the diameters and volumes of the primary particles in the aggregates. This relationship is given in Eq. (3):

$$\begin{array}{}\text{(6)}& d_i=\mathit{\xi }{v}_i^{\frac{\mathrm{1}}{D_{\mathrm{f}}}}.\end{array}$$

Here, d_i and v_i are the diameter and volume of the primary particles forming an aggregate. Costa et al. (2010), Dekkers and Friedlander (2002), and Folch et al. (2010) indicate that a fractal dimension on the order of 0.6 to 1 is sufficient for describing the geometry of volcanic ash particles and aggregates. As done in Costa et al. (2010), a unity fractal dimension factor is utilized in this study.

The sticking efficiency coefficient, α, relies heavily on the concentration of water vapor and ice (Costa et al., 2010). In order to formulate an appropriate estimate for the sticking efficiency coefficient, a new parameterization was incorporated into the WRF-Chem emissions driver that includes volcanic water vapor emissions that are specified by the user. This code adds these emissions to the ambient water vapor mass within the model environment. Van Eaton et al. (2012) demonstrated that the sticking efficiency of volcanic ash particles follows exponential curves. Using these fitted curves, the sticking efficiency coefficient, α, between two particles i and j may be calculated using a fitting coefficient, S. This coefficient varies with water vapor concentration, [H₂O], and the radius of the colliding particles, r. A lookup table was added to select sticking coefficients based on this work by utilizing the water vapor content of the model cell and the particle size (Eq. 4 and Table 2). Importantly, this equation is computationally inexpensive to solve. Although electrostatic interactions are significant enough to cause aggregation of particles, they are most likely insignificant when compared to aggregation in the presence of water (James et al., 2003; Schumacher and Schmincke, 1995). Since the modeled background water rarely approaches 0 % relative humidity, dry interactions are not parameterized in this study.

$$\begin{array}{}\text{(7)}& \mathit{\alpha }\left(ij\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]\right)=e^{-\mathrm{Sr}}\end{array}$$

The four aggregation equations (Eqs. 2 to 5) are solved for volcanic ash bins 2 to 10 (Table 3) at every time step, for every model grid cell, and account for interaction of particles between the different bins by using the total mass to calculate the available number of primary particles available for aggregation. Large particles, greater than 1 mm in diameter, are included in WRF-Chem volcanic ash bin 1, which has been designated as the “aggregate” bin. All aggregates generated by the code are moved to bin 1, and their corresponding masses are subtracted from bins 2–10. The large particles (in bin 1) assume high fall velocities and contribute to ash fallout within periods of minutes (Rose and Durant, 2011). All volcanic ash removed from the model domain is stored in the ASH_FALL variable, allowing the analysis of fallout mass and location.

Table 2 Ash aggregation coefficients based on liquid water content, w∕w, as described in Van Eaton et al. (2012). The weight percent of water (w∕w) is calculated as mass of water divided by mass of the atmosphere.

Table 3 Distribution of volcanic ash in model domain among 10 size bins corresponding to the S2 size distribution as given in Mastin et al. (2009). The percentages of mass per bin are specified in the volc_d01.asc name list and may be given any value between 0 and 100.

3 Case study and methods

The 2010 eruption of Eyjafjallajökull volcano in Iceland has been selected to test the modified WRF-Chem modeling experiment. Eyjafjallajökull erupted in April and May 2010, dispersing ash over Europe that caused numerous flight delays over the course of weeks and a resulting loss of revenue to airlines in the billions of dollars (Harris et al., 2012). Due to Eyjafjallajökull's location and the availability of observational resources, it became one of the most studied and well-documented eruptions in history, providing numerous sources of data regarding the plumes' characteristics. The German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt, DLR) took several in situ measurements of Eyjafjallajökull's ash clouds over the course of the 2 months of eruptions by flying its Falcon aircraft into forecasted plume locations. Three of these flights are used for analysis in this study from 19 April and 16 and 17 May. The flight paths corresponding to these flights are depicted using colored lines in Fig. 1. During the flights, Schumann et al. (2011) recorded particle number concentrations using a Grimm SKY-OPC 1.129 optical particle counter and a Particle Measuring Systems, Inc. (PMS) forward scattering spectrometer probe (FSSP), observing a range of particles from 0.25 and 24 µm. In addition, upper and lower mass concentration estimates were calculated using the minimum and maximum imaginary component of the refractive index, of which the FSSP was particularly sensitive. For the flight of 17 May, a medium estimate of mass concentration was calculated. From these studies, information on particle number, mass concentration, plume heights, and gas composition are available, providing one of the best in situ datasets available to study distal and proximal volcanic emissions (Schumann et al., 2011). In addition to these in situ data, Doppler radar measurements of the eruptive column and ground air sampling measurements were conducted by many groups to establish descriptive and accurate eruption source parameters (Arason et al., 2011; Devenish et al., 2012a, b; Stevenson et al., 2012). Observations of volcanic tephra fallout are also available and provide important insights into the PSD and transport of the distal Eyjafjallajökull ash clouds (Gudmundsson et al., 2012; Stevenson et al., 2012). In addition, volcanic ash aggregation was directly observed via high-speed photography near the vent, lending proof that particle aggregation occurred in the plumes Eyjafjallajökull produced (Taddeucci et al., 2011).

Figure 1 WRF-Chem model domain used for simulations in Lambert conformal projection with true latitude and longitude and center at 50^∘ N, 0^∘ E/W. Location of Eyjafjallajökull (63.62^∘ N, 19.61^∘ W) marked with red dot. DLR Falcon flight paths for flights on 19 April (red), 16 May (green), and 17 May (blue) are shown with colored lines.

3.1 Eyjafjallajökull model domain setup

The newly implemented aggregation code was applied to the April and May 2010 eruptions of Eyjafjallajökull. Additionally, sensitivity studies were conducted using a hypothetic single eruption of Eyjafjallajökull on 5 May 2010. In all studies, the model domain was centered at 50^∘ N, 0^∘ W, offsetting the Eyjafjallajökull vent (63.62^∘ N, 19.61^∘ E) to the northwest of the domain to account for the predominant southwest trajectory of the ash clouds. The model was set up with a resolution of 10 km² per grid cell and a total of 500×500 horizontal grid cells. This resolution was chosen as a compromise between the long timescale of the model study (on the order of months), the large spatial extent of the model domain required to study the physics of the distal plume, and the amount of available computational time. The domain is shown in Fig. 1 with Eyjafjallajökull marked in red. The model included 48 vertical pressure levels with the top level of the model set to 2000 Pa. The integration time step of the dynamics and chemical fields was set to 30 s.

Meteorological fields were obtained from the National Center for Environmental Prediction Final Global Operational Analysis (NCEP FNL) datasets, ds083.2, accessed through the National Center for Atmospheric Research Data Archive (NCAR, 2000). These datasets represent the final analysis of historical Global Forecast System (GFS) model output. Ingest was conducted similar to Hirtl et al. (2019), using a 9 d spinup time before the first eruption on 14 April and with meteorological initializations every 48 h. Each reinitialization of the meteorological fields required the model to idle for varying periods due to competing jobs in the supercomputer's scheduler queue. The 48 h reinitialization of the meteorological fields balanced the need for sufficient synoptic-scale time coverage and the extensive size (order of months) of the time domain. The WRF-Chem volcanic package was enabled with chemistry option 402, which includes 10 particle sizes of volcanic ash (Stuefer et al., 2013). These particle sizes are shown in Table 3. The Yonsei University Planetary Boundary Layer (YSU PBL) scheme and the Noah Land Surface Model (LSM) were included for PBL and near-ground physics (Chen and Dudhia, 2001; Hong et al., 2006).

Water was added to the model domain by multiplying the water content of Eyjafjallajökull's magma, 1.8 % (Keiding and Sigmarsson, 2012), by the total erupted mass of 400 Tg for fine and coarse ash estimated by Taddeucci et al. (2011). This 1.8 % multiplier produces water vapor emissions that agree with constraints constructed by comparing H₂O∕SO₂ emission ratios using values from Allard et al. (2011), yielding a ratio of 458 mol mol⁻¹, and SO₂ emission rates from two remote sensing studies by Boichu et al. (2013) and Thomas and Prata (2011). The code was modified to read in volcanic water vapor emissions rates into WRF-Chem as a callable Fortran module.

In addition, Hirtl et al. (2019) noted that the model topography of Eyjafjallajökull is smoothed at the 10 km² model spatial resolution, resulting in a vent height 400 m lower than the actual height of 1000 m. A 400 m height offset was applied to correct this.

3.2 Sensitivity study model setup

Multiple sensitivity studies were conducted in order to assess (1) the overall change in mass due to aggregation, (2) the effects of different fractal dimensions, D_f, on the aggregation rate, (3) the contribution of each collision kernel, A_B, A_S, and A_DS, to the decrease in domain ash mass, and (4) the effect of adding coupled water vapor emissions to the model domain on the aggregation rate. A list of these sensitivity studies, including the parameters varied and the analysis approach used, is presented in Table 4. These sensitivity studies were conducted on a smaller time slice of the parent domain, using a 9 h eruptive event on 5 May 2019, initialized at 00:00 Z with a rate of 4×10⁶ kg s⁻¹, which corresponds to an average value of Eyjafjallajökull's largest eruptions. A 72 h spinup time was included prior to the eruption initialization to allow the meteorological fields to stabilize and was then run for 6 d, ending at 00:00 Z on 11 May. The smaller model time domain allowed for new meteorological fields to be reinitialized every 24 h, as opposed to 48 h in the longer timescale study. Each volcanic ash bin was populated with 10 % of the total erupted mass in order to simplify output analysis.

Table 4 Details of the model sensitivity studies discussed, including parameters varied and the analysis methodology.

In order to assess how the aggregation code affects model output, WRF-Chem was run with and without the aggregation code enabled. Due to a lack of experimental data, a choice of fractal dimension, D_f, is difficult. Therefore, the fractal dimension, D_f, was varied to measure its effects on the overall aggregation rate. The span of fractal dimensions chosen ranges from D_f={2.5, 2.6, 2.7, 2.8, 2.9, 2.95, 2.98, 2.99, 3.0} and is based on studies by Costa et al. (2010) and a similar study of Mount St. Helens and Mount Spurr using Fall3D by Folch et al. (2010).

The contribution of each collision kernel, A_B, A_S, and A_DS to the total reduction in domain mass was also assessed by using the same domain and eruption parameters and enabling only one kernel at a time using a fractal dimension of 2.5 and 3.0. The total change in mass from each kernel was then divided by the total change in mass with all kernels enabled to find the percent contribution.

The impacts of the inclusion of water vapor on the aggregation rate were studied by running the code with and without the 1.8 % water vapor emissions included in the model domain. For the simulation run without water vapor emissions, only background water vapor from the FNL datasets were used.

3.3 Model setup for April and May 2010 eruptions of Eyjafjallajökull

WRF-Chem was also configured to simulate phase I (14–18 April 2010) and phase III (4–18 May 2010) eruptions of Eyjafjallajökull using the same model domain described above. Phase II eruptions were effusive rather than explosive and ejected tephra at much lower altitudes of 2 to 4 km a.s.l. (Gudmundsson et al., 2012) and were thus not included in this modeling case study.

Eruption source parameters (ESPs) for Eyjafjallajökull were adapted from Mastin et al. (2014) and Hirtl et al. (2019). Camera footage and C-band Doppler radar measurements were used to establish 3-hourly plume heights for the April and May 2010 eruptions (Arason et al., 2011; Mastin et al., 2009; Hirtl et al., 2019). These plume heights were used to calculate eruption rates based on the plume height–eruption rate relationship derived by Mastin et al. (2009). The total erupted mass was then scaled based on work by Gudmundsson et al. (2012) such that the total ash mass ejected over the eruptive phases agreed with the 170 Tg phase I estimate and 190 Tg phase III estimates for fine ash stated (Hirtl et al., 2019). The bimodal, silicic (S2) ESP particle size distribution (Table 3) was used to populate the 10 volcanic ash bins in the model (Mastin et al., 2009). The 3-hourly plume heights and eruption rates used in the study are presented in Fig. 2.

Figure 2 Three hourly plume heights (KM) above sea level (orange, km) and emitted mass (blue, Tg) used in the WRF-Chem modeling simulations (volc_d01.asc name list) for the eruption period 12 April until 18 May 2010. Values adapted from Hirtl et al. (2019) with dates as DD/MMM.

In this study, all aggregation collision kernels were enabled, and water vapor emissions as described previously were added to the model domain at each time step. As mentioned earlier, the choice of a fractal dimension is hindered by a lack of experimental data. Folch et al. (2010) conducted linear regression analysis of repeated model run comparisons to tephra fallout measurements from eruptions originating at Mount Spurr and Mount St. Helens. This study resulted in the use of a D_f=2.99 fractal dimension. Due to a lack of experimental data on the development of volcanic ash fractal dimensions, and the fact that aggregate fractal dimensions are not necessarily constant with time, D_f was set at the upper bound of 3.0, providing a maximum effect of particle aggregation.

4 Results

The newly implemented aggregation parameterization was first assessed with a sensitivity study of a singular eruptive event, and then by application to the entire phase I and phase III eruption periods.

4.1 Sensitivity study results

Varying the fractal dimension between 2.5 and 3.0 resulted in a range of aggregation rates. Figure 3 illustrates the change in domain mass from a single 9 h eruption on 5 May at 00:00 Z with a constant eruption rate of 4×10⁶ kg s⁻¹. As expected, higher values of D_f result in higher rates of aggregation with the largest jumps in the aggregation rate between D_f=3.0 and 2.8. The degree to which aggregation reduced the overall ash domain mass can be seen in the peak mass loadings at hour 9 in Fig. 3. Here, the peak domain mass using D_f=3.0 is 17.4 Tg. This is a 72 % reduction in peak mass compared to the non-aggregation-enabled run of 62.9 Tg. Lower values of D_f provide almost no change in the total domain mass. For example, D_f=2.5 results in a 0.7 % decrease in peak mass by about 0.5 Tg.

Figure 3 Change in total domain ash mass (teragrams) for a hypothetical eruption on 5 May, beginning at 00:00 Z and ending at 09:00 Z, for a range of fractal dimensions, D_f={3.0, 2.99, 2.98, 2.95, 2.9, 2.8, 2.7, 2.6, 2.5}. Constant eruption rate = 4×10⁶ kg s⁻¹.

To quantify the change in aggregation rate, volcanic ash lifetimes in terms of e-folding time were calculated. This analysis is presented in Fig. 4 and indicates a range of e-folding times from 72 h with no aggregation code enabled to 15 h with maximum aggregation considered (D_f=3.0). As the fractal dimension increases, the atmospheric lifetime of volcanic ash decreases due to the incorporation of more volcanic ash particles into each aggregate. When considering fractal dimensions 2.7 and lower, the total lifetime is reduced only slightly, less than 4 %. Larger decreases in lifetime become apparent with D_f=2.8 (10 % decrease) and jump thereafter to a maximum 79.5 % decrease at D_f=2.99 and D_f=3.0 (same decrease for both). Based on work by Folch et al. (2010, 2016), it is assumed that an optimal value of the fractal dimension likely lies near D_f=2.99, which corresponds to a 79.5 % difference in e-folding times. In terms of volcanic ash lifetime, on hourly timescales, there is no difference between D_f=3.0 and 2.99.

Figure 4 Volcanic ash e-folding time in hours for a hypothetical eruption on 5 May, beginning at 00:00 Z and ending at 09:00 Z, for a range of fractal dimensions, D_f={3.0, 2.99, 2.98, 2.95, 2.9, 2.8, 2.7, 2.6, 2.5}. Constant eruption rate = 4×10⁶ kg s⁻¹.

Figure 5 shows the extent to which each kernel contributed to the overall change in the model domain's ash mass by enabling each kernel independently. Two fractal dimensions were considered, D_f=2.5 and 3.0, and both affected each kernel's contribution to aggregation differently. The differential sedimentation kernel, A_DS, for example contributed to the majority of the change in domain mass over the course of the 96 h model run (≈99 %) when D_f was set to 3.0, but contributed only 5 % on average with D_f=2.5. The Brownian kernel became the major contributor to aggregation in the case of D_f=2.5, contributing to over 90 % of the aggregation. This agrees with parametric studies of varying fractal dimensions by Costa et al. (2010), who noted this trade between A_DS and A_B when considering fine ash particles (<63 µm). Overall, fluid shear interactions were the minor contributor to aggregation for both fractal dimensions. While its contribution to aggregation approaches that of A_DS for D_f=2.5, it is many orders of magnitude lower than A_B or A_DS for D_f=3.0.

Figure 5 Percentage of aggregation rate for each collision kernel (A_B: Brownian; A_S: shear; A_DS: differential sedimentation) when considering a hypothetical eruption on 5 May, beginning at 00:00 Z and ending at 09:00 Z, for two fractal dimensions, D_f={3.0, 2.5}. Constant eruption rate = 4×10⁶ kg s⁻¹.

Figure 6 illustrates the total domain mass for fine ash (bins 7–10) in Fig. 6a as well as the percentage of total domain mass in Fig. 6b, representing the PSD of the fine ash fraction. Figure 6 considers maximum aggregation with D_f=3.0. The bins with larger ash particles (1–6) were not included due to the rapid decrease in their domain mass as a result of their high settling velocities. Figure 6a depicts a decreased mass loading for each bin when aggregation is enabled, as well as a shorter lifetime, as expected. Figure 6b depicts a shift in the particle size distribution due to aggregation. The aggregation code results in less contribution from fine ash particles (bins 7–9), resulting in a shift of the PSD towards bin 10. Bin 10 in the aggregation-enabled code makes up an extra 10 % of the model domain mass upon reaching near-steady state at model hour 120. This is the result of the increased aggregation of the larger size particles since larger radii result in a larger probability cross section of collision and subsequent aggregate formation.

Figure 6 Total domain mass (a) and particle size distribution (b) of volcanic ash bins 7 to 10 when considering a hypothetical eruption on 5 May, beginning at 00:00 Z and ending at 09:00 Z, and a fractal dimension D_f=3.0. Constant eruption rate = 4×10⁶ kg s⁻¹.

Coupling water emissions resulted in a very small increase in aggregation rate, lowering the total domain mass on the order of megagrams per hour, much lower than the overall loss rate of ash due to aggregation on the order of teragrams per hour (6 orders of magnitude). The sticking efficiency, Eq. (3), is high (>90 %) for small particles (<63 µm). As the residence time of large particles is very short, the sticking efficiency is applicable to the narrow range of particle sizes that persist in the domain (bins 7–10, <32.5 µm). These particle sizes correspond to a narrow range of sticking efficiencies (0.87 to 0.97), regardless of the water vapor concentration.

4.2 Eyjafjallajökull study results

The ash cloud dynamics generated by WRF-Chem over the model period agree with other modeling studies of Eyjafjallajökull utilizing WRF-Chem (Hirtl et al., 2019; Webley et al., 2012). Figure 7 provides an example of the output from WRF-Chem for 15 and 16 April 2010. The dynamics of the ash clouds are apparent. The plume moves south and east towards the coasts of Scandinavia and northern Europe and then splits into two plumes: one residing over Sweden and Finland and the other passing through multiple northern European countries.

Figure 7 WRF-Chem-generated volcanic ash column densities for the Eyjafjallajökull eruption in April 2010 at 4 h intervals: (a) 15 April at 08:00 UTC, (b) 15 April at 12:00 UTC, (c) 15 April at 16:00 UTC, (d) 15 April at 20:00 UTC, (e) 16 April at 00:00 UTC, (f) 16 April at 04:00 UTC, (g) 16 April at 08:00 UTC, (h) 16 April at 12:00 UTC, and (i) 16 April at 16:00 UTC.

Model output also agrees with airborne in situ measurements. The DLR research aircraft conducted 13 flights on 11 different days that transected Eyjafjallajökull's ash clouds over the course of the phase I and phase III eruptions (Schumann et al., 2011). Predicted ash concentrations from WRF-Chem were compared to the in situ observational data from three of these flights: 19 April and 16 and 17 May 2010. WRF-Chem volcanic ash bins 8–10 correspond to the particle size detection limits of the Grimm OPC and PMS FSSP aboard the Falcon aircraft and were thus chosen for comparisons.

Figure 8 presents time series plots of WRF-Chem output and DLR measurements. Figure 8a, c, and e show the WRF-Chem output in mass concentration (g m⁻³). Figure 8b, d, and f show the WRF-Chem ash bin as number concentrations by using an assumed particle density of 2500 kg  m⁻³ (Brown et al., 2012) in order to make direct comparisons to the Grimm OPC and FSSP detectors.

Figure 8 Comparisons of WRF-Chem model output to in situ mass concentrations (a, c, e) and particle numbers (b, d, f) observed by DLR during 19 April (a, b), 15 May (c, d), and 17 May (e, f) 2010 flights.

Temporal changes in observed and modeled ash concentrations agreed moderately well for the 19 April flight (Fig. 8a and b). Analysis of particle number densities in Fig. 8b for 14 April shows five significant overestimations of volcanic ash by the non-aggregation-enabled code, between 50 % and 75 % at 14:55 and 15:07 UTC, between 15:15 and 15:18 UTC, between 15:35 and 15:42, and between 16:55 and 17:06 UTC. These overestimations did not occur when the aggregation code was used. One peak concentration was observed at 15:30 UTC on 19 April, which was not resolved by WRF-Chem (Fig. 8b). An analysis of the surrounding grid cells in the vertical and horizontal did not contain this peak; however the next vertical grid cell in the positive k contained higher ash concentrations (similar order of magnitude). This analysis, along with analysis of the integrated volcanic ash over the time span of the peak, leads to the conclusion that this lack of peak concentration in the model is a result of model diffusion, which is typical for all Eulerian models. Smaller domain grid cells permit better comparison with point observations, but decreases in grid cell sizes are computationally expensive and in many cases impossible to resolve completely.

Number density readings for 15 May (Fig. 8d) contained more robust data than mass concentration (Fig. 8c) and were therefore used in the analysis. Here, a large overestimation of ash is calculated by WRF-Chem when not using the aggregation code. A peak of 290 particles per cubic centimeter is observed in the unmodified code, almost 10 times higher than observed. With aggregation enabled, the WRF-Chem solution is much closer to the observed numbers at a maximum of 45 particles per cubic centimeter.

On 17 May (Fig. 8e and f), the aircraft performed a steep transect through a plume with larger ash particles. Almost no ash concentration was recorded at the lowest flight altitude reached during the middle of the flight at 16:40 UTC. At this same time, WRF-Chem predicted concentrations in excess of 400 g m⁻³. Where the plume locations do agree, there is improved agreement between the aggregation-enabled code and the airborne observations of mass concentration. For the entire time range, observations where the aggregation code produced mass readings of the same order of magnitude as those observed by DLR were counted. This total was then divided by the total flight time and resulted in an average 80 % agreement of the data (78 % for 19 April, 78 % for 15 May, and 83 % for 17 May). This fell to an average of 62 % when the code was run without aggregation, using the same methodology.

In addition to comparisons with Schumann et al. (2011) in situ measurements, WRF-Chem tephra fallout was also compared to field measurements of tephra collected by Stevenson et al. (2012) in the United Kingdom (UK). Figure 9 depicts the mass of tephra deposited in the model domain from all April 2010 eruptions in Fig. 9a and from May 2010 eruptions in Fig. 9b. Stevenson et al. (2012) report three sampling periods that overlap with the model domain times in this study. For example, Stevenson et al. (2012) counted 218 grains of tephra per square centimeter, at Benbecula in the Outer Hebrides (57.43^∘ N, 7.34^∘ W, Fig. 9a, white circle), with a mean diameter of 18±7 µm while sampling between 13 and 20 May 2010. Assuming an average density of 2500 kg m⁻³ yields a tephra concentration between 20 and 45 mg m⁻², compared to 31 mg m⁻² predicted by WRF-Chem with the aggregation code enabled during the same time range. Samples taken at Leicestershire (52.73^∘ N, 1.16^∘ W, Fig. 9b, white circle) between 25 April and 3 May 2010 estimate a range of tephra mass on the ground between 51 and 119 mg m⁻², also near the WRF-Chem estimate of 41 mg m⁻² (80 % of observed mass) between those dates. Another sample from Lincolnshire (52.74^∘ N, 0.38^∘ W, Fig. 9b, white circle) covered a period from 24 to 30 April 2010. In this case, tephra fallout between 3 and 13 mg m⁻² was measured, whereas WRF-Chem predicted a smaller value of 1.2 mg m⁻² (40 % of observed mass). The smaller estimates for the Lincolnshire and Leicestershire sites may be explained by the lack of model data covering 27 April–3 May, as the last modeled hour was 00:00 UTC on 27 April. When considering WRF-Chem run without aggregation, the modeled fallout seen in these areas is minimal, with less than 1 mg m⁻² observed.

Figure 9 Mass of tephra fallout deposited on model surface and lowest model level in WRF-Chem, for April (a) and May (b) 2010 model simulations. The white circle in (a) marks the Outer Hebrides, and the white circle in (b) marks Lincolnshire and Leicestershire, UK, corresponding to sample areas in Stevenson et al. (2012). Maximum domain fallout is 52 Mg m⁻².

The aggregation code altered the total domain mass of each volcanic ash bin. To study this change, the model domain mass was analyzed from 14 to 18 May 2010. This time frame represents the last 96 h of modeled eruptions and includes a high degree of variability in the eruption rate and plume height (see Fig. 2). The total domain mass is presented in Fig. 10 without (a) and with (c) the aggregation code enabled. To analyze the PSD, the mass of each volcanic ash bin was divided by the total model domain mass. The resulting percentages are presented in Fig. 10b and d. The top panels, Fig. 10a and b, depict WRF-Chem output without the use of the aggregation code, whereas the lower panels, Fig. 10c and d, include the aggregation code. The short atmospheric lifetime of the large particles in bins 1–3 results in small masses during this time frame compared to bins 4–10 including smaller particle sizes. As such, only bins 4–10 are depicted in Fig. 10. Major changes in the eruption rates are annotated on the time axis with red marks.

Figure 10 Total domain ash mass (a, c) and percent contribution to domain mass (b, d) for the modeled period between 14 and 18 May 2010 without (a, b) and with (c, d) aggregation code enabled. Red numbers on date/time axis denote major (>10 %) changes in the eruption rate: (1) 14:00/09:00 Z – decrease from 7.949 to 1.775 Tg every 3 h; (2) 16:00/06:00 Z – increase from 1.175 to 7.949 Tg every 3 h; (3) 17:00/00:00 Z – decrease from 7.949 to 1.056 Tg every 3 h; (4) 17:00/15:00 Z – decrease from 7.949 to 0.966 Tg every 3 h. Note there are variable increases and decreases in the eruption rate between markers 3 and 4.

Two important observations are noted when aggregation is included. First, the total domain mass in each bin is reduced, and second, the PSD shifts towards smaller-sized particles during eruptive events. For example, the initial period in Fig. 10 is eruptive until the first red mark on 14 May at 09:00 UTC. During this period, the eruption rate is 7.36×10⁵ kg s⁻¹ (7.949 Tg every 3 h). In the non-aggregation-enabled code, the dominant ash species are bins 6–8, which have peak masses of 3.7, 4.1, and 3.3 Tg, respectively. In the non-aggregation-enabled code, bins 6–8 make up the majority of the domain as mass, contributing 21.5 %, 24.1 %, and 19.3 % of the total domain mass. When the aggregation code is enabled, the total domain mass for each of the bins is reduced to 1.0, 1.5, and 1.4 Tg, respectively, which is around one-third of the original peak mass, showing an overall reduction. Additionally, their contribution to the overall domain mass changes to 14.6 %, 21.1 %, and 20.5 %. The smaller bin 8 ends up with more of the mass, with the other two contributing less to the PSD. In fact, the smaller bins 9 and 10 also contribute more to the overall domain mass, increasing from a peak of 13.1 % and 9.6 % on 14 May at 09:00 UTC without the aggregation code enabled to 15.2 % and 11.6 % with the aggregation code enabled. Overall there is a slight shift towards smaller particle bins during eruptive events.

Interestingly, this trend in the PSD is not observed during periods of decreased eruption rates, while trends in overall domain mass continue are still observed. Between marks 1 and 2, the eruption rate decreases from 7.36×10⁵ to 1.09×10⁵ kg s⁻¹. During this period of slower eruption rates, the total domain mass continues to increase; however it is much lower when aggregation is considered. The PSD, on the other hand, remains consistent, with bins 8–10 trending similarly in the non-aggregation- and aggregation-enabled cases. This suggests that the aggregation code is most effective during eruptive events when particles are in high concentration.

Without aggregation, the only sinks for volcanic ash are via settling, which is dependent on gravity and water vapor concentration, or via the plume traveling out of the model domain. For finer ash particles, removal via settling is minimal when compared to larger particles, which is evident in Fig. 10a and c. During periods of less volatile eruptions, such as between markers 1 and 2 or markers 3 and 4, the fine ash bins reach a steady state where the source of ash is almost equal to the sink, i.e., settling. This is evident in the horizontal slope of the bin domain mass. This is not true for larger particles whose settling velocities are high enough to remove them faster than they are added. Aggregation adds an additional sink that is noticed subtly during less eruptive phases as the slight dips in domain mass, as well as the more pronounced decreases in the slope of the change in domain mass during periods of higher eruption rates.

5 Summary and conclusions

A parameterization of volcanic ash particle aggregation has been implemented into the fully coupled WRF-Chem model. The new model has been tested for ash loadings and lifetimes. A simplified version of the Smoluchowski coagulation equation (Costa et al., 2010; Dekkers and Friedlander, 2002; Folch et al., 2010, 2016; Smoluchowski, 1917) was incorporated into the WRF-Chem model. This simplified method was chosen for its computational efficiency, allowing the aggregation rate to be calculated at each model time step in line with the atmospheric dynamics.

The effects of the aggregation code were assessed by applying it to a high-resolution model study of the 2010 eruptions of Eyjafjallajökull, including a single study of a 9 h test eruption. The effect of each particle collision kernel on the overall aggregation rate (Eq. 2) was studied. The degree to which each kernel affected aggregation depended on the choice of the fractal dimension, D_f. The differential sedimentation kernel provided the largest contribution by orders of magnitude when a fractal dimension of 3.0 was chosen; however the Brownian kernel dominated when a fractal dimension of 2.5 was chosen. This result suggests that vertical motion, when a fractal dimension near 3.0 is chosen, is the primary driving force behind particle interactions in the aggregation process, rather than random (Brownian) or horizontal (shear) motions. Additionally, analysis of the volcanic ash lifetime shows that varying the fractal dimension may greatly vary the lifetime, especially when considering fractal dimensions between 3.0 and 2.8.

The Eyjafjallajökull model study was assessed by comparison to aircraft in situ measurements taken by DLR as well as tephra fallout samples measured in the United Kingdom. By comparing WRF-Chem-calculated volcanic ash mass concentrations using the aggregation code to those observed by DLR, an average 80 % match in order of magnitude was observed for the three flights analyzed. Additionally, non-aggregation-enabled code calculated 20 %–50 % higher volcanic ash concentrations on numerous occasions, where the aggregation-enabled code did not. The aggregation-enabled WRF-Chem code tended not to overestimate volcanic ash, or to overestimate less than the non-aggregation-enabled version, potentially yielding more realistic ash concentrations which may benefit aircraft hazard mitigation forecasting.

As the plume was transported over the United Kingdom, WRF-Chem predicted ash fallout that compared well to field measurements. Tephra fallout generated by WRF-Chem fell within observed values at one sample location and predicted on average 60 % of the fallout at two others. This suggests that WRF-Chem may be used to model not only the atmospheric transport of ash clouds but also the deposition of ash as well.

Importantly, these observations all suggest that two factors drive volcanic ash aggregation when including aggregation in the WRF-Chem code. First, volcanic ash concentration is noted to be the primary driving factor behind aggregation rate. The majority of model domain mass decreased near the vent where concentrations of ash are high. In addition, PSD analysis indicates that bins with higher portions of the eruption PSD undergo faster rates of ash aggregation. Bins with a larger share of the eruption PSD will aggregate faster due to their increased probability of collision. Second, vertical motions of ash falling through the atmosphere also drive the aggregation process through differential sedimentation for realistic ranges of fractal dimension (between 2.95 and 3.0).

The inclusion of this aggregation scheme into WRF-Chem provides research and operational meteorological communities a second VATD model to Fall3D that includes volcanic ash aggregation and is the first to run aggregation in an inline fashion where aggregation equations are solved at each model time step (Folch et al., 2010). This inline computation of volcanic ash yields many benefits. For example, the code identifies the driving forces behind volcanic ash aggregation, i.e., ash concentration and differential sedimentation rates, and allows for the study of the effects of water vapor concentration on the aggregation rate. In addition, it allows the study of changes in particle size distributions due to enhanced ash settling as a result of aggregation processes, which are of particular importance to remote sensing communities where the effective particle size directly impacts the spectral methods used for detection. While this study focused primarily on the distal ash cloud transport and aggregation physics, the calculations integrated into WRF would also benefit a higher-resolution, nested domain over the emission source to study proximal aggregation effects. The modified code also benefits the operational volcanic ash modeling community by providing model-derived ash mass concentrations that augment existing VATD models for use in aircraft hazard mitigation. In the operational setting, first-guess, expedient model output from VATD models can be augmented by WRF-derived mass loadings as they become available. The time requirement for this is feasible in the operational setting as the modified code is computationally expedient. It ingests output from global models, such as ECMWF and GFS, and runs volcanic ash dispersion and aggregation code while simultaneously calculating mesoscale atmospheric dynamics, eliminating the need for additional, offline calculations. Additionally, this code results in another model that provides researchers a robust treatment of ash microphysical processes as they are erupted, transported, and removed from a model domain. Ultimately, this study provides another step towards the inclusion of volcanic ash aggregation, an important physical process, into VATD models.