Plume-SPH (Cao et al., 2018) is designed to describe an injection of well-mixed solid and volcanic gas from a circular vent above a flat surface into a stratified stationary atmosphere. The basic assumptions of the model are as follows:

1) Molecular viscosity and heat conduction are neglected since turbulent energy and momentum exchange are dominant;

Based on above assumptions, the governing equations of our model are given as: ∂ ρ ∂ t + ∇ ⋅ ρ v = 0 (1) ∂ ρ ξ ∂ t + ∇ ⋅ ρ ξ v = 0 (2) ∂ ρ v ∂ t + ∇ ⋅ ρ v v + p I = ρ g (3) ∂ ρ E ∂ t + ∇ ⋅ ρ E + p v = ρ g ⋅ v (4) where ρ is the density, v is the velocity, ξ is the mass fraction of ejected material, g is the gravitational acceleration, I is a unit tensor. E = e + K is the total energy which is a summation of kinetic energy K and internal energy e. An additional equation is required to close the system. In this model, the equation for closing the system is the following equation of state (EOS). p = γ m − 1 ρ e (5) where γ m = R m / C v m + 1 (6) R m = ξ g R g + ξ a R a (7) C v m = ξ s C v s + ξ g C v g + ξ a C v a (8) ξ a = 1 − ξ (9) ξ g = ξ ⋅ ξ g 0 (10) ξ s = ξ − ξ g (11) where C _v is the specific heat with constant volume, R is the gas constant. ξ is the mass fraction of erupted material. The subscript m represents the mixture of ejected material and air, s represents solid portion in the ejected material, g represents gas portion in the ejected material, a represents air, 0 represents physical properties of erupted material. ξ _g0 is the mass fraction of vapor in the erupted material.

Three different boundary conditions are applied in this model. At the vent, temperature of erupted material T, eruption velocity v, the mass fraction of vapor in erupted material ξ _g0, and mass discharge rate M ̇ are given. The pressure of erupted material p is assumed to be the same as ambient pressure for pressure-balanced eruption. The radius of the vent is determined from ρ, M ̇ and v. Nonslip wall boundary condition is applied to the flat ground, where we enforce the velocity to be zero. With further assumption that the ground is adiabatic, internal energy flux, which consists of heat flux and energy flux carried by mass flux, vanishes on the wall boundary. Pressure outlet boundary condition is applied to the surrounding atmosphere where the pressure is given. Except for the pressure, boundary values for density, velocity, and energy are determined by numerical calculation from the conservation laws. The initial condition for Plume-SPH is created based on the atmosphere profile before the eruption.

The governing equations, EOS, boundary conditions, and initial conditions establish a complete mathematical model. The model posed over the computational domain is then discretized using smoothed particle hydrodynamics (SPH) method (Gingold and Monaghan, 1977) available in the tool Plume-SPH (Cao et al., 2017, 2018) using two types of SPH particles: 1) particles of phase 1 to represent ambient air and 2) particles of phase 2 to represent erupted material. So before the eruption, the computational domain is fully occupied by particles of phase 1. During the eruption, particles of phase 2 are injected into the computational domain. The discretized model is then converted into a large computation task in the Plume-SPH tool based on a parallel data management framework (Cao et al., 2017).

The input parameters for Plume-SPH include the eruption condition at vent, the material properties, and a profile of the atmosphere. The eruption parameters, material properties, and atmosphere for the “strong plume–no wind” case in the recent comparison study on eruptive column models (Costa et al., 2016) are adopted. Eruption conditions and material properties are listed in Table 1. Note that the density of erupted material at the vent and radius of the vent can be computed from the given parameters. The eruption pressure is assumed to be the same as the atmospheric pressure at the vent; hence, it is not given in the table. The vertical profiles of atmospheric properties were based on the reanalysis data from European Centre for Medium-Range Weather Forecasts (ECMWF) for the period corresponding to the climactic phase of the Pinatubo eruption.

