An Adaptive Semi-Lagrangian Advection Model for Transport of Volcanic Emissions in the Atmosphere

Dispersion of volcanic emissions in the Earth atmosphere is of interest for climate research, air traffic control as well as human wellbeing. Current volcanic emission dispersion models rely on fixed grid structures that often are not able to resolve the fine filamented structure of volcanic emissions while being transported in the atmosphere. Here we extend an existing adaptive semi-Lagrangian advection model for volcanic emissions including the sedimentation of volcanic ash. The advection of volcanic emissions is driven by a pre-calculated wind field. For evaluation of the model, the explosive eruption of Mount 5 Pinatubo in June 1991 is chosen, which was one of the largest eruptions in the 20th Century. We compare our simulations of the climactic eruption on June 15, 1991 to satellite data of the Pinatubo ash cloud and evaluate different sets of input parameters. We could reproduce the general advection of the Pinatubo ash cloud and owing to the adaptive mesh, simulations could be performed at a high local resolution while minimizing computational cost. Differences to the observed ash cloud are attributed to uncertainties in the input parameters and the pass by of Typhoon Yunya, which is probably not completely resolved in the 10 wind data used to drive the model. Best results were achieved for simulations with multiple ash particle sizes.


Introduction
Tephra and SO 2 emissions from large volcanic eruptions have a crucial impact on short-and long-term climate variations, air traffic and the living conditions of people in the surrounding of volcanoes.Large tropical and high-latitude eruptions were primary drivers of interannual-to-decadal temperature changes in the Northern Hemisphere during the last 2500 years (e.g., Sigl et al., 2015).However, even smaller volcanic eruptions do significantly affect the living conditions on a local scale.For example, the respiration of volcanic ash and gas (e.g., Horwell and Baxter, 2006) is along with the fall out of tephra (e.g., Paladio-Melosantos et al., 1996) the most important impact on the local scale.Heavy tephra fall can lead to collapse of buildings, destruction of mechanical and electrical systems, disruption of transport systems, formation of enormous lahars, chemical and physical changes in water quality, and damage of vegetation, crops, forestry and pastures (Folch, 2012).Drifting ash clouds pose a serious threat to jet aircraft and can lead to engine failure (e.g., Casadevall, 1993).Since 1976 an average number of two damaging encounters per year between aircraft and ash clouds has been reported (Guffanti et al., 2010), andClarkson et al. (2016) lately reviewed available engine and volcanological data and proposed a new "safe-to-fly" chart with a much lower ash concentration threshold than previously recommended.
Volcanic SO 2 injected into the stratosphere has a global impact due to its conversion to sulfate aerosol which disturbs the Earth's radiation balance.Tropical volcanic erup-Published by Copernicus Publications on behalf of the European Geosciences Union.E. Gerwing et al.: An adaptive semi-Lagrangian advection model tions thereby lead to warmer winters and colder summers on the Northern Hemisphere continents through dynamical feedbacks and radiative forcing, respectively (Robock, 2000).In addition, volcanic aerosols lead to an increase in stratospheric particle surface area, enhancing the ozone destruction especially in high latitudes (Solomon, 1999).The amplitude of the diurnal cycle of the surface air temperature is reduced by volcanic tephra remaining in the atmosphere on timescales from minutes to weeks (Robock, 2000).
In order to mitigate risks and assess hazards originating from volcanic clouds, accurate observations and forecasts are needed.Advecting volcanic clouds can be tracked by satellite observations, but satellite images in the visible spectrum only result in outer contours of the cloud.Moreover, the global coverage and image frequency of satellite observations is highly inhomogeneous and satellite images only reflect the current state and cannot be used for forecasting.Therefore, numerical models predicting the advection of ash or SO 2 are necessary.
There are several models simulating the advection (and sedimentation) of ash and SO 2 clouds.They are mainly performed on a regular grid and can generally be divided into two types by their numerical framework: Eulerian models like ATHAM (Oberhuber et al., 1998), REMOTE (Regional Model with Tracer Extension) (Langmann, 2000), Fall3d (Folch et al., 2009) or Ash3d (Schwaiger et al., 2012); and Lagrangian models including Puff (Searcy et al., 1998) and NAME III (Jones et al., 2007).Additionally, there are some models using other approaches like semi-analytical tephra transport and dispersion models -HAZMAP (Macedonio et al., 2005;Pfeiffer et al., 2005) or TEPHRA (Bonadonna et al., 2005) -and the Lagrangian-Eulerian model VOL-CALPUFF (Barsotti et al., 2008) and a volcanological adaptation of HYSPLIT (Stein et al., 2015).For more details on these models the reader is referred to a recent review by Folch (2012).
In this article, we extend an existing semi-Lagrangian advection model performed on an adaptive, triangulated mesh (amatos and flash : Behrens, 1996;Behrens et al., 2000) for volcanic emissions.The semi-Lagrangian method has the advantage of a very stable and numerically efficient advection calculation and can be performed in parallel (Behrens, 1996).Adaptive mesh methods have the additional advantage of high resolution in the area where the advected cloud currently resides, while the computational cost is kept relatively low by using a coarse mesh outside the cloud.With this model, simulations forecasting the advection of an ash cloud for several days could be performed at low computational cost (CPU times of less than 1 h) in contrast to comparable simulations on a uniform mesh which required about 10 times more CPU time.We apply our new model to the advection and sedimentation of tephra, because SO 2 clouds cover a much larger area and the sedimentation of tephra occur on timescales of minutes to weeks, while SO 2 and sulfate can remain in the atmosphere for some years.We concen-trate on the advection of the ash cloud, neglecting complex eruption column dynamics and the influence of the eruption column on the surrounding atmosphere.This is a valid assumption, because far enough from the vent these effects play a minor role (Folch, 2012).Furthermore we note that our model neglects processes like ash aggregation as well as wet and dry deposition, as the main focus is on the impact of adaptive meshing as a tool to improve and accelerate ash dispersion forecast models.
In the following we first introduce the implementation of the adaptive semi-Lagrangian advection algorithm and explain how the sedimentation of particles has been implemented into this model.We then turn to the description of our case study of the climactic eruption of Mt.Pinatubo, 1991.
Here we focus first on the main advantages of our solution and then carry out a sensitivity analysis by varying different input parameters.We finish with some discussion (including a detailed performance study) and conclusion.

