Probabilistic Tsunami Hazard Analysis: High Performance Computing for Massive Scale Inundation Simulations

Probabilistic Tsunami Hazard Analysis (PTHA) quantifies the probability of exceeding a specified inundation intensity at a given location within a given time interval. PTHA provides scientific guidance for tsunami risk analysis and risk management, including coastal planning and early warning. Explicit computation of site-specific PTHA, with an adequate discretization of source scenarios combined with high-resolution numerical inundation modelling, has been out of reach with existing models and computing capabilities, with tens to hundreds of thousands of moderately intensive numerical simulations being required for exhaustive uncertainty quantification. In recent years, more efficient GPU-based High-Performance Computing (HPC) facilities, together with efficient GPU-optimized shallow water type models for simulating tsunami inundation, have now made local long-term hazard assessment feasible. A workflow has been developed with three main stages: 1) Site-specific source selection and discretization, 2) Efficient numerical inundation simulation for each scenario using the GPU-based Tsunami-HySEA numerical tsunami propagation and inundation model using a system of nested topo-bathymetric grids, and 3) Hazard aggregation. We apply this site-specific PTHA workflow here to Catania, Sicily, for tsunamigenic earthquake sources in the Mediterranean. We illustrate the workflows of the PTHA as implemented for High-Performance Computing applications, including preliminary simulations carried out on intermediate scale GPU clusters. We show how the local hazard analysis conducted here produces a more fine-grained assessment than is possible with a regional assessment. However, the new local PTHA indicates somewhat lower probabilities of exceedance for higher maximum inundation heights than the available regional PTHA. The local hazard analysis takes into account small-scale tsunami inundation features and non-linearity which the regional-scale assessment does not incorporate. However, the deterministic inundation simulations neglect some uncertainties stemming from the simplified source treatment and tsunami modelling that are embedded in the regional stochastic approach to inundation height estimation. Further research is needed to quantify the uncertainty associated with numerical inundation modelling and to properly propagate it onto the hazard results, to fully exploit the potential of site-specific hazard assessment based on massive simulations.

Probabilistic Tsunami Hazard Analysis (PTHA) quantifies the probability of exceeding a specified inundation intensity at a given location within a given time interval. PTHA provides scientific guidance for tsunami risk analysis and risk management, including coastal planning and early warning. Explicit computation of site-specific PTHA, with an adequate discretization of source scenarios combined with high-resolution numerical inundation modelling, has been out of reach with existing models and computing capabilities, with tens to hundreds of thousands of moderately intensive numerical simulations being required for exhaustive uncertainty quantification. In recent years, more efficient GPU-based High-Performance Computing (HPC) facilities, together with efficient GPU-optimized shallow water type models for simulating tsunami inundation, have now made local long-term hazard assessment feasible. A workflow has been developed with three main stages: 1) Site-specific source selection and discretization, 2) Efficient numerical inundation simulation for each scenario using the GPU-based Tsunami-HySEA numerical tsunami propagation and inundation model using a system of nested topo-bathymetric grids, and 3) Hazard aggregation. We apply this site-specific PTHA workflow here to Catania, Sicily, for tsunamigenic earthquake sources in the Mediterranean. We illustrate the workflows of the PTHA as implemented for High-Performance Computing applications, including preliminary simulations carried out on intermediate scale GPU clusters. We show how the local hazard analysis conducted here produces a more fine-grained assessment than is possible with a regional assessment. However, the new local PTHA indicates somewhat lower probabilities of exceedance for higher maximum inundation heights than the available regional PTHA. The local hazard analysis takes into account small-scale tsunami inundation features and nonlinearity which the regional-scale assessment does not incorporate. However, the deterministic inundation simulations neglect some uncertainties stemming from the simplified source