Running of Plume-SPH updates physical quantities, such as temperature, velocity, and the position of SPH particles in each time step. During Plume-SPH simulation, SPH particles of phase 2, which represent the erupted material, are injected from the eruption vent into the computation domain with an initial injection velocity. As they move upwards, these particles will get mixed with SPH particles of phase 1, which represent the air, during the whole simulation. Their physics quantities get updated as well. After the simulation, the computation domain will be filled with SPH particles of both phase 1 and phase 2. Removing all SPH particles of phase 1 from the computation domain, all of the remaining SPH particles represent the erupted material, which naturally forms a plume (see Figure 1).

Puff (Tanaka, 1991; Searcy et al., 1998) is a dynamic pollutant tracer model. The model is based on a 3D Lagrangian form of the fluid mechanics, in which the material transport is represented by the fluid motion, and diffusion is parameterized by a stochastic process of random walk. Here, the model is constructed by a sufficiently large number of Lagrangian tracer particles with a random variables R _i (t) = (x(t), y(t), z(t)), where i = 1 ∼ M, which represents position vectors of particles from the origin of the ash source at the time t. M is the total number of Lagrangian tracer particles, a sample of all the ash particles. R i t + Δ t = R i t + W t Δ t + Z t Δ t + S i t Δ t (12) Here, W accounts for local wind advection, Z is generated by Gaussian random numbers and accounts for turbulent dispersion, and S is the terminal gravitational fallout velocity or settling speed, which depends on a tracer’s size.

A collection of tracer particles can be used to start a Puff simulation. The tracer particles have three basic properties: age, size, and position. The age of each particle is the elapsed time from when it was released. Ash particles in the initial ash cloud have zero age. Initial ash size distribution is assumed to be log-normal. According to a mean and standard deviation provided by the user, Puff assigns size to each particle. Puff initializes the position of each particle according to semiempirical expressions. The height of each particle is determined according to the specified distribution from the surface (1000 mbar ≅ 0 m) to the top of the plume height, H _max , which is given by the user. Puff also supports reading predefined initial ash clouds from a file, containing the coordinates of all tracer particles.

Vertical particle distribution in Puff is usually based on the Poisson distribution. For the Poisson distribution, the vertical height of ash particles is given by Eq. 13: H = H m a x − 0.5 H w i d t h P + H w i d t h R (13) where P is an integral value drawn from a Poisson distribution of unit mean, R is a uniformly distributed random number between 0 and 1, H _max is the maximum plume height, H _width represents an approximate vertical range over which the ash will be distributed. So for Poisson distribution, the user can specify two parameters, H _max and H _width . Another commonly used vertical ash distribution in VATD simulation is Suzuki. For the Suzuki plume shape (Suzuki, 1983), the ash mass vertical distribution is assumed to follow the Eq. 14: Q z = Q m k 2 1 − z / H m a x exp k z / H m a x − 1 H m a x 1 − 1 + k exp − k (14) where Q _m is the total mass of erupted material, k is the shape factor, which is an adjustable constant that controls ash distribution with height. A low value of k gives a roughly uniform distribution of mass with elevation, while high values of k concentrate mass near the plume top. So for Suzuki distribution, besides the plume height H _max , there is another user-specified parameter, k.

Puff initializes the horizontal distribution of ash particles according to semiempirical expression as well. Puff uses a uniformly distributed random process to determine ash particle locations in a circle centered on the volcano site. The maximum radius (at plume top) at which a particle can be located is given as “horizontal spread”. The horizontal displacement from a vertical line above the volcano is a random value within a circle of which the radius equals the “horizontal spread” multiplied by the ratio of the particle height H to the maximum H _max , see Eq. 15. So the resulting shape of the particle distribution within the plume is an inverted cone in which particles are located directly over the volcano at the lowest level and extend out further horizontally with increasing plume height. r H = r m a x H / H m a x R (15) where r(H) is the radius of the horizontal circle, within which all particles at the height of H are located. r _max is the horizontal spread. H is the height, and R is an uniformly distributed random number between 0 and 1.

