Volcanic Hazard Assessment for an Eruption Hiatus, or Post-eruption Unrest Context: Modeling Continued Dome Collapse Hazards for Soufrière Hills Volcano

Effective volcanic hazard management in regions where populations live in close proximity to persistent volcanic activity involves understanding the dynamic nature of hazards, and associated risk. Emphasis until now has been placed on identification and forecasting of the escalation phase of activity, in order to provide adequate warning of what might be to come. However, understanding eruption hiatus and post-eruption unrest hazards, or how to quantify residual hazard after the end of an eruption, is also important and often key to timely post-eruption recovery. Unfortunately, in many cases when the level of activity lessens, the hazards, although reduced, do not necessarily cease altogether. This is due to both the imprecise nature of determination of the “end” of an eruptive phase as well as to the possibility that post-eruption hazardous processes may continue to occur. An example of the latter is continued dome collapse hazard from lava domes which have ceased to grow, or sector collapse of parts of volcanic edifices, including lava dome complexes. We present a new probabilistic model for forecasting pyroclastic density currents (PDCs) from lava dome collapse that takes into account the heavy-tailed distribution of the lengths of eruptive phases, the periods of quiescence, and the forecast window of interest. In the hazard analysis, we also consider probabilistic scenario models describing the flow’s volume and initial direction. Further, with the use of statistical emulators, we combine these models with physics-based simulations of PDCs at Soufrière Hills Volcano to produce a series of probabilistic hazard maps for flow inundation over 5, 10, and 20 year periods. The development and application of this assessment approach is the first of its kind for the quantification of periods of diminished volcanic activity. As such, it offers evidence-based guidance for dome collapse hazards that can be used to inform decision-making around provisions of access and reoccupation in areas around volcanoes that are becoming less active over time.

Introduction and Motivation

The ability for communities around volcanoes to recover from volcanic unrest hinge on well-judged and timely decisions about land use, as well as access to some areas around volcanoes for reasons of work, attending personal property, and/or potential reoccupation. However, identifying and defining, as such, “the end” of an eruption is notoriously difficult (Siebert et al., 2011; Ogburn et al., 2015). Furthermore, hazard assessments that account for diminishing activity, cessation in eruptive activity and/or the establishment of new post-eruption background activity levels are hard to develop given that the processes involved are poorly understood, and are necessarily characterized by high uncertainties. To our knowledge no such hazard assessments have been carried out for any eruptions so far.

This work is motivated by the need for reassessing the hazards from the Soufrière Hills Volcano (SHV), Montserrat, given the now extended state of quiescence of the volcano since 2010.

Recent Eruptive History at Soufrière Hills Volcano

The eruption has been characterized throughout by alternating periods of active lava dome extrusion (months–years), punctuated by sequences of Vulcanian explosions (up to VEI 3) (Wadge et al., 2014b). The last phase of dome growth ended on February 11, 2010 after an explosion removed about 20% of the $248\cdot {10}^6{\text{m}}^3$ dome (Stinton et al., 2014; Cole et al., 2015). Sporadic, small dome collapse PDCs and rockfalls generated by progressive collapse of the remaining dome and crater area occurred until 2013, and have subsequently phased out almost entirely. In the period 2016–2020 there were typically less than 10 rockfalls/year. The volcanic activity therefore, as observed at the surface, has been very low since 2010. This quiescence now represents the longest pause since the 1995 eruption began. Deep but low-level residual seismicity, low level degassing and slow inflation (Figure 1) indicates that the volcanic system is still in a state of unrest. The causal processes for this continued unrest are hard to elucidate. Although for a while it was interpreted as representing continued pressurization of the system, it is now considered possible that current unrest is at least in part the manifestation of ongoing visco-elastic response to the initial emplacement of magma in the crust (Scientific Advisory Committee of the Montserrat Volcano Observatory, 2019).

FIGURE 1

FIGURE 1. Seismic, Global Positioning System (GPS), and SO₂ monitoring data for the period January 1, 1995–November 7, 2019. Extrusive phases and pauses are in shown red and green, respectively. Top: Number of seismic events detected by the seismic system. Middle: Radial ground displacement of continuous GPS (cGPS) stations MVO1 (red) and GERD (blue) smoothed with 7-days running mean filter, Black: GPS Height of HARR. Bottom: Measured daily SO₂ flux, filtered with 7-days running median filter. Green: Correlation Spectrometer (COSPEC), Blue: Differential Optical Absorption Spectroscopy (DOAS), White: Traverse data, Red: new DOAS network (Adapted from MVO activity report; MVO OFR 18-02-draft).

In this work we try to address the issue of hazard assessment, specifically tuned for the current situation within this long hiatus period or potentially (still to be shown) post-eruption scenario. We focus specifically on hazards associated with lava dome collapse as these are recognized to be a key high-impact hazard that could still occur with little or no warning. We recognize that in the case of a re-start of eruptive activity, other precursory activity might be the first to manifest—such as ash venting or explosive activity. This work does not address the question of forecasting which events are most likely to occur, either during a continued hiatus period or after a re-start. Instead we address the question of how likely are dome collapse hazards to occur and, if so, where and when.

The broader rationale for this hazard assessment is provided by the following risk management questions which are currently being considered on island: To what extent is it safe to work in, visit, or re-inhabit the flanks of the Soufrière Hills Volcano and surrounding area in the south of Montserrat? And, in balancing the long term continued exposure (albeit now to very low levels of activity) with short-term development gains; for example, by supporting the extractive (sand mining) and tourist industries which might be key to the island’s economic recovery: To what extent are re-investments in this area possible and prudent? Ultimately these questions are for civil protection authorities, the Government of Montserrat, and the Montserrat Volcano Observatory to consider but they also weigh on the Scientific Advisory Committee (SAC) of the UK Government’s Foreign and Commonwealth Office Montserrat. It is under that auspice that this work is undertaken.

It is noted that the conditions and mechanisms for PDC generation are different now from those during the main eruptive period. During periods of active lava dome extrusion, the first-order control on PDC generation is by lava dome extrusion and associated dome instability and collapse (Calder et al., 2002). In the absence of dome growth during pauses, a long hiatus (such as the current one), or post-eruption unrest, other collapse mechanisms can come into play. Lava domes can collapse years, decades, or millennia after volcanic activity is over, triggered by bouts of intense rainfall, new seismic activity, and/or weakening by hydrothermal activity and structural instabilities (Ball et al., 2013; Harnett et al., 2019). Post-eruption collapses of lava domes, or more extensive parts of volcanic edifices are known as sector collapses and they generate debris avalanches. There is substantial evidence across the Lesser Antilles volcanic arc (e.g., Boudon et al., 2007; Sampler et al., 2008) and elsewhere, for the generation of modest to extensive debris avalanches sourced from collapses of lava dome complexes. At the volcano-scale these can be considered “extreme events,” and as such, for any given volcano, there is little data to constrain the frequency of occurrence and anticipated time-scale post-eruption until one of these post-eruption events might occur. One can consider the different types of collapse events on a spectrum where debris avalanches represent the larger volume, lower recurrence rate collapses, and where typical dome-collapse PDCs during eruptions are the smaller volume, higher frequency events. Our consideration of “collapses” in this work, implicitly considers this spectrum of plausible scenarios.

We use the catalog of previous collapse events at SHV, and how that has changed over time, to inform our analysis. Over the course of the SHV eruption more than 900 PDCs with runout over 1 km were produced. Most PDCs were the result of dome collapse events, but column-collapse PDCs were also common. Dome-collapse PDCs show a first-order correlation with periods of dome growth as seen in Figure 2, and there is a known association between high dome-extrusion rates and frequent dome-collapses (Calder et al., 2002; Wadge et al., 2010). Calder et al. (2002) also noted that many larger-volume collapses were associated with periods of elevated dome extrusion, intense seismicity, and/or deformation cycles. Lava extrusion mechanically destabilizes growing domes by increasing internal shear stress, increasing loading on support structures, and by over-steepening dome structures (Calder et al., 2002; Pallister et al., 2013). High extrusion rates are also frequently associated with gas pressurization (Voight and Elsworth, 2000) and explosive activity that can cause dome-collapse events or produce column-collapse PDCs. The dominant mechanisms controlling the size and frequency of dome-collapse PDCs during dome growth thus include the volume of the dome, the extrusion rate, and associated explosive activity that can cause large, violent dome-collapse events. The collapses from this period are dominantly gravitationally induced dome collapse events, but a number of other causes contribute (Calder et al., 2002). The slope of the cumulative count of PDCs in Figure 2 episodically steepens and shallows, reflecting intervals of higher and lower rates of PDC production associated with dome growth phases. Wolpert (2018) estimate the annual rate to range between a high level of roughly 260 events/year and a low level of about 30 events/year, based on evidence from the period 1995–2012. Including the quiescent period from 2013 to 2020 (in which no PDCs were recorded) reduces the estimated low-frequency PDC expected annual rate to roughly 17 event/year. Even this reduced rate, however, is inconsistent with the current behavior—the probability of a PDC gap of seven or more years is almost zero for a Poisson model with rate as high as 17/year. This indicates that a different model is required to explain the current hiatus in PDC activity. That is, with respect to dome collapse events, the activity levels during this extended pause is significantly quieter than that during the previous pauses. A Poisson model will still fit the data well, but the annual rate must be much lower (below 0.20/year or so) for a seven year PDC gap to be plausible. We now take this into account in our forward modeling for future collapse hazards.

FIGURE 2

FIGURE 2. Cumulative number of recorded, dome-collapse PDCs, beginning July 18, 1995 and ending with the last recorded PDC on March 28, 2013. Data are from FlowDat (Ogburn and Calder, 2012; Ogburn and Calder, 2017).

The work developed here presents a defensible, evidence-based approach to PDC forecasting, in a post-eruption context, which accounts for attendant uncertainties and assesses specifically the threat posed by inundation of infrequent post-eruption-unrest flows that could affect southern Montserrat. To do so, we develop a new statistical model for forecasting the probability of post-eruption-unrest flows occurring over the next s years, for a range of possible forecast periods. The model is based on considering a very low-frequency background level of activity that is balanced with the chance that activity has already completely ceased or will cease at some point in the future (Wolpert et al., 2016). More precisely, PDC events are modeled using a homogeneous Poisson process with constant expected rate λ events per year up to an uncertain time T at which the eruption ends. The absence of observed PDC events prior to the present time t suggests that either the eruption has already ended, i.e., that $T<t$, or else that $T\ge t$ but the rate λ is very small. The distribution for the end-of-eruption time T is based on observations of hundreds of dome-collapse volcano histories (Wolpert et al., 2016). This is the key, and novel contribution in this work, which we believe to be the first of its kind. The model is then combined with a simulation-based strategy to assess the impact on a specific area of interest from flows corresponding to different scenarios (e.g., flow volume, initial direction) (Spiller et al., 2014; Bayarri et al., 2015). In its own right, the new statistical method for forecasting the probability of post-eruption flows is a potentially powerful approach that could readily be applied to other volcanoes. Likewise, the combined statistical-simulation modeling approach could be undertaken at other volcanoes and used as a tool to inform decisions about land use planning for communities situated near volcanoes prone to PDCs, during or after eruptive periods or during prolonged low-level unrest.

