Tephra segregation profiles based on disdrometer observations and tephra dispersal modeling: Vulcanian eruptions of Sakurajima volcano, Japan

The profile of tephra concentration along a volcanic plume (i.e., the tephra segregation profile) is an important source parameter for the simulation of tephra transport and deposition and thus for the tephra sedimentation load. The most commonly‑used approach is to treat an eruption as a single event (i.e., with a time‑averaged mass erup‑ tion rate; MER). In this case, it is common to use pre‑determined profiles that feature most of the tephra segregate at the top of the plume. However, case studies based on observations have revealed that large concentration maxima also appear at the lower part of the plume. To investigate this discrepancy, the impact of plume height on the tem‑ poral variations in the MER is examined. To this end, we use the tephra transport and dispersion model Tephra4D with MER estimates obtained from geophysical monitoring and maximum plume height observations to calculate the spatial distribution of the tephra deposit load for 39 eruptive events that consisted of explosions and quasi‑ steady particle emission from the Sakurajima volcano, Japan. A comparison of the model results with observa‑ tions from a disdrometer network revealed that for both kinds of activity, maxima in tephra segregation can occur at heights below the reported plume height. The tephra segregation profiles of Vulcanian eruptions at Sakurajima volcano are consistent with most of the modeling studies giving profiles that feature most of the tephra segregating at the top of the plume if the temporal variation of the MER is taken into consideration to properly represent the total series of eruptive events in a sequence. This highlights that even though the activity at Sakurajima volcano is com‑ monly characterized simply as Vulcanian eruptions, in addition to the primary plume developed due to the initial instantaneous release caused by the explosion, the subsequent continuous plume that can accompany the eruption plays an important role in particle emission. Calculations could not reproduce the simultaneous deposition of parti‑ cles with a wide range of settling velocities in observations, suggesting the importance of volcanic ash fingers caused by gravitational instability in tephra transport simulations.


Introduction
The precise modeling of an eruptive column is important for simulating volcanic eruptions as it controls the way tephra is introduced into the atmosphere in a tephra transport and dispersion model (TTDM).The eruptive column is typically described by a few key parameters, such as the mass eruption rate (MER), total mass of tephra particles, column height, and total grain size distribution.Another important consideration is the vertical distribution of tephra along the column profile (Scollo et al. 2008), which we refer to here as the tephra segregation profile (TSP).
The TSP function of the plume height proposed by Suzuki (1983), hereafter referred to as the Suzuki TSP, is used as a source parameter in many TTDMs (e.g., Macedonio et al. 2005: HAZMAP;Schwaiger et al. 2012: Ash3D;Costa et al. 2006;Folch et al. 2020: FALL3D;Stein et al. 2015: HYSPRIT; Shimbori and Ishii 2021: JTA-ATM).This function is used to obtain the relation between the altitude above vent level (hereafter "avl") h seg and the mass M of the segregating particles based on a model in which the vertical velocity of the plume is a maximum at the vent and zero at the top (Suzuki 1983): where k is the deviation parameter for altitude, M tot is the total mass of the tephra, and h p is the plume height.A smaller k value leads to a larger mass fraction of particles that segregate at the lower part of the plume.The maximum segregation height can thus be controlled by fine tuning the k parameter (Fig. 1).Most of the previous studies based on the TTDM calculations, k = 4 (Hurst et al. 2017;Poulidis et al. 2017;Trancoso et al. 2022) or 8 (Schwaiger et al. 2012) are assigned so that most tephra particles segregate at the top of the plume.Pfeiffer et al. (2005) suggested a modified function from Eq. ( 1) and Folch et al. (2020) and Poulidis and Iguchi (2021) adopted the function in their modeling.The deviation parameters used in the function are also tuned so that most tephra particles segregate at the top of the plume.
Despite the ubiquitous use of top-loaded TSPs in modeling studies, a number of case studies based on the range of eruptive conditions have highlighted exceptions, implying that a large proportion of particles segregate at lower altitudes than the maximum column height, in an apparent departure from the top-concentrated TSPs applied in TTDMs using Suzuki TSP.Mannen (2014) calculated the tephra deposit load using the TTDM Tephra2 (Bonadonna et al. 2005) and estimated the TSP for the Izu-Oshima eruption in 1986, representing the observed spatial distribution of the tephra deposit load via a grid search method.It was found that most of the particles segregate at an altitude of 1-3 km, which is significantly lower than the maximum column height (10 km).Similar results have been reported for larger eruptions.Fero et al. (2008), studying the Mt.St. Helens 1980 eruption, andFero et al. (2009) and Cao et al. (2021), studying the Mt.Pinatubo 1991 eruption, showed that using an eruptive plume model to initialize tephra in a TTDM can lead to an overestimation of the main dispersal height as compared to satellite-derived data.Their results suggest that volcanic ash particles mainly segregated at an altitude lower than that of the umbrella cloud, a significant departure from the typical TSP distribution.
Tephra segregation maxima can appear below the reported plume height when the plume shape deviates from the ubiquitous mushroom shape (e.g., Fero et al. 2009) or when there are multiple maxima along the profile.When the resulting plume is modeled, the plume shape can be directly controlled by the k value of the Suzuki TSP, which is typically assigned to a specific eruption (Pfeiffer et al. 2005;Poulidis et al. 2017;Cao et al. 2021;Poulidis and Iguchi 2021) and used to describe the plume for the entire duration.The existence of secondary tephra segregation maxima is thought to arise from differences in the eruptive style and/or resulting plume heights throughout an eruptive event, suggesting the need for a time series of TSPs, which is important for large volcanic eruptions (e.g., Bruckert et al. 2022).
The aim of the present study is to construct a realistic time series of TSPs that represent the different phases (i.e.Vulcanian eruptions and continuous emission) of the eruption.We examine two factors, namely the sensitivity of the simulation results to the k value and the impact of the temporal evolution of the TSP.To create a realistic time series of TSPs that represent the different phases of an eruptive event at the Sakurajima volcano, a detailed evolution of the MER during an eruptive event is first estimated using geophysical observations (ground deformation and seismograph observations; Iguchi 2016).This method has been extensively used in case studies of eruptions at the Sakurajima volcano (Iguchi et al. 2019(Iguchi et al. , 2022)).The estimated MER is used to initialize a number of simulations using the TTDM Tephra4D (Takishita et al. 2021), which implements the wind field with high resolution to calculate the tephra deposit load.The best TSPs are then selected based on the comparison of model results with the observed tephra deposit load, which is derived from the number of tephra particles classified by their diameter and settling velocity measured by a disdrometer network using an empirical conversion formula (Takishita et al. 2022).

