A physics-based model to estimate source conditions for a tephra-dispersal model is developed. The source condition is generally expressed by a distribution of released particles along an eruption plume (referred to as “source magnitude distribution” SMD). The present model (NIKS-1D) calculates the SMD and the column height for given vent conditions (e.g. mass eruption rate and magma properties) on the basis of an eruption column model below the neutral buoyancy level (NBL), a downwind gravity current model around the NBL, and a particle sedimentation model. It quantitatively reproduces the following features of the SMD for typical explosive eruptions: (1) a significant amount of coarse particles are released from the rising eruption column, whereas most of the fine particles are carried to the NBL, (2) in a downwind gravity current, the coarse particles tend to decrease more rapidly with distance from the vent than the fine particles, (3) the SMD from the downwind gravity current decreases with distance more slowly in a strong ambient wind than that in a weak ambient wind, and (4) the SMD from the downwind gravity current for eruptions with large mass eruption rates decreases with distance more slowly than that for eruptions with small eruption rates. NIKS-1D includes a new parameter, \(\mu\), which represents the ratio of the volumetric flux at the source of the downwind gravity current to that of the eruption column model at the NBL. This parameter is determined by the physics of the entrainment process around the connection between the eruption column and the downwind gravity current, and depends on the intensity of eruptions. We propose an empirical formula to calculate the value of \(\mu\) as a function of the mass eruption rate on the basis of the observation data from two well-studied eruption events. In a real-time tephra-dispersal forecasting system, NIKS-1D estimates the mass eruption rate from the observed plume height, and calculates the SMD from the estimated mass eruption rate as a source conditions for a tephra-dispersal forecasting immediately after an eruption.

Introduction

When an explosive volcanic eruption occurs, a large amount of tephra is injected into the atmosphere. The tephra fallout covers a wide area over short periods of time, and is the most widespread volcanic hazard. Besides the impact of tephra fallout on the ground due to its mass loading, air-borne ash constitutes a serious threat to aviation safety (e.g. Casadevall 1994; Guffanti et al. 2010; Clarkson et al. 2016; Liang and Xu 2021). For the hazard mitigation, issuing information about tephra-dispersal forecast immediately after an eruption is of utmost importance.

In tephra-dispersal models for a real-time forecasting, the distribution of released particles along a plume (“source magnitude distribution (SMD)”) is used as a source term. Therefore, the uncertainties of the tephra-dispersal forecast primarily inherit the uncertainties of the SMD (e.g. Folch 2012). To estimate the SMD accurately, understanding the plume dynamics is essential. In the real-time tephra-dispersal forecasting system of the Japan Meteorological Agency (JMA), the SMD is provided on the basis of an empirical model (Suzuki 1983; Hasegawa et al. 2015). The empirical model represents the SMD from a vertical plume, and does not consider the eruption plume dynamics or the sedimentation theory. On the other hand, recent studies have shown that the SMD critically depends on the dynamics of the eruption plume and the sedimentation processes from a turbulently mixed plume (e.g. Bursik et al. 1992; Bursik 2001; Koyaguchi et al. 2009; Girault et al. 2014, 2016). The amount of particles released from an eruption plume is governed by its shape and rising velocity, which is in turn controlled by the ambient wind and the vent conditions such as the mass eruption rate. For an improvement of the tephra-dispersal forecast for various eruption conditions, it is essential to evaluate the physical processes in the plume dynamics and the sedimentation processes.

In this paper, we develop a physics-based SMD model (referred to as “NIKS-1D” hereafter) for the tephra-dispersal forecasting that considers the plume dynamics (e.g. Woods 1988; Bursik et al. 1992) and the sedimentation theory (e.g. Bursik et al. 1992; Bursik 2001; Koyaguchi et al. 2009). To estimate the SMD immediately after an eruption, NIKS-1D also includes an inversion algorithm to estimate the mass eruption rate from the observed plume height. We expect that NIKS-1D improves the accuracy of the SMD used in JMA’s current real-time tephra-dispersal forecasting system.

Methodology

Framework of NIKS-1D

We are concerned with the amount of tephra-particles released from eruption plumes per unit length and unit time (i.e. the SMD). The SMD is controlled by the shapes of eruption plumes as well as the column dynamics (e.g. rising velocity). Generally, an eruption plume is composed of two parts: an eruption column below the neutral buoyancy level (NBL) and a downwind advection current around the NBL. The height of the eruption column is governed by the vent condition (e.g. mass eruption rate, Woods 1988), and its shape is influenced by the wind field (e.g. Bursik 2001). The downwind advection current around the NBL is described as the crosswind-direction spreading due to the gravity current (e.g. Bursik et al. 1992). The SMD can be expressed as particles released from these two parts of the plume.

The NIKS-1D model is defined as a combination of the three models: an eruption column model below the NBL, a downwind gravity current model around the NBL, and a particle sedimentation model (Fig. 1). The interaction between the plume and an ambient wind expresses the plume bending (Bursik 2001). The downwind gravity current model includes dynamics that express crosswind spreading of a plume as a gravity current coupled with downwind advection (Bursik et al. 1992). In addition, the particle sedimentation model describes the process of particles segregation from the side of the eruption column and the base of the downwind gravity current (e.g. Bursik et al. 1992; Bursik 2001; Koyaguchi et al. 2009). The combination of these three models makes it possible to express the plume dynamics and the SMD for various eruption conditions. See Table 1 for the notations and the preset values used in the later sections.