In summary, particle distributions in the initial ash cloud are controlled by several parameters, for example, H _max , H _width , and r _max if the user chooses to use semiempirical expressions, Eqs 13, 15. Users can optimize or calibrate these parameters to adjust the initial condition for Puff so that the simulated results match better with observations. Besides the initial ash cloud, other input parameters for Puff are diffusivity in the vertical and horizontal directions, start and end time of the eruption, and eruption duration. When creating initial conditions from the output of Plume-SPH, the total number of Lagrangian tracers is the count of all SPH particles of phase 2 in the plume. The same total number of Lagrangian tracers is used when creating the initial ash cloud based on semiempirical expressions. All input parameters for Puff are listed in Table 2.

In this study, we convert the output of Plume-SPH into an initial ash cloud which serves as the initial condition for Puff. The method proposed consists of generating the initial ash cloud directly from Plume-SPH, foregoing assumptions and estimates, or inverse modeling, regarding ash injection height and timing. The steps to create an initial ash cloud based on the raw output of Plume-SPH are shown in Figure 1. The initial ash cloud is created from SPH particles of phase 2, which represents the erupted material in the model. After reaching the maximum rise height and starting to spread horizontally, particles of phase 2 form an initial umbrella cloud (Figure 2). The 3D plume simulation is considered complete once the umbrella cloud begins to form. Parcels that will be transported by the ambient wind are those above the “corner” region, where mean plume motion is horizontal rather than vertical. With such consideration, we introduce an elevation threshold, which is the lower elevation limit of the ash that will be transported by the VATD. All SPH particles with elevation lower than the threshold are excluded when creating the initial ash cloud. The inflection point from vertical raising to horizontal spreading happens around 15 km according to the averaged vertical velocity [(d) in Figure 2] and horizontal velocity [(e) in Figure 2]. Below this inflection point, particle trajectories are primarily vertical in the stalk-like eruption column. Above this level, particle trajectories are primarily horizontal, as they flow into the umbrella cloud gravity current. So we choose 15 km to be the elevation threshold in this study.

Considering that SPH particles are only discretization points, each is assigned a grain size according to a given total grain size distribution (TGSD) (Paladio-Melosantos et al., 1996), and a concentration according to the mass and volumetric eruption rate. The Plume-SPH discretization points are thus switched to Puff Lagrangian tracer particles having grain sizes and concentrations. The coordinates of these tracer particles, which are initially in the local Cartesian coordinate system of Plume-SPH, are converted into Puff’s global coordinate system, which is given in terms of (longitude, latitude, height). Puff takes the initial ash cloud, consisting of the collection of Lagrangian tracer particles with grain size and concentration, and propagates from time t to time t + Δt via solution to an advection/diffusion equation (Eq. 12).

To summarize, there are four steps to create an initial ash cloud from the raw output of Plume-SPH:

1) filter by SPH particle type to select SPH particles that represent erupted material (phase 2);

The features of the volcanic plume and resulting initial ash cloud used in the case study are shown in Figure 2. It is important to point out that since both Plume-SPH and Puff are based on the Lagrangian method, there is no extra step of conversion between an Eulerian grid and Lagrangian particles.

The plume and ash transport models are run at different time scales and length scales. The spatial and temporal resolutions of the plume simulations are much finer than those of the ash transport model. It takes tens of minutes (600 s in this case) for the Pinatubo plume to reach a steady height. However, the eruption persisted for a few hours (9 h for the climactic phase of Pinatubo eruption), and it may be necessary to track ash transport for days following an eruption. At present, it is too computationally expensive to run 3D plume simulations for several hours in real time. In order to handle the difference in time scale, we mimic a continuing eruption with intermittent pulses releasing ash particles. In particular, we restart Puff at an interval of 600 s, i.e., the physical time of the plume simulation to reach a steady height. At every Puff restart, we integrate the output of the last Puff simulation and Plume-SPH into a new ash cloud. This new ash cloud serves as a new initial condition with which to restart a Puff simulation. A sketch demonstrating the overall restart process is shown in Figure 3. The total number of Lagrangian tracer particles used in Puff thus equals the summed number of particles in all releases. The total number of tracer particles is therefore no longer a user-selected parameter. Fero et al. (2008) proposed using more realistic time-dependent plume heights. We do not adopt that strategy here for simplicity, although the idea would be straightforward in execution, given time-dependent eruption conditions.

