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)

The San Salvador volcanic complex (El Salvador) and Nejapa-Chiltepe volcanic complex (Nicaragua) have been characterized by a significant variability in eruption style and vent location. Densely inhabited cities are built on them and their surroundings, including the metropolitan areas of San Salvador (∼ 2.4 million people) and Managua (∼ 1.4 million people), respectively. In this study we present novel vent opening probability maps for these volcanic complexes, which are based on a multi-model approach that relies on kernel density estimators. In particular, we present thematic vent opening maps, i.e., we consider different hazardous phenomena separately, including lava emission, small-scale pyroclastic density currents, ejection of ballistic projectiles, and low-intensity pyroclastic fallout. Our volcanological dataset includes: (1) the location of past vents, (2) the mapping of the main fault structures, and (3) the eruption styles of past events, obtained from critical analysis of the literature and/or inferred from volcanic deposits and morphological features observed remotely and in the field. To illustrate the effects of considering the expected eruption style in the construction of vent opening maps, we focus on the analysis of small-scale pyroclastic density currents derived from phreatomagmatic activity or from low-intensity magmatic volcanism. For the numerical simulation of these phenomena we adopted the recently developed branching energy cone model by using the program ECMapProb. Our results show that the implementation of thematic vent opening maps can produce significantly different hazard levels from those estimated with traditional, non-thematic maps.

hazard is based on the determination and analysis of a number of expectable eruptive scenarios (e.g., Cioni et al., 2008;Neri et al., 2008;Martí et al., 2012;Ferrés et al., 2013;Wright et al., 2019). These scenarios are typically based on the characteristics of past eruptions, while their potential effects in the surrounding area can be evaluated by coupling field observations, uncertainty quantification and numerical simulations (e.g., Bayarri et al., 2015;Neri et al., 2015;Bevilacqua et al., 2017b;35 Rutarindwa et al., 2019). In addition to the expected magnitude and eruption style, other sources of uncertainty for volcanic hazard assessment may derive from the unknown position of the vent during future eruptions, especially for volcanic fields associated with large calderas or controlled by regional tectonics (e.g., Connor et al., 2019). Accordingly, vent opening probability maps have become fundamental tools for the assessment of hazard and risk in different volcanic systems (e.g. Alberico et al., 2002;Bevilacqua et al., 2015;Bevilacqua et al., 2018;Poland and Anderson, 2020). The procedures adopted 40 to construct these maps typically consider the data associated with vent position of past eruptions and, sometimes, major fault structures. Thus, they require the availability of a detailed data set of geological and volcanological information (Marzocchi and Bebbington, 2012;Connor et al., 2015;Németh and Kereszturi, 2015).
Vent opening maps have been developed for a growing number of volcanoes, such as Campi Flegrei (Alberico et al., 2002;Selva et al., 2012;Bevilacqua et al., 2015;Bevilacqua et al., 2017b), Somma-Vesuvius (Tadini et al., 2017a;Tadini et al., 45 2017b), Etna (Cappello et al., 2012;Del Negro et al., 2020), Yucca Mountain Region (Connor and Hill, 1995;Connor et al., 2000), Main Ethiopian Rift (Mazzarini et al., 2013;Mazzarini et al., 2016), Auckland Volcanic Field (Bebbington, 2013;2015;Ang et al., 2020), Long Valley Volcanic Region Bevilacqua et al., 2018), El Hierro Island (Becerril et al., 2013), Pacaya Volcanic Complex (Rose et al., 2013), Harrat-Rahat Volcanic Field (Runge et al., 2014), Snake River Plain (Gallant et al., 2018), Taupo Volcanic Zone (Kósik et al., 2020), among others. Furthermore, the expected style of activity 50 and resulting volcanic phenomena are likely influenced by vent position. In other words, vent position is not only important for volcanic risk assessment because it controls the final dispersal of volcanic products, but also because it can be correlated with eruption style and intensity. A number of different volcanic processes may produce this correlation, such as the involvement of different portions of the plumbing system, the presence of groundwater or surface water in specific zones of the volcanic field, or variations in the mechanical characteristics of the country rocks (e.g., Andronico et al., 2005;Coppola et 55 al., 2009;Aravena et al., 2020b). This is evident for example in partially submerged calderas, in which the style of activity of vents in the submerged zones is clearly influenced by their position. Thus, the spatial distribution of the expected style of activity could be strongly relevant when vent opening maps are used to produce long-term hazard maps for different eruption scenarios.
In this context, in order to deal simultaneously with the uncertainty associated with vent position, eruptive style and their 60 interdependence, here we present thematic vent opening maps that consider separately the occurrence of four selected volcanic phenomena: lava emission, small-scale pyroclastic density currents, ejection of ballistic projectiles, and low-intensity pyroclastic fallout. The input geological information adopted for the construction of these maps includes the mapping of the main fault structures, the position of past eruptive vents, and various information on the eruption styles and processes. This is deduced from a critical analysis of literature, the study of volcanic deposits, and from remote analysis of their morphological characteristics. The resulting thematic probability maps include uncertainty quantification both on past vents locations and on model parameters. To illustrate the effects of using thematic vent opening maps instead of ordinary non-thematic maps, we show a systematic application of the recently developed branching energy cone model (Aravena et al., 2020a), which is able to describe the invasion area of pyroclastic density currents considering channelization processes.
In this investigation, done in the ambience of the RIESCA project, funded by the Italian Development Cooperation and aimed 70 at promoting applied training in risk scenarios in Central America, we construct thematic vent opening maps for two volcanic complexes of the Central America Volcanic Arc: San Salvador Volcanic Complex (SSVC, El Salvador; Figure 1a) and Nejapa-Chiltepe Volcanic Complex (NCVC, Nicaragua; Figure 1b). SSVC and NCVC are both associated with the subduction of the Cocos Plate beneath the Caribbean Plate at an average rate of 70-90 mm/year (Barckhausen et al., 2001;DeMets, 2001)  This implies that any future volcanic activity could threaten the life of hundreds of thousands of people. Thus, in addition to be an interesting application of our new methodology, the availability of vent opening maps becomes relevant for improving hazard assessment in these regions.
In this paper, after the description of the geological framework of SSVC and NCVC (Sect. 2), we present new thematic vent 80 opening maps for the two studied volcanic complexes (Sect. 3). Then, we present invasion probability maps of small-scale pyroclastic density currents (PDC), which are derived from the application of the branching energy cone model (program ECMapProb) and the adoption of the thematic vent opening maps described previously (Sect. 4). Finally, Sect. 5 and Sect. 6 present the discussion and conclusions of this study.