Eruption column model

Eruption column models calculate steady distributions of physical quantities such as density, velocity, and temperature with height by integrating conservation laws of mass, momentum, and specific enthalpy fluxes (e.g. Woods 1988). We use a steady 1-D model for a bending eruption column in a wind field developed by Bursik (2001), where the effects of the particle fallout are considered. This model takes into account an interaction between the eruption column and the ambient wind, and calculates the bending of the plume. In this model, the conservation laws along the plume are formulated as follows,

where s is the local coordinate for the vent location along the centerline of the plume, \(\varvec{v}\) is the average velocity vector of the plume which has a direction of the local coordinate s, \(\rho\) is the bulk plume density, r is the plume radius, \(\rho _{\mathrm {a}}\) is the ambient air density, g is the gravitational acceleration (\(\varvec{g}=(0,0,g)\)), \(\theta\) is the local angle of the plume (i.e. bending of the plume), \(\varvec{w}\) is the horizontal wind velocity vector, T is the temperature of the bulk mixture, \(T_{\mathrm {a}}\) is the ambient air temperature, \(C_{\mathrm {p}}\) is the specific heat of the bulk mixture, \(C_{\mathrm {pa}}\) is the specific heat of the air, and \(C_{\mathrm {ps}}\) is the specific heat of the pyroclast. \(-\frac{dQ_\phi }{ds}\) is the released mass along the plume per unit time and unit length (i.e. [kg/(s m)]), and we call it SMD (i.e. SMD\(\equiv -\frac{dQ_\phi }{ds}\)).

The set of the conservation laws of Eqs. (1)-(3) is complemented by the equation of state (Woods 1988),

where n is the gas mass fraction, \(R_{\mathrm {g}}\) is the gas constant of the bulk mixture, p is the pressure, \(\rho _{\mathrm {p}}\) is the density of the pyroclasts. The thermodynamical properties of the mixture are expressed by

where \(n_{\mathrm {a}}\) is the air mass fraction, \(n_{\mathrm {v}}\) is the volcanic gas mass fraction (\(n = n_{\mathrm {a}} + n_{\mathrm {v}}\)), \(R_{\mathrm {a}}\) is the gas constant for the air, \(R_{\mathrm {v}}\) is the gas constant for the volcanic gas, and \(C_{\mathrm {pv}}\) is the specific heat of the volcanic gas.

In the steady 1D eruption column model, the rate of the turbulent entrainment of the ambient air into the plume \(U_{\mathrm {e}}\) is parameterized as follows (Morton et al. 1956; Hoult et al. 1969; Hewett et al. 1971; Woods 1988; Bursik 2001; Devenish et al. 2010),

$$\begin{aligned} U_{\mathrm {e}} = \alpha \left| v - w \cos \theta \right| + \beta \left| w \sin \theta \right| \end{aligned}$$

(7)

where \(v=|\varvec{v}|\), \(w=|\varvec{w}|\), \(\alpha\) and \(\beta\) are the radial entrainment parameter and the wind-entrainment parameter, respectively. This formulation of \(U_{\mathrm {e}}\) is based on the assumption that the entrainment process for the bending eruption column depends on the two mechanics: the velocity differences parallel to the plume axis (the first term of Eq. (7)) and the velocity differences normal to the plume axis (the second term of Eq. (7)).

NIKS-1D solves the set of the conservation laws (1)-(3) using the fourth-order Runge-Kutta algorithm with a step size of ds, and provides the shape of the eruption plume (e.g. the heights of the plume top and the NBL) as well as the SMD from the input parameters related to the magma properties, the atmospheric condition, and the eruption condition at the vent. For given magma properties (see Table 1), the eruption condition at the vent is specified by giving the exit velocity \(v_0\) and the mass eruption rate \(\dot{M}\); the vent radius is calculated from \(\dot{M}=\pi r^2 \rho _0 v_0\). In this study, we assume that the exit velocity \(v_0\) equals the sonic velocity (i.e. Mach 1; \(v_0=\sqrt{nR_{\mathrm {g}}T}\)).

In a real-time tephra-dispersal forecasting, the mass eruption rate \(\dot{M}\) is required to execute NIKS-1D immediately after an eruption. Therefore, NIKS-1D includes an inversion to estimate the mass eruption rate \(\dot{M}\) from the observed plume height.

Downwind gravity current model

After the bending plume reaches the NBL, the eruption plume is laterally advected by the ambient wind while spreading orthogonally to the wind direction. This behavior of the laterally advected part of the bending plume is understood as a gravity current coupled with the ambient wind (referred to as the “downwind gravity current”). Under the assumption that the volumetric flux remains constant with distance, the downwind gravity current is calculated by a simple model as follows (Simpson 1987; Bursik et al. 1992),

$$\begin{aligned} \frac{dL}{ds} = \frac{\lambda N V_1}{w_{\text {nbl}}^2 L} \end{aligned}$$

(8)

where L and \(V_1\) are the width and the volumetric flux of the downwind gravity current, respectively, N and \(w_{\text {nbl}}\) are the Brunt-Väisälä frequency and the ambient wind velocity at the NBL, respectively, \(\lambda\) is the spreading constant of order unity that depends on a flow geometry (Bursik et al. 1992). The Brunt-Väisälä frequency N has the order of \(10^{-2}\) [1/s] for a typical atmosphere. Equation (8) expresses how the spreading of the downwind gravity current depends on the volumetric flux \(V_1\) and ambient wind speed \(w_{\text {nbl}}\). It has analytical solutions as follows,