INTRODUCTION
Tsunamis are infrequent hazards that can potentially lead to devastating consequences. Earthquakes are the most common source of tsunamis (about 80% of tsunamis worldwide, see e.g. the NCEI global tsunami database: https://www.ngdc.noaa. gov/hazard/tsu_db.shtml) and we restrict the analysis here to seismic sources (coined seismic PTHA or S-PTHA, see e.g., Lorito et al., 2015). Their induced tsunamis pose a risk toward the global coastal population, both related to human casualties (Løvholt et al., 2012), direct economic losses (e.g., Løvholt et al., 2015), to critical infrastructures useful for crisis management (e.g., harbors: Argyroudis et al., 2020), or through secondary cascading events such as in the Fukushima event (e.g., Synolakis and Kânoglu, 2015). The uncertainty in tsunami hazard models is great, resulting from the infrequency of events (and consequent relatively small datasets of past events), from the vast number of potential sources and tsunami-generating mechanisms (e.g., Grezio et al., 2017;Davies et al., 2018), from the accuracy of high resolution topo-bathymetry and friction models for inundation calculations (e.g., Park et al., 2014;Bricker et al., 2015;Griffin et al., 2015;Sepúlveda et al., 2020), and from the approximations of numerical simulations (e.g., Behrens and Dias, 2015). Among these, the hazard is largely controlled by the source probability of occurrence which is highly uncertain, in particular for large events like megathrust subduction earthquakes of the scale of the 2004 Indian Ocean and 2011 Tohoku event (Lay et al., 2005;Kagan and Jackson, 2013) or large crustal events (Basili et al., 2013). Furthermore, the hazard and related uncertainty stem from the complexity and variety of the various types of earthquake mechanisms such as tsunami earthquakes (e.g., Newman et al., 2011), outer rise events (e.g., the 2009 Samoa tsunami: Fritz et al., 2011), other significant unknown or only partially known crustal sources (Basili et al., 2013;Selva et al., 2016), variable slip (e.g., Geist, 2002;Scala et al., 2020), or generally due to unexpected source mechanisms such as revealed for the Palu tsunami in 2018 (Ulrich et al., 2019). Moreover, tsunamis often happen simultaneously with other hazards, and may interact with them (e.g. earthquakes, landslides, or volcanoes) in a complex manner (e.g., Goda and De Risi, 2018;Pitilakis et al., 2019;Argyroudis et al., 2020). This was demonstrated by the 2018 Palu earthquake and tsunami, where the earthquake (Bao et al., 2019), liquefaction (Cummins, 2019) and tsunami (e.g., Omira et al., 2019;Ulrich et al., 2019) impacted almost simultaneously. Clearly, it is important to have well established methods that can capture these complexities to represent the hazard.
In recent years, Probabilistic Tsunami Hazard Analysis (PTHA: Geist and Parsons, 2006;Grezio et al., 2017) has become the standard way of estimating this complex tsunami hazard. PTHA estimates the probability of exceeding a specified tsunami metric (e.g. flow depth, tsunami height, or momentum flux) at a given location within a given time interval, such as the probability of exceeding a specified inundation height within the next 50 years. Tsunami observations are usually not sufficient for constraining a PTHA. Computation-based PTHA considers a discretization of the total hazard into many potential source scenarios, together with the estimated probability of occurrence of each scenario, for as many scenarios as necessary to represent the expected natural variability in a probabilistic source model. To resolve tsunami source uncertainty adequately, many thousands or sometimes even millions of scenarios need to be simulated (e.g., Selva et al., 2016). To be feasible, PTHA applications have therefore often resorted to estimating the hazard offshore and extrapolating it onshore, and sometimes applying stochastic inundation modeling (e.g. Power et al., 2007;Burbidge et al., 2008;Horspool et al., 2014;Davies et al., 2018;Glimsdal et al., 2019). It has so far not been possible to conduct tsunami hazard analysis running this broad range of scenarios without renouncing some details needed for practical applications. For instance, the first widely accepted probabilistic tsunami hazard map for Europe was developed within the TSUMAPS-NEAM project (http://www.tsumaps-neam.eu/), which covers the hazard in the North-eastern Atlantic, the Mediterranean and connected seas (NEAM). This so-called NEAMTHM18 assessment (Basili et al., 2018;Basili et al., 2019) includes millions of sources, but estimates inundation probability at regional scales based on offshore analysis of tsunami simulations (Glimsdal et al., 2019), and so lacks the highresolution inundation simulation typical of sitespecific PTHA.
Local scale applications to date have needed to reduce the number of simulations dramatically in order to be feasible. González et al. (2009) provided the first local PTHA, generating inundation maps for a location in Oregon using the MOST simulation software (Titov and Gonzalez, 1997) on a system of nested grids. However, this study was limited to a small number of megathrust earthquake scenarios deemed to dominate the hazard at the target region. Lorito et al. (2015) address the feasibility of inundation maps for Seismic PTHA (S-PTHA) for the Mediterranean region, with strategies for reducing the number of simulations required to optimize accuracy in the hazard. This work was taken further by Volpe et al. (2019) who emphasized the need to differentiate between near-field and far-field sources due to the alteration of the coastal height as a consequence of co-seismic displacement. However, to date, there exists no available benchmark that covers the source variability sufficiently, which can be used to test the degree to which the above simplifications are viable. A necessary element that has been lacking moving forward on this aspect, is the availability of computation resources and related workflows that allow effective use of PTHA on the major computational facilities, namely Tier-0 High Performance Computing (HPC) systems .
We here attempt to bridge the gap between regional-scale PTHA and scenario-specific local inundation simulation, by developing and prototyping a new workflow for site-specific high-resolution PTHA using HPC. This local PTHA workflow is a so-called pilot demonstrator of the H2020 funded ChEESE Center of Excellence (https://cheese-coe.eu/) for addressing geophysical problems related to solid earth processes using (future) pre-Exascale and Exascale supercomputers. A comprehensive PTHA with high resolution inundation calculations at local scale is a problem that is only tractable given such computational resources. In this paper, we used as a starting point for our site-specific analysis the existing regional-scale NEAMTHM18 tsunami hazard assessment (Basili et al., 2019). Hence, we will not assess source probabilities from scratch here, but rather use those probabilities derived by Basili et al. (2019) as input to our analysis. Our primary objective is to extend this regional analysis to local hazard combining high resolution topobathymetric data with nonlinear shallow water (NLSW) inundation modeling using the multi-GPU finite volumes Tsunami-HySEA model (de la Asunción et al., 2013;Macías et al., 2017;Macías et al., 2020a;Macías et al., 2020b), restricting the simulations to the sources deemed relevant by NEAMTHM18 at the specific site considered. To this end, we present a new workflow embedding Tsunami-HySEA into PTHA and demonstrate the suitability for HPC usage. We also present comparison with the previous NEAMTHM18 analysis for the offshore tsunami hazard, as well as new inundation hazard curves and maps. This paper is organized as follows: In The Seismic Probabilistic Tsunami Hazard Assessment in the Mediterranean Region: The Regional Model NEAMTHM18, we describe briefly the NEAMTHM18 analysis relevant for creating the input to the local hazard. In Implementation of a High-Performance Computing Oriented Seismic Probabilistic Tsunami Hazard Analysis Workflow, the PTHA workflow is described, including the source disaggregation from NEAMTHM18 and the inundation simulations with Tsunami-HySEA on Tier-0 GPU clusters. The setup for the example case presented here, Catania harbor, is described in Setup for Hazard Analysis Toward the City of Catania. In Results, the inundation calculations and hazard aggregation are discussed in the context of previous results (Setup for Hazard Analysis Toward the City of Catania). We finally provide future perspectives.

THE SEISMIC PROBABILISTIC TSUNAMI HAZARD ASSESSMENT IN THE MEDITERRANEAN REGION: THE REGIONAL MODEL NEAMTHM18
The NEAMTHM18 tsunami hazard model provided a rigorous analysis of the annual rates of possible earthquake events and of the tsunami hazard curves for the coastlines of the NEAM region using many millions of scenarios. It is also the first community and consensus based regional tsunami hazard assessment in the NEAM, where the quantification of epistemic uncertainties was heavily based on expert opinion distilled through formal elicitation processes. For details related to the establishment of the source probabilities, we refer to Basili et al. (2019).
The NEAMTHM18 PTHA considered two types of earthquake sources (Selva et al., 2016), coined Predominant Seismicity (PS) and Background Seismicity (BS). PS consists of earthquakes associated with subduction interfaces and major fault systems, where the fault geometries and mechanisms are relatively well understood. The second (BS) class comprises (crustal) seismicity associated with other fault systems, the knowledge and geometry of which may be more incomplete. The BS earthquakes can, in principle, occur anywhere. This distinction was adopted for the sources established in TSUMAPS-NEAM and Basili et al. (2019) provide a comprehensive description of the employed source discretization. In the Mediterranean, PS comprises three main subduction interfaces: the Calabrian, Hellenic and Cyprian arcs. Here, PS scenarios are defined by a specified slip on each element of triangular meshes, modelling 3D fault geometries. BS scenarios are defined over a regular grid covering the entire Mediterranean Sea and nearby lands, with sea-bottom deformations modelled by considering uniform slip on rectangular faults (Okada, 1992). The set of tsunamigenic scenarios is defined by systematic discretizations of the earthquake parameters describing these two classes of seismicity.
The NEAMTHM18 tsunami scenarios were produced by linear combinations of previously simulated elementary Gaussian sources . These simulations involved approximately 200,000 Tsunami-HySEA simulations carried out offshore for the entire NEAM region. This was possible since the offshore tsunami simulations were linear, and linear combinations could be employed to provide a much higher number of scenarios. Amplification factors translated the offshore wave characteristics to inundation height depending on the local bathymetry, and on the polarity and dominant period of the incident wave (Glimsdal et al., 2019). Amplification factors simplify the assessment for a regional hazard quantification. However, they do not capture the nonlinearity, nor include the detailed dynamics of inundation on local topography. This quantification is associated with very large uncertainties, mainly epistemic. These uncertainties are quantified by hundreds of NLSW inundation simulations carried out at different sites, and stem largely from the topographic variability onshore. These uncertainties are then treated by means of conditional probabilities as a function of the amplified inundation heights (Glimsdal et al., 2019). However, these conditional probabilities are not site-specific, as they were estimated by aggregating inundation simulations from a variety of different coastal sites. The application of local high-resolution topography and inundation simulations here is mainly to try to reduce this epistemic uncertainty.
The number of scenarios considered in NEAMTHM18, restricted to the Mediterranean area only, is in the order of 10 6 (Basili et al., 2019). Even with massive Tier-0 HPC resources available, the number of source scenarios in NEAMTHM18 must be dramatically trimmed down to be feasible for local tsunami hazard computation. For this purpose, the most important sources can be determined through hazard disaggregation (e.g., Bazzurro et al., 1999;for tsunamis: Selva et al., 2016;Power et al., 2018). In the next section, we explain the entire workflow for the local PTHA, including details of the disaggregation in the context of NEAMTHM18.

Overview of the local Probabilistic Tsunami Hazard Analysis Workflow
The local PTHA workflow consists of the following main components, shown in Figure 1, all elaborated in separate subsections below: • Provision of user specifications, including: definition of hazard metrics, thresholds for the specified metrics, availability of computational resources, and physical input parameters for the hydrodynamic simulations (e.g., topo-bathymetric grids, Manning friction coefficient, CFL number, dry land threshold value, etc.). • Source selection of scenarios, here using the NEAMTHM18 disaggregation as input (Disaggregation and Source Selection). In this step, it is also possible to refine the sources (to provide more sources to cover more broadly the source variability, for example by more finely sampling location and slip distribution of local sources). This source refinement is not used in the example studies provided herein. • A micro-workflow for HPC inundation simulation (High-Performance Computing Inundation Simulations and Workflow), using the NLSW model Tsunami-HySEA as the computational engine for the inundation simulations, capable of managing large ensembles of hundreds of thousands multi-GPU simulations; • A hazard aggregation step, combining the different model runs to provide hazard curves and maps for potentially inundated areas (Hazard Aggregation). This step also manages the epistemic uncertainty by considering the alternative modelling and/or parametrizations for the seismic sources, their recurrence rates, and amplification models, providing the hazard's uncertainty statistics (returning a mean and percentiles). This quantification of epistemic uncertainty is not implemented here, and only the mean hazard curves and inundation maps are presented for the sake of conciseness.

Disaggregation and Source Selection
Local tsunami hazard analysis utilizes non-linear models. Hence it cannot exploit superposition of unit sources as in conventional regional-scale PTHA (Burbidge et al., 2008;Basili et al., 2019). As explained above, it may not be feasible to simulate millions of scenario simulations for local inundation analysis. However, for a specific site, only a limited set of these sources contributes significantly to the hazard. To identify the most significant scenarios, a disaggregation analysis is carried out on the regional hazard estimated in NEAMTHM18 as the first step of the local PTHA workflow. The disaggregation algorithm first ranks all the scenarios contributing to the target site (i.e. to one or more offshore points close to the target site) for a given intensity, or intensity interval, according to their relative importance for the site of interest, measured as their relative contribution to the local offshore hazard curve in terms of mean annual rates. Then, a desired "degree of accuracy" to which the local hazard should be approximated can be defined in terms of the resemblance of the original offshore hazard curves and of the ones calculated only with the selected scenarios, corresponding to the given intensity value or interval. This degree of accuracy formally corresponds to the probability that the occurrence of the target event (a tsunami in the selected interval) is caused by one of the selected sources, as computed from standard disaggregation (e.g., Bazzurro and Cornell, 1999). An example of an application of this procedure using the 1-4 m interval is given in Scenario Selection and Representation of Probabilities. The degree of accuracy controls the number of simulations required, and it should be selected based on the available resources and the computational cost, which mainly depends on the size and the spatial resolution of the computational grids. The input parameters to the disaggregation procedure applied here consist of 1) the selection of the offshore Points of Interest (PoIs) to consider (should be in the vicinity of one or more sites for modelling inundation), 2) the mean annual rate of the set of scenarios retrieved from regional hazard, and 3) a specified value or range of values of interest, provided in terms of Maximum Inundation Height (MIH). We here define the MIH as the largest inundation height above mean sea level at any onshore point in the computational domain. However, when MIH is derived from the regional assessment using amplification factors, it represents a larger area, and must hence be interpreted as a stochastic quantity inheriting the spatial variability of the local MIH.
From NEAMTHM18, we obtain a first order source discretization and corresponding probabilities. The disaggregation of the TSUMAPS-NEAM assessment provides, for a predefined PoI, a list of those earthquake scenarios which dominate the tsunami hazard locally as estimated from offshore simulation results. When the available computational resources allow, we can subsequently perform a refinement of sources and corresponding splitting of probability, to improve the existing Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 source discretization for the local hazard. We stress that this source splitting step is not carried out in the example cases presented in this paper, but it nevertheless is an essential element of the local PTHA workflow. The source refinement procedure depends on the source type and, in this case, should be closely tied to the PS and BS source definition in NEAMTHM18. The PS sources may take larger Moment Magnitudes (M w ) that are represented by stochastic fields of heterogeneous slip, embedding the possibility of shallow slip amplification controlled by rigidity and coupling variations with depth (Scala et al., 2020). For the largest magnitudes, the PS sources were modeled by five different stochastic realizations for each given source region. In the local PTHA workflow, the source refinement may allow the user to extend this to an arbitrary increase of heterogeneous slip realizations (to e.g., 20, 30, 40 etc.) depending on computational resources available. Similarly, each individual BS source can be refined with greater resolution in location, fault orientation, and focal mechanism variability. For both the PS and the BS sources, source probabilities from NEAMTHM18 are split from the "parent" scenarios to the "children" scenarios which constitute the refinement of the parent sources, redistributing the total mean annual rates of parent scenarios. This procedure, with or without refinements, ends with the definition of a list of scenarios to be modelled through HPC. To ensure the feasibility of PTHA, an iterative step in addition to the disaggregation and source refinement procedures is required. This step considers the computational resources available for the tsunami simulations and assesses how many scenarios can be simulated with these constraints. If necessary, a lower disaggregation level is selected, and the scenario list is consequently reduced. Similarly, a reduction of refinement level can be considered.

