Tephra4D: A Python-Based Model for High-Resolution Tephra Transport and Deposition Simulations—Applications at Sakurajima Volcano, Japan

: Vulcanian eruptions (short-lived explosions consisting of a rising thermal) occur daily in volcanoes around the world. Such small-scale eruptions represent a challenge in numerical modeling due to local-scale effects, such as the volcano’s topography impact on atmospheric circulation and near-vent plume dynamics, that need to be accounted for. In an effort to improve the applicability of Tephra2, a commonly-used advection-diffusion model, in the case of vulcanian eruptions, a number of key modiﬁcations were carried out: (i) the ability to solve the equations over bending plume, (ii) temporally-evolving three-dimensional meteorological ﬁelds, (iii) the replacement of the particle diameter distribution with observed particle terminal velocity distribution which provides a simple way to account for the settling velocity variation due to particle shape and density. We veriﬁed the advantage of our modiﬁed model (Tephra4D) in the tephra dispersion from vulcanian eruptions by comparing the calculations and disdrometer observations of tephra sedimentation from four eruptions at Sakurajima volcano, Japan. The simulations of the eruptions show that Tephra4D is useful for eruptions in which small-scale movement contributes signiﬁcantly to ash transport mainly due to the consideration for orographic winds in advection.


Introduction
Volcanic eruptions release tephra to the atmosphere causing a hazard that can adversely affect the lives and livelihood of surrounding populations, leading to damages in sectors such as agriculture, forestry, fisheries. Furthermore, airborne tephra can be widely dispersed and have a serious impact on aircraft operations, as shown in the 2010 eruption of Eyjafjallajökull volcano [1]. For these hazard assessments, numerical simulations using advection-diffusion models are used to predict the dispersion of tephra in the atmosphere and the mass of tephra deposits [2].
After an eruption and the formation of a plume, tephra is advected by the wind, diffused due to turbulent processes, and descends due to its mass. In order to simulate the transport and deposition of tephra, numerical models require input parameters such as the plume height, total mass of tephra, wind field, and particle size, all of which are either based on observations, estimated, or simulated. Based on the meteorological conditions and the size of the eruptions, tephra transport and deposition can be simulated on scales of tens to hundreds of kilometers (e.g., [3][4][5]).
Tephra deposition prediction models that can reproduce small-scale eruptions are necessary for the study of volcanic plume dynamics. Up until now, most studies have focused on large-scale eruptions such as Plinian eruptions, but they have not been able to establish parameters such as the distribution of the segregation height of tephra in the plume (tephra segregation profile (TSP); e.g., [5,6]) and the interaction between particles and between particles and fluid (e.g., [7,8]). This is partly because large eruptions are infrequent, leading to the shortage of usable observations. Vulcanian eruptions present an opportunity to cover this lack of data due to their high frequency and relative safety of data acquisition.
For tephra transport numerical modeling to be suitable for vulcanian eruptions, it is first necessary to improve the consideration of the advection process. When tephra segregates from the plume at a low altitude from the surface of a volcano, it is important to take into account the bending of the ascending plume [9] and the heterogeneity of the wind field caused by the presence of the mountain (e.g., [10,11]).
When initializing a tephra transport simulation, it is common to assign grain size distribution to characterize the range of released particles. However, laboratory experiments have shown that the settling velocity of a particle of the same size depends on its effective density, shape, and whether or not it is bound to other particles (such as [6,12]), suggesting the importance of constraining the descent process not by the particle size distribution but by the settling velocity distribution. Due to the limitation of the observation, the particle size distribution has been assigned as the input parameter of the previous tephra deposit prediction model, and the settling velocity distribution was approximated from the particle size distribution based on the simple parameterizations. Particles with irregular effective density were calculated with an option (e.g., [13]) or neglected in these models. If the model handles a settling velocity distribution as an input parameter, such considerations can be made by simply changing the mass distribution without any special options.
In recent years, optical disdrometers have been increasingly used to observe tephra sedimentation [8,14,15] because they can measure the settling velocity of falling particles with a high temporal resolution. This allows for an alternative strategy: the characterization of the released tephra through a settling velocity distribution.
As such, in this study, we develop an advection-diffusion model Tephra4D based on the advection-diffusion model Tephra2 [3]. Tephra4D can take into account the spatial distribution of the time series of three-component wind velocity field and atmospheric density in the advection process and constrain the descent process by the settling velocity. Using disdrometer data from four eruptions from Sakurajima volcano, Japan, two sets of simulations are carried out. A simulation with the simple wind field is performed with a reimplemented model of Tephra2 constrained by the settling velocity distribution and simulations using the modified Tephra4D model. By comparing the results with the observations using a disdrometer network, we examine how the specification changes made in Tephra4D contribute to the improvement of the calculation and which points need further improvement.
The paper is structured as follows. Initially, the detailed calculation of the tephra deposit by Tephra4D is described in Section 2. Then the eruptions studied, model setup, and the results of the calculation are detailed in Section 3. In Section 4, we discuss the advantages of Tephra4D and important factors. Finally, Section 5 summarizes the main conclusions of the study.

