We used the eruptive events listed in the IVESPA (Independent Volcanic Eruption Source Parameters) database from Aubry et al. (2021) to test the modelling approaches. These entries include not only “complete” eruptions, but also individual phases to reflect the transitions between eruptive styles or column regimes in more complex eruptions. From the 134 eruptive events listed, 130 entries provide plume top heights, vent level altitude, duration of eruption and independent estimates (i.e., not inferred by reversed plume height modelling) of total erupted mass (M _IVESPA ). 64 of these entries provide information on the mass uncertainties, while no such information is given for the other 66 eruptive events. For this study, all top heights were converted to plume heights above vent level.

We used three different subsets of IVESPA for our analyses:

• “All” data (sample size N = 130, see Supplementary Table S1): this dataset includes all eruptive events listed in IVESPA with given plume heights, event durations and total erupted mass. For cases in which no information on the uncertainties of the total erupted mass was provided, an error of ±50% was assumed. This value represents the approximate average of mass uncertainties listed in IVESPA.

The values for wind speeds V ¯ , wind shear W s ∼ , atmospheric temperatures at the vent T _ao , the height-averaged buoyancy frequencies N ¯ , the duration of the eruption t and the magmatic temperature T ₀ used in this study were taken directly from the IVESPA database (Aubry et al., 2021). All weather parameters taken from IVESPA are based on ERA5 re-analysis data (Hersbach et al., 2018; 2020). When no temperature data were available for T ₀ , a value between 950°C and 1,100°C was assumed, depending on the composition of the magma (see Supplementary Tables S1–S3), with an uncertainty of ±50°C. Following the considerations in Dürig et al. (2023), we used (1,250 ± 100) J kg⁻¹ K⁻¹ for the specific heat capacity C ₀ of the ash and (998± 20) J kg⁻¹ K⁻¹ for the specific heat capacity of the atmosphere C _ao .

All statements on statistical significance are based on the findings from two-tailed Student’s (Student, 1908) and Welch’s (Welch, 1947) t-tests. Which type of t-test was applied depends on the homogeneity of the datasets’ variances (Davis, 2002; Dürig et al., 2021). To find out if the tested datasets were homogeneous or heterogeneous, we used Levene tests (Levene, 1960).

The six models tested in this study are named after their first authors and are described below.

The Wilson model is a simple numerical 0D model based on the buoyant plume theory by Morton et al. (1956), which estimates the mass eruption rate MER by: M E R W i l = h / c 4 (1) with h, the height of the plume top (in m), as the only input parameter. The constant c is 236 m (s/kg)^1/4.

The Sparks (Sparks et al., 1997) and Mastin (Mastin et al., 2009) models are 0D models, which calibrated both c and the exponent within Eq. 1 against data from historic eruptions to link plume top height h with the MER. The relationship found by Sparks et al. (1997) is: M E R S p a = ρ ∙ h / c 3.86 (2) where ρ is the dense-rock equivalent (DRE) density of the tephra within the plume and c is a constant of 1,670 m (s/kg)^1/3.86.

The Mastin model is given by: M E R M a s = ρ ∙ h / c 4.15 (3) with c being 2000 m (s/kg)^1/4.15 (Mastin et al., 2009).

Another empirical model studied is the Gudmundsson model, which was calibrated for the different eruptive phases of the 2010 Eyjafjallajökull eruption, Iceland (Gudmundsson et al., 2012). M E R G u d = ρ ∙ a ∙ k I ∙ h a v g + h max / c 4.15 (4)

This model uses the same constant c as the Mastin model. In addition, it introduces a dimensionless constant a which is calibrated to be 0.0564. In contrast to other models, the Gudmundsson model introduces an additional scale factor k _I , which is dependent on the eruptive conditions, such as wind speed and/or eruptive style (Gudmundsson et al., 2012). For simplification, we used a value of 1.6 for all eruptions. This is the suggested optimal value for modelling the complete 2010 Eyjafjallajökull eruption (Gudmundsson et al., 2012). The Gudmundsson model is also the only model from those analysed in this study, which uses two different measures of plume height as input parameters: the average plume top height h _avg and the maximum plume top height h _max .