High-Performance Computing Inundation Simulations and Workflow
Tsunami-HySEA (de la Asunción et al., 2013) is an NLSW model implemented in CUDA for NVIDIA GPU computations and parallelized with MPI for running in multi-GPU architectures. Tsunami-HySEA models both open-ocean tsunami propagation and nested grid inundation using progressively finer grid resolution of the coastal areas in a single code. The code has undergone an intensive process of model validation and verification following, in particular, the benchmarking standards of the NTHMP, the National Tsunami Hazard and Mitigation Program, USA (Macías et al., 2017;Macías et al., 2020a;Macías et al., 2020b). The nested grid meshes are fixed with an arbitrary number of levels satisfying power-of-two refinement ratios. The nested grid algorithm updates the nested grids at coarser grid levels through spatial projection of the mesh values at the finer grids levels. In addition, values at nodes of the coarser grids along the boundaries of each finer grid level are used to drive the simulation on the finer grid. In this way, a two-way update between the fields on each grid level is performed. A new nested grid algorithm has been implemented for the PTHA simulations and adopted thereafter. Specifically, in the new algorithm, the values of the FIGURE 1 | Overview of the Probabilistic Tsunami Hazard Analysis (PTHA) workflow. The step involving HPC resources is shown in the green box. Preprocessing and hazard aggregation is performed outside the HPC system. As large amounts of data can be produced, attention is paid to continuous data transfer during the entire HPC project.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 free surface and water height are projected on the coarse grids and the bathymetry on the coarse grids is modified accordingly to preserve the consistency in the relation h (x,t) H (x,t)+s(x,t), where H is the bathymetry, h the total water depth, and s the sea surface elevation (from a fixed reference level based on the bathymetric depth value). The Tsunami-HySEA code is formulated using Finite Volumes and is implemented for multi-GPUs utilizing 2D domain decomposition with load balancing. The efficient use of GPUs makes it one of the most efficient NLSW models available . It includes methods for conveying seabed deformations into tsunami generation based on full potential hydrodynamic theory (Kajiura, 1963;Nosov and Kolesov, 2007). For treating the overland flow, a quadratic Manning friction term is implemented. Moreover, a minimum dry land threshold depth h m is adopted for stability purposes.
The Tsunami-HySEA code has recently been optimized for improved HPC efficiency. The output is written into NetCDF files that are generated synchronously, or asynchronously using C++11 threads for efficient input-output. The output, one file for each of the two finest sub-meshes and one for storing time series, consists of data that has been compressed using the algorithm described in Tolkova (2007) to further reduce the size of the output files. Tsunami-HySEA has been tested on four Tier-0 supercomputers: the CTE-POWER cluster in the Barcelona Supercomputing Center (BSC), where it has already been used intensively, the DAVIDE cluster and the new MARCONI100 machine at the CINECA consortium in Italy, and Piz-Daint GPU at the Swiss National Supercomputer Center. Relevant to the present application is the synchronous simulations of many independent Tsunami-HySEA simulations in parallel, naturally achieving perfect parallelism (sometimes coined "embarrassing parallelism"). To this end, Tsunami-HySEA provides good so-called weak scalability on the different architectures it has been tested on, meaning that the computational speed is only marginally reduced when many simulations are run synchronously using embarrassing parallelism. For example, the weak scalability obtained for the present set of scenarios was 99.94% when using four GPUs and 98.96% for 64 GPUs on the Marconi-100 supercomputer. Hence, Tsunami-HySEA is well placed to utilize Tier-0 resources to conduct the high number of scenario simulations necessary for PTHA.
Tsunami-HySEA is the engine of the inner PTHA workflow that is carried out within the HPC Tier-0 environment, shown as the green box in Figure 1. This inner workflow consists of Faster Than Real Time (FTRT) numerical simulations for each of the selected earthquake scenarios from the disaggregation step on the provided system of nested grids. The nested grids allow for high resolution inundation calculations at the coastal region of interest while keeping a coarser spatial resolution for the open sea wavepropagation, appropriate for the temporal and spatial scale of the wave. It is the system of nested grids, and in particular the grids with the finest resolution, which dominates the simulation time. A more detailed look at the inner-workflow for the HPC resources is provided in Figure 2. This provides a closer look at the nested grid structure and illustrates the calculation of inundation maps exemplified for Catania for four example scenarios derived from the disaggregation and source splitting step explained above: two of the BS class, crustal seismicity, and two of the PS class, subduction earthquakes. Moreover, the innerworkflow also includes procedures for filtering spurious simulation results. In the cases where such simulations are discovered, scenarios can be removed from the assessment, and all other probabilities can be normalized. In the analysis presented here, such spurious simulations were not detected, and the renormalization option was not invoked.

