Alvaro Aravena

Olivier Roche

Influence of the topography of stratovolcanoes on the propagation and channelization of dense pyroclastic density currents analyzed through numerical simulations

). The relevance of this classification lies on the possibility of defining volcanic analogues (defined here as volcanoes that share a suite of topographic characteristics and may be considered comparable to a certain extent) and identifying the main processes that may affect PDC propagation in specific topographic contexts. These aspects are useful for studying poorly documented volcanic edifices and for volcanic hazard assessment. Additionally, we compare this classification with published morphometric characteristics of volcanoes, showing that morphometric parameters such as mean slope of the low flank, irregularity index, ratio of volcano height and basal width, and ratio of crater width and basal width, are useful variables for recognizing the groups we defined. These parameters can be used as rough indicators of the expected interaction patterns between the topography of a given volcano and dense PDCs.

Introduction

The propagation dynamics of pyroclastic density currents (PDCs) and the resulting run-out distance is controlled by the eruption source parameters (e.g., mass flow rate, volume, concentration of solid particles, temperature, grain size distribution, initial velocity; Esposti [START_REF] Ongaro | Transient 3D numerical simulations of column collapse and pyroclastic density current scenarios at Vesuvius[END_REF][START_REF] Roche | Discharge rate of explosive volcanic eruption controls runout distance of pyroclastic density currents[END_REF][START_REF] Shimizu | The run-out distance of large-scale pyroclastic density currents: A two-layer depth-averaged model[END_REF] and by the topography that the pyroclastic mixture encounters during propagation and emplacement, which is the result of the complex interplay of constructive and destructive geological processes (e.g., [START_REF] Grosse | Morphometry and evolution of arc volcanoes[END_REF][START_REF] Germa | Construction and destruction of Mont Pelée volcano: Volumes and rates constrained from a geomorphological model of evolution[END_REF][START_REF] Castruccio | The Influence of Plumbing System Structure on Volcano Dimensions and Topography[END_REF]. In fact, topographic features frequently observed in volcanic areas such as craters, calderas and high-slope radial valleys control the propagation of the dense basal part of PDCs (e.g., [START_REF] Douillet | Sedimentology and geomorphology of the deposits from the August 2006 pyroclastic density currents at Tungurahua volcano, Ecuador[END_REF][START_REF] Martí | Topographical controls on small-volume pyroclastic flows[END_REF][START_REF] Doronzo | Energy facies: A global view of pyroclastic currents from vent to deposit[END_REF]), which we call dense PDCs hereafter. Interaction with topography can affect flow rheology through a series of complex mechanisms, including excess pore pressure due to reduced basal/wall friction [START_REF] Breard | The Basal Friction Coefficient of Granular Flows With and Without Excess Pore Pressure: Implications for Pyroclastic Density Currents, Water-Rich Debris Flows, and Rock and Submarine Avalanches[END_REF]) and bulking processes [START_REF] Bernard | Pyroclastic flow erosion and bulking processes: comparing field-based vs. modeling results at Tungurahua volcano, Ecuador[END_REF]). Thus, volcano topography influences critically the hazard zonation of volcanoes [START_REF] Itoh | Hazard estimation of the possible pyroclastic flow disasters using numerical simulation related to the 1994 activity at Merapi Volcano[END_REF][START_REF] Macías | Hazard map of El Chichón volcano, Chiapas, México: Constraints posed by eruptive history and computer simulations[END_REF][START_REF] Charbonnier | Numerical simulations of block-and-ash flows using the Titan2D flow model: Examples from the 2006 eruption of Merapi Volcano, Java, Indonesia[END_REF][START_REF] Neri | Quantifying volcanic hazard at Campi Flegrei caldera (Italy) with uncertainty assessment: 2. Pyroclastic density current invasion maps[END_REF][START_REF] Charbonnier | New Insights Into the 2070 cal yr BP Pyroclastic Currents at El Misti Volcano (Peru) From Field Investigations, Satellite Imagery and Probabilistic Modeling[END_REF][START_REF] Bevilacqua | Thematic vent opening probability maps and hazard assessment of small-scale pyroclastic density currents in the San Salvador volcanic complex (El Salvador) and Nejapa-Chiltepe volcanic complex (Nicaragua)[END_REF]. For example, earlier works have shown the critical influence of Mt.

Somma and Posillipo Hill on the propagation dynamics of PDCs at Vesuvius (Gurioli et al. 2010) and Campi Flegrei [START_REF] Rossano | Numerical simulation of pyroclastic density currents on Campi Flegrei topography: a tool for statistical hazard estimation[END_REF][START_REF] Neri | Quantifying volcanic hazard at Campi Flegrei caldera (Italy) with uncertainty assessment: 2. Pyroclastic density current invasion maps[END_REF], respectively, as well as the effect of the asymmetric crater configuration of Merapi on PDCs generated by dome collapses [START_REF] Thouret | Toward a revised hazard assessment at Merapi volcano, Central Java[END_REF][START_REF] Procter | Mapping block-and-ash flow hazards based on Titan 2D simulations: a case study from Mt. Taranaki, NZ[END_REF][START_REF] Charbonnier | Evaluation of geophysical mass flow models using the 2006 block-and-ash flows of Merapi Volcano, Java, Indonesia: Towards a short-term hazard assessment tool[END_REF]. Flank collapse scars such as those of Tungurahua or Reventador are also considered as major topographic features in controlling the propagation direction of PDCs [START_REF] Hall | Tungurahua Volcano, Ecuador: structure, eruptive history and hazards[END_REF][START_REF] Pennec | Combining historical and 14C data to assess pyroclastic density current hazards in Baños city near Tungurahua volcano (Ecuador)[END_REF].