The Advection-Diffusion Model Tephra2 and Its Modifications
Tephra2 is an advection-diffusion model improved from the analytical model of [6], which numerically solves the equation of terminal velocity and dispersion [16]. The migration path of the center of gravity of the particle cluster as they segregate from the plume and fall into the atmosphere (hereinafter referred to as the "trajectory") is calculated and then the tephra deposit load is obtained considering the diffusion process ( Figure 1). Tephra2 employs a horizontally homogeneous steady state, does not account for plume bending, and particles are classified by diameter. These aspects of the model need to be improved in order to successfully model vulcanian eruptions. We do this by employing temporally evolving wind fields obtained from the Weather and Research and Forecasting (WRF) model, including a mechanism for plume bending, and classifying the particles by terminal velocities instead of diameter. We refer to this improved model as Tephra4D.
A model that implements plume bending in Tephra2 was developed by [9], while in this study we implement another simple assumption. The vertical velocity when the plume rises is calculated in the same way as the particle falls. The ascending velocity and the horizontal coordinate of the segregation point are independent of the settling velocity.
Tephra4D takes into account wind heterogeneity in advection, but not in diffusion, and follows the assumption in Tephra2 that particles diffuse concentrically only in the horizontal direction. This is because the horizontal size of the plume is about 1.5 km, even for the largest plume in this study, and the horizontal wind is approximately uniform within this region. We reimplemented Tephra2 to use arbitrary TSP and settling velocity distributions as input parameters for control experiments. This reimplemented version is referred to as Tephra2 PY in this study. The difference of specifications among Tephra2, Tephra2 PY , Tephra4D are described in Table 1. We prepare a representative particle corresponding to each settling velocity class. The movement of the center of gravity of a group of particles is calculated as the movement of the representative particle. The effective density of the particle is assumed to be 2640 kg/m 3 [18]. The diameter and shape parameter of the particle are adjusted to take the same settling velocity as the median of the settling velocity class of the disdrometer at 0 m asl.
Let L A (h seg , d) be the areal concentration of particles of settling velocity v t arriving at point A on the ground from a plume at the segregation height h seg , then the tephra deposit load at point A S A can be expressed as the sum of the areal concentrations of particles as a function of the segregation height and terminal velocity: Assuming that a falling particle cluster diffuses only horizontally and that the diffusion equation is expressed by a two-dimensional Gaussian distribution (variance σ) of the horizontal deviation (x, y) from the center of gravity of the particle cluster, the areal concentration L A (h seg , v t ) at point A, (x A , y A ) away from the center when the center reaches an altitude of point A, is obtained from the following equation: where M is the mass of the particle cluster with settling velocity v t segregated from the altitude h seg and σ A is the dispersion of the cluster at the altitude at which point A is located. Dispersion σ 2 is expressed as a function of dispersion time t and empirically given by an expression proportional to the first power of t when t does not exceed a threshold FTT (Fall Time Threshold) or an expression proportional to the 2.5 power of t when t exceeds FTT [6]. In this study, it is assumed that the horizontal spread of a plume is also caused by the diffusion of pyroclastic particles when the plume rises blowing in the wind. Assuming t = t A (the time between the departure at the vent and the arrival at point A), σ 2 is calculated as follows: where K and C are the diffusion coefficient. FTT is defined to be 3600 s, the setting value in Tephra2, and considering the continuity of σ at t = FTT, which is not considered in Tephra2 and considered in Tephra2 PY , C is defined as follows: The relationship between the radius of the plume r and the elevation from the vent h is empirically given as follows [19]: and in Tephra2 and Tephra2 PY the following relationship is also given [3]: This relationship is also followed in Tephra4D. Using Equations (3)-(6), the time between the departure at the vent and the segregation at altitude h seg , t up (h seg ), is obtained as follows: Applying Equations (3), (5) and (6) to a vulcanian eruption at Sakurajima, K was estimated from the time evolution of the plume top height. The height was read from the live camera image [20] every minute as shown in Figure 2. The optimal value is K = 100 m 2 /s although there is a large variation. The vertical velocity when particles rise through the plume is calculated based on Equation (7). We do not take into account the temporal variation of the mass eruption rate.