Hazard Aggregation
The results from all the simulations are combined in the hazard aggregation step. This consists also of the postprocessing and visualization that are performed outside of the HPC resources since these tasks do not require the same level of performance as the simulation calculations themselves. In the aggregation step, the different inundation simulation results are combined to provide hazard curves (probability of exceeding a given MIH during a given exposure time). For a local inundation site, a large number (typically hundreds of thousands) of such curves will be available, with one curve for each inundated point.
These curves are obtained in the following way (for a more complete description we refer to Basili et al., 2019); the individual scenario list and related mean annual rates are derived from the disaggregation and source refinement step, representing the most relevant sources for the local site of interest. For each of these scenarios, we can have different estimates of the mean annual rates and, consequently, of probabilities in the reference exposure time, depending on the epistemic uncertainty. In the aggregation, we assume a time-independent Poisson process implying that the exceedance probabilities can be computed using: where P(I > I C ,ΔT) represents the probability of the tsunami metric I exceeding a threshold value I C during an exposure time period ΔT. λ(I > I C ) represents the mean annual rates for an individual source giving rise to the tsunami exceeding the given impact metric. Due to their independence, all these rates are summed for deriving the probability of at least one exceedance (PoE) in ΔT for a given threshold I C .
In this study, for simplicity, we assume that λ(I>I C ) can be evaluated, for each individual source, as the product of the mean annual rate of the source and of an identity function which equals 1 if the simulated tsunami exceeds I C , and 0 otherwise (e.g., Grezio et al., 2017). This procedure neglects the potential uncertainty stemming from the tsunami generation, propagation and inundation model, as well as the one due to the oversimplifications of the source model (Choi et al., 2002). Since this study has the main purpose of illustrating the new workflow, rather than an assessment for operational purposes, we considered our approach sufficiently elaborated. For a specific application, uncertainties can be applied in the form of a conditional probability (Glimsdal et al., 2019) but starting from inundation simulations rather than from offshore results. Ideally, the parameters of the distributions (usually assumed lognormal, Glimsdal et al., 2019 and reference therein) should be calibrated with run-up observations. Because we based all the present results on a subset of the potential sources selected through disaggregation, the mean annual rates are finally normalized by the degree of accuracy chosen in the disaggregation step. Assuming that the subset can be considered an unbiased sample of all the sources that may generate local tsunamis, we may compensate the removal of potentially impacting sources by re-normalizing the mean annual rates of the selected sources. For example, for a 99% disaggregation level, λ(I>I C ) can be simply multiplied with a normalization factor 1/0.99. The PoE can also be converted into average return periods T through the expression T ΔT/ abs(ln(1 -P(I > I C ,ΔT)). In the application presented herein, we compute PoEs assuming a ΔT 50-years exposure time. For example, a 2% PoE in 50 years then yields T ≈ 2,475 years, and a 10% PoE in 50 years yields T ≈ 475 years.
In this paper, we consider the following impact metrics, the MIH, the maximum flow depth H m , and the maximum depth averaged momentum flux defined as the maximum of the instantaneous product of the square velocity and flow depth (U 2 H) m .
When multiple models for source probabilities and/or for tsunami generation, propagation and inundation are implemented, the results are represented as a family of different hazard curves, where each curve represents one realization of the epistemic uncertainty (Grezio et al., 2017). The epistemic uncertainty in annual rates may be inherited from the regional hazard model (here, NEAMTHM18). The epistemic uncertainty in tsunami generation, propagation, and inundation should be evaluated by analyzing the impact of source simplifications and limitations of the numerical tsunami simulations. However, in this paper, only the average of the epistemic uncertainty from the source model is presented: the weighted average of the alternative mean annual rates from NEAMTHM18. The alternatives are weighted according to the relative credibility assigned to the different considered source models in NEAMTHM18. The coarser intermediate grids (40 m and 160 m) were derived from the 10 m grid by using a bilinear resampling algorithm to assure depth compatibility among all the nesting involved in the simulations. To ensure stability, Tsunami-HySEA simulations use a threshold for water height at inundation (here set to 0.001 m), and a maximum allowed current velocity as a simulation stopping criteria (here set to 40 m/s). A uniform Manning coefficient n 0.03 is employed for representation of the bed friction. A CFL number of 0.5 is used for the time stepping.

Scenario Selection and Representation of Probabilities
We carried out the disaggregation analysis for PoIs close to the East coast of Sicily, Italy, offshore Catania, based on the NEAMTHM18 regional hazard assessment. To make use of an example of potential practical interest, in this specific application, we selected the interval of MIH derived from the regional assessment between 1 m and 4 m. This should allow the inclusion, up to the 98th percentile, all the scenarios which may generate a run-up comparable to the design value used for the construction of the tsunami evacuation maps in the Catania plain (Dipartimento della Protezione Civile, 2018; Basili et al., 2019). Locally, in fact, the design MIH value used for the subsequent preparation of the evacuation maps is 1.2 m. Since the disaggregation procedure (Disaggregation and Source Selection) attempts to approximate the total hazard by ranking the most probable scenarios, and since these highest ranked scenarios generate a smaller MIH in the considered interval, the scenarios in proximity of 1.2 m are almost surely included in the selection when using the 1-4 m interval. Note that the corresponding run-up design value is chosen as three times the design MIH. Based on the reanalysis by Basili et al. (2019) of the simulations performed by Glimsdal et al. (2019), the 98% of the run-up values in the vicinity of the POI where the MIH is estimated would not exceed three times the design value. Figure 3 (panel A) displays a trade-off curve of the number of scenarios required to quantify the hazard at the Catania PoI as a function of the level of approximation of the total hazard in the 1-4 m interval. In this study, we chose to disaggregate to retain those sources that collectively have a 99% probability of causing hazards in the target interval. Adding up their individual contributions allows the reproduction of 99% of the total hazard. We indicate the number of scenarios required to reproduce 90, 95, and 99% of the hazard. We performed a number of disaggregations over different and broader intervals and note that very few additional scenarios result from extending the upper limit of the interval, for example, from 4 m to 8 m. The scenarios that result in the highest inundations are associated with exceptionally low probabilities and few of these scenarios contribute to the PTHA for the time-intervals of interest. This is confirmed in panel B) where the hazard curves calculated from the full NEAMTHM18 model and approximated at the 99% level are nearly the same, providing confidence that neglected scenarios do not contribute significantly in the considered MIH interval. Altogether 32,363 scenarios were obtained from the disaggregation. Of these scenarios, 11,120 were of type BS and 21,243 were of type PS. We display these scenarios in Figure 4, with BS and PS sources covered in panels (A) and (B) respectively. The colors indicate the cumulative rates of the sources, with the darker symbols and shapes indicating a higher probability of seismic slip at the location indicated within the relevant time interval. Both panels display simplifications of the earthquake generated tsunami sources. It can be qualitatively assessed that the hazard is driven by subduction earthquakes and by local crustal sources, consistently with the findings of Selva et al. (2016). For each symbol in panel (A), many earthquake scenarios defined by Okada source parameters have their hypocenter within the region indicated. The scenarios vary in all the other parameters: depth, rupture width, rupture length, slip, and the angles of strike, dip, and rake. The lateral dimensions of the rupture for any one earthquake scenario may exceed significantly the size of the pixel displayed. Similarly, for the PS sources in panel (B), each scenario consists of a slip distribution represented across many elements of a triangular mesh. Since variation in the slip distribution can have significant consequences for inundation, NEAMTHM18 therefore considered several stochastic realizations of slip for a given fault geometry (see Scala et al., 2020). The color of each triangular cell in panel (b) indicates the number of PS scenarios which have a non-zero slip on that particular cell.
For epistemic uncertainty quantification, we adopt the ensemble of mean annual rates describing the epistemic uncertainties on seismic rates in NEAMTHM18. This implies that each scenario simulation is associated with 1,000 alternative estimates of the mean annual rates λ to represent the epistemic uncertainty, that are sampled from the epistemic uncertainty alternative tree (Selva et al., 2016). These epistemic uncertainties represent uncertainty related to the seismicity models such as the Magnitude Frequency Distributions (MFDs), scaling relation used for the earthquake scenario, crustal rigidity model etc., as well as to the simplified inundation model (Glimsdal et al., 2019). All these uncertainties have been unraveled in NEAMTHM18 through a thorough expert elicitation process. As mentioned above, we here present only results based only on averaging the epistemic uncertainty.