In this section, we discuss the vertical distribution of ash particles in the initial ash cloud. The majority of volcanic ash particles are usually injected at an elevation lower than the plume height. For instance, Holasek et al. (1996a,b) reported the maximum Pinatubo volcanic plume, i.e., source height, as ∼ 39 km while the distal volcanic cloud heights were estimated at ∼ 20–25 km. Self et al. (1996) reported that the maximum plume height could have been >35 km, but that cloud height was ∼23– 28 km after ∼ 15–16 h. The neutral buoyancy height of the Pinatubo aerosol cloud was estimated with different methods at ∼ 17–26 km (lidar) by DeFoor et al. (1992), ∼ 20–23 km (balloon) by Deshler et al. (1992), ∼ 17–28 km (lidar) by Jäger (1992), and ∼ 17–25 km (lidar) by Avdyushin et al. (1993). Based on a comparison between simulated clouds with early infrared satellite imagery, Fero et al. (2008) reported that the majority of ash was transported between 16 and 18 km. These observations make good physical sense, as particles are concentrated or centered around the intrusion height of the umbrella cloud, not near the plume top, because the plume top is due to momentum overshoot. However, the empirical expressions for the height-MER relation, which are commonly adopted to create initial conditions for VATD simulations, tend to place the majority of ash particles closer to the top if one uses observed plume height in the empirical expressions.

Here, we investigate two commonly used plume shapes, the Poisson (see Eq. 13) and Suzuki (see Eq. 14). Particle distributions (in terms of mass percentage or particle number percentage) in the vertical direction in the initial ash cloud are shown in Figure 7. In that figure, the vertical particle distribution based on Plume-SPH output is compared with the vertical particle distribution based on semiempirical shape expressions. Both Poisson and Suzuki distributions in Figure 7 take H _max = 40 km, which is close to the reported observed plume height. When adopting the Poisson distribution [(c) in Figure 7], the majority of the particles are between 30 and ∼ 40 km. Obviously, the Poisson function distributes the majority of ash at a higher elevation than was observed (e.g. Fero et al., 2008). As for the Suzuki distribution, (D) in Figure 7, the majority of ash particles also occur in a range that is significantly higher than 25 km. Note that in the plot (d), the Suzuki constant k is set to 4, which is commonly used for Sub-Plinian and Plinian eruption columns (Pfeiffer et al., 2005). As for initial ash clouds in Plume-SPH simulations, most ash particles are distributed between ∼ 17–28 km, which matches well with observations. The plume height is also consistent with observation.

For the Poisson distributions, the ash particles cannot be lower without changing the plume height. To distribute the majority of ash particles at a lower elevation, the plume height must be reduced to a value smaller than the observed plume height. Adjusting parameters such as plume height in the empirical expression is actually the traditional source term of calibration method. A set of initial ash clouds using different plume heights based on the Poisson distribution is shown in Figure 8. The plume heights adopted in plume shape expressions are not obtained from any plume model or observation of plume height, but by a posteriori calibration to later-observed ash cloud transport heights. For Suzuki distribution, adjusting the Suzuki constant can adjust the distribution of ash particles in vertical direction. As shown in Figure 7, when k is equal to 1 [see (e)], the majority of ash particles are at a lower elevation than observation. With k = 3 and k = 6 [figure (g) and (h)], the majority of ash particles are at a higher elevation than observation. When k is set to 2 [see (f)], we can see that the majority of ash particles are roughly distributed in the range 17–28 km. But the shape does not look like a typical plume, as particles are more uniformly distributed in the vertical direction. In addition, the “best fit” Suzuki constant is different from the typical value, which is 4 (Pfeiffer et al., 2005), for Sub-Plinian and Plinian eruptions, meaning that we can not apply previous experiences into the semiempirical expression for this eruption.