where \(s_0\) denotes the position of the beginning of the downwind gravity current. Here, the thickness of the downwind gravity current h(s) is obtained from Eq. (9) and the relationship \(V_1=Lhw_{\text {nbl}}\). A source condition for the downwind gravity current model (the mass flux \(Q_\phi\), the volumetric flux \(V_1\), and the width of the plume at the beginning of the downwind gravity current \(L_0\)) is provided from the output of the eruption column model at the NBL (the mass flux \(Q_\phi\) and the volumetric flux \(V_0\)). We assume that \(Q_\phi\) is continuous at the connection of the eruption column and the downwind gravity current. On the other hand, \(V_1\) is considered to be substantially greater than \(V_0\) because of the entrainment of the ambient air around the top of the eruption column (e.g. Suzuki et al. 2005). To express this effect, we introduce a new parameter, \(\mu =V_1/V_0 (\ge 1)\); we evaluate the value of this parameter in a later section. The width at the beginning of the downwind gravity current (\(L_0\)) is determined such that it satisfies the continuity of the SMD (i.e. \(-\frac{dQ_\phi }{ds}\)) between the eruption column and the downwind gravity current at the NBL as follows,

We consider that the above method to connect the eruption column and the downwind gravity current models is a tentative one. This method does not fully describe the behavior around the top of a strong plume generated by an eruption with a large mass eruption rate, where the formation of an umbrella cloud plays a role in the tephra dispersal (Koyaguchi and Ohno 2001a, b; Bonadonna and Phillips 2003; Costa et al. 2013). A more sophisticated method that takes the dynamics of umbrella clouds into account needs to be developed in future studies.

Sedimentation model

The pyroclasts in the plumes are released from the side of the eruption column and the base of the downwind gravity current into the surroundings atomosphere (e.g. Woods and Bursik 1991; Bursik et al. 1992; Koyaguchi et al. 2009). To express a distribution of released particles along an eruption plume (i.e. SMD), we adopt the settling model for a turbulent fluid as follows (Martin and Nokes 1988),

Here, \(v_\phi\) is the terminal velocity for particles of grainsize \(\phi\), v is the velocity of the plume tangential to the plume axis (i.e. the direction of local coordinate s). \(\delta V\) is the control volume, and \(\delta S\) is the “sedimentation area” for the control volume \(\delta V\). For the terminal velocity \(v_\phi\), we use a formula proposed by Klawonn (2014), in which \(v_\phi\) is calculated on the basis of the formula of Cheng (2009) with a reduction of 5\(\%\) to account for the deviation from sphericity for actual volcanic particles.

The configurations of \(\delta V\) and \(\delta S\) are different between the eruption column model and downwind gravity current model. \(\delta V\) is expressed as a truncated cone with the height of ds for the eruption column model (Fig. 2a), whereas it is expressed as a rectangular box with the thickness h, the width L, and the length ds for the downwind gravity current model (Fig. 2b). The values of \(\delta V\) are calculated as

On the other hand, \(\delta S\) is defined as the projection of the lateral area of the control volume onto the ground for the eruption column model (Fig. 2a), whereas it is defined as a bottom area of the rectangular box for the downwind gravity current model (Fig. 2b). For vertical eruption columns, the projection of the lateral area of the control volume onto the ground is simply given by \(\pi (r_2^2-r_1^2)\) (Bursik et al. 1992), where \(r_1\) and \(r_2\) are radii of the bottom and top of the truncated cone. In this study, the effects of a bending eruption column are considered as is shown in Fig. 2a. In this case, \(\delta S\) are approximately calculated as,

All these geometrical parameters (i.e. \(\Theta\), \(\psi\) and dl) are determined from the results of the eruption column model.

Model parameters in NIKS-1D

NIKS-1D consists of the combination of the eruption column model, the downwind gravity current model, and the sedimentation model. Therefore, NIKS-1D inherits uncertainties associated with parameters included in each model (Fig. 1); the eruption column model includes the parameters about the entrainment process (i.e. \(\alpha\), \(\beta\)), and the gravity current model includes an uncertain parameters about flow geometry (i.e. \(\lambda\)). In this study, the values of these parameters are following the previous studies: \(\alpha\)=0.06 and \(\beta\)=0.2 (Suzuki and Koyaguchi 2015), \(\lambda\)=0.83 (Bursik et al. 1992). In addition to these parameters, NIKS-1D includes a new parameter at the connection of the eruption column model and the gravity current model; the parameter to parameterize the entrainment process at the top of the eruption column (i.e. \(\mu =V_1/V_0\)). In a later section, we discuss the effects of \(\mu\) on the output of NIKS-1D ("Run 2: the effects of μ" section). In addition, a formula to provide the recommended value of \(\mu\) is proposed on the basis of the observation data of the real eruption cases (The estimation of μ section).

Result

In this section, we present typical results for NIKS-1D defined in the previous section. An outline of the settings for four Runs is shown in Table 2. Firstly, we present the result of NIKS-1D with the typical condition as a reference to the other Runs ("Run 1: a reference" section). Secondly, we focus on the effect of \(\mu\) which is the new parameter introduced in NIKS-1D ("Run 2: the effects of μ" section). Finally, in order to describe the model results for the typical conditions, we perform parametric studies for an ambient wind ("Run 3: the effects of ambient wind" section) and a mass eruption rate ("Run 4: the effects of the mass eruption rate" section).