Being an explicitly wind-affected model, the Degruyter model (Degruyter and Bonadonna, 2012) was developed by combining the buoyant plume theory of Morton et al. (1956) and its modification by Hewett et al. (1971). It estimates the mass eruption rate by: M E R Deg = π ρ a 0 g ′ 2 5 / 2 α 2 N ¯ 3 z 1 4 H C 4 + β 2 N ¯ 2 V ¯ 6 H C 3 (5)

The subscripts a and 0 refer respectively to the atmosphere and volcanic source vent height. Therefore, the density of the atmosphere at vent level is given by ρ _a0 . The plume height-averaged wind speed is defined as V ¯ , and the height-averaged buoyancy frequency as N ¯ . The radial and wind entrainment coefficients are denoted α and β. The maximum non-dimensional height resulting from numerical integration of the governing equations of Morton et al. (1956) is given by z₁ = 2.8 (Degruyter and Bonadonna, 2012). The reduced gravity at the source g’ is defined as: g ′ = g C 0 T 0 − C a 0 T a 0 C a 0 T a 0 (6) with T ₀ being the magma temperature at the source, C ₀ the specific heat capacity of the magma, T _ao the atmospheric temperature at the vent and C _ao the specific heat capacity of the atmosphere at the vent.

The assessment of the plume height H _c used as model input parameter depends on the type of plume featured in the eruption and can in many cases not be simply assumed to be identical with the top height h (Mastin, 2014; Devenish, 2016; Dürig et al., 2023). This is due to the fact that the wind-affected plume models relate MER to the plume maximum centreline height H _c . If the plume trajectory is purely vertical, then H _c coincides with the top plume height h, otherwise the two heights do not coincide and can potentially lead to significantly different estimates of MER (Mastin, 2014; Devenish, 2016; Dürig et al., 2023).

The sixth and last model we tested was the Woodhouse model, which is given by a relationship derived from the application of the 1D model developed by Woodhouse et al. (2013) with β = 0.9. The Woodhouse model estimates the MER using: M E R W o o = 1 0.318 H C 1 + 4.266 W s ∼ + 0.3527 W s ∼ 2 1 + 1.373 W s ∼ 3.953 (7) where the wind shear from the ground to a reference height H ₁ is described by the parameter W s ∼ , according to: W s ∼ = 1.44 W s = 1.44 V 1 N ¯ H 1 (8)

V ₁ is the wind speed at reference height H ₁ = H _c (Woodhouse et al., 2013). Like the Degruyter model, Woodhouse uses H _c as main model input parameter.

To distinguish the plume types, we use the scaled parameter Π, which describes the relative influence of buoyancy and cross-wind of windspeed V ¯ on the plume dynamics and is characterized by (Degruyter and Bonadonna, 2012; Bonadonna et al., 2015): Π = N ¯ h 1.8 ∙ V ¯ ∙ α β 2 (9)

For consistency with published results in literature (Degruyter and Bonadonna, 2012; Bonadonna et al., 2015; Scollo et al., 2019; Aubry et al., 2023; Dürig et al., 2023), we used the entrainment coefficients α = 0.1 and β = 0.5 for the computation of Π. Following the approach applied in previous studies on paroxysms at Etna (Scollo et al., 2019) and on the ash plume of Eyjafjallajökull 2010 (Dürig et al., 2023), Eq. 9 can be used to discriminate three plume types (see Figure 1 in Mastin (2014)):

(i) buoyancy-dominated plumes, which rise vertically and are often referred to as “strong plumes” (Carey and Sparks, 1986; Mastin, 2014). They are characterized by Π > 0.5, and since the maximum elevation of the plume’s centreline coincides in such a case with the top height h, we can use this value for H _c .

For simplification, we followed the approach introduced by Scollo et al. (2019) and applied plume height corrections only for wind-dominated plumes, using an approximation for the plume radius R suggested by Devenish (2016). Thus H _c was obtained by: H C = h ∙ 1 − β β + 1 ;   i f   Π   <   0.1 (10) H C = h ;   i f   Π ≥ 0.1 (11)

To test the effect of the correction on the outcome, we also computed the “uncorrected” model outcome for Degruyter and Woodhouse, by assuming H _c = h for all cases. The uncorrected model results are marked by the suffix “uncor.”.