Model description
In the model we solve the advection equation in conservation form, where C(x, y, z, t) = C(x, t) is a scalar concentration, u = (u x , u y , u z ) is a given wind field and R is a right-hand side which can contain sources and sinks.As usual, we will denote with (x, y, z) = x and t being the spatial and time coordinate, respectively.Using the integral form of this conservation law, d dt where V (t) is a (time-dependent) reference volume and R * is the right-hand side integrated over this volume, we can derive a semi-Lagrangian discretization: Following the description in (Behrens, 2006) the integral expressions can be discretized by a midpoint rule, multiplying the point value of concentration C(x, t) with the area of the dual cell V (x, t) corresponding to grid point x.We obtain where α is the semi-Lagrangian upstream displacement from a given grid point and R * is evaluated at the upstream position.Note that the cell area corresponding to the concentration point changes in time and -depending on the velocity field -may be distorted, such that mass is conserved even in non-divergent flow fields.This scheme deviates from a pointwise semi-Lagrangian method as described in Staniforth and Cote (1990) by the correction factor . It is unconditionally stable as long as V (x, t) does not degenerate.
In order to compute the areas |V (x − α, t − t)| and |V (x, t)| on a nonuniform tetrahedral mesh simple computational geometry methods are used.The upstream value C(x − α, t − t) is interpolated by a linear interpolation within the upstream tetrahedron.While linear interpolation is theoretically very dissipative, it is positivity preserving and monotonous, and in combination with the local mesh refinement described below only small numerical smoothing is observed.
The advection of the cloud is driven by a precalculated wind field u(x, y, z) from the regional scale atmospheric chemistry and climate model REMOTE (for details see Langmann, 2000).Initial meteorological data are taken from the European Centre for Medium-Range Weather Forecasts and boundary conditions are updated every 6 h.The horizontal resolution of the wind field is 0.5 • × 0.5 • (approximately 55 km×55 km).In the vertical direction a sigma-pressure coordinate system subdivides the model atmosphere into 31 layers of increasing thickness between the Earth surface and the 10 hPa pressure level.Here we utilized all vertical layers and interpolated the wind in the x, y and z direction to the grid resolution.This interpolation is linear in space and time to maintain monotonicity and preserve the shape of the vector field.
The semi-Lagrangian method described above employs an adaptive mesh following Behrens (1996): in regions where a high spatial resolution is required, the mesh is refined, whereas the mesh size in other parts of the model domain is kept relatively coarse.Thereby, memory requirements can potentially be decreased by orders of magnitude without losing accuracy.In our case the refinement criterion is based on the concentration gradient ∇C| τ i in a mesh element τ i ∈ T , where the triangulation T denotes the complete set of tetrahedra representing the computational domain.A mesh element is refined if with being the maximum of all local concentration gradients.Accordingly, a mesh element is coarsened if The parameters θ ref and θ crs (with 0 < θ ref ≤ 1 and θ crs < θ ref ) define the relative tolerances for refinement and coarsening, respectively.The a posteriori adaptation strategy computes local gradients in each element after each time step.This computation is cheap and requires just a few operations per element.Those elements fulfilling Eqs. ( 4) or ( 6) are marked for refinement or coarsening, respectively.The mesh is only refined or coarsened, if a certain fraction w ref/crs of grid cells is marked, in order to balance accuracy requirements with computational cost.The cell refinement or coarsening follows a bisection strategy, described in Bänsch (1991).
In our scenario computations we used the following parameters for controlling the mesh: This means that a cell is flagged for refinement if its gradient is larger than 2 % of the maximum local gradient.And the mesh is changed if 0.1 % of the cells are marked.