Post-Eruptive Pyroclastic Density Currents and Debris Avalanches

Currently, neither lava dome growth nor explosive activity has occurred at SHV since February 11, 2010, and no PDCs have been produced since 2013. However, the large lava dome remaining at the summit means that dome-collapse PDCs can still be generated. While frequency-volume statistics exist for dome-collapse PDCs during eruptions, only very sparse information exists for PDC-generating collapses of inactive lava domes (e.g., Ball et al., 2015; Harnett et al., 2019). However, collapses during low to no eruptive activity do occur, both from the collapse of still hot lava domes to produce PDCs and from collapse of unstable, older domes and volcanic flanks to produce colder debris avalanches and rockfalls.

Statistical analysis of rockfall activity at SHV led Calder et al. (2002) to suggest that, during periods of low dome extrusion or pauses in dome growth, external forces, such as extreme rainfall, seismic activity, and talus apron erosion more strongly control failure of lava domes. For example, following heavy rainfall, the July 3, 1998 dome collapse removed 20% of a $110\cdot {10}^6{\text{m}}^3$ meta-stable dome that developed between November 1995 and February 1998, but which was no longer actively extruding (Norton et al., 2002; Elsworth et al., 2004). Shallow or intense seismic activity at SHV was shown to destabilize unconsolidated talus aprons and the dome itself via direct shaking of the edifice and to contribute to collapses (Calder et al., 2002).

Dome growth also contributes to the instability of adjacent hydrothermally altered volcanic edifices, which may remain unstable and prone to failure in the form of debris avalanches for decades after dome growth has ceased. Hydrothermal alteration and partial flank and talus collapse triggered the December 26, 1997 debris avalanche at SHV. Though this collapse did occur during active dome-growth, evidence suggests (Ball et al., 2013; Ball et al., 2015; Harnett et al., 2019) that dome and flank collapses caused by instability due to hydrothermal alteration pose a significant, and difficult to anticipate hazard. For example, debris avalanches of lava domes have occurred at La Soufrière, Guadeloupe (Salaün et al., 2011); Cerro Chascon-Runtu Jarita Complex, Bolivia (Watts et al., 1999); and Las Derrumbadas Dome Complex, Mexico (Capra et al., 2002).

Several mechanisms are known to generate post-eruptive PDCs and/or debris avalanches:

Rainfall-Triggered Collapses of Lava Domes

Especially heavy rainfall occurred before and during the SHV collapse of July 3, 1998 (Norton et al., 2002), and this has been recognized elsewhere as an important collapse-triggering mechanism (Mastin, 1994; Yamasato et al., 1998; Ratdomopurbo and Poupinet, 2000; Voight, 2000; Matthews et al., 2002; Barclay et al., 2006). Although rainfall-triggered collapses may also occur during active growth periods, they are likely to represent a more important forcing mechanism during the months to years following cessation of extrusion. Barclay et al. (2006) examined rainfall and PDC data between 1998 and 2003 and found that, given a rainfall event of more than 20 mm, the probability of producing a multiple dome-collapse PDC episode increased by up to 9.2%. In fact, several of the largest dome-collapses at SHV were associated with extreme rainfall events (Matthews et al., 2002; Barclay et al., 2006).

Sector Collapse Associated With Hydrothermally Altered Systems

Lava domes and lava dome complex volcanoes, with mature hydrothermal systems are frequently associated with sector collapses (Vallance et al., 1995; Capra et al., 2002; Clavero et al., 2004). Lava dome and/or flank failures are often attributed to hydrothermally altered and weakened edifices [e.g., Mount Rainier in the United States (Reid et al., 2001); Casita volcano in Nicaragua (Opfergelt et al., 2006); and Santa María volcano (Ball et al., 2013)].

Two pertinent examples are the 1998 debris flow at Casita in Nicaragua (Sheridan et al., 1998; Voight et al., 2002; Opfergelt et al., 2006) and the 1997 debris avalanche at Soufrière Hills, Montserrat (Sparks et al., 2002; Voight et al., 2002). At Casita, an ∼8 ka dacite lava dome complex, a $1.6\cdot {10}^6{\text{m}}^3$ collapse generating a debris flow on October 30, 1998 was triggered by intense rainfall associated with Hurricane Mitch, resulting in more than 2,500 fatalities (Sheridan et al., 1998; Kerle, 2002). It was found that intense hydrothermal activity had altered the rocks to smectite clays (Opfergelt et al., 2006). For the December 26, 1997 debris avalanche at Soufrière Hills Volcano, Montserrat, hydrothermal alteration of the retaining crater wall was implicated as a major contributor to the destabilization and subsequent collapse and depressurization of the lava dome (Sparks et al., 2002; Voight et al., 2002). In both these cases low-strength, low-permeability, fines-rich alteration materials were thought to have helped to concentrate water along and lubricate structural discontinuities, thereby reducing the shear strength of the rocks and ultimately leading to destabilization of the edifice (Voight et al., 2002; Opfergelt et al., 2006).

Volcanic Island Flank Collapses

Large scale flank collapses are recurring processes in the Lesser Antilles volcanic arc (Deplus et al., 2001; Le Friant et al., 2004; Boudon et al., 2007; Sampler et al., 2008). Horseshoe-shaped structures within Dominica, Saint Lucia and Martinique consistently open toward the Grenada basin. The recurrence rate of flank collapse events in the Lesser Antilles may be linked to tectonic activity, and/or magma production rate, with intervals as high as ${10}^4$ years to several ${10}^5$ years (Sampler et al., 2008). There have been 15 flank collapses in the last 12,000 years in the Lesser Antilles (Boudon et al., 2007). In the northern arc, the potential for future flank collapse is high. Boudon et al. (Boudon et al., 2007) stress that edifice collapse can trigger laterally directed explosions and are often associated with tsunamis (Boudon et al., 2007). The same authors state that collapse of La Soufrière, Guadeloupe, is a “highly plausible geological scenario, whose risks must be taken into account.”

Downward Propagation of Cooling Fractures

At Mount St. Helens, this appeared to be the most likely mechanism that caused dome collapse-explosions in 1989–1991 that were associated with rainfall (Mastin, 1994).

Few forecasting models exist for eruption durations (Sparks and Aspinall, 2004; Wolpert et al., 2016); this in part is difficult because the definition of “end” for the catologue of long-duration eruptions is not very precise. It could refer to the end of any activity, such as dome collapse PDCs or post-eruption debris avalanches [which we can broadly consider as a spectrum from small to large mass flow events which can defensibly be modeled by flow simulators like TITAN2D (Patra et al., 2005)]. Likewise, as we define level of activity we do not distinguish among different types of collapse events, but merely focus on their average frequency or recurrence rate. As we will see, precise knowledge of the activity rate is not needed, as for any forecast window we can find the most conservative frequency (i.e., the one leading to the highest forecast probability) to find an upper bound on the probability of a collapse event whatever the actual frequency might be.

Methodology

Probabilistic hazard forecasting and analysis as we present here hinge on probabilistic modeling of aleatory scenarios, and combining those scenario models with physical models. A road block for that combination is the computational burden of exercising physical models for Monte Carlo calculations. As such, we advocate utilizing statistical emulators of computationally expensive simulators to overcome that burden, an approach that we have developed since 2009. In the section titled Forecasting and Planning With Dynamic Probabilistic Hazard Maps, we present a summary of our previously published statistical emulator methodology for efficient hazard calculations and walk through the process we advocate for creating probabilistic hazard maps. [Detailed descriptions of emulator-based probabilistic hazard methodology are available in Spiller et al. (2014) and Bayarri et al. (2015)]. In the section titled Forecasting With a Chance of Ending, we expand upon the previous approaches for the specific application to the current activity level, a protracted hiatus in eruptive activity, or a potentially post-eruption scenario. To do so, we present a new probabilistic model that balances a low-level of potential volcanic activity with the possibility that the eruptive phase has actually ended. In the section titled Quantifying Uncertainty in Hazard Forecasts, we describe how to quantify epistemic uncertainties associated with probabilistic (aleatoric) scenario descriptions.

As before (Spiller et al., 2014; Bayarri et al., 2015), we explore potential inundation footprints by future PDCs as modeled by TITAN2D. The TITAN2D computational environment simulates large, dry mass flows over terrain. TITAN2D solves depth-averaged mass and momentum conservation equations (Patra et al., 2005). A digital elevation model (DEM) provides a local coordinate system, yielding normal, down-slope, and cross-slope directions. A fixed mass of material is released with sampled dissipation parameters (internal and basal friction angles). Simulations were performed on a 1 m 2010 LiDAR-based digital elevation model with a vertical resolution of 15 cm m⁻² (Cole et al., 2010; Scientific Advisory Committee of the Montserrat Volcano Observatory, 2011) obtained from the MVO. The original DEM was smoothed to 10 m grid spacing to eliminate spurious features at high spatial frequency which can degrade the flow simulations.

Forecasting with a Chance of Ending

Sparks and Aspinall (2004) fit a generalized Pareto probability distribution to data consisting of the recorded durations of a subset of 137 dome-building eruptions taken from the Smithsonian Institution Global Volcanism Program Volcanoes of the World database (Simkin and Siebert, 1994) (note that all the distributions used in this work are described in Appendix 1). In this model the probability that the duration τ of an eruption exceeds any time $t>0$ is given by

$\mathrm{P}\left[\tau >t\right]={\left(1+t/\beta \right)}^{-\alpha }(1)$

for some constants $\alpha >0$ and $\beta >0$, to be estimated from the data (we used maximum likelihood methods to estimate $\alpha \approx 0.58$ and $\beta \approx 0.60$ years). Wolpert et al. (2016) applied an extension of that model to a more recent collection of 177 dome-building eruptions taken from the DomeHaz database (Ogburn et al., 2012; Ogburn et al., 2015); in this extension, the parameters $\alpha ,\beta$ (and hence the survival distribution) may depend on the magma composition. (Note, data used to fit this model and extension, including estimation of α and β, is available within (Wolpert et al., 2016).) This led in Wolpert et al. (2016) to a median forecast of 35–47 years of additional activity at SHV, under the explicit assumption of stationarity—i.e., assuming that by 2016 the volcano remained in the same active state it had begun in 1995. The present work is intended to remove that assumption, and make forecasts of future activity that reflect the likely possibility that the eruption has already subsided, or will soon do so.

Wolpert (2018) constructed and fit a dynamic model for volcanic activity in which the activity level λ (measured in expected PDCs per year) and preferred (i.e., most likely) direction could change over time. The best fit activity level function alternated between two values: a high activity of ${\lambda }_{\text{ hi}}=261.3$ events/year, and a low level of ${\lambda }_{\text{ lo}}=30.9$ events/year (see slope of Figure 2).

During the seven-year period from March 2013 through March 2020 there have been no PDC events large enough to be reflected in these models (the threshold was a volume of $0.15\cdot {10}^6{\text{m}}^3$). The probability that a volcano remaining active at a rate as high as $\lambda \ge {\lambda }_{\text{ lo}}=30.9$ events/year would exhibit no PDC in ten years is negligible (below ${10}^{-100}$), so it is statistically untenable that SHV has remained at the same activity level up to the present.