Moreover, the channelization of concentrated pyroclastic material allows reducing the energy dissipation rate, permitting the flows to reach larger distances than their non-channelized counterparts, and also enhances thermal insulation and thus promote hot overspills in unconfined areas at valleys bends [START_REF] Hutchison | Generation of Overspill Pyroclastic Density Currents in Sinuous Channels[END_REF]. Different strategies have been adopted for the morphometric characterization of volcanoes (e.g., [START_REF] Pike | Volcanoes on the inner planets -Some preliminary comparisons of gross topography[END_REF][START_REF] Pike | Revised classification of terrestrial volcanoes and catalog of topographic dimensions, with new results on edifice[END_REF][START_REF] Grosse | Morphometry and evolution of arc volcanoes[END_REF][START_REF] Grosse | Systematic morphometric characterization of volcanic edifices using digital elevation models[END_REF][START_REF] Grosse | A global database of composite volcano morphometry[END_REF] and a robust dataset is currently available in the literature [START_REF] Grosse | A global database of composite volcano morphometry[END_REF]. However, although morphometric data have been interpreted in terms of the growth history and evolution of volcanoes, the influence of topographic features on the dispersion of volcanic products, including PDCs, has yet not been addressed systematically.

The complexity and variability of the topography of stratovolcanoes, as well as the incompleteness of the volcanological record, have hampered the development of field databased studies on the effect of stratovolcanoes topography in PDC propagation. Alternatively, in this work we use an approach based on numerical modelling, which allows studying different volcanic systems using a common input dataset and methodology. Different models allow simulating the propagation of PDCs [START_REF] Dufek | Pyroclastic Density Currents: Processes and Models[END_REF]. Simple formulations such as the energy cone model [START_REF] Malin | Computer-Assisted Mapping of Pyroclastic Surges[END_REF][START_REF] Sheridan | Application of computer-assisted mapping to volcanic hazard evaluation of surge eruptions: Vulcano, lipari, and vesuvius[END_REF] do not describe properly the effect of topography and the occurrence of channelization. More complex models such as depthaveraged or multi-phase formulations [START_REF] Ongaro | Transient 3D numerical simulations of column collapse and pyroclastic density current scenarios at Vesuvius[END_REF][START_REF] Charbonnier | Numerical simulations of block-and-ash flows using the Titan2D flow model: Examples from the 2006 eruption of Merapi Volcano, Java, Indonesia[END_REF][START_REF] Procter | Mapping block-and-ash flow hazards based on Titan 2D simulations: a case study from Mt. Taranaki, NZ[END_REF][START_REF] Kelfoun | A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents[END_REF][START_REF] De' Michieli | IMEX_SfloW2D 1.0: a depth-averaged numerical flow model for pyroclastic avalanches[END_REF]) are limited by their computational cost and thus they cannot be applied systematically on a large set of volcanoes.

As a compromise solution, in this work we adopt the branching energy cone model [START_REF] Aravena | Tree-Branching-Based Enhancement of Kinetic Energy Models for Reproducing Channelization Processes of Pyroclastic Density Currents[END_REF], which is a recently developed reformulation of the traditional energy cone model. This model, which suits better for simulating the dense basal part of PDCs irrespective of their source mechanisms [START_REF] Cole | Deposits from dome-collapse and fountain-collapse pyroclastic flows at Soufrière Hills Volcano, Montserrat[END_REF][START_REF] Gueugneau | Investigation of surge-derived pyroclastic flow formation by numerical modelling of the 25 June 1997 dome collapse at Soufrière Hills Volcano[END_REF], allows describing PDC channelization processes with a limited computational cost. The systematic application of the branching energy cone model on a large set of stratovolcanoes allows identification of how the main topographic features of volcanoes (recognizable in a 30-m resolution DEM; e.g., summit crater, decametric or larger valleys, and proximal barriers) are able to affect the propagation of PDCs. Moreover, this permits us to classify stratovolcanoes in terms of the expected interaction pattern between their topography and dense PDCs, and to compare them with published morphometric data [START_REF] Grosse | A global database of composite volcano morphometry[END_REF]). Thus, this approach offers the possibility of identifying the main processes that may affect PDC propagation in specific topographic contexts and recognizing eventual volcanic analogues, defined by [START_REF] Tierz | VOLCANS: an objective, structured and reproducible method for identifying sets of analogue volcanoes[END_REF] as volcanoes that share enough characteristics to be considered comparable to a certain extent. This is particularly useful for volcanic hazard assessment and for studying poorly documented volcanic systems.

We remark that our study is exclusively devoted to the analysis of the potential effect of stratovolcanoes topographies on the propagation of dense PDCs, i.e. with no consideration on the occurrence probability or the expected size of PDCs in the volcanic systems studied. This paper is organized in five sections. We first describe the methods, with emphasis on the type of specific results extracted from each numerical simulation. Then we present the results, including the classification of stratovolcanoes based on the interaction pattern between their topographies and the simulated PDCs. Finally, we describe the comparison of our results with morphometric data of the studied volcanoes, and we present the discussion and concluding remarks.

Methods