Run 1: a reference

As a reference, we set the input parameters of Run 1 such that its mass eruption rate is approximately equal to that of the Mt. St. Helens 1980 eruption. In this eruption, a downwind gravity current with an elongated parabolic outline was established after the initial formation of an umbrella cloud (Sparks et al. 1986; Bursik et al. 1992). Because NIKS-1D does not fully describe the effect of the formation of an umbrella cloud, it is considered that this eruption is roughly maximum eruption case to which NIKS-1D is applicable.

Figure 3a shows the centerline of an eruption column (red line). The eruption column is bent by an ambient wind, and connects to the downwind gravity current at the NBL (10309 [m], blue line in Fig. 3a). Figure 3b shows the top and side views of the downwind gravity current; the width of the downwind gravity current increases (red line), whereas the thickness decreases with distance from the vent (blue line).

Figure 4a shows the proportion of particles to reach the NBL for different grainsizes. Most of the pyroclasts finer than \(\phi\)=0 reach the NBL, whereas considerable parts of coarse particles tend to separate from the rising eruption column. Figure 4b shows the SMD from the eruption column as a function of altitude. The SMD for the coarse particles tends to be larger than that of fine particles, which accounts for the tendency that coarse particles separate from the rising eruption column in Fig. 4a. Figure 4b also shows that the SMD remarkably decreases with height below the altitude of 1km for any grainsize, whereas it is approximately constant above 1km except for \(\phi =-6, -5\). The remarkable decrease in the SMD below the altitude of 1km results from the fact that the radius of the eruption column is so small that \(\delta S/\delta V\) in Eq. (12) has large values at low altitudes.

Figure 4c shows the SMD from the downwind gravity current; the SMD for the coarse particles (e.g. \(\phi =-6, -4\)) decreases rapidly with distance from the vent, whereas changes of the trend in the SMD for the fine particles (e.g. \(\phi\)=4, 6) are moderate. The slight increase in the SMD with distance for fine particles results from the fact that the width of the gravity current increases with distance (see the red curve in Fig. 3b); the supply rate of particles per unit area from the downwind gravity current decreases with distance even for the fine particles.

The quantitative features observed in the results of Run 1 depend on the eruption and atmospheric conditions as well as model parameters. In the rest of this section, we investigate the effects of the new parameter (i.e. \(\mu\)), ambient wind, and the mass eruption rate.

Run 2: the effects of \(\mu\)

The volumetric flux at the beginning of the downwind gravity current (\(V_1\)) is generally greater than that of the NBL of the eruption column (\(V_0\)) due to entrainment at the top of the eruption column (e.g. \(\mu =V_1/V_0=3.6\) in Run 1 as recomended value; see Table 1). To evaluate this effect on NIKS-1D, we show the results for \(\mu =1\) (i.e. the case of no entrainment at the top of the eruption column) in the second columns of Figs. 3 and 4. Because the eruption column model does not include \(\mu\), the results concerning the eruption column model are identical between Run 1 (Figs. 3a, 4a,b) and Run 2 (Figs. 3c, 4d,e).

The parameter \(\mu\) significantly influences the results of the downwind gravity current (Figs. 3d and 4f). Figure 3d shows that the larger \(\mu\) leads to the larger width and thickness of the downwind gravity current; the solutions of the downwind gravity current model indicate that width L(s) and thickness h(s) are proportional to \(\sqrt{\mu s}\) and \(\sqrt{\mu /s}\), respectively. Figure 4f, on the other hand, shows that the SMD (i.e. \(-dQ_\phi/ds\) in Eq. (12)) for Run 2 decreases more rapidly than that for Run 1. This tendency is explained by the fact that the decay constant of \(Q_\phi\) in Eq. (12) is proportional to 1/h.

The above result indicates that \(\mu\) governs the width of the downwind gravity current. That implies that \(\mu\) can be estimated from the observed width of the downwind gravity current such as satellite observations. In a later section ("The estimation of μ" section), we will attempt to estimate \(\mu\) from the actual eruption cases, and propose recommended values \(\mu\) as a function of the mass eruption rate.

Run 3: the effects of ambient wind

The third columns of Figs. 3 and 4 show the results for the strong ambient wind (see Table 2). As pointed out in previous studies (e.g. Bursik 2001), the strong ambient wind leads to a bending plume and low plume height (Fig. 3e). The major causes of these trends are a supply of large horizontal momentum to the plume from the strong ambient wind and an increase in the dilution rate of the buoyancy of the plume due to large entrainment associated with the strong ambient wind. Figure 4g shows that the strong ambient wind slightly decreases the proportion of particles to reach the NBL. This tendency is explained by the fact that the SMD of the eruption column (Fig. 4h) slightly increases by the increase of \(\frac{\delta S}{\delta V}\) due to the plume bending.

The strong ambient wind also affects the shape of the downwind gravity current. Figure 3f shows that the strong ambient wind leads to the small width of the downwind gravity current, whereas the thickness slightly depends on the ambient wind speed. The solution of the gravity current model (Eqs. (9)(10)) indicates that the width L(s) is proportional to \(1/w_{\text {nbl}}\), whereas the thickness h(s) is almost independent of \(w_{\text {nbl}}\) except for the region around the source (i.e. \(s \sim s_0\) in Eq. (10)). The slight difference in the thicknesses between Run 1 and Run 3 results from slight changes in the Brunt-Väisälä frequencies N and volumetric fluxes \(V_1\) due to the change of the NBL. Figure 4i shows that the SMD from the downwind gravity current for Run 3 decreases with distance more slowly than that for Run 1. This tendency is explained by the fact that the decay constant of \(Q_\phi\) in Eq. (12) is proportional to \(1/w_{\text {nbl}}\). As a result, the strong wind has an effect to carry the pyroclast far from the vent, and tends to keep the mass flux in the downwind gravity current.