The ash clouds created by the Poisson distribution with different plume heights are used as initial conditions in Puff simulations, whose results are shown in Figure 9. Except for the plume height, all other parameters for creating an initial ash cloud are the same as those in Table 2. Of course, the range over which the majority of ash particles are located is lower when using lower plume heights. Figure 9 thus shows that the plume height has a significant influence on the ash transport simulation. The maximum heights of the simulated ash cloud are completely different when using different H _max values in the Poisson expression. When the plume height is 10 km, the ash lags behind that observed and its FMS is 0.055, which is very close to FMS when H _max is 40 km. For the cases that H _max is 20 and 30 km, the FMS values are 0.121 and 0.142, respectively. Taking 20 km as the plume height better represents the lower elevation portion of the ash cloud while taking 30 km as the plume height better represents the higher elevation portion of the ash cloud.

Simulation results based on a calibrated plume height of 30 km show a footprint similar to those of “Plume-SPH + Puff” although smaller in terms of area. However, the initial ash cloud created by a Poisson distribution with a plume height around 35 km generates the best match with observation in terms of FMS metric, with the FMS value reaching 0.227. That is to say, a plume height lower than the real plume height is required by the Poisson plume shape to distribute ash particles at elevations comparable to the “true” ash distribution. Even for the best-matched results, the high concentration area does not match with an observation well.

It is clear that the initial condition of vertical ash distribution has a dominant effect on VATD simulation, so it is critical for the forecast capability of VATD simulations to explore more accurate and adaptive ways for establishing the initial ash distribution, especially methods that do not rely on a posteriori parameter calibration or inversion.

In the previous section, we explored the effects of adjusting the plume height to change the vertical ash distribution at the source. In this section, we investigate the importance of another parameter in the semiempirical Poisson expression (Eq. 13). We vary the “vertical spread”, H _width , in the range ∼ 3–10 km. A set of initial ash clouds with different vertical spreads are shown in Figure 8. Except for vertical spread, all other parameters for creating an initial ash cloud are the same as those in Table 2. The vertical range within which the majority of ash particles are located becomes narrower when a smaller value for the vertical spread parameter is used. The ash clouds based on different vertical spread parameters are then used as initial conditions in Puff simulations.

The VATD results are shown in Figure 9. Adjusting the vertical spread changes particle distribution in the vertical direction, and thus, not surprisingly, affects the VATD simulation results. None of the VATD simulations based on initial ash clouds with vertical spreads equal to 3 km or 5 km yield better results than do VATD simulations based on initial conditions created by Plume-SPH (see Figure 9). But when we take 10 km as the vertical spread, we get a FMS that is very close to Plume-SPH, even though the shape of the ash cloud footprint and the maximum height of the ash cloud are completely different.

The calibration tests on vertical spread, carried out here, are certainly not exhaustive. One could do a more comprehensive calibration throughout the multidimensional parameter space (for Poisson distribution, the parameter space is two-dimensional) and find better results. In addition, with a more complicated semiempirical plume shape expression, one could have more control over plume shape and might be able to get an initial condition that yields a more accurate ash transport forecast. However, more complicated and adaptable plume shape expressions imply a higher-dimensional parameter space, which requires more effort in calibration, even though the degrees of freedom to adjust plume shape are still limited. Creating initial conditions based on 3D plume simulations avoids such parameter calibration.