Two possibilities remain: that SHV remains active, but at a much lower rate λ (perhaps once in a decade, or once in a century); or that PDC generation from the 1995 eruption has ended. In this work we explore variations on those two possibilities.

Results

Three Possibilities

probability of no PDC occurring over the time frame including t years since the last recorded PDC and s years of forecast can be broken down into three cases:

(1) Eruptive activity may have ceased sometime between the last recorded PDC and time t, and hence we will see no PDCs through the forecast window.

(2) Eruptive activity may cease at some point during the forecast period, and until that time we get “lucky” that no PDC happens to occur.

(3) Eruptive activity may not cease until after the forecast window, and we get “lucky” that no PDC happens to occur during the forecast window.

We explore several cases for different activity levels, λ, and visualize the probability of PDCs occurring (or not) through a series of figures.

The thick solid lines in Figures 3A–C (see Eq. 3 for their derivation) indicate the forecast probability of at least one PDC at SHV in the next s years, as a function of s, for $0\le s\le 20$ years, assuming that the volcano remains active but at a much lower rate than in the decade 1995–2005: a rate of once annually in Figure 3A, once a decade in Figure 3B, and once a century in Figure 3C. The thin dashed and dotted lines (see figure legend) show the probability of no PDC in that same period (dashed black line), broken down into three possible cases: that the 1995 SHV eruption has already subsided (red dotted line); that it has not yet subsided, but will before s years have passed (blue dash-dot line); or that it has not yet subsided and will not during the next s years (green dashed line). Mathematical expressions for these curves are given in Eq. 3 below. Evidently the forecast PDC probability is lower for both a rate as high as $\lambda \approx 1.0\gg 0.10$, because with such a high annual rate λ it is overwhelmingly likely that a PDC would have occurred during these $t=7$ years, unless the eruption had paused; and also for $\lambda \approx 0.01\ll 0.10$, in which case the rate is so low that any PDC would be unlikely, even if the eruption is continuing. Figure 4 shows that the highest possible forecast probability of a PDC in the next $s=5$ years would be $P=0.0947$, achieved at an annual rate of $\lambda =0.085$, about one event per decade. Thus even without any knowledge of the rate λ, the forecast probability of a PDC at SHV in the next $s=5$ years is bounded above by about 9.5%. Similar results are available for any forecast period s.

FIGURE 3

FIGURE 3. Solid black curve shows forecast probability (nearly zero) of at least one PDC in s years, for $0\le s\le 20$, starting $t=7$ years after the most recent eruption. Assumed activity rates are (A)$\lambda =1$ event per year, (B)$\lambda =0.1$ yr⁻¹ (once per decade), and (C)$\lambda =0.01$ yr⁻¹ (once per century), on average, until the uncertain time T the eruption ends. Thin dashed black curve shows probability of the complimentary event, zero PDCs in s years, as the sum of three parts: the probability the eruption has already paused by time t (red curve); and the probabilities that the eruption pauses during the next s years (blue dash-dot curve) or after the next s years (green dashed curve) without an intervening PDC.

FIGURE 4

FIGURE 4. Solid black curve shows forecast probability of at least one PDC in $s=5$ years, starting $t=7$ years after the most recent eruption and assuming an activity rate of λ events per year (on average) until the uncertain time T when the eruption ends, plotted for $0\le \lambda \le 0.5$. Peak probability of $P=0.0947$ is achieved at $\lambda =0.085$, about one event per decade. Results are similar for other choices of s in the range $0\le s\le 20$ years. Other curves have same meaning as in Figure 3.

The Model Behind the Plots

Let us begin by measuring time starting at the onset of eruption, and take the recorded event times of PDC events with volume exceeding $0.15\cdot {10}^6{\text{m}}^3$ to be $\mathrm{t}=\left\{t_1,\dots ,t_n\right\}$ with $0\le t_1\le \cdots \le t_n$ for some positive integer n (the number of recorded PDCs), all in the interval $\left[0,t\right]$ from the eruption onset 0 to the present time t. Our data (Ogburn and Calder, 2012; Ogburn and Calder, 2017) currently include $n=931$ PDCs with runout of at least 1 km, ending $t=7$ years ago. Denote by T the uncertain time of the eruption’s end; necessarily $T>t_n$ but we don’t know whether or not T exceeds the present time t, i.e., whether or not the eruption remains active. We suppose that PDC events arrive at a Poisson rate of $\lambda >0$ during $\left(0,T\right]$, and zero outside that interval. Then the likelihood function for T is

$f\left(T|\mathrm{t},\lambda \right)\propto {\lambda }^n{e}^{-\lambda T}{1}_{\left\{T>t_n\right\}}$

(where the notation $1_E$ denotes the indicator function of an event E, the function that is one if E occurs—here, if $T>t_n$—and otherwise zero.) Utilizing the generalized Pareto $\mathsf{GPa}\left(\alpha ,\beta \right)$ prior distribution of (1), the posterior probability density function for T is

$p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)=c{\lambda }^n{e}^{-\lambda T}{\left(1+T/\beta \right)}^{-\alpha -1}{1}_{\left\{T>t_n\right\}}$

with normalizing constant (to ensure that $p\left(T|\cdot \right)$ integrates to one, as any pdf must) given by

$c={\beta }^{-\alpha -1}{\lambda }^{-\alpha -n}{e}^{-\beta \lambda }/\text{Γ}\left[-\alpha ,\lambda \left(\beta +t_n\right)\right].$

Here the function $\text{Γ}\left(z,x\right):={\displaystyle {\int }_x^{\infty }{t}^{z-1}{e}^{-t}dt}$ of two variables $z\in \mathrm{ℝ}$ and $x>0$ denotes the incomplete Gamma function [Abramowitz and Stegun, 1964, §6.5.3]. Simplifying a little, the posterior density function for T, given $\mathrm{t}=\left\{t_j\right\}$ (the vector of observed PDC times), $\left(\alpha ,\beta \right)$ (parameters of the $\mathsf{GPa}\left(\alpha ,\beta \right)$ distribution), and $\lambda >0$ (the PDC mean annual event rate), is:

$p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)=\frac{e^{-\lambda \left(\beta +T\right)}{\left(1+T/\beta \right)}^{-\alpha -1}}{{\lambda }^{\alpha }{\beta }^{\alpha +1}\text{Γ}\left[-\alpha ,\lambda \left(\beta +t_n\right)\right]}{1}_{\left\{T>t_n\right\}}.(2)$

The values of $\alpha ,\beta$ are estimated in Wolpert et al. (2016) from the entire DomeHaz data set (Ogburn and Calder, 2012; Ogburn et al., 2012; Ogburn and Calder, 2017).

At any time $t>t_n$ the probability of no PDC in the next $s>0$ years under this model, given none in $\left(t_n,t\right]$, is given by a sum of three simple integrals over different time intervals, each available in closed form:

$\begin{array}{l}\mathrm{P}\left[\text{No PDC in first}s\text{years after present time}t|\mathrm{t},\alpha ,\beta ,\lambda \right]=\\ {\displaystyle {\int }_{t_n}^{t}1p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT}+{\displaystyle {\int }_t^{t+s}{e}^{-\lambda \left(T-t\right)}p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT}+{\displaystyle {\int }_{t+s}^{\infty }{e}^{-\lambda s}p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT},\end{array}(3)$

where each integrand is the product of the conditional probability of no PDC in s years (given the uncertain value of T) and the conditional pdf for T. The first of these three terms, shown as a red dotted curve in Figures 3A–C, covers the range $T\le t$, in which the eruption has completed (or paused) before the current time t (and so the conditional probability of no PDC in s years is one). The second, shown as a blue dash-dot line, covers the range $t<T\le t+s$, in which the eruption is still active at the present time t but ends before an additional s years transpire. The third, shown as a green dashed line, covers the range in which $T>t+s$, so the eruption remains active throughout the period of interest. The dashed black line shows their sum (see Eq. 3), the probability $S\left(s,t\right)$ of no PDC in the first s years after the present. Finally, the solid black line shows the probability $1-S\left(s,t\right)$ of the complimentary event, i.e., the forecast probability of at least one PDC event in the next s years, given none in the t years since the last recorded PDC.

As time passes without additional PDCs, the forecast probability that the 1995 SHV eruption has paused or ended will rise, and the forecast probabilities of future events will fall. Figure 5 shows the forecast probability that SHV has paused or ended as a function of the time t since the last PDC, for $7\le t\le 20$ (note $t=7$ in 2020). Forecast probability that the eruption has already paused or ended (shown as a dotted red curve) rises from 67% at $t=7$ years up to 95% at $t=20$, with the forecast probability of a PDC in the next century (solid black curve) dropping to below 2%.

FIGURE 5

FIGURE 5. Solid black curve shows forecast probability of at least one PDC in $s=100$ years, assuming an activity rate of $\lambda =0.1$ event per year, starting t years after the most recent eruption for $7\le t\le 20$. Dotted red curve shows forecast probability that the eruption has already ended by time t, rising from 67% at $t=7$ years (i.e., in 2020) up to 95% at $t=20$ years.

Forecasting and Planning With Dynamic Probabilistic Hazard Maps

Now that we have probabilistic forecasts of future flow events occurring at SHV, one may ask:

If a flow event happens at SHV, what locations might be affected?

It is still useful to think of this question probabilistically, as the flow scenario—flow volume, valleys affected, hot vs. cold flows etc. — are governed by aleatoric uncertainty.

Many locations surrounding SHV have not been inundated by any flows to date, or by just one or two of the several hundred on record since 1995. This does not, however, imply that these locations are immune from inundation by future flows. To study the effect of future flows we will rely on flow simulations [in this work using TITAN2D, see (Patra et al., 2005)] at scenarios consistent with activity at SHV in a probabilistic sense. That is, we can probe scenarios that have some reasonable probability of occurring and exercise TITAN2D at those scenarios to approximate which locations would be inundated by such flows. The approach taken in this work follows (Bayarri et al., 2009; Spiller et al., 2014; Bayarri et al., 2015) and for clarity we summarize important aspects of those contributions in the next subsection.

Probabilities of Pyroclastic Density Current Inundations and Emulator-Based Calculations

Bayarri et al. (2009) explored this question of potential PDC inundation under the assumption that PDC event times follow a stationary Poisson point process [Note, one can consider non-stationary models of event timing in a similar fashion, but for purposes of elucidating our approach, we will follow (Bayarri et al., 2009)]. Let us further consider the case where the aleatory uncertainty in volume and initial angle is described by a probability density function (pdf), $p\left(V,\varphi \right)$. (Note, we intentionally leave $p\left(V,\varphi \right)$ unspecified at this point, so that we can consider different functions and explore the impacts of uncertainty in these functions). Under these assumptions, the expected number of flows in s years that inundate a specific location, ${\left(x,y\right)}_k$, is given by