Geological framework and volcanological dataset 85
We performed a detailed literature analysis in order to collect volcanological data and structural information on SSVC and NCVC (Figure 1a-b).
First, we collected information about past vent locations and their uncertainty areas, which were defined using polygons ( Figure   1a-b). Inside these polygons the probability is distributed uniformly. Small polygons are associated with well constrained vent locations (e.g., by the existence of a well-preserved crater), while large polygons are related to significant uncertainty in vent 90 position. Figure 1c-d presents the eruptive sequence of the two volcanic complexes addressed in this work, where blue symbols represent their well-documented eruptions (see Tables 1 and 2) and black symbols represent regional stratigraphic markers. In this figure, eruption numbers refer to the source vents (see Fig. 1a-b and Tables 1 and 2) while the letters are used to identify different eruptions from a single vent. In both the case studies there are many monogenetic vents and only one polygenetic vent, i.e. San Salvador -Boqueron volcano and Apoyeque volcano, respectively.
Then, we considered the major fault structures present in the studied volcanic areas (Figure 1a-b), which likely had a significant influence in the vent position of past events because of their effect in magma ascent dynamics (Gaffney and Damjanac, 2006;Valentine and Krogh, 2006;Ferrés et al., 2011;Avellán et al., 2012;Le Corvec et al., 2013).
Finally, we considered the range of eruptive phenomena associated with past events, mainly based on the nature of the different volcanic edifices and deposits, and on the available volcanological literature (see Tables 1 and 2). In this way, for each vent, 100 we defined one of four probability levels, i.e. Y = 1, PY = 0.75, PN = 0.25, and N = 0, associated with the occurrence of four different eruptive phenomena: (1) lava emission, (2) small-scale pyroclastic density currents associated with phreatomagmatic activity or low-intensity magmatic volcanism, (3) ejection of ballistic projectiles, and (4) low-intensity pyroclastic fallout due to Strombolian, violent Strombolian or Vulcanian activity (see Tables 1 and 2 and their footnotes). We remark that, although sub-Plinian to Plinian eruptions able to produce large-scale pyroclastic flows and fallout deposits are also present in the 105 geological record of SSVC and NCVC, this eruption type was only associated with the two main central edifices of San Salvador volcano and Apoyeque volcano, respectively. Thus, the probability of vent opening for these phenomena is largely concentrated on the apical part of these edifices and, consequently, the construction of vent opening maps for these events is not considered in this work.

