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

^{1}The Norwegian Geotechnical Institute (NGI), Oslo, Norway^{2}Istituto Nazionale di Geofisica e Vulcanologia (INGV), Rome, Italy^{3}Departamento de Análisis Matemático, Estadística e Investigacíon Operativa y Matemática Aplicada, Facultad de Ciencias, Universidad de Málaga, Málaga, Spain^{4}Istituto Nazionale di Geofisica e Vulcanologia (INGV), Bologna, Italy^{5}GFZ German Research Centre for Geosciences, Potsdam, Germany^{6}CINECA SuperComputing Applications and Innovation, Rome, Italy^{7}Department of Physics “Ettore Pancini”, University of Naples “Federico II”, Naples, Italy

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.

## 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ânoğlu, 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 high-resolution inundation simulation typical of site-specific 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 (Løvholt et al., 2019).

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 topo-bathymetric 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 (Molinari et al., 2016). 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.

## Implementation of a High-Performance Computing Oriented Seismic Probabilistic Tsunami Hazard Analysis Workflow

### 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.

**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.

### 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 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 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)+ŋ(x,t), where H is the bathymetry, h the total water depth, and ŋ 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 (Løvholt et al., 2019). 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 wave-propagation, 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 inner-workflow 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.

**FIGURE 2**. Overview of simulation process, from source models and seafloor displacements **(left)** for BS and PS sources, the nested grid procedure **(center)**, and tsunami inundation **(right)** for different scenarios. The Tsunami-HySEA code simulates the flow depth on the co-seismically deformed bathymetry which is different from scenario to scenario. The aggregated hazard takes this into account. See Scala et al. (2020) for a detailed description of the stochastic slip distributions for the Predominant Seismicity sources.

### Hazard Aggregation

The results from all the simulations are combined in the hazard aggregation step. This consists also of the post-processing 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 log-normal, 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.

## Setup for Hazard Analysis Toward the City of Catania

### Computational Grids and Hydrodynamic Parameters