Using topographic information derived from the database SRTM 30 m [START_REF] Rabus | The shuttle radar topography mission-a new class of digital elevation models acquired by spaceborne radar[END_REF], we performed different sets of simulations using the branching energy cone model [START_REF] Aravena | Tree-Branching-Based Enhancement of Kinetic Energy Models for Reproducing Channelization Processes of Pyroclastic Density Currents[END_REF][START_REF] Aravena | Calibration strategies of PDC kinetic energy models and their application to the construction of hazard maps[END_REF] considering 50 stratovolcanoes (Table 1 andFigure 1). This model is a reformulation of the traditional energy cone model [START_REF] Malin | Computer-Assisted Mapping of Pyroclastic Surges[END_REF][START_REF] Sheridan | Application of computer-assisted mapping to volcanic hazard evaluation of surge eruptions: Vulcano, lipari, and vesuvius[END_REF]. It allows consideration of flow channelization and thus captures the effect of topography on PDC propagation [START_REF] Aravena | Tree-Branching-Based Enhancement of Kinetic Energy Models for Reproducing Channelization Processes of Pyroclastic Density Currents[END_REF][START_REF] Aravena | Calibration strategies of PDC kinetic energy models and their application to the construction of hazard maps[END_REF][START_REF] Bevilacqua | Thematic vent opening probability maps and hazard assessment of small-scale pyroclastic density currents in the San Salvador volcanic complex (El Salvador) and Nejapa-Chiltepe volcanic complex (Nicaragua)[END_REF]. In this formulation, a root energy cone is complemented with branch energy cones along the directions of preferential channelization. Each branch energy cone is defined considering a collapse height controlled by the residual potential energy computed in its channelization zone. The branch energy cones are organized in a tree-like structure whose construction is stopped when the new energy cones do not add pixels to the resulting inundation area. Note that the inputs of the branching energy cone model are exactly the same as that of the traditional formulation, i.e. initial height of the root energy cone (𝐻 0,0 ), energy cone slope (tan (𝜑)) and location of collapse. Each of the 50 sets of simulations comprises 1,000 runs with variable values of 𝐻 0,0 (from 100 m to 1000 m), tan (𝜑) (from 0.2 to 1.0) and collapse location, which was sampled uniformly within a 500 m-radius circle centred on the summit or crater area of each volcano (Figure 2). The values adopted for 𝐻 0,0 and tan (𝜑) are within ranges expected for dense PDCs sourced from collapsing domes or eruptive columns from low to moderate height (up to a few kilometres, note that the interpretation of H 0,0 as equal to the collapse height may be misleading in PDCs derived from column collapse; [START_REF] Aravena | Calibration strategies of PDC kinetic energy models and their application to the construction of hazard maps[END_REF]. These input ranges allow the simulation of run-out distances from <1 km to a few tens of kilometres, as we show below.

Note that column collapse from greater heights would be dominated by the generation of voluminous dilute PDCs that are better described using other formulations such as the box model (Esposti [START_REF] Ongaro | A fast, calibrated model for pyroclastic density currents kinematics and hazard[END_REF]. We stress that, in the branching energy cone model, the initial collapse of pyroclastic material is described as axisymmetric, and thus this formulation is not able to simulate directional flows. Because the input parameters were not calibrated using the volcanological record of each volcano, we did not analyse the results in terms of the simulated inundation zones but rather in terms of the statistical distribution of model outputs and the relationships among them (in other words, the resulting probability maps of PDC inundation are not considered relevant in terms of hazard evaluation). In particular, for each numerical simulation, we extracted the following parameters from the inundation polygon: [START_REF] Gueugneau | Investigation of surge-derived pyroclastic flow formation by numerical modelling of the 25 June 1997 dome collapse at Soufrière Hills Volcano[END_REF] maximum run-out distance (𝑅 𝑚𝑎𝑥 ), (2) minimum run-out distance (𝑅 𝑚𝑖𝑛 , minimum distance between the source and a point belonging to the inundation area contour), (3) inundation area (IA), ( 4) perimeter (𝑃), and ( 5) solidity (𝑆), the latter defined as the inundation area divided by the area of the smallest convex polygon containing the invasion zone. From these parameters, we also computed (1) IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), (2) 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 , and (3) 𝐶 𝐹 = 2√𝜋 • 𝐼𝐴/𝑃. These parameters, as well as 𝑆, range between zero and one, and their combination allows understanding the degree of channelization of the simulated PDCs. For instance, a perfectly circular inundation area would produce a value of 1 for all these parameters, while the concomitance of channelization zones in different directions would generate values close to 0 for all the described parameters (note that a single well-developed channelization zone would translate into IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 , and 𝐶 𝐹 close to 0 and 𝑆 close to 1). We compared our numerical results with published morphometric information of volcanoes [START_REF] Grosse | A global database of composite volcano morphometry[END_REF], including volcano size parameters, profile shape parameters, plan shape parameters, and slope parameters (see Section 4).

We highlight that the use of a ~30 m-resolution DEM [START_REF] Rabus | The shuttle radar topography mission-a new class of digital elevation models acquired by spaceborne radar[END_REF] does not permit us to consider small-scale channels and therefore represents a limitation of our approach. However, with such a resolution we can apply a common methodology for the complete set of volcanoes.

We stress also that the morphology in the summit zone of some volcanoes, such as Merapi and Sangay, changed significantly during the last decade, which is not considered in the DEMs adopted. Note, however, that the simulations performed for these volcanoes (see Section 3.1)

do not include PDCs that stopped in the summit area and thus the effect of summit topography modifications on the resulting inundation areas is expected to be limited.

Results

In this Section, we use two approaches to address the effect of topography on PDC propagation.

In Section 3.1 we describe the main topographic features that are recognizable from numerical results, while in Section 3.2 we classify the studied stratovolcanoes based on the statistical distributions of IA/(π • 𝑅 𝑚𝑎𝑥 2 ), 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 and 𝑆.

Main topographic features

Our results (Supplementary Material) show that the volcano topography has significant effects on the simulated inundation polygons. Here we describe the main topographic features (TF) of volcanoes whose effects on PDC propagation are clearly recognizable from numerical results:

(a) Steep slopes in proximal zones (TF1).

In some cases, there is a gap in the simulated run-out distances in very proximal domains, or even the absence of simulations with small run-out distances. Some examples are Fuego, Guallatiri, Merapi and Sangay (see Table 1, Supplementary

Material and the case of Tungurahua in Figure 3). This is a consequence of the presence of particularly steep slopes in proximal zones, which inhibit flow stopping near the source. This is confirmed by the comparison of our results with the morphometric parameters presented by [START_REF] Grosse | A global database of composite volcano morphometry[END_REF]. In fact, a two-sample t-test, which allows us evaluating the hypothesis that the morphometric parameters of both sets of volcanoes (i.e. with and without gaps in the simulated run-out distances in very proximal domains) come from independent random samples from normal distributions with equal means and equal but unknown variances, shows that volcanoes with a significant gap in the simulated run-out distances present larger maximum average slopes than the rest of the analysed volcanoes (mean value of 33.2° and standard deviation of 3.6° compared to 24.6° and 4.7°, with p-value much lower than 0.01). Moreover, other morphometric variables for which both sets of volcanoes present significantly different mean values (i.e. with p-values lower than 0.05) include the ratio of height and basal width (𝐻/𝑊 𝐵 ), the ratio of summit width and basal width (𝑊 𝑆 /𝑊 𝐵 ), mean slope angle of the main flank, and summit mean slope angle, among others. We speculate that, for these volcanoes, documented small run-out distance PDCs that stopped on steep slopes were probably limited by their volume, which cannot be taken into account in kinetic energy models (see the analysis for Merapi in [START_REF] Aravena | Calibration strategies of PDC kinetic energy models and their application to the construction of hazard maps[END_REF]).

(b) Summit crater (TF2).

Some volcanoes present a cluster of simulations with particularly small run-out distances, such as San Salvador, Chaitén, and Kelut (see Table 1, Supplementary

Material and the case of Chaitén in Figure 3). This behaviour is related to the presence of a summit crater deep/wide enough to limit the propagation of the smallest PDCs (i.e.

those characterized by low values of 𝐻 0,0 and high values of tan (𝜑)), which remain confined in the summit area. In this case, the comparison with the morphometric parameters of [START_REF] Grosse | A global database of composite volcano morphometry[END_REF] shows that volcanoes exhibiting the above-described effect of the summit crater tend to present smaller values of 𝐻/𝑊 𝐵 (0.09 ± 0.03 versus 0.15 ± 0.04), where 𝐻 is volcano height and 𝑊 𝐵 is basal width, with a p-value lower than 0.01.

(c) Proximal topographic obstacles (TF3).

The expected positive correlation between IA/(π • 𝑅 𝑚𝑎𝑥 2 ) and 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 , for some volcanoes, is partially masked by the presence of a set of simulations with very low values of 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 (less than 0.2) and variable results of IA/(π • 𝑅 𝑚𝑎𝑥 2 ), typically between 0.1 and 0.5 (e.g. Tungurahua, Merapi and Fuego, Table 1 and Figure 3). This reflects the presence of proximal topographic obstacles (e.g. an asymmetrical crater configuration such as those observed at Merapi and Tungurahua) influencing the preferential propagation direction of PDCs during early transport phases. This process may increase significantly the run-out distance because it allows reduction of the early energy dissipation rate and prevents the spreading of pyroclastic material over a larger area (Kubo Hutchison and Dufek 2021). Consistently, the volcanoes exhibiting proximal topographic obstacles tend to present a proximal gap in the simulated run-out distances (cf. TF1; Table 1).

(d) Radial valleys with slope breaks (TF4).

In some cases, the distribution of simulated run-out distances is clearly multimodal (e.g.

Chillán, Peteroa and El Misti; see Table 1, Supplementary Material and the cases of Galeras and Teide in Figure 3). This indicates that one or more valleys control the propagation of PDCs, and these valleys are characterized by one or more zones of slope break that generate a set of peaks in the resulting distribution of run-out distance. All the examples recognized with a clear multimodal distribution of run-out distance (Table 1) present well-developed channelization zones (see Supplementary Material). In fact, the average values of IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) for volcanoes with and without multimodal distributions of run-out distance are 0.22 ± 0.05 and 0.38 ± 0.12, respectively; while the average values of 𝐶 𝐹 are 0.48 ± 0.06 and 0.62 ± 0.10, respectively.