Run 4: the effects of the mass eruption rate

The fourth columns of Figs. 3 and 4 show the result for the small mass eruption rate (see Table 2). The small eruption rate leads to a low plume height (Fig. 3g) and small width and thickness of the downwind gravity current (Fig. 3h) when compared to the previous model Runs. These features of the downwind gravity current are accounted for by a small volumetric flux due to the small mass eruption rate; the solutions of the downwind gravity current model indicate that the width L(s) and thickness h(s) are proportional to \(\sqrt{V_1s}\) and \(\sqrt{V_1/s}\), respectively.

According to the column dynamics model (e.g. Morton et al. 1956; Woods 1988; Bursik 2001), the small mass eruption rate leads to a decrease in a rising velocity of the plume v. The small velocity of the plume (i.e. Run 4) results in the large SMD from the eruption column (see Eq. (12)). Due to this effect, the mass of coarse particles (\(\phi<\)3) that reaches the NBL is remarkably reduced for the small mass eruption rate compared with Run 1 (Fig. 4j). Figure 4l shows that the SMD for the small mass eruption rate decreases more rapidly than that of the large mass eruption rate. This tendency is explained by the fact that the decay constant of \(Q_\phi\) in Eq. (12) is proportional to 1/h.

Discussion

In this section, we firstly attempt to estimate the values of \(\mu\), which is the new parameter to connect the eruption column model and the downwind gravity current model. As shown in the previous section, \(\mu\) significantly influences the behavior of the SMD from the downwind gravity current. Therefore, an accurate estimation of \(\mu\) is indispensable for NIKS-1D to be applied to an operation for the real-time tephra-dispersal forecasting. After the estimation of \(\mu\), we explain the outline of a real-time tephra-dispersal forecasting system including NIKS-1D under development ("Outline of the real-time tephra-dispersal forecasting system" section).

The estimation of \(\mu\)

Previous studies indicate that the entrainment process around the top of the eruption column is characterized by complex interaction between the turbulence inside the eruption column and the upward/downward motion around the NBL, and that it strongly depends on the mass eruption rate (Bonadonna and Phillips 2003; Suzuki and Koyaguchi 2009; Suzuki et al. 2016; Devenish and Cerminara 2018). In this study, we attempt to parameterize \(\mu\) as a function of the mass eruption rate on the basis of observation data.

The analytical solution for the downwind gravity current model indicates that the volumetric flux \(V_1\) can be estimated from an observed width of the downwind gravity current, because the Brunt-Väisälä frequency N and wind velocity \(w_{{\rm nbl}}\) are known from the atmospheric conditions (see Eq. (9)). On the other hand, the volumetric flux \(V_0\) at the NBL is obtained from the eruption column model. Consequently, we can determine the value of \(\mu\) using the estimated \(V_1\) and \(V_0\). Here, we first estimate \(\mu\) from the observed width of the downwind gravity current for well-studied eruption events with different mass eruption rates (the Mt. St. Helens 1980 eruption (\(\dot{M} \sim 1\times 10^7\)[kg/s]; e.g., Sarna-Wojcicki et al. (1981); Sparks et al. (1986)) and the Ruapehu 1996 eruption (\(\dot{M}\sim 2\times 10^5\)[kg/s]; e.g., Mastin et al. (2009)), and show that the estimated \(\mu\) leads the SMD to be consistent with the SMD estimated from the tephra-fall deposits for those eruption events ("St. Helens eruption on May 18th 1980" and "Ruapehu eruption on June 17th 1996" sections). On the basis of these results, we attempt to parameterize \(\mu\) as a function of the mass eruption rate (Relationship between parameter μ and mass eruption rate Ṁ section).

St. Helens eruption on May 18th 1980

For this eruption case, the plume established an elongated parabolic outline at 10:45 PDT on the 18th May 1980; the width of the plume at the distance x was estimated as \(L(x)\sim 588 \sqrt{x}\); the effect of the source position of the downwind gravity current is assumed to be sufficiently small compared with the whole shape of the downwind gravity current (Bursik et al. 1992). By substituting this estimation of L(x) to Eq. (9), we can obtain \(\mu\)=3.6 with the Brunt-Väisälä frequency N=4.2\(\times 10^{-2}\) [1/s], the wind velocity \(w_{\text {nbl}}\)=25 [m/s] and the volumetric flux \(V_0\)=0.87 [km\(^3\)/s] at the NBL (\(\sim\)12.6 [km] asl) calculated by the eruption column model in NIKS-1D with the settings of the Table 1 and an atmospheric condition at the time of the eruption. Here, the atmospheric condition is obtained by interpolating the gridded data of the global atmospheric reanalysis (JRA-55C; Kobayashi et al. 2014) to the location of St. Helens for 11:00 PDT on the 18th May 1980.

Using the estimated \(\mu\) (=3.6), the SMD can be calculated from NIKS-1D. On the other hand, the SMD can be estimated independently from the tephra-fall deposits (see Fig. 3b in Bursik et al. 1992). Figure 5 shows that the SMD from NIKS-1D is roughly consistent with the SMD based on the tephra-fall deposits.

Ruapehu eruption on June 17th 1996

For this eruption case, the plume established an elongated parabolic outline at 08:30 NZT on the 17th June 1996 (Bonadonna et al. 2005); the width of the plume at the distance x can be approximately estimated as \(L(x) \sim 158\sqrt{x}\) from satellite imagery (Bonadonna et al. 2005). By substituting this estimated L(x) to Eq. (9), we can obtain \(\mu\)=4.7 with N=3.2\(\times 10^{-2}\) [1/s]), \(w_{\text {nbl}}\)=28 [m/s], and \(V_0\)=0.079 [km\(^3\)/s] at the NBL (\(\sim\)6.4 [km] asl) calculated by the eruption column model in NIKS-1D with the settings of Table 1 and an atmospheric condition at the time of the eruption. Here, the atmospheric condition is obtained by interpolating the gridded data of the global atmospheric reanalysis (JRA-55C; Kobayashi et al. 2014) to the location of Mt. Ruapheu for 06:00 NZT on the 17th June 1996.