San Salvador Volcanic Complex (SSVC) 110
The Pleistocene-Holocene San Salvador Volcanic Complex (SSVC, El Salvador; Figure 1a) sits in a EW trending graben and is formed by an active stratovolcano called Boqueron Volcano (BV) located inside an ancient volcanic edifice (San Salvador Volcano, SSV), and by at least 25 monogenetic volcanic edifices located at the SE, NW, and N of the main edifice ( Figure   1a). The monogenetic edifices include scoria cones, maars, tuff rings, tuff cones, and explosion craters. Their positions are strongly influenced by two NW-trending normal faults ( Figure 1a). In this work, following Ferrés et al. (2013) and Ferrés 115 (2014), these structures are referred as faults A (i.e., the N40W-trending fault) and B (i.e., the N65W-trending fault; see Figure   1a). Ferrés et al. (2013) and Ferrés (2014) identified three main periods in the eruptive history of SSVC: (a) Stage I (>72 ka -36 ka) corresponds to the construction and collapse of SSV. It is represented by pyroclastic deposits intercalated with andesitic and basaltic andesitic lavas, and fallout sequences associated with the Coatepeque caldera (Kutterolf 120 et al., 2008;Ferrés et al., 2011). El Picacho (~1960 m a.s.l.), El Jabali (~1400 m a.s.l.) and the SW portion of the current volcanic edifice are the present remnants of the collapsed and deeply eroded SSV ( Fig. 1a; Ferrés, 2014). No monogenetic centers have been identified during this stage (Ferrés, 2014). The collapse of the SSV was a consequence of a phreato-Plinian eruption responsible of the emplacement of the pyroclastic sequence G1 (event 25a in the nomenclature adopted in Fig. 1c; Sofield, 1998). 125 (b) Stage II (36 ka -3 ka) represents the construction of the active Boqueron Volcano (volume of circa 8.5 km 3 ; Ferrés, 2014). It includes both effusive and explosive volcanism associated with the emission of basaltic andesitic and andesitic magmas, intercalated with pyroclastic deposits from the Ilopango caldera (Kutterolf et al., 2008;Ferrés et al., 2011;Smith et https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License. al., 2020). During this period, lava flows were issued from the SSV (in particular, fourteen lava flows were recognized by Fairbrothers et al., 1978), while the explosive activity included at least six events after the pyroclastic sequence G1 (events 130 25b-g in the nomenclature of Fig. 1c, the last one marked the limit between Stages II and III; Ferrés et al., 2013).
(c) Stage III (3 kapresent) includes eruptive activity mainly concentrated in the flanks of the BV and in the plain nearby. Sofield (1998) suggests that this shift was related to the reaching of the critical height of the younger edifice (i.e., BV). Several monogenetic edifices were formed during this period as a consequence of both explosive and effusive volcanism, which frequently involved magma interaction with external water, such as Crater La Escondida and Loma Caldera (Table 1). In any 135 case, no hydromagmatic explosions have occurred in the area in the last 1000 years (after Talpetate I, see Fig. 1c; Sofield, 1998), possibly as a consequence of a reduction in the level of groundwater in the flanks of the volcano. Ferrés (2014) suggests that the enlargement of the volcanic field during Stage III is associated with a general decrease of activity in the Central Crater. It is important to note that, although flank activity started approximately 10 ka ago, it has been dominant only during the last 3 ka, exhibiting a strong structural control (Fig. 1a). On the other hand, during this period, BV 140 volcanic activity has included at least three volcanic events: Talpetate I (San Andres Tuff, which includes fallout and surge deposits widely distributed toward SW, event 25h in the nomenclature of Fig. 1c), Talpetate II (event 25i in the nomenclature of Fig. 1c), and the last eruption (i.e., A.D. 1917 event, vent 29; Fig. 1c and Table 1), when minor explosive activity in the central crater, associated with the construction of the intra-crateric Boqueroncito tuff ring, coexisted with an important lava effusion from parasitic vents of the N flank (vents 27 and 28 in our nomenclature, Fig. 1c and Table 1). 145 Sofield (2004) proposed five different eruption scenarios: (1) hydromagmatic flank eruption of VEI 1-3, (2) monogenetic magmatic eruption of VEI 1-3, (3) small-scale eruption of VEI 1-3 within Boqueron crater, (4) sub-Plinian eruption from the central vent (VEI 4-5) and (5) Plinian eruption from the central vent (VEI 6). Major et al. (2001) also described plausible scenarios associated with lahar occurrence. A first hazard assessment exercise for BV was done by Ferrés et al. (2013), considering ash fall, ballistic projectiles, and pyroclastic density currents issued from the central vent. However, since a 150 recurrence period of 85±50 years was proposed by Sofield (1998) and Sofield (2004) for flank eruptions during the Stage III, the most probable future event is associated with monogenetic volcanism, and thus the analysis of the hazard associated with flank eruptions is of paramount importance in SSVC in order to complement the recent literature.
In this study, following Sofield (1998), Ferrés et al. (2011), and Ferrés (2014, we considered 29 vents in the SSVC ( Fig. 1a and Table 1), most of them with a monogenetic character (the only exception is vent 25, i.e., the central crater). Mainly 155 considering the edifice type, the deposits characteristics, and the historical activity, for each volcanic vent we defined the probability of occurrence of the four volcanic phenomena studied here (i.e., lava emission, small-scale PDCs, emission of ballistic projectiles, and low-intensity pyroclastic fallout; Table 1). In general terms, small-scale PDCs were not associated with the activity of scoria cones, while all the four volcanic phenomena considered here have been linked, although with different probability of occurrence, to maars, explosion craters, tuff cones, and tuff rings as well as to the central activity of 160 BV. It is important to highlight that the products of Plinian and sub-Plinian volcanism, which have been produced at the central https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License. vent of SSVC, were not considered in the construction of our vent opening maps due to the small uncertainty in vent position for these events.
NCVC deposits are intercalated with pyroclastic deposits derived from Apoyo and Masaya calderas, which represent useful stratigraphic markers (Figure 1d; Kutterolf et al., 2008). Apoyeque products are mainly dacitic to rhyolitic in composition, while the products of monogenetic vents range from basaltic to andesitic basaltic (Avellán et al., 2012). Most of the eruptions associated with NCVC involved hydromagmatic activity, possibly triggered by the interaction between rising magma and 175 shallow aquifers (Avellán et al., 2012). The possible presence of sources of surface water is so of primary importance for volcanic hazard assessment in this area.
No clear trends are recognized in the temporal evolution of vent position and eruption style, even if important climate variations have been documented in this zone during the Holocene (Freundt et al., 2010). In any case, it is important to highlight that phreatomagmatic activity in central and southern NCVC has been dominantly observed relatively near the Nejapa fault (e.g., 180 El Plomo craters, Ticomo craters, Refineria crater), while scoria cones are preferentially observed in the peripheral zone (e.g.,  Table 2). Also in this case, for each volcanic vent we defined the probability of occurrence of the four 190 volcanic phenomena studied in this work, mainly based on the edifice type, the characteristics of the documented deposits, and the current presence of surface water (Table 2). In general terms, scoria cones were mainly related to the production of smallvolume fallout deposits, ejection of ballistic projectiles and the emission of lavas. On the other hand, the hydromagmatic activity typical of maars, tuff cones, and tuff rings has been linked to the emission of ballistic projectiles and generation of https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License.
relatively small PDCs and fallout deposits; while lava flow activity is generally absent at this type of edifices. Finally,in 195 Apoyeque volcano (i.e., vent 31, Fig. 1d and Table 2), we consider that the typical volcanism is characterized by the emission of ballistic projectiles and the eventual generation of relatively small volumes of fallout deposits and small-scale PDCs. We remark that the analysis of large-scale explosive eruptions, which have characterized the activity of Apoyeque volcano, is beyond the objective of this work and thus they were not considered in the construction of our vent opening maps.

Methods
In both the case studies of SSVC and NCVC, we adopt the kernel-based multi-model approach developed in Bevilacqua et al. (2015), Bevilacqua et al. (2017a) and Tadini et al. (2017b). These are long-term assessments based on the record of past eruptions and mapped faults, which lie on the assumption that a new vent will likely open close to previous vents and to geological structures that, in the past, have favoured the ascent of magma to the surface. 205 We linearly integrate two models named Model 1 and Model 2, with weights affected by uncertainty (Fig. 2). Model 1 considers the faults and the position of past vents following a Bayesian approach, while Model 2 adopts a Gaussian kernel density estimator applied on a uniform distribution within the uncertainty areas enclosing the past vents. We remark that in both the models past vents do not comprise simple points, but areas of uncertainty of different extent. Each area can cover several cells of our computational grid, some of them completely, others only partially. Therefore, for each cell, it was taken into account 210 the fraction of each uncertain area that it contains Tadini et al., 2017a). Both models are reviewed in the Appendix A. A key novelty of this study is the fact that we weight the past vents differently according to the chance that the new vents related to them will be the source of a specified hazardous phenomenon, based on the past eruptions locally occurred. The thematic weights, which leads to the production of thematic maps, are reported in Tables 1 and 2.
The location of the next eruptive vent is modelled as a random variable X with a continuous probability density function (PDF) 215 in the spatial domain. The vent opening map is then displayed as the probability density of the variable X per km 2 . We remark that the probabilities of vent opening that we obtain are conditioned on the occurrence of a new eruption and there is no associated temporal window. We did not count the polygenetic vents multiple times in these models. In Sect. 5 we discuss this assumption and we test the opposite choice.
Moreover, our estimates are doubly stochastic (e.g., Cox and Isham, 1980;Daley and Vere-Jones, 2003;Jaquet et al., 220 2008;Jaquet et al., 2012;, which means that the statistical distribution of the location of the next eruptive vent is represented using ill-constrained parameters that are treated as uniformly distributed random variables (e.g., Bevilacqua, 2016).
These parameters are two distance values 1 and 2 tuning the kernel bandwidth of the two models, and two probability values 1 and 2 tuning the probabilistic relevance of Model 1 compared to Model 2 and of the mapped faults compared to unknown structures, respectively. As a consequence of this approach, the PDF values will have their own confidence intervals. 225 https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License.
In particular, the two map layers associated with different conceptual models (Model 1 and Model 2) are linearly combined with specific weights, 1 and (1 − 1 ), for the development of a multi-model vent opening map (Fig. 2). A similar approach is adopted inside Model 1 to define the prior fault map as the combination of two layers, weighted 2 and (1 − 2 ). One layer, 1 , is related to the mapped faults (i.e. fault outcrops), and the other layer, 2 , is a uniform PDF representing the unknown (i.e. buried) faults. Figure 2a illustrates the logic of this vent opening model combination. Figure 2b shows the two models and 230 their differences when applied to the test example of a past vent near a fault.
The input parameter 1 = Unif(0.75,1) was constrained through a straightforward expert judgment after the assumption that Model 1 should have a greater relevance than Model 2. This is because both SSVC and NCVC are characterized by robust geological information which suggests significant relationships between past vents, fault structures and volcanism, with several aligned vents along the main faults. Similarly, 2 = Unif(0.75,1) was chosen through expert judgment due to the greater 235 relevance of mapped faults compared to unknown structures that are uniformly distributed in the region. Thus, the vent opening maps produced here are dominated by Model 1 and by the information of mapped faults, but the possibility of a new vent not influenced by structural alignment or by an unknown fault is considered. Further research focused on the deeper understanding of the two case studies of SSVC and NCVC could improve these constraints, either through more structured expert judgment techniques Tadini et al., 2017a) or likelihood based techniques Bevilacqua 240 et al., 2018).
On the other hand, our multi-model approach depends on two additional parameters: 1 and 2 (see Appendix A). The parameter 1 in Model 1 is the average distance between sub-parallel regional faults. Different distances are measured in different sub-regions, thus defining the uncertainty range of 1 . Specifically, 1 = Unif(5,10) km for both SSVC and NCVC.
The parameter 2 in Model 2 is the average of the distance of the -th nearest neighbor of each past vent. The number is 245 varied in [1,3] to consider the presence of spatial clusters, thus defining the uncertainty range of 2 . Specifically, 2 = Unif(1,2.5) km for SSVC and 2 = Unif(1,2) km for NCVC.
In summary, if ( ) is the PDF of Model 1 and ( ) is the PDF of Model 2, then the multi-model "integrated" vent opening map is given by: 250 that is: where we also expressed the dependence of and on the mapped faults 1 and the uniformly distributed map of unknown faults 2 , and the distance parameters 1 and 2 .
https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License. percentile and the 95 th percentile, expressed in terms of probability density per km 2 , in percentage. For comparison purposes, Figure 5a presents the results associated with the mean value of a non-thematic vent opening map constructed using the same procedure but with all the vents adopting the same weight. Please note that, in practice, this map is equivalent to the vent opening map of ballistic projectiles (Fig. S1).