The Degruyter model (Eq. 5) and the plume height correction (Eq. 10) require knowledge of the entrainment coefficients. Based on theoretical considerations (e.g., Turner, 1986; Carazzo et al., 2006; Papanicolaou et al., 2008) and experimental findings (Dellino et al., 2014), the radial entrainment coefficient α is known to be ∼0.1. In contrast, the wind entrainment coefficient β is more challenging. The majority of suggested values for β gravitate around 0.5 (Huq and Stewart, 1996; Contini et al., 2011; Michaud-Dubuy et al., 2020), and this is also the value used in the studies of Degruyter and Bonadonna (2012) and Scollo et al. (2019). There are, however, also other studies that suggest different values for β, ranging between 0.1 and 1.0 (Bursik, 2001; Suzuki and Koyaguchi, 2012; Woodhouse et al., 2013). For example, a recent study on the ash plume for Eyjafjallajökull 2010 found that the prediction accuracy of the Degruyter model can be significantly improved when using a value of β = 0.3 (Dürig et al., 2023). Here, we therefore tested two values for β: 0.3 and 0.5 when computing Equations 5, 10. We note that the choice of β in Eq. 10 also influences the results of the Woodhouse model and is therefore always reported along with the results.

To find out which of the tested approaches is most accurate in estimating the MER, we computed the predicted total erupted mass M _model , according to: M m o d e l = M E R m o d e l ∙ t (12) where t is the duration of the eruptive event. The uncertainties of M _model were estimated by using the endmember values (i.e., minimum and maximum values) for h and t, when computing MER _model and Eq. 12. This results in a range of predicted mass eruption rate MER _model and predicted erupted mass M _model , constrained by the endmembers MER _min and MER _max , and M _min and M _max , respectively. Since such an approach is not applicable for the Gudmundsson model (Eq. 4), only the best estimate for MER and M _model was computed for this model.

For each eruptive event, the predicted mass is compared to the independently assessed total erupted mass (M _IVESPA ) from the IVESPA database (Aubry et al., 2021), and the number of matches between M _model and M _IVESPA is counted. The number of “successful” predictions is then referred to the total number of eruptive events tested, and this ratio is termed “success rate” S. The outcome of this approach is strongly dependent on the mathematical definition of the conditions, under which a prediction provides a match and therefore qualifies as “successful”. In this study we tested two different matching criteria (see Figure 1):

In addition to testing the six models individually, we also examined linear combinations of these models. The strategy to combine the models via model weight factors W _model is, for example, used by the software REFIR for real-time prediction of the MER (Dürig et al., 2018; Dioguardi et al., 2020). According to this “mixing” strategy, the mass eruption rate is assessed by: M E R m i x = ∑ m o d e l W m o d e l ∙ M E R m o d e l (16) where the sum of the six model weight factors is 1. The range of predicted mass M _mix was computed with Eq. 12, analogue to the procedure for individual models.

A Matlab^® script (see Supplementary File S1) was applied that iteratively computes M _mix for all possible combinations of model weight factors with step size 0.01. In a second step, for the selected matching criterion, the weight factor combinations with the highest success rate S were found. In cases where several weight factor solutions exist, we report the combination that i) involves the lowest number of models and ii) prioritizes the simple plume models (Wilson, Sparks, Mastin and, where applicable, Gudmundsson). This prioritization is chosen to minimize the uncertainties that result from combined model uncertainties and from the errors introduced by the additional input parameters required for the explicitly wind-affected models (Degruyter and Woodhouse).

We also examined if S can be further increased by “tailoring” the choice of the model and/or model weight factors to the eruptive conditions. For this purpose, we split the dataset introduced above as “all” into smaller subsets, grouped by:

• Magmatic composition: “basalt”, “andesitic basalt”, “andesite” and “dacite/rhyolite”;

All grouping variables were provided by, or computed on the basis of the IVESPA database (Aubry et al., 2021). As illustrated in Figure 2, this results in the following types of model solutions for optimized S:

i. Optimal global model: representing the model that provides the highest success rate for a tested dataset (“all”, “precise” and “control”);

