Topographic Controls on Pyroclastic Density Current Hazard at Aluto Volcano (Ethiopia) Identified Using a Novel Zero‐Censored Gaussian Process Emulator

Aluto volcano (Central Ethiopia) displays a complex, hybrid topography, combining elements typical of caldera systems (e.g., a central, flat caldera floor) and stratovolcanoes (e.g., relatively high and steep, radial flanks, related to eruptions occurring clustered in space). The most recent known eruptions at Aluto have commonly generated column‐collapse pyroclastic density currents (PDCs), a hazardous phenomenon that can pose a significant risk to inhabited areas on and around the volcano. In order to analyze and quantify the role that Aluto's complex topography has on PDC hazard, we apply a versatile probabilistic strategy, which merges the TITAN2D model for PDCs with a novel zero‐censored Gaussian Process (zGP) emulator, enabling robust uncertainty quantification at tractable computational costs. Results from our analyses indicate a critical role of the eruptive vent location, but also highlight a complex interplay between the topography and PDC volume and mobility. The relative importance of each factor reciprocally depends on the other factors. Thus, large PDCs (≥0.1–0.5 km3) can diminish the influence of topography over proximal regions of flow propagation, but PDCs respond to large‐ and small‐scale topographic features over medial to distal areas, and the zGP captures processes like PDC channelization and overbanking. The novel zGP can be applied to other PDC models and can enable specific investigations of PDC dynamics, topographic interactions, and PDC hazard at many volcanic systems worldwide. Potentially, it could also be used during volcanic crises, when time constraints usually only permit computation of scenario‐based hazard assessments.


Introduction
Pyroclastic density currents (PDCs) are hot, gravity-driven mixtures of volcanic particles and gas that are generated by a variety of explosive (e.g., eruption-column collapse) and non-explosive (e.g., lava-dome collapse)

10.1029/2023JB028645
Key Points: • A novel zero-censored Gaussian Process emulator is systematically applied to explore pyroclastic density current hazard at Aluto (Ethiopia) • A complex, multi-scale interplay between flow volume, flow mobility and Aluto's varied topography arises from our quantitative analyses • Vent location is the primary control on hazard but flow-source forcing and flow mobility can also outweigh the other controlling factors
In this study, we propose an innovative and versatile probabilistic strategy to flexibly investigate the links between topography and PDC hazard.The novelty of the strategy is twofold.First, we apply a multi-model framework (after Clarke, 2020;Clarke et al., 2020) to parameterize the TITAN2D model (Patra et al., 2005;Pitman et al., 2003) using flux sources, which better approximate a column-collapse mechanism for the generation of PDCs, relevant for Aluto and many other volcanic systems (e.g., Rutarindwa et al., 2019;Tennant et al., 2021;Tierz et al., 2018).Second, we systematically apply a novel zero-censored Gaussian Process (zGP) emulator (Spiller et al., 2023) to explore the links between topography and PDC hazard.The zGP emulators are built using sets of TITAN2D simulations and act as a substitute (or surrogate model, see Section 2.2) for the simulator.They enable efficient quantification of the aleatory and epistemic uncertainties relevant for volcanic hazard assessment in a way that would be computationally unfeasible using TITAN2D alone (e.g., Bayarri et al., 2009;Spiller et al., 2014).The zGP emulators also represent an important advance in dealing with the "zero problem" in hazard assessment: a collection of issues including the fact that model outputs tend to be strictly positive or that "near-hit" zero values are indistinguishable from those zero values far away from the zero-/ positive-output envelope, within the model parameter space (Spiller et al., 2023).We use our novel TITAN2D-zGP strategy to quantify uncertainty and explore the influence of topography on PDC hazard.We focus on characterizing the aleatory uncertainty (e.g., C. B. Connor et al., 2003;Marzocchi et al., 2004), linked with the natural variability in PDC events (e.g., vent location, PDC volume, mobility, etc.).Nonetheless, we also incorporate some epistemic uncertainty (e.g., Marzocchi et al., 2004;Marzocchi et al., 2008), related to incomplete knowledge, into our analysis, similarly to previous, doubly-stochastic probabilistic hazard assessments for PDCs (e.g., Bevilacqua, Neri, et al., 2017;Neri et al., 2015;Rutarindwa et al., 2019;Sandri et al., 2012;Sandri et al., 2018;Spiller et al., 2014;Tierz et al., 2018;Tierz, Sandri, Costa, Sulpizio, et al., 2016).
We present an application of our strategy to Aluto volcano, in Central Ethiopia (Figure 1).The volcano hosts a geothermal power station, and a diffuse but significant rural population (Central Statistical Agency, 2007;Clarke, 2020;Vye-Brown et al., 2016), as well as important agricultural and economic resources for the region.In volcanological terms, Aluto is a peralkaline rhyolitic caldera whose most recent eruptive history (at least since 60 ka) is associated with around 100 known eruptive vents, most of which are interpreted to relate to pumice-conebuilding eruptions.These enigmatic eruptions commonly display both explosive (including PDCs) and effusive phenomena (e.g., Clarke, 2020;Clarke et al., 2019Clarke et al., , 2020;;Di Paola, 1972;Fontijn et al., 2018;Hutchison et al., 2016;Kebede et al., 1985).The PDCs known from the deposits of these eruptions (e.g., Clarke, 2020;Clarke et al., 2020;Fontijn et al., 2018;Hutchison et al., 2016): (a) developed concentrated (dense) lower-flow boundary zones (Branney & Kokelaar, 2002); (b) were generated by the collapse of unsteady eruption columns, and (c) generally, tended to be valley-confined.That is, the flows tended to be strongly influenced by the topography of the volcano.Aluto possesses a complicated and varied topography (Figure 1), as it combines elements of caldera systems (e.g., a topographic depression, or caldera floor, formed during the latest caldera-forming eruption(s), 300-320 Ka, Hutchison et al., 2016;Fontijn et al., 2018) and stratovolcanoes (e.g., relatively long and steep flanks, which dip both outward and inward, toward the caldera floor, Figure 1).These positive reliefs at Aluto are linked to the eruptive products from multitude of overlapping eruptive centers that form cones, ramparts and silicic lavas, and collectively obscure the underlying caldera-rim faults (Clarke, 2020;Clarke et al., 2020;Hunt et al., 2019;Hutchison et al., 2015Hutchison et al., , 2016)).The volcano is also variously incised by narrow gorges and wider valleys, forming drainages of many scales and in all directions.The co-existence of these features in one volcanic system, though not unique, (a) provides an excellent test-bed to explore the influence of topography on PDC propagation, and (b) necessitates particular consideration for PDC hazard, which is so strongly influenced by topography.Previous probabilistic hazard analyses for PDCs at Aluto have focused on hazard domains covering the entire volcanic edifice and surrounding areas (Clarke et al., 2020).In this study, we seek a complementary analysis and focus on specific points along the main roads that traverse the volcanic edifice, from the lower sectors of the volcano flanks up to the caldera floor and several points on the caldera rim (Figure 1).Given their critical role for communication, transportation, and as potential evacuation routes, roads and bridges are key elements to daily life at/around volcanoes.They can also dynamically affect exposure to volcanic hazard (e.g., if rendered non-functional by PDC impact, erosion or deposition), hence controlling changes in volcanic risk with time (e.g., Charbonnier et al., 2023;Hayes et al., 2022;Loughlin, Calder, et al., 2002;Mossoux et al., 2019;Osman et al., 2019).
We describe the PDC model (TITAN2D) and the uncertainty-quantification technique (zGP emulator), as well as new metrics to assess the influence of topography on PDC hazard, in Section 2. We then present our findings for Aluto, focusing on different spatial scales as a function of vent proximity, in Section 3. Finally, we discuss the links between topography and PDC inundation in a general context, and the use of emulators in such assessments, in Sections 4 and 5.