The Trajectory Calculation
In the trajectory calculation, we consider a Cartesian coordinate system with the eastward, northward, and upward as the positive directions of x-, y-, and z-axis, respectively. Since the position of the meteorological field data is based on the UTM coordinate system in the horizontal direction and the altitude asl in the vertical direction, as shown in Figure 3, the particle exists in a hexahedron. The meteorological field inside this hexahedron is assumed to be a spatially homogeneous field represented by the value given to point C in Table 2. The entire meteorological field will be composed of a large number of hexahedra, each with a homogeneous wind field and atmospheric density.  Table 2. The symbols to be used in trajectory calculation. The index n indicates nth step.

Symbol
Definition The lower and the higher surface of the hexahedron.
(x n , y n , z n ) The coordinate of a particle inside a hexahedron.
(u n , v n , w n ) The wind vector in the hexahedron.

v tn
The terminal velocity of the particle.
t n Time to start. ∆t x , ∆t y , ∆t z Time for the particle to reach the interface almost vertical to the x-, y-, and z-axis.
In each step, after the particle starts to move, the particle goes straight in the hexahedron, reaches the interface of the hexahedron, and moves to the next step in the neighboring hexahedron. The movement of the particle at nth step is expressed as the following recursion formula: Each symbol is described in Table 2. The recursion formula about time is as follows: The same calculation is performed inside the hexahedron which the particle entered. If the particle moves to a neighboring hexahedron and then immediately returns to the original hexahedron, i.e., the hexahedron on which the particle resides at a given step is the same as the two steps earlier, the coordinate shift is calculated with the velocity component in the direction of boundary travel set to zero and the other velocity components set to the average of the velocity fields on the two hexahedra. Such a sequential calculation of trajectory is repeated until the particle reaches the ground or the horizontal limit of the calculation domain every indefinite spatial step. The calculation of ∆t x , ∆t y , ∆t z is detailed in Supplementary Material S1.

Wind Field and Atmospheric Density
Tephra4D uses temporally evolving three-component wind and atmospheric density data. The data were computed using the WRF model [21] with a horizontal grid spacing of 300 m, 58 layers in the vertical, and a temporal output of 10 min. The results from WRF are interpolated to a vertical grid to be used by Tephra4D. The detailed calculation of the wind field by WRF is described in Appendix C.
In Tephra2 and Tephra2 PY , the wind field is assumed to be horizontally uniform and does not have a vertical component. The wind field above the crater is interpolated from the wind speed and direction given by the user for each section. The wind speed at sea level is assumed to be zero and the wind speed v(h) at an altitude h asl under the crater is assumed as follows: where h vent is the altitude of the vent asl. Assuming the density of fluid at 0 m asl as ρ a0 , the density of fluid ρ a at the altitude h is calculated following the hypsometric formula: where the unit of h and 8200 is m. The simulated wind field for the eruption on 16 June 2018 is shown as an illustrative example. Tephra transport occurred to the west of the vent, agreeing with the resolved wind direction (Figure 4a). The atmospheric flow over the volcano is dominated by orographic waves, indicated by the oscillation seen in the west side (leeside) of the volcano (repeated pattern of downwards and upwards wind extending outwards) [22]. Such orographic activity has been shown to play an important role in the transport of tephra over complex topography [10].