San Salvador Volcanic Complex
Results exhibit significant differences between the various thematic maps. While the highest vent opening probability 265 associated with small-scale PDCs follows the N65E-trending fault B (Figure 4), in the other thematic maps the highest vent opening probabilities tend to locate along the northern portion of the N40E-trending fault A (Figs. 3, S1 and S2).
In particular, in the case of events able to produce lava flows (Fig. 3), the maximum probability density is located on the NNW flank of San Salvador volcano along the fault A, reaching probabilities of up to 1.0% per km 2 in the mean value map, with 90% confidence interval [0.7%, 1.6%] per km 2 . Considering the mean value map, the total probability of vent opening in a 4 270 km wide belt across the northern portion of fault A is 39.3%, excluding the pixels that are closer to fault B. On the other hand, the total probability of vent opening close to fault B is 13.9%. This value was also computed by considering a 4 km wide belt across the fault under examination. Please note that this criterion is also used in the sequel.
Instead, as mentioned above, the region of maximum probability of vent opening conditioned on the occurrence of an eruption able to produce small-scale PDCs (Fig. 4) is located along the N65E-trending fault B, with a peak probability in the mean 275 value map of 1.1% per km 2 , with 90% confidence interval [0.8%, 1.6%] per km 2 . In this case, in fact, the total probability of vent opening near the fault B (i.e. 23.5%) is significantly higher than the results observed in the thematic map of lava emission.
Instead, the probabilities of vent opening along the northern and southern portions of fault A are 26.8% and 6.1%, respectively.
The thematic maps associated with ballistic projectiles (Fig. S1) are quite similar to the case of lava emission, with the maximum values observed along fault A. These maps show a peak probability density of 0.9% per km 2 in the mean value map, 280 with 90% confidence interval [0.7%, 1.3%] per km 2 . The probabilities of vent opening near the northern portion of fault A, the southern portion of fault A and fault B are 35.0%, 7.1% and 16.7%, respectively.
Finally, the vent opening maps related to eruptions able to produce small-scale fallout deposits (Fig. S2) are similar to those presented for lava emission and ballistic projectiles. The maximum probability of vent opening, located at the NW flank of the volcano, is 0.9% per km 2 , with 90% confidence interval [0.7%, 1.3%] per km 2 . In this case, the probabilities of vent opening 285 close to the northern portion of fault A, the southern part of fault A and fault B are 35.1%, 6.5% and 16.3%, respectively. https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License.