Classification of volcanoes based on the interaction between their topographies and dense PDCs

According to the distributions of IA/(π • 𝑅 𝑚𝑎𝑥 2 ), 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 and 𝑆 (Table 1 andTables S1-S2 in the Supplementary Material), we classified the studied volcanoes in five groups (Figure 3). Note that we define the proximal and distal domains according to the simulated range of run-out distances for each volcano.

(a) Group A: strong channelization in different valleys up to distal domains.

The topography of these volcanoes (e.g. Colima and Peteroa; Table 1 andFigure 3) is able to induce intense channelization through different radial valleys, causing positively skewed distributions ofIA/(π • 𝑅 𝑚𝑎𝑥 2 ) (skewness higher than 0.85) and nearly symmetric to positively skewed distributions of 𝑆 (skewness higher than -0.3). The combined effect of propagation valleys in different directions is also manifested in multimodal distributions of run-out distance (i.e. TF4), with values of run-out distance typically higher than that observed for the other groups. Most of these volcanoes present an inverse relationship between run-out distance and 1, Supplementary Material and Figure 3), which we considered to define this group (see caption of Table 1 for the details; [START_REF] Hartigan | The Dip Test of Unimodality[END_REF].

𝐶
(e) Group E: weak channelization in proximal domains.

The topography of this group of volcanoes (e.g. Kelut and Akagi; see Supplementary Material and Figure 3) is not able to induce efficient channelization in proximal domains, due to the presence of a large crater (e.g. Pinatubo and San Salvador) or the absence of proximal ravines able to control significantly PDC propagation. In fact, most of the volcanoes presenting a cluster of simulations with particularly small run-out distances due to the effect of the summit crater (i.e. TF2, see Section 3.1) are part of Group E (Table 1). The simulated flows able to overcome the proximal domain of limited channelization eventually propagate through radial valleys causing efficient channelization (e.g. Kelut and Pinatubo), but in any case the significant proximal energy dissipation is typically manifested in run-out distances much smaller than those simulated for the other groups.

Comparison with morphometric parameters

In this Section, the groups identified in Section 3.2 are discussed according to the morphometric parameters presented by [START_REF] Grosse | A global database of composite volcano morphometry[END_REF], allowing to recognize the main features of volcanic edifices that determine the groups they belong.

Our results indicate that the volcanoes able to induce intense channelization through different valleys up to distal domains (i.e. Group A) present high values of low flank mean slope angle (17°-25°) and relatively high outline irregularity indexes (>1.22; Figure 4a). On the other hand, high values of low flank mean slope angle in concomitance with small outline irregularity indexes are typically related to type C volcanoes (i.e. intense channelization near the source and moderate distal channelization; Fig. 4a). Groups A and C overlap in the plots of the morphometric parameters 𝐻/𝑊 𝐵 and 𝑊 𝑆 /𝑊 𝐵 as functions of the low flank mean slope angle (Fig. 4c-d, where 𝑊 𝑆 is the summit width), and they partially overlap when the average irregularity index is considered (Fig. 4b). Instead, volcanoes with the potential to induce intense reduced PDC run-out distance. Finally, Group B represents a sort of intermediate member between the four categories described above (Figure 4).

Discussion and concluding remarks

In this study, we have shown that the systematic application of the branching energy cone model to avoid subjective considerations, we defined specific numerical thresholds to set the different groups (see caption of Table 1) from the statistical distributions of different parameters extracted from the simulated inundation polygons (see Table 1 and Tables S1-S2 in the Supplementary Material). Importantly, these groups permit us to identify the expected topographical effect on the propagation of PDCs, which is useful for defining hazard assessment strategies, for studying poorly documented volcanoes, and eventually for defining volcanic analogues.

We have shown that some of the morphometric parameters defined by [START_REF] Grosse | A global database of composite volcano morphometry[END_REF] (in particular, mean slope of the low flank, outline irregularity index, average irregularity index, ratio of volcano height and basal width, and ratio of crater width and basal width) can be used to recognize the five groups we defined. We recall that our study, which is based on the application of the branching energy cone model, addresses exclusively the influence of the topography of stratovolcanoes on the propagation dynamics of dense PDCs (irrespective of their origin, from dome/column collapse-derived to surge-derived pyroclastic flows; [START_REF] Cole | Deposits from dome-collapse and fountain-collapse pyroclastic flows at Soufrière Hills Volcano, Montserrat[END_REF][START_REF] Druitt | Small-volume, highly mobile pyroclastic flows formed by rapid sedimentation from pyroclastic surges at Soufrière Hills Volcano, Montserrat: an important volcanic hazard[END_REF][START_REF] Kelfoun | Suitability of simple rheological laws for the numerical simulation of dense pyroclastic flows and longrunout volcanic avalanches[END_REF][START_REF] Gueugneau | Investigation of surge-derived pyroclastic flow formation by numerical modelling of the 25 June 1997 dome collapse at Soufrière Hills Volcano[END_REF]), while we have not considered the probability or expected volume of PDCs for the stratovolcanoes studied, nor the possible presence of structural specificities able to control vent position. An additional, relevant process that should be considered in the study of dense PDCs is the possible detachment of an upper, dilute portion of the PDC [START_REF] Druitt | Small-volume, highly mobile pyroclastic flows formed by rapid sedimentation from pyroclastic surges at Soufrière Hills Volcano, Montserrat: an important volcanic hazard[END_REF][START_REF] Jenkins | The Merapi 2010 eruption: An interdisciplinary impact assessment methodology for studying pyroclastic density current dynamics[END_REF][START_REF] Wibowo | Emplacement Temperature of the Overbank and Dilute-Detached Pyroclastic Density Currents of Merapi 5 November 2010 Events using Reflectance Analysis of Associated Charcoal[END_REF], able to propagate independently from the dense basal part. Taking into account these volcanological considerations (e.g. expected magnitude, uncertainty in vent position, eruption mechanism), as well as using DEMs with finer resolution, is in fact required for refining the definition of the volcanic analogues presented here and for the development of studies devoted to volcanic hazard assessment.

Table 1. Stratovolcanoes considered in this work and main characteristics of numerical results.

Volcano Location Main topographic features clearly recognizable from numerical results

Group of volcanoes

1 TF1 TF2 TF3 TF4 Akagi (A) Japan X E Asakusa (Ask) 2 Japan X X E Asama (Asm) Japan X E Bandaisai (B) Japan X E Calbuco (Ca) Chile X X E Ceboruco (Ce) Mexico X E Chaiten (Cht) 2 Chile X E Chichon, El (Chc) 2 Mexico X E Chillán, Nevados de (NCh) Chile X A Chimborazo (Chm) Ecuador X X X E Chokai (Chk) Japan X B Colima, Nevado de (NCo) Mexico X X A Cotopaxi (Co) Ecuador X X D Fuego (F) Guatemala X X X C Galeras (Ga) Colombia X A Guallatiri (Gu) Chile X D Haku, Mount (H) Japan X A Kelut (K) Indonesia X E Lascar (L) Chile C Machin, Cerro (CM) Colombia X X C Mayon (Ma) Indonesia X X D Merapi (Mp) Indonesia X X D Meru (Mr) Tanzania X X B Misti, El (EM) Peru X X X C Momotombo (Mo) Nicaragua X X D Ngauruhoe (N) 2 New Zealand X X D Orizaba, Pico de (O) Mexico X X X A Peteroa (Pe) Chile X A Pinatubo (Pi) Philippines X E Quizapu (Q) Chile X X X A Reventador (Re) Ecuador X X B Ruapehu (Ru) New Zealand X X X E Ruiz, Nevado del (Rz) Colombia X A Sangay (Sg) Ecuador X X C San Miguel (SMg) El Salvador X X D San Salvador (SS) El Salvador X E Santa María (SM) Guatemala X X C Semeru (Se) Indonesia X X C Sinabung (Si) Indonesia X X D Socompa (So) Chile X X B Soufrière, La (SG) Guadeloupe E Soufrière Hills (SHi) Montserrat E Spurr (Sp) 2 USA X X X A St. Helens (SHe) USA X X B Taranaki (Ta) New Zealand X X D Teide (Te) Spain X B Tolima (To) Colombia X X A Tungurahua (Tn) Ecuador X X X C Tutupaca (Tt) Ecuador X E Vesuvius (V)
Italy E 1 : Classification based on the distributions of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), 𝑆 and 𝐶 𝐹 . The conditions were tested in the following order (see Tables S1 andS2 in the Supplementary Material):