Particle sedimentation
The sedimentation of tephra from an advecting ash cloud is mainly dependent on the grain size, the density of the particles and the properties (viscosity and density) of the surrounding air.In order to account for sedimentation, the terminal settling velocity v t (balance between drag force and gravitational force) is calculated for atmospheric conditions at every mesh point and every time step.The terminal settling velocity is given by with ρ p being the density of an ash particle with a diameter D p , ρ the density of the surrounding fluid (calculated here from REMOTE simulation results), g the gravitational acceleration and C d the drag coefficient.Settling of particles is then accounted for in the advection equation (Eq. 1) by modifying the vertical advection term as follows: The terminal settling velocity is dependent on the drag coefficient C d which is a function of the Reynolds number (Re) which in turn depends on the settling velocity.Empirical formulations of the drag coefficient for different regimes of the Reynolds number have been suggested by several authors (Seinfeld and Pandis, 2006;Dellino et al., 2005;Bonadonna et al., 1998;Herzog et al., 1998;Ganser, 1993;Arastoopour et al., 1982;Wilson and Huang, 1979).Here we use the model introduced by Ganser (1993), which gave the best results for our conditions (look at the Supplement for more details).
For calculating Re the dynamic viscosity of the air, depending on the air temperature, is required (Pruppacher and Klett, 1997): Particles require a certain time (a so-called relaxation time) to reach their terminal settling velocity (Seinfeld and Pandis, 2006): Here, M p is the mass of the ash particle and C s is a slip correction factor defined by Seinfeld and Pandis (2006).
Using reasonable values for the parameters in Eq. ( 10), it is obvious that the maximum time required by particles to reach their terminal settling velocity is relatively short.Even for particles with a diameter of 1 mm (φ = 0) the maximum relaxation time is only about 9 s.Compared to the default simulation time step of 10 min used here and a simulation period of about 5 days, the relaxation time is negligible.Therefore, we assume that the ash particles are falling directly with their terminal settling velocity.
With the assumption that in dilute clouds ash particles of different sizes do not affect each other but each particle settles individually (i.e., we neglect particle aggregation as well as particle-particle interaction), we apply the model for individual particle diameters and then combine the results of different runs to predict the sedimentation of the complete grain size distribution.Furthermore, as we model the fallout from the umbrella cloud, we assume that the ash particles already reached their maximum injection height at the start of the model.Complex eruption column dynamics are neglected and we suppose no interaction with and no re-entrainment into the eruption column.In addition, the settling of ash particles is strongly affected by rainfall and particle aggregation (see e.g., Brown et al., 2012).Since this work is a first case study of the modeling of sedimentation of ash particles on an adaptive mesh, the impact of rain on the sedimentation and aggregation of ash particles is neglected.Finally, we treat the ash particles as passive tracers and do not monitor the thickness of the ash deposited on the ground.The 1991 Mt. Pinatubo eruption on the Luzon island in the Philippines (see Fig. 1) was one of the largest explosive eruptions in the 20th century.The amount of erupted SO 2 induced a global cooling of at least 0.5 • C in the 2 years following the eruption (Self et al., 1996).Pyroclastic flows, lahars (thick volcanic mudflows) and ash fall made more than 50 000 people homeless, affected the lives of more than a million people and caused 200 to 300 deaths (Punongbayan et al., 1996;Bautista, 1996).
Following a quiet period of about 500 years (Newhall et al., 1996) activity at Mt. Pinatubo started in July 1990 with a surface-wave magnitude 7.8 earthquake along the Philippine fault, about 100 km northeast of the summit (Punongbayan et al., 1991).Many smaller earthquakes were recorded in the following months, and in April 1991 the first eruptions with column heights between 1 and 8 km took place.The explosive phase began on 12 June 1991 and lasted until 16 June with sub-Plinian to Plinian eruptions and column heights between 19 and 40 km.The most violent eruptions occurred between 13:40 and 22:40 PHT (Philippine Time; all following times are given in PHT) on 15 June with more or less continuous high-output activity.The intensity of this eruption period began to decrease after about 3 h at 16:40 on 15 June (Wolfe and Hoblitt, 1996).In the climactic eruption phase lasting 9 h, 80 % of the total erupted volume was ejected and the highest eruption columns were reached (Holasek et al., 1996).
During the first phase of the climactic eruption, the ash expanded radially and formed a huge umbrella cloud.Koyaguchi and Tokuno (1993) analyzed the hourly multispectral images of the Geostationary Meteorological Satellite (GMS) on 15 June 1991 and showed that, following the onset of the climactic eruption at 13:41, the erupted material expanded radially for about 5 to 6 h hours in a giant umbrella cloud.At 14:40 Koyaguchi and Tokuno (1993) identified an ash cloud of 280 km diameter in the satellite images and at 15:40 the umbrella cloud covered an area with a diameter of 400 km.Similar studies using infrared (Lynch and Stephens, 1996) and GMS-4 visible band data were used to determine the radial expansion and advection of the ash cloud in the Nat.Hazards Earth Syst.Sci., 18,[1517][1518][1519][1520][1521][1522][1523][1524][1525][1526][1527][1528][1529][1530][1531][1532][1533][1534]2018 www.nat-hazards-earth-syst-sci.net/18/1517/2018/west-southwest direction.The southwestward advection of the umbrella cloud mainly reflects the wind direction in the stratosphere (Koyaguchi and Tokuno, 1993).Light to moderate tephra was displaced southward and moderate to heavy tephra northeastward by Typhoon Yunya which passed at a distance of about 75 km northeast of the erupting volcano at around 14:00 on 15 June 1991 (Oswalt et al., 1996).This atypical wind in the lower and middle troposphere caused by the passing typhoon lead to the wide distribution of tephra in nearly all directions around the volcano (Wolfe and Hoblitt, 1996).The heaviest tephra falls occurred during the climactic eruption on 15 June, producing tephra fall deposits with up to 33 cm thickness (Paladio-Melosantos et al., 1996).An area of around 7500 km 2 on Luzon was cov- ered by more than 1 cm thick tephra deposits and the entire island obtained at least a trace of ash (Paladio-Melosantos et al., 1996).Paladio-Melosantos et al. (1996) examined the grain sizes of the Pinatubo 1991 tephra fall deposits on the Luzon island relatively close to the vent (≤ 30 km distance), while Wiesner et al. (1995) recorded the fallout of tephra following the climactic eruption by two sediment traps moored at 14.60 • N and 115.10 • at a water depth of 1190 and 3730 m in the South China Sea.