The Tsunami-HySEA simulations use three levels of nested local grids, as well as one global 0-grid for the open ocean propagation covering the Mediterranean Sea. The finest grid resolution is 10 m, and a spatial refinement ratio of four is used here. Hence, the other grids have resolutions 40 m, 160 m, and 640 m respectively. Their extents are displayed in Figure 2. The grid with the coarsest resolution (640 m) uses the open GEBCO topo-bathymetry model (https://www.gebco.net/), which was resampled. Regarding the finest grid, a 10 m DEM model was built interpolating the following datasets:

• LiDAR inland points with 2 m resolution. The Geology and Geotechnologies Laboratory (INGV, http://istituto.ingv.it/it/36-laboratori/1656-laboratorio-geologia-e-geotecnologie.html) provided the LiDAR dataset in the framework of an agreement with the Ministry of the Environment, Earth and Sea (http://www.pcn.minambiente.it) - Italian National Geoportal, owner of the data).

• EMODnet Digital Bathymetry (DTM from EMODnet Bathymetry Consortium (2018) http://doi.org/10.12770/18ff0d48-b203-4a65-94a9-5fd8b0ec35f6, which has a resolution of nearly 112 m.

• EU-DEM, EU-DEM-4258: 1 arcsec - five arcsec, EU-DEM-3035: 25 m, Color shaded DEM derived from the EU-DEM-3035: 25 m “Produced using Copernicus data and information funded by the European Union - EU-DEM layers.”

• MaGIC data, foglio 32 e foglio 33, data from the MaGIC project, Dipartimento della Protezione Civile (http://dati.protezionecivile.it/geoportalDPC/catalog/main/home.page).

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.

**FIGURE 3**. Overview of the hazard disaggregation/scenario selection procedure **(A)** trade-off curve between number of scenarios contributing to the hazard at the site and percentage of hazard reproduction for the 1–4 m interval. The scenarios are ranked such that those contributing the most come first **(B)** The 1–4 m interval hazard curve reproduced with the selected scenarios as well as the complete one.

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.

**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.

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.

## Results

### 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 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 (Molinari et al., 2016), 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.

**FIGURE 5**. Locations of PoIs offshore Catania **(A)** from the current study (blue squares, with 2 km spacing) and NEAMTHM18 (the red and black circles only, with 20 km spacing). Comparison of Probability of Exceedance curves **(B)** for the PoIs indicated on the map. The blue lines in panel **(B)** are the aggregated hazard curves from the new simulations at each of the PoIs marked with blue squares in panel **(A)**. The red and black curves in panel **(B)** are offshore PoE derived from the corresponding NEAMTHM18 hazard curves with the offshore height calculated from the MIH using Green’s Law. The PoIs labeled POI1 (black) and POI2 (red) have geographical coordinates 37.537^{o}N 15.133^{o}E and 37.355^{o}N 15.158^{o}E respectively.

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.

**FIGURE 6**. Comparison of offshore PoE curves for the location 37.355^{o}N 15.158^{o}E (50 m depth). Panel **(A)** shows the location of this PoI with respect to the coastline and panel **(B)** displays the new mean offshore hazard curve (blue), the second, 50th, and 98th percentiles of offshore hazard estimated from the NEAMTHM18 MIH curves using Green’s Law (red), and the corresponding offshore PoE curves calculated from the offshore point values (black/green).

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.

**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.

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 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.

**FIGURE 8**. Probability of exceedance within 50 years of flow depth of **(A)** 1 cm **(B)** 1 m **(C)** 3 m, and **(D)** 5 m, for Catania and neighboring areas. The zero-elevation contour line for the undeformed topography is shown in blue.

**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.

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-50-years 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 pre-event coastline, likely caused by near-field sources which limited the inundation. Such an effect could not be appreciated by the 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 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 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: the northern ones (A), the ones in the central part of the domain (B–D), and the southern ones (E,F).

**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.

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 high-resolution 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 25 min on the MARCONI-100 supercomputer at CINECA, meaning that approximately 13,600 GPU hours are required to complete the set of scenarios selected. With 1,024 simulations able to run in parallel, around 14 h of clock-time are needed to perform these calculations.

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), Geist and Lynett (2014), SipkinBeck et al. (2016), Løvholt et al. (2020) (landslides), Paris et al. (2019) (volcanoes), Geist et al. (2014) (meteotsunamis), and Grezio et al. (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 post-processing 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.

## Funding

This work is partially funded by the European Union’s Horizon 2020 Research and Innovation Program under grant agreement No 823844 (ChEESE Center of Excellence, www.cheese-coe.eu). Computing resources at the Barcelona Supercomputer Center (BSC), Spain, were used under grant AECT-2020-1-0009. Optimization and testing of the Tsunami-HySEA computational engine was partly performed at Piz Daint at CSCS (The Swiss National Computer Center) under the PRACE preparatory access N. 2010PA4869. Some of the authors also benefited from the financial support from the project “Assessment of Cascading Events triggered by the Interaction of Natural Hazards and Technological Scenarios involving the release of Hazardous Substances”, funded by the Italian Ministry MIUR PRIN (Progetti di Ricerca di Rilevante Interesse Nazionale) 2017 – Grant 2017CEYPS8. Initial tests done in common with PD8 for the PRACE award 2019215169 TSU-CAST. Additional prototyping and testing was performed on resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer (AA) declared a past co-authorship with several of the authors (SL, FR, MV) to the handling editor.

## Acknowledgments

We are grateful for assistance from BSC (Spain), CINECA (Italy), CSCS (Switzerland), and UNINETT Sigma2 (Norway), in implementing and executing the workflows and codes on the different HPC resources. Many of the figures in this paper are generated using GMT software (Wessel et al., 2019).

## References

Argyroudis, S. A., Fotopoulou, S., Karafagka, S., Pitilakis, K., Selva, J., Salzano, E., et al. (2020). A risk-based multi-level stress test methodology: application to six critical non-nuclear infrastructures in Europe. *Nat. Hazards* 100, 595–633. doi:10.1007/s11069-019-03828-5

Bao, H., Ampuero, J.-P., Meng, L., Fielding, E. J., Liang, C., Milliner, C. W. D., et al. (2019). Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake. *Nat. Geosci.* 12, 200–205. doi:10.1038/s41561-018-0297-z

Basili, R., Brizuela, B., Herrero, A., Iqbal, S., Lorito, S., Maesano, F. E., et al. (2018). NEAM Tsunami Hazard Model 2018 (NEAMTHM18): online data of the Probabilistic Tsunami Hazard Model for the NEAM Region from the TSUMAPS-NEAM project. Istituto Nazionale di Geofisica e Vulcanologia (INGV). Available at: http://doi.org/10.13127/tsunami/neamthm18

Basili, R., Brizuela, B., Herrero, A., Iqbal, S., Lorito, S., Maesano, F. E., et al. (2019). NEAMTHM18 documentation: the making of the TSUMAPS-NEAM tsunami hazard model 2018. Istituto Nazionale di Geofisica e Vulcanologia (INGV). doi:10.5281/ZENODO.3406625

Basili, R., Tiberti, M. M., Kastelic, V., Romano, F., Piatanesi, A., Selva, J., et al. (2013). Integrating geologic fault data into tsunami hazard studies. *Nat. Hazards Earth Syst. Sci.* 13, 1025–1050. doi:10.5194/nhess-13-1025-2013

Bazzurro, P., and Cornell, C. A. (1999). Disaggregation of seismic hazard. *Bull. Seismol. Soc. Am.* 89, 501–520.

Behrens, J., and Dias, F. (2015). New computational methods in tsunami science. *Philos Trans A Math Phys Eng. Sci.* 373, 20140382. doi:10.1098/rsta.2014.0382

Bommer, J. J., and Abrahamson, N. A. (2006). Why do modern probabilistic seismic-hazard analyses often lead to increased hazard estimates?. *Bull. Seismol. Soc. Am.* 96, 1967–1977. doi:10.1785/0120060043

Bricker, J. D., Gibson, S., Takagi, H., and Imamura, F. (2015). On the need for larger Manning’s roughness coefficients in depth-integrated tsunami inundation models. *Coast Eng. J.* 57, 1550005–1–1550005–13. doi:10.1142/S0578563415500059

Burbidge, D., Cummins, P. R., Mleczko, R., and Thio, H. K. (2008). “A probabilistic tsunami hazard assessment for western australia,” in *tsunami science four Years after the 2004 Indian ocean tsunami*. basel: birkhäuser basel, 2059–2088. doi:10.1007/978-3-0346-0057-6

Choi, B. H., Pelinovsky, E., Ryabov, I., and Hong, J. S. (2002). Distribution functions of Tsunami wave heights. *Nat. Hazards*. 25, 1–21. doi:10.1023/A:1013379705323

Cummins, P. R. (2019). Irrigation and the Palu landslides. *Nat. Geosci.* 12, 881–882. doi:10.1038/s41561-019-0467-7

Davies, G., Griffin, J., Løvholt, F., Glimsdal, S., Harbitz, C., Thio, H. K., et al. (2018). A global probabilistic tsunami hazard assessment from earthquake sources. *Geol. Soc. London, Spec. Publ.* 456, 219–244. doi:10.1144/SP456.5

de la Asunción, M., Castro, M. J., Fernández-Nieto, E. D., Mantas, J. M., Acosta, S. O., and González-Vida, J. M. (2013). Efficient GPU implementation of a two waves TVD-WAF method for the two-dimensional one layer shallow water system on structured meshes. *Comput. Fluids* 80, 441–452. doi:10.1016/j.compfluid.2012.01.012

Dipartimento della Protezione Civile (2018). “Presidenza del Consiglio dei Ministri - Dipartimento della Protezione Civile - decreto 2 ottobre 2018alle Strutture operative del Servizio nazionale di protezione civile per l’aggiornamento delle pianificazioni di protezione civile per il rischio maremoto, (18A07309) (GU Serie Generale n.266 del 15-11-2018)”. Available at: http://www.protezionecivile.gov.it/amministrazione-trasparente/provvedimenti/dettaglio/-/asset_publisher/default/content/indicazioni-alle-componenti-ed-alle-strutture-operative-del-servizio-nazionale-di-protezione-civile-per-l-aggiornamento-delle-pianificazioni-di-prot-1

Fritz, H. M., Borrero, J. C., Synolakis, C. E., Okal, E. A., Weiss, R., Titov, V. V., et al. (2011). Insights on the 2009 South Pacific tsunami in Samoa and Tonga from field surveys and numerical simulations. *Earth Sci. Rev.* 107, 66–75. doi:10.1016/j.earscirev.2011.03.004

Geist, E. L. (2002). Complex earthquake rupture and local tsunamis. *J. Geophys. Res.* 107, 2086. doi:10.1029/2000JB000139

Geist, E. L., and Parsons, T. (2006). Probabilistic analysis of tsunami hazards. Nat. *Hazards* 37, 277–314. doi:10.1007/s11069-005-4646-z

Geist, E. L., ten Brink, U. S., and Gove, M. (2014). A framework for the probabilistic analysis of meteotsunamis. *Nat. Hazards* 74 (1), 123–142. doi:10.1007/s11069-014-1294-1

Geist, E., and Lynett, P. (2014). Source processes for the probabilistic assessment of tsunami hazards. *Oceanography* 27 (2), 86–93. doi:10.5670/oceanog.2014.43

Glimsdal, S., Løvholt, F., Harbitz, C. B., Romano, F., Lorito, S., Orefice, S., et al. (2019). A new approximate method for quantifying tsunami maximum inundation height probability. *Pure Appl. Geophys.* 176, 3227–3246. doi:10.1007/s00024-019-02091-w

Goda, K., and De Risi, R. (2018). Multi-hazard loss estimation for shaking and tsunami using stochastic rupture sources. *Int. J. Disaster Risk Reduction.* 28, 539–554. doi:10.1016/j.ijdrr.2018.01.002

González, F. I., Geist, E. L., Jaffe, B., Kânoğlu, U., Mofjeld, H., Synolakis, C. E., et al. (2009). Probabilistic tsunami hazard assessment at Seaside, Oregon, for near- and far-field seismic sources. *J. Geophys. Res.* 114, C11023. doi:10.1029/2008JC005132

González-Vida, J. M., Macías, J., Castro, M. J., Sánchez-Linares, C., de la Asunción, M., Ortega, S., et al. (2019). The Lituya Bay landslide-generated mega-tsunami. Numerical simulation and sensitivity analysis,. *Nat. Hazards Earth Syst. Sci.* 19, 369–388. doi:10.5194/nhess-19-369-2019

Grezio, A., Babeyko, A., Baptista, M. A., Behrens, J., Costa, A., Davies, G., et al. (2017). Probabilistic tsunami hazard analysis: multiple sources and global applications. *Rev. Geophys.* 55, 1158–1198. doi:10.1002/2017RG000579

Grezio, A., Cinti, F. R., Costa, A., Faenza, L., Perfetti, P., Pierdominici, S., et al. (2020). Multisource bayesian probabilistic tsunami hazard analysis for the gulf of naples (Italy). *J. Geophys. Res. Oceans* 125 (2), e2019JC015373. doi:10.1029/2019JC015373

Grezio, A., Sandri, L., Marzocchi, W., Argnani, A., Gasparini, P., and Selva, J. (2012). “Probabilistic tsunami hazard assessment for Messina Strait area (Sicily, Italy)”. *Nat. Hazards* 64 (1), 329–358. doi:10.1007/s11069-012-0246-x

Griffin, J., Latief, H., Kongko, W., Harig, S., Horspool, N., Hanung, R., et al. (2015). An evaluation of onshore digital elevation models for modeling tsunami inundation zones. *Front. Earth Sci.* 3, 32. doi:10.3389/feart.2015.00032

Horspool, N., Pranantyo, I., Griffin, J., Latief, H., Natawidjaja, D. H., Kongko, W., et al. (2014). A probabilistic tsunami hazard assessment for Indonesia. *Nat. Hazards Earth Syst. Sci.* 14, 3105–3122. doi:10.5194/nhess-14-3105-2014

Kagan, Y. Y., and Jackson, D. D. (2013). Tohoku earthquake: a surprise?. *Bull. Seismol. Soc. Am.* 103, 1181–1194. doi:10.1785/0120120110

Kim, J., Løvholt, F., Issler, D., and Forsberg, C. F. (2019). Landslide material control on tsunami genesis—the Storegga Slide and tsunami (8,100 years BP). *J. Geophys. Res. Oceans* 124 (6), 3607–3627. doi:10.1029/2018JC014893

Lay, T., Kanamori, H., Ammon, C. J., Nettles, M., Ward, S. N., Aster, R. C., et al. (2005). The great sumatra-andaman earthquake of 26 december 2004. *Science.* 308, 1127–1133. doi:10.1126/science.1112250

Lorito, S., Selva, J., Basili, R., Romano, F., Tiberti, M. M., and Piatanesi, A. (2015). Probabilistic hazard for seismically induced tsunamis: accuracy and feasibility of inundation maps. *Geophys. J. Int.* 200, 574–588. doi:10.1093/gji/ggu408

Løvholt, F., Bondevik, S., Laberg, J. S., Kim, J., and Boylan, N. (2017). Some giant submarine landslides do not produce large tsunamis. *Geophys. Res. Lett.* 44 (16), 8463–8472. doi:10.1002/2017GL074062

Løvholt, F., Glimsdal, S., and Harbitz, C. B. (2020). On the landslide tsunami uncertainty and hazard. *Landslides* 17, 2301–2315. doi:10.1007/s10346-020-01429-z

Løvholt, F., Glimsdal, S., Harbitz, C. B., Zamora, N., Nadim, F., Peduzzi, P., et al. (2012). Tsunami hazard and exposure on the global scale. *Earth Sci. Rev.* 110, 58–73. doi:10.1016/j.earscirev.2011.10.002

Løvholt, F., Griffin, J., and Salgado-Gálvez, M. A. (2015). “Tsunami hazard and risk assessment on the global scale,” in *Encyclopedia of Complexity and systems science*. Berlin, Heidelberg: Springer, 1–34. doi:10.1007/978-3-642-27737-5_642-1

Løvholt, F., Lorito, S., Macias, J., Volpe, M., Selva, J., and Gibbons, S. (2019). “Urgent tsunami computing” in IEEE/ACM HPC for urgent decision making (UrgentHPC) Denver, CO, November 17, 2019 (IEEE), 45–50. doi:10.1109/UrgentHPC49580.2019.00011

Macías, J., Castro, M. J., and Escalante, C. (2020a). Performance assessment of the Tsunami-HySEA model for NTHMP tsunami currents benchmarking Laboratory data. *Coast. Eng.* 158, 103667. doi:10.1016/j.coastaleng.2020.103667

Macías, J., Castro, M. J., Ortega, S., and González-Vida, J. M. (2020b). Performance assessment of Tsunami-HySEA model for NTHMP tsunami currents benchmarking. Field cases. *Ocean Model.* 152, 101645. doi:10.1016/j.ocemod.2020.101645

Macías, J., Escalante, C., and Castro, M. J. (2020c). “Multilayer-HySEA model validation for landslide generated tsunamis. Part I Rigid slides” Development of efficient hydrodynamic and morphodynamic simulators for risk assessment and forecasting (SIMURISK), May 22–August 27, 2000. doi:10.5194/nhess-2020-171

Macías, J., Escalante, C., and Castro, M. J. (2020d). Multilayer-HySEA model validation for landslide generated tsunamis Part II Granular slides. *Nat. Hazards Earth Syst. Sci.* doi:10.5194/nhess-2020-172

Macías, J., Castro, M. J., Ortega, S., Escalante, C., and González-Vida, J. M. (2017). Performance benchmarking of tsunami-HySEA model for NTHMP’s inundation mapping activities. *Pure Appl. Geophys.* 174, 3147–3183. doi:10.1007/s00024-017-1583-1

Macías, J., Vázquez, J. T., Fernández-Salas, L. M., González-Vida, J. M., Bárcenas, P., Castro, M. J., et al. (2015). The Al-Boraní submarine landslide and associated tsunami. A modelling approach. *Mar. Geol.* 361, 79–95. doi:10.1016/j.margeo.2014.12.006

Molinari, I., Tonini, R., Lorito, S., Piatanesi, A., Romano, F., Melini, D., et al. (2016). Fast evaluation of tsunami scenarios: uncertainty assessment for a Mediterranean Sea database. *Nat. Hazards Earth Syst. Sci.* 16, 2593–2602. doi:10.5194/nhess-16-2593-2016

Newman, A. V., Hayes, G., Wei, Y., and Convers, J. (2011). The 25 October 2010 Mentawai tsunami earthquake, from real-time discriminants, finite-fault rupture, and tsunami excitation.*Geophys. Res. Lett.* 38, L05302. doi:10.1029/2010GL046498

Nosov, M. A., and Kolesov, S. V. (2007). Elastic oscillations of water column in the 2003 Tokachi-oki tsunami source: *in-situ* measurements and 3-D numerical modelling. *Nat. Hazards Earth Syst. Sci.* 7, 243–249. doi:10.5194/nhess-7-243-2007

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. *Bull. Seismol. Soc. Am.* 82, 1018–1040.

Omira, R., Dogan, G. G., Hidayat, R., Husrin, S., Prasetya, G., Annunziato, A., et al. (2019). The september 28th, 2018, tsunami in palu-sulawesi, Indonesia: a post-event field survey. *Pure Appl. Geophys.* 176, 1379–1395. doi:10.1007/s00024-019-02145-z

Paris, R., Ulvrova, M., Selva, J., Brizuela, B., Costa, A., Grezio, A., et al. (2019). Probabilistic hazard analysis for tsunamis generated by subaqueous volcanic explosions in the Campi Flegrei caldera Italy. *J. Volc. Geothermal Res.* 379, 106–116. doi:10.1016/j.jvolgeores.2019.05.010

Park, H., Wiebe, D. M., Cox, D. T., and Cox, K. (2014). Tsunami inundation modeling: sensitivity of velocity and momentum flux to bottom friction with application to building damage at seaside. *Oregon. Coast. Eng. Proc.* 1, 1. doi:10.9753/icce.v34.currents.1

Pitilakis, K., Argyroudis, S., Fotopoulou, S., Karafagka, S., Kakderi, K., and Selva, J. (2019). Application of stress test concepts for port infrastructures against natural hazards. The case of Thessaloniki port in Greece. Reliab. Eng. Syst. *Saf. Now.* 184, 240–257. doi:10.1016/j.ress.2018.07.005

Power, W., Downes, G., and Stirling, M. (2007). “Estimation of tsunami hazard in New Zealand due to south American earthquakes,” in *Tsunami and its Hazards in the Indian and pacific oceans*. basel: birkhäuser basel, 547–564. doi:10.1007/978-3-7643-8364-0-15

Power, W., Wang, X., Wallace, L., Clark, K., and Mueller, C. (2018). The New Zealand Probabilistic Tsunami Hazard Model: development and implementation of a methodology for estimating tsunami hazard nationwide. *Geol. Soc. London, Spec. Publ.* 456, 199–217. doi:10.1144/SP456.6

Scala, A., Lorito, S., Romano, F., Murphy, S., Selva, J., Basili, R., et al. (2020). Effect of shallow slip amplification uncertainty on probabilistic tsunami hazard analysis in subduction zones: use of long-term balanced stochastic slip models. *Pure Appl. Geophys.* 177, 1497–1520. doi:10.1007/s00024-019-02260-x

Selva, J., Tonini, R., Molinari, I., Tiberti, M. M., Romano, F., Grezio, A., et al. (2016). Quantification of source uncertainties in seismic probabilistic tsunami hazard analysis (SPTHA). *Geophys. J. Int.* 205, 1780–1803. doi:10.1093/gji/ggw107

Sepúlveda, I., Tozer, B., Haase, J. S., Liu, P. L. F., and Grigoriu, M. (2020). Modeling uncertainties of bathymetry predicted with satellite altimetry data and application to tsunami hazard assessments. *J. Geophys. Res. Solid Earth* 125, 1–25. doi:10.1029/2020JB019735

SipkinBeck, E. M., Mountjoy, J. J., Power, W. L., and Mueller, C. (2016). “Probabilistic hazard of tsunamis generated by submarine landslides in the cook strait canyon (New Zealand)” in *Global tsunami science: past and future, volume I. Pageoph topical volumes*. Editors E. L. Geist, H. M. Fritz, A. B. Rabinovich, and Y. Tanioka (Cham: Birkhäuser). doi:10.1007/978-3-319-55480-8_6

Synolakis, C., and Kânoğlu, U. (2015). The Fukushima accident was preventable. *Philos. Trans. A Math. Phys. Eng. Sci.* 373, 20140379. doi:10.1098/rsta.2014.0379

Takagi, H., Mikami, T., Fujii, D., Esteban, M., and Kurobe, S. (2016). Mangrove forest against dyke-break-induced tsunami on rapidly subsiding coasts. *Nat. Hazards Earth Syst. Sci.* 16, 1629–1638. doi:10.5194/nhess-16-1629-2016

Titov, V. V., and Gonzalez, F. I. (1997). Implementation and testing of the method of splitting tsunami (MOST) model. NOAA technical memorandum. Available at: https://repository.library.noaa.gov/view/noaa/10979/noaa_10979_DS1.pdf

Tolkova, E. (2007). Compression of MOST propagation database (NOAA technical memorandum OAR PMEL-134). Available at: https://repository.library.noaa.gov/view/noaa/13838/noaa_13838_DS1.pdf

Ulrich, T., Vater, S., Madden, E. H., Behrens, J., van Dinther, Y., van Zelst, I., et al. (2019). Coupled, physics-based modeling reveals earthquake displacements are critical to the 2018 Palu, sulawesi tsunami. *Pure Appl. Geophys.* 176, 4069–4109. doi:10.1007/s00024-019-02290-5

Volpe, M., Lorito, S., Selva, J., Tonini, R., Romano, F., and Brizuela, B. (2019). From regional to local SPTHA: efficient computation of probabilistic tsunami inundation maps addressing near-field sources. *Nat. Hazards Earth Syst. Sci.* 19, 455–469. doi:10.5194/nhess-19-455-2019

Keywords: tsunami, hazard, probabilistic tsunami hazard analysis, high-performance computing, gpu, inundation, earthquakes

Citation: Gibbons SJ, Lorito S, Macías J, Løvholt F, Selva J, Volpe M, Sánchez-Linares C, Babeyko A, Brizuela B, Cirella A, Castro MJ, de la Asunción M, Lanucara P, Glimsdal S, Lorenzino MC, Nazaria M, Pizzimenti L, Romano F, Scala A, Tonini R, Manuel González Vida J and Vöge M (2020) Probabilistic Tsunami Hazard Analysis: High Performance Computing for Massive Scale Inundation Simulations. *Front. Earth Sci.* 8:591549. doi: 10.3389/feart.2020.591549

Received: 04 August 2020; Accepted: 17 November 2020;

Published: 11 December 2020.

Edited by:

Joanna Faure Walker, University College London, United KingdomReviewed by:

Alberto Armigliato, University of Bologna, ItalyHyoungsu Park, University of Hawaii at Manoa, United States

Copyright © 2020 Gibbons, Lorito, Macías, Løvholt, Selva, Volpe, Sánchez-Linares, Babeyko, Brizuela, Cirella, Castro, De La Asunción, Lanucara, Glimsdal, Lorenzino, Nazaria, Pizzimenti, Romano, Scala, Tonini, Manuel González Vida and Vöge. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Steven J. Gibbons, steven.gibbons@ngi.no