$\mathrm{E}\left[\#\text{ flows that inundate }{\left(x,y\right)}_k|\lambda \right]=\lambda s{\displaystyle {\int }_{\left(-\mathrm{180,180}\right]\times {\mathrm{ℝ}}_{+}}{1}_{\{h_k>h_{\text{ crit}}\}}\hspace{0.17em}p\left(V,\varphi \right)dVd\varphi }(4)$

where $h_k$ is the maximum flow depth at map location ${\left(x,y\right)}_k$, $h_{\text{ crit}}$ is the critical threshold whose exceedance indicates “inundation,” and $1_E$ is again the indicator function that takes on the value one if the event E occurs (here, that location ${\left(x,y\right)}_k$ is inundated) and otherwise zero. The probability of inundation is then given by one minus the probability that no flow leads to inundation at ${\left(x,y\right)}_k$ in s years, namely

$\begin{array}{l}\mathrm{P}\left[\text{location }{\left(x,y\right)}_k\text{ is inundated over }s\text{years }|\lambda \right]\\ =1-\text{exp}\left[-\lambda s{\displaystyle {\int }_{\left(-\mathrm{180,180}\right]\times {\mathrm{ℝ}}_{+}}{1}_{\left\{h_k>h_{\text{ crit}}\right\}}\hspace{0.17em}p\left(V,\varphi \right)dVd\varphi }\right].\end{array}(5)$

From a practical standpoint, evaluating the indicator function is a significant barrier to calculating Eq. 5. Maximum simulated flow height $h_k$ is not available in closed form (for TITAN2D or other partial differential equation based models solved over digital elevation models), so such integrals must typically be approximated using Monte Carlo (MC) simulation. A straight-forward MC simulation may require thousands to millions of unique samples and hence indicator function evaluations (e.g., $h_k\left(V_j,{\varphi }_j\right)$ with $V_j,{\varphi }_j\sim p\left(V,\varphi \right)$). This approach is a prohibitive computational expense, as a single TITAN2D run takes O(mins–hours) on a super-computer.

To overcome this obstacle, Spiller et al. (2014) and Bayarri et al. (2015) built statistical models of the physical model of flowing mass. Specifically they fit Gaussian process (GP) emulators to TITAN2D output at each location ${\left(x,y\right)}_k$. This approach requires an initial set of N training runs, typically several hundred to several thousand. (Note that all map locations share this common set of N training runs, and the following approach can be applied in parallel at each location). Input values $\left\{\mathrm{V},\varphi \right\}$ for these TITAN2D runs, called design points, are sampled in a space-filling manner to cover the support of any potential pdf $p\left(V,\varphi \right)$. The responses, $h_k(\cdot )$ of the design points are collected into an $N\times 1$ response vector, ${\mathrm{y}}_k$. Then a Gaussian process is fit to these training runs, yielding

${\tilde{h}}_k\left(V^{\text{*}},{\varphi }^{\text{*}}\right)\sim \mathrm{MVN}\left[\mu \left(V^{\text{*}},{\varphi }^{\text{*}}\right),\mathrm{R}|{\mathrm{y}}_k,\mathrm{V},\varphi \right],$

where $\mathrm{MVN}$ denotes the multivariate-normal distribution, $\left(V^{\text{*}},{\varphi }^{\text{*}}\right)$ are an untested input pair, $\mu (\cdot )$ is the mean of the GP, and $\mathrm{R}$ is an $N\times N$ covariance matrix for the GP at the design points. “Fitting” a GP entails finding good choices of the mean function, correlation structure, and correlation parameters; further details can be found in the references we cite (additionally, and in a more general context, (Rasmussen and Williams, 2006; Santner et al., 2018; Welch et al., 1992) provide an excellent background on this topic). It is worth mentioning that ${\tilde{h}}_k$ acts as an interpolator, so at all design points ${\tilde{h}}_k\left(\mathrm{V},\varphi \right)=h_k\left(\mathrm{V},\varphi \right)$. For an untested input scenario, $\left(V^{\text{*}},{\varphi }^{\text{*}}\right)$, ${\tilde{h}}_k\left(V^{\text{*}},{\varphi }^{\text{*}}\right)$ is a draw from a GP whose variance reflects the uncertainty we introduce by emulation, i.e., by taking ${\tilde{h}}_k(\cdot )$ as a surrogate for $h_k(\cdot )$. The benefit of utilizing this approximation in Eq. 5 is that each indicator function evaluation in a MC sample is now a computationally “free” function evaluation available in milliseconds, instead of a O(min–hours) computation. Additionally, once the GP has been fit, no additional runs of TITAN2D are necessary.

Dynamic Probabilistic Pyroclastic Density Current Inundation Forecast Maps

To address the motivating question—of the extent to which the area surrounding the flanks of SHV in southern Montserrat could be impacted by PDC hazards—we now seek not just a probability of a PDC occurring in the next s years, but the probability of a PDC occurring that is big enough, and oriented close enough toward a specific location ${\left(x,y\right)}_k$, to inundate that location in the next s years. Effectively this means combining Eqs 3 and 4. In particular, the expected number of PDCs in s years with annual rate λ is $s\lambda$, and the expected number of PDCs in s years that lead to inundation at location ${\left(x,y\right)}_k$ is given by Eq. 4. Thus, the arguments of the exponential functions in the second and third integrals in Eq. 3 get replaced by Eq. 4, and we have

$\begin{array}{c}\mathrm{P}\left[\text{A flow inundates}{\left(x,y\right)}_k\text{in the s years after present time }t\right|\mathrm{t},\alpha ,\beta ,\lambda ]=\\ 1-[\begin{array}{l}{{{\displaystyle \int }}^{}}_{t_n}^{t}1p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT\\ +{{{\displaystyle \int }}^{}}_t^{t+s}\text{exp}(-\lambda \left(T-t\right){\displaystyle {\int }_{\left(-\mathrm{180,180}\right]\times {\mathbb{R}}_{+}}{1}_{\{h_k>h_{\text{crit}}\}}p\left(V,\varphi \right)dVd\varphi })\times p\left(T|\mathbf{t},\alpha ,\beta ,\lambda \right)dT\\ +{{{\displaystyle \int }}^{}}_{t+s}^{\infty }\text{exp}(-\lambda s{\displaystyle {\int }_{\left(-\mathrm{180,180}\right]\times {\mathbb{R}}_{+}}{1}_{\{h_k>h_{\text{crit}}\}}\times p\left(V,\varphi \right)dVd\varphi })p\left(T|\mathbf{t},\alpha ,\beta ,\lambda \right)dT\end{array}].\end{array}(6)$

We now consider various probabilistic scenarios for flow volume and initial orientation. In Wolpert (2018), motivated by the near-linearity of a log-log plot of the number of observed PDCs of volume exceeding v vs. v, we fit a Pareto probability distribution to PDC volumes at SHV, $\mathrm{V}=\left\{V_1,\dots ,V_n\right\}$ at or exceeding volume $5\cdot {10}^6{\text{m}}^3$. We further assume that this holds for flows down to the smallest considered, $\mathit{\epsilon }=0.15\cdot {10}^6{\text{m}}^3$—a reasonable approximation, since small volume flows are only hazardous to areas very near the source. Now we consider the aleatory scenario model for the volume of flows at SHV, given by

$p\left(V|{\alpha }_V\right)=\frac{{\alpha }_V}{\mathit{\epsilon }}{\left(\frac{V}{\mathit{\epsilon }}\right)}^{-\left({\alpha }_V+1\right)}{1}_{\left\{V>\mathit{\epsilon }\right\}}(7)$

where the uncertain “shape” parameter ${\alpha }_V>0$ must be estimated from the data (see below). Further, in Wolpert et al. (2018), we show that there is not strong evidence against independence for volumes and initial directions of flow events. Thus we take $p\left(V,\varphi \right)=p\left(V\right)p\left(\varphi \right)$ to be of product form and here we consider three cases for $p\left(\varphi \right)$: that flows are distributed uniformly over all possible angles; that flows are initiated exactly to the east; or that they are initiated exactly to the northwest, i.e., we take $p\left(\varphi \right)$ to be uniform on $\left(-\mathrm{180,180}\right]$, or a point mass at $0^{\circ }$, or one at ${135}^{\circ }$, respectively, for angles ϕ increasing counter-clockwize from due east. These three choices represent 1) no preference for a particular collapse direction, or high uncertainty about the directionality, here represented by a uniform distribution of initial angles; 2) collapse direction and flow due east, down the Tar River Valley to the sea, which has been the dominant flow direction throughout the eruption (in part due to the substrate topography slope); and 3) collapse direction and flow toward the northwest, in order to understand the worst-case scenarios for the Belham Valley where there is proximity to population and to sand mining activities.

Note that Eq. 6 depends on the choice of the shape parameter ${\alpha }_V$ through the Pareto probability density function. Likewise it depends on α, β, and λ from the posterior distribution of eruption duration. Although there is very little data (from SHV or other volcanoes) to constrain the low-level frequency λ, we see the maximum probability of an event occurs for $\lambda \approx 0.10$ flows per year for any combination of current time, t, and forecast length, s (see Figure 4 for $t=7$ years and $s=5$ years). Thus we use $\lambda =0.10$ flows per year for all forecasts as probabilities calculated with that value will offer a conservative bound for all forecasts even if the “true” frequency is as high as 1 flow per year or as low as 1 flow per 50 years. This leaves choices for ${\alpha }_V$, α, and β which, respectively, have maximum likelihood estimates of ${\alpha }_V=1.01$, $\alpha =0.58$, and $\beta =0.60$ years. One can imagine plugging these numbers into Eq. 6, and proceeding with a Monte Carlo simulation to estimate Eq. 6; this will reflect the aleatoric uncertainty associated with the randomness of volcanic processes, but not the epistemic uncertainty arising from our lack of certainty about model parameter values. Alternately, one can account for aleatoric and epistemic uncertainty in the models for T and V by treating the parameters in the associated probability distribution functions as random variables. Following the Bayesian paradigm we use probability theory to describe both kinds of uncertainty. To help distinguish them, we systematically use “$p(\cdot )$” to denote probability density functions (pdfs) describing aleatoric uncertainty associated with the inherent randomness of the governing physical processes, and “$\pi (\cdot )$” to denote pdfs describing epistemic uncertainty associated with our imperfect knowledge of model parameters describing these processes. Details of this approach are explored in the next subsection.

We create dynamic probabilistic hazard maps by repeating the calculation in Eq. 6 for every location on the map (in the maps presented here, we do this for each point on a $200\text{m}$ grid covering roughly the southern half of Montserrat). We then repeat it for different choices of angle distribution $p\left(\varphi \right)$ and forecast time length s. Again, the key to making these calculations tractable is that we rely on a single initial design set of roughly 2,000 TITAN2D runs at different scenarios covering the support of $p\left(V,\theta \right)$ and build emulators for each location, e.g., ${\tilde{h}}_k\left(V,\varphi \right)\approx h_k\left(V,\varphi \right)$. With these emulators in hand and by running MC simulations for Eq. 6 at each location in parallel, the computational cost of our MC simulations is very low (roughly 1 min per map on a laptop computer).

Results