Methods
Our general modeling strategy (Figure 2) aims to obtain probabilistic estimates of PDC inundation; a quantification of PDC hazard, at a relatively low computational expense (e.g., Bayarri et al., 2009;Dalbey et al., 2008).To achieve our goal, we use an appropriate combination of a PDC model, herein also referred to as a simulator (TITAN2D, Patra et al., 2005) (see Section 2.1) and an uncertainty-quantification technique (zero-censored Gaussian Process, zGP, emulators, Spiller et al., 2023) (see Section 2.2).Several other strategies, differing in the choice of PDC model and uncertainty-quantification technique, have been previously proposed (e.g., Bayarri et al., 2009;Bevilacqua et al., 2021;Bevilacqua, Neri, et al., 2017;Charbonnier et al., 2020;Dalbey et al., 2008;Neri et al., 2015;Rutarindwa et al., 2019;Sandri et al., 2018;Spiller et al., 2014Spiller et al., , 2020;;Tierz et al., 2018;Tierz, Sandri, Costa, Sulpizio, et al., 2016;Tierz, Sandri, Costa, Zaccarelli, et al., 2016), albeit so far without the emergence of a standard approach (Tierz et al., 2021).This is partly related to the fact that each PDC model may be best suited to specific flow conditions and/or types (e.g., "dense" vs. "dilute" flows).Additionally, the most complex PDC models are too costly, in computational terms, to enable uncertainty quantification, even with the most-efficient techniques.zGP emulators can be seen as one of such highly-efficient uncertainty-quantification techniques.Our study represents the first-ever systematic application of zGP emulators to probabilistic hazard assessment of PDCs.This is an important step forward as key hazard variables related to PDCs (flow thickness and speed, area of PDC invasion, maximum runout, etc.) are bounded to be strictly positive, which represents a zero-problem that previous approaches had to deal with following more ad-hoc procedures (e.g., Rutarindwa et al., 2019;Spiller et al., 2014Spiller et al., , 2020)).
In this work, we focus on the three main road tracks that go from the lower northwest, northeast and southeast flanks of Aluto's edifice into the caldera floor, in the central area of the volcano (Figure 1).The roads provide the main foot and vehicular access and egress to settlements, water resources, and geothermal infrastructure on the volcano.They, therefore, represent great social and economic importance, as well as potential evacuation routes during a volcanic crisis.We discretize the roads into individual "road points" spaced 550 m apart, along the road.At each of these road points, plus the Aluto-Langano geothermal power plant (45 points in total), one zGP emulator is built, based on sets of TITAN2D simulations (see Sections 2.1 and 2.2).  , c) selected examples of TITAN2D simulations run at Aluto, which inform about some of the topographic effects discussed in this study.Please note that inset (a) is an approximate representation of the source conditions, and flow propagation, as simulated by TITAN2D, where only the dense, basal part of PDCs is actually modeled (e.g., Tierz et al., 2018).Flow volume and bed friction angle used in the simulations shown in (b, c) are indicated on the upper part of the maps.Blue lines correspond to 0.1 m flowthickness contours.The other contours correspond to 0.2 intervals in ln(h + 1), starting in 2.4, where h is the flow thickness.Labels show flow thickness in meters.