Regional setting
The current study focuses on the Sakurajima volcano (Fig. 2a, b), one of the most active and closely monitored volcanoes in Japan (Iguchi et al. 2019).It has two active craters, namely Minami−Dake and Showa (Fig. 2c). Figure 3a shows the annual number of explosions at the Sakurajima volcano from October 1955 to December 2021 and the annual ejecta from 1979 to 2021.Since 1955, the typical eruptive style for the volcano has been ash-rich Vulcanian eruptions accompanied by impulsive sounds with shockwaves that eject volcanic bombs and ash, forming mushroom-shaped volcanic clouds.These eruptions are triggered by the destabilization of the brittle plug over the vent caused by a pressure increase (Iguchi et al. 2008).Although the activity of the Sakurajima volcano is commonly treated as typical Vulcanian explosions (i.e., instantaneous explosions), in reality, the tephra emission at a lower plume height that follows these explosions lasts for several hours to a day (Poulidis et al. 2019a) and is accompanied by long-lasting tephrafall.Some plumes reach a maximum height of 9.5 km avl (e.g. a Vulcanian eruption on 2:59 JST, Jun 4, 2020; Meteorological Research Institute 2020).Generally, the plume height and the frequency of explosions are inversely correlated up to 5 km avl.
Figure 3b shows the spatial distribution of cumulative tephra load at stations within 20 km of the Minami− Dake crater from 1991 to 2020 based on the data collected by the Kagoshima Prefectural Government.In total, over the last 30 years of activity, the load exceeded 30 kg m −2 at all sites on Sakurajima.At the sites closest to the Minami−Dake crater (Arimura, Futamata-ue, Kurokami, and Yunohira), the load exceeded 300 kg m −2 .Away from the volcano, tephra-fall maxima are located at sites southeast of the vent.This trend is consistent with the frequency distribution of wind direction and wind speed (Fig. 3c) based on radiosonde observations by the Kagoshima Local Meteorological Observatory, highlighting the impact of wind on tephra transport.

TSP time series TSP profile function
In this study, the TSP is derived from Eq. ( 1).The Suzuki TSP has been used in numerous settings in modeling studies, either coupled with an advection-diffusion model for sub-Plinian to Plinian eruptions (k = 4; HAZ-MAP; Pfeiffer et al. 2005), such as the Vulcanian eruptions at Sakurajima (k = 4; WRF-chem and FALL3D; Poulidis et al. 2017;Poulidis et al. 2021) or as a control test along a three-dimensional plume model based on conservation laws (k are between 1 and 6; Cao et al. 2021) for ultra-Plinian eruptions.In order to cover an exhaustive parameter space based on the previous uses of the Suzuki TSP, we here examine four values (1, 2, 4, and 8) for the k parameter.

TSP time series and TTDM calculations
The TSP is based on the observed initial plume height reported by the Japan Meteorological Agency (JMA), which is the agency in charge of monitoring activity from the volcano, and a geophysically estimated MER signal.Specifically, the JMA assigns a plume height based on the maximum visually observed plume height, but it does not monitor temporal changes in the plume height.Owing to this lack of observational data, we convert the temporal change of the MER into the temporal change of plume height based on the relationship between the MER and plume height (e.g., Morton et al. 1956;Mastin et al. 2009;Aubry et al. 2023).
The time series of the MER (10 3 kg min −1 ) is estimated at 1 min intervals based on geophysical observations (ground deformation data and the seismograph spectrum) using the linear-regression-based empirical method proposed by Iguchi (2016) where A is the sum of the seismograph spectrum between 2 and 3 Hz (m s −1 ), V is the pressure source volume change sum (m 3 ), c 1 and c 2 are the conversion factors, and c 3 is the correction term.The values of the conversion (2) MER estimated using this method has been successfully used in case studies (Poulidis et al. 2019a, b).We use the 1 min seismic and ground deformation data measured at the Arimura Volcanological Observation Tunnel located approximately 2.2 km southeast of the Minami−Dake summit crater.
A time series of TSPs is made by discretizing the MER signal from 15 min before the eruption onset up to 117 min after the onset in 3 min intervals.The upper time limit is set so that the frequency of the repeated explosions is sufficiently low.The plume height at time t is estimated from the relationship with MER (e.g., Morton et al. 1956;Mastin et al. 2009;Aubry et al. 2023): where h 0 is the plume height reported by the JMA (km avl) and MER 0 is a maximum MER at the time step of the eruption onset.For events during which the plume entered a cloud, the value measured by an X-band MP radar (Iguchi et al. 2019) is used for h 0 .Although X-band radars are reported to underestimate the plume height as compared to the visual observations by approximately 25% (Kobori et al. 2022), we do not make any calibration because the apparent underestimation is caused by fine ash which is not observable by radars.Note that the residence time in the plume is considered in the TTDM calculation meaning that the temporal variation of plume height is not simultaneous with that of MER.The explosive plumes in this study are unsteady; the duration is shorter than the time for the plume to rise.Although variations in MER influence the initial phases of plume rise, mixing and collapse (Chojnicki et al. 2015), as it is suggested that the classical entrainment closure (Morton et al. 1956) applies to unsteady plumes (Woodhouse et al. 2016), Eq. (3) was used for the plumes studied here as well.Since plume height is known to be reduced by the wind (e.g., Degruyter and Bonadonna 2012;Woodhouse et al. 2013), constraining the relationship between plume height and the MER based on the observation of the explosive column formed at the eruption onset reduces uncertainty by implicitly including the reduction of plume height caused by wind.Equation ( 3) is made under an assumption that the effect of wind-driven reduction in plume height is constant from the eruption onset to the end, and such an assumption is suggested to be valid in 2 h, the considered period of eruption.
By including the temporal change of plume height in the studied activities, we consider TSPs not only for the instantaneous explosion, but also for three other types of characteristic activity, namely (i) quickly repeated impulsive plumes with intervals of several minutes or tens of minutes, (ii) long-lasting and fading explosions, the duration of each varying from several minutes to tens of minutes and the intensity weakening with time (hereafter referred to as MER decay), and (iii) quasi-steady plumes with a lower height than that of the initial impulsive plume produced after an impulsive plume.The ability to reproduce different styles of activity allows the total TSP (i.e., the TSP based on all phases of the eruption) to provide substantial tephra that segregate at an altitude lower than the maximum height of the event, even if a fixed value of k is used.For accurate tephra deposit load calculations of smallscale Vulcanian eruptions, we used the TTDM Tephra4D (Takishita et al. 2021).Analytic reanalysis data (Mes-oScale Model; MSM provided by JMA) with a horizontal resolution of 5 km are downscaled into approx.0.5 km (0.005°) using the Weather Research and Forecasting (WRF) model (version 4; Skamarock et al. 2019).The detailed procedure is described in Iguchi et al. (2022).Atmospheric field data were interpolated horizontally into 0.25 km (0.0025°) and vertically from several isobaric surfaces with inconstant intervals of altitude to every 200 m from the sea level.Other input parameters used in this study are described as part of the Supplementary Material.