The estimates resulting from all tested models are shown in Figure 3 for all tested data subsets. In Figure 3A, the modelled erupted mass is plotted over the measured erupted mass M _IVESPA for all 130 events (dataset “all”). This includes cases with unknown M _IVESPA errors, for which a 50% uncertainty was assumed. An example of the uncertainties on the observed and simulated total erupted mass for the Wilson model is shown as error bars in Figure 3B. Figures 3C,D are analogue to Figure 3A, but they show the model results for the data subsets “precise” and “control”, respectively. The latter excludes events that were used for the calibration of the models Mastin and Sparks, while the former only includes data points with known mass errors.

Table 1 presents the success rates S resulting from applying the “best fit” matching criterion to the tested datasets. The values for S range from 7.8% (uncorrected Degruyter model with β = 0.5 and the “precise” subset) to 39.1% (Mastin model with the “precise” subset). When using the FMER definition of Eq. 15 for the “best fit” criterion, the success rates for Sparks, Mastin, uncorrected Woodhouse, plume-height corrected Woodhouse with β = 0.5 and plume-height corrected Degruyter with β = 0.3 are higher than when using the CMER definition of Eq. 14. Using the FMER definition when applying Mastin always leads to a significant increase of S. This finding is in agreement with the optimal MER prediction strategy reported for the different eruptive phases of Eyjafjallajökul 2010 (Dürig et al., 2022) and is further supported by similar results for Mastin with the data subsets “precise” and “control”, where the use of FMER instead of CMER leads to an increase in S from 26.6% to 39.1% and from 28.6% to 30.5%, respectively.

The success rates S for the “best range” criterion are listed in Table 2. Since the “best range” criterion is less strict than the “best fit” criterion, the success rates are larger, ranging from 54.3% (uncorrected Degruyter with β = 0.5 for dataset “control”) to 84.4% (plume-height corrected Degruyter with β = 0.3 for dataset “precise”). This large variation in success rates demonstrates the importance of applying plume-height corrections, when using the Degruyter model. For both tested matching criteria, the highest success rates are achieved, when a plume height correction is applied, and a wind entrainment coefficient β = 0.3 is assumed. This is in agreement with the findings of Dürig et al. (2023).

Also, the success rates of the other wind-affected plume model, Woodhouse, are improved when using corrected plume heights (Table 2). An important difference to the optimal correction for the Degruyter model, however, is that Woodhouse (which does not explicitly depend on β but is implicitly based on the assumption β = 0.9, i.e., for this model the variation of β only affects the top plume height correction) yields higher success rates when using β = 0.5 for the plume height correction with Eq. 10. These outcomes were confirmed for all tested datasets and matching criteria. For the rest of this study, we therefore exclusively focus on model versions with highest success rates, which are the plume height corrected models Degruyter with β = 0.3 and Woodhouse with β = 0.5.

The optimal global model (Figure 2A) depends on the choice of the matching criterion and partly on the tested dataset. When using the “best fit” criterion, Mastin yields the highest success rates for all tested datasets (accompanied by the other empirical 0D model Sparks for dataset “control”). When applying the “best range” criterion, Mastin is the optimal global model for dataset “all” (S = 80.0%), but for the dataset “precise”, Mastin is equally successful to plume-height corrected Woodhouse with β = 0.5 (S = 82.8%), and even excelled by (plume-height corrected) Degruyter with β = 0.3 (S = 84.4%). For dataset “control”, the optimal global models are Mastin and (plume-height corrected) Degruyter with β = 0.3 (S = 77.1%).

Table 3 shows the success rates of models tested for subsets of the dataset “all” grouped by specified eruptive conditions. For example, when the eruptive events from the IVESPA database are grouped by composition, we find that for basaltic eruptions the optimal situation-specific plume model (i.e., the model with highest S) is Mastin, whereas for eruptions producing andesitic basalts it is Woodhouse. Both models also show the largest success rates in predicting mass eruption rates of dacitic/rhyolitic eruptions. For andesitic eruptions, the most successful plume models are Wilson, Sparks, Mastin and Degruyter with success rates of 75.0%.