Comparison of Offshore Hazard Curves
Before assessing onshore hazard curves, we first verify that the offshore tsunami hazard obtained with the new simulations is consistent with the regional scale hazard estimated in the NEAMTHM18 assessment. This is a general consistency test. Even if all parts of the study are implemented correctly, differences may arise because of several specific reasons. Firstly, we use here only the subset of sources selected after disaggregation and we renormalize the source mean annual rates accordingly to better approximate the hazard curves, even if we FIGURE 4 | Cumulative annual rates in NEAMTHM18 source discretization resulting in a maximum inundation height (MIH) in the range 1-4 m for Catania (white square) based upon the offshore tsunami simulation result, displayed together with topo-bathymetry of the Mediterranean. The outlines of the PS-source meshes are displayed in Panel (A), demonstrating how crustal (BS) sources are also defined in the same locations as PS sources. "Remote" BS sources, geographically separated from the contiguous source clusters, may appear in the disaggregation given eg marginally more efficient tsunami wave propagation. However, their contribution to the total hazard at Catania is likely to be small given the total number of scenarios considered.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 have seen in the previous section that this factor should not be very important. Secondly, a linear combination of elementary sources was used for NEAMTHM18 , while in the current study we use direct simulations for each scenario considered. The third reason is that the bathymetric and topographic data and applied grid resolution within the system of nested grids are very different. It should be noted that the NEAMTHM18 hazard curves include log-normal uncertainty (Glimsdal et al., 2019). This treatment is meant to be a pragmatic approach for dealing with uncertainty stemming from potentially inaccurate source and inundation modeling, including DEM inaccuracy. This uncertainty treatment tends to increase the hazard relative to the point-value hazard curves which would be obtained without it, because both a bias and a dispersion are taken into account (see also Davies et al., 2018, their equation 15). We note that a similar observation is found for probabilistic seismic hazard, when the sigma of the ground motion prediction equation is increased (e.g., Bommer and Abrahamson, 2006). Figure 5A shows the locations of the two offshore Points of Interest (PoIs) from the NEAMTHM18 assessment closest to the town of Catania. Separated by approximately 20 km along the 50 m depth isobath, the two PoIs have very different locations with respect to the coastline. The northernmost point (black), located on steep seafloor, is less than 1 km from the shoreline. The southernmost point (red) is approximately 6 km offshore on a gentler slope facing the coastline. The topo-bathymetric contrast between these two sites alone provides motivation to improve the resolution of the hazard assessment. The curves in Figure 5B display offshore maximum elevation Probability of Exceedance (PoE) for a 50-years interval. The red and black curves display the NEAMTHM18 model uncertainty where the 2 nd , 50th, and 98th percentiles are estimated for 50 m water depth from the onshore MIH curves using Green's Law. Consistently higher offshore surface elevations are obtained for the northernmost PoI than for the southernmost PoI.
From the HySEA simulations carried out here, we stored time-series of wave height at each of the blue squared points indicated in Figure 5A, which are spaced at 2 km intervals along the 50 m isobath. The blue lines in Figure 5B show PoE curves for this refined set of PoIs. The new offshore PoE curves lie within the uncertainties estimated from NEAMTHM18 MIH curves up to a maximum offshore wave height of approximately 4 m. Above this level, the newly estimated PoE curves lie below the NEAMTHM18 estimates yet partly within the second percentile of the hazard corresponding to PoI 2. We stress that the (blue) hazard curves in Figure 5B, obtained from the local analysis, are not subject to the same uncertainty as the (red and black) curves derived from normalized NEAMTHM18 MIH percentiles. From this comparison, it is evident that the difference grows with growing intensity, as expected from the log-normal uncertainty. We then suggest that the divergence of the local hazard curves is strongly controlled by the uncertainty treatment. To better address our claim, we provide a more direct comparison, limited to the single PoI indicated in Figure 6. Here, the single mean local hazard curve (blue) is displayed together with both the NEAMTHM18 curves derived from the onshore amplified MIH (red) and NEAMTHM18 percentiles based on offshore point values (black/green). These point-value curves do use the NEAMTHM18 sources and linear combinations, but do not use the amplification factors and associated uncertainty treatment. Whereas the red and the blue curves diverge for the higher tsunami surface elevations (lower probabilities), with the NEAMTHM18 hazard curves indicating a higher PoE for higher offshore surface elevations, the local hazard curve and the point-value NEAMTHM18 curves follow the same trend. We conclude that the NEAMTHM18 model and the local hazard analysis are relatively consistent, despite clear discrepancies on the mean curve that can possibly stem from differences in bathymetric data and spatial resolution as well as the linear combination vs. distinct specification of sources.
While the uncertainty in MIH curves from the amplification factors is neglected, the use of NLSW models reduces the overall uncertainty. However, a significant uncertainty remains due to factors such as NLSW approximations, uncertainties associated with topo-bathymetric data, assumptions in friction modeling, etc. This uncertainty likely has a big impact in the tail of the distribution but, to the best of our knowledge, a method for treating it consistently has not been yet developed. This goes beyond the goals of this paper. However, if site-specific hazard for higher intensities is required, for example due to the presence of a critical infrastructure, then an extended analysis of this uncertainty would be required, as clearly demonstrated in Figure 5B.