The Settling Velocity of Pyroclastic Particles
We consider that the difference in settling velocity for the same diameter corresponds to the difference in effective density due to the shape of the particles and the degree of agglomeration. We first obtain the relationship among the effective density of the pyroclastic particle ρ p (kg/m 3 ), the terminal velocity v t (m/s), and particle diameter d (mm). Assuming that the shape of the pyroclastic particle is a spheroid, the relationship between ρ p , v t and d is represented by the following set of equations connected by the drag coefficient C D and the Reynolds number R a [6]: where g is the gravitational acceleration, η a , ρ a are the viscosity and density of the surrounding fluid (i.e., the atmosphere in this study), and F is the shape parameter of the particles. Substituting Equations (13) and (14) into Equation (12) we get: Based on the diagram of the diameter and shape parameter F of the pyroclastic particles at Stromboli volcano [23], F is calculated using the following equation: Assuming that v t and d are equal to the settling velocity and particle size observed by the disdrometer, respectively, η a , ρ a are the values for the atmosphere at 20 • C, 1.8 × 10 −5 Pa·s, 1.205 kg/m 3 , the settling velocity, effective particle density is calculated as shown in Figure 5. Tephra2 and Tephra2 PY calculates the terminal velocity from [17], while Tephra4D calculates it from Equation (15), which is used to calculate the observed tephra deposit load in Appendix B. As shown in Supplementary Material S2, the distance to reach the terminal velocity is negligible concerning the height of the plume top, so the particles are considered to have reached their terminal velocity when they are segregated from the plume. The spatial distribution of the terminal velocity of a pyroclastic particle with a terminal velocity of 2.2 m/s on the ground is shown in Figure 6a for the eruption on 16 June 2018. The terminal velocity decreased by about 20% between 5 and 0 km asl. Due to the influence of downward orographic wind to 3 km asl (Figure 6b), at the point where the downwind is strongest, the particle falls with 8.0 m/s in the settling velocity, which is the sum of the terminal velocity of 2.3 m/s plus the downward wind of 5.7 m/s. Thus, when strong mountain waves are formed, the change in settling velocity is more strongly affected by the change in the downward wind than by the change in atmospheric density. Solid gray lines are topography used to calculate trajectories, and dashed gray lines are topography used to calculate tephra deposit distribution. The vertical direction is emphasized twice as much as the horizontal direction. Shading indicates the total settling velocity of a particle (red upwards and blue downwards motion) in (b).

Tephra Segregation Profile
Since TSP M(h seg , v t ), which is substituted in Equation (2), is highly arbitrary, it was constrained using the model of the distribution function. At an altitude interval of ∆h, the mass of tephra with terminal velocity v t at altitude h seg discretized between the vent altitude h vent to the plume altitude h vent + h plume is given by the following equation defining the function of TSP as m(h): where M total (v t ) is the settling velocity distribution of the total mass fraction of tephra. In this study, the observed settling velocity distribution in each eruption is applied to M(v t ). We applied three distribution functions m(h) in this study.
1. Logarithmic Gauss distribution with top concentration;

Logarithmic Gauss distribution with bottom concentration;
3. Uniform distribution.
These distributions are shown in Figure 7. TSP 1 and 2 have m(h) > 0 below the vent and above the plume top, respectively, and h seg is restricted to the height from the vent to the plume top. As such, the sum of M in Equation (17) is smaller than M total .

The Eruptions Studied
In this study, we used eruptions from Sakurajima volcano, one of the most active and closely monitored volcanoes in Japan [24] to evaluate Tephra2 PY and Tephra4D. The eruptive activity of Sakurajima volcano is detailed in Appendix A. We analyzed four eruptions shown in Table 3 that occurred at the Minamidake crater, one large eruption (L), two moderate eruptions (M1, M2), and one small eruption (S). A disdrometer network was installed on the island of Sakurajima to observe the tephra deposit load for error evaluation. The detailed observation is described in Appendix B. The height and direction of the plume tops are announced by JMA (eruptions L and M1), and in case the plume was obscured by clouds, they were estimated using the X-band MP radar image (eruptions M2 and S) [24,25]. Eruption volumes were calculated based on [26] from the strain changes associated with the eruption at Arimura (AVOT in Figure A4). The range of the plume height covers the typical height in Sakurajima ( Figure A2a). The eruptions were chosen as they represent a range of typical eruptions from Sakurajima that led to tephra sedimentation over the volcano at a large number of disdrometers (≥4). In the two eruptions where the direction of the flow was above the crater (T in the table), tephra deposit was detected at stations from north to west as seen from the crater.
The temporal variation of the tephra deposit load at six sites associated with this eruption is shown in Figure 8a. Since vulcanian eruptions and subsequent continuous tephra emissions have different parameters such as plume height and TSP, we extract only the mass of pyroclastic particles ejected immediately after the eruption. To do this we consider that sedimentation from the initial explosion ends if there is a 10-min period of nondetection for each disdrometer location. For the analyzed sample, the settling velocity distributions of the tephra deposit load for the four eruptions were calculated according to Appendix B (Figure 8b) and applied to M total (v t ) in Equation (17). The detected settling velocity range was widest and the heaviest tephra deposit was detected during eruption M2. The spatial distributions are shown in Figure 9. The azimuths of tephra deposit were distributed with the range of about 90 • in all the eruptions. Eruptions L to S produced the largest tephra deposit in the west, northwest, west-northwest, and north, respectively. More than 1 kg/m 2 of tephra deposit was detected at a total of three sites in eruptions L and M2.  Table 3.