Nejapa-Chiltepe Volcanic Complex
The thematic vent opening maps associated with NCVC are presented in Figures 6 (lava flows), 7 (small-scale PDCs), S3 (ballistic projectiles), and S4 (small-scale fallout pyroclastic deposits). Again, each figure presents the mean value of probability per km 2 , in percentage, including the results of the 5 th percentile and the 95 th percentile as well. Instead, Figure 5b  290 shows the mean value of a non-thematic vent opening map, i.e. with no differences in the weight of the different vents. Also in this case, this map is equivalent to the vent opening map of ballistic projectiles (Fig. S3).
All the maps show the maximum probability near Asososca maar (vent number 13 in Fig. 1b,d and Table 2), and the only significant differences are related to the N-S extent of the high-probability zone, which tends to be longer in the maps associated with the occurrence of small-scale PDCs (Fig. 7). A secondary peak is observed near Miraflores scoria cone. 295 In particular, the peak value of vent opening probability density in the mean value map related to the emission of lavas is 1.6% per km 2 , with 90% confidence interval [1.0%, 2.3%] per km 2 (Fig. 6). The maximum values in the other maps are: 1.7% per km 2 for events able to produce small-scale PDCs, with 90% confidence interval [1.2%, 2.5%] per km 2 (Fig. 7); 1.6% per km 2 for ballistic projectiles, with 90% confidence interval [1.1%, 2.2%] per km 2 (Fig. S3); and 1.6% per km 2 for eruptions that produce small-scale fallout deposits, with 90% confidence interval [1.1%, 2.3%] per km 2 (Fig. S4). The resulting probabilities 300 of vent opening inside the limits of Managua are 28.9%, 35.1%, 31.2% and 32.0% for events able to produce lava flows, smallscale PDCs, ballistic projectiles and small-scale fallout deposits, respectively.
The previous hazard assessment exercise in Connor et al., (2019) produced spatial density maps which differ from our results.
As said above, they did not associate any uncertainty area to the past vent locations, they considered slightly less events, and they did not use the map of the main fault structures. Moreover, their method is not doubly stochastic, thus they did not quantify 305 the uncertainty affecting their results. In other words, their spatial map is not a PDF but a spatial density. However, once divided by the number of past vents considered, it becomes a vent opening map. Their maximum probability values are ~2.8-3.2% per km 2 and shifted of ~1 km towards west compared to ours. Their region above 1% probability per km 2 has a similar areal extent to ours, but is significantly stretched in the N-S direction, due to the use of an elliptical kernel. Further comparison would be needed to fully evaluate the strengths and weaknesses of both approaches. 310

Numerical simulation of small-scale PDC invasion hazard maps
For both SSVC and NCVC we extracted 100 sets of 1024 samples of vent position by using the thematic vent opening maps associated with small-scale PDCs. Each set of 1024 vent positions is associated with one of 100 different arrays of the model parameters ( 1 , 2 , 1 , 2 ) adopted in the construction of our thematic vent opening maps. The sampling method is based on a Latin Hypercube Sampling scheme, generalized with Orthogonal Arrays to increase space-filling properties (e.g., Bevilacqua 315 et al., 2019b;Patra et al., 2020). We modified the method to work under the assumption of a non-uniform PDF, as detailed in Appendix A.
Through the use of the branching energy cone model, implemented in the program ECMapProb and described in Aravena et al. (2020a), we performed these 100 x 1024 simulations (i.e., 102,400 simulations) for each volcanic complex. We fixed the other PDC initial conditions. In particular, we assumed = 500 m and tan( ) = 0.25. These input conditions are consistent 320 with the runout distances observed in the small phreatomagmatic deposits recognized in the field, e.g., Loma Caldera in SSVC, and are assumed to be representative of the studied PDCs. We decided not to use variable input conditions for initial PDC characteristics, so our hazard assessment is only valid in this specific scenario of PDC size and friction angle. A more complete PDC hazard assessment considering variable size and friction properties would require additional information to properly calibrate these input parameters (e.g., Cioni et al., 2020). Further analysis could investigate the sensibility of numerical results 325 on these parameters, and increase the number of inputs that are sampled in the LHS scheme, also considering correlations between PDC size and friction properties, but this is beyond the purpose of this study (e.g., Spiller et al., 2014;Ogburn et al., 2016;Ogburn and Calder, 2017;Rutarindwa et al., 2019;Patra et al., 2020). with a maximum value of 23%, with 90% confidence interval [19%, 29%] and high invasion probabilities near the cities of Nueva Ciudad del Niño and Lourdes (Fig. 8). Modelled invasion probability at the SE flank of San Salvador volcano and in the city of San Salvador tends to be low, where a significant effect is exerted by the topographic barrier of Cerro El Picacho (Fig. 8). The highest small-scale PDC invasion probability calculated at the metropolitan area of San Salvador, ~2-3%, is located in its western sector, i.e. Santa Tecla. This is favoured by the high slope of the main edifice and the absence of 340 significant topographic barriers in this direction.
Instead, in the case of NCVC, numerical results indicate the highest invasion probabilities along the Nejapa fault, near the western border of Managua city and inside it. The maximum values of invasion probability, i.e. 22% and 23% in the mean value map, are associated with the low-elevation areas of Asososca maar and Nejapa maar, respectively (Fig. 9). Nonnegligible invasion probabilities have been calculated at the eastern portion of Ciudad Sandino as well, while the invasion 345 probability of small-scale PDCs at the northern part of the studied zone, i.e. near Apoyeque volcano, tends to be low (Fig. 9).
Finally, for the two volcanic complexes addressed here, Figure 10 presents the difference between the mean PDC invasion probability percentages obtained by adopting thematic vent opening maps (i.e. Figs. 8a and 9a) and the equivalent results that would derive from the application of non-thematic vent opening maps, i.e. by using the maps displayed in where the application of thematic vent opening maps implies an increase in the calculated PDC invasion probability, while negative values (i.e. blue zones) are the sectors where the use of thematic vent opening maps instead of the traditional ones is manifested in a reduction of the computed PDC invasion probability. SSVC exhibits differences ranging from -7.2%, on the N-NNE flank of Boqueron Volcano (BV), near Lourdes, to +5.8% on the NW flank of BV, south of Quezaltepeque (Fig. 10a). 355 These results involve dramatic modifications in the computed PDC invasion probability. In fact, at the NW flank of BV, the PDC invasion probability is halved when thematic vent opening maps are considered, while this parameter increases in a 70% in some portions of the NW flank of the BV. On the other hand, the computed differences for NCVC are moderate and vary from -2.9% in the northern portion of NCVC to +3.2% in the southern portion of this volcanic field, including a significant area that belongs to the city of Managua. In any case, in relative terms the differences between the results derived from the 360 application of thematic and non-thematic maps are small both in Managua and in its surroundings.