Volcanic ash emissions
In order to properly model the source term R in Eq. ( 8), mass eruption rates need to be defined as model inputs.Our estimate of the mass eruption rates is based on observations by Holasek et al. (1996) in visible and infrared satellite images.For some satellite data Holasek et al. (1996) could not determine the altitude of the eruption plume and we completed values with data from Self et al. (1996).A complete list of eruption heights used to estimate mass eruption rates is given in Table A1.After 16 June 10:41, secondary explosions were induced by the interaction of the hot ignimbrite with water, but these secondary eruption plumes were of lower intensity and are not considered in this study.Holasek et al. (1996) calculated an average eruption rate of 1.4 × 10 5 m 3 s −1 for the climactic phase lasting 9 h from 15 June 13:41 to 22:41.Accordingly, we estimated mass eruption rates for the other eruption phases from the data of Holasek et al. (1996) and Self et al. (1996).The mass eruption rates used in this study are listed in Table 1.Ash is not released evenly along the eruption column into the atmosphere but mostly close to the neutral buoyancy level.Fero et al. (2009) simulated the Pinatubo eruption with different tephra dispersal models and determined that most of the ash was advected in an umbrella cloud at the level of the tropopause at around 17 km -significantly below the maximum column heights of 40 km and below the main trans- port level of SO 2 at around 25 km.We follow their approach and release ash over the time intervals indicated in Table 1 in a 4 km long cylinder with a radius between 2 and 5 • and a medium height of 17 km centered above the Pinatubo volcano (see also Table 2).The size of the cylinder is varied during the sensitivity study as well as its center location in the atmosphere.The mass eruption rates (listed in Table 1) were divided by the cylindrical ash injection volume and multiplied by the simulation time step.Therefore, the particle concentration C at each mesh point and each time step during the release of ash was obtained.Ash settles out of the eruption cloud with a settling velocity depending on the grain size.Summarizing the results of Paladio-Melosantos et al. (1996) and Wiesner et al. (1995) in an area of about 600 km around the volcano, ash fallout ranges from −4 to 9 φ (16 to 0.00195 mm).Since large particles below 1 φ would settle too fast and would require very small time steps and very small particles above 8 φ sediment outside the simulation domain, we neglected these particle sizes.For the sake of simplicity, we used a Gaussian distribution around a mean grain size of 4.5 φ with a standard deviation of σ φ = 2.5 which corresponds well with the estimated bulk mean grain size of 4 φ of the Pinatubo tephra deposit (Fero et al., 2009).The grain size categories used for the es- timation of the settling velocities are listed together with the particle diameter and the particle density in Table 3.