Choice of appropriate forecast timescale should be linked directly to the types of decisions that need to be made. Decisions around short-term access to sites of interest—for example, tourist access or sand mining, both of which have significant economic and livelihood implications—could be based on short-term, annual, forecasts. Such forecasts can be updated to reflect the extending time frame of relative quiescence, and should provide a defensible quantification of background hazard level for visitors. We would recommend this approach over using a longer-term forecast map in this context. Alternatively, for decisions such as investment in infrastructure and longer term development, a time scale consistent with a return on the investment might be more appropriate. In these circumstances it might be more appropriate to use 5 or even 20 years forecast maps.

Figure 6 assumes a uniform probability density function for $p\left(\varphi \right)$, $t=7$ years since the most recent PDC. The three panels show forecast snapshots in time, with $s=1,s=5$ and, $s=20$ years, respectively, in panels (A–C). Directionality of a potential flow event is random in an aleatoric sense and should be described by a probability density function (pdf). It is a useful exercise to think about different possible probabilistic descriptions of directionality depending on available data, expert elicitation, or interest in worst-case scenarios. To that end, we have illustrated this process for three separate choices of initiation angle pdfs: uniform (Figure 6), east (Figure 7), and northwest (Figure 8).

FIGURE 6

FIGURE 6. Probability of PDC inundation forecasts plotted on a ${\text{log}}_{10}$ color scale, under the assumption of a uniform distribution for $p\left(\varphi \right)$ over the entire range $-{180}^{\circ }<\varphi \le {180}^{\circ }$. Inundation probability contours are plotted for visual reference. Panel (A) is a forecast for $s=1$ year, (B) for $s=5$ years, and (C) for $s=20$ years into the future. Location of Plymouth is indicated by small triangle due west of SHV.

FIGURE 7

FIGURE 7. Probability of PDC inundation forecasts plotted on a ${\text{log}}_{10}$ color scale, under the assumption of a distribution for $p\left(\varphi \right)$ that is a point mass at $\varphi =0^{\circ }$, i.e., due east via the Tar River Valley to the sea. Inundation probability contours are plotted for visual reference. Panel (A) is a forecast for $s=1$ year, (B) for $s=5$ years, and (C) for $s=20$ years into the future.

FIGURE 8

FIGURE 8. Probability of PDC inundation forecasts plotted on a ${\text{log}}_{10}$ color scale, under the assumption of a distribution for $p\left(\varphi \right)$ that is a point mass at $\varphi ={135}^{\circ }$, i.e., northwest toward the Belham Valley. Inundation probability contours are plotted for visual reference. Panel (A) is a forecast for $s=1$ year, (B) for $s=5$ years, and (C) for $s=20$ years into the future.

Quantifying Uncertainty in Hazard Forecasts

Equation 3 presents the probability of no PDC in the next s years starting t years after the eruption onset, for a specific value of the uncertain parameter vector $\theta$ that determines the prior distribution of T (here $\theta =\left(\alpha ,\beta \right)$). To reflect epistemic uncertainty in hazard forecasts, we integrate this with respect to a density function $\pi \left(\theta \right)$ governing prior uncertainty about $\mathit{\theta }$, then subtract from one, to get:

$\begin{array}{c}\mathrm{P}\left[\text{At least one PDC in first}s\text{years after present time t}|t,\lambda \right]\\ =1-{\displaystyle {\int }_{{\mathrm{ℝ}}_{+}^2}\left\{\begin{array}{l}{\displaystyle {\int }_{t_n}^{t}1p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT}\\ +{\displaystyle {\int }_t^{t+s}{e}^{-\lambda \left(T-t\right)}p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT}\\ +{\displaystyle {\int }_{t+s}^{\infty }{e}^{-\lambda s}p\left(T|\mathrm{t},\alpha ,\beta ,\lambda \right)dT}\end{array}\right\}\pi \left(\mathit{\theta }\right)d\mathit{\theta }.}\end{array}(8)$

Results

Here we explore the effects of epistemic uncertainty on both forecasting the probability of a PDC occurring and on the probability of a PDC inundating a particular location. Figure 9A shows a scatter plot of samples from $\pi \left(\mathit{\theta }\right)$. Using these samples in a MC simulation of Eq. 8 yields an estimate for the probability of a future flow as a function of time, as is plotted by the pink curve in Figure 9B. This curve is analogous to the thick black curve in Figure 3B. Another approach to explore the impacts of epistemic uncertainties on forecasts is to again sample $\mathit{\theta }\sim \pi \left(\mathit{\theta }\right)$, but instead of averaging over the probabilities of a PDC corresponding to each $\mathit{\theta }$ (e.g., the term calculated within the large brackets on the right hand side of Eq. 8) as one would for MC, collect those probabilities and visualize them as a histogram. To explore how these histograms of probabilities evolve in time, we calculate one for each forecast year and fit a kernel density estimate to each histogram. These posterior probabilities reflect the epistemic uncertainty in eruption duration on PDC forecasts and are plotted as a color map in Figure 9; the mode of this posterior is also plotted in time as a white curve. For the first five years, the spread of PDC forecasts is rather narrow and the mean of this posterior distribution lines up its mode. This behavior then transitions to a steady state forecast after about 20 years. At that point, the spread of PDC forecasts reflecting epistemic uncertainties in the duration model goes from a probability of roughly 0.06–0.16, and this posterior distribution is skewed to the left (see Figure 9C), resulting in slightly larger forecasts when using the posterior mode rather than when using the posterior mean.

FIGURE 9

FIGURE 9. Top Left (A): a scatter plot of posterior Monte Carlo samples from $\pi \left(\alpha ,\beta \right)$ reflecting uncertainty in the duration model, i.e., samples from the posterior probability distribution for T (Wolpert et al., 2016). Each sample (plotted point) of $\left(\alpha ,\beta \right)$ represents a single possible shape for the duration distribution. Top Right (B): Probability density function estimates (in color) of flow forecast probabilities (left axis) as they evolve in forecast time s (bottom axis). Uncertainty in these forecasts reflects uncertainty in the duration model. Bottom (C): Histogram of (normalized) pdf for probability of at least one PDC in next $s=5$ years, following $s=7$ years of quiescence—a slice from panel (B) at abscissa $s=5$. Note left skewness [also visible in (B)].

To understand how epistemic uncertainty affects the forecast of inundation probabilities, we will focus on forecasts for the Plymouth area, take $\mathit{\theta }=\left\{\alpha ,\beta ,{\alpha }_V\right\}$, and integrate Eq. 6 against $\pi \left(\mathit{\theta }\right)$. This approach yields:

$\begin{array}{l}\mathrm{P}\left[\text{A flow inundates Plymouth in the s years after present time t }\right|\mathrm{t},\lambda ]=1-{\displaystyle {\int }_{{\mathbb{R}}_{+}^3}[{\displaystyle {\int }_{t_n}^{t}1} p(T\mid \mathrm{t},\alpha ,\beta ,\lambda )dT}\\ {\displaystyle +{\int }_t^{t+s}\exp }(-\lambda (T-t){\displaystyle \underset{(-\mathrm{180,180}]\times {\mathbb{R}}_{+}}{\int }{\mathrm{1}}_{\{h_k>h_{\text{crit}}\}}p(V\mid {\alpha }_V)p(\varphi )dV\hspace{0.17em}d\varphi })p(T\mid \mathrm{t},\alpha ,\beta ,\lambda )dT\\ {\displaystyle +{\int }_{t+s}^{\infty }\exp }(-\lambda s{\displaystyle \underset{(-\mathrm{180,180}]\times {\mathbb{R}}_{+}}{\int }{\mathrm{1}}_{\{h_k>h_{\text{crit}}\}}p(V\mid {\alpha }_V)p(\varphi )dVd\varphi })p(T\mid \mathrm{t},\alpha ,\beta ,\lambda )dT ]\pi (\theta )d\theta .\end{array}(9)$

Just as in the flow forecasting case, we calculate the probability of a flow inundating Plymouth by sampling $\theta \sim \pi \left(\theta \right)$, computing the probability in the bracket on the right hand side of Eq. 9 for each sample, and averaging over those probability forecasts in a Monte Carlo simulation. This approach was followed at Plymouth and all locations on the $200\text{m}$ grid of South Montserrat to produce Figures 6–8 which account for both aleatoric and epistemic uncertainty. Location of Plymouth is indicated by small triangle, due west of SHV. Again, as in the flow forecasting case, we can calculate histograms and kernel approximations to the posterior density of inundation probabilities—these are plotted in a color map in Figure 10C. The mean inundation forecast posterior probability is plotted against time as a pink curve, and the mode posterior is plotted as a white curve. The mean curve gives inundation forecasts more than twice as high as those corresponding to the mode posterior curve indicating that the posterior probability distributions are strongly skewed to the right.

FIGURE 10

FIGURE 10. Probability density function estimates (in color) of forecast inundation probabilities at Plymouth (left axis) as they evolve in forecast time s (bottom axis). In (A) uncertainty in these estimates reflects only uncertainty in the posterior probability distribution for the flow volume V. In (B) uncertainty in these estimates reflects only uncertainty in the posterior probability distribution for T, the eruption duration. In (C) uncertainty in these estimates reflects uncertainty in both the posterior probability distribution for the eruption duration T and the flow volume V. Evidently the dominant uncertainty is that in V. In all cases, the pink curve represents the mean of the probability density function estimates, and the white curves represent the mode of the mode of the posterior distribution of the inundation forecasts. In Panel (C) the modes from (A) and (B) are plotted as black dashed curves for reference (Note, color scales are truncated for visualization). Panel (D) shows histogram of posterior pdf for probability of inundation of Plymouth within next $s=5$ years, reflecting uncertainty in both volume V and end-of-eruption time T—a slice from panel (C) at abscissa $s=5$. Mean and 95% credible appear in bold on Table 1.

TABLE 1

TABLE 1. Numbers in the table (note division by one million on the left) are mean probabilities of inundating Plymouth within s years for flows with initial direction $\varphi =0^{\circ }$ (due east), $\varphi ={135}^{\circ }$ (northwest), or ϕ drawn uniformly from $\left(-{180}^{\circ },{180}^{\circ }\right)$, along with 95% posterior credible intervals. Middle column of bottom row (in bold) corresponds to the histogram shown in Figure 10D below, with mean $5.8\cdot {10}^{-4}$ and 95% credible interval $\left(3.8\cdot {10}^{-5},2.4\cdot {10}^{-2}\right)$. Note intervals are asymmetric due to skewness of the distribution. Probabilities increase with interval duration s and are highest for flows directed NW toward Plymouth and lowest for flows directed due east down the Tar River Valley toward the sea.