Discussion
Several probabilistic assessments of vent opening at different volcanoes have been presented during the last decade (e.g., Marzocchi and Bebbington, 2012;Connor et al., 2015;Poland and Anderson, 2020), improving the methodologies implemented to construct these maps and sensibly increasing the number of volcanic fields for which a probabilistic assessment 365 of vent opening is available. It is remarkable, for example, the inclusion of structured expert judgment procedures to constrain the input parameters of the models adopted to construct these maps (Chapman et al., 2012;Bevilacqua et al., 2015;Tadini et al., 2017a;Bebbington et al., 2018), and their coupling with models aimed at describing the dispersal of volcanic products (e.g., Neri et al., 2015;Gallant et al., 2018;Hyman et al., 2019;Rutarindwa et al., 2019). A growing effort is aimed at the production of short-term vent opening maps which can modify the long-term estimates after plugin-in the monitoring 370 information that progressively evolves during volcanic unrest (e.g., Bevilacqua et al., 2019c;Patra et al., 2019;Bevilacqua et al., 2020b;Sandri et al., 2020).
In this context, in this work we adopted a novel approach where we have considered separately different volcanic hazardous phenomena (i.e. emission of lava flows, small-scale PDCs, ballistic projectiles, and small-scale pyroclastic fallout). This led to the construction of thematic vent opening maps which allow to assess the hazard that different volcanic phenomena involve 375 separately. We remark that previous studies already produced vent opening maps devoted to specific types of eruptions (e.g. Plinian and sub-Plinian eruptions in Tadini et al., 2017b), to eruptions inside selected sub-regions (Bevilacqua et al., 2017b) or to a suite of pre-imposed eruptive scenarios (Ang et al., 2020).
To do this, in addition to the mapped faults and the record of past eruptions, we considered the eruption style and the likely eruptive phenomena that the past eruptions could have produced, on the basis of deposits and morphological features of 380 volcanic edifices. This is justified by the suggested relation of vent position and eruption style and intensity, which has been observed at different volcanic fields (e.g., Andronico et al., 2005;Coppola et al., 2009). In particular for the studied volcanic fields, at SSVC, small-scale PDCs have been preferentially produced by vents located along the N40W-trending fault B (Fig. https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License.
1a and Table 1), while lava flows have been commonly observed in vents located along the N65W-trending fault A (Fig. 1a and Table 1). Apparently, no small-scale PDCs have been issued in the past from the central vent of the present Boqueron 385 volcano. On the other hand, at NCVC, maars and tuff rings have been mainly produced in events whose vent was located near the Nejapa fault, while scoria cones are predominant in the peripheral zones. Additional evidence is reported also for the possible presence of hydromagmatic eruptive centers located in the south-western corner of Managua Lake.
As expected, the integration of this information in the construction of thematic vent opening maps is manifested in their results (Figures 3, 4, 6, 7 and S1-S4), which in some cases present significant differences as a function of the considered hazardous 390 phenomenon. We also estimated how these changes are sensible to some of the main sources of uncertainty. In fact, the 90% range of the probability per km 2 can strongly enhance the differences between thematic maps as can be observed, for example, in Figures 3c and 4c.
The design and implementation of mitigation strategies of volcanic risk is inescapably associated with the characteristics of the volcanic process under consideration and vent position. In this sense, the importance of using thematic maps when the 395 studied volcanic complex shows a vent position-controlled eruption style has been illustrated by modeling the propagation dynamics of small-scale PDCs. Our results, which were obtained by adopting the branching energy cone model (Aravena et al., 2020a), show the coupled effect of vent opening maps and volcano topography in determining the PDC invasion hazard.
Results stress that the adoption of thematic vent opening maps is able to produce significant effects in the construction of hazard maps, which are particularly relevant for SSVC. 400 Finally, an important aspect in our analysis is that we did not count multiple times the polygenetic vents (San Salvador -Boqueron volcano in SSVC, and Apoyeque volcano in NCVC), i.e. every vent was weighted only in consideration of its past eruptive phenomena and not the number of eruptions. Similarly, in Connor et al. (2019) the Apoyeque volcano counted one in the production of their spatial density maps. In contrast, a straightforward event-counting approach could lead to different results, and this fact deserves some additional comments. 405 In the case of SSVC we counted nine explosive events after 36 ka, including the collapse of SSV at the end of Stage I, six eruptions in Stage II, and two in Stage III (Figure 1c). Moreover, in Stage II at least fourteen lava flows were issued from the SSV (Fairbrothers et al., 1978). In the NCVC we counted three explosive eruptions from Apoyeque, and no lava flows after 30 ka (Fig. 1d). In our results we did not follow an event counting approach because of several reasons explained in the sequel.
Firstly, an event counting approach would rely on the assumption that the volcanic system is stationary over the time period 410 considered in the statistics. In contrast, for example at 3 ka SSVC apparently shifted from a central volcanism to the development of the monogenetic field on the flanks of the BV and in the plain nearby. Most the twenty-eight known monogenetic events are after 3 ka, while only two of the BV eruptions are after 3 ka. In addition, how much the Central Crater should count depends on the time interval considered. In a similar situation, Bevilacqua et al. (2017a) treated the uncertainty related to a shift of volcanic activity from Long Valley Caldera to Mono-Inyo chain by weighting the events in the older site 415 from 0% to 50% of the more recent events. Bevilacqua et al. (2015) weighted differently the vents in the three past epochs of https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License.
Campi Flegrei post-collapse activity, after an expert elicitation. Again, the older vents tended to weigh less than the more recent vents.
Secondly, in the two case studies of SSVC and NCVC the specific knowledge of many eruptive centers is variable and sometimes poor, under-recording is likely, and a robust event counting is not possible. However, the possible errors in event 420 counting are not the same in the various thematic maps. For example, in the case of small-scale PDCs at SSVC the capability of producing low energy events in the central edifice should not be substantially different from monogenetic vents. In contrast, BV is a stratovolcano built after a large number of lava flows, while monogenetic vents typically produce a single flow. Figure 11 shows two illustrative examples of modified vent opening maps of SSVC after following an event counting-based approach. We only present the results of SSVC because the potential effect of this modification would be minor in NCVC, 425 due to the relatively small number of known polygenetic events. In particular, Figure 11a shows a vent opening map for lava flow emission that counts the central crater 16.25 times, i.e. 14 times for the past lava flows in Stage II, plus 9 x 0.25 times for the past explosive eruptions after 36 ka. Figure 11b shows a vent opening map for ballistics that counts the central crater 9 times, i.e. the past explosive eruptions. The significant concentration of vent opening probability towards the polygenetic center is evident in both figures. Figure S5 and S6 in the Supporting Information show the 5 th and the 95 th percentile values of 430 these maps. Conversely, in terms of hazard assessment, our choice (i.e. with no consideration of a higher weight for polygenetic vents) spreads the vent opening probability further from the central volcano and potentially closer to highly inhabited areas.
In summary, to improve this aspect of the analysis, separate maps for the central volcano and the surrounding volcanic field could be produced, and their differential weight better constrained after additional research.