Mesh generation and refinement settings
The initial mesh consists of three cubes, each extending from 4.5 to 24.5 • N and from 0 to 23 km height.In the east-west direction, the cubes stretch from 69.5 to 89.5 • E (cube 1), 89.5 to 109.5 • E (cube 2) and 109.5 to 129.5 • E (cube 3).Each cube contains 6 tetrahedra, so initially there are 18 tetrahedra in the domain.Due to the minimum refinement level of 8 (coarse-mesh level in Table 2) each of these tetrahedra is refined seven times, resulting in a total of 2304 tetrahedra in the initial model domain before local refinement.

E. Gerwing et al.: An adaptive semi-Lagrangian advection model
Figure 5. Ash concentration of the ash cloud on 14 June at 14:00 and on 15 June 22:00 for particles with a grain size of 4.5 φ.A cross section of isosurfaces of the ash concentration is displayed in kg m −3 .The minimum displayed ash concentration is 0.001 kg m −3 .The coarse-mesh level gives the level of global and uniform refinement at the initialization of the grid and the maximum level to which an element is coarsened, whereas the finemesh level defines the maximum level to which an element can be refined.After the initial refinement of the mesh (corresponding to a mesh level of 8), the minimum length of the edges of the tetrahedra in the horizontal direction is 4.5 • .
With the maximum refinement level of 17, minimum edge lengths of 0.6 • in the horizontal direction and 0.7 km in the vertical direction are achieved (compare Table 4).In Fig. 2, the structure of the three-dimensional adaptive tetrahedral mesh of the ash cloud on 15 June 22:00 PHT is displayed.The cloud is seen form the side.The mesh is composed of tetrahedra with variable size.Please note that particle settling is not considered in this simulation.On the right side -where the mesh resolution is higher -the ash emissions are inserted into the model (see above), leading to a larger gradient in the ash concentration which in turn starts the mesh refinement depending on the refinement and coarsening tolerances (see Table 2 and Sect. 2).An animation of the advection of the three-dimensional cloud can be found in the video Supplement (Gerwing, 2017a).
The impact of the maximum level of refinement (finemesh level) is demonstrated in Fig. 3 where a horizontal cross section at the height of the mean transport level (17 km) is shown.Between the ash clouds simulated with a fine-mesh level of 14 (Fig. 3a) and a fine-mesh level of 17 (Fig. 3b), significant differences in shape and ash concentration of the clouds can be observed.When comparing the ash concentration for different mesh resolutions, it is import to consider that the total amount of ash inserted initially is slightly dif-ferent for the different refinement levels, because the discrete volume of the cylinder in which the ash is inserted is dependent on the cell size.Using a fine-mesh level of 14 the discrete volume amounts to 97.6 % of the analytical volume of the cylinder (compare Sect. 3.2.1),while with a fine-mesh level of 17 the inserted ash mass already accounts for 99.5 % of the original ash mass.
For a fine-mesh level larger than 16, we found that the results do not change significantly any more apart from small differences in the ash concentration, but the general behavior of the ash cloud is preserved.Hence, from here on we utilize a fine-mesh level of 17 as maximum refinement level allowing for fast computation of the ash spreading.

