Uncertainty in eruption parameters
Quantifying uncertainty in the eruption parameters is challenging, particularly in the first hours of an eruption when reliable quantitative observations are often not available (Aubry et al. 2021). Three key eruption source parameters (ESPs) for modelling volcanic ashfall transport are the eruption plume height (H), mass eruption rate (MER), and eruption duration (D). Plume height is not trivial to reliably measure, can be variable on short timescales, and estimates may have large uncertainties. The mass eruption rate is not directly measurable but usually derived from empirical correlations based on the observable plume height, in relationships that have confidence intervals spanning orders of magnitude (Aubry et al. 2021; Mastin et al. 2009). Finally, eruption duration is related to mass eruption rate and total volume erupted, and difficult to forecast while an eruption is ongoing. Existing data consists of a mix of directly observed eruptions and derived estimates based on total ash erupted. What constitutes the ‘end’ of an eruption can also be difficult to ascertain.
Volcano database and ESP PDFs
We undertook a comprehensive desktop study to develop, for the first time, PDFs for eruption plume height, the mass eruption rate, and eruption duration for the 12 volcanic centres monitored by GeoNet (2022). The intention is to have a priori distributions ready to use if there is confirmation of an eruption. To achieve this, we compiled a large dataset, combining datasets from Deligne (2021) and Aubry et al. (2021). We briefly describe both below.
Deligne (2021) compiled measured, estimated, and calculated estimates for volcanic eruption mass eruption rates, column heights, and durations from 213 events, a third of which are from New Zealand volcanoes. These include both prehistoric (estimates based on geologic record) and historic eruptions. Eruptions were categorized according to magma type (mafic, intermediate, silicic) and eruption style (steam-driven, magmatic small to moderate, magmatic large). Every entry was coded according to data type (e.g. instrumentally measured, observed, derived from geologic study) and assumptions made in determining values. There was no attempt to evaluate the quality of the studies providing individual parameter values, nor an exclusion of derived values (e.g. mass eruption rate is often calculated based on eruption column height, duration, and/or other parameters, but is rarely independently measured). For this study, values compiled in Deligne (2021) derived from empirical or numerical models were removed.
Aubry et al. (2021) independently compiled estimates of total erupted mass of fall deposits, duration, eruption column height, and atmospheric conditions from 134 eruptions from around the world from 1902–2016. Eruptions had to have independently measured/available values for all four parameters. Each eruption was independently evaluated by at least two experts, and uncertainty associated with each parameter was systematically evaluated, with estimates requiring a high degree of interpretation flagged.
For this study, we attributed magma type and eruption style categories using the criteria documented in Deligne (2021) to eruptions that solely appear in Aubry et al. (2021). For 22 eruptions, there was no indication of eruption style apart from being classified as VEI 4 or smaller by the Smithsonian Global Volcanism Program: we assigned these to the ‘Magmatic small to moderate’ eruption style category. Table 1 summarizes the collected data and Fig. 3 shows the PDFs of ESPs by magma type.
Table 1
Summary of eruptions data compiled to estimate ESP distributions.
|
Steam-driven
|
Magmatic small to moderate
|
Magmatic large
|
Total
|
Mafic
|
18 (7%)
|
78 (30%)
|
4 (2%)
|
100 (38%)
|
Intermediate
|
27 (10%)
|
85 (33%)
|
22 (8%)
|
134 (51%)
|
Silicic
|
1 (< 1%)
|
15 (6%)
|
10 (4%)
|
26 (10%)
|
Total
|
46 (17%)
|
178 (69%)
|
36 (14%)
|
260 (100%)
|
Uncertainty in weather
Uncertainty in weather is nowadays quantified with NWP ensemble datasets (Cheung 2001). These datasets seek to capture forecast error growth by running many model realisations, with perturbed initial conditions that reflect uncertainty in the current atmospheric state, and perturbed model parameters that reflect uncertainty in physical processes. The non-linear nature of the equations governing atmospheric motion means that even with perfectly known initial conditions, the tiniest perturbations in atmospheric variables in the most accurate of NWP models may result in different outcomes after a finite time (Lorenz 1969; Zhang et al. 2019). The output of these ensembles is probabilistic, which more accurately reflects what is knowable about the future atmospheric state than the deterministic picture provided by a single model.
NWP ensembles
The NWP ensemble datasets used in this study are from the Global Ensemble Forecast System version 12 (GEFSv12) running at the National Oceanic and Atmospheric Administration (NOAA). The operational dataset is updated 4 times per day (cycles starting at 0, 6, 12 and 18:00 UTC) with 31 members (30 perturbed and 1 unperturbed/control) with approximately 25 km horizontal resolution and 64 vertical hybrid levels for atmospheric components, and out to 16 days of forecast at each cycle, except for 35 days at 0000 UTC (Stajner et al. 2020). The model output is public and available via NOMADS Server (2022) at a 0.5° resolution grid, 3 hours temporal resolution, with surface and 12 pressure level fields up to 10 hPa, for the past 30 days.
We also made use of the historical dataset spanning from 2000–2019 consisting of GEFSv12 reforecasts (GENS-3), available via Open Data on AWS (2022) but with only 21 members (20 perturbed and 1 unperturbed/control) available at a lower resolution (1° resolution grid and 6 hours temporal resolution) and with less variables and less vertical levels (surface and 10 pressure level fields up to 10 hPa).
Combining uncertainties and sampling with Latin Hypercube technique
Multiple approaches can be taken to quantify uncertainty both from the source parameters as well as from the weather conditions. In our rapid response context, the biggest challenge of probabilistic forecasting is to optimize resources while ensuring the spread and natural variability of our multidimensional parameter space are well represented.
Monte Carlo sampling methods have been previously used in ash probabilistic forecasting (Hurst and Smith 2004; Magill et al. 2006) but they are only computationally feasible with simplifying assumptions that limit the sample size. For instance, only using a few eruption sizes, and/or only using wind profiles at the vent location. While these assumptions may be reasonable for long-term hazard studies, they would likely result in inaccurate forecasts for a single eruption and near real-time forecasting.
LHS is a stratified sampling approach that draws representative samples using a smaller sample size than Monte Carlo sampling, ensuring every section of the parameter space is sampled once (Mckay et al. 2000). Sigg et al. (2018) and Prata et al. (2019) both used LHS techniques as an alternative to Monte Carlo in dispersion applications. Sigg et al. (2018) considered a known amount of pollutant and the other LHS parameters were associated with microscale boundary layer meteorology important for their specific small-scale application. Prata et al. (2019) considered large scale transport of airborne ash and considered various eruption parameters as independent parameters in their LHS.
For our purposes, we considered the GEFSv12 ensemble dataset to be reliable in a statistical sense (i.e., over time, ensemble member counts of some categorical weather event are consistent with their actual probability). The atmospheric evolution can be sampled by selecting an ensemble member, and as we considered the eruption and atmospheric state to be independent, the ensemble member label can be considered an independent parameter with uniform distribution (Stein et al. 2015b).
LHS with Dependence
The nature of the sampling problem for ESPs is different because they are not independent of each other – an underlying assumption of standard LHS technique. To introduce correlation between ESPs, we assume that they approximately follow a multivariate Gaussian distribution (Packham and Schmidt 2008). We then first draw samples for each parameter from a standard Gaussian distribution \(x\sim N\left(\text{0,1}\right)\) and then match the mean \(\mu\) and covariance matrix \(\varSigma\) of the prior eruption parameter distribution. Because we have three ESPs, \(x\) and \(\mu\) are 3x1 vectors and \(\varSigma\) is a 3x3 matrix. A single LHS draw \(y\) can then be expressed as:
\(y = Lx + \mu with \varSigma =L{L}^{T}\)Equation 1
where \(L\) is the Cholesky factorisation of \(\varSigma\) (Murphy 2022). It follows that \(y\) has the same mean as the original distribution. To see that \(y\) also has the same covariance matrix we can calculate the covariance of \(y\):
\(Cov\left(y\right) = LCov\left(x\right) {L}^{T}= L I {L}^{T} = \varSigma\) Equation 2
with \(I\) being the Identity matrix.
Figure 4 compares Monte Carlo and LHS techniques, demonstrating that LHS allows us to generate a representative sample of a given distribution with considerably fewer samples than Monte-Carlo sampling requires.
LHS modelling runs
As previously mentioned, we considered each perturbed atmospheric state to be independent of each other and of the ESPs. To optimize computational resources and maximize the weather variability information available, we constrained our ESP sample size to the same as the number of NWP ensembles, (30 for GEFSv12 forecasts and 20 for GENS-3 reforecasts). Figure 5 shows that even a sample number of 20 suffices to have a good approximation to the PDFs derived from historic eruptions (Fig. 3).
We then used these “eruption scenarios”, i.e., each point with a sampled MER, H, D and NWP ensemble member, to force independent ashfall dispersion model runs, constraining H to a maximum of 20 km above mean sea level to fit within the weather forcing data, and D to a maximum of 24 hours which is the forecast length. Table 2 summarizes the run parameters used as input to the dispersion model.
Table 2
HYSPLIT dispersion run parameters selected for each sample taken with LHS
Parameter
|
Sampling range
|
Plume height (H)
|
0–20 km above vent
|
Mass eruption rate (MER)
|
Unconstrained
|
Eruption duration (D)
|
0–24 h
|
Meteorological fields (Operational case study)
|
GEFSv12 members 1–30, 0.5° resolution, 3 hourly
|
Meteorological fields (Historical case study)
|
GENS-3 members 1–20, 1° resolution, 6 hourly
|
HYSPLIT configuration
Ashfall dispersion was computed with the HYSPLIT model version v5.0.1 (April 2020), a state-of-the art hybrid Eulerian and Lagrangian dispersion model that is extensively used for volcanic ash transport and deposition by the international community (Rolph et al. 2017; Stein et al. 2015a).
The real and hypothetical eruptions studied here were located at Tongariro volcano (39.130 °S, 175.642 °E) with the main vent at 1978 m above sea level. It typically erupts andesite lavas (Leonard et al. 2021), which we classified as “intermediate magma” in our eruption database (Table 1 and Fig. 3). The configuration options used were implemented in Hurst and Davis (2017) and a brief description follows. The particle size distribution used for this volcano was andesite with a density of 1300 kg/m3, divided into 37 size bins (see Fig. 6). The total mass erupted was distributed along the eruption column according to a Suzuki distribution with constant 4 and discretized in 10 levels. This distribution adds more mass to the upper levels of the eruption plume, as is typical for most eruptions (Suzuki 1983). Figure 7 shows the vertical plume mass profile for a sample of ESPs. Model results were integrated to a grid of 1 km horizontal resolution spanning 14 degrees latitude and longitude, from 48.5 °S to 33.5 °S and 165.5 °E to 180.5 °E.