Using estimated \(\mu\) (=4.7), the SMD can be calculated from NIKS-1D. On the other hand, the SMD can be estimated independently from the tephra-fall deposits (Klawonn 2014). Figure 6 same as Fig. 5 shows that the SMD from NIKS-1D is roughly consistent with the SMD based on the tephra-fall deposits.

Relationship between parameter \(\mu\) and mass eruption rate \(\dot{M}\)

The estimation of \(\mu\) for the two eruptions suggests that the parameter \(\mu\) decreases with the mass eruption rate. Figure 7 shows that the tendency of the variation in \(\mu\) estimated from the two examples is crudely parameterized by a relationship between \(\mu\) and \(\log _{10}(\dot{M}({\text {kg/s}}))\) as

$$\begin{aligned} \mu = a - b \log _{10} (\dot{M}) \end{aligned}$$

(18)

where a=8.1 and b=0.65. Although the strong plume is out of scope in NIKS-1D, this parameterization seems to be consistent with the strong plume generating an umbrella cloud; the 3-D simulation of the umbrella cloud for the 1991 eruption of Mt. Pinatubo indicate that the value of \(\mu (=V_1/V_0)\) can be as small as 2 when the mass eruption rate is as large as \(10^9\) [kg/s] (Suzuki and Koyaguchi 2009).

In this section, \(\mu\) which is needed to execute NIKS-1D was parameterized by the relationship between \(\mu\) and \(\dot{M}\) in Eq. (18). We emphasize here that this relationship is tentative one and that further observations, as well as theoretical studies about the entrainment process around the top of the plume, will improve the parameterization of \(\mu\) in the future.

Outline of the real-time tephra-dispersal forecasting system

In this section, we briefly explain the outline of the real-time tephra-dispersal forecasting system in the JMA using NIKS-1D. Figure 8 shows a flowchart of the real-time tephra-dispersal forecasting system under development. This system is composed of three steps as follows. The first step is the estimation of the mass eruption rate \(\dot{M}\) by the inversion based on NIKS-1D from the observed plume height with various sources (e.g. satellite analysis and pilot reports, etc.). The second step is a calculation of the SMD by NIKS-1D from the mass eruption rate \(\dot{M}\) estimated in the first step. The SMD calculated in this step is used as an initial condition for a tephra-dispersal model in the next step; the effect of variation in sedimentaion area (\(\delta S\)) is taken into consideration as necessary. The third step is a calculation of the tephra-dispersal forecasting by a tephra-dispersal model (JMA-ATM, Shimbori and Ishii 2021). After these steps, the information about the tephra-dispersal forecast is sent to local governments and JMA branch offices around the volcano, and posted on the JMA website. Here, the first step and second step have computational costs of a few seconds, each. The third step has computational costs of a few tens of seconds. These computational costs are much shorter than the time-scale of the acquisition of the column-height data. In this sense, we consider that the information about the tephra-dispersal forecast is expected to be issued in real-time.

The model developed here (i.e. NIKS-1D) is concerned with the first and second steps. In the system used in JMA at present, the mass eruption rate and the SMD are estimated on the basis of two empirical models: the power-law for the mass eruption rate (Morton et al. 1956; Mastin et al. 2009) and the empirical function for the SMD (Suzuki 1983), respectively. We propose that the mass eruption rate and the SMD are estimated by physics-based models. The accuracy of NIKS-1D is expected to be improved through further observations as well as future theoretical studies about plume dynamics, which will improve both estimations of the mass eruption rate and the SMD. For example, NIKS-1D uses petrological and geophysical data at the vent (e.g. initial temperature and exit velocity) as model parameters. This means that constraints of these model parameters by additional petrological and geophysical observations will lead to improvements of the tephra-dispersal forecasting based on NIKS-1D.

Summary and future work

In this study, we proposed the one-dimensional plume model (NIKS-1D) for the real-time tephra-dispersal forecasting system. This model was implemented as a combination of the eruption column model, the downwind gravity current model, and the sedimentation model. The combination of these three models makes it possible to calculate the SMD from the eruption plume which is used as an initial condition for the tephra-dispersal model. To comprehend the characteristics of the SMD from NIKS-1D for typical conditions, we studied the effects of ambient wind and mass eruption rate. On the other hand, the combination of the eruption column model and the downwind gravity current model requires a new parameter \(\mu\) which is related to the entrainment around their connection. We estimated \(\mu\) for real eruption cases, and propose recommended values \(\mu\) as a function of the mass eruption rate. We also discussed the role of NIKS-1D in the real-time tephra-dispersal forecasting system under development; in this system, NIKS-1D is used to estimate a mass eruption rate from an observed plume height, and provide the SMD for the tephra-dispersal model.