Results for the standard model setup
Figure 4 shows the evolution of the settling ash cloud for a particle size of 4.5 φ (0.0469 mm) (animation in the Supplement).At 10:00 on 13 June, the ash cloud is centered above the volcano at a mean height of around 16 km (Fig. 4a).In the following 4 h, the ash cloud settles down to a medium height of around 13 km and is slightly advected to the west (Fig. 4b).On 14 June 14:00, ash particles of the eruption on 13 June have sedimented down to the ground, while the eruption cloud from the second eruption phase starting on 14 June 13:09 (compare Table 1) is still close to its initial position (Fig. 4c).One day later on 15 June at 14:00, shortly after the onset of the climactic eruption, the ash column from the second eruption phase is centered above the South China Sea.While settling, ash particles are advected to the southwest, especially between heights of 8 to 10 km (Fig. 4d).In the following hours, the ash cloud drifts further in the southwest direction (Fig. 4e and f).After the last eruption phase ended (on 16 June 10:41), the ash cloud sinks down and is advected to the west-southwest (Fig. 4g and h).An animation of the settling ash cloud for a particle size of 4.5 φ can be found in the video Supplement (Gerwing, 2017b).
Figure 4 only shows the ash concentration on the surface of the ash cloud.Figure 5 allows a look into the ash cloud where cross sections of the isosurfaces of the ash concentra-  3 and for a simulation neglecting the settling of particles.The concentration on the surface of the ash cloud is displayed in kg m −3 .The minimum displayed ash concentration is 0.001 kg m −3 .tion on 14 June at 14:00 and 15 June at 22:00 are displayed.When ash is inserted, the ash concentration inside the cloud is initially homogeneous (Fig. 5a).While the ash is advected and sediments, ash particles are dispersed and the ash concentration decreases.Since ash is inserted continuously between a height of 14 to 20 km, the highest ash concentrations are obtained in the upper part of the cloud (Fig. 5b).After the onset of the climactic eruption, the ash concentration inside the cloud significantly increases.
The results of the advected and sedimented ash cloud on 15 June at 22:00 are displayed for model runs with different grain sizes and for a simulation without the settling of particles in Fig. 6.In the simulation without the settling of particles (Fig. 6a), most of the ash is advected in the stratosphere www.nat-hazards-earth-syst-sci.net/18/1517/2018/Nat.Hazards Earth Syst.Sci., 18, 1517-1534, 2018 Figure 7. Advected ash cloud on 16 June 04:30 for particle sizes of 4, 4.5 and 5 φ and for a simulation without the settling of particles.In the background, the outer edge of the observed ash cloud is displayed in 3 h intervals (observed contours from Lynch and Stephens, 1996).The outer edge of the simulated ash cloud was defined as the isosurface with an ash concentration of 0.05 kg m −3 , producing the closest fit to the observed data.Note that the observational contour that corresponds to our model time (16 June 04:30) is marked in red.
and upper troposphere in the west-and southwestward direction.The larger the ash particles are, the higher the settling velocity is and the more the ash is advected to the south due to the changing wind patterns in the atmosphere.For particles with grain sizes of 7.5 φ and 6.5 φ, advection dominates and the ash cloud drifts more or less horizontally to the west-southwest (Fig. 6b and c).The effect of sedimentation becomes visible for particle sizes larger or equal to 5.5 φ (Fig. 6d).For simulations with particle sizes between 4.5 and 3.5 φ, a certain amount of the ash particles sedimented to the ground on 15 June 22:00, but advection still had an impact on the motion of the cloud (Fig. 6e and f).For even larger particles, sedimentation dominates the advection of ash particles and the particles settle down in a nearly vertical ash column (Fig. 6g and h).
Projecting the extent of the calculated ash cloud onto satellite observations made during the Pinatubo eruption (see Sect. 3.1) is also quite instructive for determining which particle size or sizes are the best to achieve a good fit with observations.In Fig. 7a, the result of a simulation without the settling of particles is shown.The winds in the lower stratosphere advect most of the ash to the west and only a very small amount of ash travels to the southwest.For particles with a grain size of 5 φ (0.0313 mm), the ash settles slowly to lower altitudes, where the wind is directed southwestwardly.The particles are more or less evenly distributed between west and southwest (Fig. 7b).With increasing particle size, ash is advected in a heart-shaped cloud to the west and the southwest (Fig. 7c and d).The simulated ash cloud with particles of a grain size of 4 φ (0.0625 mm) is of a much smaller  Comparing all simulations with different settling velocities, the best fit to the data of Lynch and Stephens (1996) is obtained by simulations with a particle size of 5 φ.The agreement between the simulated ash cloud and the outline of the umbrella cloud identified by Lynch and Stephens (1996) increases when the results for particle sizes of 4, 4.5 and 5 φ are combined (Fig. 8).An animation of the compared results can be found in the video Supplement (Gerwing, 2017c).On 15 June at 22:30, the simulated ash cloud matches the data of Lynch and Stephens (1996) very well (Fig. 8a).At 04:30 on 16 June, the modeled ash cloud and the observed contour of the ash cloud generally agreed with each other, but the southern extension of the umbrella cloud could not be reproduced (Fig. 8b).

Discussion
In this study we adapted an existing adaptive semi-Lagrangian advection model to model the dispersion of volcanic ash emissions including the sedimentation of particles from the ash cloud.We matched our results with published satellite data of the umbrella cloud during the climactic phase of the Pinatubo eruption and found that the simulation of the advection and the sedimentation matched the observations quite well.The best fits between modeled and observed data were obtained by combining results from simulations with multiple particle sizes (Fig. 8b).In the following we will first test the sensitivity of our results to some of the main input parameters before we turn to a discussion of the performance advantage of our model compared to fixed-grid models.