The success rates listed in Table 3 under “optimal selection” for each group variable were obtained by selecting the optimal situation-specific models and dividing the sum of the accurate predictions by 130, which is the total number of events tested. For example, when selecting Degruyter for wind-dominated plumes, Wilson for intermediate plumes and Woodhouse for buoyancy-dominated plumes, the success rate for this composition-specific optimal model selection strategy is 81.5%, hence slightly higher than the 80% of the global optimal model. Higher gains in S can be achieved, if one of the other tested grouping variables is used as discriminatory eruptive condition for selecting the optimal model. The grouping variable with the highest success rates for the optimal situation-specific model approach turns out to be eruptive style, with a gain in S of 4.6%, compared to the optimal global model predictions. The optimal style-specific models are Woodhouse and Degruyter for magmatic or phreatomagmatic eruptions, respectively, while it is Wilson for eruptions of unknown style. For modelling phreatic eruptions, the success rates appear to be independent from the choice of the model, but the sample size of 5 is too low for any further general inference from this result.

In Table 4, model weight factors W are presented that lead, with Eq. 16, to the highest success rates S for the tested datasets and matching criteria. These weight factors characterize solutions for the optimal global mix (see Figure 2C). For example, for the dataset “all” with “best range” matching criterion, a success rate of 82.3% is obtained for mixing the models with W _Wilson = 0.20, W _Sparks = 0.20, W _Degruyter = 0.15 and W _Woodhouse = 0.45. We note that in contrast to this particular case, the solutions presented are generally not unique, which means that a large numbers of different weight factor solutions can lead to the same given maximum S. The number of solutions is listed as Y in Table 4. For “best fit”, the number of equivalent optimal global mix solutions Y ranges from 3 to 10, while for “best range”, Y ranges from 1 to 97.

Similar to the approach for the optimal situation-specific models (Table 3), it was examined for which grouping variable an optimal selection of weight factors would lead to the highest success rates. Testing dataset “all”, we found that if we distinguish by eruptive style and always select the optimal style-specific combination of weight factors, this results in a success rate of 38.5% with the “best fit” criterion and in 84.6% with the “best range” criterion. The latter value is the same as found for the optimal selection of style-specific models. In other words, within the analyzed step size of 0.01 no weight factor combination was found that would further improve S compared to the solutions listed for optimal style-specific model selection in Table 3. For other grouping variables, however, the combination of various models can lead to higher success rates (all solutions can be found in Supplementary Table S4).

When sorting the dataset by plume type, composition and duration, we find with “best range” optimal situation-specific mix solutions with success rates of 84.6%, 83.8% and 84.6%, respectively. With 85.4%, the highest success rates are obtained when sorting the dataset by eruptive strength, leading to optimal MER-specific mix solutions. An example of optimal composition-specific weight factor combinations is presented in Table 5.

The maximum success rates for all tested model solutions, datasets and matching criteria are presented in Figure 4. From the datasets tested, the highest success rates for all model solution types are found for the dataset labelled “precise”. The same observation is made for the success rates of all models (see Tables 1, 2) or mix solutions (see Table 4). Assuming an error of 50% for cases of unknown mass uncertainty increases the sample size, but also introduces an additional source of error, which explains the slightly lower success rates for dataset “all”. The “control” data subset was specifically created to remove any potential bias on the empirical models caused by using the same datasets for testing as was used for calibration. When being tested with the “control” subset for both matching criteria, Mastin remains among the models with largest S. We therefore infer that Mastin can be seen as the “optimal global model” without falling into the trap of circular argumentation. If the aim is to obtain MER predictions without the requirement or knowledge of additional eruptive conditions, applying the Mastin model is a reasonable strategy. This approach is currently applied, e.g., at the London VAAC (Beckett et al., 2020). Our findings are also complementary to (and in good agreement with) findings from a recent study which tested models with the IVESPA data and found empirical models to be more accurate than more complex wind-affected models, including the (in their case uncorrected) Degruyter model (Aubry et al., 2023). While the exact values for S depend on the test dataset used, the general trends of S remain unaffected by the choice of the dataset.

The lower panel of Figure 4 displays the gains in S if optimized model solutions other than the optimal global model is applied. For example, with the “best range” matching criterion, S is increased by 1.2%–2.9% if, instead of simply using Mastin, an optimal combination of models (such as given by Table 4) is globally applied. The “best fit” matching criterion reveals similar gains in S. It is worth noting that for this matching criterion the same weight factor solutions are found for all tested datasets (Table 4): the highest success rates are achieved by using a model combination of 90% Mastin and 10% Sparks. By adding 10% of Sparks to the Mastin model predictions, the success rates of the optimal global mix slightly increased by 0.7%–1.9% (see also Table 1).