Disdrometer observations
As Vulcanian eruptions are small-scale and eject relatively small amounts of tephra, deposits can be swept away by wind and rain or can be mixed with deposits released during preceding or later events (Poulidis et al. 2019b).As such, an automated in-situ high temporal resolution observation is required for observing tephra deposit loads in real time.Optical disdrometers were selected for this.
Optical disdrometers are laser-based devices that count the precipitation particles and classify them into diameter and settling velocity classes.Disdrometers have increasingly been used to observe tephra sedimentation in these years (e.g.Kozono et al. 2019;Freret-Lorgeril et al. 2019;2022;Marchetti et al. 2022).A disdrometer network has been established to monitor the tephra sedimentation around the Sakurajima volcano (Iguchi et al. 2022).Takishita et al. (2022) investigated the characteristics of diameter-settling velocity distribution of tephra sedimentation measured by the network, in comparison with the collected tephra deposit beside them.
Specifically, the disdrometers installed are OTT's Parsivel 2 model, which is designed to detect particles larger than 0.25 mm in diameter.However, Takishita et al. (2022) noted that tephra-fall was also detected when the sedimentation contained few particles larger than 0.25 mm and seemed not to contain aggregates.This result indicates that disdrometers also detect particles smaller than 0.25 mm when the number density is sufficiently high.Importantly, this also suggests that disdrometers do not measure diameter accurately in such situations.In addition, results suggested that the disdrometers have a sedimentation rate threshold: most events where detected when the tephra sedimentation rate exceeded 10 g m −2 min −1 , while they were rarely detected when the tephra sedimentation rate was below 1 g m −2 min −1 .In addition to this analysis, we performed 0.125 mm sieving of the same set of samples collected by Takishita et al. (2022).The fraction of the particles smaller than 0.125 mm exceeded 50 wt% in 28 samples among 44 (the detail of the result is described in Supplements).Hence, about two-thirds of the tephra sedimentations detected by the disdrometer network in Sakurajima volcano consist mainly of particles smaller than 0.125 mm.
Based on these observations, Takishita et al. ( 2022) obtained an empirical formula to convert the number of detected particles in each diameter-settling velocity class into a tephra deposit load.In this formula, the inner product of the vector of the number of detected particles and the vector of the load per detected particle, for each diameter and settling velocity class respectively, are calculated.
One of the advantages of tephra-fall observations using disdrometers is that they can measure settling velocity.Only diameter is generally assigned to classify deposits in the sample collection method, the most generally adopted tephra sedimentation observation.Following this, the particle size distribution, not the settling velocity distribution itself, is commonly assigned as the input parameter of advection-diffusion models.However, the terminal velocity of particles with the same diameter varies because of the particle density and shape parameters.Here, we mitigated this issue by simply handling a settling velocity distribution without any special options (e.g.Bonadonna et al. 2002).Note that we assume that the settling velocity measured by disdrometers is the terminal velocity of each particle.In the current study, to estimate the tephra deposit load in each settling velocity class every minute, we updated the formula obtained by Takishita et al. (2022) by adding samples and disdrometer data and filtering the diameter-settling velocity class to be calculated based on mass fraction instead of number fraction.The settling velocity range of the updated formula is 0.6-7.2m s −1 .

Criteria for chosen eruptions
A total of 39 of the 668 eruptions that occurred between May 2018 and November 2019 were selected as the ground truth for the simulations (Table 1).These eruptions were selected because they had disdrometer observations.The tephra deposit load was estimated at 1 min intervals from the particle number in each diametersettling velocity class measured by the disdrometer at Sakurajima (Takishita et al. 2022).Specifically, the selection criteria were the detection of tephra at 3 or more of the 17 disdrometer sites (Fig. 4), in order to get a variation of distance or direction from the vent in the observed data, and the lack of rainfall during the eruption.As the eruptions occurred at the Minami−Dake summit crater, their coordinates were set to 31.5806°N 130.6580°E and 1000 m above sea level.The data collected within 3 h after the eruption onset were used in the analysis.For the cases where rain or a tephra deposit from another explosion was detected between 2 and 3 h after the onset, the end of the analysis period was brought forward 1 h (i.e., the analysis period was shortened by 1 h).

Table 1 Eruptions used for analysis
The ID number assigned to each eruption was created as follows.The first two digits correspond to the last two digits of the year when the eruption occurred and the last three digits correspond to the number in the annual JMA catalogs.The eruption onset is presented in Japan Standard Time (JST; UTC + 9).For the eruptions for which the time is marked with an asterisk, the end of the analysis period was set to two hours after the onset.A direction of "T" indicates that the plume moved steadily upward.The ejecta is the sum derived using the method mentioned in Sect."TSP time series and TTDM calculations".The wind speed is the mean value of the wind speed between the altitude of the vent and the top of the plume at the horizontally nearest grid to the vent."Working site" is the number of the sites that were working during the eruption and "Detected site" is the number of the sites that detected tephra-fall

Eruptive activity characterization
Here, we establish criteria to distinguish an impulsive plume from a quasi-steady plume based on the MER.Note that we distinguish "events" as the series of an eruption from 15 min before the eruption onset to 117 min after the onset, consisting of impulsive plumes and quasisteady plumes, with "explosions" that last for a few minutes.A histogram of the 3 min means of MER among 39 events is shown in Based on the above threshold, the total number of identified explosions increased to 82.Note that 20 of the 39 events included multiple explosions.For 22 events, 25 explosions continued to emit plumes until the next 3 min time step or later.For 12 events, quasi-steady plumes with MER values of more than 3 × 10 3 kg min −1 formed for more than 60 min in total.
The total TSPs are classified into four types based on the following criteria: • The total TSP contains a large mass fraction of the particles that segregate from a fading explosive plume or repeated explosive plume.There is a local maximum of the mass fraction in the initial TSP (i.e., the TSP of the initial plume) between the maximum altitude of the quasi-steady plume and 0.8 times the altitude at which the maximum mass fraction is.• The total TSP contains a large mass fraction of the particles segregating from a quasi-steady plume.The maximum concentration in the altitude range of the quasi-steady plume exceeded the maximum of that in the range of the explosive plume.
The events that meet none of the criteria, only the former criterion, only the latter criterion, and both criteria are classified as types A, B, C, and D, respectively.The numbers of events for these types are 18, 11, 5, and 5, respectively.More than half of the events are classified into types other than A (i.e., total TSPs with a significant fraction of particles segregating to an altitude below the maximum indicated by the initial TSP).

Tephra deposit characteristics
In order to simulate tephra particle transport and deposition using the TTDM Tephra4D (Takishita et al. 2021), the total grain size distribution is required (Bonadonna  and Houghton 2005; in addition to the temporally evolving TSPs).However, as the aim here is to compare the calculated load with the load estimated using the disdrometer network data, we adopt the settling velocity as the representative parameter for particles because the disdrometer registers a cluster of particles smaller than the lower limit as a single particle with a larger diameter (Takishita et al. 2022).This approach comes with the added benefit of implicitly accounting for aggregation in the simulations.For each eruption, an optimal k is assigned for each settling velocity class.
The settling velocity classes are based on the mass fraction distribution for settling velocities of tephra deposits.Figure 6 shows the accumulated results for all 39 eruptions.The mass fractions of particles with velocities of 0.9-1.8m s −1 exceed 5 wt%.There is a prominent mass fraction peak for particles with velocities of 1.0-1.2m s −1 .The mass fractions for particles with velocities of 1.8-3.6 m s −1 are between 1 and 4 wt%.There is a moderate mass fraction peak for particles with velocities of 2.0-2.4 m s −1 .As individual settling velocity distributions for each eruption can be biased due to the small number of sampling locations, the cumulative distribution to the tephra-fall calculation for all eruptions.We set six settling velocity classes based on the two peaks, as shown by the dashed vertical lines in Fig. 6.We also show the individual distributions in the Additional file 1: Fig. S1.