The advantage of NIKS-1D is that it is based on the plume dynamics. Therefore, further theoretical understanding and observations of plume dynamics can lead to improvements of NIKS-1D while used in an operation. For example, the basic idea of NIKS-1D is naturally extended to the umbrella clouds which spread radially at the NBL just above the vent as a turbulent gravity current (e.g. Koyaguchi and Ohno 2001a, b; Bonadonna and Phillips 2003; Suzuki and Koyaguchi 2009; Costa et al. 2013; Johnson et al. 2015); substituting the umbrella cloud model for the downwind gravity current model is expected to improve of NIKS-1D, particularly for large-scale eruptions.

Finally, we point out the limitations of the real-time tephra-dispersal forecasting system including NIKS-1D. NIKS-1D does not consider some processes that can affect the plume dynamics and the tephra-dispersal such as particle aggregation (e.g. Brown et al. 2012; Macedonio et al. 2016) and water phase change (e.g. Mastin 2007). Some of the hypotheses used in the eruption column model (e.g., the entrainment hypothesis) are still inconclusive in a quantitative sense (e.g. Costa et al. 2016). These limitations of NIKS-1D can introduce bias or uncertainties in the result of the tephra-dispersal model via the SMD. The estimation of the mass eruption rate from observed plume heights is another possible source of uncertainty. Preliminary reports of plume heights based on monitoring cameras and pilot reports often have a large uncertainty, which causes an uncertainty in the input parameter for NIKS-1D. In an operational environment, the information about these uncertainties will be as vital as the forecast itself.

Availability of data and materials

Not applicable.

Abbreviations

JMA:

Japan Meteorological Agency

NBL:

Neutral Buoyancy Level

SMD:

Source Magnitude Distribution

References

Bonadonna C, Phillips JC (2003) Sedimentation from strong volcanic plumes. J Geophys Res Solid Earth 108(B7):2340

Bonadonna C, Phillips JC, Houghton BF (2005) Modeling tephra sedimentation from a Ruapehu weak plume eruption. Geophys Res Solid Earth 110(B8). https://doi.org/10.1029/2004JB003515

Brown RJ, Bonadonna C, Durant AJ (2012) A review of volcanic ash aggregation. Phys Chem Earth A/B/C 45–46:65–78

Bursik MI, Sparks RSJ, Gilbert JS, Carey SN (1992) Sedimentation of tephra by volcanic plumes: I. Theory and its comparison with a study of the Fogo A plinian deposit, Sao Miguel (Azores). Bull Volcanol 54:329–344

Casadevall JT (1994) Volcanic ash and aviation safety; proceedings of the first international symposium on volcanic ash and aviation safety. US Geol Surv Bull 2047:107–145

Cheng NS (2009) Comparison of formulas for drag coefficient and settling velocity of spherical particles. Powder Technol 189:395–398

Clarkson RJ, Majewicz EJ, Mack P (2016) A re-evaluation of the 2010 quantitative understanding of the effects volcanic ash has on gas turbine engines. Proc Inst Mech Eng G J Aerosp Eng 230(12):2274–2291

Costa A, Folch A, Macedonio G (2013) Density-driven transport in the umbrella region of volcanic clouds: Implications for tephra dispersion models. Geophysical Research Letters 40:4823–4827

Costa A, Suzuki YJ, Cerminara M, Devenish BJ, Ongaro TE, Herzog M, Van Eaton AR, Denby LC, Bursik M, de’ Michieli Vitturi M, Engwell S, Neri A, Barsotti S, Folch A, Macedonio G, Girault F, Carazzo G, Tait S, Kaminski E, Mastin LG, Woodhouse MJ, Phillips JC, Hogg AJ, Degruyter W, Bonadonna C (2016) Results of the eruptive column model inter-comparison study. J Volcanol Geotherm Res 326:2–25

Folch A (2012) A review of tephra transport and dispersal models: Evolution, current status, and future perspectives. J Volcanol Geotherm Res 235–236:96–115

Girault F, Carazzo G, Tait S, Ferrucci F, Kaminski E (2014) The effect of total grain-size distribution on the dynamics of turbulent volcanic plumes. Earth Planet Sci Lett 394:124–134

Girault F, Carazzo G, Tait S, Kaminski E (2016) Combined effects of total grain-size distribution and crosswind on the rise of eruptive volcanic columns. J Volcanol Geotherm Res 326:103–113

Guffanti M, Schneider DJ, Wallace KL, Hall T, Bensimon DR, Salinas LJ (2010) Aviation response to a widely dispersed volcanic ash and gas cloud from the August 2008 eruption of Kasatochi, Alaska, USA. J Geophys Res Atmos 115(D2). https://doi.org/10.1029/2010JD013868

Hasegawa Y, Sugai A, Yo H, Yu H, Saito S, Shimbori T (2015) Improvements of volcanic ash fall forecasts issued by the Japan Meteorological Agency. J Appl Volcanol 4(2)

Hewett TA, Fay JA, Hoult DP (1971) Laboratory experiments of smokestack plumes in a stable atmosphere. Atmos Environ 5:767–789