Sensitivity study
Table 2 lists various input parameters used in the model calculations.Since the initial radial expansion of the Pinatubo cloud is not included in our model, the initial radius of the ash cloud was chosen to be relatively large.In the models of Fero et al. (2009), an initial radius of around 400 km was found to be necessary to account for the radial expansion.We varied the radius between 2 and 5 • (approximately 222 to 555 km).Following Fero et al. (2009), most of the ash was advected at a height of around 17 km, but some amount of the ash was injected at much higher altitudes (Holasek et al., 1996).We therefore tested medium cloud heights between 15 and 21 km.Koyaguchi (1996) reported that the thickness of the umbrella cloud was between 3 and 5 km, while Self et al. (1996) mentioned a cloud thickness of 10 to 15 km.We therefore varied the umbrella cloud thickness between 2 and 8 km.Varying the cloud thickness between 2 and 8 km height did not impact the ash dispersion significantly so this is not discussed in further detail below.In the following sensitivity analysis we use a medium grain size of 4.5 φ (0.0469 mm) in order to not obscure the results due to the differences in settling velocities.
The impact of the initial cloud radius is shown in Fig. 9.The ash cloud was simulated with an initial cloud thickness of 6 km, a medium height of 17 km, a grain size of 4.5 φ, and an initial radius of 2, 3, 4 and 5 • .The area covered by the cloud decreases with a decreasing radius and the shape of the modeled ash cloud becomes more circular for larger radii (compare Fig. 9a and d).At 22:30 on 15 June, the best fit to the outline of the observed ash cloud (identified by Lynch and Stephens, 1996) is obtained with an initial cloud radius of 4 • .
Figure 10 shows the effect of variation in the initial mean cloud height on the extent of the simulated ash cloud.The higher the ash is inserted, the more the cloud is advected to the west by stratospheric winds (Fig. 10d).At lower altitudes, wind mainly carries the ash in the southwest direction and the ash cloud covers a heart-shaped area (Fig. 10a).The modeled ash clouds with an initial height of 17 and 19 km (Fig. 10b  and c) better matches the contour of the ash cloud identified by Lynch and Stephens (1996), but the expansion of the umbrella cloud to the south is not reproduced.
Since none of our simulations, even the best fit one shown in Fig. 9, match the observations completely -in particular in the south -we attribute the remaining differences between the modeled and the observed ash dispersion to the following facts: 1.The initial conditions of the eruption cloud are not well known, including the vertical mass distribution in the plume.
2. The radial expansion of the umbrella cloud is accounted for by a larger initial cloud radius, which might induce deviations in cloud shape.
3. The precalculated wind field might not reproduce correctly the conditions during the eruption, especially the course of Typhoon Yunya.
4. The heavy rain caused by Typhoon Yunya is neglected in these model simulations, so washout and vertical transport will be underestimated.
5. Uncertainties in the outline of the umbrella cloud identified by Lynch and Stephens (1996).Satellite observations only reflect the outer contours of the ash cloud in the uppermost layers.In addition, ash clouds are hard to distinguish from meteorological clouds.
6. Negligence of ash aggregation in the model calculation.

Performance due to adaptive meshing
The main advantage of this model is the use of an adaptive tetrahedral mesh leading to significantly reduced computational costs while tracking the volcanic emissions with a very high local mesh resolution.In order to determine the advantage of our adaptive meshing approach compared to fixedgrid calculations we carried out a series of model runs with a fixed grid (i.e., we set the fine-grid level equal to the coarsegrid level to achieve a fixed grid) and compared them to our standard run with an adaptive mesh.In this case study, all simulations are carried out without the sedimentation of particles.In Fig. 11, results of simulations on a uniform grid with refinement levels of 13, 15 and 17 are compared to a result of a calculation on the adaptive mesh.The shape of the ash cloud is recovered quite well in all calculations, but small patterns in the shape and the ash concentrations are much better recovered with the finer grid structure.As mentioned above (Sect.3.2.2), the initial ash mass varies slightly for different mesh resolutions.In Fig. 12 we compare the computational costs of the different model calculations.The simulation on the adaptive mesh needed only about 50 min, while the calculation on a uniform mesh with a fine-grid level of 17 (i.e., the same maximum local resolution as in our adaptive mesh calculation) required already around 9 h.Those significantly reduced computational cost would allow for ensemble runs with varying meteorological boundary conditions as well as different ash injection assumptions to better constrain and forecast probable dispersion patterns and directions (see e.g., Madankan et al., 2014).