Elimination of outlier events
After the simulations of the 39 events were carried out, in order to focus the analysis on cases of realistically replicated deposition, criteria for eliminating events were established using a threat score (TS), which is defined as the ratio of the number of sites where tephra-fall is detected in both the calculation and observation to the number of sites where tephra-fall is detected in either the calculation or observation: where L obs and L cal,k are the observed and calculated loads, respectively, for each eruption and N is the number of sites.Hereafter, for k values of 1, 2, 4, and 8, the calculated load is denoted as L cal1 , L cal2 , L cal4 , and L cal8 respectively.Illustrative examples of the relationship between TS and the spatial distribution of tephra for six cases are shown in Fig. 7.A TS value of less than or equal to 0.5 represents cases where tephra sedimentation occurs in areas very different from observations, and a ( 4) TS value greater than 0.5 corresponds to good agreement between the model and observations.The eruptions in which the maximum TS for the six settling velocity classes exceeds 0.5 was used for k evaluations; 33 of the 39 eruptions (~ 85%; all except #18149, #18456, #18477, #19145, #19905, and #19365) meet this criterion.
Note that the evaluation performed here does not aim to be a general model evaluation (in which case eliminating events would lead to artificially improved sensitivity of k values).The aim of event exclusion here is to eliminate events whose simulations led to unphysical results for a number of possible reasons (e.g., incorrect plume height data or local bias in the meteorological data; Poulidis and Iguchi 2021), since including such results would add uncertainty to the main analysis of this study.Furthermore, only the presence or absence of tephra is used, not the calculated or observed load.

Selection of an appropriate k value
The root mean square ratio (RMSR) is adopted to evaluate the agreement between observations and simulation results: For each eruption, the k value that minimizes the RMSR was selected.For cases where the difference between the maximum and minimum RMSR was 0.1 or smaller, an optimal k was not selected.
As the number of events is too large to show all individual results in detail, the spatial distribution of simulated sedimentation and observations of the total tephra deposit load is shown for a specific case (#19147) as an illustrative example of the selection process (Fig. 8).This case was chosen because the eruption contained repeated explosive plume emission, MER decay after an explosion, and quasi-steady plume formation (i.e., all eruption styles discussed in Sect."TSP time series and TTDM calculations").Similar plots for all other cases are shown in the Supplementary Material (Additional file 1: Figs.S2-S34).As expected from the TSP profile characteristics, L cal1 indicates a larger tephra deposit in the proximal area (within 2 km from the vent; e.g., Fig. 8t), and L cal8 generally indicates more tephra distributed towards the distal area (farther than 3 km from the vent; e.g., Fig. 8b).For In the subplots, "NA" indicates that both the calculation and observation results were zero at all sites the specific case shown here, for particles with v t = 0.6-0.9 and 0.9-1.2m s −1 , L cal8 underestimated the observed load by two orders of magnitude at proximal site A, overestimated one by one order of magnitude at site B, and comparable to one at site C (Fig. 8a, b).The underestimations of L cal1 and L cal2 were less than that of L cal8 at the proximal site; however, L cal1 and L cal2 also underestimated the observed load at site C (Fig. 8m, n, s, t).For particles with a higher settling velocity, the calculated loads were zero at all sites, as shown in the map.Based on the RMSRs (shown above each subplot), an optimal value of k = 4 was selected for settling velocities 0.6-0.9 and 1.8-2.4m s −1 , and an optimal value of k = 8 was selected for settling velocities 0.9-1.2 and 1.2-1.8m s −1 .No selection was made for settling velocities larger than 2.4 m s −1 .

Overall dataset evaluation
Figure 9a shows the RMSR distributions between L cal1 , L cal2 , L cal4 , andL cal8 and L obs .In the two settling velocity classes 1.8-2.4 and 2.4-3.6 m s −1 , there is little difference among the RMSR distributions calculated with the four k values, indicating that the TSP has a negligible effect on the level of agreement between tephra deposit load calculation and observation.For the particles with the three smallest and largest settling velocity classes, the medians of the RMSR tend to decrease with increasing k. Figure 9b shows the frequency of the optimal k based on the 33 cases for each settling velocity class.For the three lowest settling velocity classes, the optimal k was determined to be 8 for ~ 50% of the events and 1 for 15-25% of the events.For the larger classes, an optimal k was Fig. 9 a RMSR distributions between calculations with k = 1, 2, 4, and 8 and observations.The dashed lines separate the settling velocity classes.The box-and-whisker diagrams are based on the median (orange line), the 25-75th percentile (boxes), and the range of 1.5 times the box (whiskers), with outliers circled.b Frequency of optimal k among all cases for each settling velocity class.The numbers in the legend indicate the optimal k and NA (not applicable) indicates that the optimal k was not selected because the difference in RMSRs between k = 1, 2, 4, and 8 did not exceed the threshold.c The relationship between optimal k and frequency distribution of MR distribution for the settling velocity class 0.6-2.4m s −1 .N vt is the number of total settling velocity classes among 39 eruptions not selected for ~ 50% of the events.For the remaining events, k values of 1 and 8 led to the optimal solutions in most cases (approximately 65-75% for each value).
For the four lowest settling velocity classes (v t < 2.4 m s −1 ), we introduce the mean ratio (MR) as an evaluation index to evaluate the bias of overestimation or underestimation of calculations as compared to observations at all sites: The effect of the optimal k value on the frequency distribution of all settling velocity classes classified by MR is shown in Fig. 9c.For the events with an optimal k value of 8, the majority of MRs were within ± 0.1 or exceeded 0.1, indicating good agreement and overestimation, respectively.For the events with an optimal k value of 1, ~ 50% of the MRs were below − 0.2, indicating major underestimation.As the most frequent optimal k values were 1 and 8, only results for these values are discussed below. (6

Characteristics of TSP time series
An example of the total and initial TSPs for eruption #19147 for k = 1 and 8 are shown in Fig. 10.The time series of the geophysically estimated MER is also shown.Four explosions occurred within 30 min of the eruption onset, with MER decay after the first two explosions.
For k = 8, particles segregate at the top of the plume in the initial profiles, but the total TSP also includes a fraction of the particles that segregate at the bottom of the plume, reflecting the changes in the MER.For k = 1, the mass is negatively correlated with the altitude or release for both the initial and total TSPs, with the overall impact of the quasi-steady particle emission being an increase in mass close to the surface.For both k values, the results suggest that not only the plume emission at the onset of the eruption but also the particle emission that follows is important in the total profile.It can be seen that the profile in each time step changes drastically, indicating that the temporal change of the plume height should not be neglected when obtaining the time series of the tephra deposit load.

Evaluation of sedimentation temporal evolution
Here, the time series of observed and calculated loads are compared.The time series of L obs , L cal1 , and L cal8 for eruption #19147 at the sites where the total L obs , L cal1 , and L cal8 all exceed 0 are shown in Fig. 11a-d.The locations of the sites are shown in Fig. 11f.The results differ among the disdrometer sites.At site HKU, the largest amount of tephra fell in L obs , L cal1 , and L cal8 in terms of the maximum 5 min load and duration.Within 60 min after the eruption onset, the periods when the 5 min load exceeded 10 g m −2 almost overlapped for L obs , L cal1 , and L cal8 .The maximum 5-min load was underestimated for L cal1 and L cal8 with v t < 2.4 m s −1 .The tephra-fall over 10 g m −2 later than 60 min after the eruption onset was not reproduced for L cal1 or L cal8 .
At site HKD, the trend differs between the period until one hour from the eruption onset and later.One hour after the eruption onset, the 5-min L obs exceeded 10 g m −2 for 15 min for velocities of 0.9-1.2m s −1 and 5 min for velocities of 0.6-0.9m s −1 .L cal1 exceeded 10 g m −2 for more than 20 min for velocities of 0.9-1.2 and 1.2-1.8m s −1 and the maximum exceeded 100 g m −2 , leading to overestimation.For L cal8 , the 5 min value exceeded 10 g m −2 within 10 min, which was slightly shorter than that for L obs .The observed tephra-fall that exceeded 10 g m −2 later than 100 min after the eruption onset was not reproduced.
At site AKA, the periods when the 5-min tephra deposit load exceeded 10 g m −2 (e.g., between 40 and 75 min from the onset) were reproduced for both L cal1 andL cal8 , with a temporal difference of around 10 min.For v t < 0.9 m s −1 , only L cal1 was overestimated.For particles with 0.9 < v t < 1.2 m s −1 , the maximum value of L obs was 10-30 g m −2 and those of L cal1 andL cal8 were 30-100 g m −2 , slightly exceeding that of L obs .For the v t bin of 1.2-1.8m s −1 , only the maximum value of L obs exceeded 10 g m −2 .
Finally, at site SBT,L cal8 was closer to L obs than L cal1 for both the 5 min load and duration.In summary, the reproducibility was generally high within 60 min of eruption initiation and low thereafter.The arrival periods of particles in L cal1 and L cal8 that exceeded 0 overlapped to a large degree and the 5 min load seems to have affected the total load.For L obs at sites HKD, HKU, and AKA, the length of the arrival period significantly exceeded 30 min.From the TSP in Fig. 10, the repeated explosions finished 30 min after the onset and the quasi-steady plume followed.Therefore, the particles reaching each site during the later period might be segregated from the quasi-steady plume.
The arrival periods for L obs , L cal1 , and L cal8 at each site are summarized in Table 2.At sites HKD and HKU, the particles from the quasi-steady plume continued to arrive after the arrival of the particles from the explosive plumes in both L cal1 and L cal8 .The particles from the explosive plumes finished reaching each site 68-70 min from the eruption onset.This time is included in the periods in which no particles were detected in L obs , that is, 60-75 min at HKU and 50-95 min at HKD and from the eruption onset, respectively (Fig. 13a, b).In contrast, the particles from quasi-steady plumes finished reaching each site 108-139 min from the eruption onset.The particles detected after the detection break at sites HKU and HKD are thus likely to have sedimented from the quasisteady plume.In addition, the length of the detection period of L cal8 from the quasi-steady plume is closer to that of L obs than to that of L cal1 .
In summary, the particles from the quasi-steady plume were likely detected in the observations and reproduced in the calculations.However, the estimation accuracy is poor for the period when the particles from the quasisteady plume reach a site (more than one hour from the eruption onset).

Relationship between settling velocity and arrival time
Here we define t arr as the time between the eruption onset and the arrival of a particle with each settling velocity class at first.t arr at each site for eruption #19147 is shown in Fig. 11e.For the observations, particles in all settling velocity classes arrived first at HKU, which is the closest to the vent and the highest in altitude, and arrived last at AKA, which is the farthest from the vent and lowest in altitude.On the other hand, for the calculations, particles also arrived first at HKU, but the order of the other three sites was reversed.For L cal1 and L cal8 , t arr at HKU for all settling velocity classes and t arr at HKD and SBT for settling velocity classes above 0.9 m s −1 are the same.For other combinations of sites and settling velocity classes, the differences in t arr are generally small.For all the combinations of sites and settling velocity classes except particles at HKU and SBT with settling velocities of 0.6-0.9m s −1 , the calculated t arr were sooner than the observed times.Regarding the relationship between settling velocity and t arr , for the calculations, particles with a lower settling velocity tend to arrive later than those with a higher settling velocity, but for the observations, particles in more than half of the settling velocity classes arrived simultaneously at all sites.

Fig. 12
The probability distribution of a particle with the indicated settling velocity class included in the particles that arrived at each site in the first minute.For the case where the period is longer than 15 min when L obs , L cal1 , and L cal8 are also zero, the events are separated and those during which the particles in four or more settling velocity classes arrive are extracted.In each region (separated by dashed lines), the lower plots are the loads with lower settling velocities.The y-axis is the settling velocity of the particles (m s −1 ) The difference in t arr for L obs , L cal1 , and L cal8 , shown in Fig. 11e, is examined here.Figure 12 shows the probability distribution of a particle with the indicated settling velocity class included in the particles that arrived at each site in the first minute.A comparison between simulation results and observations reveals significant deviations.For the simulation results, particles with settling velocities of 0.6-0.9m s −1 tend not to be included in the particles that arrive in the first minute (< 40%) and the highest probability among all classes is only slightly higher than 60%.For the settling velocity class between 0.9 and 1.2 m s −1 , there was a similar trend to the class between 0.6 and 0.9 m s −1 in L cal8 and to the classes between 1.2-1.8 and 1.8-2.4m s −1 in L cal1 .For the observations, the particles in all the settling velocity classes simultaneously arrive for most events (≈80% at lowest).For the calculations, particles with a higher settling velocity generally arrive sooner, whereas, for the observations, particles with a wide range of settling velocities generally arrive simultaneously.

Characteristics of total TSP
In the analysis presented in the previous sections, for most of the events studied here, rather than a skewed profile (i.e., k = 1), a top-loaded (k = 8) TSP with multiple maxima due to the temporal evolution of the MER is the most likely reason for sedimentation from the lower parts of the plume.This creates distinct total TSPs for the four eruptive types characterized by the MER temporal evolution (types A-D; Sect."Criteria for chosen eruptions").
A comparison of the typical total TSP with k = 8 with the TSP at the eruption onset (dashed line) is shown in Fig. 13.Type A eruptions show one main concentration maximum, meaning that the total TSP is almost identical to the initial profile (46% of events; first row in Fig. 13).

Fig. 13
Comparison between mass distribution profiles for total TSPs (solid lines) and initial TSPs (dashed lines) for k = 8 for a-e total TSP almost identical to initial TSP (type A), f-j total TSP with large mass fraction of particles segregating from a quasi-steady plume (type B), k-o total TSP with large mass fraction of particles segregating from a fading explosive plume or repeated explosive plumes, and p-t total TSP with characteristics of types B and C. The y-axis is the segregation height with intervals of 0.5 km Type B eruptions show two maxima at altitudes similar to those in the original TSP (28% of events; second row in Fig. 13).Type C eruptions show an additional distinctive maximum near the surface with lower or equal magnitude to the initial distribution (13% of events; third row in Fig. 13).Finally, for type D eruptions, particles segregate at a wide range of altitudes between the maximum plume height and the surface (13% of events; fourth row in Fig. 13).Type A is the typical profile of an explosive phase of the eruption.The profiles for types B, C, and D include particle segregation at an altitude lower than the maximum plume height, indicating the existence of quasi-steady particle emission with a lower plume height and long duration, or multiple explosions with different MER and plume heights.
Note that despite the use of a single value for k for all of the profiles shown here (k = 8), the effective maximum for the total TSP is reduced for type B and D eruptions (e.g., Fig. 13h, i, r, s, t).This shows that even for a fixed value of k, accounting for the temporal evolution of the MER can replicate the impact of a lower k.

Parameter for constraining TSP
The main objective of this study was to understand tephra sedimentation that occurs at a height below the maximum plume height.Two factors were examined, namely the sensitivity of the simulation results to the k value and the impact of the temporal evolution of the TSP.Accounting for changes in the MER for a given k value was found to lead to realistic total TSPs with multiple maxima.
Allowing for a temporally evolving TSP enabled us to characterize the total TSP based on the underlying activity, namely single explosion, quickly repeated explosions, smooth reduction of the MER during the eruption, or the formation of a quasi-steady plume.Of the four types, a typical TSP with tephra segregation at the top of the plume (corresponding to a single explosion) was common but represented less than half of the events under study.All other TSPs featured a primary or secondary maximum beneath the plume height.This is similar to the situation observed for Plinian eruptions, which have larger magnitudes and less temporal variation in plume height, and perhaps top-concentrated TSPs can reproduce TSPs estimated from observations (Fero et al. 2008(Fero et al. , 2009;;Mannen 2014;Cao et al. 2021) by considering the temporal changes in plume height with a high temporal resolution.
A k value of 8 was optimal for 33 of the 39 events studied, suggesting that for each eruption phase, almost all particles segregate at the top.This profile is consistent with ones given in previous studies (Schwaiger et al. 2012;Hurst and Davis 2017;Trancoso et al. 2022).Suzuki (1983) suggested that larger particles segregate from the lower part of the plume than smaller particles, but such a trend is not seen in our results.One possible reason for this may be that the range of k is shorter than the interval of the choices of k in this study.Such suggestions have also been made in previous studies based on 1D plume modeling.Ernst et al. (1996) and Girault et al. (2014) gave a fundamental assumption that the mass of the segregated particles per height is proportional to the terminal velocity of the particle (Eq.( 1) and Eq. ( 9) respectively), but the resulting TSP was such that most of the particles segregate from the top of the plume.
The results here provide evidence that commonly occurring activity at the Sakurajima volcano differs from the textbook definition of Vulcanian eruptions as discrete instantaneous events (Clarke et al. 2015).Instead, some eruptive events at Sakurajima can be seen as the result of repeated Vulcanian explosions of varying intensities, often accompanied by the passive release of material before a new plug starts being formed.Results here suggest that rather than analyzing isolated transient eruptions, certain eruptive events at Sakurajima can be approached as quasi-steady eruptions, mirroring the ideas of Wilson et al. (1978) which postulated similar findings for repeated strombolian eruptions during the 1973 eruption of Eldfell, Iceland.The repetition of transient events may explain why the top-concentrated plume accurately describes the activity of this volcano.
This has important ramifications for the operational management of the volcano.When the plume reaches its maximum height, the MER is also maximum in the series of eruptions and the total TSP tends to be such that most of the particles segregate at the maximum plume height.For Vulcanian eruptions, the focus is on the initial explosion, which produces a plume with the maximum height.In JMA operations, the maximum plume height is thus recorded as the observation that characterizes the eruption.However, among the eruptions of the Sakurajima volcano identified as Vulcanian eruptions, sometimes there are quickly repeated explosions of the same or gradually reduced intensity that may be followed by longlasting quasi-steady particle emission.In these cases, the plume height becomes lower than the maximum height reached by the initial explosion.Then, the mass fraction of particles that segregate at the lower height increases.Considering the variation of the plume height, even for k = 8, which makes a TSP in which most of the particles segregate at the top, the total TSP has a local maximum at an altitude lower than the maximum for more than half of the events.In the Vulcanian eruptions and the series of following eruptive activities, the plume height drastically changes, indicating the importance of considering the variation of plume height, which is also suggested by the total TSP.As the plume height is estimated based on the MER of the total activity, the MER can be used to control the TSP.Some inconsistencies do persist even after the application of a temporally evolving MER.The underestimation of the tephra deposit load, during the arrival of the particles from the quasi-steady plume, is likely caused by the MER or plume height of the quasi-steady plume being inaccurate.The factors and terms in Eq. ( 2) are estimated based on the monthly tephra deposit loads.Equation ( 3) is only constrained by the parameters at the onset.The factors and terms in Eq. ( 2) should thus be separately estimated for explosive particle emission and quasisteady particle emission.

Re-examination of the relationship between MER and plume height
For a non-negligible fraction of the events studied, the optimal value for k was 1, indicating particle segregation at the lower part of the plume.In these cases, calculations with k = 8 tend to produce underestimation (Fig. 9c).Possible causes of this underestimation are the small number of observation sites in the regions reached by the tephra particles, the discrepancy between the assumed and actual total settling velocity distributions, and the high observed plume height relative to the estimated MER.To verify this, we re-examine the relationship between the MER and the plume height.
The plume height at the eruption onset was constrained by the observations and that after the eruption onset was estimated using the quarter-power-law between the MER and plume height (e.g., Morton et al. 1956).As the entrainment of the ambient atmosphere is known to contribute to plume ascent, Degruyter and Bonadonna (2012) proposed the following relationship between plume height h p (km avl) and the buoyancy flow rate at the vent Ḟb : where α is the entrainment coefficient under a calm wind (determined to be 0.1 through observations and experiments; Morton et al. 1956;Carazzo et al. 2008;Devenish et al. 2010), f BV is the buoyancy frequency, and z 1 is the maximum non-dimensional height (determined to be 2.8 through numerical integration of the non-dimensional governing equations; Morton et al. 1956).( 7) Degruyter and Bonadonna (2012) also proposed a relationship between Ḟb and h p ′ considering the deteriora- tion of the plume height by the wind: where β and w h are the wind entrainment coefficient and the horizontal wind speed, respectively.We set β to 0.5, an empirically obtained value (Briggs 1972;Devenish et al. 2010).The mean value between the vent and the top of the plume above the vent is applied to w h .
Ḟb is obtained from the following formula with MER Ṁ (kg s −1 ) (Morton et al. 1956;Wilson et al. 1980;Woods 1988;Glaze et al. 1997;Degruyter and Bonadonna 2012) where g is the gravitational acceleration.ρ a0 , c a0 , andθ a0 are the density, heat capacity, and temperature of the ambient air, respectively, and c 0 andθ 0 are the specific heat capacity and temperature of the source, respectively.Then, the relationship between plume height h p , h p ' , and MER Ṁ considering the entrainment can be expressed by Eq. ( 11) with no wind and Eq. ( 12) with wind: Here, ρ a0 is 1.23 kg m −3 , c a0 and c a are 998 and 1250 J kg −1 K, respectively, and θ a0 and θ 0 are 288 and 1200 K, respectively, after Woods (1988) and Degruyter and Bonadonna (2012).The mean value between the vent and the top of the plume above the vent is applied to f BV .
Figure 14 shows the relationship between h p (Fig. 14a), h p ' (Fig. 14b), and h 0 , colored with the optimal k for the settling velocity class 0.9-1.2m s −1 .MER 0 is substituted into Ṁ .When the wind effect is ignored, all the observed plume heights were lower than the estimated values (Fig. 14a; RMSE = 1.7 km among eruptions with an optimal k value of 8).When the wind effect is considered, the heights match (Fig. 14b; RMSE = 0.7 km among the same eruptions).From the limited cases studied here, we suggest that when the wind is considered, for eruptions with an optimal k value of 1, a higher proportion of the plume height is underestimated when compared with that for eruptions with an optimal k value of 8.The segregation (8) height is defined as the altitude where particles end their ascent in the plume and start to descend, leaving the plume.For a plume that grows upward, tephra particles leave the plume from the top horizontally.For a plume that is distorted by the wind, particles are released from the bottom of the plume.Therefore, for plumes whose growth is highly affected by the wind, the gap between their maximum altitude and the maximum segregation altitude will correspond to their width (e.g., Fig. 15).A larger width leads to a larger gap between h p ′ and h 0 .For the case where h 0 is higher than h p ′ , the TSP with k = 8 overestimates the segregation height of many particles, leading to an overall underestimation of the load (Fig. 8c).
In such cases, compensating for the gap, the profile with k = 1, lower segregation, and providing relatively heavier tephra load in the proximal area where the observation network is located, are suggested.The focus of the present work has been on the TSPmost other model parameters have been assumed as set.Of course, the modeling of tephra transport and deposition is complex and critically relies on a number of other input parameters.As Scollo et al. (2008) investigated, the total mass of tephra and plume height have the most significant impact to tephra deposit load calculation.Macedonio et al. (2016) also reported that the total erupted mass has a first-order effect on the extension of the hazard zone.Meteorological input data and aggregation scheme also affect the calculation (Poulidis and Iguchi 2021).A lot of work has been made to make the gap between input parameters and natural conditions both in theoretically and empirically in MER (e.g., Iguchi 2016, specialized to Sakurajima Vulcanian eruptions), plume height (e.g., Morton et al. 1956;Mastin et al. 2009;Aubry et al. 2023, derivation from MER of steady plume; Degruyter and Bonadonna 2012;Woodhouse et al. 2013, vending plume model), total grain size distribution (e.g., Costa et al. 2016), and the meteorological field (e.g., Poulidis et al. 2017;Takishita et al. 2021, consideration of orographic effect).To maximize the accuracy of the simulations presented here, an effort was made to use literature-based parameter values and methodologies appropriate for Sakurajima volcano.Naturally, inconsistences with observations remain and the contribution of each parameter used in Tephra4D to the inaccuracy of calculations needs to be evaluated in an appropriatelydesigned comprehensive parametric study.1) and calculated plume height (a) h p obtained from Eq. ( 11) and (b) h' p obtained from Eq. ( 12) Fig. 15 Ash fingers formed from a plume during the eruption at 12:21 JST on June 27, 2022.The photograph was taken 5 min after the eruption started

Inconsistencies in modeled tephra arrival times
Particles with a broad range of settling velocities arrive at observation sites simultaneously, as shown in Fig. 12.To explain this, we consider the case where particles with a high settling velocity leave the plume from the higher part and particles with a low settling velocity leave from the lower part; the two types of particle reach the observation sites at the same time.This trend appears in the relationship between the particle settling velocity and the initial segregation height of the first particle reaching each site in the simulations (Fig. 16), but the simultaneous arrival is not reproduced.
In the case where the particles mostly follow similar trajectories from the segregation from the plume to the ground, simultaneous arrival makes sense.This can be explained using the concept of ash fingers (e.g., Carazzo and Jellinek 2012;Scollo et al. 2017;Freret-Lorgeril et al. 2020;Fries et al. 2021;Lemus et al. 2021;Fig. 15).The overturning and stirring of the ash-air mixture enhance the production of particle boundary layers (PBLs) and the formation of ash fingers (Freret-Lorgeril et al. 2020).In the case where the settling velocity of the ash-air mixture is greater than that of individual particles, the mixture behaves as a continuum, and ash fingers form (Hoyal et al. 1999).The volume of fingers increases as they fall and the settling velocity quickly decreases.The fingers then evolve into sediment thermals (Freret-Lorgeril et al. 2020).When the settling velocity of a sediment thermal becomes lower than the velocity of each component, the component leaves the thermal and starts to fall at its individual terminal velocity.In the case where the period from the departure of a particle from the thermal to the deposit is short, particles with various settling velocities are detected on the ground almost simultaneously.In the case where the plume height is low and fingers form at a sufficiently low altitude, particles reach the ground as a continuum.As particles have various settling velocities due to turbulence in the continuum, they are simultaneously detected with various settling velocities.The plume height for the eruption shown in Fig. 15 is 1 km, which is sufficiently low for the fingers to keep their constitution close to the ground, suggesting that the particles reached the ground as a continuum.
Considering that ash fingers can aid with the interpretation of the results regarding settling velocity in this study, the velocity of the fingers is estimated from the disdrometer observations.It has been reported that even for tephra-fall in which particles smaller than 0.25 mm are dominant, the detected settling velocity is up to 3.2 m s −1 (Takishita et al. 2022, Fig. 7b).Focusing on the settling velocity classes under 3.2 m s −1 , we compared the total number of sites where the particles are detected in each settling velocity class between L cal8 and L obs among 39 eruptions (Fig. 17).The sites were classified based on the duration of tephra-fall.The total number of sites is almost the same among the four settling velocity classes between 0.6 and 2.4 m s −1 , suggesting that particles in these settling velocity classes form the finger.From the consistency of the total number of sites and the duration distribution, the mean settling velocity of the finger is estimated to be 0.9-1.2m s −1 .The estimated mean velocity overlaps the value estimated from visible image analysis of the Eyjafjallajökull volcano eruption in 2010 (Manzella et al. 2015), namely 1 ± 0.5 m s −1 around the median.Under the assumption that the observed range  Next, the validity of the source parameters of the Vulcanian eruptions of Sakurajima volcano is examined to explain the temporal characteristics of the observed tephra-fall within an ash finger.Regarding the MER threshold, Scollo et al. (2017) conducted laboratory experiments using water, sugared water, glass beads, and volcanic ash and showed that the threshold of the MER for forming the finger is 10 5 kg s −1 .Since the Sakurajima eruption at 12:21 JST on June 27, 2022, for which the plume height was 1500 m and the MER was 1.56 × 10 4 kg s −1 , formed ash fingers (Fig. 15); the MER threshold for finger deposition can be reduced by approximately one order.Most of the events (30 of 39) analyzed in this study exceeded the reduced threshold, suggesting that these events meet the conditions for which fingers are formed.Regarding the size of the constituents of the finger, Lemus et al. (2021) suggested a formulation for the threshold based on the finger velocity and the terminal velocity of the individual particle.We estimated the threshold in the same way using natural conditions (particle volume fraction X p = 1-4 × 10  Suzuki (1983) and shape parameter empirically estimated by Suh et al. (2019).As shown in Fig. 18, the threshold of particle diameter is estimated to be 0.19-0.41mm.When compared with our observations, as mentioned in Sect."Disdrometer observations"., the majority of the particles detected by our disdrometer network are smaller than 0.125 mm.Disdrometers should thus be able to detect ash finger constituents.While the maximum particle size may be affected by wind (Freret-Lorgeril et al. 2020), the size of the constituents and the MER in this study are suggested to satisfy the sufficient condition in a windless environment.

Conclusions
Considering the temporal variations in the volcanic plume height, the tephra deposit load was calculated using a large number of Vulcanian eruptions and continuous emissions at the Sakurajima volcano.The time series of the TSP in which most of the particles segregate at the top of the plume makes the calculations correspond to the spatial distributions of the tephra deposit load measured by the disdrometer network.The results show that the TSPs given in most modeling studies are indeed appropriate for individual events.Consideration of the temporal variations in the volcanic plume height allows the total TSP, even one that consists of plumes in which most of the particles segregate at the top of the plume, to lead to not only top-concentrated TSP but also feature a large fraction of particles segregating at an altitude lower than the maximum height of the event, suggested from the previous studies based on observations.As plume height depends on MER, the essential constraint of TSP is suggested to be MER.
The comparison of the calculations and observations of tephra deposit loads, made in this study is innovative in terms of applying a temporal constraint to the discussion of tephra transport.As the time series of the tephra deposit load obtained by the disdrometer network contains more information than the cumulative load, accuracy can be separately evaluated for each phase of the tephra-fall.Observations indicate that particles with a wide range of settling velocities reach the ground simultaneously; however, this was not reproduced by the calculations.A possible reason for this is the model's inability to resolve volcanic ash fingers caused by gravitational instabilities in the plume.A reanalysis of the results under this assumption revealed that the observed characteristics can broadly be explained using literature values associated with ash fingers, suggesting the importance of this phenomenon in tephra transport simulations with accurate calculations of the arrival time.

Fig. 1
Fig. 1 Effect of k parameter value in Suzuki TSP for plume height of 2.5 km avl factors and correction term (c 1 = 3.8 × 10 −5 kg min −1 (m s −1 ) −1 , c 2 = 2.6 kg min −1 m −3 , c 3 = − 1.03 × 10 −5 kg min −1 ) were based on a minimizing the Root Mean Square-Error (RMSE) against a training dataset based on the monthly ashfall data collected by the Kagoshima Prefectural Government at 62 stations between 2009 and 2013.

Fig. 2 a
Fig. 2 a Location of Sakurajima volcano (red triangle) and large cities in Kagoshima Prefecture.b General view of Sakurajima volcano (origin at Minami−Dake crater A) and c magnified view of a rectangle in b showing locations of Minami−Dake craters A and B and Showa crater around the summit of Minami−Dake

Fig. 3 a
Fig. 3 a Annual number of explosions (N ex ) from October 1955 to December 2021 and annual ejecta (Ej an ) from 1979 to 2021.M and S above the plot indicate the active periods of the Minami−Dake and Showa craters, respectively.b Spatial distribution of cumulative tephra load (L cum ) at stations within 20 km of Minami−Dake crater (red triangle) from 1991 to 2020.c Frequency distribution of wind direction and wind velocity w h (m s −1 ) at an isobaric surface of 850 hPa (about 1.5 km above sea level) from 1991 to 2020 Fig. 5.It consists of the major maximum for low MER values (quasi-steady plumes) and the secondary maximum for high MER values (explosive phase).To separate the two types of activity, based on the minimal value in Fig. 5, we used a threshold of 2 × 10 5 kg min −1 .Values over the threshold indicate an explosive eruption and the others indicate a quasi-steady plume.

Fig. 4
Fig. 4 Disdrometer observation network.AVOT (square marker) is the station where ground deformation and seismographs are measured

Fig. 5
Fig. 5 Histogram of 3 min means of MER for all eruptions

Fig. 7
Fig. 7 Examples of tephra sedimentation pattern for TS values of a 0, b 0.11, c 0.4, d 0.5, e 0.6, and f 0.8.In all cases, L cal8 is compared to L obs for the settling velocity class 0.9-1.2m s −1

Fig. 8
Fig. 8 Spatial distribution of simulated observed total tephra deposit load for eruption #19147.Different rows show simulation results for different k values and different columns show different sedimentation velocity bins.The intervals of the x and y axes are 2 km.In the subplots, "NA" indicates that both the calculation and observation results were zero at all sites

Fig. 10
Fig. 10 Time series of TSPs (left) and total (solid lines) and initial TSPs (dashed lines) (right) for a k = 8 and b k = 1 for event #19147.EXP and QS indicate explosions and quasi-steady plumes, respectively, and shading indicates the mass fraction.c Time series of geophysically estimated MER.Note that the time indicates the emission time of the particles segregating at each altitude.The residence time in the plume is not shown

Fig. 11
Fig. 11 Time series of L obs , L cal1 , and L cal8 for eruption #19147 at sites a HKU, b HKD, c AKA, and d SBT.In each region (separated by dashed lines), the lower plots are the loads with lower settling velocities.The y-axis is the settling velocity of the particles (m s −1 ).e Relationship between particle settling velocity and arrival time (the time when the first particle reaches each site) from the eruption onset.Colors indicate the sites corresponding to the marker colors in f. f Location of each site.The intervals for the x and y axes are 2 km

Fig. 16
Fig.16Relationship between particle settling velocity and initial segregation height of the first particle reaching each site in simulations.The box-and-whisker plots are drawn in the same way as done in Fig.8a

Fig. 17
Fig. 17 Settling velocity dependence of total number distribution of sites classified in tephra deposit duration.The left (right) side of the bar indicates the result for L obs (L cal8 ) in each v t class.Events were extracted in the same way as in Fig. 12 −6 ; PBL thickness δ PBL = 70-100 m, Manzella et al. 2015; atmospheric density ρ a0 = 1.23 kg m −3 , Woods 1988 and Degruyter and Bonadonna 2012; particle density ρ p = 2640 kg m −3 , Haruyama et al. 1977; dynamic viscosity μ = 3 × 10 −5 Pa s −1 , Carazzo and Jellinek 2012).The terminal settling velocity of the individual particle is calculated based on the formula suggested by

Fig. 18
Abbreviations avlAbove vent level MER Mass eruption rate

Table 2
Arrival period of eruption #19147 at sites shown in Fig.11The values indicate the time in minutes after the eruption onset