A further increase in S can be obtained by using situation-specific model solutions. When applying the “best fit” matching criterion with FMERs, the success rates of optimal situation-specific models range from 16.7% (Wilson, Sparks, Mastin for MER < 10⁴ kg/s) to 65.0% (Sparks for buoyancy-dominated plumes). Mastin qualifies as the optimal situation-specific model for 14 of 21 tested eruptive conditions, Sparks for 6, Woodhouse for 4, and Gudmundsson and Wilson each for 2. If the “best range” matching criterion is used instead (see Table 3), both Mastin and Degruyter are optimal situation-specific models for 8 of 21 tested eruptive conditions, Woodhouse for 6, Wilson for 5, and Sparks for 2. The success rates range from 60% (for phreatic eruptions) to 100% (Degruyter for eruptions with more than 72 h duration). These findings demonstrate that when compared with each other, the tested plume models should not be regarded as generally “worse” or “better”, but rather as being more likely or less likely to provide an accurate MER prediction under given eruptive conditions.

The “best range” matching criterion is a rather lenient one, which is in principle already fulfilled by the slightest overlap between measured and modelled mass. The type of overlap can be further statistically analyzed by distinguishing between three cases, which describes where the modelled mass lies compared to the measured, see Figure 1 for illustration. These are i) cases where the modelled mass is situated in the lower sector of M _IVESPA ; ii) cases where modelled mass and M _IVESPA completely overlap; iii) cases where the modelled mass is situated in the upper sector of the M _IVESPA , and therefore is higher than the measured range. The pie diagrams of Figure 5 are an example of such a statistical analysis. The sum of presented success rates for the lower sectors (green), upper sectors (dark blue) and for complete overlap (bright blue) coincide with the success rates for the “best range” criterion for different plume types, as presented in Table 3. For wind-affected plumes, the “best range” success rates of Mastin and Degruyter coincide, but the closer examination of Figure 5 reveals that the matches by Degruyter are more frequently characterized by a complete overlap than for Mastin (49% vs. 38%). Based on this result, it could be inferred that from the two models, Degruyter might have a higher likelihood for successful MER prediction for wind-affected plumes. We note, however, that a large percentage of complete overlap is only meaningful for cases where the error bars of the predicted MER are of comparable size to ΔM _IVESPA . In cases where the model uncertainties are far larger, a high overlap is not necessarily an indicator for preciseness and should therefore not be over-interpretated.

When comparing the success rates for matches with the lower sector with those for matches with the upper sector, Mastin shows with 28% versus 10% a clear prevalence for the former. This indicates a tendency to underestimate the MER for wind-dominated plumes. The same bias is also visible for intermediate plumes (21% vs. 8%) but vanishes for buoyancy dominated plumes (15% vs. 15%). This dependency of prediction quality on windspeed is a phenomenon for simple empirical plume models that is well known from previous studies (Bursik, 2001; Mastin, 2014; Dürig et al., 2022) and was the main reason for the development of wind-affected plume models in the first place (e.g., Bursik, 2001; Degruyter and Bonadonna, 2012; Devenish, 2013; Woodhouse et al., 2013; Mastin, 2014; de’Michieli Vitturi et al., 2015; Folch et al., 2016). It is therefore an interesting side note that for intermediate plumes, Wilson, the simplest (and oldest) of all tested models, is the plume-type specific model with the highest success rates. The complexities of modelling wind-affected plumes are reflected by the fact that the success rates for wind-dominated plumes are significantly lower than for buoyancy-dominated plumes (see Table 3; Figure 5).

The Gudmundsson model was developed to model the 2010 Eyjafjallajökull eruption, which produced basaltic andesites (Gudmundsson et al., 2012). Surprisingly, for Gudmundsson the “best fit” success rates for basaltic andesitic eruptions are the lowest (16.2%), while it shows its highest success rates for dacites/rhyolites (30%). This might indicate that the model scale factor k _I is less dependent on composition than on other parameters, such as wind conditions and eruptive style.

The other two empirical models (Sparks and Mastin) are trained on datasets where basaltic eruptions are arguably underrepresented. However, we do not observe any significantly lower success rates for eruptive events of basaltic composition (Table 3). This finding speaks against a significant influence of a compositional bias in the underlying datasets.