The differences between the semiempirical plume particle distribution and actual (or simulated by the 3D plume model) are not only in the vertical direction. The importance of the horizontal distance of each initial ash particle from a line extending upward from the volcano is investigated in this section. Puff uses a uniformly distributed random process to determine ash particle locations in a circle centered on the volcano site as described in Puff and Initial Ash Cloud for the output of Plume-SPH, and an effective (maximum) radius is determined according to a given threshold of ash concentration, following Cerminara et al. (2016b). A time-averaged, spatial integration of the dynamic 3D flow field is conducted to remove significant fluctuations in time and space. Figure 10 compares the radius of the initial ash clouds created by 3D plume simulations with that assumed in the semiempirical plume shape expression adopted in Puff. It is impossible for the simple, assumed plume shapes to capture the complex and more realistic shapes developed by Plume-SPH. Additional parameterization may generate more reasonable shapes, but these would continue to be ad hoc; none would likely have the potential fidelity of the 3D simulation to reality, and adding a temporally changing distribution would be difficult.

Comparison between cross-sectional views of the initial ash clouds is shown in Figure 11. The cross-sectional view of horizontal particle distribution using the semiempirical method (last figure in Figure 11) is similar to a cross-sectional view of a simulated 3D plume, in a general sense. However, for simulated 3D plumes, the ash particle distribution in cross section varies with height, which factor would become increasingly important with increasing wind speed, where wind speed to be included in the estimate of initial plume shape. It is difficult for the semiempirical expressions to accommodate such a complex distribution.

Despite the obvious difficulty of correctly estimating ash distribution near the vent, or for short propagation times, assigning different values for the horizontal spread has a negligible effect on VATD simulation results at large time. We investigated horizontal spread values between 50 and 1600 km to create initial ash clouds; all of them generated similar results at large propagation times (> 1 day). Figure 9 shows two different simulation results based on initial ash clouds with horizontal spread equal to 50 and 600 km, respectively. No visible differences are apparent between them. The FMS values, 0.073 and 0.074, respectively, are also very close. This implies that horizontal distribution has a less significant influence on VATD simulation results than does vertical distribution for long distance or large time. Perhaps, the most important ramification of this result is that it means the time at which the “handshake” is made between Plume-SPH and the VATD does not affect results significantly for relatively large distances and times.

Besides the initial ash cloud, other parameters for Puff simulations are horizontal diffusivity, vertical diffusivity, mean grain size, grain size standard deviation, and total number of tracers. We present in this subsection informal sensitivity studies on these parameters. We also investigate the influence of eruption duration. The sensitivity analyses will serve as the basis for identifying possible sources of disparities between simulation and observation.

Fero et al. (2008) simulated the volcanic ash transport in the Pinatubo eruption in 1991. He carried out systematic sensitivity analysis with respect to input parameters of Puff and found that all other parameters except for the plume height have negligible effect on long-term ash transportation of Pinatubo. Inspired by Fero et al. (2008), we carried out similar informal sensitivity analysis with much fewer sample points in the parameter space and got similar results. Among the parameters explored, the eruption duration and beginning time show the most obvious influence on simulated ash distribution although the effect is still small. To show the differences in an intuitive way, (a) - (c) in Figure 12 shows simulated ash distribution corresponding to 4.9 h duration, 9 h duration, and 11 h duration, respectively. After 72 h, relative to the simulation starting time, these three cases generate very similar results with tiny visible differences. Daniele et al. (2009) did sensitivity analysis with respect to the input parameters of Puff on different volcanoes and found that for eruptive eruptions, the most dominant factors are the wind field and plume height, while all other input parameters are relatively less important. The significance of the wind field has been confirmed by other researchers (e.g. Stefanescu et al., 2014).

We conducted several simulations with eruption duration varying in the range of 5–11 h with slightly different starting time of climactic phase. Table 3 lists all these simulations. However, only slight visible differences are observed among the simulated ash transport outputs. We can see that the eruption duration has negligible effects on long-term ash transport.