Klawonn M, Frazer LN, Wolfe CJ, Houghton BF, Rosenberg MD (2014) Constraining particle size-dependent plume sedimentation from the 17 June 1996 eruption of Ruapehu Volcano, New Zealand, using geophysical inversions. J Geophys Res Solid Earth 119:1749–1763

Kobayashi C, Endo H, Ota Y, Kobayashi S, Onoda H, Harada Y, Onogi K, Kamahori H (2014) Preliminary Results of the JRA-55C, an Atmospheric Reanalysis Assimilating Conventional Observations Only. SOLA 10:78–82

Koyaguchi T, Ohno M (2001) Reconstruction of eruption column dynamics on the basis of grain size of tephra fall deposits: 1. Methods. J Geophys Res Solid Earth 106:6499–6512

Koyaguchi T, Ohno M (2001) Reconstruction of eruption column dynamics on the basis of grain size of tephra fall deposits: 2. Application to the Pinatubo 1991 eruption. J Geophys Res Solid Earth 106:6513–6533

Koyaguchi T, Ochiai K, Suzuki YJ (2009) The effect of intensity of turbulence in umbrella cloud on tephra dispersion during explosive volcanic eruptions: Experimental and numerical approaches. J Volcanol Geotherm Res 186:68–78

Liang Y, Xu J (2021) The impact of volcanic ash on the safety of aviation industry: review of china’s current situation. Geol Soc London Spec Publ 510:253–262

Macedonio G, Costa A, Scollo S, Neri A (2016) Effects of eruption source parameter variation and meteorological dataset on tephra fallout hazard assessment: example from vesuvius (italy). J Appl Volcanol 5:5

Mastin LG (2007) A user-friendly one-dimensional model for wet volcanic plumes. Geochem Geophys Geosyst 8(3)

Mastin LG, Guffanti M, Servranckx R, Webley P, Barsotti S, Dean K, Durant A, Ewert JW, Neri A, Rose WI, Schneider D, Siebert L, Stunder B, Swanson G, Tupper A, Volentik A, Waythomas CF (2009) A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. J Volcanol Geotherm Res 186:10–21

Morton BR, Taylor GI, Turner JS (1956) Turbulent gravitational convection from maintained and instantaneous sources. Proc R Soc Lond Ser A Math Phys Sci 234:1–23

Sarna-Wojcicki AM, Shipley S, Waitt R, Dzurisin D, Wood S (1981) Areal distribution, thickness, mass, volume, and grain size of air-fall ash from the six major eruptions of 1980. U.S. Geological Survey Professional Paper 1250:577–600. https://doi.org/10.3133/pp1250

Shimbori T, Ishii K (2021) Design of the Japan Meteorological Agency Atmospheric Transport Model. Technical Reports of the Meteorological Research Institute 84

Simpson JE (1987) Gravity Currents: In the Environment and the Laboratory. Ellis Horwood Limited, Chichester

Sparks RSJ, Moore JG, Rice CJ (1986) The initial giant umbrella cloud of the may 18th, 1980, explosive eruption of mount st. helens. J Volcanol Geotherm Res 28:257–274

Suzuki, T (1983) A theoretical model for dispersion of tephra. In: Shimozuru D Yokoyama I (eds)Arc Volcanism: Physics and Tectonics, 93–113. Terra Scientific Publishing Company (TERRAPUB), Tokyo. https://doi.org/10.1186/s13617-016-0045-2

Suzuki YJ, Koyaguchi T (2009) A three-dimensional numerical simulation of spreading umbrella clouds. J Geophys Res Solid Earth 114:B03209. https://doi.org/10.1029/2007JB005369

Suzuki YJ, Koyaguchi T, Ogawa M, Hachisu I (2005) A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model. J Geophys Res Solid Earth 110

Suzuki YJ, Koyaguchi T (2015) Effects of wind on entrainment efficiency in volcanic plumes. J Geophys Res Solid Earth 120:6122–6140

Suzuki YJ, Costa A, Cerminara M, Esposti Ongaro T, Herzog M, Van Eaton AR, Denby LC (2016) Inter-comparison of three-dimensional models of volcanic plumes. J Volcanol Geotherm Res 326:26–42

This work was supported by the “Integrated Program for Next Generation Volcano Research and Human Resource Development” by the Ministry of Education Culture, Sports, Science and Technology-Japan (MEXT).

Funding

Not applicable.

Author information

Authors and Affiliations

Meteorological Research Institute, Japan Meteorological Agency, Ibaraki, Japan

Kensuke Ishii

Earthquake Research Institute, the University of Tokyo, Tokyo, Japan

Akira Nishijo, Takehiro Koyaguchi & Yujiro J. Suzuki

Current address: Seismological and Volcanological Department, Japan Meteorological Agency, Tokyo, Japan

NIKS-1D was developed with the cooperation of the Japan Meteorological Agency and Earthquake Research Institute, the University of Tokyo; all the authors contributed to the establishment of the basic idea of the present model. An early version of the source code was developed as a master thesis project by AN under the supervision of TK and YJS. KI developed the source code for this paper, and implemented the inversion algorithm to NIKS-1D.

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Ishii, K., Nishijo, A., Koyaguchi, T. et al. A physics-based source model for real-time tephra-dispersal forecasting for weak eruption plumes.
J Appl. Volcanol.11, 15 (2022). https://doi.org/10.1186/s13617-022-00127-w