Inundation and Coastal Zone Hazard Results
Figure 7 displays, as an example, the MIH for a single simulation out of the many considered for the full hazard assessment, for the entire region of the inner grid (left panel). In addition, close-up views for near Catania harbor and on the plains south of the city are also shown (right panels). Clearly, the 10 m grid in the simulation resolves the distribution of inundation heights at a scale of detail of relevance to onshore infrastructures. For each scenario, Tsunami-HySEA stores the maximum inundation height over the co-seismically deformed topography achieved at each grid location. Additional metrics such as maximum current velocities and momentum fluxes, the latter measuring the maximum of the instantaneous products of H and U, are also stored. In the simulations conducted here, it took about 25 min to run a single simulation. However, up to 1,024 simulations could be run synchronously in the Tier-0 system Marconi-100, which meant that we were able to produce more than 2000 scenario simulations during just 1 h. Figures 8-10 shows the aggregated probability maps using the local PTHA workflow for the city of Catania using three different metrics, namely the flow depth, maximum surface elevation, and the maximum instantaneous momentum flux. These give a first rough overview of the hazard products that can be obtained using the workflow. In each probability map, we assume a 50-years exposure time, and the average hazard with respect to epistemic Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 uncertainty (percentile values are presented in a companion paper). Evidently, there is a significant local variability in the flow depth pattern, a clear indication that the local-scale hazard assessment adds much more information, more fine-grained through the utilization of details in the topography, than any regional assessment (using just the two offshore points displayed in Figure 5 to present the overall hazard). The maps not only clearly differentiate between the more hazardous areas close to the shoreline, and the less hazardous areas located inland, but also resolve interesting differences along the shoreline. We show the probability of exceeding flow depths of 1 cm, 1 m, 3 m, and 5 m in Figure 8. In the south, the PoE-50-years map for 1 m flow depth covers a rather large area extending up to 3-4 km inland, which is largely due to a rather flat and low-lying topography. The inner points have the smallest mapped PoE-50years of the order of 10 −5 , corresponding to an average return period T of about half a million years. If dealing with critical infrastructures, the need for looking at even longer return periods may arise.
For evacuation mapping design purposes, a flow depth below 1 m should be considered, given that a flow only a few tens of centimeters deep is considered sufficient in certain circumstances for dragging people away (Takagi et al., 2016). In the north, the flow depths do not extend as far inland due to steeper topography. PoE-50-years values exceeding 0.01 (T ∼ 5000 years) are mostly confined within 500-1,000 m distances from the shoreline. Moreover, we see clearly that the highest likelihood of exceeding 3 m flow depth is confined to a small strip along the coastline. Note that the flow depth displayed in Figure 8 is evaluated at locations for which the deformed topography is greater than zero. Hence, co-seismic displacement was dominantly positive, that is an uplift with respect to the preevent coastline, likely caused by near-field sources which limited the inundation. Such an effect could not be appreciated by the FIGURE 7 | Inundation visualizations at different scales at and close to Catania, Sicily, for a single scenario HySEA simulation to display flow depth in relation to infrastructure and communication routes. The maximum wave height and flow depth are displayed offshore and onshore, respectively. All imagery from and available via Google Earth.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 regional assessment not involving local scale inundation simulations for each considered seismic scenario.
In Figure 9, we zoom in on the Catania harbor area, and display PoE-50-years maps for exceeding maximum elevations of 1 m, 3 m, and 5 m, including both the maximum surface elevations (offshore) and the MIH (onshore). The different panels show again the relatively large heterogeneity in PoE depending on the location. Interesting results are obtained for the harbor and its vicinity. By taking a closer look to the panels of Figure 9, we note that the spatial distribution of the hazard depends, to first order, on the intensity threshold considered. The hazard is higher inside the harbor and lower outside it for the Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 13 lowest threshold considered of 1 m ( Figure 9A). This might be an indication of harbor resonance under the forcing of tsunami oscillations related to relatively small magnitude and likely local earthquake sources, characterized in turn by relatively small characteristic periods. Conversely, for the largest intensity of 5 m ( Figure 9C), the largest hazard is found offshore and north of the breakwaters, and south of the breakwaters. The breakwater and harbor orientation likely have in this case a clear positive effect in reducing the hazard inside the harbor, as PoEs are smaller here relative to the surrounding areas. Local extrema close to breakwater tips may be the indication of a tendency for the creation of horizontal eddies (e.g., Volpe et al., 2019). An intermediate more homogeneous situation is observed for intermediate intensity ( Figure 9B).
The PoE-50-years map for the maximum momentum flux is shown in Figure 10. This impact metric is more sensitive to the position than the two other metrics displayed as it involves the current speed non-linearly. It seems to differ with respect to the spatial distribution too, having higher values in the northern part of the study region. Moreover, we see that the momentum fluxes have larger values inside the harbor, which may have implications for boats and other offshore objects.
To better investigate the differences between the present local-scale and the previous NEAMTHM18 regional-scale FIGURE 9 | Probability of exceedance within 50 years of maximum surface elevation (offshore) and MIH (onshore) of (A) 1 m (B) 3 m, and (C) 5 m, for Catania and neighboring areas. The zero-elevation contour line for the undeformed topography is shown in blue.
FIGURE 10 | Probability of exceedance within 50 years of maximum momentum flux of (A) 50 m 3 /s 2 (B) 120 m 3 /s 2 , and (C) 300 m 3 /s 2 , for Catania and neighboring areas. The zero-elevation contour line for the undeformed topography is shown in blue.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 assessment, and to highlight the improvement in the spatial representation allowed by the site-specific analysis, Figure 11 compares the hazard curves obtained from the two analyses. Figure 11 displays the mean PoE-50-years hazard curves for MIH and offshore surface elevations at all the grid points along six different transects (A-F) in the innermost computational domain, compared with the NEAMTHM18 hazard curves derived using MIH values (related to the red curves in Figures 5 and 6) at the two PoI closest to Catania (coined PoI1 and PoI2, see also Figure 5). For both of these PoIs, we show the mean hazard curves as well as their 2nd and 98th percentile values. The local hazard curves are colored according to their topographic and bathymetric elevation/ depth values. The regional hazard curves for PoI1 are displayed using red curves, with the thick solid line representing the mean value and the thin lines the percentiles. For PoI2, the black dashed curve shows the mean value and the dash-dotted curves the percentiles. By investigating the different transects, we see that it may be possible to separate the results roughly into families of curves: FIGURE 11 | Comparison of surface elevation (MIH onshore) between local and regional assessments. The red and black lines indicate percentiles as indicated for MIH estimated from offshore heights using amplification factors (Glimsdal et al., 2019) and the local hazard curves are displayed for every grid node along the transects indicated. Colors indicate the topographic elevation of the point as indicated: blue lines are offshore, green lines are low elevation onshore, and brown lines indicate higher elevations. Since we display maximum surface elevation (relative to the baseline sea level), and not flow depth, the hazard curve is flat until it reaches the local topographic height and only starts to be meaningful at values greater than this. Further details are provided in the text.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 591549 the northern ones (A), the ones in the central part of the domain (B-D), and the southern ones (E,F). Transect A intersects the harbor and hence the curves have rather different characteristics than the other transects. Up to at least 5 m, the nearshore hazard values (blue lines) resemble the mean value of PoI1, while they are clearly larger than for PoI2. For the largest MIHs (above 5-10 m), the local hazard curves drop much faster than in the regional assessment, as already commented above. The distribution of offshore maximum surface elevation gives a narrower height distribution than the regional 98th and 2nd percentile values. Onshore (green to brown lines), however, the variability is much larger. The probability for high inundation heights in the vicinity of the harbor infrastructure is also seen in Figure 9, with (dark green) onshore curves sometimes exceeding the 98th percentile. Figure 10 illustrates the high momentum flux into the harbor; the resulting influx of water is trapped by the harbor wall and inundates the town. Such an effect would not be predicted in the regional assessment. The sensitivity of the hazard as a function of location, in relation to topographic details and coastal infrastructure, is only resolved in the local assessment. On the other hand, the regional hazard is averaged over coastal features, and so exceedances in few critical points are expected.
For transects B-D, we see that the nearshore hazard curves follow the same trend as for the transect A curves, closely following the mean value for PoI1 from NEAMTHM18 up to 5 m height (MIH or maximum surface elevation). Along these transects, the MIH hazard curves for the onshore points close to the shore do not exceed the offshore values, but rather seem to be in the same order of magnitude (the onshore curves are plotted on top of the offshore ones). For the largest MIHs, the discrepancies between the regional and local hazard curves are even larger than for transect A, with no MIH value exceeding 10 m even for the smallest PoE-50-years investigated (10 −6 ). However, based on the comparisons shown in Figure 6, this discrepancy is expected, as the uncertainty based NEAMTHM18 hazard curves provide larger probabilities at the highest intensities also for the offshore values. As we move southwards where the topography flattens, we see that the offshore curves are distributed over a larger probability range, and the probability depends mainly on the distance from the shoreline.
For curves E,F, we see a similar trend as for curves B-D, but here the offshore hazard curves are clearly higher than the mean hazard curves from NEAMTHM18 up to 3-5 m. In fact, the largest mean hazard curves from the local assessment along these transects correspond to the 98th percentile. However, there seems to be an opposite tendency onshore, where the tsunami hazard in the tail of the curves (low probability/high intensity) seems to be clearly lower than in the regional assessment. Again, we stress that this is expected due to offsets in the offshore hazard curves shown in Figure 6.
In summary, there is clear spatial variability in the hazard. Moreover, the local hazard seems to provide lower hazard estimates than regional assessment for the highest hazard intensities. We interpret the higher hazard for the highest intensities in the regional analysis as being due to the inclusion of uncertainties, which are not included in the local analysis. The spatial representation in the local hazard curves is superior to that in the regional assessment. This is both due to the spatial variability induced by applying a proper inundation model over realistic topo-bathymetric data, but also due to the more sophisticated model used to simulate the inundation, which includes in this case the co-seismic displacement associated with each scenario. The latter may compound the uncertainty treatment to increase the differences between the regional and site-specific assessments. The comparison with the NEAMTHM18 hazard curves show that the hazard can be locally even higher than the 98th percentiles of the regional assessment. This is not unexpected, as the amplification factors developed for NEAMTHM18 averaged the uncertainty treatment over a wide range of inundation sites within the NEAM region. Hence, very localized characteristics for a single site could not be resolved.
Topographic effects and friction have dramatic effects on the local inundation and overlooking the uncertainty on modelling them may reduce the hazard dramatically onshore. This is underpinned by the fact that the trend in the offshore hazard is opposite in the local analysis (with larger offshore hazard towards the south) than the NEAMTHM18 PoIs with larger hazard for PoI1 than for PoI2. The areas in the southern part of the domain have a gentler slope, which should exaggerate shoaling and hence increase the local hazard. On the other hand, the friction may act to reduce the tsunami amplitude during inundation, which may explain why the hazard for the inundated points are clearly lower for the largest MIH (and smallest probabilities). A second explanation for this discrepancy is the already observed prevalence of the coseismic uplift from local sources. Altogether, more sensitivity studies would be necessary to quantify the most influential factors of inundation variability.

