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 · 10 6 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).
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 . In the absence of dome growth during pauses, a long hiatus (such as the current one), or posteruption 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) 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 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 . 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.
The work developed here presents a defensible, evidencebased 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 ≥ t but the rate λ is very small. The distribution for the end-oferuption time T is based on observations of hundreds of domecollapse 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 simulationbased 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 statisticalsimulation 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 domecollapse 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 · 10 6 m 3 metastable dome that developed between November 1995 and February 1998, but which was no longer actively extruding 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 .
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 ; 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 , 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 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, 2017).
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 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 domecollapse 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).
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 Voight et al., 2002). At Casita, an ∼8 ka dacite lava dome complex, a 1.6 · 10 6 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 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 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 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 for some constants α > 0 and β > 0, to be estimated from the data (we used maximum likelihood methods to estimate α ≈ 0.58 and β ≈ 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., 2015); in this extension, the parameters α, β (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 λ hi 261.3 events/year, and a low level of λ 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 · 10 6 m 3 ). The probability that a volcano remaining active at a rate as high as λ ≥ λ 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.

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 ≤ s ≤ 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 λ ≈ 1.0 ≫ 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 λ ≈ 0.01 ≪ 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 λ 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.

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 · 10 6 m 3 to be t {t 1 , . . . , t n } with 0 ≤ t 1 ≤ / ≤ t n for some positive integer n (the number of recorded PDCs), all in Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 the interval [0, t] from the eruption onset 0 to the present time t. Our data 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 λ > 0 during (0, T], and zero outside that interval. Then the likelihood function for T is f (T|t, λ) ∝ λ n e −λT 1 {T> tn} (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 GPa(α, β) prior distribution of (1), the posterior probability density function for T is p T|t, α, β, λ cλ n e −λT 1 + T/β −α−1 1 {T> tn} with normalizing constant (to ensure that p(T|·) integrates to one, as any pdf must) given by Here the function Γ(z, x) : ∞ x t z−1 e −t dt of two variables z ∈ R 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 t {t j } (the vector of observed PDC times), (α, β) (parameters of the GPa(α, β) distribution), and λ > 0 (the PDC mean annual event rate), is: The values of α, β are estimated in Wolpert et al. (2016) from the entire DomeHaz data set 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 (t n , t], is given by a sum of three simple integrals over different time intervals, each available in closed form: 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 ≤ 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 ≤ 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(s, t) of no PDC in the first s years after the present. Finally, the solid black line shows the probability 1 − S(s, t) 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 ≤ t ≤ 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%.

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. 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 nonstationary models of event timing in a similar fashion, but for 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 ≤ λ ≤ 0.5. Peak probability of P 0.0947 is achieved at λ 0.085, about one event per decade. Results are similar for other choices of s in the range 0 ≤ s ≤ 20 years. Other curves have same meaning as in Figure 3.

Probabilities of Pyroclastic Density Current Inundations and Emulator-Based Calculations
FIGURE 5 | Solid black curve shows forecast probability of at least one PDC in s 100 years, assuming an activity rate of λ 0.1 event per year, starting t years after the most recent eruption for 7 ≤ t ≤ 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.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 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(V, ϕ). (Note, we intentionally leave p(V, ϕ) 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, (x, y) k , is given by E # flows that inundate x, y k |λ where h k is the maximum flow depth at map location (x, y) k , h 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 (x, y) 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 (x, y) k in s years, namely P location x, y k is inundated over s years λ 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 straightforward MC simulation may require thousands to millions of unique samples and hence indicator function evaluations (e.g., h k (V j , ϕ j ) with V j , ϕ j ∼ p(V, ϕ)). 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 (x, y) 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 {V, ϕ} for these TITAN2D runs, called design points, are sampled in a space-filling manner to cover the support of any potential pdf p(V, ϕ). The responses, h k (·) of the design points are collected into an N × 1 response vector, y k . Then a Gaussian process is fit to these training runs, yielding where MVN denotes the multivariate-normal distribution, (V * , ϕ * ) are an untested input pair, μ(·) is the mean of the GP, and R is an N × 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 thath k acts as an interpolator, so at all design pointsh k (V, ϕ) h k (V, ϕ). For an untested input scenario, (V * , ϕ * ),h k (V * , ϕ * ) is a draw from a GP whose variance reflects the uncertainty we introduce by emulation, i.e., by taking h k (·) as a surrogate for h k (·). 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 (x, y) 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λ, and the expected number of PDCs in s years that lead to inundation at location (x, y) 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 P A flow inundates x, y k in the s years after present time t t, α, β, λ 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, V {V 1 , . . . , V n } at or exceeding volume 5 · 10 6 m 3 . We further assume that this holds for flows down to the smallest considered, ε 0.15 · 10 6 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 where the uncertain "shape" parameter α 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(V, ϕ) p(V) p(ϕ) to be of product form and here we Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 consider three cases for p(ϕ): 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(ϕ) to be uniform on (−180, 180], or a point mass at 0 + , or one at 135 + , 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 α 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 λ ≈ 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 λ 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 α V , α, and β which, respectively, have maximum likelihood estimates of α V 1.01, α 0.58, and β 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(·)" to denote probability density functions (pdfs) describing aleatoric uncertainty associated with the inherent randomness of the governing physical processes, and "π(·)" 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 m grid covering roughly the southern half of Montserrat). We then repeat it for different choices of angle distribution p(ϕ) 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(V, θ) and build emulators for each location, e.g.,h k (V, ϕ) ≈ h k (V, ϕ). 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 shortterm 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(ϕ), 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).

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 θ that determines the prior distribution of T (here θ (α, β)). To reflect epistemic uncertainty in hazard forecasts, we integrate this with respect to a density function π(θ) governing prior uncertainty about θ, then subtract from one, to get: P At least one PDC in first s years after present time t|t, λ

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 π(θ). 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 θ ∼ π(θ), but instead of averaging over the probabilities of a PDC corresponding to each θ (e.g., the term Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 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.
To understand how epistemic uncertainty affects the forecast of inundation probabilities, we will focus on forecasts for the Plymouth area, take θ {α, β, α V }, and integrate Eq. 6 against π(θ). This approach yields: P A flow inundates Plymouth in the s years after present time t t, λ Just as in the flow forecasting case, we calculate the probability of a flow inundating Plymouth by sampling θ ∼ π(θ), 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 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.
To explore the effect of epistemic uncertainties in the volume model vs. those in the duration model, we take π(θ) to be of product form π(θ) π(α, β) π(α V ) and explore uncertainty in α V . We find posterior histograms reflecting uncertain in π(α V ) by repeating the process of calculating probability histograms from samples of π(α V ), but marginalizing Eq. 9 over π(α, β). 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 π(α, β) with Eq. 9 marginalized over π(α V ). 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 π(α V ) and a secondary mode dominated by the epistemic uncertainty in the duration model π(α, β).

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  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 FIGURE 9 | Top Left (A): a scatter plot of posterior Monte Carlo samples from π(α, β) 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 (α, β) 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)].
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 535567 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 lowfrequency model can readily be folded into a simulation-based probabilistic hazard mapping approach as presented here. With 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 | 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 ϕ 0 + (due east), ϕ 135 + (northwest), or ϕ drawn uniformly from (−180 + , 180 + ), 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 · 10 − 4 and 95% credible interval (3.8 · 10 − 5 , 2.4 · 10 − 2 ). 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. 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 statisticalsimulation 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 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 . DomeHaz: Dome-forming eruptions database v2.4 (https:// vhub.org/groups/domedatabase). 2) FlowDat data are available in spreadsheet snapshot format at: . 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: highperformance 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.