Model Setup
In this study, the settings in Table 4 were applied. ∆h is the vertical resolution of the wind field and the slice of plume model, h vent is the altitude of the vent, ∆t is the temporal resolution of the wind field, ρ a (h) is the atmospheric density at the altitude h, K and C are the diffusion coefficients in Equation (3), and η a is atmospheric viscosity.

Results Comparison
The trajectories of tephra calculated in Tephra2 PY and Tephra4D are shown in Figure 10. The horizontal position of the segregation point at the top of the plume was calculated almost the same as the horizontal coordinate of the vent in eruptions M1, M2, S, and in the same direction where the deposit was detected in eruption L. The maximum horizontal travel distance during ascent was calculated to be about 1.5 km, and the consideration of plume bending in this study does not seem to have a significant impact. Trajectories in Tephra4D change direction more frequently than those in Tephra2 PY , reflecting the consideration of horizontal heterogeneity in wind fields. The point reached on the ground is generally nearer in Tephra4D than in Tephra2 PY when particles segregate from low altitude and farther when particles segregate from high altitude. This is because the downward wind blows near the ground and the upward wind blows far from the ground in the region just leeward from the vent in the cases in this study (Figure 4c). Based on these trajectories, the tephra deposit load distributions are calculated as shown in Figure 11. Overall the change in the way the wind field is used in the cal-culations can have a considerable impact. As with the trajectories, the most significant impact is seen close to the vent as enhanced deposition leads to larger deposit loads in the Tephra4D calculations. We compared the calculation accuracy of Tephra2 PY and Tephra4D by the root mean square error (RMSE) and mean absolute percent error (MAPE) between calculated tephra deposit load S cal and observed tephra deposit load by disdrometer S obs . RMSE and MAPE are calculated using the following equations, respectively: S cal (er, v t , TSP, n) − S obs (er, v t , TSP, n) S obs (er, v t , TSP, n) where er, n are eruption and observation site, respectively. RMSE is shown in units of kg/m 2 while MAPE at a percentage. The sites where S obs or S cal is bigger than 0 are considered in Equation (21) and the sites where S obs is bigger than 0 are considered in Equation (22). N is the number of considered sites. RMSE and MAPE are classified by the settling velocity class. The number of combinations of eruption, TSP, and settling velocity class where the number of sites where S cal > 0 and S obs > 0 N', are shown in Table 5. Based on the disdrometer observations, calculations using Tephra4D consistently improve results for particles with settling velocities less than 0.8 m/s but accuracy significantly drops for larger particles, particularly for settling velocities over 7.2 m/s. The settling velocity groups where more than half of the combinations where N' > 0 are 0.8-2.4 and 2.4-7.2 m/s. RMSE and MAPE distributions in these groups are shown in Figure 12. The ranges of RMSE are generally similar with a small but notable reduction in Tephra4D, while MAPE is reduced in Tephra4D especially with the settling velocity of 2.4-7.2 m/s. RMSE and MAPE distribution for the individual eruption is shown in Figure 13. Both RMSE and MAPE distributions improved significantly for the Tephra4D simulation in eruption L. In eruption M1, both the maximum and minimum of RMSE and MAPE reduced but the median increased with settling velocity 0.8-2.4 m/s. Accuracy was decreased based on both RMSE and MAPE and for both settling velocity groups in eruption M2. In eruption S, the maximum of RMSE and MAPE decreased but the minimum increased with settling velocity of 0.8-2.4 m/s.