CONCLUDING REMARKS
We have designed and implemented a workflow for site-specific Probabilistic Tsunami Hazard Assessment (PTHA) using High Performance Computing (HPC) to conduct tens of thousands of numerical tsunami simulations. To the best of our knowledge, this is the first occasion on which a PTHA has been performed for both a comprehensive discretization of seismic sources and highresolution inundation calculations. To limit the number of scenarios calculated to a manageable number, we used the NEAMTHM18 source discretization and regional offshore PTHA (part of the TSUMAPS-NEAM project: Basili et al., 2019). A hazard disaggregation was performed and selecting those sources expected to constitute 99% of the total hazard for generating a maximum inundation height in the range 1-4 m for Catania, Sicily, resulted in a total of 32,363 scenarios. This range was selected as that of most likely relevance to evacuation planning. We used the Tsunami-HySEA program to model tsunamigenesis, open sea propagation, and coastal inundation using a system of nested grids with 10 m-spacing at the finest resolution. A complete Tsunami-HySEA tsunami simulation for Catania with the grid system presented here took approximately We compare directly offshore maximum surface elevation in the new calculations with the regional NEAMTHM18 assessment and find that the Probability of Exceedance curves for multiple locations offshore Catania fall comfortably within the uncertainty range indicated in the regional assessment for maximum inundation heights up to about 5 m. Above 5 m, the PoE curves are clearly lower than the NEAMTHM18 curves, likely related to the fact that the regional hazard methodology inherits uncertainties that are not present in our local analysis. This leads to increased probabilities for the highest MIH values in the regional assessment in comparison with the local hazard model presented here. However, existing uncertainty in the local analysis is neglected, preventing an effective evaluation of the robustness of this drop in the hazard tails. Therefore, to better interpret this drop, future studies should be focused on the development of methods to efficiently treat the uncertainty in tsunami generation, propagation and inundation also in local studies.
We display hazard maps for flow depth, maximum inundation height, and momentum flux for the Catania region, illustrating the higher resolution possible in the local PTHA. We find the local hazard curves, aggregated from the Tsunami-HySEA inundation maps, to correspond well with the uncertainty range indicated by the offshore based NEAMTHM18 hazard curves. However, the details in the local PTHA are far better resolved than in the regional PTHA and there are noteworthy deviations from the uncertainties predicted from the regional analysis. For many parts of the low-lying coastal region, the PoE curves are lower than the NEAMTHM18 curves, likely due to the more advanced inundation model. In the vicinity of the harbor and downtown Catania, the local hazard curves significantly exceed the regional curves, due to the presence of the harbor wall: effects that could not be resolved in the regional assessment.
This study has provided a comprehensive demonstration of a local PTHA workflow and analysis for a single coastal region with simulation of the vast number of scenarios made possible through parallel GPU computations on the pre-exascale MARCONI-100 machine. Improved accuracy in the tsunami hazard can be achieved both with regards to the inundation modeling and the discretization of sources. Increasing the resolution and/or dimensions at either end will increase the computational demands.
With regards to inundation, the Tsunami-HySEA program can handle arbitrarily many systems of nested grids in the same simulation. For example, we have run calculations with the nested grid system displayed here for Catania together with a corresponding system of nested grids for Siracusa (Sicily), further down the coast. To constitute 99% of the hazard for Catania for the selected hazard thresholds then 32,363 scenarios are deemed. A similar analysis for Siracusa found 32,514 scenarios which were deemed to constitute 99% of the hazard for that region. However, many of these scenarios are common to the hazard at both sites, and a total of 42,720 scenarios are needed to cover this specification of the hazard for both locations. However, were a third, fourth, or fifth site, further afield, to be added to the PTHA target region, the overlap of scenarios would likely be smaller and there is likely a trade-off in the efficiency of covering multiple stretches of coastline in the same simulations (the open sea tsunami propagation accounts for a small fraction of the total simulation time; the greatest part is spent calculating the inundation).
With regard to the accuracy of inundation, we stress that uncertainty in the inundation simulation is not accounted for here. Significant uncertainty may arise from the topo-bathymetric model, the spatially uniform approach to model friction, as well as from the nonlinear shallow water approximation. These uncertainty sources may significantly affect the entire hazard curve.
Regarding the accuracy of the hazard with regards to the source discretization, we can refine the source specifications of both the so-called PS and BS seismic scenarios. The BS scenarios have discrete sets of strike, dip, and rake angles, in addition to the hypocenter fault dimensions, and slip. This is a high-dimensional parameter space with a refinement factor in each parameter combining multiplicatively to the total number of BS scenarios required. We see that the vast majority of the BS scenarios dominating the local hazard are very concentrated close to the shore where the inundation is considered. For the PS scenarios, with the significance attached to the pattern of heterogeneous slip, five slip realizations were employed for each geometrical earthquake hypothesis. Would 10, or 25, or 50 provide a useful refinement of the source discretization? A sensitivity analysis would be required to assess the effect of all changes in the source.
A truly comprehensive PTHA for a given site would also need to consider other source mechanisms for tsunamis, including landslide and meteorological sources, together with realistic estimates of source probabilities. For non-seismic sources, e.g., landslides, volcanoes and meteotsunamis, PTHA is still in its infancy with just a handful of applications. Examples include for instance Grezio et al. (2012) (2020) (multiple hazards). The methods applied for these sources lack the methodological complexity that earthquake PTHA has, and often diversity and lack of source frequency data make estimation of annual probability more uncertain and less constrained. Nevertheless, models for landslide tsunamis exist, including models such as BingClaw (Løvholt et al., 2017;Kim et al., 2019). Moreover, two codes of the HySEA family, Landslide-HySEA and Multilayer-HySEA (Macías et al., 2015;González-Vida et al., 2019, Macías et al., 2020c, 2020d, are also available. Related to HySEA, the numerical simulation of meteotsunamis will require the implementation of the pressure term that takes into account the key generating mechanism for these events.
This study is a demonstration of a workflow and is not an operational hazard assessment. The necessity of further investigating uncertainty quantification related to numerical modelling has been clearly highlighted. Nevertheless, we have first and foremost demonstrated that PTHA with high resolution inundation calculations is now within reach using modern HPC resources.

DATA AVAILIABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SG prepared the manuscript, assisted with the experiment design, performed analysis of the results, and generated illustrations and plots. SL contributed to conceiving the workflow, to design the experiment, to the post-processing of the simulation results, and to the interpretation of the hazard results. JM contributed to the workflow and in particular to the tsunami code development, to the numerical simulations and to code improvements for PTHA applications. FL contributed to the experiment design, analysis and interpretation of the results, and preparation of the manuscript. JS contributed to conceiving the workflow, to design the source selection, the post-processing, and the hazard aggregation, as well as to design of the case study and to the interpretation of the hazard results. MVol contributed to conceiving and implementing the workflow, performed the disaggregation, contributed to the tsunami numerical simulations setup, to the numerical simulations, to the postprocessing of the simulation results and to the hazard aggregation calculations. CL performed the numerical simulations and sanity checks, 3D visualization and data management. AB contributed with code development for tsunami generation. MC contributed to the development of the tsunami code, the implementation of the new nested mesh algorithm and output compression. MA contributed to the development of the tsunami code, and to the implementation of the new nested mesh algorithm and output compression. FR contributed the tsunami numerical simulations setup, to the post-processing of the simulation results and to the hazard aggregation calculations. SGl and MVög prepared codes and tested performance and implementation of the HPC workflow. AC and AS contributed to the slip distribution simulation setup and computation for the PS scenarios. RT contributed to implement the workflow and to conceive the disaggregation step. ML and MN contributed to implementation of the workflow. BB and LP prepared the Digital Elevation Model and numerical grids used for numerical modeling. All the authors contributed to the article and approved the submitted version.