Conclusions
In this study we have demonstrated the versatility of adaptive meshing algorithms for modeling the dispersion of volcanic emissions.More specifically, the high performance of this code would allow, if implemented into operational ash dispersion models, a significant improvement of disper-  Comparison of the computation times of simulations on a uniform mesh with refinement levels of 13, 14, 15, 16 and 17 and on the adaptive mesh with a coarse-mesh level of 8 and a fine-mesh level of 17. Calculations were carried out on a Lenovo ThinkPad with an i3-2310M processor and 8 GB of main memory.sion predictions as model runs could be carried out significantly faster compared to codes using a fixed grid.The research community benefits from such a faster code by being able to resolve the fine filamented structure of volcanic emissions during their transport as well as test more boundary conditions, newly developed sedimentation models (e.g., Bagheri and Bonadonna, 2016) and complex chemical reactions which could occur between different trace gases in the atmosphere while being transported (e.g., Hoshyaripour et al., 2015).
In our sensitivity study we have shown that the initial conditions of the ash cloud significantly influence the region impacted by the ash cloud.Even in cases where meteorological predictions, the initial height, the extent of the ash cloud and the mean grain size of the erupted particles are not very well constrained, our model could be used for forecasting the advection and the sedimentation of ash after a volcanic eruption through ensemble runs and thereby contribute to assessment and mitigation of risks posed by drifting ash clouds.
A further application of the model is for the prediction of the ash loading on the Earth's surface from tephra fallout which only needs an additional two-dimensional array to sum up the deposited ash.
In order to enable the abovementioned application fields, we envision several extensions of our model.First, a more realistic tephra reaction and aggregation model could be implemented.Methodologically, this is a relatively straight forward extension of the right-hand side of Equation (1).Secondly, a multicomponent tephra simulation could be implemented.In order to support adaptive mesh refinement for multiple grain sizes with their own sedimentation rates, a combined criterion would be necessary, which could be an additive combination of the individual component's refinement criteria.Finally, higher-order interpolation for the semi-Lagrangian time stepping could further increase accuracy.

Figure 2 .
Figure2.Triangulated mesh on the surface of the ash cloud.The ash cloud is represented on 15 June at 22:00 in a side view (from south).The concentration on the isosurfaces of the ash cloud is displayed in kg m −3 .The minimum displayed ash concentration is 0.001 kg m −3 .

Figure 3 .
Figure 3. Horizontal cross section of the ash cloud at a height of 17 km, modeled with a fine-mesh level of 14, 17, 20 and 23.Results for 14 June at 14:00 PHT are shown.The colors indicate ash concentration in kg m −3 .The minimum displayed ash concentration is 0.005 kg m −3 .

Figure 4 .
Figure 4. Settling of the ash cloud over time for particles with a grain size of 4.5 φ.The concentration on the surface of the ash cloud is displayed in kg m −3 .Note that different color bars are used.The minimum displayed ash concentration is 0.001 kg m −3 .

Figure 6 .
Figure 6.Modeled ash cloud for the grain size categories listed in Table3and for a simulation neglecting the settling of particles.The concentration on the surface of the ash cloud is displayed in kg m −3 .The minimum displayed ash concentration is 0.001 kg m −3 .

Figure 8 .
Figure 8. Combined results for particle sizes of 4, 4.5 and 5 φ and the outline of the umbrella cloud identified byLynch and Stephens (1996).

Figure 9 .
Figure9.The modeled ash cloud on 15 June 22:30 with an initial radius of 2, 3, 4 and 5 • and the outer edge of the observed umbrella cloud in 3 h intervals (observed outlines fromLynch and Stephens, 1996).The modeled ash cloud is displayed for 15 June 22:30.

Figure 10 .
Figure10.Advected ash cloud on 16 June 04:30 inserted at a mean height of 15, 17, 19 and 21 km in combination with the extent of the observed ash cloud analyzed byLynch and Stephens (1996).

Figure 11 .
Figure11.Cross sections at a height of 18 km on 17 June at 13:20.Comparison between simulations on a uniform mesh with a refinement level of 13 (a), 15 (c) and 17 (b) and a simulation on the adaptive mesh with a coarse-mesh level of 8 and a fine-mesh level of 17 (d).The same figure including the mesh structure can be found in the online Supplement.

Figure
Figure12.Comparison of the computation times of simulations on a uniform mesh with refinement levels of 13, 14, 15, 16 and 17 and on the adaptive mesh with a coarse-mesh level of 8 and a fine-mesh level of 17. Calculations were carried out on a Lenovo ThinkPad with an i3-2310M processor and 8 GB of main memory.

Table 1 .
Eruption phases implemented in the right-hand side.For time periods not listed in this Table, the right-hand side was set to zero.The values given in italics are those for the climactic phase of the eruption.For converting m 3 s −1 to kg s −1 we used an average density of 1500 kg m −3 .

Table 3 .
Macedonio et al. (1988)densities utilized in the sedimentation simulations.The particle densities as a function of the grain size listed here are extracted fromMacedonio et al. (1988).

Table 4 .
Minimum edge length of the tetrahedra in the horizontal and vertical direction.