The optimal style-specific models found were Woodhouse for magmatic and Degruyter for phreatomagmatic eruptions. Phreatomagmatic eruptions are characterized by a much larger water content in the plume, which significantly changes the heat capacity of the mixture (Sparks et al., 1997; Degruyter and Bonadonna, 2012). Unlike magmatic eruptions, for which the source of energy is provided by expanding gas, phreatomagmatic eruptions are driven by explosive molten fuel-coolant interaction (MFCI) (Zimanowski et al., 2015; White and Valentine, 2016). MFCI is a thermohydraulic mechanism that converts the heat of the magmatic melt into shock waves and mechanical stress, which drives the fragmentation of the melt and leads to the generation of ash particles (Dürig et al., 2012b; Dürig and Zimanowski, 2012; Dürig et al., 2020). The energy expended for both vaporization of the water and fragmentation of the melt results in a significant cooling of the fragmented melt (Spitznagel et al., 2013; Moitra et al., 2020). Consequently, the temperature of tephra in phreatomagmatic plumes is expected to be considerably lower than for ash columns resulting from magmatic explosions. Another important difference between phreatomagmatic and magmatic eruptions is the abundance of fine “interactive” ash particles in phreatomagmatic plumes that together with agglomeration effects, particle aggregation processes and the presence of fragmented host rock lead to differences in grain size distributions and the plume’s bulk density (Zimanowski et al., 2003; Dürig et al., 2012a; Dürig et al., 2020; Andrews et al., 2014; White and Valentine, 2016). All these factors influence the plume dynamics, which might not only explain why the models’ S are so different for both plume types, but also why the success rates for magmatic eruptions are generally higher (with a maximum “best range” S of 88.7%, see Table 3) than for the more varied phreatomagmatic eruptions (with a maximum S of 80.6%).

The results in Table 3 also show that a correct identification of the eruptive condition is crucial for an effective use of the situation-specific model strategies, because a misclassification would lead to significantly lower success rates. For example, let’s assume that a phreatomagmatic eruption is wrongly identified as magmatic. If the optimal style-specific model for magmatic eruptions (Woodhouse) is applied to phreatomagmatic eruptions, the success rates are not as expected (88.7%) but only 67.7%, which is lower than for any other of the tested models. Conversely, Degruyter has, with 77.4%, the lowest success rates of all models for magmatic eruptions. In both cases, the optimal global model (Mastin) yields higher success rates, which means that the situation-specific modelling strategy was not helpful, but in fact harmful for optimizing the forecast quality. This phenomenon is observed for all of the other grouping variables and datasets tested. We therefore conclude that the strategies of optimal situation-specific model (or optimal situation-specific mix) solutions should only be applied if the eruptive condition used for tailoring the modelling strategy is well known. Otherwise, it is highly recommended to fall back on the optimal global mix or optimal global model, which are both easier and more accurate than the situation-specific solutions.

The only situation-specific mix solutions that have shown to be more successful than the optimal eruptive style-specific model solution are the MER-specific solutions (see example in Table 5). Using the MER for the decision of choosing the weight factors is, however, not expected to always be a very practical procedure, given that MER is the unknown parameter which we aim to simulate with the models. The only viable way for the MER-specific mix strategy to succeed is by applying an iterative procedure. Based on these considerations, we designed the procedure displayed in the flow chart in Figure 6 to forecast the MER in a real-time plume monitoring scenario with highest accuracy. First it is checked if meteorological data (required for the wind-affected plume models Woodhouse and Degruyter) are available. If this is not the case, we suggest using the optimal global model for MER forecast (i.e., Mastin). If meteorological data are available, in a next step the MER category of the current eruption is classified based on a preliminary MER estimated with the optimal global mix (Table 4).

Examples for eruptions of unambiguous strengths are eruptions with a MER << 10⁴ kg/s, MER >>10⁷ kg/s or not significantly wind-affected eruptions with a mass eruption rate that is well limited within 10⁴ kg/s and 10⁶ kg/s. Based on the MER classification, the appropriate weight factors should be selected according to Table 5. Interestingly, for eruptions with MER <10⁴ kg/s, Mastin has shown to be the optimal model, despite the fact that such weak eruptions were underrepresented in the dataset underlying this model (Mastin et al., 2009; Mastin, 2014). When it is not possible to constrain the category of eruptive strength with sufficient certainty, it is better to choose another grouping variable, despite the slightly lower success rates.

