The Longevity of Lava Dome Eruptions

The Longevit y of La v a Dome Eruptions Rob ert L. W olpert 1 , Sarah E. Ogburn 2 , 3 , Eliza S. Calder 4 1 Departmen t of Statistical Science, Duke Univ ersit y , Durham NC, USA, 2 Departmen t of Geology , Univ ersit y at Buffalo, Buffa lo NY, USA, 3 USGS and USAID V olcano D isaster Assistance Progra m, V ancouv er W A, USA, 4 Sc ho ol o f Geoscience s, Univ ersit y of Edin burgh, Edin burg h, UK F eb 08, 2016 Key Poin ts • The durations of la v a dome eru ptions are hea vy-tailed and dep end on comp osition. • Ob jectiv e Ba ye s ian statistical mo dels can describ e the dep endence on comp osition. • Mo del-based forecasts are made for the con tinued d u ration of activ e dome erup tions. Index T erms and Key W or ds 8488 V OLCANOLO GY: V olcanic hazards and risks; 3245 MA THEMA TICAL GEOPHYSICS: P r obabilistic forecasting. Abstract Understanding the duration of past, on-going and future v olcanic eruptions is an impor tan t scientific goal and a key so cietal need. W e present a new metho dology for for e c asting the du- ration of on-going a nd future lav a dome eruptions based o n a database (DomeHaz) recently compiled by the authors. The data base includes dur ation and co mposition for 17 7 such er up- tions, with “eruption” defined as the p erio d enco mpa ssing individual episo des of do me growth along with asso ciated quiesce nt p erio ds during which extrusion pauses but unrest contin ues. In a key finding we show that pr obabilit y distributions for do me eruption dura tions a re b oth heavy-tailed a nd comp osition-dep enden t. W e cons truct o b jective Bay e sian statis tica l mo dels featuring heavy-tailed Generalized Pareto distributio ns with comp osition-sp ecific para meters to make forecasts abo ut the durations of new and o n-going eruptions that dep end on b oth eruption dura tio n-to-date a nd co mposition. O ur Bay esia n predictive distributions reflect b oth uncertaint y ab out mo del para meter v alues (epistemic uncer tain ty) and the na tural v aria bilit y of the geologic processes (aleatoric uncertaint y). The r esults are illustr ated by presenting likely tra jectories for fourteen dome-building eruptions on-g oing in 2015. F ull r e presen tation o f the uncertaint y is presented for tw o key eruptions, Soufri´ ere Hills V olcano in Montserrat (10–139 years, median 35yr ) and Sinabung, Indonesia (1 –17 years, median 4 y r). Uncertainties ar e high, but, imp ortantly , qua n tifiable. This w ork provides for the first time a quan titative and transfer- able metho d and rationa le on which to base long-term planning decisio ns fo r lav a dome forming volcanoes , with wide p otential use and tra nsferabilit y to foreca sts o f o ther types o f eruptions and other a dv er se even ts acro ss the g eohazard s pectrum. 1 1 In tro duction The lik ely duration of ongoing v olcanic eru ptions is a topic of great in terest to vo lcanologi sts, v olcano observ atories, and comm unities aroun d volc ano es. Ho we v er, few stud ies ha v e in vesti- gated eruption durations ( Bebbi ngton 2007, Gunn et al. 2014, Mastin et al. 2009, Simkin 1993, Sp arks and Aspinal l 2004), in p art b ecause the d ata are sparse and distrib u ted. This work uti- lizes a new database of la v a dome eruptions to analyze the d urations of dome b u ilding erup tions using an ob jectiv e Ba ye s statistical mo del, and inv estigates p ossible c haracteristics ( e.g. , magma comp osition) that affect th ose eruption du r atio n s. La v a dome forming eruptions can b e long-l iv ed, and can pro duce violent and difficult-to-forecast activit y includ in g p lin ian and vulcanian explosiv e activit y and menacing p yro clastic density cur- ren ts. The erup tiv e p erio ds asso ciated with domes are notorious for their tendency to cease ex- trusiv e activit y and th en to start up again w eeks, m on ths or y ears later. P erio ds of activ e dome extrusion and gro wth are in tersp er s ed with p er io ds of relativ e qu iesce nce, durin g whic h extrusion ma y slo w or ev en pause altogether, but where p ersistent v olcanic un r est con tin ues and the v olcano do es not return to a long-term state of dormancy . This contribution fo cuses on th e du rations of these longer-term unrest p hases, hereafter termed eruptions , whic h include p erio ds of b oth la v a extrusion and in tervening quiescence. 2 The La v a Dome Database When studying p opulation c haracteristics of vol canic erup tions, it is critical to hav e clear and consisten t criteria for s electing wh ic h d ata to include ( R o dado et al. 2011). Our study includes all eruptions in a new la v a dome d ata base DomeHaz ( Ogburn et al. 2012, v2.2) for whic h d uration is recorded. T h is d atabase con tains information from 419 dome-forming episo des that comprise 228 eruptions at 127 vo lcano es. F or most eruptions the information includ es dur atio n of eru ptions, p erio ds and pauses of d ome gro wth; extrusion rates; and the timing and magnitude (VEI) of an y asso ciated large explosions. Entries were collected systematically from p eer-review ed sources, v olcano observ atory data sources, the S mithsonian Institution Global V olcanism Pr ogram (GVP) database ( Glo b al V olc anism Pr o gr am 2013), the Bullet in of the Global V olcanism Net work (B GVN) ( V enzke et al. 2002), the Large Magnitude Explosiv e V olcanic Eru p tions (LaMEVE, version 2) database ( Cr oswel ler et al. 2012), and N ewha l l and Melson (1983). There is alw a ys am biguit y in w hat p recisely constitutes a single “eruption”. Simkin et al. (1981) treated inactiv e p erio ds of three mont hs or less as a pause in an eruption, b ut to ok longer p erio ds without activit y to b e gaps b etw een eru ptions. Dome-forming eru ptions, w hic h can b e long-liv ed and often feature cyclical dome gro wth episo des, are particularly difficult to c haracterize consisten tly . F ollo wing Ogbu rn et al. (2015), w e treat as a single erup tion an y p er io d du ring which 1. The vo lcano is describ ed as “con tinuously activ e” in the literature, and/or 2. Dome quiescence lasts less th an 2 y ears, and/or 3. F r equ en t or con tinuous signs of unrest are rep orted throughout dome quiescence, without an apparen t retur n to bac kgroun d lev els. As an example, this analysis will classify a p ersisten tly activ e system lik e Merapi to ha ve b een con tinuously activ e since 19 August, 1768. V ery few cases r ely on only one of these criteria, rather, 2 there are usually m ultiple lines of evidence th at supp ort classifying discrete episodes of dome gro wth as a single eruption. W e id en tify 177 dome-forming eruptions (of the 228 tota l) with rep orted d urations in DomeHaz, 14 of wh ich were still ongoing eruptions at the time of pu b licat ion. Most completed DomeH az eruptions (89%, 145 / 163 ) lasted less than 6 ye ars. Of the remainin g 18 longer-li v ed completed eruptions, ho wev er, 14 lasted longer th an 10 years, and 9 lasted longer than 20 yea rs. All but one (Sinabung) of the 14 ongoing eruptions has lasted 10 y ears or more, half of them ov er 20 y ears, and 5 o ver 50 yea r s. V er y long-liv ed dome-forming eruptions ha v e o ccurred at S an ta Mar ´ ıa (San tiaguito) V olc ano, Guatemala (92 yea rs, and ongoing); Sangay V olcano, Ecu ador (187 ye ars, ending in 1916) ; and Merapi V olcano, Ind onesia (246 y ears, and ongoing). Whereas the San ta Mar ´ ıa la v a dome eruption has b een extru ding almost con tinuously since 1922 , most other long- liv ed eruptions are c h aracterized by frequ en t p auses lasting up to sev eral y ears. 3 A Statistical Mo del F ollo w ing Sp arks and Aspinal l (2004), we mo del the dur ati on of la v a dome er u ptions as rand om v ariables with the generalized Pareto distribution. This h ea v y -tailed distrib ution is supp orted empirically by th e near linearit y of the log-log plot in Figure (1) of frequ ency-v ersus-dur ati on for the eruptions considered here. A f ormal go o dness-of-fit test for the P areto distribution ga ve a test statistic v alue of χ 2 = 14 . 4 on 10 degrees of freedom for a P - v alue of 0 . 1557 , far from rejecting th is h yp othesis. In con tr ast, a test of the Exp onentia l distrib ution gav e χ 2 = 151 . 3 on 11 degrees of freedom, for P = 8 . 25 · 10 − 27 sho win g this more commonly-used sur viv al distrib ution is completely inconsisten t with our d ata . Th e Paret o mod el is also sup p orted theoretically by the statistical theory of extreme v alues, whic h asserts that the exceedances of sufficient ly high thresh olds of indep endent replicate s from a wide r ange of pr obabilit y d istributions will ha ve a generalized P areto distribution (see the d iscussion of th e “p eaks o ver thr eshold” approac h to extreme statistic s in Coles 2001, Ch. 4). W e p resen t Ba yesian mo dels for forecasting in Section (5.2), to ac hieve forecasts that reflect all sour ces of mo del uncertain ty . Compu tat ional details are presente d in App endix A. 3.1 Generalized Pareto Mo dels Sp arks and Aspinal l (2004) fit a generalized P areto distribu tion to data consisting of a selection of 137 dome-building eruptions tak en from the Sm ith s onian Institution database ( Simkin and Sieb ert 1994). W e apply an extension of this mo del to a more recen t collection of 177 dome-building erup- tions ta k en fr om DomeHaz ( Ogburn et al. 201 2 ); the extension r eflects v ariation among the erup tiv e durations for vol cano es of differing comp ositions (d aci tic/rh yo litic, andesitic, basaltic). Note th at the Sp arks and Aspinal l (2004) study al so differs b eca use it consid ers du rations of (shorter) discrete dome-building episo des, while w e study durations of (lo n ger) en tire erup tiv e p eriod s (see Section 2), so b oth the d riving questions and the results are n ot directly comparable. In our p aramete rization a random v ariable (here, the duration T of an eruptiv e ph ase) with the generalized P areto distrib ution can take on an y p ositive num b er, with sur viv al p robabilit y function and pr ob ab ility densit y function giv en f or t > 0 by P [ T > t ] =  1 + t/β  − α , f ( t ) = ( α/β )  1 + t/β  − α − 1 (1) 3 1 2 5 10 20 50 100 200 1 2 5 10 20 50 100 Duration T (y) Number lasting ≥ T y ears MERAPI 1768 → SANGA Y 1728 → SANT A MARIA [SANTIAGUITO] 1922 → SANGA Y 1934 → SEMERU 1946 → BEZYMIANNY 1955 → COLIMA 1957 → Figure 1: F requency of lav a dome eruptions lasting longer than T y ears, v ersu s T , on log-log scale. Near-linearit y su ggests a Pareto distribution for erup tio n du rations. The seven eruptions lasting o ve r 50 y ears, w ith their resp ectiv e start dates, are sp ecified on the fi gure and listed in T able (3) in the App endix. for some dimensionless “shap e” p arameter α > 0 and a “scale” p aramete r β > 0 measured in the same time units as T — h ere, years. W e denote this d istribution b y GPa ( α, β ). F or 0 < α < 1 (as suggested b y our data) the surv iv al pr ob ab ility and densit y function d ecrease so slo wly with increasing t that the mean sur v iv al time is in finite, so that the a ve rage of a gro wing list of eruption durations sh ou ld b e exp ected to grow without b ound. F or that reason we presen t the med ian and quartiles for these distrib utions, which are well-defined and finite for all α , rather than means and v ariances, whic h are n ot. Eruptions whic h hav e already lasted for some considerable time are more lik ely to last a longer additional time than are n ew eruptions. F or that reason it is of particular in terest to find the c onditional distribution of the remaining p erio d of activit y T , for an erup tion that has already lasted for some p erio d s . This to o has a generalized P areto distribution, with up dated p aramete rs: T ∼ GP a ( α, β + s ). The most p opular metho d for estimating u n certain statistica l parameters from d ata is Maxim um Lik eliho o d estimation ( Casel la and Ber ger 1990, § 7.2.2 ), in whic h parameters su c h as α and β are estimated b y the v alues ˆ α and ˆ β f or whic h the Lik eliho od F unction (the j oin t p robabilit y densit y function for all th e observ ed qu an tities, regarded as a function of the mo del parameters) attains its maxim u m v alue. In Ap p endix A we p r esen t the lik eliho od fu nction for b oth fu lly-observed d ata (when b oth the start and end times for an eruption are known) and for censored data (where only a low er b ound is kno wn for th e duration of an eru ption, usu ally b ecause it is still ongoing at the time of the analysis). Bo th fully-observ ed and censored data are present in our dataset. Here we present the resulting forecasts of the remaining duration of eac h of the fourteen vo lcano es in our data s et that are still activ e, as well as f orecasts of total d uration for p ossible new eru p tions of eac h comp osition t yp e. 4 W e exp lore in tw o different w ays the p ossibilit y that the remainin g duration of activit y T for a la v a dome eruption may dep end on observ able co v ariates suc h as the magma comp ositio n in terms of silica con tent X (in p er cent). First, we classify into three grou p s b y the eruption comp osition or, where un a v ailable, charact eristic comp osition for that v olcano, and fit mo dels s eparate ly to eac h class; second, w e introdu ce a log-linear regression mo del in whic h b oth th e shap e parameter α and scale p aramete r β for the P areto du ration distribu tion m a y dep end on the silica con tent X . F ull details are given in Ap p endix A. 4 Results In this section w e fit thr ee v ariations on the mo del of Section 3 to the DomeHaz data: an Aggregate mo del fi tting a sin gle generalized P areto m odel to all 177 eruptions; a Group ed mo del, fitting separate P areto m odels to eac h of thr ee classes of eruption, based on comp osition; and a log- linear Regression mo del in wh ic h th e P areto parameters dep end exp licit ly on silica level s. In eac h case mo del parameters are estimated b y n umerical maximization of the lik eliho od fu nction to find Maxim u m Likelihoo d Estimators (MLEs), and Standard Err ors (SEs) are estimated fr om the in verse Hessian matrix at the MLE. 4.1 Aggregate Mo del MLEs and SEs for the parameters of the aggregate GPa mo del without comp ositional dep enden ce are ˆ α = 0 . 6487 ± 0 . 0132 and ˆ β = 0 . 7018 ± 0 . 055 1 yr . P anel (a) of Figure (2) sh o ws an empirical plot of f ractio n vs. duration on a log-log scale, along with the b est m odel fit, for the 163 v olcanic eruptions in the d ata s et whose erup tiv e ep iso de h ad end ed by 15 Marc h 2015 when the dataset w as constructed. Panel (b) shows a similar plot for all 177 eru ptions, including the 14 ones then ongoing, sh ifted forw ard by the median pr o jected remaining duration ∆ = ( ˆ β + s )[2 1 / ˆ α − 1] (und er the mo del) for an eru ption already lasting dur atio n s years (mark ed with an emp t y diamond), to s + ∆ (filled diamond). In eac h plot b oth d uration t and the fraction with d u ration exceeding t are displa yed on logarithmic scales, so for large t the mo del fit will b e appr o ximately a straigh t line with slop e − α . 4.2 Group ed by Comp osition Bottom ro w of Figure (2 ) shows fits to data separated by magma comp osition in to three classes: mafic, typica lly b asaltic v olcano es; intermediate , includ in g andesitic and basaltic and esitic volc a- no es; and evo lv ed systems, typica lly dacitic or r h yo litic vol cano es. Go od n ess-of-fit tests again sh ow that Pareto mo dels fit eac h of these w ell, while Exp onentia l mo dels (with one exception) do not. Exp onen tial mo dels offer a go o d fit to the t wen t y-nin e completed dacitic/rh y olitic erup tions, all of whic h lasted less than ten y ears and t wen t y-six of wh ic h lasted less than five— suggesting falsely that the hea v y -tailed P areto distrib ution ma y b e unnecessary . T he thirtieth is the on-going eruption at the San ta Mar ´ ıa/San tiaguito vo lcano complex in Guatemala, whose 92-y ear-old duration to date confirms the hea vy-tailed n ature of du r atio n distributions. F or the 29 completed Dacite/Rh y olite eruptions w e fit the simpler Exp onenti al Ex ( λ ) mo del, a limiting case of the GPa ( α, β ) f or large α and β with α/β ≈ λ . Mo del parameters f or the GPa ( α, β ) and Ex ( λ ) distrib utions are estimated separately within eac h class. Both th e plots and th e parameter estimates presented in T able 1 show that the duration distribution d iffers mark edly by comp osition. 5 Duration t Fraction > t 3d 1m 6m 2y 10y 50y 0.005 0.05 0.2 1 Duration of Completed Lav a Dome Eruptions Duration t Fraction > t 3d 1m 6m 2y 10y 50y 500y 0.005 0.05 0.2 1 SINABUNG MERAPI SEMERU KARANGET ANG [API SIAU] IBU KARYMSKY BEZYMIANNY SHIVELUCH COLIMA POPOCA TEPETL SANT A MARIA [SANTIAGUITO] FUEGO REVENT ADOR SOUFRIERE HILLS Duration of Completed & Ongoing Lav a Dome Eruptions (a) (b) Duration t Fraction > t 3d 1m 6m 2y 10y 50y 0.01 0.05 0.2 0.5 Duration of Completed Lav a Dome Eruptions Grouped b y Composition Basaltic Andesitic Dacitic, Rh yolitic Duration t Fraction > t 3d 1m 6m 2y 10y 50y 500y 0.01 0.05 0.2 0.5 Duration of Completed & Ongoing Lav a Dome Eruptions Grouped b y Composition (c) (d) Figure 2: Empirical fraction vs. duration for 163 completed eruptions (left, panels a and c) and for all 177 eruptions (right , panels b an d d), completed eru ptions (soli d) and ongoing eruptions (op en), on log-lo g scale. T op ro w (a,b) shows aggregate mo del without explicit comp osition dep end ence, b ottom ro w (c,d) displays fr act ions (within comp ositional classes) and mo del forecasts f or eac h comp ositional class. Horizonta l bars ind icat e m edian pro jected remaining duration und er our mo del for the 14 sp ecified ongoing eru p tions. 6 Completed by 15 Marc h 2015 Completed & Ongoing Comp osition # α ± SE β ± SE # α ± SE β ± SE Basalt: 40 0.6712 ± 0.0540 0.3698 ± 0.1272 42 0.5 440 ± 0.0390 0.293 2 ± 0.1002 Andesite: 94 1.0900 ± 0.0540 1.2140 ± 0.1738 1 05 0.5769 ± 0.0180 0.5993 ± 0.0746 Dacite/ Rhy olite: 29 Ex ( λ ) , ˆ λ = 0 . 3390 ± 0 . 6 758 30 1. 8920 ± 0.3969 5.044 0 ± 1.9380 ± T able 1: P arameter estimates for generalized Pareto and exp onen tial mo dels, for d ata group ed b y comp osition. Mean d uration is E [ T ] = β / ( α − 1) if α > 1, and infinite for α ≤ 1 for Pa reto and 1 /λ for exp onential . Exp onent ial Ex ( λ ) is limiting case of generalized Pareto GP a ( α, β ) for large α, β with α/β ≈ λ . 4.3 A Log-linear R egression Mo del for Comp ositional Dep endence The generalized P areto p aramete r estimates ˆ α and ˆ β in T able 1 v ary strikingly and consisten tly across comp osition classes— eac h parameter increases with increasing silica conten t. This feature can b e captured in the regression mo del T | X ∼ GP a  αe γ α ( X − 60) , β e γ β ( X − 60)  (2) expressing log-linear dep end ence of the generalized Pareto distribu tio n on the s ilic a con tent X , f or new regression parameters γ α , γ β ∈ R . Suc h a mo del offers the adv ant age o v er th e class-sp ecific mo del of Section (4.2) th at eviden ce from all 177 eru ptions can b e us ed in generating forecasts f or eac h v olcano, ev en those wh ose comp osition class is rare (or absent ) in the DomeHaz database. P arameter estimates for this m odel are ˆ α = 0 . 6923, ˆ β = 0 . 7915 yr , ˆ γ α = 0 . 0447, and ˆ γ β = 0 . 1302, again sh o wing that the generalize d P areto p aramet ers dep end on comp osition (b ecause ˆ γ α and ˆ γ β are s ev eral SEs fr om zero). Because ˆ γ α and ˆ γ β are b oth p ositive , b oth the sh ap e and scale parameters increase w ith increasing silica, α by abou t 4 . 5% and β b y ab out 13% for eac h additional p ercen t silica. 4.4 Whic h mo del is b est? The three mo dels pr esented in Sections 4.1–4.3 h a v e different degrees of complexit y , with t wo free parameters for the Aggrega te mo del, six for the Group ed mo del, and four for the Regression mo del. Tw o traditional ap p roac hes to m o del comparison are the Ak aik e information cr iterion, or “AIC” ( Aka ike 1974), and the Bay esian in f ormatio n criterion, or “BIC” ( Schwarz 1978). Eac h fa vors mo dels that fit the data b etter on a log-lik eliho o d scale, and eac h p enalizes mo dels f or complexit y . The Gr ou p ed and Regression mo dels were comparable usin g the AIC criterion (Group ed AIC= 751 . 61 , Regression AIC= 751 . 97) , and eac h out-p erformed significan tly the Aggregate mo del (AIC= 756 . 52 ). The BIC crite rion, whose co mplexit y p enalties are more sev ere, fav ors the Aggregate model for its s im p licit y (Aggregate BIC= 762 . 87, Regression BIC= 764 . 68, Group ed BIC= 770 . 67). 7 5 F orecasting: Ho w m uc h longer for Soufri` ere Hills and S inabu ng? In this section w e present forecasts for the remaining length of t w o ongoing erup tions, the 20-y ear- long erup tion at Soufri` ere Hills, Montserrat and the more recent on e at Sin abung, In d onesia. Under the generalized P areto mod el GP a ( α, β ), the q th quan tile of the remaining p erio d of activit y for an eruption already activ e for s y ears is P [ T > ∆ q | t , α, β ] = q ∆ q = ( β + s ) h (1 − q ) − 1 /α − 1 i . (3) The er u ption at S ou f ri ` ere Hills V olcano commenced in 1995 and had lasted 7189 days, ab out 19 . 7 y ears, as of 15 Marc h 2015, while that at Sinabun g b egan in 2013 and had then lasted 546 days. F r om Eqn (3) forecasts can b e made for the m ed ian ( q = 0 . 50) and quartiles ( q = 0 . 25 , 0 . 75 ) of the distributions for their con tinued eruption durations, using s = 19 . 7 and s = 1 . 49 yea rs and suitable v alues f or the u ncertain parameters α , β for these eru ptions, based on the historical record t of observ ed eruption du rations of all 177 v olcano es. 5.1 Maxim um Lik eliho o d E stimates T able 2 pr esen ts MLEs for the quartiles and median r emaining activ e time for the S oufri ` ere Hills and Sinabu ng V olcano es as of 15 Marc h 2015, often called “plug-in” estimates b ecause they are tak en from (3) w ith α, β replaced by their MLEs ˆ α, ˆ β . Resu lts are presen ted for three v ariations on our generalize d P areto m o del: one f rom the aggreg ate mo del of S ect ion (4.1) with no comp ositional dep endence, one f or the group ed mo del of Section (4.2) basing eac h forecast on only the eruptions with similar comp osition, and one based on the four-p arameter r egression mo del of Section (4.3) with log-linear silica dep endence. Data in cluded b oth the 163 completed eruptions (whose ent ire duration length is kno wn ) and the 14 ongoing erup tions (for whic h only a lo wer b ound on the duration length is kn o w n). Omitting long-lasting ongoing eruptions w ould distort the evidence by in tro du cing a strong down wa r d bias for estimates based on ly on completed eruptions. Soufri` ere Hills Sinabun g Mo del q 25 q 50 q 75 q 25 q 50 q 75 MLE Aggregate: 11.36 38. 91 152.18 1. 23 4.20 16.42 MLE Group ed: 13.10 47.1 0 203.68 1.3 5 4.87 21.06 MLE Log Linear: 10.54 35. 21 131.02 1. 18 3.94 14.65 Ba y es Log Linear: 10.2 9 35.01 138.56 1.21 4.19 16.70 T able 2: Estimated quartiles for pro j ect ed remaining activ e time (in y ears) at SHV (left) and at Sinabun g (right) as of 15 Marc h 2015, based on all 177 eruptions, includ ing the 14 ongoing. T op three r o ws are plug-in estimates based on MLEs ˆ α, ˆ β ; b ottom ro w sho ws ob jectiv e Ba yes p osterior quartiles for log-linear regression m odel. Pro jected remaining dur ati on of all four teen current ly-activ e dome forming volc ano es, colo r - co ded b y comp osition, are giv en in Figure 3. Estimates are based on the regression mo del of Section (4.3), u sing the MLEs f or the four m o del parameters. 8 1 5 10 50 500 0.0 0.2 0.4 0.6 0.8 1.0 Remaining Activ e Time (y) Probability 1 2 3 4 5 6 7 8 9 10 11 12 13 14 A B−A D A A A BA B A A A BA−A A BA−A MERAPI SEMER U SANT A MARIA BEZYMIANNY COLIMA SHIVELUCH KARANGET ANG FUEGO SOUFRIÈRE HILLS IBU KAR YMSKY POPOCA TEPETL REVENT ADOR SINABUNG Figure 3: Pro jected remaining d u ration of eruptions at all f ourteen cur ren tly-activ e dorm forming v olcano es, based on plug-in p arameter estimates for the log-linear regression mo del of S ectio n (4.3) . Colors indicate comp ositional class. Thick blue line is Soufri` ere Hills V olcano, lo wermost green line is Sinabun g. 9 5.2 Ba yesia n Posterior & Predictiv e Distributions Plug-in forecasts b ased only on p oin t estimates ˆ α , ˆ β of the marginal pr obabilit y P [ T > t ] that a new v olcanic eru ption will last longer than t y ears or the conditional probabilit y P [ T > t + s | T > s ] that an s -y ear old v olcanic er u ption will con tinue at least t m ore ye ars, ma y b e distorted b ecause they ignore uncertain ty ab out the parameters α and β . A Bay esian appr oac h is more satisfactory here b eca use it reflects fully b oth epistemic uncertain t y ab out mo del p aramete r v alues and aleatoric uncertain ty from the natural v ariabilit y of the geologic pro cesses. Ou r obje ctive Bay es appr oac h, using r eference prior d istributions ( Bernar do 1979, Ber ger et al. 2009) rather than s ub jectiv e priors, is describ ed in App endix A.4 . Ob jectiv e Ba yesia n p osterior quartiles for the remaining duration of erup tions at Souf r i ` ere Hills V olca no and Sin ab u ng are present ed in the b ottom ro w of T able 2. P osterior predictiv e probabilities that these erup tions will contin ue for at least t add itio nal ye ars are sho wn in Figure 4, for 0 ≤ t ≤ 50 yr. T he solid red line indicates o verall p osterior probability of contin uing t more y ears, while the width of the 90% prediction in terv al (blue lines) indicates ho w uncertain that forecast is on the basis of av ailable evidence. Th e plug-in estimate, shown as a dashed blac k line, is close to the mean b ut obscures the un certain t y . F orecast probabilit y that remaining duration at Sin abung exceeds ten years has a mean of 33.7%, with a wide 90% range of 17.8%–49 .8%, while forecast pr obabilit y that remaining duration at Soufr i ` ere Hills V olca no exceeds ten y ears has a m ean of 75.5%, w ith a narrow er 90% range of 64.8%–8 4.2%. The difference in width is attributable principally to the longer duration of the current erup tion at Soufri´ ere Hills. The median pro j ect ed remaining du rations are 35 . 01 y ears for S oufri ` ere Hills V olcano and 4 . 19 y ears for Sinabu ng. Thus, there is a 50:50 c han ce th at eac h of these eruptions con tinue m ore (or less) than th ose resp ectiv e spans. Th ese forecasts ma y includ e extended quiescent p erio ds of up to t wo y ears or more (see Section 2 for the criteria) and, in th e case of Soufri` ere Hills V olcano, it is stressed that the forecast is conditional on th e current (2015) quiescence b eing a pause and not already the onset of pr olonge d dormancy . 6 Discussion: Implicat ions for ph ysical eruptiv e pro cesses Although th e statistical analysis undertake n h ere cannot pro vid e direct in dicatio n of sp ecific erup - tiv e pro cesses, the asso ciation found b etw een comp osition and longevit y can aid in constraining some asp ects of th ose p ro cesses or indicating w hic h kind of pro cesses ma y b e most imp ortan t. The relativ e con tr ib utions of shallo w, cond uit-lev el, pro cesses v ersu s deep er, m agma c hamb er lev el, pro cesses in the regulation of b oth extrus ion p erio d icit y and erup tio n d u ration for la v a d ome eruptions is tied to time-scales. S hort time-scale patterns ha ve b een mo deled as fun ctions of shallo w nonlinear cond u it dynamics ( Costa et al. 2007, Denlinger and Hoblitt 1999) whereas long time -scale patterns and er u ption ov erall dur ati ons h a v e generally b een regarded as relating to magma rheolog y and deep er magma c hamber conditions ( Barmin et al. 2002, Melnik and Sp arks 2005). This study is fo cused on taking a longer-term view of eruptions, our definition of whic h may in- clude non-extrusiv e p hases where unrest contin ues. Consisten t with earlie r work ( Sp arks and Aspinal l 2004) on individu al dome extrusions, w e fi nd that eru p tion durations are hea vy -tailed, w ith sub - stan tial d rop-off in duration b et ween 1–5 yea r s. Sp arks and A spina l l (2 004) suggested that five y ears may b e sufficien t for the ma jorit y of conduits to mature and stabilize, and that conduit dynamics ma y co n trol dome eru ption longevit y rather than gradual fr eez ing of shallo w magma. 10 0 10 20 30 40 50 0.0 0.2 0.4 0.6 0.8 1.0 Time T (years) P[ Eruption Lasts T More Y ears ] Remaining Duration of 1995 SOUFRIERE HILLS Erupti on 0 10 20 30 40 50 0.0 0.2 0.4 0.6 0.8 1.0 Time T (years) P[ Eruption Lasts T More Y ears ] Remaining Duration of 2013 SINABUNG Eruption (a) (b) Figure 4: Pr o jected remaining duration of eru ptions at Soufri` ere Hills V olcano a) and Sinabun g b), reflecting u n certain t y ab out mo del parameters. Eac h estimate is based on log-linear regression mo del reflecting evidence from all 177 volca n oes. On e hundred d ra ws from p osterior distr ib ution are sho wn (faint pin k lines), along w ith 90% credible inte r v al (medium blue lines) and p osterior mean (thick red line). Dashed b lac k line sho ws plug-in forecast using estimated parameter v alues, ignoring parameter uncertaint y . Poste rior median p ro jected du ration is 35 . 01 y ears at Sou f ri ` ere Hills V olcano, 4 . 19 y ears at Sinab u ng. F or short-term eru ptions, their results and ours are comparable and thus w e concur with their suggestion that cond uit dynamics ma y con trol longevit y f or these. T he longer duration eru ptions w e stud y , in cluding unr est p erio ds with ou t extru s ion, reflect longer-term s ystem dyn amics that are more like ly regulated by the d eeper magmatic system. F u rther, this work demonstrates that er u ption d urations also v ary by system comp osition— and, by extension, the combined effects of magma rh eolo gy , temp erature, crystal conte n t, etc . Figure 2 (c,d) sh o ws that d ome-forming eruptions p redominan tly o ccur at andesitic systems, and of those 76% do not con tinue for more th an 5 y ears (blac k squares). Evolv ed systems s h o w an abrupt d rop-off in eruption duration after 3 y ears (red d iamonds) whereas more mafic systems sho w earlier-onset deca ys, but slow er duration d rop-off rates (blue triangles). F or ev olve d systems, this may reflect that these dome eruptions are mo dulated b y shallo w-lev el, conduit, p rocesses bu t ma y also, imp ortant ly , reflect th e fact that in some of these cases these ev olv ed domes repre- sen t residual vola tile-depleted melts as late-stage squeeze-ups f ollo wing ma jor explosive erup tio ns ( Ogburn et al. 2015). In essence, many ev olv ed la v a domes are asso cia ted with already-depleted magma c hamb ers— conditions which are less lik ely to b e dr iv ers of pr olo nged activit y . T he long- liv ed Sant a Mar ´ ıa domes are an interesting exception to th is, ho wev er. It has b een su gge sted that the 1902 explosiv e erup tion (comp ositionally related to th e dome la v as) ma y ha ve b een shut down prematurely by collapse of edifice m ate rial in to the ven t region ( Andr ews 2014). The su bsequen t long-liv ed d ome episo des are therefore not sourced from a depleted reservoir as is often the case with other ev olv ed dome eruptions. The slo w er deca y of eruption duration for inte rmediate to mafic systems could reflect that more of these eru ptions are mo dulated by deep er pro cesses, a combination of d eep and shallo w pro cesses, or that there is more v ariability within th ese systems. 11 F orecasting of lik ely remaining eruption dur atio ns is extremely v aluable f or hazard mitigation and risk managemen t, esp ecially for ev acuation an d longer term land-use planning. The present analysis has already allo w ed the USGS V olcano Disaster Assistan t Program (VD AP) to p r o vide information to the Indonesian Geolog ical Agency that the Sin abung eruption w ould most lik ely con tinue f or sev eral more y ears, and th u s has pla yed a r ole in supp ortin g d ecisio ns regarding long term managemen t of ev acuations ( Pal lister 2015). 7 Conclusions La v a dome eruptions sho w con tinuous but h ea vy-tailed distributions that v ary with comp osition. Eruption d uration d riv ers ma y b e d ominan tly mo du lat ed by th e shallo w pro cesses (condu it) for ev olv ed systems, wh ereas intermediate to mafic systems ma y ha ve b oth shallo w and deep er (c h am- b er) signatures to them. The metho dology devel op ed here allo ws for comparison of su rviv al cur ves for differen t ongoing eruptions as w ell as pro jected probabilit y distribu tions for remaining dur atio n s for ongoing eruptions. T he pro jections are based on evid en ce pr ovided from a suite of 177 dome eruptions and synthesize d in an ob jectiv e Ba yesian s tatistical mo del. This analysis illustrates the uncertain ties, whic h are h igh but, most imp ortantl y , quant ifiable. Th is w ork pro vid es, for th e first time, qu an titativ e and transferable metho ds and r ati onale on whic h to base long-term planning decisions for la v a d ome formin g vol cano es. A App endix: Generali zed P areto Mo del Calculations This App endix in cludes d eriv ation of the lik eliho o d functions for all three v ariations of the Ba y esian generalized P areto statistical mo del present ed in Sectio n 3, along with a description of our ob jectiv e Ba y esian pr ior distr ib utions and a p roof that the corresp onding p osterior distrib u tions are prop er. Also includ ed as T able 3 is an ordered list of all dome-building v olcanic erup tions in the DomeHaz database ( Ogburn et al. 2012, v2.2) with du ratio n five y ears or longer. Under the generalized Pareto distribution GPa ( α, β ), the surviv al fun ctio n for th e er u ption duration T of a new vo lcano is P [ T > t ] = (1 + t/β ) − α , (4a) dep ending on t w o parameters: the dimensionless “shap e” parameter α > 0, and the “scale” pa- rameter β > 0, measured in ye ars. T h e m ean is E [ T ] = β / ( α − 1) if α > 1, or infinit y if α ≤ 1; the m edian is β [2 1 /α − 1] for any α > 0. The remaining dur atio n of a v olcanic eruption that has already b een in progress for s y ears also has a generalized Pareto distribu tion, P [ T > t | Er uption s ye ars old] = (1 + t/ ( β + s )) − α , (4b) with the same shap e parameter b ut an up dated scale parameter: T ∼ G P a ( α, β + s ). A.1 Lik eliho o d I: No Comp ositional Dep endence The lik eliho od function for α, β u p on observing the exact (u ncensored) erup tio n durations { t i : i ∈ I 1 } of some n umber n 1 of completed eruptions, and also the r igh t-censored observ ations { t i : i ∈ I 0 } 12 of th e form “ T > t i ” fr om some num b er n 0 of eru ptions that we re still con tinuing at the time of data collection, is f ( t | x , α, β ) = Y i ∈ I 1 h ( α/β )  1 + t i /β  − α − 1 i Y i ∈ I 0  1 + t i /β  − α where I 0 and I 1 are sets indexing the censored and non-censored eru ptions, resp ectiv ely . The negativ e log like liho od is ℓ ( t | x , α, β ) = X ( α + δ i ) log (1 + t i /β ) + n 1 log( β /α ) , (5) where δ i = 1 for uncensored and δ i = 0 for censored observ ations, and where n 1 = | I 1 | is the n u m b er of u ncensored observ ations and n = | I 0 ∪ I 1 | is the total num b er of observ ations. Maxim um Lik eliho o d Es timates ˆ α , ˆ β can b e found by a t wo-dimensional searc h to m in imize Eqn (5) or, more efficien tly , by a one-dimensional searc h to minimize ℓ  t | x , ˆ α ( β ) , β  where ˆ α ( β ) = n 1 / P log(1 + t i /β ) is the conditional MLE for α , giv en β . Ba y esian estimation is discussed in Section (A.4). A.2 Lik eliho o d I I: Group ed Comp ositional D ep endence Equation (5) can also b e u sed for p aramete r inference for the stratified mo del of Section (4.2), by limiting I 0 , I 1 to those eruptions featuring a comp osition in a sp ecified group of t y p es. P arameters are estimated separately for eac h comp osition group. A.3 Lik eliho o d I I I: Regression Mo deling of Composit ional Dep endence The negativ e log like lih oo d f unction for th e mo del of Eqn (2) in Section (4.3), in wh ic h the gener- alized Pareto parameters α , β are eac h log-linear functions of silica conten t for eac h eruption, can b e expressed as: ℓ ( t | x , α, β , γ α , γ β ) = X ( α i + δ i ) log (1 + t i /β i ) + X δ i log( β i /α i ) (6) where again δ i is zero for censored and o ne for uncensored observ ations, and wh ere α i = α exp  γ α ( x i − 60)  and β i = β exp  γ α ( x i − 60)  are the comp osition-sp ecific parameters go verning the dur atio n t i of the i th eruption, with silica con tent x i . Th e MLEs can now b e found with a four -d imensional searc h o ver α, β ∈ R + and γ α , γ β ∈ R or b y a three-dimensional search o v er β , γ α , γ β using th e conditional MLE ˆ α ( β , γ α , γ β ) = n 1 . X i e γ α ( x i − 60) log  1 + t i e − γ β ( x i − 60) /β  where n 1 = | I 1 | = P δ i is the num b er of uncensored observ ations. 13 A.4 Ob jectiv e Ba yesian Estimates and F orecasts Ob jectiv e Bay esian indep end en t reference p rior distr ibutions ( Bernar do 1979, Ber ge r e t al. 2009) α ∼ 1 /α and β ∼ 1 /β were us ed f or the mo del parameters, b oth improp er scale-in v ariant distri- butions on R + . Post erior distributions are prop er and hav e fin ite means and v ariances so long as n 1 ≥ 3 (see b elo w). Results w ere insens itiv e to these c hoices. Ba y esian p osterior estimates of parameter v alues and d uration forecasts based on this prior and the negativ e log likeli ho o ds of Eqn s (5, 6) are ev aluated using the Metrop olis-Ha stings v ariation ( Hastings 1970, Metr op olis et al. 1953) of th e Marko v chain Mon te Carlo simulat ion-based com- putational metho d ( Besag et al. 1995, Gelfand and Smith 1990, Gilks et al. 1996, Tierney 199 4). After an initial burn-in p erio d of 10 4 steps a fur ther 10 6 MCMC iterations were p erformed. MCMC samples we re th inned at r ate 1 / 10 3 to eliminate measurable auto correlation, lea ving a sample of 10 3 essen tially in d ep enden t and id en tically-distributed (iid) ob s erv ations { ( α j , β j , γ α j , γ β j ) } from the joint p osterior distribu tion. Sample quant iles and momen ts from this sample are used to find in terv al and p oint estimates for the parameters, while ev aluating Eq n (3) and Eqns(4a, 4b) along the MCMC s ample give s the forecasts u sed to generate Figure (4) and the b ottom row of T able 2. Pro of of Posterior P r opriet y Fix a, b, c, d ≥ 0 and let α ∼ Ga ( a, b ) and β ∼ Ga ( c, d ) hav e indep enden t Gamma prior distribu tions. Set X β := P log(1 + t i /β ) and Y β := P δ i log( β + t i ), eac h a fu n ctio n of β . If the data in clude at least one u ncensored ob s erv ation t i > 0 ( i.e. , n 1 ≥ 1) then X β > 0 and Y β > n 1 log β . As β → 0, X β ≍ n log(1 /β ) → ∞ and Y β → Y 0 := P δ i log t i ∈ R ; as β → ∞ , X β → 0 and Y β ≍ n 1 log β → ∞ . The joint p osterior p robabilit y d istribution for ( α, β ) has a den sit y pr op ortional to π ( α, β | t ) ∝ α a − 1 e − bα β c − 1 e − dβ exp n − X ( α + δ i ) log (1 + t i /β ) − n 1 log β + n 1 log α o = α a + n 1 − 1 exp {− α [ b + X β ] } β c − 1 exp {− dβ − Y β } In tegrating wrt α o ver R + giv es the marginal p osterior for β : π ( β | t ) ∝ [ b + X β ] − a − n 1 β c − 1 exp {− dβ − Y β } Asymptotically this is π ( β | t ) ≍ β c − 1 / log (1 /β ) a + n 1 as β → 0 and π ( β | t ) ≍ β c − n 1 − 1 e − dβ as β → ∞ so p osterior p ropriet y will follo w if either c > 0 or a + n 1 > 1, for integrabilit y near β ≈ 0, and either d > 0 or n 1 > c , for in tegrabilit y near β ≈ ∞ . Po sterior means and v ariances are fi nite if in addition d > 0 or n 1 > c + 2. In our app lica tion n 1 = 163 ≥ 3, so p osteriors are p rop er and hav e finite means and v ariances eve n f or our reference priors with a = b = c = d = 0. Ac kno wledgmen ts All d ata used for this analysis were tak en from the DomeHaz database ( Ogburn et al. 2012). T h e authors would lik e to thank R. S. J . Sparks and W. P . Aspin all for helpful conv ersations, the staff of the Mon tserrat V olcano Ob serv atory for hospitalit y and encour agement, and t wo anon ymous review ers for helpful suggestions. This wo r k was sup p orted in part b y US National Science F oun - dation grants DMS–122 8317, DMS–1228 217, and EAR–08095 43. 14 Duration Start V o lcano Duration Start V o lcano (years) Y ear Name (years) Y ear Name 5.0 13 10 OKA T AINA 16.2 199 8 IBU ∗ 5.4 19 70 KARANGET ANG [AP I SIAU] 1 8.5 19 13 COLIMA 5.4 18 70 CEBORUCO, V OLCAN 19.7 199 5 SOUFRI ` ERE HILLS ∗ 5.4 19 91 SOPUT AN 23.0 197 2 BAGA NA 5.4 19 44 SHIVELUCH 23.7 199 1 KARANGET ANG [API SIAU] ∗ 5.5 19 51 LAMINGTON 27.0 188 3 BOGOSLOF 6.0 18 72 SINARKA 27.1 179 6 BOGOSLOF 6.6 19 80 ST. HELE NS 27.6 197 3 LANGILA 7.1 19 94 ETNA 34.6 198 0 SHIVELUCH ∗ 8.6 19 84 LASCAR 40.0 186 9 COLIMA 8.7 18 97 DONA JUANA 42.5 196 8 ARENAL 10.2 200 5 POPO CA TE PETL ∗ 45.0 189 0 VICTOR Y 10.3 200 4 REVENT ADOR ∗ 57.8 195 7 COLIMA ∗ 11.3 200 0 SOPUT AN 59.4 195 5 BEZYMIANNY ∗ 12.4 197 0 KAR YMSKY 68.4 194 6 SEMERU ∗ 13.0 197 3 CHILLAN, NEV ADOS DE 78.8 193 4 SANGA Y 13.2 200 2 FUEGO ∗ 92.7 192 2 SANT A MARIA [SANTIAGUITO] ∗ 13.3 200 1 KAR YMSKY ∗ 187.7 17 28 SANGA Y 15.4 199 9 MA YON 246.6 17 68 MERAPI ∗ T able 3: Durations (in yea rs) of all eruptions lasting five yea r s or more in DomeHaz ( Ogburn et al. 2012). Th ose ongoing at pu b licat ion date are mark ed “*”. 15 References Ak aike , H. (1974) , A new lo ok at the statistical mo del identificati on, IEEE T r ansactions on Auto- matic Contr ol , 19 (6), 716–723, doi:10.110 9/T AC.19 74.1100705. Andrews, B. J. (20 14), Magmatic storage conditions, decompression rate, and incipient caldera collapse of th e 1902 eruption of San ta Mar ´ ıa Volcano, Guatemala, Journal of V olc anolo gy and Ge othermal R ese ar ch , 282 , 201–114, doi:10.101 6/j.jv olgeores.2014.06.009. Barmin, A., O. E. Melnik, and R. S. J. Sparks (2002 ), Periodic b eha viour in lav a dome erup tions, Earth and Planetary Sci enc e L etters , 199 (1–2), 173–1 84, doi:10.1016 /S0012-821X(02)00 55 7- 5 . Bebbington, M. S. (2007), Identifying v olcanic regimes using h idden Mark ov mo dels, Ge ophysic al Journal International , 171 (2), 921–9 42, doi:10.1111 /j.1365- 246X.2007.03559.x. Berger, J. O., J. M. Bernardo, and D. Sun (2009 ), Th e formal defin ition of reference priors, Ann. Stat. , 37 (2), 905–938 , d oi:1 0.1214 /07- A OS587. Bernardo, J. M. (1979), Referen ce p osterior distributions for Ba yesia n inferen ce (with discuss ion), J. R oy. Stat. So c. B, 41 (2), 113–147 . Besag, J., P . J. Green, D. Higdon, and K. Mengersen (1995), Ba yesia n computation and sto c h astic systems (with discussion), Stat. Sci. , 10 (1), 3–66, d oi:10.1214/ ss/11770 10123. Casella, G., an d R. L. Berger (1990), Statistic al Infer enc e , Duxbur y Pr ess, Belmont , CA. Coles, S. G. (2001), An Intr o duction to Statistic al Mo deling of Extr eme V alues , Springer-V erlag, New Y ork, NY. Costa, A., O. E. Melnik, R. S. J. Sp arks, and B. V oigh t (20 07), Con trol of magma fl o w in dyk es on cyclic la v a dome extrusion, Ge ophysic al R ese ar ch L etters , 34 (L02303), 1–5, doi: 10.102 9/2006GL027466. Crosw eller, H. S ., B. A rora, S . K. Bro wn , E. Cottrell, N. I. Deligne, N. O. Guerrero, L. Hobb s, K. Kiy osugi, S. C . Loughlin, J. Lo w n des, M. Na yem bil, L. Sieb ert, R. S . J . Sparks, S. T ak arada, and E. V enzk e (2012), Global database on large magnitude explosiv e vol canic erup - tions (LaMEVE), Journal of Applie d V olc anolo gy , 1 (1), doi:10 .1186/2191 - 5040- 1- 4, on-line at http://d x.doi.org/10. 1186/2191- 5040- 1- 4 . Denlinger, R. P ., and R. P . Hoblitt (1999), Cyclic eru ptiv e b eha vior of silicic volc ano es, Ge olo gy , 27 (5), 459–46 2, doi:10.1130/ 0091-7613. Gelfand, A. E., and A. F. M. Smith (1990 ), Sampling-based app roac h es to calculating marginal densities, J. Am. Stat. Asso c. , 85 (410) , 398–409, doi:10.230 7/228 9776. Gilks, W. R., S. Ric hard son, and D. J. Sp iege lhalter (Ed s.) (1996), Markov Chain Monte Carlo i n Pr actic e , Chapman & Hall, New Y ork, NY. Global V olcanism Program (2013), V olc ano es of the World , v ol. 4.4.1, Venzk e, E. (ed.). Smithsonian Institution, d oi:1 0.5479/si.GVP .V O T W4-2 013. 16 Gunn, L. S., S . Blak e, M. C . Jones, an d H. Rymer (201 4), F orecasting the d uration of v olcanic eruptions: an empirical probabilistic mo del, Bul letin of V olc anolo g y , 76 (780), 1–18, doi:10.10 07/ s00445 - 013- 0780 - 8. Hastings, W. K. (197 0), Monte Carlo sampling metho ds u sing Mark ov chains and their applications, Bometrika , 57 (1), 97–109, doi:10.23 07/2334940. Mastin, L. G., M. Guffan ti, R. Servr anc kx, P . W ebley , S. Barsotti, K. Dean, A. Duran t, J. W. Ewert, A. Neri, W. I. Rose, D. S c hneider, L. S iebert, B. Stu nder, G. Sw anson, A. T upp er, A. V olen tik, and C. F. W a ythomas (2009), A m u ltidisciplinary effort to assign realistic source p aramete rs to mo dels of volca n ic ash-cloud trans p ort and disp ersion during eruptions, Journal of V olc anolo gy and Ge othermal R ese ar ch , 186 (1-2), 10–21 , doi:10.101 6/j.jv olgeores.2009.01.008. Melnik, O. E., and R. S. J. S p arks (2005), Cont rols on conduit magma fl o w dynamics du ring la v a dome building eruptions, Journal o f Ge ophysic al R ese ar ch: Solid Earth , 110 (B2), 1–21, doi:10.10 29/20 04JB003183. Metrop olis, N. C., A. W. Rosenbluth, M. N. Rosenbluth, A. H. T eller, and E. T eller (1953 ), Equ a- tions of state calculatio ns b y fast compu ting mac h ines, Journal of Chemic al Physics , 21 (6), 1087– 1091, doi:10.1063 /1.169 9114. Newhall, C. G., and W. G. Melson (1983), Explosive activit y asso ciate d with the growth of v olcanic domes, Journal of V olc anolo gy an d Ge othermal R ese ar ch , 17 (1-4), 111–131 , doi: 10.101 6/0377- 0273(83)90064- 1. Ogburn, S. E., S . C. Loughlin, and E. S. Calder (2012), DomeHaz: Dome-forming eruptions datab ase , on-line at https:/ /vhub.org/re sources/1742 . Ogburn, S. E., S. C. Loughlin, and E. S . Calder (2015), The asso ciati on of la v a d ome-g ro wth with ma jor explosive activit y (VEI ≥ 4): DomeHaz, a global dataset, Bu l letin of V olc anolo gy , 77 (40), doi:10.10 07/s00 445- 015- 0919- x. P allister, J. S. (2015), p ersonal communicatio n. Ro dado, A., M. S. Bebbington, A. Noble, S. Cronin, and G. Jolly (2011), On selection of analog v olcano es, M ath ematic al Ge oscienc es , 43 (5), 505– 519, doi:10.1007 /s11004- 011- 9345- 6. Sc hw arz, G. E. (1978 ), Estimating the dimension of a mo del, Ann. Stat. , 6 (2), 46 1–464, doi: 10.121 4/aos/1176344136. Simkin, T. (1993), T errestrial v olcanism in s pace and time, Annual R eviews of Earth and Planetary Scienc e , 21 , 427–452, d oi:1 0.1146 /annurev.ea.21 .050193.002235. Simkin, T., and L. Sieb ert (1994), V olc ano es of the World: a R e g iona l Dir e ctory, Gazette er, and Chr onolo gy of V olc anism D uring the L ast 10,000 Y e ars , 2nd ed ., Geoscience Press, T ucson, AZ. Simkin, T., L. Sieb ert, L. McClelland, and D. Bridge (19 81), V olc ano es of the World: a R e gional D i- r e ctory, Gazette er, and Chr onolo gy of V olc anism D uring the L ast 10,000 Y e ars , 1st ed., Hutchison Ross, Strou d sburg, P A, USA. 17 Sparks, R. S. J ., and W. P . Aspinall (2004), V olcanic activit y: F ron tiers and c hallenges in fore- casting, pr edictio n and r isk assessment, in The State of the Planet: F r ontiers and Chal lenges in Ge ophysics , IUGG Mono gr aph , v ol. 19, edited by R. S. J. Sparks and C. J. Ha wke sw orth , p p. 359–3 73, A GU, New Y ork, NY, d oi:10.1029/ 150GM2 8, geoph ysical Monograph 150. Tierney , L. (1994 ), Marko v chai ns for exploring p osterior distributions (with discussion), Ann. Stat. , 22 (4), 1701–17 62, d oi:1 0.1214 /aos/1176325750. V en zke, E., R. W. W underman, L. McClelland, T. Simkin , J. F. Luhr, L. Sieb ert, G. Ma yb erry , and S. K. Sennert (Eds.) (2002), Glob al V olc anism, 1968 to the P r esent , Sm ith s onian Institute, Global V olcanism Program, global V olcanism Pr ogram Dig ital Information Series, GVP-4, on-line at http: //www.volcano .si.edu/reports/ . 18