Conclusions 435
In this study, we have presented vent opening probability maps for two Central American volcanic fields: San Salvador Volcanic Complex and Nejapa-Chiltepe Volcanic Complex. These volcanic fields present some features that make critical the availability of this tool for the design and implementation of volcanic risk mitigation procedures: (1) they are next to highly inhabited cities, i.e. San Salvador (2.4M people) and Managua (1.4M people), respectively; (2) they present a significant variability in vent position; and (3) they produced a number of different types of eruption style. 440 We implemented quantitative vent opening probability estimates which are based on main fault structures and past vent positions, and include uncertainty quantification and the type of activity of the documented eruptions. We produced different thematic maps related to emission of lavas, small-scale pyroclastic density currents, ballistic projectiles and small-scale fallout deposits.
The main findings derived from the new vent opening maps are: 445 (a) At SSVC, the maps show their maximum values on the NW flank of San Salvador volcano and in the northern region of this volcanic complex. In particular, results indicate that the zone of high probability of vent opening (>0.5% mean probability per km 2 ) of events able to produce small-scale PDCs follows the N40W-trending fault B, while the high-probability zone https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License. related to lava emission, ballistic projectiles and the production of small-volume fallout deposits follows the N65W-trending fault A. This shows that the separate consideration of different hazardous phenomena is able to produce relevant differences 450 in vent opening probability maps.
(b) At NCVC, results show that the maximum vent opening probability is located near Asososca maar. The only differences between the different thematic maps is the N-S extent of the high-probability zone, which tends to be larger in the maps associated with small-scale PDCs than in the other thematic maps. Importantly, a significant portion of the vent opening probability distribution is located inside the limits of Managua City, which implies major challenges for dealing with the 455 volcanic hazard associated with NCVC.
We remark that we did not consider the occurrence of large volume pyroclastic deposits derived from sub-Plinian or Plinian eruptions, e.g., large-volume fallout or PDC deposits. Indeed, in the past they were only emitted from the polygenetic central vents of both volcanic complexes, i.e. San Salvador and Apoyeque, respectively. Thus in these cases the construction of vent opening probability maps is not meaningful at the scale of the other thematic maps, i.e. tens of kilometers. 460 Finally, we have adopted the branching energy cone model described in Aravena et al. (2020a) to illustrate the use of the thematic vent opening maps presented here. In particular, we performed a Monte Carlo simulation of ~10 5 small-scale PDCs, representative of those derived in the past from phreatomagmatic activity or low-intensity magmatic volcanism. We varied the vent location following a doubly stochastic approach that enabled us to estimate the effects of the uncertainty affecting the vent opening PDF. We did not change the size and the friction properties of the flow to better analyse the effects of the different 465 vent opening maps.  The findings presented in this paper represent relevant information for the management of volcanic hazards in these volcanic complexes, and complement previous works focused on the hazard assessment of activity in the central vents (e.g., Avellán et al., 2012;Ferrés et al., 2013). Coupling of thematic vent opening probability maps with numerical models able to describe other eruptive phenomena, e.g. lava flows and ballistics, would provide further relevant data for the assessment of volcanic 485 hazard, as shown here for PDCs.

Appendix A
In the following we review the details Model 1 and Model 2, and the non-uniform Latin Hypercube Sampling scheme adopted in this work.