The optimal situation-specific mix solutions grouped by plume type, eruptive style and duration all have identical success rates of 84.6% (Supplementary Table S4). From these parameters, duration is realistically the least practical one for real-time monitoring purposes, due to its low predictability. Since a plume-type classification has to be conducted anyway to apply the necessary plume height corrections for the wind-affected plume models, we suggest using optimal situation-specific mix solutions tailored for plume-types. Suitable model weight factor settings are presented in Table 5.

The plume heights used for this study are taken from the IVESPA database, which were acquired by a multitude of different methods, including visual observations and measurements with radar, lidar or satellites (Aubry et al., 2021). This bears a potential of errors, since it has been shown that the MER prediction quality is dependent on the method of plume height source (Aubry et al., 2023). Future studies could therefore focus on subsets of dataset “all” that only contain events with known plume heights from one measurement technique. A similar approach was followed, for example, in recent studies on Eyjafjallajökull 2010 (Dürig et al., 2022) and Etna (Mereu et al., 2023). As demonstrated for the 2010 Eyjafjallajökull eruption, the temporal resolution of the plume height average has a significant impact on the outcome of the modelled MER (Dürig et al., 2022). All plume height values from IVESPA are averaged over the total duration of the eruptions, a parameter which itself varies considerably. While a single average value for plume height can be sufficient to constrain a rather short-lived event of, for example, less than 6 h duration, it is less useful for quantifying the strengths of long-lived eruptions that feature changing source and atmospheric conditions. The fact that we found different results for different durations reflects this effect.

The weather data from IVESPA used in our study are based on ERA5, which is probably the best global product available, but due to effect of necessary interpolation from grid points of the global model to the location of the volcano, there can be significant errors in, e.g., windspeed. However, it is expected that this error is small in comparison to the larger uncertainties introduced from plume height and duration.

Seven out of the 130 eruptive events from dataset “all” were identified as “outlier events” (see Supplementary Table S5). These are events for which the simulated mass never meets the “best range” matching criterion with the range of the measured erupted mass (M _IVESPA ), regardless of the type of optimized model solutions. All identified events are also part of the “precise” data subset.

In the following, the outlier events are listed and potential aggravating factors and circumstances are discussed that might have complicated the predictability of these outlier events:

(i) Redoubt Alaska 1989 December 19: this was a very short-lived phreatomagmatic eruption, which lasted only 9 min (Miller and Chouet, 1994). All tested models tested assume a steady plume, which is fed by a stream of ash and gas with a constant MER (Morton et al., 1956; Sparks et al., 1997). This assumption is not suitable for eruptions of very short duration (for example, when the ascent time of the ash column is in the same order of magnitude as the total duration of the eruption), or for pulsatile behavior with sufficiently long repose times (Dürig et al., 2015b). Therefore, all models underestimate the MER.

The six models examined in this study have the advantage that they use with plume height h as main input parameter, which is relatively easy to obtain, despite the important limitation that the value can be dependent on the method of measurement itself (Dürig et al., 2022; Aubry et al., 2023). Relatively small errors in h, however, result in considerable uncertainties in MER, often covering one or more orders of magnitude (see Figure 3), due to the exponential relationship between MER and h. Applying the discussed strategies of mixing different plume height-based models cannot compensate for this effect, because it is intrinsic to all models tested.

Approaches for real-time MER forecast that do not use plume-height exist, but they require i) more sophisticated instruments and ii) knowledge of additional parameters that are often difficult to acquire. For example, MER prediction by infra-sound requires detailed information on the vent geometry (Johnson and Ripepe, 2011; Ripepe et al., 2013), without which the uncertainties are considerable (Dürig et al., 2015b). Other methods, like electric field measurements are still in development and need further calibration and testing under different eruptive conditions (Büttner et al., 2000; Calvari et al., 2012).

If instruments and necessary input data are available, the most promising strategy to better constrain the MER is to combine several independent prediction methods (e.g., estimates from infra-sound with the presented plume-height based 0D model forecasts).