Discussion
We estimated tephra deposits in four eruptions with various plume heights and wind fields in Tephra4D and Tephra2 PY . The tephra deposit load calculation is improved by Tephra4D in eruptions L and S compared with Tephra2 PY , were mostly unaffected in eruption M1, and worsened in eruption M2. The order of plume height is L, M1, M2, and S (M1 and M2 are the same), and there is no correlation between plume height and Tephra4D's advantage. Since L was the eruption with the most dynamically moving trajectories during ascent and descent, and S was the eruption with the smallest plume height, Tephra4D is suggested to be useful for eruptions in which small-scale movement contributes significantly to ash transport.
The most significant reason for this advantage of Tephra4D is suggested to be the orographic wind consideration in the advection process, whose importance in tephra dispersion has been suggested in [10,27]. As shown in Figure 4b, strong downslope winds are present within 3 km of the crater, and the tephra deposit load near the crater was heavier in Tephra4D with vertical winds than in Tephra2 PY without vertical winds. Such difference in consideration into orographic winds affects the settling velocity range of tephra to deposit in the calculation. The settling velocity range where N'>0 is shown in Table 6. The lower limits of the velocity were more accurate in Tephra4D, while the upper limits of the settling velocity were more accurate in Tephra2 PY , suggesting that Tephra4D is useful for eruptions in which many particles with low settling velocity are ejected. The wind in the cases M1 and M2 was relatively weaker than that in the cases L and S leading to weaker orographic wave activity (Supplementary material S3). This can explain the lack of improvement for these two cases. Despite the reproducibility of orographic winds in Tephra4D, Tephra2PY sometimes produces calculations more consistent with observations. This may be because Tephra4D does not take into account the orographic flows in the diffusion process. The movement of the particle cluster is represented by a single trajectory for each settling velocity and segregation height, and the diffusion of the particle cluster is considered to be concentric without considering the heterogeneity of the wind field. The wind heterogeneity in the cases L and S might be strong enough that the advantage of considering it in advection outweighed the disadvantage of not considering it in diffusion. Implementing the effect of orographic wind on diffusion, Tephra4D will contribute to the estimation of unknown parameters such as TSP and the interaction between particles and between particles and fluid.
The accuracy of the calculation of tephra dispersion affects the diffusion coefficient. When the diffusion coefficient is small, the accuracy of the calculation of the tephra deposit distribution is greatly affected by the accuracy of the particle movement estimation, and it is difficult to perform a highly accurate calculation in substituting a small diffusion coefficient. In this study, we applied 100 m 2 /s as K from the time evolution of the plume top height and the tangential width of tephra deposit are properly reproduced. In previous studies, K was estimated as shown in Table 7 and is longer than K in the present study. In the studies for the eruptions from Ruapehu [28] and Etna [29], diffusion during the plume rising was neglected and K on descent might be overestimated, but K in Tephra4D was small even taking it into account. Since these previous studies gave a horizontally homogeneous wind field, the diffusion coefficient likely compensated for the three-dimensional heterogeneity of the wind field. In Tephra4D, the heterogeneity in advection and the mountain wind contributed to the decrease in the diffusion coefficient compared to previous studies. By considering the heterogeneity in the diffusion as well as in the advection, we might calculate more accurate dispersion and reproduce the local variation of the tephra deposit load distribution more accurately.

Conclusions
The study of tephra dispersion from small-scale eruptions is required because of its high frequency. We created a new advection-diffusion model modified from Tephra2 PY in order to study small-scale eruptions, which need to consider local deviations of wind fields. Two advection-diffusion models with different wind fields and ascending plume trajectories were used to calculate the tephra deposit load distribution for the vulcanian eruptions from Sakurajima volcano, Japan.
Results from Tephra4D, the model developed in this study, show a shift in the deposition of tephra closer to the vent compared to the conventional model and enhanced lateral diffusion of the tephra cloud leading to wider deposition patterns. The model is expected to be useful for eruptions in which small-scale movement contributes significantly to ash transport. The main factor that led to these differences was the consideration of orographic winds. The comparison against results using Tephra2 PY shows reduced errors by considering the heterogeneity of the wind field, potentially leading to improvement in our study to establish unknown parameters in plume dynamics. 2006 and 2016 [11]. Minamidake is further divided into two nearby craters (A and B) and has been bearing the main activity since October 2017. The majority of explosive activity from the volcano since 1955 has been ash-rich vulcanian eruptions, occurring after increasing pressure causes the brittle plug over the vent to destabilize [30]. The duration of such eruptions is commonly limited to a few minutes, up to an hour [11,31]. Plumes from Sakurajima can reach up to 10,500 m above ground level (agl) in height [32]; however, plume heights are commonly constrained within 500-5000 m [11] with the number of eruptions decreasing exponentially against plume height ( Figure A2a). The cumulative spatial distribution of the total deposit load from 2009 to 2019 around the volcano is shown in Figure A2b based on data collected by Kagoshima prefectural government [33]. The accumulated mass ranges between 30 kg/m 2 near the northwestern coast and up to 250 kg/m 2 to the south of the vent.