A.1 Model 1: Kernel-based Bayesian update of fault map 490
Model 1 merges the faults and past vents locations following a Bayesian approach. This model is also detailed in Bevilacqua et al. (2017a), but differs in the method to select the distance parameter 1 . We assume that the mapped faults provide information in the forecast of future vent locations (e.g., Felpeto et al., 2007;Cappello et al., 2012;Cappello et al., 2015). In fact, we expect the new vents to more likely occur near mapped faults that previously interacted with rising magma (Martin et al., 2004;Jaquet et al., 2012). The key idea in Model 1 is that the vent opening map depends on an additional spatial parameter 495 = ( 1 , 2 ) representing the outcrop of a fault interacting with a future potential rising dike. We call the prior PDF of , i.e. the prior fault map. Namely, = 1 is the distribution of mapped faults, and = 2 is the distribution of unknown faults.
Considering that linear averaging is not commutative with the Bayes theorem, we average the posterior probability maps because, in this way, the linear weights directly impact the final results, and they are not altered by the multiplicative step.
For each past vent location ( , ), we use Bayes theorem to calculate the posterior probability density values ( 1 , 2 │ , ), 500 up to a multiplicative normalization constant detailed below. This represents the fault locations that may have interacted with the rise of magma that led to the past vent ( , ). We integrate past vent locations ( = 1, … , ) on their uncertainty areas ( ), according to the expression: where (• │ •) is a Gaussian likelihood defined below (Eq. 4). We calculate different posterior probabilities, one for each past 505 vent considered ( = 1, … , ), and we define the global posterior of as the weighted average of them, with weights = , where is related to the specified hazardous phenomenon (see Tables 1 and 2). Different hazards will lead to different posterior maps.
Then, conditioned on the value of parameter , we define the likelihood for a vent location x as a symmetrical, two dimensional Gaussian function of mean and covariance matrix 2 : This kernel function creates a link between the past vents and the portions of the faults next to them, i.e. . The value of the standard deviation is = 1 /2, where 1 is the average distance between sub-parallel regional faults. In other words, we assume that about 95% of the vents that are related to a fault will open closer than the next sub-parallel fault. This is a different strategy from that adopted by Bevilacqua et al. (2017a), where the distance bound was obtained through a geometric argument 515 on the maximum depth of the fault-dike interaction and the dip angle of the faults.
Finally, we calculate the vent opening map applying again the kernel function (• │ •) to the global posterior probability distribution of . So, the kernel function creates a link between and the new vents. Thus, the vent opening probability density according to Model 1 is defined by: where with = 1,2 are the two posterior probability density functions of , 1 = 2 , and 2 = 1 − 2 .

A.2 Model 2: Kernel density estimator
Model 2 is a Gaussian kernel density estimation applied on a uniform distribution within the uncertainty areas enclosing the past vents (Bebbington and Cronin, 2011;Connor et al., 2012;Bevilacqua et al., 2015). This method is detailed in Tadini et al. (2017b), but differs in the selection of the distance parameter 2 . The kernel function describes the expected distance of 525 propagation of new eruptive vents from pre-existing vents. The kernel function can be any positive function that integrates to one (Weller et al., 2006), and in general, given a finite sample = ( , ), a kernel density estimator can be defined as follows: where ℎ is the bandwidth, and are the vent weights related to the specified hazardous phenomenon. We assume equal to 530 a two-dimensional radially symmetric Gaussian function, and ℎ is the standard deviation parameter. In this study we consider ℎ = 2 2 / , where 2 is the average of the distance of the -th nearest neighbor of each past vent, formally represented by a Rayleigh distribution. This is a different procedure from that adopted by Tadini et al. (2017b), where = 1 and ℎ = 2 .

A.3 Non-uniform Latin Hypercube Sampling
Our sampling technique of the input variables in the PDC simulations is based on the Latin Hypercube Sampling (LHS) idea, 535 and in particular, on the improved space-filling properties of the orthogonal array-based Latin Hypercubes.
The LHS is a well-established procedure for defining pseudo-random designs of samples with good properties with respect to the uniform probability distribution on an hypercube [0, 1] d . In particular, if compared to a random sampling, LHS: (1) https://doi.org/10.5194/nhess-2020-382 Preprint. Discussion started: 21 November 2020 c Author(s) 2020. CC BY 4.0 License. enhances the capability to fill the space; (2) avoids the overlapping of point locations projections; and (3) reduces the dependence on dimensionality (Stein, 1987;McKay et al., 2000;Ranjan and Spencer, 2014). 540 Once the desired number of samples is selected, [0,1] is divided in equal bins, then each bin will contain one and only one projection of the samples over every coordinate. There is a large number of possible designs, i.e. the number of permutations of the bins in the -projections. If these permutations are randomly sampled, clusters of points or regions of void space may be observed. For this reason, we base our design on the orthogonal arrays, which constrain the samples to fill the space with respect to a regular grid at a coarser scale than the LHS bins (Owen, 1992;Tang, 1993;Ai et al., 2016;Patra et 545 al., 2018;Bevilacqua et al., 2019b;Patra et al., 2020).
In this study we further modified a two-dimensional LHS to work under the assumption of a non-uniform PDF, i.e. the vent opening probability map (·). This approach was also implemented in Bevilacqua et al. (2019a). In particular, for any given y we defined the conditional PDF along x-direction: where is the appropriate normalizing constant.
Then we calculated the marginal PDF along y-direction: where ( ) and ( ) are the limits of the mapping region along the x-direction for any given .
Thus, once sampled a uniformly distributed LHS, we apply the transformation (·) to obtain a set of input variables that are distributed according to the vent opening map (·), but preserves the good space-filling properties of the original LHS.

Author contribution 560
All co-authors contributed in the compilation and critical revision of volcanological data of SSVC and NCVC. AB, RC and AN contributed in the choice of the different strategies and assumptions adopted to construct vent opening probability maps.
AB constructed the presented vent opening probability maps. AA performed the branching energy cone simulations and then processed the numerical results. AB, AA and RC prepared the manuscript with contributions from all co-authors.

Competing interests 565
The co-authors do not have any competing interests.