Group A: skewness of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) higher than 0.85, 95 th percentile of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) lower than 0.7, and skewness of 𝑆 higher than -0.3.

Group B: skewness of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) higher than 0.85, 95 th percentile of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) lower than 0.7, and skewness of 𝑆 lower than -0.3.

Group C: 95 th percentile of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) lower than 0.6.

Group D: at least one of the distributions of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) or 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 is not unimodal. This was tested by computing the Hartigan's dip statistic for unimodality [START_REF] Hartigan | The Dip Test of Unimodality[END_REF]. When the value of dip is less than 0.035, we consider that the distribution is clearly multimodal.

Group E: clear unimodal distributions of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) and 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 .

2 : Not included in the analysed dataset of [START_REF] Grosse | A global database of composite volcano morphometry[END_REF].

TF1: efficient PDC propagation in proximal zones.

TF2: summit crater.

TF3: proximal topographic obstacles.

TF4: radial valleys with slope breaks. We exclude bimodal distributions of run-out distance when one of the peaks is related to the summit crater effect. Supplementary Tables Table S1 5 Dip statistic, derived from the application of the Hartigan's test for unimodality [START_REF] Hartigan | The Dip Test of Unimodality[END_REF]. Values greater than 0.035 imply that the distribution is clearly multimodal.

  𝐹 over almost the entire range of run-out distances (e.g. Peteroa and Chillán, with respective values of 𝐶 𝐹 as small as ~0.25 and ~0.3 for high values of run-out distance), while the associated values of IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), typically lower than 0.4, tend to be poorly correlated with run-out distance (see Supplementary Material and the case of Galeras in Figure3). The regular decreasing trend of 𝐶 𝐹 with run-out distance and the resulting inundation maps (see Supplementary Material) suggest that long run-out distance, channelized flows necessarily involve the presence of proximal catchments of pyroclastic material (i.e. the valleys heads) and long, uninterrupted ravines able to reduce efficiently the rate of energy dissipation during a significant fraction of the PDC propagation. (b) Group B: intense channelization through a single dominant valley up to distal domains. These volcanic systems (e.g. Teide, Reventador and Mt. St. Helens) present positively skewed distributions of IA/(π • 𝑅 𝑚𝑎𝑥 2 ) (skewness higher than 0.85) and negatively skewed distributions of 𝑆 (skewness lower than -0.3). While the low values of IA/(π • 𝑅 𝑚𝑎𝑥 2 ) are associated with intense channelization, their concomitance with high values of 𝑆 is typically related to the presence of only one dominant channelization valley, as observed in the resulting inundation maps (see Figure 3 and Supplementary Material). Preferential channelization directions are caused by asymmetric crater configurations and/or proximal topographic obstacles (i.e. TF3). (c) Group C: intense channelization near the source and moderate distal channelization. These volcanoes (e.g. Fuego and El Misti; see Figure 3 and Supplementary Material) present well-defined proximal ravines producing intense channelization (95 th percentile of 𝐼𝐴/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ) lower than 0.6), often with proximal topographic obstacles (i.e. TF3). This topography causes frequently a clear preferential propagation direction and hinders the simulation of small run-out distance PDCs (i.e. TF1) due to the inefficient energy dissipation during early propagation phases. The combined effect of several radial valleys gives rise to a poor dependency between run-out distance and channelization efficiency in proximal domains. At longer distances from the source, channelization decreases moderately, being poorly correlated with run-out distance. (d) Group D: potentially intense channelization only near the source. These volcanoes (e.g. Sinabung and Mayon) are able to induce flow channelization only in proximal domains, while the presence of flat topographies downstream reduces flow channelization. Note that DEM resolution limitations may accentuate the reduction of channelization efficiency in case of relatively narrow valleys. The combination of wellchannelized flows with small run-out distance and poorly channelized flows with long run-out distance gives rise to bimodal distributions of the parameters describing channelization efficiency (see Table

  channelization only near the source (i.e. Group D) are typically related to low values of irregularity index (outline and average) and of low flank mean slope. These volcanoes also present high values of 𝐻/𝑊 𝐵 and low values of 𝑊 𝑆 /𝑊 𝐵 . These characteristics are consistent with the presence of flat areas in the volcano surroundings, which inhibit channelization in distal domains. On the other hand, Group E volcanoes (i.e. weak channelization in proximal domains) present low values of 𝐻/𝑊 𝐵 and of low flank mean slope angle, and high values of irregularity index (outline and average) and 𝑊 𝑆 /𝑊 𝐵 . These characteristics are consistent with the presence of a relatively extended summit zone where the flow propagates radially (i.e. in absence of channelization zones), resulting in efficient, early energy dissipation and thus

  on a large set of stratovolcanoes allows us to recognize the potential effect of different topographic features of volcanoes on PDC propagation (e.g. steep proximal slopes >~30°, summit crater, topographic obstacles, and radial valleys with slope breaks). Despite possible limitations due to the use of 30-m resolution DEMs, we have shown that these topographic features critically affect the hazard zonation of PDCs and related parameters such as the runout distance and inundation area. Note that[START_REF] Doronzo | Energy facies: A global view of pyroclastic currents from vent to deposit[END_REF] discussed also the interaction between PDCs and volcano topographic features, which were defined in four categories: open topography, channelled topography, topographic barrier and steep slope. Interestingly, our simulations show that volcanic topographies are frequently able to induce under-and overrepresentation of specific run-out distances giving rise to multimodal distributions of this parameter (i.e. TF4), which is due to the presence of significant slope breaks along the channelization valleys. The latter translates into frequent flow stopping in specific zones and improbable flow stopping in other sectors, and should not be interpreted necessarily as the result of multimodal distributions of eruption source parameters caused by the concomitance of different collapse/eruption mechanisms. The relationships between run-out distance and parameters describing the properties of the simulated inundation areas (i.e. IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 , 𝐶 𝐹 and 𝑆) indicate that the channelization efficiency is strongly influenced by volcano topography and PDC volume. While most of the volcanoes are able to induce strong PDC channelization in proximal areas (typically, <~5 km), strong channelization at larger distances from the source (typically, >10-20 km) is possible for stratovolcanoes with steep flanks, with long, uninterrupted radial valleys whose heads (i.e. the zone from which pyroclastic flows can become channelized) are located near the vent, being able to reduce efficiently energy dissipation during a significant portion of the PDC propagation. We defined five groups of stratovolcanoes in terms of the mode of interaction between their topographies and dense PDCs: (1) intense channelization through different valleys up to distal domains; (2) intense channelization through a single dominant valley up to distal domains; (3) intense channelization near the source and moderate distal channelization; (4) potentially intense channelization only near the source; and (5) weak channelization in proximal domains, manifested in efficient early energy dissipation. In order

Figure 1 .

 1 Figure 1. Location of the stratovolcanoes considered in this study. SeeTable 1 for abbreviations.

Figure 2 .

 2 Figure 2. Illustrative example of the input parameters of the simulations performed on a specific volcano (Fuego volcano, Guatemala). (a) Collapse height (sampled uniformly from 100 m to 1000 m). (b) Energy cone slope (sampled uniformly from 0.2 to 1.0). (c) Collapse position, sampled uniformly within a 500 m-radius circle centred on the summit or crater area of the volcano.

Figure 3 .

 3 Figure 3. Illustrative examples of the output parameters for the different groups of volcanoes recognized in this work. From top to bottom: map showing the fraction of simulations that reach each pixel of the map (dark grey zones indicate pixels inundated by most of the simulations, while light grey zones are associated with low inundation probabilities), histogram of 𝑅 𝑚𝑎𝑥 , 𝐶 𝐹 as a function of 𝑅 𝑚𝑎𝑥 , histogram of IA/(𝜋 • 𝑅 𝑚𝑎𝑥 2 ), histogram of 𝑅 𝑚𝑖𝑛 /𝑅 𝑚𝑎𝑥 , and histogram of solidity (𝑆).

Figure 4 .

 4 Figure 4. Relationship between the morphometric parameters of Grosse et al. (2014) with colours indicating the different groups (A-E) of volcanoes recognized in this study. (a) Outline irregularity index versus low flank mean slope. (b) Average irregularity index versus low flank mean slope. (c) Ratio of height and basal width versus low flank mean slope. (d) Ratio of summit width and basal width versus low flank mean slope. Note that the irregularity indexes quantify the irregularity or complexity of the elevation isolines (Grosse et al. 2014).

.

  Main statistical parameters of the results describing PDC channelization for each stratovolcano (1/2).

	Volcano			𝐈𝐀/(𝝅 • 𝑹 𝒎𝒂𝒙	𝟐 )					𝑹 𝒎𝒊𝒏 /𝑹 𝒎𝒂𝒙
		P5 1	Mean	P50 2	P95 3	Sk 4	HD 5	P5 1	Mean	P50 2	P95 3	Sk 4	HD 5
	Akagi	0.07	0.35	0.33	0.69	0.24	0.01	0.05	0.36	0.35	0.67	0.03	0.01
	Asakusa	0.09	0.33	0.32	0.66	0.42	0.01	0.05	0.30	0.30	0.63	0.29	0.02
	Asama	0.07	0.30	0.26	0.61	0.47	0.01	0.04	0.26	0.26	0.51	0.26	0.01
	Bandaisai	0.09	0.40	0.41	0.70	-0.05	0.02	0.04	0.33	0.34	0.61	-0.09	0.03
	Calbuco	0.05	0.30	0.24	0.72	0.72	0.01	0.04	0.29	0.23	0.64	0.49	0.02
	Ceboruco	0.09	0.44	0.47	0.76	-0.11	0.03	0.10	0.43	0.44	0.73	-0.19	0.01
	Chaiten	0.16	0.43	0.42	0.77	0.31	0.01	0.15	0.40	0.36	0.75	0.61	0.01
	Chichon, El	0.08	0.42	0.44	0.71	-0.03	0.01	0.13	0.44	0.46	0.69	-0.22	0.01
	Chillán, Nevados de	0.04	0.21	0.18	0.51	0.98	0.01	0.03	0.22	0.23	0.50	0.50	0.03
	Chimborazo	0.07	0.33	0.31	0.67	0.33	0.01	0.02	0.29	0.22	0.69	0.45	0.02
	Chokai	0.06	0.25	0.21	0.67	1.19	0.01	0.03	0.27	0.26	0.66	0.66	0.03
	Colima, Nevado de	0.03	0.19	0.12	0.54	1.24	0.01	0.01	0.18	0.15	0.51	1.05	0.01
	Cotopaxi	0.07	0.38	0.35	0.75	0.17	0.05	0.03	0.37	0.38	0.74	0.03	0.09
	Fuego	0.09	0.27	0.24	0.54	0.48	0.02	0.01	0.13	0.13	0.27	0.40	0.04
	Galeras	0.03	0.14	0.12	0.32	3.01	0.01	0.02	0.16	0.15	0.36	1.42	0.02
	Guallatiri	0.07	0.37	0.36	0.69	0.06	0.05	0.03	0.32	0.30	0.68	0.27	0.01
	Haku, Mount	0.05	0.19	0.16	0.40	1.34	0.01	0.03	0.20	0.19	0.41	0.60	0.01
	Kelut	0.07	0.40	0.40	0.73	0.10	0.01	0.08	0.37	0.37	0.68	0.15	0.01
	Lascar	0.07	0.29	0.25	0.57	0.65	0.01	0.05	0.25	0.25	0.50	0.58	0.01
	Machin, Cerro	0.08	0.26	0.24	0.50	0.83	0.01	0.05	0.28	0.28	0.51	0.17	0.01
	Mayon	0.11	0.56	0.58	0.93	-0.14	0.08	0.02	0.48	0.53	0.93	-0.06	0.14
	Merapi	0.08	0.37	0.32	0.71	0.29	0.04	0.02	0.30	0.28	0.68	0.31	0.02
	Meru	0.08	0.16	0.15	0.30	1.71	0.01	0.03	0.14	0.11	0.29	1.00	0.01
	Misti, El	0.07	0.29	0.25	0.59	0.46	0.02	0.02	0.27	0.19	0.60	0.32	0.06
	Momotombo	0.11	0.60	0.73	0.95	-0.50	0.04	0.03	0.52	0.70	0.90	-0.41	0.11
	Ngauruhoe	0.12	0.47	0.53	0.73	-0.41	0.01	0.05	0.40	0.46	0.70	-0.24	0.06
	Orizaba, Pico de	0.06	0.22	0.18	0.50	0.90	0.01	0.01	0.23	0.25	0.50	0.20	0.05
	Peteroa	0.04	0.16	0.13	0.35	1.37	0.01	0.01	0.17	0.17	0.35	0.63	0.02
	Pinatubo	0.10	0.67	0.78	0.96	-0.90	0.02	0.16	0.63	0.69	0.91	-0.76	0.01
	Quizapu	0.06	0.21	0.17	0.49	0.93	0.01	0.02	0.18	0.15	0.44	0.53	0.01
	Reventador	0.07	0.22	0.18	0.51	1.06	0.01	0.02	0.19	0.16	0.48	0.78	0.01
	Ruapehu	0.04	0.24	0.16	0.71	1.31	0.01	0.04	0.24	0.19	0.69	1.38	0.01
	Ruiz, Nevado del	0.03	0.16	0.12	0.49	2.08	0.01	0.03	0.21	0.19	0.47	1.08	0.01
	Sangay	0.08	0.27	0.26	0.51	0.26	0.02	0.01	0.25	0.27	0.51	0.01	0.09
	San Miguel	0.07	0.38	0.33	0.74	0.20	0.05	0.04	0.33	0.27	0.67	0.18	0.08
	San Salvador	0.12	0.41	0.40	0.71	0.08	0.01	0.11	0.37	0.35	0.66	0.29	0.01
	Santa María	0.09	0.27	0.28	0.46	0.03	0.03	0.02	0.24	0.29	0.47	-0.12	0.10
	Semeru	0.08	0.29	0.25	0.57	0.33	0.03	0.02	0.21	0.23	0.43	0.06	0.06
	Sinabung	0.11	0.44	0.50	0.76	-0.20	0.04	0.03	0.41	0.52	0.75	-0.22	0.09
	Socompa	0.05	0.18	0.14	0.42	1.32	0.01	0.02	0.13	0.12	0.29	1.09	0.01
	Soufrière, La	0.06	0.32	0.32	0.62	0.24	0.02	0.03	0.30	0.32	0.56	-0.02	0.02
	Soufrière Hills	0.08	0.42	0.40	0.87	0.30	0.02	0.05	0.36	0.36	0.71	0.10	0.02
	Spurr	0.05	0.20	0.16	0.46	0.91	0.01	0.01	0.17	0.16	0.46	1.26	0.05
	St. Helens	0.06	0.19	0.14	0.54	1.82	0.01	0.05	0.21	0.17	0.58	1.46	0.01
	Taranaki	0.09	0.46	0.39	0.87	0.16	0.08	0.02	0.40	0.37	0.85	0.14	0.06
	Teide	0.06	0.18	0.16	0.37	1.62	0.01	0.01	0.18	0.17	0.39	0.30	0.01
	Tolima	0.05	0.16	0.15	0.34	0.93	0.01	0.02	0.20	0.20	0.43	0.30	0.04
	Tungurahua	0.06	0.26	0.23	0.51	0.30	0.02	0.01	0.23	0.23	0.48	0.12	0.06
	Tutupaca	0.10	0.38	0.38	0.67	0.04	0.02	0.04	0.30	0.31	0.62	0.17	0.02
	Vesuvius	0.10	0.48	0.44	0.92	0.34	0.03	0.07	0.43	0.40	0.93	0.48	0.03
	1 Percentile 5.

2 Percentile 50.

3 Percentile 95.

Table S2 .

 S2 Main statistical parameters of the results describing channelization for each stratovolcano (2/2).

	Volcano			𝐒						𝐂 𝐅
		P5 1	Mean	P50 2	P95 3	Sk 4	HD 5	P5 1	Mean	P50 2	P95 3	Sk 4	HD 5
	Akagi	0.34	0.72	0.77	0.93	-0.76	0.01	0.23	0.58	0.61	0.84	-0.36	0.01
	Asakusa	0.44	0.73	0.77	0.91	-0.81	0.01	0.30	0.54	0.54	0.78	-0.03	0.01
	Asama	0.52	0.76	0.78	0.93	-0.63	0.01	0.39	0.60	0.61	0.80	-0.04	0.01
	Bandaisai	0.50	0.79	0.83	0.96	-0.93	0.01	0.38	0.64	0.66	0.84	-0.43	0.01
	Calbuco	0.28	0.60	0.61	0.91	-0.06	0.01	0.20	0.45	0.45	0.70	0.41	0.01
	Ceboruco	0.45	0.82	0.90	0.96	-1.17	0.01	0.32	0.69	0.74	0.95	-0.51	0.01
	Chaiten	0.50	0.82	0.88	0.97	-0.94	0.01	0.39	0.70	0.77	0.91	-0.51	0.03
	Chichon, El	0.45	0.80	0.87	0.96	-1.20	0.01	0.34	0.69	0.71	0.96	-0.34	0.02
	Chillán, Nevados de	0.33	0.61	0.60	0.91	0.07	0.01	0.26	0.46	0.42	0.76	0.74	0.01
	Chimborazo	0.38	0.69	0.72	0.95	-0.25	0.01	0.27	0.51	0.49	0.78	0.21	0.01
	Chokai	0.39	0.69	0.72	0.94	-0.39	0.01	0.27	0.51	0.51	0.79	0.13	0.01
	Colima, Nevado de	0.29	0.56	0.52	0.89	0.34	0.02	0.18	0.40	0.34	0.76	0.77	0.01
	Cotopaxi	0.39	0.74	0.80	0.95	-0.57	0.02	0.27	0.57	0.60	0.83	-0.13	0.03
	Fuego	0.36	0.66	0.68	0.90	-0.22	0.03	0.23	0.46	0.45	0.70	0.04	0.05
	Galeras	0.35	0.60	0.62	0.86	0.00	0.01	0.23	0.43	0.40	0.71	1.50	0.01
	Guallatiri	0.42	0.74	0.77	0.94	-0.50	0.01	0.34	0.59	0.58	0.82	-0.04	0.02
	Haku, Mount	0.33	0.57	0.54	0.89	0.49	0.01	0.23	0.41	0.36	0.79	1.16	0.01
	Kelut	0.28	0.74	0.84	0.95	-0.97	0.01	0.21	0.64	0.73	0.94	-0.50	0.02
	Lascar	0.50	0.76	0.79	0.94	-0.68	0.01	0.37	0.61	0.61	0.83	0.09	0.01
	Machin, Cerro	0.45	0.75	0.78	0.94	-0.76	0.01	0.36	0.60	0.59	0.86	0.09	0.01
	Mayon	0.46	0.82	0.88	0.99	-0.75	0.02	0.31	0.68	0.71	0.92	-0.41	0.04
	Merapi	0.31	0.71	0.73	0.96	-0.51	0.02	0.24	0.54	0.53	0.81	-0.06	0.03
	Meru	0.61	0.81	0.82	0.93	-1.08	0.01	0.36	0.61	0.62	0.81	-0.35	0.01
	Misti, El	0.35	0.69	0.71	0.94	-0.24	0.04	0.26	0.53	0.50	0.79	0.22	0.06
	Momotombo	0.56	0.87	0.96	0.99	-1.36	0.01	0.45	0.77	0.86	0.92	-0.90	0.02
	Ngauruhoe	0.65	0.86	0.89	0.96	-1.30	0.01	0.48	0.70	0.72	0.86	-0.57	0.01
	Orizaba, Pico de	0.31	0.60	0.60	0.86	-0.07	0.02	0.21	0.40	0.39	0.64	0.60	0.01
	Peteroa	0.32	0.59	0.58	0.89	0.20	0.01	0.25	0.44	0.42	0.75	0.83	0.02
	Pinatubo	0.37	0.86	0.95	0.97	-1.85	0.01	0.28	0.81	0.91	0.99	-1.45	0.02
	Quizapu	0.42	0.66	0.67	0.86	-0.26	0.01	0.32	0.47	0.46	0.64	0.90	0.01
	Reventador	0.46	0.70	0.73	0.88	-0.56	0.01	0.27	0.50	0.51	0.74	0.00	0.01
	Ruapehu	0.35	0.64	0.62	0.95	0.14	0.03	0.22	0.47	0.41	0.90	0.81	0.01
	Ruiz, Nevado del	0.31	0.64	0.65	0.95	-0.13	0.02	0.17	0.49	0.47	0.93	0.43	0.03
	Sangay	0.42	0.69	0.73	0.90	-0.48	0.01	0.26	0.44	0.45	0.62	-0.03	0.01
	San Miguel	0.44	0.79	0.82	0.97	-0.69	0.02	0.34	0.66	0.66	0.89	-0.27	0.04
	San Salvador	0.55	0.88	0.94	0.97	-1.95	0.01	0.35	0.80	0.90	0.95	-1.46	0.01
	Santa María	0.45	0.75	0.81	0.93	-0.77	0.01	0.32	0.56	0.59	0.73	-0.48	0.01
	Semeru	0.38	0.68	0.71	0.91	-0.27	0.03	0.24	0.46	0.48	0.67	-0.09	0.03
	Sinabung	0.53	0.84	0.93	0.98	-0.98	0.02	0.40	0.70	0.77	0.90	-0.58	0.03
	Socompa	0.40	0.65	0.67	0.84	-0.40	0.01	0.27	0.42	0.42	0.60	0.55	0.01
	Soufrière, La	0.40	0.74	0.79	0.92	-0.84	0.01	0.28	0.57	0.59	0.80	-0.26	0.01
	Soufrière Hills	0.40	0.76	0.81	0.97	-0.91	0.02	0.31	0.64	0.67	0.87	-0.53	0.01
	Spurr	0.44	0.66	0.65	0.90	0.11	0.01	0.26	0.43	0.39	0.74	0.97	0.01
	St. Helens	0.48	0.76	0.78	0.93	-0.83	0.01	0.31	0.59	0.59	0.85	-0.08	0.02
	Taranaki	0.34	0.73	0.77	0.98	-0.40	0.01	0.24	0.57	0.57	0.86	-0.10	0.04
	Teide	0.49	0.73	0.74	0.91	-0.46	0.01	0.31	0.51	0.50	0.73	0.17	0.01
	Tolima	0.39	0.64	0.64	0.88	-0.01	0.01	0.26	0.46	0.44	0.69	0.40	0.01
	Tungurahua	0.31	0.61	0.61	0.88	-0.09	0.03	0.22	0.42	0.41	0.66	0.24	0.01
	Tutupaca	0.51	0.80	0.85	0.95	-0.99	0.01	0.39	0.64	0.66	0.85	-0.33	0.01
	Vesuvius	0.70	0.88	0.90	0.99	-1.33	0.01	0.53	0.76	0.77	0.92	-0.61	0.01
	1 Percentile 5.

Percentile 50.

Percentile 95.

Skewness.

Dip statistic, derived from the application of the Hartigan's test for unimodality[START_REF] Hartigan | The Dip Test of Unimodality[END_REF]. Values greater than 0.035 imply that the distribution is clearly multimodal.

Acknowledgements

We thank Domenico Doronzo and an anonymous reviewer for useful comments. This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001 . This is Laboratory of Excellence ClerVolc contribution number XXX.