TITAN2D: Model and Parameterization
The TITAN2D model (Patra et al., 2005;Pitman et al., 2003) simulates dry granular flows by using a "shallowwater" approximation (i.e., flow vertical extent is significantly smaller than its lateral extent and vertical velocity component can be considered negligible) for a continuum that propagates over realistic topography, using a Digital Elevation Model (DEM), gaining momentum due to the action of gravity and losing it mainly due to frictional forces.For a complete description of TITAN2D governing equations, see Pitman et al. (2003); Patra et al. (2005).For a comparison with the equations of other "shallow-water" PDC models, see Gueugneau et al. (2021).Two advantages of using TITAN2D are (a) the open-source nature of the numerical code, and (b) the relatively short runtimes (minutes to tens of hours) given the complexity of the physical model, which are achieved through code parallelization and the use of an adaptive computational grid (Patra et al., 2005).In terms of frictional forces, several rheological models are available to describe how these affect flow propagation in TITAN2D v4.0.0, the version that we use in this study (https://vhub.org/resources/4057/download/Titan2D_User_Guide.pdf).We use the Coulomb rheological model which has been successfully applied to several other volcanoes worldwide (e.g., Charbonnier & Gertisser, 2012;Spiller et al., 2014;Tennant et al., 2021;Tierz et al., 2018;Vázquez et al., 2019).
One important element of our modeling strategy (Figure 2a) is the use of flux sources to simulate PDC initiation, which better approximates the inferred PDC generation mechanism (column collapse) during pumice-conebuilding eruptions at Aluto volcano (Clarke, 2020;Clarke et al., 2020).Flux sources have been seldom used in TITAN2D (e.g., Charbonnier & Gertisser, 2012;Ogburn & Calder, 2017;Tennant et al., 2021) and never in the context of PVHA of PDCs (e.g., Bayarri et al., 2009;Rutarindwa et al., 2019;Spiller et al., 2014;Tierz et al., 2018).We adopt an innovative, multi-model approach to parameterize these flux sources (after Clarke, 2020;Clarke et al., 2020), based on three TITAN2D parameters: flux (vent) radius, flux rate per unit area (hereinafter named just "flux rate"), and flux duration.In order to constrain appropriate values for these parameters at Aluto volcano, we follow the modeling strategy described in Clarke et al. (2020).In summary, we use the results from running the eruption-column-dynamics model PlumeRise (Woodhouse et al., 2013) to build a multi-variate regression model that provides estimates for the Mass Eruption Rate (MER) at the onset of total column collapse.This model is built using two independent variables only: vent radius and magma gas-mass fraction (after Clarke, 2020;Clarke et al., 2020).Using appropriate values for the eruptive-mixture density, controlled by the magma gas-mass fraction (Clarke, 2020), we can derive values of flux rate, the TITAN2D parameter, which are compatible with the MER at the onset of column collapse.Our parameterization explores vent radii between 1 and 150 m, based on observations at Aluto and a "precautionary principle" of expanding the range toward larger events (Clarke, 2020;Clarke et al., 2020), and compatible flux rates (per unit area) of 20-180 m/s, according to the procedure explained above.This implies a calculated, approximate range of PDC volumes between 0.004 and 763 Mm 3 , consistent with the available field constraints (Clarke, 2020).
To parameterize the flux-source duration, we choose a somewhat arbitrary but informed value of 240 s.Three aspects motivate this choice: (a) shorter flux durations would eventually approximate the instantaneous release of pile sources and could generate unrealistically-high collapse heights (Ogburn & Calder, 2017;Rutarindwa et al., 2019); (b) longer flux durations could result in unrealistically-low collapse heights (relative to their mass eruption rates, Clarke et al., 2020;Clarke, 2020), and would also result in overly long runtimes, as the stopping time increases with the duration of the flux source (e.g., Charbonnier & Gertisser, 2012;Tierz et al., 2018;Yu et al., 2009); (c) the chosen value of flux duration generates collapse heights that are compatible across PDC models, in this case TITAN2D and the Energy Cone (EC) model (after Clarke, 2020;Clarke et al., 2020).Finally, in terms of the TITAN2D stopping time (e.g., Dalbey, 2009;Yu et al., 2009), we choose a value of 400 s, based on the selected flux duration of 240 s, and on the observed dependence of the stopping time on both flux duration and bed friction angle (e.g., Charbonnier & Gertisser, 2012;Sulpizio et al., 2010;Tierz et al., 2018).
Apart from vent radius and flux rate, another three TITAN2D parameters are modeled as uncertain in our simulation-emulation strategy (Figure 2a): bed friction angle, and UTMx-y coordinates of the simulated vent location (flux source).We explore a range between 6 and 30°for the bed friction angle, which is compatible with: (a) previous studies combining TITAN2D and GP emulators (e.g., Rutarindwa et al., 2019;Spiller et al., 2014); and (b) previous studies analyzing PDC mobility at Aluto volcano (Clarke et al., 2020).Finally, in terms of vent locations, we explore a large area covering the following coordinates (EPSG:32637-WGS 84/UTM zone 37N): UTMx [468,000, 484,000] m; and UTMy [852,490, 872,000] m.These limits are also fully compatible with what is known in terms of vent-opening patterns and potential future vent locations at Aluto volcano (Clarke, 2020;Clarke et al., 2020;Hutchison et al., 2015;Tierz et al., 2020).

Zero-Censored Gaussian Process Emulator: Model and Parameterization
Surrogate modeling can be defined as the ensemble of statistical and computational procedures aimed at obtaining approximations, or representations, of a computationally-expensive computer model (simulator), that are reliable enough as to be used as substitutes, or surrogates, of this simulator.One realm of surrogate models, also termed statistical emulators, are those that use stochastic processes to quantitatively capture key relationships between the computer model parameters (inputs) and model predictions (outputs) (e.g., O'Hagan, 1978;Razavi et al., 2012;Sacks et al., 1989;Smith, 2013).Through the emulator, users can obtain appropriate estimates of output values, together with a measure of epistemic uncertainty, at untried locations within the simulator input parameter space.This helps save a significant amount of computational time and resources, hence enabling procedures such as forward uncertainty quantification or Monte Carlo (MC) calculations of hazard probabilities that would otherwise be intractable using the simulator alone, through classic forward simulation (e.g., Dalbey, 2009;Smith, 2013;Spiller et al., 2014).
Gaussian Process (GP) emulators are one kind of surrogate model that has become increasingly relevant in, and widely applied to, a large variety of scientific problems in different fields and areas of research (e.g., Baker et al., 2022;Bayarri et al., 2009;Beck & Guillas, 2016).In our case, we consider an input to TITAN2D to be a set of values x = {x 1 , x 2 , x 3 , x 4 , x 5 }, where x 1 = flux (vent) radius, x 2 = flux rate, x 3 = bed friction angle, x 4 = UTMx (Easting), x 5 = UTMy (Northing), all normalized to the ranges indicated in Section 2.1.Note that x 4 , x 5 define the location of the eruptive vent.For each input x, the output y is taken to be the maximum flow thickness (h) at a particular road point.More precisely, acknowledging that maximum flow thicknesses can vary form hundreds of meters down to less that a meter, we take y as the ln(h + 1).Considering one road point, we will let T (D for design) be a vector of simulated output corresponding to each of N TITAN2D runs at scenarios in the design X D (NB.We provide details on the selected design for Aluto later in this section).That is, Then, the GP emulator is given by: where Z(•) is a constant-variance, mean-zero Gaussian process and μ(•) is a user-specified trend function, often chosen to be constant or linear.In this application, we take the trend function to be linear in the vent radius, flux rate and distance between road point and vent location, and constant in the bed friction angle.As with the trend function, the choice of correlation functions is user-specified and here we use the Matèrn 5-2 (Stein, 1999), given by: where ρ k is the correlation length in the kth dimension of input space and x ik , x jk are the ith, jth design point values in the kth dimension.The . The predictive mean and predictive variance at an untested scenario in input space, x, are given by: where 1 is an N-vector of ones and r ) is an N-vector of the correlation between the untested scenario and each of the design points.The parameters in these definitions-the scalar variance, σ 2 z , parameters in μ(⋅), and parameters in the correlation structure, ρ 1 , …, ρ 5 -need to be estimated and we do so using the RobustGaSP-in-Matlab package (https://github.com/MengyangGu/RobustGaSP-in-Matlab;Gu & Berger, 2016;Gu et al., 2018;Gu, 2019).The semi-binary nature (many design scenarios lead to no-flow) of the TITAN2D data poses a significant obstacle to fitting GPs as the Gaussian processes employed are: (a) assumed to be stationary, and (b) can take on any value, positive or negative.
The zero-censored GP, or zGP (Spiller et al., 2023), overcomes this challenge by imputing negative y values for design points that resulted in no flow.The basic idea is to impute negative y values by sampling a GP fit to the subset of positive-response values, and a few selected zero-response values, and keeping only negative samples.In practice, the algorithm is effectively a Gibbs sampler that cycles through each design point leading to a zero response, and samples a negative replacement value from the GP fit to all of the other design points, including the other design points with imputed negative responses (Spiller et al., 2023).Once the process is complete, the vector y D z of the original positive and newly imputed negative responses is used in place of the original y D for the predictive mean in Equation 3 above.Then, the predicted output at an untested point, x, is taken to be the max (0, h(x)).Now, by evaluating h( ⋅ ) at scenarios where TITAN2D was not run, for example, x ∈ D, but x ∉ X D , we have faithful estimates that adhere to the semi-binary nature of TITAN2D.Further, s 2 (•) quantifies how much uncertainty we introduce by replacing h(•), the true mean of the Gaussian process, with h(•), its estimator.We stress that, even though the zGP is able to emulate the TITAN2D outputs, in situations where the simulator performance is an inadequate model, the surrogate will inherit the simulator's shortcomings.
For our TITAN2D-zGP design at Aluto, we run a total of 4,200 TITAN2D simulations, following Latin Hypercube Sampling designs aimed at effectively covering the model parameter space.The total runtime of such a sample size of simulations, using High-Performance-Computing (HPC) resources, might be around 2-4 weeks, depending on the availability of the HPC resources.To build emulators specific of each of the road points (cf.Rutarindwa et al., 2019;Spiller et al., 2014), we sub-sample: (a) all simulations that result in PDC inundation at the road point (i.e.positive output here taken as h > 0.1 m); and (b) 250 simulations that result in zero output at the road point.The number of zero-output design points used is similar to the number of positive-output design points, and chosen to strike a balance between: (a) the zGP emulators having enough design points at both sides of the zero-flow isoline as to capture it successfully (Figures 3b and 3c); and (b) the zGP construction not being overburdened, computationally, by having too many design points, which can compromise the inversion of large matrices required in the building of the emulators (Spiller et al., 2014(Spiller et al., , 2023)).Regarding (a), it is important to note that, depending on the selection of design points, as well as on the use of either zGP emulators or classical GP emulators, the emulation results may be different (Figures 3b and 3c).One recurring issue that can be observed is that of emulator surfaces that cannot properly capture the flat nature of the response across areas of the parameter space dominated by zero-flow outputs (Figure 3b).At Aluto, emulator surfaces analyzed over the UTMx-y parameters (having fixed the other three parameters) tend to show either large areas of positive predictions over zero-flow domains (Figure 3c, top) and/or "rebounds" from zero or negative predictions to positive ones away from the zero-flow isoline (Figure 3c, middle).The use of the zGP emulator, and a careful selection of the 250 zero-output design points (cf.Spiller et al., 2014;Rutarindwa et al., 2019, see also Supporting Information S1), help address the aforementioned issues.The zGP models in this study (after Spiller et al., 2023;Tierz, 2023) are built in Matlab (MATLAB, 2022), and include the use of the Robust GaSP Matlab package (https://github.com/MengyangGu/RobustGaSP-in-Matlab;Gu & Berger, 2016;Gu et al., 2018;Gu, 2019).

Exploring Topographic Effects on PDC Hazard
We focus on three spatial scales of analysis depending on where the emulated eruptive vents are located with respect to the road point analyzed: (a) proximal (i.e., vents located about 0-2 km away from the road point, Figure 4a), (b) medial (i.e., vents located about 2-5 km away from the road point, Figures 4a, 4c, and 4d), and (c) distal (i.e., vents located up to about 10 km away from the road point).
To analyze the proximal topographic effects, we design vent-opening square regions centered at each of the road points, with increasing side lengths of 1, 2 and 3 km (maximum distance vent-road point ≈2.1 km).The emulated PDCs are sourced from eruptive vents sampled uniformly within these vent-opening squares (Figure 4a).We emulate 10 5 PDCs for each vent-opening square size and road point (excluding the geothermal power plant), and calculate the probability of PDC inundation as the frequency of events that result in PDC inundation at the road point.To explore the influence of PDC size (volume) and mobility (bed friction angle), we carry out our analyses using three sets of emulated PDC events: (a) any volume and bed friction; (b) small PDC events (i.e., lower-thanaverage volume, higher-than-average bed friction); and (c) large PDC events (i.e., higher-than-average volume, lower-than-average bed friction).Here, "average events" refer to V PDC = 0.01 km 3 and ϕ bed = 15°, based on data from PDCs analog to Aluto's, as presented in Clarke et al. (2020).
Further, to assess the influence of proximal topographic features on PDC inundation (thus, hazard), we design two simple DEM-based metrics (Figure 4b): (a) the difference between the mean altitude across the vent-opening square and the altitude at the road point; and (b) a metric (Γ jk ) informing about the degree to which the topography within the kth vent-opening square deviates from a funnel shape centered at the jth road point (Figure 4b).This metric is calculated as follows: (5) where N jk is he total number of DEM points within the kth vent-opening square, and γ i is the smallest angle between the vector of terrain aspect (QGIS Development Team, 2023), at each of the ith DEM points inside the kth vent-opening square, and a line connecting each of such points with the jth the road point (e.g., road point 7 in Figure 4b).This angle takes on values between 0°(terrain aspect directly points toward the road point) and 180°( terrain aspect points in the opposite direction of the road point).If we calculate Γ jk /(180°• N jk ), we obtain a normalized metric, Γ n jk , indicating how much the topography within the kth vent-opening square deviates from a funnel shape centered at the jth road point.
To analyze the medial topographic effects, we select two vent-opening areas located close to or along the southern caldera rim of Aluto volcano (Figures 4a, 4c, and 4d).We sample vent locations uniformly from within these areas (10 5 times).These vent-opening areas correspond to high vent-opening probability values according to the model presented in Clarke et al. (2020), based on the known locations of past eruptive vents at Aluto volcano.The road points that we analyze are located in medial regions away from these vent-opening areas.In particular, the road points are located over the western sector of the caldera floor (medial zone 1, Figure 4c) and over the eastern sector of the caldera floor, including the Aluto-Langano geothermal power plant, as well as the lower southeastern flank of Aluto volcano (medial zone 2, Figure 4d).We note that, in the design of the zGP emulators at each of the road points, we include all zero-flow simulations whose vent locations fall within the limits of the corresponding vent-opening area.This aims at providing the zGP emulators with further constraints to successfully carry out the detailed analyses that we present for the medial topographic effects.
Differently to the case of proximal topographic effects, one focus of the medial topographic analyses is on emulated PDC events with vent locations inside the same vent-opening area and that may inundate several road points simultaneously.A main goal is to explore, in probabilistic terms, the effects of processes such as PDC channelization and overbanking, which can be critical for PDC hazard in medial to distal areas (e.g., Charbonnier et al., 2023;Charbonnier & Gertisser, 2012;Gueugneau et al., 2021;Ogburn & Calder, 2017;Ogburn et al., 2014).Given the key influence of the bed friction angle on medial to distal PDC propagation using the TITAN2D model (e.g., Charbonnier & Gertisser, 2012;Sheridan et al., 2005;Spiller et al., 2014;Sulpizio et al., 2010;Tierz et al., 2018), we focus on this variable while exploring the medial topographic effects at Aluto volcano.In particular, we randomly divide the 10 5 sampled vent locations into 25 equal-size groups of bed friction angles (at intervals of 1°from ϕ bed = 6°to ϕ bed = 30°) and evaluate the zGP emulators at those fixed values of bed friction.In terms of PDC volume, we fix the values of vent radius and flux rate (hence, fixing the volume) for all 10 5 emulator evaluations.We repeat the same procedure using three different PDC volumes around 0.1, 0.3 and 0.6 km 3 .We focus on flows with volumes above 0.1 km 3 , as small-size PDCs infrequently propagate across the relatively long distances between vent-opening areas and road points that are explored in the medial topographic analyses.
To analyze the distal topographic effects, we explore vent locations within the ≈16 × 18 km vent-opening domain in Clarke et al. (2020) (Figure 4a).For the other TITAN2D parameters, we explore their full ranges (Section 2.1).We define probability density functions (pdfs) for all five uncertain parameters, sample these distributions (10 5 times) and use those samples to evaluate the zGP emulators built at each of the road points at Aluto, propagating the aleatory uncertainty from model parameters into model outputs.In other words, our distal topographic analyses are, in practice, an illustrative example of PVHA for dense PDCs at Aluto volcano, using the TITAN2D-zGP strategy.Similar to previous studies using GP emulators, we incorporate an inverse relationship between the PDC volume and the bed friction angle into our model of aleatory uncertainty (e.g., Ogburn et al., 2016;Rutarindwa et al., 2019;Spiller et al., 2014).Our simple model is based on fitting a log 10 -log 10 linear regression to the aforementioned data set of analog PDCs (pumice flows) to those known at Aluto volcano (Clarke et al., 2020): where V PDC denotes PDC volume (in m 3 ), and ΔH/L denotes the effective friction coefficient, a proxy for flow mobility, which we relate to the bed friction angle as follows: ϕ bed ≈ tan 1 (ΔH/L) (e.g., Charbonnier & Gertisser, 2012;Ogburn et al., 2016;Spiller et al., 2014;Tierz et al., 2018).We also explore the influence of epistemic uncertainty by using alternative aleatory-uncertainty models.one but without incorporating the relationship between PDC volume and bed friction (cf.Dalbey et al., 2008;Tierz et al., 2018).In that case, bed friction angles are sampled from a truncated Gaussian pdf with the same mean and variance as the one used in Clarke et al. (2020), but truncated at a higher minimum bed friction angle of 6°(instead of 2°), as TITAN2D cannot deal with such very-high-mobility flows (e.g., Ogburn et al., 2014;Procter et al., 2010;Tierz et al., 2018).The last model (all-Uniform model) is a purely non-informative model, where all five uncertain TITAN2D parameters are modeled using Uniform pdfs (please note that the PDC volume is not distributed uniformly, as it is a non-linear transformation of vent radius and flux rate).The all-Uniform model allows us to better understand the impact of having a full-aleatory-uncertainty model to quantify PDC hazard (mainly based on Clarke, 2020;Clarke et al., 2020) versus the situation of not having any such model available (see Sections 3 and 4).

Results
Figures 2b and 2c show four examples of TITAN2D simulations run in this study.They illustrate PDCs of different sizes and mobility, sourced from different vent locations.Simulation 2b1 shows a large flow (≈0.5 km 3 ), with intermediate mobility (ϕ bed = 15°), that initiates on the upper areas of the catchment explored in medial zone 1 and propagates, more or less radially, but mainly toward the northeast.Down the main SW-NE valley, the flow produces extensive overbanking, with some channelization observed over proximal areas only.Simulation 2b2 corresponds to a moderate-volume flow (≈0.15 km 3 ), with relatively high mobility (ϕ bed = 11°), that originates from the E-W ridge that marks the SE sector of the caldera rim and mainly propagates southward, following the wide valley shown in profile c-c' (Figure 1c).The flow tends to follow the bend located on lower sectors of the valley (profile d-d', Figure 1c), although it also produces some overbanking over the road points on the left bank of the valley (Figure 2b).Simulation 2c1 shows a moderate-volume flow (≈0.1 km 3 ), with high mobility (ϕ bed = 9°), sourced from a similar location to simulation 2b1 but that mainly follows the catchment area and main SW-NE valley toward the caldera floor.Overbanking of the flow does occur, particularly after the break-in-slope, which coincides with a narrowing and deepening portion of the valley, but it is relatively limited.Finally, simulation 2c2 displays a moderate-to-large flow (≈0.4 km 3 ), with high mobility (ϕ bed = 8°) and a vent location slightly west of simulation 2b2, that propagates both toward the southern, outer flank of the volcano, but also northward, toward the inner flank of the volcano, onto the caldera floor (Figure 2c).Toward the south, the flow follows the main valley, and responds to the bend in the valley, but it also produces considerably greater overbanking that in the case of simulation 2b2 (Figures 2b and 2c).Propagation over the flat area beyond the break-inslope of the outer flank, toward the W-SW, is also enhanced.Toward the north, the flow inundates three main valleys: one going NW, another going NE (also inundated in simulation 2b2), and the central main valley going N, toward the geothermal power plant (profile e-e'-e", Figure 1c).PDCs that propagate down these three valleys converge over the caldera floor on its central to eastern sector, and cause extensive PDC inundation over the flat terrain (Figure 2c).
Figure 3d displays the values of correlation lengths (ρ k ) for each TITAN2D parameter and each zGP emulator built at each of the road points analyzed in our study.We note that the correlation lengths are solely based on the ensemble of TITAN2D simulations run according to the design described in Section 2.2, and do not change depending on the choice of parameter values used to carry out the zGP evaluations described across Sections 3.1-3.3.The correlation lengths describe the influence of each TITAN2D input parameter on the model output (maximum flow thickness): the more influential the parameter, the higher the correlation length.Also, changes over higher values are more relevant (please notice the logarithmic scale).The results highlight how the vent location (UTMx-y parameters) is critical to quantify PDC inundation (hence, hazard) at Aluto volcano, mainly due to the spatially-varied topography of the volcano (Figure 1).These results further justify the focus of our analysis on the role of topography.It is also worth noting the greater influence of the UTMy (north-south) parameter, in general, and in particular around road points 18-30, located the closest to the main E-W ridge on the S-SE sector of the caldera rim, the highest area of the volcanic edifice (Figure 1).The other three TITAN2D parameters are less influential than the vent location, with decreasing average influence from vent radius to bed friction to flux rate.PDC volume depends on the square of the vent radius, and the latter also marks the lateral extent of the source term, hence influencing the very-proximal regions of flow propagation.We observe that the correlation lengths for vent radius and flux rate have more homogeneous values across road points compared to those for the bed friction angle.This may hint to the complexity of the relationship between bed friction and flow thickness (or invasion area), which again justifies the use of statistical emulators to explore such influence in detail, without having to run impractical numbers of TITAN2D simulations (e.g., Rutarindwa et al., 2019;Spiller et al., 2014).

Proximal Topographic Effects
We analyze the data for all the road points together, focusing on the two simple morphological metrics described in Section 2.3: (a) the difference between the mean altitude across the vent-opening area (1, 2, 3 km squares centered at each road point) and the altitude at the road point; and (b) the normalized morphological deviation of the vent-opening area from a funnel shape centered at each road point (Figure 4b).We compare these two metrics with the probability of PDC inundation at each road point for small, large and all PDCs (see Section 2). Figure 5 shows that, over very proximal spatial scales (i.e., vent-opening squares of 1 km side, Figures 5a and 5b), large PDCs are associated with probabilities of PDC inundation close to 1, independently of the specific proximal topography surrounding each road point (Figure 5a).For small PDCs, the probabilities of PDC inundation show random patterns (influenced by specific combinations of vent radius, flux rate, bed friction and vent location), when compared to the morphological metrics (Figure 5b).On the other hand, for vent-opening squares of 3 km side, we observe that for both large and all PDCs, the greater the altitude difference between the vent-opening area and road point, and the more funnel-shaped the vent-opening area, the higher the probability of PDC inundation computed at the specific road point (Figures 5c and 5d).

Medial Topographic Effects
In Figure 6, we display a summary of the changes observed in the mean emulated flow thicknesses across the three sets of road points (columns in Figure 6) analyzed for the medial topographic analyses.We also explore PDCs of different volumes (rows in Figure 6) and bed friction angles.
Under a morphological point of view (Figures 1 and 4), medial zone 1 is meant to approximately cover vent locations whose PDCs will tend to be channelized down one main valley with SW-NE direction, connecting the elevated areas on the SW sector of Aluto's volcanic edifice with the western sector of the caldera floor (Figures 1b  and 4c): from east to west, road points 35-39.A clear break-in-slope in the topography marks the junction between the elevated areas and the caldera floor, with only a narrow valley continuing toward point 36, and 35 eastward (Figure 4c).Road points 37 and 38 are "overbank locations" north of the narrow valley and are located in front of the main outlet of the catchment at the break-in-slope (Figure 4c).Road point 39 is located beyond a small topographic obstacle formed by past eruptive products (Figure 4c).Medial zone 2 corresponds with the highestaltitude area along the southern caldera rim of Aluto.PDC events sourced from this zone are expected to impact road points located on the eastern sector of the caldera floor (inner flank) as well as road points located on the lower SE flank of the volcano (outer flank, Figure 4d).The slopes connecting these highest-altitude areas with the low-altitude areas can be relatively steep on both sides (Figure 1c, profile b-b').Down the outer flank, one main, wide valley runs N-S for about 2 km before a bend marks a change to an ENE-WSW orientation of the valley (Figures 1c and 4d).From east to west, road point 26 is at an "overbank location" right out the valley bend, while road points 27 to 30 are at the valley's margin, running approximately parallel to it (Figure 4d).Down the inner flank, two main valleys connect the high, caldera-rim areas with the topographic low of the caldera floor (Figures 1c and 4d): from east to west, road point 15 is located beyond the exit of the E valley (Figure 1c), and road points 14, 2 (the geothermal power plant) and 32 are roughly located in front of the exit of the W valley (Figures 1c and 4).Road point 14 is slightly north of road point 2, and over a gently up-sloping sector of the caldera floor.
We observe that flow thicknesses are considerably lower for road points linked to medial zone 1 compared to those linked to medial zone 2. Among the road points down the inner and outer flanks, flow thicknesses are generally similar, apart from the case of PDC volumes ∼0.6 km 3 for which flow thicknesses are greater at the outer-flank road points (Figures 6h and 6i).In medial zone 1, the general pattern is that of decreasing flow thicknesses with increasing bed friction angles, for all road points and PDC volumes explored.At road point 36, flow thicknesses are largest for moderate PDC volumes (Figures 6a and 6d).At road point 37, flow thicknesses increase, comparatively to other road points, as PDC volume increases.For PDC volumes ∼0.6 km 3 , and values of bed friction angle ≥17°, flow thicknesses are greater at road points 37 and 38, on the overbanking regions, compared to those at road point 36 (Figure 6g).Elsewhere, flow thicknesses are the smallest at road points 35 and 39, respectively (Figures 6a, 6d and  6g).In medial zone 2, the road points located at the base of the inner flank, show relatively similar patterns to those linked to the medial zone 1.At road point 2, flow thicknesses are largest, independently of the flow volume considered (Figures 6c, 6f and 6i).At road point 32, flow thicknesses are generally second-largest, particularly for low bed friction angles.For PDC volumes ∼0.6 km 3 , and bed friction angles ≥15°, flow thicknesses at road point 15 are greater than at road point 32 (Figure 6i).At road point 14, flow thicknesses are lowest.Similar to what observed at road points down the medial zone 1, a change in slope in the flow-thickness curves tends to occur at bed friction angles around 15°, with the curves flattening out toward larger values (i.e., lower PDC mobility, Figure 6).On the road points down the outer flank, the patterns of the flow-thickness curves are slightly more complex.In general, flow thicknesses are largest at road point 26.At road point 27, flow thicknesses are second-largest, and show some relationship with those at road point 26, particularly for flow volumes ≥0.3 km 3 and bed friction angles ≥15°, approximately (Figures 6e and 6h).Hence, flow-thickness curves display a "hump" at road point 26, while they display a "trough" at road point 27, particularly for PDC volumes ∼0.6 km 3 .At road points 28 to 30, flow-thickness curves show similar values and patterns, with the lowest flow thicknesses generally happening at road 30, furthest away from the vent-opening area (Figures 6b, 6e and 6h).

Distal Topographic Effects
Figures 7a and 7b show two profiles of the probability of PDC inundation (given distal vent locations) at different road points across the caldera floor and the NE and SE inner/outer flanks of the volcano.A measure of epistemic uncertainty in PDC inundation (hence, hazard) is given by the length of the vertical lines, highlighting the differences in probability of PDC inundation between the full-aleatory-uncertainty and all-Uniform models (see Section 2).On the E-W profile across the caldera floor, the epistemic uncertainty is slightly higher over the edges of the caldera floor, compared to its central area.Probabilities of PDC inundation are generally higher for the fullaleatory-uncertainty model (blue lines) apart from the very central point of the caldera floor, where the all-Uniform model has a slightly higher probability (red line, Figure 7a).The presence of the Artu Jawe fault scarp across the caldera floor (Figure 4a) can be noted in the difference in probability between road points 3 and 15.On the N-S profile, the largest extent of the epistemic uncertainty occurs around road points 19-24 (Figure 7b).These points run between two main reliefs marking the caldera rim of Aluto (Figure 1), and these two elevated areas are associated with high vent-opening probabilities in the Clarke et al. ( 2020) vent-opening model.In some cases (red lines), the probability of PDC inundation is higher in the all-Uniform model than in the full-aleatory-uncertainty model.These all occur toward the lower sections of the outer flanks of the volcano on both the north and south sides (Figure 7b).In Figures 7a and 7b, we also display the probabilities of PDC inundation, at the nearest grid point, calculated using the Energy Cone-Monte Carlo (EC-MC) strategy of Clarke et al. (2020) (after Tierz, Sandri, Costa, Sulpizio, et al., 2016;Tierz, Sandri, Costa, Zaccarelli, et al., 2016).We observe that, in general, and especially in the N-S profile (Figure 7b), the patterns in these probabilities are compatible with what observed using the TITAN2D-zGP strategy.However, the values are always higher in the EC-MC strategy.The changes in probability of PDC inundation around the Artu Jawe fault scarp are considerably greater when using the latter strategy (Figures 7a and 7b).
In Figure 7c, we report curves showing the probability of PDC inundation (conditional on PDC volume) at the Aluto-Langano geothermal power plant, for the three aleatory-uncertainty models considered.We observe how the full-and partial-aleatory-uncertainty models yield practically identical curves, while the all-Uniform model displays lower probability values, particularly as the PDC volume increases (Figure 7c).In Figure 7d, we show a spatial comparison between the probabilities of PDC inundation using the EC-MC and TITAN2D-zGP strategies.
As shown by Figures 7a and 7b, the differences are strictly positive.However, we observe that these differences are higher over the caldera floor and lower sectors of the SE outer flank, and lower along the valley traversing the caldera rim on the SE sector of the volcano, and on the road points down the NW and NE flanks of the volcano (Figures 7a,7b,and 7d).

Discussion
The ensemble of results derived from our proximal, medial and distal analyses of the role of topography on PDC inundation at Aluto volcano (Figures 5-7) indicate a complex interplay between flow size (volume), flow mobility (bed friction angle) and the scale of the main topographic features that the flows encounter during their propagation (Figure 8).On the proximal scale, large flows are able to "force" high probabilities of PDC inundation (i.e., high hazard), independently of the general topography across their area of propagation (Figure 5a).Small PDCs are dependent on specific combinations of flow size, mobility and vent location (albeit proximal), and this also leads to a negligible influence of topography on PDC hazard (Figure 5b).As the modeled eruptive vents are located further away from the road point (Figures 5c and 5d), the influence of the PDC source conditions decreases in favor of a higher influence of topography, and the probabilities of PDC inundation can be linked with simple topographic parameters (Figures 5c and 5d).Our observations thus indicate that the "PDC-source forcing" (Figure 8) depends on the distance between vent location and target area (e.g., road point), as well as on the flow volume, at the proximal scale of analysis.The use of flux sources in TITAN2D can also influence the lateral extent of this "forced regime," meaning that conditions at the source, or very proximal areas (e.g., steep slopes), strongly influence flow propagation (Charbonnier & Gertisser, 2012;Doronzo, 2012;Ogburn & Calder, 2017).Pile sources, on the contrary, given their pseudo-instantaneous release of material, would create flows that transition into an "inertial regime," meaning flows that are minimally influenced by (forcing) conditions at the source (cf.Breard et al., 2023;Doronzo, 2012), more rapidly, and over shorter distances, for a given combination of flow size and mobility.Investigating such flow transitions following a systematic, ideally multi-model, approach (e.g., Gueugneau et al., 2021;Ogburn & Calder, 2017;Patra et al., 2020), and using probabilities of PDC inundation to fully capture uncertainty, could be a fruitful topic for future research, facilitated by the use of zGP emulators (Spiller et al., 2023).
Over medial scales, the links between PDC inundation, flow characteristics and topography become more complex, as flow mobility also starts playing a role (shaded area, Figure 8).PDC channelization and overbanking processes are key to analyze PDC hazard at Aluto.In medial zone 1, PDC channelization occurs down the main valley in the area (e.g., simulation 2c1, Figure 2c), but only if PDC volume is not large enough as for flows to completely fill and overflow that topographical feature, hence decreasing its role on the aforementioned processes.The latter happens for large (e.g., V PDC > 0.4-0.5 km 3 ) flows (e.g., simulation 2b1, Figure 2b), for which road point 37, and even 38, can show greater mean emulated flow thicknesses, especially for low-mobility flows (ϕ bed > 15-20°).Such flows tend to show more radial footprints, less influenced by topography (Figures 6g and  8).In medial zone 2, down the inner flank toward the caldera floor, PDC channelization down the main valley (western valley in Figure 1b) implies that mean emulated flow thicknesses are largest at road points 2 (geothermal power plant) and 32.PDC channelization down the eastern valley in Figure 1b (e.g., simulation 2b2, Figure 2b) contributes to greater mean emulated thicknesses at road point 15 compared to 14, or even compared to road point 32 for high-volume flows (Figures 6c,6f and 6i).At the road points down the outer flank, a competition between PDC channelization and overbanking can be clearly observed, in particular, at road points 26 and 27 (Figure 4d).Intermediate-to high-mobility flows (ϕ bed < 15°) tend to inundate both road points more readily, as the overbanking process is promoted and affects all road points analyzed (e.g., simulation 2c2, Figure 2c).However, for transitional flows (15°< ϕ bed < 25°, approximately, Figure 8), greater mean emulated flow thicknesses are observed at the overbanking location (road point 26, Figure 1c).Approximately in correspondence with this (Figures 6b, 6e and 6h), lower mean emulated flow thicknesses are observed at the valley bend, just on the edge of the valley (road point 27, Figure 1c; e.g., simulation 2b2, Figure 2b).These results highlight the complex interplay between flow size, mobility and the topographical features of the volcano.In general, the forcing caused by the PDC source decreases as flow size decreases and distance from the vent increases (gray/red arrow, Figure 8).This general effect is compounded by processes like PDC channelization and overbanking, which also depend on flow size and mobility, and topography surrounding the vent location.All this is compatible with previous observations made at several volcanoes worldwide (e.g., Aravena & Roche, 2022;Breard et al., 2023;Charbonnier et al., 2020;Charbonnier & Gertisser, 2012;Gueugneau et al., 2020;Ogburn & Calder, 2017;Ogburn et al., 2014;Tennant et al., 2021;Tierz et al., 2018).Notwithstanding, it is also important to place the PDC modeling developed using TITAN2D into the wider context of what is currently known in terms of PDC dynamics.For example, benchtop and large-scale analog experiments of PDCs have demonstrated that pumice flows are fluidized through porepressure feedback, which translates into a reduction in the effective friction coefficient (e.g., Lube et al., 2019;Penlou et al., 2023).This is possible due to the low permeability of the granular mixture (Breard et al., 2019).In this context, friction is actually space-and time-variant as pore pressure diffuses over time during flow propagation.Flows of larger volume are typically thicker and defluidize slower, given that the diffusion duration (in seconds) scales with h 2 /D, where h is flow thickness (in meters) and D is the hydraulic diffusion coefficient (in m 2 /s) (Roche, 2012).Such processes contribute to the inverse relationship between flow volume and flow mobility, as observed by Hayashi and Self (1992); Toyos et al. (2007); Ogburn et al. (2016) and demonstrated in Breard et al. (2018).This relationship has been previously incorporated into TITAN2D modeling strategies (Spiller et al., 2014;Rutarindwa et al., 2019, this study).However, the complexity of space-and-time-variant friction is here, and commonly elsewhere, simplified and encapsulated into single internal and bed friction parameters that are kept constant in each simulation.Finally, as introduced earlier, the TITAN2D model assumes that the pumice flows are emplaced as dense currents.This is compatible with the main body of knowledge on PDCs at Aluto volcano (e.g., Clarke, 2020;Clarke et al., 2020).Nevertheless, for dilute flows to be incorporated into the hazard assessment, one would need to use a different PDC model, which accounted for settling of particles and entrainment of air into the flow (e.g., Esposti Ongaro et al., 2007;de' Michieli Vitturi et al., 2023).This said, the results obtained in this study illustrate the great potential and high versatility of the use of statistical emulators (Bayarri et al., 2009;Rutarindwa et al., 2019;Spiller et al., 2014Spiller et al., , 2020)), particularly the zGP emulator (Spiller et al., 2023), to analyze, in a high and customizable level of detail, the particularities of PDC inundation and hazard at a given volcanic system of interest.Such analyses can expand on previous work based on selected scenarios, case studies and/or test runs (Charbonnier & Gertisser, 2012;Gueugneau et al., 2021;Ogburn & Calder, 2017;Ogburn et al., 2014).Importantly, given that zGP emulators can be built from any PDC model, the scope of these analyses can be very wide (Spiller et al., 2023).
On the distal scale, results mainly highlight the key role of vent location on PDC hazard at Aluto volcano, as it is also demonstrated by the high values of zGP correlation lengths (i.e., high influence on model outputs) of the UTMx-y TITAN2D parameters (Figure 3d).The greatest extent of the epistemic uncertainty (vertical lines in Figures 7a and  7b) corresponds with those road points (e.g., 18-24, 38-40) located closer to the areas of highest vent-opening probability (Figure 4 in Clarke et al., 2020), which tend to coincide with topographic highs (W and S-SE sectors of the caldera rim, Figure 1).The curves of probability of PDC inundation, cumulative over PDC volume (Figure 7c), confirm the key role of vent location.For example, the differences between curves are greatest when vent locations are sampled from the all-Uniform model (Figure 7c), and these differences are negative, even though mean PDC volume is slightly higher for this aleatory-uncertainty model.Given the paramount importance of vent location on PDC hazard at Aluto, future research efforts could be directed toward improving the data and models available to assess vent-opening probabilities (Clarke et al., 2020;Tierz et al., 2020).One key step forward, could be the dating of eruptive products, so specific vents could be associated with specific eruption dates, and spatial and/or temporal trends in eruption occurrence might be better constrained (e.g., Bebbington, 2013;Bebbington & Cronin, 2011;Bevilacqua et al., 2016).Challenges to find useable material to date recent eruptive products (e.g., charcoal to date pyroclastic deposits using radiocarbon) are known across the Main Ethiopian Rift (e.g., Fontijn et al., 2018).Nevertheless, any further knowledge gained about the local to regional stratigraphic/chronological framework (e.g., time ordering of eruptive vents; knowledge on the existence of "missing vents"; bracketing of eruption dates, even within large epistemic-uncertainty ranges, etc.) could enable the application of spatial and/or temporal models that improve over the existing hazard assessments (e.g., Bevilacqua, Bursik, et al., 2017;Bevilacqua et al., 2015;Selva et al., 2022).It is also important to note that emulators can play a highly significant role, not only for long-term, but also for short-term hazard assessments.Hence, they enable very rapid, yet detailed, quantifications of the probabilities of PDC inundation as estimates on key variables, such as vent location, change over time.For example, interpretation of monitoring signals during a volcanic crisis can imply that vent-opening probabilities become higher over a reduced area compared to the "pre-crisis" vent-opening domain (e.g., Sandri et al., 2012;Selva et al., 2014).Emulators can readily, and rapidly, accommodate such changes into updated, dynamic probabilistic volcanic hazard assessments (e.g., Rutarindwa et al., 2019;Spiller et al., 2020).
In terms of the comparison between the EC-MC and TITAN2D-zGP strategies, the highest differences in the probability of PDC inundation occur over the caldera floor, and on the lower SE and NE outer flanks (Figure 7d).This may relate to the EC outputs responding more markedly to the net difference in altitude between PDC-source and runout areas (e.g., Malin & Sheridan, 1982;Tierz, Sandri, Costa, Sulpizio, et al., 2016;Tierz, Sandri, Costa, Zaccarelli, et al., 2016), while the TITAN2D outputs are capturing additional topographic features and their influence on PDC propagation: namely, PDC channelization and overbanking (Figure 6; Charbonnier & Gertisser, 2012;Ogburn & Calder, 2017).Under a broader perspective, PDC model selection is still a hot topic of debate, both regarding the investigation of PDC dynamics as well as the quantification of volcanic hazard (e.g., Aravena et al., 2020;Bevilacqua, Neri, et al., 2017;Esposti Ongaro et al., 2016;Esposti Ongaro et al., 2020;Gueugneau et al., 2021;Kelfoun, 2017;Neglia et al., 2021;Ogburn & Calder, 2017;Tierz et al., 2018;Tierz, Sandri, Costa, Zaccarelli, et al., 2016).For almost a couple of decades, TITAN2D has arguably been the PDC model most widely applied to assess volcanic hazard associated with PDCs, with approaches spanning much of the spectrum between scenario-based to fully-probabilistic hazard assessments (Tierz et al., 2021): for example, Macías et al. (2008), Bayarri et al. (2009), Oramas-Dorta et al. (2012), Charbonnier and Gertisser (2012), Sandri et al. (2014), Spiller et al. (2014, 2020), Tierz et al. (2018), Rutarindwa et al. (2019), Charbonnier et al. (2020), Tennant et al. (2021), Constantinescu et al. (2022).Other PDC models with similar levels of complexity and that have also been commonly applied are VolcFlow (Kelfoun & Druitt, 2005); with a two-layer, dense-dilute version presented by Kelfoun (2017);and, recently, IMEX_SfloW2D (de' Michieli Vitturi et al., 2019).Their runtimes are expected to be shorter than TITAN2D (Gueugneau et al., 2021), thus providing, in principle, the chance to carry out forward-modeling, MC integration to estimate probabilities of PDC inundation, albeit with modest sample sizes (e.g., Aravena et al., 2020;Charbonnier et al., 2020).If the use of efficient uncertainty-quantification techniques, such as statistical emulators, is not required to compute PVHA of PDCs using these models, we envision two main implications for the future of probabilistic hazard assessment of PDCs: (a) simpler models, such as the classical EC (e.g., Tierz, Sandri, Costa, Sulpizio, et al., 2016;Tierz, Sandri, Costa, Zaccarelli, et al., 2016), the tree-branching EC (Aravena et al., 2020), or the box model (e.g., Neri et al., 2015): might lose their niche in PVHA of PDCs; and (b) statistical emulation could begin to seek application to more complex PDC models, for example, to multi-phase, multi-component models such as MFIX (e.g., Dufek & Bergantz, 2007) or PDAC (Esposti Ongaro et al., 2007;Neri et al., 2003).Nevertheless, as counter-arguments to these two implications, we also suggest that: (a) the large sample sizes required to obtain high-precision probability estimates using MC integration, especially at volcanoes where vent-opening variability is large (e.g., Bevilacqua, Bursik, et al., 2017;Rutarindwa et al., 2019;Sandri et al., 2018), could still pose issues to such an approach using VolcFlow or IMEX_SfloW2D (NB.Vent-opening variability might be smaller for stratovolcanoes, e.g., Tadini et al. (2017); Tierz et al. (2018); Tennant et al. (2021); and/or within shorter timescales, for example, during volcanic crises, e.g., Sandri et al. (2012); (b) depending on the specific conditions of each hazard analysis, the number of simulations required to successfully build statistical emulators (ranging from several hundreds up to 1,000 simulations, approximately) might still have prohibitive computational costs when using very complex multi-phase models; (c) both in relation to (a) and (b), even if tractable, the computational cost of the volcanic hazard assessment may be too high in the context of many volcano observatories, public institutions and universities worldwide.Therefore, strategies that enable the sharing of computational and other resources (including knowledge transfer), and frameworks that operate on the basis of open data, codes and workflows (e.g., Pallister et al., 2019;Tierz, 2020), will be much needed to help develop future quantitative hazard assessments of PDCs.Finally, we argue that multi-model approaches to PDC hazard assessment would also be very beneficial, as they would provide more comprehensive quantifications of the epistemic uncertainty associated with the hazard (e.g., Ogburn & Calder, 2017;Ogburn et al., 2016;Sandri et al., 2018;Selva et al., 2010), and they would also enable probabilistic model comparisons and benchmarking.In this broad context, stochastic-process emulators, and in particular the novel zero-censored Gaussian Process emulator (Spiller et al., 2023) used in this study, can play a key role, given that such statistical models can be applied to any PDC model of interest.

Conclusions
Aluto volcano, in Central Ethiopia, is known to have produced pyroclastic density currents (PDCs) during its recent-most, post-caldera eruptions (at least last 60 ka).Such volcanic flows are primarily influenced by topography and pose a significant risk to areas located on and around the volcanic edifice.Previous and recent research looking at better understanding the PDCs generated at Aluto volcano, and their associated hazard, has allowed us to develop a novel probabilistic framework to analyze PDC inundation in detail.This framework merges the TITAN2D model (Patra et al., 2005) with a novel zero-censored Gaussian Process, zGP, emulator (Spiller et al., 2023).The obtained results highlight a complex interplay between flow size (volume), mobility (bed friction angle) and topography, which also depends on the spatial scale considered, in terms of the distance (proximal, medial, distal) between the flow source and the target point analyzed.The forcing of the PDC source region can diminish the influence of the other factors and/or tune the scale of the topographic features that condition PDC propagation and inundation downstream.Low-mobility flows tend to both: (a) propagate more readily down valleys, and (b) produce enhanced flow overbanking at valley bends.Finally, we note how the high versatility of the zGP emulator is what enables great, and customizable, levels of detail to analyze PDC inundation at the target volcano.The goals of such detailed analyses, which are not achievable through classical forward modeling combined with MC integration for complex simulators, can range from investigating PDC dynamics to quantifying PDC hazard.In addition, the zGP emulator is "model-agnostic."That is, it could be applied to other PDC models, such as VolcFlow, hence contributing to: (a) improved, multi-model quantifications of epistemic uncertainty (cf.Selva et al., 2010;Tierz, Sandri, Costa, Zaccarelli, et al., 2016); and (b) detailed, probabilistic model comparisons that would complement recent model-comparison and benchmarking efforts (e.g., Esposti Ongaro et al., 2020;Gueugneau et al., 2021).

Figure 1 .
Figure 1.Aluto volcano (Central Ethiopia) and its particular topography: (a) geographical setting; (b) combined Digital Elevation Model (DEM) of Aluto volcano (2 m LiDAR DEM, Hutchison et al., 2014, over the central edifice of Aluto, 30 m ASTER DEM over the edges of the edifice onto the surrounding plain) on which three main branches of roads traversing the edifice are shown.Orange points indicate locations along the roads where zero-censored Gaussian Process (zGP) emulators are built.Topographic profiles displayed in (c) are indicated in different colors; (c) topographic profiles on a 10 m DEM, interpolated from the 30 m DEM, used to run the TITAN2D simulations in this study.The profiles highlight the general morphology of Aluto's edifice, as well as the cross-sections of some valleys analyzed in this study.

Figure 2 .
Figure 2. Simulation and emulation strategy adopted at Aluto: (a) schematic of the general strategy, which combines the TITAN2D model to simulate pyroclastic density currents (PDCs) and the zero-censored Gaussian Process (zGP) emulator to assist in the quantification of uncertainty related to PDC hazard; (b, c) selected examples of TITAN2D simulations run at Aluto, which inform about some of the topographic effects discussed in this study.Please note that inset (a) is an approximate representation of the source conditions, and flow propagation, as simulated by TITAN2D, where only the dense, basal part of PDCs is actually modeled (e.g., Tierz et al., 2018).Flow volume and bed friction angle used in the simulations shown in (b, c) are indicated on the upper part of the maps.Blue lines correspond to 0.1 m flowthickness contours.The other contours correspond to 0.2 intervals in ln(h + 1), starting in 2.4, where h is the flow thickness.Labels show flow thickness in meters.

Figure 3 .
Figure 3. Zero-censored Gaussian Process (zGP) emulator and its application to Aluto: (a) illustrative example of a zGP built to predict values from a known toy function, blue symbols correspond to positive-output and red symbols denote zero-output design points, which become negative after imputation.Colored surface represents the mean zGP predictions, which are censored to be zero on imputed, negative design points (adapted from Spiller et al., 2023); (b) plain view of the 2D toy function in (a), where the white line represents the true zero isoline and the red line denotes the zero-output isoline modeled by a classic GP emulator (colors show its mean predictions), based on the set of design points displayed in the figure ("+" symbols are positive-output and "o" symbols are zero-output design points); (c) comparison of PDC (maximum) flow thicknesses emulated using classical GP emulators (fitted to positive-output points only, top row, or to positive-and zeroresponse points, middle row) and the innovative zGP emulator (fitted to positive-and zero-response points, bottom row), at road point 35.Dots show zero-output points only, for flows of either lower (black) or higher (red) PDC volume than the one corresponding to the emulator predictions; (d) values of zGP correlation lengths (the inverse of the range parameters fitted using the Robust GaSP package, Gu et al., 2018) for each of the TITAN2D model parameters used to build the zGP emulators at each of the road points considered in this study.

Figure 4 .
Figure 4. Quantifying topographic effects on pyroclastic density current (PDC) hazard at Aluto: (a) general view of vent-opening areas used to quantify proximal (1-3 km-side black squares centered at the road point), medial (colored rectangles) and distal (dashed black box, after Clarke et al., 2020) topographic effects; (b) details on the procedure to quantify proximal effects, Δz denotes the difference between the mean altitude across the vent-opening area and the altitude at the road point, and γ i corresponds to the variable in Equation 5; (c, d) zoomed-in maps to zones 1 and 2, respectively, in the context of analyzing medial topographic effects on PDC hazard.
One such model (full-aleatory-uncertainty model) incorporates all information available for Aluto for each of the five uncertain TITAN2D parameters considered: vent radius and flux rate pdfs calculated afterClarke et al. (2020), using the joint-parameterization strategy to obtain appropriate values for the TITAN2D parameters; bed friction angles derived from the aforementioned log 10 -log 10 linear regression; and a bi-variate, UTMx-y pdf taken from the kernel density estimate model presented inClarke et al. (2020), based on previous work by C. B.Connor and Hill (1995);Weller et al. (2006); L. J.Connor et al. (2012).One alternative model (partial-aleatory-uncertainty model) is similar to the full-aleatory-uncertainty

Figure 5 .
Figure 5. Summary results from the proximal topographic analyses: plots of total angular deviations (γ i , normalized) of the shape of the vent-opening area from a funnel centered at each road point, versus the difference (in meters) between the average altitude across the vent-opening area and the altitude at each road point (Δz), considering (a) high-volume, low-bed-friction pyroclastic density currents (PDCs) and (b) low-volume, high-bed-friction PDCs sourced within a 1 km-side ventopening area centered at each road point; and (c) high-volume, low-bed-friction PDCs and (d) PDCs of any volume and bed friction sourced within a 3 km-side ventopening area centered at each road point.

Figure 6 .
Figure 6.Summary results from the medial topographic analyses.Bed friction angles versus the mean emulated flow thicknesses when the zero-censored Gaussian Process (zGP) emulators are evaluated at those values of bed friction and vent locations are uniformly sampled from within vent-opening areas defined as medial zones 1 and 2 in Figure4.All road points are sorted, within each column, from those easternmost, E, to those westernmost, W (please see Figures1 and 4for topographical and morphological context for these road points).Rows show zGP evaluations at different pyroclastic density current volumes: ∼0.1 km 3 (a-c), ∼0.3 km 3 (d-f), and ∼0.6 km 3 (g-i).

Figure 7 .
Figure 7. Summary results from the distal topographic effects: (a, b) probabilities of pyroclastic density current (PDC) inundation, P(PDC | distal) considering two aleatory-uncertainty models (full-aleatory-uncertainty and all-Uniform models, see text for more details), and across two spatial profiles: (a) E-W profile over the caldera floor of Aluto, and (b) N-S profile traversing the NE and SE flanks, as well as the eastern sector of the caldera floor.Blue colors indicate the road points where P (PDC | distal), given the full-aleatory-uncertainty model, is higher than P(PDC | distal) given the all-Uniform model.Red colors indicate the opposite case.Black dots denote the probabilities of PDC inundation based on the EC-MC analysis of Clarke et al. (2020) (afterTierz, Sandri, Costa, Sulpizio, et al., 2016; Tierz, Sandri, Costa,  Zaccarelli, et al., 2016).(c) PDC volume plotted against the conditional probability of PDC inundation given zGP evaluations at PDC volumes equal or smaller than each point along the x-axis, for a zGP emulator built at the location of the Aluto-Langano geothermal power station (road point 2).The area between the different curves, generated from different aleatory-uncertainty (AU) models (see text for more details), is a measure of the extent of the epistemic uncertainty (cf.Rutarindwa et al., 2019;  Tierz, Sandri, Costa, Sulpizio, et al., 2016); (d) map showing the differences in probability of PDC inundation at each road point between the values reported inClarke et al. (2020) [grid points closest to each road point] using an EC-MC strategy, and those of the partial-aleatory-uncertainty model, computed using the TITAN2D-zGP strategy presented in this study (see text for more details).

Figure 8 .
Figure 8. Summary of the influence of topographic effects on pyroclastic density current (PDC) inundation at Aluto volcano.A complex interplay between flow size (volume), flow mobility (bed friction angle) and the scale and location of topographic features is observed.For example, *α: a high-volume PDC in proximal areas responds only to large-scale topographic features; bed friction has a minimal influence in this situation.*γ: a low-volume PDC at medial-to-distal distances responds to both large-and small-scale topographic features, as well as to the bed frictional properties of the flow.See text for more details.