The new methodology for generating initial ash clouds introduces a new parameter: elevation threshold, which was specified based on averaged vertical velocity and horizontal velocity. We carry out a separate, informal sensitivity analysis on this parameter by varying the elevation threshold from 1.5 km (the height of the vent) to 25 km. The simulated ash distributions show obvious differences, especially when the elevation threshold is either very high or very low. However, varying the elevation threshold in the range of 12–18 km generates relatively small differences in ash transport simulation results. Figure 12 (d) and (f) compare the simulated ash distributions corresponding to elevation thresholds of 1.5 and 15 km. Compared with the ash distribution for a threshold of 15 km, an extra-long tail appears when using an elevation threshold of 1.5 km. The maximum height of the tail is around 10 km. Adopting lower elevation thresholds adds more tracer particles at lower elevation. As the winds at different elevations are different, the tracers at lower elevations propagate in different directions. The HYSPLIT forward trajectory tracking indicates that the wind between elevations of 10 and 15 km blew from north-east to south-west, while winds of higher elevation blew from east to west (see Figure 6).

The full range of research issues raised by numerical forecasting of volcanic clouds is diverse. We focused on the effect of initial conditions in this study. During the plume modeling, secondary factors, such as microphysical processes, even though they play lesser roles, likely need to be included to improve accuracy for a particular eruption. Wind fields are not considered in the current version of Plume-SPH, but for weak plumes, wind plays an important enough role that it has to be considered in the plume model. In addition, eruption conditions are subject to change with time, even during the climactic phase of an eruption. For example, ash just west of Pinatubo observed in satellite images does not show up in “Plume-SPH + Puff” simulation results. This disparity is likely due to the fact that Pinatubo continued erupting (with smaller plume height) after the climactic phase, while we only simulate the climactic phase. In the future, time-dependent initial conditions for VATDs can be created from 3D plume simulations based on time-dependent eruption conditions. It is worth mentioning that the eruption conditions at the vent are usually inferred from observable information based on 1D plume models. Using a 3D plume model will not reduce uncertainties from the eruption conditions.

Additional assumptions made during computations in each VATD model or even measurements may also generate additional disparity. Analysis of the results (see large decrease in FMS shown in (f) in Figure 5) indicates that Puff underestimates the fallout of ash particles, which together with satellite pictures not capturing low-level ash clouds can explain the FMS decrease.

One implicit assumption in the current method is that ash transportation is dominated by wind advection (the passive dispersion approximation). However, during the growth of the volcanic umbrella, the dominant factors are various in different regimes (Pouget et al., 2016a) depending on the characteristics of a particular eruption. Webster et al. (2020) suggested that the lateral spread by the intrusive gravity current dominates the transport of the ash cloud in this stage. Studies by Mastin (Mastin et al., 2014; Mastin and Van Eaton, 2020) also showed that neglecting the umbrella cloud formation for larger eruptions led to significantly different footprints for the resulting VATD fallout maps. Their studies imply that including mapped velocities of the plume as a perturbation on the winds can better capture the radial spreading of an umbrella. In the current method, the 3D plume model generated initial ash cloud has a radius of around 25 km. For the Pinatubo 1991 eruption, the passive dispersion approximation can be reasonably applied when radius is greater than 450 km and can be fully valid only when the radius is greater than 1800 km (Costa et al., 2013). So the umbrella stage during the ash transportation is very likely oversimplified in the current simulation. It is computationally too expensive for the Plume-SPH model to continue simulation until the plume radius reaches, at least, for example, 450 km. An additional umbrella model, with a much coarser resolution and simplified physics, in between the plume model and the VATD model would presumably better model the whole ash transportation process.

Besides the errors from assumptions in the model, errors are also introduced from the reanalysis wind field data and the satellite observations, which are retrievals, with their associated errors, rather than the “truth”. In addition, metrics based on footprint cannot account for the disparities at different heights and ash concentrations. Comparing the simulation and observation purely based on footprint-based metric sometimes is biased.