To explore the effect of epistemic uncertainties in the volume model vs. those in the duration model, we take $\pi \left(\theta \right)$ to be of product form $\pi \left(\theta \right)=\pi \left(\alpha ,\beta \right)\pi \left({\alpha }_V\right)$ and explore uncertainty in ${\alpha }_V$. We find posterior histograms reflecting uncertain in $\pi \left({\alpha }_V\right)$ by repeating the process of calculating probability histograms from samples of $\pi \left({\alpha }_V\right)$, but marginalizing Eq. 9 over $\pi \left(\alpha ,\beta \right)$. The resulting marginal posterior density is plotted as a color map in Figure 10A along with the mode posterior for the marginal posterior in white and the mean in pink. We also look at the other case where we calculate histograms of inundation probabilities based on samples of $\pi \left(\alpha ,\beta \right)$ with Eq. 9 marginalized over $\pi \left({\alpha }_V\right)$. For this case, the resulting marginal posterior is plotted as a color map in Figure 10B along with the corresponding mode posterior in white and the mean in pink. Note, the mean for each of the three cases is identical, but it is insightful to plot the marginalized mode posterior curves on the full inundation posterior color map. We do this with the two black dashed curves in Figure 10C: the upper one corresponding to the mode posterior from 10B, and the lower from 10A. Note over the first 10 or so years, the color map of the full posterior suggests that the posterior is bimodal with a dominant lower mode corresponding to the epistemic uncertainty in the volume model $\pi \left({\alpha }_V\right)$ and a secondary mode dominated by the epistemic uncertainty in the duration model $\pi \left(\alpha ,\beta \right)$.

Discussion

In this work, we present an evidence-based approach for probabilistic hazard assessment that combines available data on flow events from SHV, statistical modeling, and flow simulations while also accounting for attendant uncertainties. We develop and apply this methodology to assess threats posed by inundation from infrequent, possibly post-eruption unrest flows that could affect the flanks of SHV and parts of southern Montserrat. To do so, we consider a very low-frequency background level of activity that is balanced with the chance that the activity has completely ceased or will cease at some point in the future. Under these assumptions, we can forecast the probability of a post-eruption unrest flow occurring over the next s years, for a range of possible values of s. We then combine this model with simulation-based strategies to glean understanding about how different probabilistic scenarios (e.g., flow volume, initial direction) would impact hazard assessment at specific locations of interest on the flanks of SHV.

This work provides direct information on which relevant stakeholders and/or decision-makers could base practical decisions about managing risk including decisions about livelihoods, and to some extent reoccupation and development around the SHV Montserrat. The context for doing this is that lava domes can remain unstable, or become unstable by various mechanisms during eruptive pauses, or even post eruption. Here we assess collapse hazards related to both rare pyroclastic density currents from an inactive lava dome to debris avalanches from flank or sector collapse. Wisdom gained from other Lesser Antilles lava dome complexes indicate that such events are infrequent but are “highly plausible geological scenarios, whose risks must be taken into account” (Boudon et al., 2007). Our approach provides a means for decision-support by helping to address questions regarding the risk associated with short-term accessibility as well as long-term land-use planning, and infrastructure investments. This method is being used by MVO and SAC to determine the residual level of hazard and associated risk around the volcano and to develop a new micro-zonation hazard map for the island, but it could also be used to address questions of risk around still active volcanoes elsewhere.

Of note, the probabilities of inundation for 1 year forecast windows as determined by this method are also in line with those determined independently through the expert estimates and the expert elicitation process undertaken during SAC meetings (Scientific Advisory Committee of the Montserrat Volcano Observatory, 2011; Scientific Advisory Committee of the Montserrat Volcano Observatory, 2019). Although this cannot be considered a validation as the two assessment methods are so different, it is reassuring that the differing methodologies offer similar forecast probabilities.

In its own right, this new statistical method for forecasting the probability of a post-eruption-unrest flows is a potentially powerful approach that could readily be applied to other volcanoes. This model provides a robust and defensible approach for long term hazard assessment as it treats both the eruption duration and flow frequency as uncertain. The model combines a heavy-tailed duration model introduced in Wolpert et al. (2016), based on an extensive record of dome collapse PDC producing volcanoes worldwide, with a low-level background frequency of activity. Results of the model indicate that the background frequency need not be constrained or inferred by volcano or episode specific data, but instead that a conservative “worst-case” frequency can be found, and that rate only varies with the period of quiescence (t in the model) for the volcano under study and the forecast period s. For the current situation at SHV, the conservative frequency is roughly one PDC event in 10 years.

Further, with the use of statistical emulators, this low-frequency model can readily be folded into a simulation-based probabilistic hazard mapping approach as presented here. With this approach, one can efficiently develop a suite of probabilistic hazard maps under different aleatoric scenario models as there is no need for further time-consuming flow simulations. Additionally, the emulator-based strategy for developing probabilistic hazard maps allows one to explore the effect of epistemic uncertainties. Specifically, at SHV, this has allowed us to compare the uncertain “tail” effects for both the eruption duration and flow volumes, each of which follows a heavy-tailed distribution. To this end, we compared histograms of forecast probabilities whose spread occurs as a result of uncertainty in fitting those probabilistic models.

Although the probabilistic models for flow volume, as well as DEMs, are volcano specific, the general combined statistical-simulation modeling approach is entirely transferable to PDC hazard assessment at other volcanoes worldwide. It is also entirely adaptable for use and application to other geohazards (tephra, tsunamis, etc.).

Conclusion

In summary, we presented a new statistical model for assessing pyroclastic density current (PDC) hazards during conditions of post-eruption-unrest and have combined that with flow simulations to develop a probabilistic hazard assessment and uncertainty quantification for rare, but plausible dome and edifice collapses during periods of quiescence at Soufrière Hills Volcano (SHV), Montserrat. The primary takeaways include:

(1) We introduced a new statistical model for forecasting the probability of post-eruption flow events that is a potentially powerful approach that could readily be applied to other volcanoes.

(2) We combined multiple probabilistic models describing aleatory uncertainty with flow simulations by utilizing statistical emulators. As such, we can rapidly create probabilistic hazard maps for a number of scenarios of interest.

(3) This approach enables efficient uncertainty quantification of both aleatoric and epistemic uncertainties associated with a probabilistic hazard forecast.

(4) We focused on PDC hazard assessment over timescales of 1, 5, and 20 years at SHV. Assuming a uniform distribution of a flow’s initial direction, over 1 year the probability of inundation at Plymouth is $150/{10}^6$. Over 5 years the probability of inundation nearly quadruples to $580/{10}^6$. Yet in 20 years, the probability of inundation is $\mathrm{1,100}/{10}^6$—not quite double the 5 years forecast window level.