Appendix B Tephra Load Observation Using the Disdrometer
A disdrometer is an instrument used to measure the particle size and settling velocity of each particle that passes through the laser beam irradiated area. Here, tephra deposit load was measured using a network of one-dimensional disdrometer Parsivel 2 (OTT; Figure A3), which can automatically measure the diameter and settling velocity of pyroclastic particles with a temporal resolution sufficient to measure a time series of tephra deposit load (e.g., [31]). The disdrometer emits a laser beam with a wavelength of 650 nm at an area 30 mm wide from the transmitter and detects the range and duration of the laser beam drop when the fallen particle blocks the laser beam at the receiver unit 180 mm away from the transmitter unit as the voltage change, and then measures the size and settling velocity of each particle. The number of particles in an effective total of 960 classes is recorded by combining the particle size (30 classes from 0.25 to 26 mm) and settling velocity (32 classes from 0 to 22.4 m/s). To detect tephra deposits in any direction, disdrometers were installed at 17 points in all directions from the crater of Minamidake ( Figure A4). The disdrometer can record the number of particles detected at user-selected intervals from the previous recording time, and we set the recording interval to 1 min to take into account the duration of the vulcanian eruption, the accuracy of the arrival time calculated by the advection-diffusion model, and the data volume.
To estimate the deposit load from the disdrometer observations the following empirical equation is used: where ρ ij is the effective density calculated as shown in Figure 5. Note that ρ ij ' in Equation (A2) is the lower limit of the effective density for each particle size and settling velocity class. The correlation between the tephra deposit load calculated by Equation (A1) and the load determined from the tephra deposit load in 59 events collected during the same period, from May 2017 to October 2019, is shown in Figure A5. The calculated values are estimated to be about 0.1-2.3 times larger than the measured values. Figure A4. The disdrometer station network. Red circles are the disdrometer stations and the gray circle AVOT is the strain station. Figure A5. Correlation between tephra deposit load obtained from deposit sampling and the disdrometer observations and Equation (A2).

Appendix C WRF Setup
Version 4 of the Weather Research and Forecasting (WRF) model [21] was used to provide the initialization meteorological data for Tephra4D. Although Tephra4D calculations are carried out in an offline manner, a high temporal output in WRF was used to capture the resolved variability of the wind field, as this has been shown to lead to an improvement in tephra transport modeling [34].
The simulations were initialized using the Medium-Range Weather Forecasts (ECMWF) reanalysis dataset (ERA5 [35]) (∆x 31 km, 137 vertical levels, hourly output [35]). In total three one-way nested domains were used with horizontal grid spacing decreasing from 7.5 km (151 × 151 grid points; ∆t = 30 s) to 1.5 km (301 × 301 grid points; ∆t = 6 s) and finally to 0.3 km (376 × 426 grid points; ∆t = 0.167 s) for the innermost domain ( Figure A6). The same grid spacing is used for the Tephra4D simulations. The Shin-Hong scale-aware Planetary Boundary Layer (PBL) scheme [36] was used to parametrize turbulent motion in the innermost domain as at ∆x = 300 m it falls within the "turbulence gray zone" [37]. Other than the PBL scheme, a full physics-parameterization set was used: hydrometeor microphysics [38], long-and short-wave radiation [39], a surface layer scheme [40], and a land surface model [41]. The simulations were initialized 18 h before each eruption to allow for model spin up time. In the vertical, 58 vertical levels were used, with vertical grid spacing at 50 m from the surface up to 1 km asl, increasing to a spacing of 1 km near the top of the domain at 50 hPa (≈21 km). Figure A6. Placement of the WRF domains.