(5) This approach has the potential to be transferable to other locations, other hazards, and different times scales. Of course details—data collection, aleatoric modeling, hazard process simulations, etc.—are specific to the hazard and location, but building a flexible framework for rapid probabilistic hazard assessment utilizing statistical emulators is a compelling and widely applicable approach to hazard assessment that can help inform decision making.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher. Additionally, the data are available at the following locations: 1) DomeHaz data are available in spreadsheet snapshot format at Ogburn et al. (2012). DomeHaz: Dome-forming eruptions database v2.4 (https://vhub.org/groups/domedatabase). 2) FlowDat data are available in spreadsheet snapshot format at: Ogburn and Calder (2012). FlowDat: Mass flow database v2.2 (https://vhub.org/groups/massflowdatabase).

Author Contributions

ES: overall project management, co-developer of our dynamic modeling framework, led code development efforts. RW: stochastic modeling of PDC events and eruption durations, co-developer of our dynamic modeling framework, assisted in code development. SO: volcanology and geology expertise, data acquisition and curation, figure generation. EC: motivated project; volcanology expertise, deep familiarity with SHV eruption, interaction with stakeholders (MVO). JB: statistical modeling and uncertainty quantification expertise. AP: high-performance computer expertise, co-developer of TITAN2D flow simulator. EP: granular flow expert, co-developer of TITAN2D flow simulator.

Funding

US National Science Foundation grants DMS-1821289-1821311-1821338, DMS-1621853-1622403-1622467, DMS-1638521, and SES-1521855 covered travel expenses, PhD student support, and some researchers’ compensation. UK NRC grants EP/K032208/1 and NERC Standard Grant NE/R011001/1 covered some travel and subsistence. Duke University will contribute some Open Access publication fees through their CODA program.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF. ES would like to thank the Isaac Newton Institute for Mathematical Studies, Cambridge, United Kingdom, for support and hospitality during the program Mathematical and statistical challenges in landscape decision making where work on this paper was undertaken. The authors would like to thank MVO scientists, in particular Paul Cole, Victoria Miller, and Graham Ryan, for discussions which have helped to direct this work. The authors would further like to thank the reviewers including Larry Mastin who performed a U.S. Geological Survey internal review.

Appendix 1: Distributions

This short appendix describes the probability distributions used in this work.

Generalized Pareto
The $\mathrm{GPa}\left(\alpha ,\beta \right)$ is a continuous distribution for positive random variables T, with probability density function (pdf), survival function (or inverse cumulative distribution function (CDF)), and mean given by:

$\begin{array}{c}f\left(t\right)=\left(\alpha /\beta \right){\left(1+t/\beta \right)}^{-\alpha -1},\hspace{1em}\hspace{1em}t>0\\ \overline{F}\left(t\right)=\mathrm{P}\left[T>t\right]={\left(1+t/\beta \right)}^{-\alpha },\hspace{1em}t>0\\ \mathrm{E}\left[T\right]=\{\begin{array}{cc}\beta /\left(\alpha -1\right)& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\text{if }\alpha >1\\ \infty & \hspace{1em}\hspace{1em}\hspace{1em}\text{if }\alpha \le 1\end{array}\end{array}$

for parameters $\alpha >0$, $\beta >0$. For large values of $\alpha ,\beta$ these are close to those of an Exponential distribution $\mathrm{Ex}\left(\lambda \right)$ with

$\begin{array}{c}f\left(t\right)=\lambda \exp \left(-\lambda t\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}t>0\\ \overline{F}\left(t\right)=\mathrm{P}\left[T>t\right]=\exp \left(-\lambda t\right),\hspace{1em}t>0\\ \mathrm{E}\left[T\right]=1/\lambda \end{array}$

with parameter $\lambda =\left(\alpha /\beta \right)$, but for modest values of $\alpha ,\beta$ the Generalized Pareto has much heavier tails than the exponential—i.e., $\mathrm{P}\left[T>t\right]$ falls off much more slowly with increasing t, so very large values of T are more likely. The unitless “shape” parameter α governs the weight of the tails—for large α (say, larger than 10 or 20) the distribution is almost indistinguishable from the light-tailed exponential distribution, while for values smaller than $\alpha \le 2$ the tails are so heavy that the distribution has infinite variance, and for $\alpha \le 1$ even the mean is infinite. This distribution is commonly used to model incomes, long durations, and other phenomena exhibiting heavy tails. We first use it in this paper in the subsection titled Forecasting with a Chance of Ending. Evidence is very strong that heavy-tailed distributions like the $\mathrm{GPa}\left(\alpha ,\beta \right)$ fit observed eruption durations much better than “light tailed” distributions (exponential, gamma, etc.) do. It was used by Sparks and Aspinall (2004) for a similar purpose. Frequently the Generalized Pareto is introduced as a three parameter distribution family, $\mathrm{GPa}\left(\alpha ,\beta ;\mathit{\epsilon }\right)$, taking values in $\left[\mathit{\epsilon },\infty \right)$ for some $\mathit{\epsilon }\in \mathrm{ℝ}$, with survival function

$\mathrm{P}[T>t]={\left(1+\left(t-\mathit{\epsilon }\right)/\beta \right)}^{-\alpha },\mathsf{t}>\mathit{\epsilon };$

this is simply our $\mathrm{GPa}\left(\alpha ,\beta \right)$ distribution plus a constant offset ϵ.

Pareto
A close relative of the Generalized Pareto is the Pareto distribution $\mathrm{Pa}\left(\alpha ,\beta \right)$ with pdf, survival function, and mean given by:

$\begin{array}{c}f\left(t\right)=\left(\alpha /\beta \right){\left(t/\beta \right)}^{-\alpha -1},t>\beta \\ \overline{F}\left(t\right)=\mathrm{P}\left[T>t\right]={\left(t/\beta \right)}^{-\alpha },t>\beta \\ \mathrm{E}\left[T\right]=\{\begin{array}{cc}\alpha \beta /\left(\alpha -1\right)& \text{if }\alpha >1\\ \infty & \text{if }\alpha \le 1\end{array}\end{array}$

It is a special case of the three-parameter $\mathrm{GPa}\left(\alpha ,\beta ;\beta \right)$, or simply the $\mathrm{GPa}\left(\alpha ,\beta \right)$ plus a constant offset of β. In this work it is used to model PDC volumes V, beginning in the subsection Dynamic probabilistic Pyroclastic Density Current Inundation Forecast Maps.

Poisson
The $\mathrm{Po}\left(\mu \right)$ is a discrete distribution for integer-valued count data $N\in \left[\mathrm{0,1,2},\dots \right]$, with probability mass function (pmf) and mean given by

$\begin{array}{c}\mathrm{p}\left(k\right)=\mathrm{P}\left[N=k\right]=\frac{{\mu }^k}{k!}{e}^{-\mu },k\in \left\{\mathrm{0,1,2},\dots \right\}\\ \mathrm{E}\left[N\right]=\mu \end{array}$

for a single parameter $\mu >0$. It is commonly used for modeling the number of events occurring in some specified time interval or spatial area. In that case one often considers a Poisson process $N_t$, the number of events occurring in time interval $\left[0,t\right]$, under an assumption that event counts in disjoint time intervals $\left(t_{j-1},t_j\right]$ are independent. In that case the mean is an increasing function of t,

${\mu }_t:=\mathrm{E}\left[N_t\right].$

The process is called “homogeneous” if ${\mu }_t=\lambda t$ for some constant average event rate $\lambda >0$. In this case the joint pmf for counts $\left\{N_j\right\}$ of events in J intervals $\left(t_{j-1},t_j\right]$ has the simple form

$\mathrm{P}\left[N_j=n_j,1\le j\le J\right]={\lambda }^{n_{+}}{e}^{-\lambda \left(t_J-t_0\right)}{\displaystyle \prod _{j=1}^J{\left(t_j-t_{j-1}\right)}^{n_j}/n_j!}$

where $n_{+}:={{\displaystyle \sum }}^{\text{}}{n}_j$ is the total number of events in the entire time range $\left(t_0,t_J\right]$.In this work Poisson models are first mentioned in the Recent Eruptive History at Soufrière Hills Volcano section, but are used more extensively in the section The Model Behind the Plots.

Circular Distributions
We consider two possible distributions for the initial angle ϕ (measured in degrees counter-clockwize from due east, from the interval $\left(-{180}^{\circ },{180}^{\circ }\right]$):

• Uniform on $\left(-{180}^{\circ },{180}^{\circ }\right]$, with pdf

$f\left(\varphi \right)=\{\begin{array}{cc}1/360& -180<\varphi \le 180\\ 0& \text{else}.\end{array}$

This indicates that flows in all directions are equally likely.

• Point mass at a specific value ${\varphi }_0$ in $\left(-{180}^{\circ },{180}^{\circ }\right]$. This indicates certainty that $\varphi ={\varphi }_0$ (specifically, we consider ${\varphi }_0=0^{\circ }$, due east toward the Tar River Valley, and ${\varphi }_0={135}^{\circ }$, northwest toward the Belham valley).

These are both limiting special cases of von Mizes distribution, a flexible family of circular distributions used to model angles and defined in Wolpert et al., 2018. We discuss these angular distributions in the subsection entitled Dynamic probabilistic Pyroclastic Density Current Inundation Forecast Maps.

Multivariate Normal
The $\mathrm{MVN}\left(\mu ,\text{Σ}\right)$ is the probability distribution of a random vector $Y\in {\mathrm{ℝ}}^d$ for some integer dimension d, with ${\mathrm{ℝ}}^d$-valued mean $\mu =\mathsf{E}\left[Y\right]$ and $d\times d$-dimensional variance matrix $\text{Σ}=\mathsf{E}\left(Y-\mu \right){\left(Y-\mu \right)}^{\mathsf{ T}}$. It is ubiquitous in statistical modeling of multi-variate data. In particular, it is the joint probability distribution for the values $y_k=Y\left[x_k\right]$ of a Gaussian Process (GP) Y, evaluated at any finite collection of design points. We use it in that context, beginning in the subsection Probabilities of Pyroclastic Density Current Inundations and Emulator-Based Calculations.

References

Abramowitz, M., and Stegun, I. A. (Editors) (1964). Handbook of mathematical functions with formulas, graphs, and mathematical tables, Washington, DC: National Bureau of Standards, 1046.

Google Scholar

Ball, J. L., Calder, E. S., Hubbard, B. E., and Bernstein, M. L. (2013). An assessment of hydrothermal alteration in the Santiaguito lava dome complex, Guatemala: implications for dome collapse hazards. Bull. Volcanol. 75, 676. doi:10.1007/s00445-012-0676-z

CrossRef Full Text | Google Scholar

Ball, J. L., Stauffer, P. H., Calder, E. S., and Valentine, G. A. (2015). The hydrothermal alteration of cooling lava domes. Bull. Volcanol. 77, 102. doi:10.1007/s00445-015-0986-z

CrossRef Full Text | Google Scholar

Barclay, J., Johnstone, J. E., and Matthews, A. J. (2006). Meteorological monitoring of an active volcano: implications for eruption prediction. J. Volcanol. Geoth. Res. 150, 339–358. doi:10.1016/j.jvolgeores.2005.07.020

CrossRef Full Text | Google Scholar

Bayarri, M. J., Berger, J. O., Calder, E. S., Dalbey, K., Lunagómez, S., Patra, A. K., et al. (2009). Using statistical and computer models to quantify volcanic hazards. Technometrics 51, 402–413. doi:10.1198/TECH.2009.08018

CrossRef Full Text | Google Scholar

Bayarri, M. J., Berger, J. O., Calder, E. S., Patra, A. K., Pitman, E. B., Spiller, E. T., et al. (2015). Probabilistic quantification of hazards: a methodology using small ensembles of physics-based simulations and statistical surrogates. Int. J. Uncertainty Quantif. 5, 297–325. doi:10.1615/Int.J.UncertaintyQuantification.2015011451

CrossRef Full Text | Google Scholar

Boudon, G., Le Friant, A., Komorowski, J.-C., Deplus, C., and Semet, M. P. (2007). Volcano flank instability in the Lesser Antilles Arc: diversity of scale, processes, and temporal recurrence. J. Geophys. Res. 112, B08205. doi:10.1029/2006JB004674

CrossRef Full Text | Google Scholar

Calder, E. S., Luckett, R., Sparks, R. S. J., and Voight, B. (2002). Mechanisms of lava dome instability and generation of rockfalls and pyroclastic flows at Soufrière Hills Volcano, Montserrat. Geol. Soc. London Mem. 21, 173–190. doi:10.1144/GSL.MEM.2002.021.01.08

CrossRef Full Text | Google Scholar

Capra, L., Macı́as, J. L., Scott, K. M., Abrams, M., and Garduño-Monroy, V. H. (2002). Debris avalanches and debris flows transformed from collapses in the Trans-Mexican Volcanic Belt, Mexico - behavior, and implications for hazard assessment. J. Volcanol. Geoth. Res. 113, 81–110. doi:10.1016/S0377-0273(01)00252-9

CrossRef Full Text | Google Scholar

Clavero, J. E., Sparks, R. S. J., Pringle, M. S., Polanco, E., and Gardeweg, M. C. (2004). Evolution and volcanic hazards of Taapaca Volcanic Complex, Central Andes of Northern Chile. J. Geol. Soc. London. 161, 603–618. doi:10.1144/0016-764902-065

CrossRef Full Text | Google Scholar

Cole, P. D., Bass, V. A., Christopher, T. E., Eligon, C., Murrell, C., Odbert, H. M., et al. (2010). Report No.: OFR 10-02a. Report to the scientific advisory committee on volcanic sctivity at Soufrière Hills Volcano, Montserrat: Report on activity between 28 February 2010 and 31 October 2010, 50.

Google Scholar

Cole, P. D., Neri, A., and Baxter, P. J. (2015). “Hazards from pyroclastic density currents,” in The encyclopedia of volcanoes. New York, NY: Elsevier, 943–956. doi:10.1016/B978-0-12-385938-9.00054-7.

CrossRef Full Text | Google Scholar

Deplus, C., Le Friant, A., Boudon, G., Komorowski, J.-C., Villemant, B., Harford, C., et al. (2001). Submarine evidence for large-scale debris avalanches in the Lesser Antilles Arc. Earth Planet Sci. Lett. 192, 145–157. doi:10.1016/S0012-821X(01)00444-7

CrossRef Full Text | Google Scholar

Druitt, T. H., and Kokelaar, B. P. (Editors) (2002). The eruption of Soufrière Hills Volcano, Montserrat from 1995 to 1999. Geol. Soc. London Mem. 21, 71–91, doi:10.1144/GSL.MEM.2002.021

CrossRef Full Text | Google Scholar

Elsworth, D., Voight, B., Thompson, G., and Young, S. R. (2004). Thermal-hydrologic mechanism for rainfall-triggered collapse of lava domes. Geology. 32, 969. doi:10.1130/G20730.1.

CrossRef Full Text | Google Scholar

Harnett, C. E., Thomas, M. E., Calder, E. S., Ebmeier, S. K., Telford, A., Murphy, W., et al. (2019). Presentation and analysis of a worldwide database for lava dome collapse events: the Global Archive of Dome Instabilities (GLADIS). Bull. Volcanol., 81 (3), 1–17. doi:10.1007/s00445-019-1276-y

CrossRef Full Text | Google Scholar

Kerle, N. (2002). Volume estimation of the 1998 flank collapse at Casita volcano, Nicaragua: a comparison of photogrammetric and conventional techniques. Earth Surf. Process. Landforms. 27, 759–772. doi:10.1002/esp.351

CrossRef Full Text | Google Scholar

Le Friant, A., Harford, C. L., Deplus, C., Boudon, G., Sparks, R. S. J., Herd, R. A., et al. (2004). Geomorphological evolution of Montserrat (West Indies): importance of flank collapse and erosional processes. J. Geol. Soc. London. 161, 147–160. doi:10.1144/0016-764903-017

CrossRef Full Text | Google Scholar

Mastin, L. G. (1994). Explosive tephra emissions at Mount St. Helens, 1989-1991: the violent escape of magmatic gas following storms? Geol. Soc. Am. Bull. 106, 175–185. doi:10.1130/0016-7606(1994)106<0175:eteams>2.3.co;2

CrossRef Full Text | Google Scholar

Matthews, A. J., Barclay, J., Carn, S., Thompson, G., Alexander, J., Herd, R., et al. (2002). Rainfall-induced volcanic activity on Montserrat. Geophys. Res. Lett. 29, 1644. doi:10.1029/2002GL014863

CrossRef Full Text | Google Scholar

Norton, G. E., Watts, R. B., Voight, B., Mattioli, G. S., Herd, R. A., Young, S. R., et al. (2002). Pyroclastic flow and explosive activity at Soufrière Hills Volcano, Montserrat, during a period of virtually no magma extrusion (March 1998 to November 1999). Geol. Soc. London Mem. 14, 467–481. doi:10.1144/GSL.MEM.2002.021.01.21

CrossRef Full Text | Google Scholar

Ogburn, S. E., and Calder, E. S. (2012). FlowDat–Mass flow database. Available at: https://vhub.org/groups/massflowdatabase.

Google Scholar

Ogburn, S. E., and Calder, E. S. (2017). The relative effectiveness of empirical and physical models for simulating the dense undercurrent of pyroclastic flows under different emplacement conditions. Front. Earth Sci. 5, 83. doi:10.3389/feart.2017.00083

CrossRef Full Text | Google Scholar

Ogburn, S. E., Loughlin, S. C., and Calder, E. S. (2015). The association of lava dome growth with major explosive activity (VEI ≥ 4): DomeHaz, a global dataset. Bull. Volcanol. 77, 40. doi:10.1007/s00445-015-0919-x

CrossRef Full Text | Google Scholar

Ogburn, S. E., Loughlin, S. C., and Calder, E. S. (2012). DomeHaz: Dome-forming eruptions database. Available at: https://vhub.org/groups/domedatabase.

Google Scholar

Opfergelt, S., Delmelle, P., Boivin, P., and Delvaux, B. (2006). The 1998 debris avalanche at Casita volcano, Nicaragua: investigation of the role of hydrothermal smectite in promoting slope instability. Geophys. Res. Lett. 33, L15305. doi:10.1029/2006gl026661

CrossRef Full Text | Google Scholar

Pallister, J. S., Schneider, D. J., Griswold, J. P., Keeler, R. H., Burton, W. C., Noyles, C., et al. (2013). Merapi 2010 eruption-chronology and extrusion rates monitored with satellite radar and used in eruption forecasting. J. Volcanol. Geoth. Res. 261, 144–152. doi:10.1016/j.jvolgeores.2012.07.012

CrossRef Full Text | Google Scholar

Patra, A. K., Bauer, A. C., Nichita, C. C., Pitman, E. B., Sheridan, M. F., Bursik, M., et al. (2005). Parallel adaptive numerical simulation of dry avalanches over natural terrain. J. Volcanol. Geoth. Res. 139 (1–2), 1–21. doi:10.1016/j.jvolgeores.2004.06.014

CrossRef Full Text | Google Scholar

Rasmussen, C. E., and Williams, C. K. I. (2006). Gaussian processes for machine learning. Adaptative computation and machine learning series. Boston, MA: The MIT Press, 266.

Google Scholar

Ratdomopurbo, A., and Poupinet, G. (2000). An overview of the seismicity of Merapi volcano (Java, Indonesia), 1983-1994. J. Volcanol. Geoth. Res. 100, 193–214. doi:10.1016/S0377-0273(00)00137-2

CrossRef Full Text | Google Scholar

Reid, M. E., Sisson, T. W., and Brien, D. L. (2001). Volcano collapse promoted by hydrothermal alteration and edifice shape, Mount Rainier, Washington. Geology 29, 779–782. doi:10.1130/0091-7613(2001)029<0779:vcpbha>2.0.co;2

CrossRef Full Text | Google Scholar

Salaün, A., Villemant, B., Gérard, M., Komorowski, J.-C., and Michel, A. (2011). Hydrothermal alteration in andesitic volcanoes: trace element redistribution in active and ancient hydrothermal systems of Guadeloupe (lesser Antilles). J. Geochem. Explor. 111, 59–83. doi:10.1016/j.gexplo.2011.06.004

CrossRef Full Text | Google Scholar

Sampler, A., Quidelleur, X., Boudon, G., Le Friant, A., and Komorowski, J.-C. (2008). Radiometric dating of three large volume flank collapses in the lesser antilles arc. J. Volcanol. Geoth. Res. 176, 485–492. doi:10.1016/j.jvolgeores.2008.04.018

CrossRef Full Text | Google Scholar

Santner, T. J., Williams, B. J., and Notz, W. I. (2018). The design and analysis of computer experiments. 2nd Edn. New York, NY: Springer-Verlag, 436, doi:10.1007/978-1-4939-8847-1

CrossRef Full Text | Google Scholar

Scientific Advisory Committee of the Montserrat Volcano Observatory (2011). Assessment of the hazards and risks associated with the Soufrière Hills Volcano, Montserrat: 16th report of the scientific advisory committee on Montserrat volcanic activity. Technical report, MVO. November 14–16, 2011, Montserrat Volcano Observatory.

Google Scholar

Scientific Advisory Committee of the Montserrat Volcano Observatory (2019). Assessment of the hazards and risks associated with the Soufrière Hills Volcano, Montserrat: 24th report of the scientific advisory committee on Montserrat volcanic activity. Technical report, MVO.

Google Scholar

Sheridan, M. F., Bonnard, C., Carreno, R., Siebe, C., Strauch, W., Navarro, M., et al. (1998). Report on the 30 October 1998 rock fall/avalanche and breakout flow of Casita Volcano, Nicaragua, triggered by Hurricane Mitch. Landslide News., 12, 2–4.

Google Scholar

Siebert, L., Simkin, T., and Kimberly, P. (2011). Volcanoes of the world. 3rd Edn. Oakland, CA: University of California Press, 568.

Google Scholar

Simkin, T., and Siebert, L. (1994). Volcanoes of the world: a regional directory, gazetteer, and chronology of volcanism during the last 10,000 years. 2nd Edn. Tucson, AZ: Geoscience Press, 349.

Google Scholar

Sparks, R. S. J., and Aspinall, W. P. (2004). “Volcanic activity: frontiers and challenges in forecasting, prediction and risk assessment,” in The state of the planet: frontiers and challenges in geophysics, volume 19 of IUGG monograph. Editors R. S. J. Sparks, and C. J. Hawkesworth, New York, NY: AGU, 359–373. doi:10.1029/150GM28. Geophysical Monograph 150

CrossRef Full Text | Google Scholar

Sparks, R. S. J., Barclay, J., Calder, E. S., Herd, R. A., Komorowski, J.-C., et al. (2002). Generation of a debris avalanche and violent pyroclastic density current on 26 December (Boxing Day) 1997 at Soufrière Hills Volcano, Montserrat. Geol. Soc. London Mem. 21, 409–434. doi:10.1144/gsl.mem.2002.021.01.18

CrossRef Full Text | Google Scholar

Spiller, E. T., Bayarri, M. J., Berger, J. O., Calder, E. S., Patra, A. K., Pitman, E. B., et al. (2014). Automating emulator construction for geophysical hazard maps. SIAM/ASA J. Uncertain. Quantif. 2, 126–152. doi:10.1137/120899285

CrossRef Full Text | Google Scholar

Stinton, A. J., Cole, P. D., Stewart, R. C., Odbert, H. M., and Smith, P. (2014). The 11 February 2010 partial dome collapse at Soufrière Hills Volcano, Montserrat. Geol. Soc. London Mem. 29, 133–152. doi:10.1144/M39.7

CrossRef Full Text | Google Scholar

Vallance, J. W., Siebert, L., Rose, W. I., Girón, J. R., and Banks, N. G. (1995). Edifice collapse and related hazards in Guatemala. J. Volcanol. Geoth. Res. 66, 337–355. doi:10.1016/0377-0273(94)00076-S

CrossRef Full Text | Google Scholar

Voight, B. (2000). Structural stability of andesite volcanoes and lava domes. Philos. Trans. R. Soc. London Ser. A 358, 1663–1703. doi:10.1098/rsta.2000.0609

CrossRef Full Text | Google Scholar

Voight, B., and Elsworth, D. (2000). Instability and collapse of hazardous gas-pressurized lava domes. Geophys. Res. Lett. 27(1), 1–4. doi:10.1029/1999GL008389

CrossRef Full Text | Google Scholar

Voight, B., Komorowski, J.-C., Norton, G. E., Belousov, A. B., Belousova, M., Boudon, G., et al. (2002). The 26 December (Boxing Day) 1997 sector collapse and debris avalanche at Soufrière Hills Volcano, Montserrat, Geol. Soc. London Mem. 21, 363–407. doi:10.1144/GSL.MEM.2002.021.01.17

CrossRef Full Text | Google Scholar

Wadge, G., Herd, R., Ryan, G., Calder, E. S., and Komorowski, J.-C. (2010). Lava production at Soufrière Hills Volcano, Montserrat: 1995-2009. Geophys. Res. Lett., 37, L00E03. doi:10.1029/2009GL041466

CrossRef Full Text | Google Scholar

Wadge, G., Robertson, R. E. A., and Voight, B. (Editors) (2014a). The eruption of Soufri`ere Hills Volcano, Montserrat from 2000 to 2010. London, UK: Geological Society, 501.

Google Scholar

Wadge, G., Voight, B., Cole, P. D., Loughlin, S. C., and Robertson, R. E. A. (2014b). The eruption of Soufrière Hills Volcano, Montserrat from 2000 to 2010. Geol. Soc. London Mem. 39, 1–40. doi:10.1144/M39.1

CrossRef Full Text | Google Scholar

Watts, R. B., De Silva, S. L., Jimenez de Rios, G., and Croudace, I. (1999). Effusive eruption of viscous silicic magma triggered and driven by recharge: a case study of the Cerro Chascon-Runtu Jarita Dome Complex in Southwest Bolivia. Bull. Volcanol. 61, 241–264. doi:10.1007/s004450050274

CrossRef Full Text | Google Scholar

Welch, W. J., Buck, R. J., Sacks, J., Wynn, H. P., Mitchell, T. J., and Morris, M. D. (1992). Screening, predicting, and computer experiments. Technometrics 34, 15–25. doi:10.2307/1269548

CrossRef Full Text | Google Scholar

Wolpert, R. L., Ogburn, S. E., and Calder, E. S. (2016). The longevity of lava dome eruptions. J. Geophys. Res. Solid Earth. 121, 676–686. doi:10.1002/2015JB012435

CrossRef Full Text | Google Scholar

Wolpert, R. L., Spiller, E. T., and Calder, E. S. (2018). Dynamic statistical models for pyroclastic density current generation at Soufrière Hills Volcano. Front. Earth Sci. 6, 55. doi:10.3389/feart.2018.00055

CrossRef Full Text | Google Scholar

Yamasato, H., Kitagawa, S., and Komiya, M. (1998). Effect of rainfall on dacitic lava dome collapse at Unzen volcano, Japan. Pap. Met. Geophys. 48, 73–78. doi:10.2467/mripapers.48.73

CrossRef Full Text | Google Scholar