Testing Tsunami Inundation Maps for Evacuation Planning in Italy

Inundation maps are a fundamental tool for coastal risk management and in particular for designing evacuation maps and evacuation planning. These in turn are a necessary component of the tsunami warning systems’ last-mile. In Italy inundation maps are informed by a probabilistic tsunami hazard model. Based on a given level of acceptable risk, Italian authorities in charge for this task recommended to consider, as design hazard intensity, the average return period of 2500 years and the 84th percentile of the hazard model uncertainty. An available, regional-scale tsunami hazard model was used that covers the entire Italian coastline. Safety factors based on analysis of run-up variability and an empirical coastal dissipation law on a digital terrain model (DTM) were applied to convert the regional hazard into the design run-up and the corresponding evacuation maps with a GIS-based approach. Since the regional hazard cannot fully capture the local-scale variability, this simplified and conservative approach is considered a viable and feasible practice to inform local coastal risk management in the absence of high-resolution hazard models. The present work is a first attempt to quantify the uncertainty stemming from such procedure. We compare the GIS-based inundation maps informed by a regional model with those obtained from a local high-resolution hazard model. Two locations on the coast of eastern Sicily were considered, and the local hazard was addressed with the same seismic model as the regional one, but using a higher-resolution DTM and massive numerical inundation calculations with the GPU-based Tsunami-HySEA nonlinear shallow water code. This study shows that the GIS-based inundation maps used for planning deal conservatively with potential hazard underestimation at the local scale, stemming from typically unmodeled uncertainties in the numerical source and tsunami evolution models. The GIS-based maps used for planning fall within the estimated “error-bar” due to such uncertainties. The analysis also demonstrates the need to develop local assessments to serve very specific risk mitigation actions to reduce the uncertainty. More in general, the presented case-studies highlight the importance to explore ways of dealing with uncertainty hidden within the high-resolution numerical inundation models, e.g., related to the crude parameterization of the bottom friction, or the inaccuracy of the DTM.


INTRODUCTION
Tsunami early warning systems (TEWS) play a significant role in protecting coastal areas from the tsunami threat. The first TEWS was developed for the Pacific Ocean after the tsunami following the Mw 8.6 Aleutian Islands earthquake on the April 1, 1946. Since then, the development of increasingly accurate and reliable TEWS was strongly fostered under the lead of United Nations international agencies, in collaboration with the Member States, in particular after the devastating December 26, 2004 Indian Ocean tsunami. These political efforts, together with rapid technological developments during recent decades, increased the number of TEWS worldwide and their monitoring capabilities, extending the coverage of tsunami monitoring activities to virtually all the oceans worldwide (Bernard and Titov, 2015;Angove et al., 2019;Mulia and Satake, 2020).
Within the framework of the Intergovernmental Coordination Group for the Tsunami Early Warning and Mitigation System in the North-eastern Atlantic, the Mediterranean and connected seas (ICG/NEAMTWS), the Italian Tsunami Alert Center at the Istituto Nazionale di Geofisica e Vulcanologia (CAT-INGV) started operating 24/7 from October 1, 2014 as Candidate Tsunami Service Providers (TSPs) for the NEAM region (North-eastern Atlantic, the Mediterranean and connected seas, Bernardi et al., 2015;Amato, 2020). CAT-INGV is the upstream component of the Italian Tsunami Warning System (SiAM; see also the national implementing decrees of SiAM directive, DPCM, 2017; Dipartimento della Protezione Civile, 2018), constituted also by the Italian Civil Protection Department (DPC) and the Italian Institute for Environmental Protection and Research (ISPRA). As the other NEAM TSPs (CENALT, France; NOA, Greece; KOERI, Turkey; IPMA, Portugal), CAT-INGV is committed to deliver messages for potential earthquakegenerated tsunamis, containing the alert levels and the estimated time of arrival at preselected forecast points along the coasts. The two alert levels established by the ICG/ NEAMTWS agency are advisory/orange (tsunami runup ≤ 1 m) and watch/red (tsunami runup >1 m). Once an alert message has been generated, it is disseminated by CAT-INGV to the IOC/UNESCO subscriber Member States, and by DPC over Italian territory. DPC also carries out specific actions related to the forecasted alert levels (Dipartimento della Protezione Civile, 2018).
An upper limit for the expected runup is not specified by the NEAMTWS protocols and alert messages, which only forecast a state of significant inundation with runup >1 m. Coastal risk managers must then define this limit by identifying the inundation-prone areas with evacuation maps. Since it is impossible to define a deterministic maximum inundation, an upper boundary must be set as the evacuation distance which has a certain probability of being exceeded in a given time frame, in principle corresponding to a given risk reduction level for the same vulnerability and exposure levels. The probabilistic hazard frameworks are ideal in this respect, as they provide estimates of the exceedance probability for given intensity measures and for specific time frames. Such an approach, which we can refer to as a "uniform hazard", is commonly adopted for seismic building codes. For example, the design seismic action is often set as the one having a probability of being exceeded of 0.1 (10%) or 0.02 (2%) in 50 years (corresponding respectively to 475 and 2,475 years Average Return Period, hereinafter ARP), mostly within a "multi-tier" approach related to different performance levels, e.g., ranging from life-safety to building collapse as a target (ASCE/SEI 7-16, 2017;NTC, 2018). The use of a homogeneous risk target is sometimes advocated: it would be more consistent with the final use of seismic design maps to adopt a "uniform risk" assumption, in which the design ground motions are defined to provide the same level of risk everywhere, e.g., the annual probability of collapse (e.g., Luco et al., 2007). More recently, an approach based on probabilistic tsunami hazard analysis has been adopted in tsunami building codes, considering different performance levels, similarly to what is done for seismic design, and assuming a design probability of 0.02 (2%) in 50 years corresponding to the longer ARP of 2,475 years for the collapse prevention performance level (Chock et al., 2016). Authorities and stakeholders can use the same principles for land-use design (e.g., Geist and Lynett, 2014).
A specific difficulty with Probabilistic Tsunami Hazard Analysis (PTHA) is that representing the seismic source variability can require up to millions of scenarios, especially in areas where the crustal intermediate-magnitude seismicity provides a relevant contribution (Selva et al., 2016). Although the use of numerical simulations is by far the most accredited tool to model tsunami propagation and inundation (Behrens and Dias, 2015), a massive use of tsunami simulations on highresolution topo-bathymetric grids still represents a challenge in several practical applications because of the computational cost. Inundation simulations are also subject to the availability of very accurate and detailed DTMs both on-and off-shore of the sites of interest which, even when available, may represent an important uncertainty source (e.g., Griffin et al., 2015;Song and Goda, 2019), also related to the choice of appropriate roughness values (e.g., Kaiser et al., 2011;Griffin et al., 2015). Therefore, the combination of numerical simulations offshore and rapid approximated methods onshore is a common practice as alternative to explicit numerical inundation simulations (Grezio et al., 2017, Grezio et al., 2020. The method proposed by the Italian DPC to local authorities to define evacuation maps (Dipartimento della Protezione Civile, 2018) embodies this philosophy: 1) the design probability for the watch/red alert is fixed through the selection of the 2,500 years ARP and 84th percentile of the epistemic uncertainty from a regional and relatively coarse-scale Seismic PTHA (SPTHA, in the current implementation is the NEAMTHM18, NEAM Tsunami Hazard Model 2018; Basili et al., 2018;Basili et al., 2019;Basili et al., 2021), evaluated off-shore and amplified to the coast; 2) the extension of the evacuation area is the outcome of a chain of simplified and conservative approximations, including some "safety factors". The details of the procedure are subsequently provided in this paper. Very similar approaches that inspired the Italian authorities have been adopted by New Zealand (Leonard et al., 2008;Fraser and Power, 2013;MCDEM, 2016) and Samoa (Wright et al., 2011). Defining inundation maps is not only important for the population to know where people have to evacuate in case of a tsunami alert, but it is also important for the local contingency planning, e.g., to identify tsunami-safe emergency management locations (local emergency centers and areas), further than for coastal planning in general.
The aim of the present work is to evaluate the performance of the Italian evacuation maps, as defined at the national level, in comparison with two more detailed case-studies at local scale. Starting from the same seismic model defined in the regional NEAMTHM18, a site-specific SPTHA, including massive high-resolution numerical inundation simulations (Gibbons et al., 2020), is here used to obtain the inundation maps for several different design levels, for two selected target sites in South-Eastern Sicily, Italy. These sites are the coastal zones in the vicinity and around the cities of Catania and Siracusa, two touristic, commercial and industrial urbanized areas, comprising harbour infrastructures and heterogeneous morphological terrain features (e.g., flat beaches, cliffs, jagged inlets), representing a test for a broad range of different coastal settings. This area has experienced destructive tsunamis in the past, as documented by both historical (Maramai et al., 2014;Maramai et al., 2019) and geological (De Martini et al., 2012;Polonia et al., 2017) evidence, with a rich and still open scientific debate on the causative sources of the most important events, in particular regarding the earthquakes occurred in 1693 and 1908 (Piatanesi and Tinti, 1998;Gutscher et al., 2006;Gerardi et al., 2008;Favalli et al., 2009;Pino et al., 2009;Billi et al., 2010;Tonini et al., 2011;Aloisi et al., 2013;Convertito and Pino, 2014;Ridente et al., 2014;Meschis et al., 2019).
The local SPTHA quantifies the exceedance probability as a function of the inland tsunami intensity (here, the maximum inundation height) on a regular grid of points at the simulation resolution (here, 10 m at the two target sites). This allows the extraction of high-resolution inundation lines corresponding to different ARPs and epistemic uncertainty percentiles, including the ones recommended for planning by the Italian DPC, for comparison with the evacuation maps derived from the simplified method. Uncertainties not explicitly modeled in the local SPTHA are also subsequently considered and estimated to first order. To discriminate between different uncertainty sources, the procedure used for the official inundation zones, originally applied to all the Italian coasts using a coarser DTM, is here applied to the same highresolution DTM used for the numerical inundation simulations, yet based on the same hazard input. This updated version of the inundation maps can be compared with those obtained through local SPTHA to deal with the method approximations, and with the official maps to directly address the impact of DTM uncertainties.

DEVELOPMENT OF INUNDATION MAPS
In this section, the basic elements defining the inundation maps are introduced. We describe the procedure developed to design the maps adopted for tsunami risk-management in Italy from regional SPTHA, as well as the methodology for developing inundation maps starting from site-specific high-resolution SPTHA. In the following, for convenience, any quantity derived from the Italian authorities' indications to define the inundation maps (i.e., the 2,500 years ARP and the 84th percentile) are labeled as "design" quantities.

The Tsunami Hazard Regional Model
The indications issued by the Italian DPC for evacuation maps (Dipartimento della Protezione Civile, 2018) are based on a regional SPTHA model. In their current implementation, the reference SPTHA is the NEAMTHM18 (Basili et al., 2018). The NEAMTHM18 was produced in the frame of the TSUMAPS-NEAM project (http://www.tsumaps-neam.eu/; Basili et al., 2019). The model is the first long-term SPTHA for the entire NEAM region, covering a very high number of seismic sources and estimating inundation probability at a relatively coarse regional scale.
In NEAMTHM18, well-constrained seismic sources (as subduction interfaces or major fault systems) and diffuse crustal seismicity are treated separately and the epistemic uncertainty is estimated using an ensemble modeling technique (Selva et al., 2016;Basili et al., 2021). Offshore sea surface elevations are obtained by linear combinations of precomputed elementary source scenarios , which are propagated with the 2D shallow-water tsunami numerical code Tsunami-HySEA (de la Asunción et al., 2013), up to water depths of 50 m along a set of points of interest (POIs) spaced at about 20 km from each other (a subset of these points in the Mediterranean Sea is shown in Figure 1A). Offshore sea surface elevations are then transformed into maximum inundation height (MIH) through a simplified amplification factor (AF) method accounting for amplitudes, periods and polarities of the incident waves (Glimsdal et al., 2019).
The final results consist of a set of hazard curves (probability of exceedance in 50 years vs. MIH) at each POI. This set describes the epistemic uncertainty on the hazard, and it is summarized through the mean, median, and different percentiles ( Figures  1B,C). NEAMTHM18 includes several sources of uncertainty, among which those related to the seismic recurrence modeling, to the simplified modeling of the seismic source, as well as to the tsunami modeling (generation, propagation and inundation). Noteworthy, several rather strong approximations exist, which can be considered appropriate for a regional-scale SPTHA. Full details can be found in Basili et al. (2021), in the TSUMAPS-NEAM project documentation (Basili et al., 2019), and in several studies dealing with specific aspects of the SPTHA methodology used therein (e.g., Molinari et al., 2016;Selva et al., 2016;Davies et al., 2017;Glimsdal et al., 2019;Scala et al., 2020;Taroni and Selva, 2020;Tonini et al., 2020).

Design of Italian Inundation Maps
Inundation maps in Italy feature two zones, corresponding to the advisory/orange and watch/red alert levels, respectively (Dipartimento della Protezione Civile, 2018). All the inundation maps are available through the Tsunami Map Viewer (http://sgi2.isprambiente.it/tsunamimap).
The extent of the advisory/orange zone is nominally set to 1 m elevation everywhere, to be consistent with the definition of advisory alert level (runup ≤1 m). Currently, to conservatively accommodate uncertainties of the DTM, the boundary of the advisory inundation zone is fixed where the DTM features a topography of 2 m above the average sea level.
The watch/red zone should define the inland extension of the potentially inundated area after an alert message forecasting a runup >1 m. To define an upper limit for the inundation distance the following procedure is adopted: 1. From NEAMTHM18, the MIH with the 2,500 years ARP at the 84th percentile of the epistemic uncertainty distribution is selected for each POI along the Italian coasts (among those in Figure 1A).   Figure 2. There are several exceptions to this algorithm, described in the Annexes to the "Indicazioni" (Dipartimento della Protezione Civile, 2018), due to specific local geomorphological features: the very shallow waters in the northern Adriatic Sea where no POI lies within 40 km from the coast; the need to avoid associating POIs offshore the opposite coasts of the narrow peninsula of southern Calabria; the need to associate more than one POI to small islands, for which a search radius of 80 km is chosen, to better characterize SPTHA variability. 3. The design MIH d value is then converted into a design maximum run-up value R p d 3pMIH d . This factor of three is based on the analysis of the ratio between the maximum simulated runup and the mean of the simulated MIH at six different coastal locations in the Mediterranean Sea, including a variety of coastal morphologies (Glimsdal et al., 2019). The inundation simulations were carried out taking into account tsunamis generated at different distances by earthquakes having different fault orientations and various magnitudes. Figure 3 shows that the 98% of such ratios are smaller than 3. 4. The design runup R d defines the design inundation distance with the function D d (R d ) as the distance from the coast where the topography height over the sea level is R d ( Figure 4A). An exception is made when the coast is relatively flat (e.g., a coastal plain), which would result in a too large inundation given the dissipative character of inland propagation. An empirical dissipation factor is introduced to limit the maximum inundation distance (e.g., Leonard et al., 2008;Løvholt et al., 2012;Davies et al., 2017). The dissipation factor is derived from observed data collected during post tsunami field surveys (Fraser and Power, 2013, and references therein): an empirical rule featuring a 1 m reduction of the R d every 200 m away from the coastline (1 m every 400 m along and surrounding a river mouth) is suggested. This rule is applied to the expected resulting maximum inundation distance D d from the coast, which cannot exceed D MAX d 200pR d ( Figure 4B). D MAX d is obtained using a GIS-based approach to convert R d into inundation lines over a DTM, through following a landward oblique line tilted down of an angle α defined by the attenuation coefficient as tg(α) 1/200 or 1/400 ( Figure 4). Summarizing, the .

The local inundation zone is obtained for
One meter is added for taking into account the uncertainty on the DTMs and the variability in the accuracy of the coast line mapping, both of which may be difficult to quantify. 6. To find the inundation zone corresponding to R p d , a discrete set of possible inundation zones along the Italian coasts is first defined according to the classes D p that is by constructing with the GIS-based procedure described above the inundation lines corresponding to the set of discrete input run-up values R d {2, 5, 10, 15, 20, 25} expressed in meters. Then, the R p d is rounded to the closest R i ≥ R p d . This defines the inundation zone. This first generation of the Italian tsunami inundation maps for emergency planning was elaborated using a composite DTM combining elevation data from different sources. High-resolution LIDAR data (1 m and 2 m resolution with a vertical accuracy as good as 25 cm), provided by MATTM (Ministry of the Environment) for the coastal zone only, have been integrated with the best available mid-resolution DTMs (2, 5, 10 m), provided by single administrative Regions of Italy, in order to FIGURE 3 | Empirical probability density function (EPDF, blue bars) and empirical cumulative density function (ECDF, orange bars) of the ratio between the maximum Run-up (R d ) and the mean of the maximum inundation height (MIH ) at each POI.
FIGURE 4 | Sketch of the method used to drawn inundation maps: R p d is the maximum estimated run-up value for a given stretch of coast, D p d is the inundation distance and, in the example, the dissipation factor is tgα 0.005 (no rivers). The cartoon illustrates the two different behaviours when R p d crosses the topography before than the projected ray on the sea level (panel A) or vice versa (panel B).
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 628061 6 obtain the complete coverage of the potential inundation areas. Other input data are the trace of the main waterways (scale 1: 10.000) as an update of the national river networks of ISPRA-SINAnet (http://www.sinanet.isprambiente.it/it/sia-ispra/ download-mais/reticolo-idrografico/view), and the coast line derived from satellite and aerial images analysis available for the year 2009 at the Geoportale Nazionale (http://www.pcn. minambiente.it/mattm/), used only when applying the GIS-based procedure to achieve greater precision (item 6. above), while it was not deemed necessary to replace the less precise ISTAT coastline of (item 2.), since the latter is used in conjunction with search radii of 40 km. For both the segments of coast here considered as case-studies in south-eastern Sicily we used the hazard curves of the two nearest POIs to the locations ( Figure 1) and the run-up design value approximated for excess to the next discrete values is R p d 10 m.  To test the performance of the method described in the previous section, in this study we produce several new inundation map versions.
A first version, coined GIS-based version, is obtained by repeating the procedure described in Design of Italian Inundation Maps section, but adopting a more accurate local DTM. More specifically, we used the same high-resolution DTM (10 m) implemented in the tsunami simulations to compute local SPTHA (see next section and also Gibbons et al., 2020). This DTM was produced by mosaicking rasters with different horizontal cell resolutions and then resampling them to a 10 m horizontal cell. For the land portion the used datasets were a 1 m cell resolution raster near the coastal zone (The Geology and Geotechnologies Laboratory, INGV, http://www.ingv.it/it/ monitoraggio-e-infrastrutture-per-la-ricerca/laboratori/ laboratorio-geologia-e-geotecnologia, 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) with the 25 m horizontal cell resolution EU-DEM Digital Elevation Model over Europe (EU-DEM-4258: 1 arcsec -5 arcsec, EU-DEM-3035: 25 m, color shaded DEM derived from the EU-DEM-3035: 25m, produced using Copernicus data and information funded by the European Union -EU-DEM layers); for the offshore portion the DTM was obtained by interpolating the data available from the MaGIC project (foglio 32 e foglio 33, data from the MaGIC project, Dipartimento della Protezione Civile, http://dati.protezionecivile.it/) and then mosaicking the result to the EMODnet Digital Bathymetry raster (EMODnet Bathymetry Consortium, 2018). For both case-studies, we used the R p d values of 2, 5, 10 m, as defined for the current national administrative classification. The comparison between the inundation zones obtained with this DTM and the official ones shows some discrepancies resulting from the different characteristics of the DTM (resolution; ground filtering; different datasets used in the mosaic; different sources and data processing techniques for building the DTM from raw data) and from the different number and extension of the rivers considered in the analysis. While differences in the DTMs are responsible for relatively minor changes in the inundation maps, neglecting or including one or more rivers can impact significantly on the shape and the extent of the inundation line, since the assumed dissipation rate is halved in the presence of a river ( Figure 5 and Figure 6).
A second version of inundation maps is estimated adopting the same procedure, but using specific values of R p d , derived in this study from the NEAMTHM18 hazard curves for the two closest target POIs (Figure 1), and without approximating for excess to the next bigger discrete value. The obtained values are R p d,POI 3.6 m and R p d,POI 6.9 m, respectively for Catania and Siracusa (Figure 7). We anticipate that a third version is obtained for the values of R p d resulting by considering a rough quantification of the tsunami intensity uncertainty superposed onto the Non Linear Shallow Water (NLSW) results, for which a hybrid version of the inundation maps will be coined. We obtain R p

Inundation Maps From Local High-Resolution NLSW-Based SPTHA for Catania and Siracusa
Another version of inundation maps, coined SPTHA-based inundation maps, is obtained using a classical approach (e.g., González et al., 2009) directly from the local SPTHA, based on explicit, high-resolution (NLSW) modeling of the tsunami inundation (Gibbons et al., 2020).
The use of a regional SPTHA as input for a higher-resolution local SPTHA allows for a pre-screening of the full variability of seismic sources, selecting only those which contribute significantly to the hazard at the site of interest. Reducing the number, or the computational cost, of high-resolution tsunami Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 628061 inundation simulations is often essential. In this respect, techniques based on cluster analysis Volpe et al., 2019;Williamson et al., 2020) or statistical emulators (Sarri et al., 2012) have been proposed. In the present work, given the availability of massive High-Performance Computing resources, a significant reduction of simulations is obtained with a hazard disaggregation approach (Gibbons et al., 2020). The number of scenarios defined in NEAMTHM18 for the Mediterranean basin is of the order of 10 7 , making the computational task potentially extremely demanding. However, the tsunamis from many of these scenarios are negligible at the sites of interest. A first significant reduction is obtained by selecting only the scenarios which produce a tsunami elevation greater than 0.01 m at the selected target POIs (red dots in Figure 8), resulting in little more than 10 6 scenarios, most of which in common between Catania and Siracusa, see Table 1. Then, hazard disaggregation is performed to select only the scenarios which reproduce a prescribed fraction of the total hazard (99%) simultaneously at the two POIs in the MIH interval between 1 m and 4 m (Figure 9). The reason for choosing this interval is that it should guarantee the inclusion of the scenarios which significantly contribute to the hazard curves around the design 2,500 years ARP and 84th percentile of the model uncertainty selected for planning, (Figures 1B,C). This is exactly the same approach followed by Gibbons et al. (2020), and most of the numerical simulations are the same as well. As already discussed by Gibbons et al. (2020), given that the selection and the hazard reproduction are operated at two offshore locations, we cannot Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 628061 expect to perfectly reproduce the hazard everywhere onshore, because the difference between two tsunami waveforms may change during inundation: for example, this can be due to nonlinear effects such as wave breaking, energy dissipation during inundation, or because of co-seismic coastal uplift or subsidence. It is also worth noting that the regional hazard curve for Siracusa is significantly higher than that for Catania for the considered ARP ( Figure 1). For this reason, while the MIH range between 1 and 4 m turns out to be suitable to produce the onshore hazard curves for Catania (Gibbons et al., 2020), to better reproduce the hazard for Siracusa, we also modeled the scenarios obtained from the disaggregation in the interval 4-8 m, in addition to the 1-4 m range. We may conclude that the disaggregation procedure and applied ranges turned out to be quite accurate while remaining remarkably efficient, since the scenarios to be modeled were reduced, in total, to about 46,000 (i.e., the amount of scenarios is reduced by ∼96%, see Table 1).
For the sake of completeness, we recall that the tsunami scenarios are simulated with Tsunami-HySEA, a non-linear hydrostatic shallow-water multi-GPU code (GPU, Graphics Processing Unit; de la Asunción et al., 2013;Macías et al., 2017), exploiting a nested grid algorithm. Tsunami-HySEA models in a single code tsunami generation, open-ocean propagation and nested grid inundation using progressively finer grid resolution of the coastal areas (Macías et al., 2016). 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, United States (Macías et al., 2017;Macías et al., 2020a;Macías et al., 2020b). Four-level nested topobathymetric grids were used here (Figure 8). Three of them are local grids with progressively higher resolution (160, 40, and 10 m) and one is a global grid with the coarsest resolution (640 m), ensuring the analysis of the source-to-target propagation in the open sea. The 640 m grid was derived from the GEBCO topo-bathymetry model. The finest local 10 m grid is the same presented in the previous section for the construction of the inundation maps. To guarantee depth compatibility between all the nested local grids, a bilinear resampling algorithm was then applied to the 10 m grid, producing the local intermediate grids with 40 m and 160 m cell resolution (Gibbons et al., 2020). Bottom friction is treated with a Manning's coefficient of 0.03. The initial sea-level anomaly is obtained for each modeled seismic scenario through an elastic instantaneous dislocation model (Okada, 1992). The single-scenario Tsunami-HySEA simulation takes approximately 25 min on the MARCONI-100 supercomputer at CINECA, meaning that approximately a Total is calculated excluding scenarios whose simulation should be repeated due to the vicinity of the two target areas.
FIGURE 9 | Number of scenarios needed to represent the 99% ( 19,000 GPU hours were required to complete the set of selected scenarios. With 1024 simulations able to run in parallel, around 18 h of clock-time are needed to perform these calculations. The simulation results obtained in the highest-resolution innermost grids at the two locations are then aggregated to calculate the hazard curves in each grid cell on land if reached by a minimum threshold of water (here, 0.01 m). Specifically, the hazard curves are obtained in each point x as where dT is the exposure time (50 years as in NEAMTHM18), h is the threshold inundation height, s i are individual scenarios with mean annual rate λ(s i ), H I(si,x) > h is an indicator function (sometimes called Heaviside step function) that is 1 if the intensity I(s i , x) computed by NLSW simulations for the scenario s i in x is > h, and 0 otherwise, and λ h is the overall mean annual rate of exceedances of h. P(I > h; x, 50yr) represents the Probability of Exceedance (hereinafter PoE), that is the probability that h has been exceeded at least once during the time dT 50 years. The use of indicator functions implies that no uncertainty is modeled into the propagation factor of the hazard, that includes the tsunami generation, propagation and inundation process (Selva et al., 2016;Grezio et al., 2017). The same seismic model [s i , λ(s i )] used to create NEAMTHM18, which is at the base of the inundation maps described in Design of Italian Inundation Maps and Inundation Maps for the Case-Studies in South-Eastern Sicily sections, is also considered for calculating this local SPTHA. This source model includes epistemic uncertainty (Basili et al., 2021), represented through an ensemble of 1000 equally plausible models with alternative definitions of the scenarios and of their rates of occurrence. Collectively, this ensemble represents the seismic source variability as a function of many alternative but scientifically acceptable modeling choices (Selva et al., 2016;Basili et al., 2021). This epistemic variability propagates to hazard curves. Ideally, it should be combined with the epistemic uncertainty in the propagation factor (not modeled here). The resulting ensemble of hazard curves describes the epistemic uncertainty on the hazard model through different statistics such as mean or percentiles.
Inundation maps are then derived by selecting the tsunami intensity corresponding to a given value of the ARP, which is the reciprocal of the mean annual rate λ h in Eq. 1 and thus can be evaluated as a function of the PoE. Graphically, this means to intersect the hazard curves, at each grid point, with an horizontal plane having a constant PoE, corresponding to the selected ARP, and to draw a contour line corresponding to the wet area (maximum flow depth greater than 0.01 m, since this is the first threshold for which we calculated SPTHA; the tolerance for simulations with Tsunami-HySEA was instead set at 0.001 m). If needed, it is possible to draw different inundation maps for different ARPs, in the time span defined in the frame of the NEAMTHM18 model (up to about 100,000 years). Inundation maps can be obtained either from the mean hazard curve, or from any other hazard curve calculated for a different percentile.
Finally, it is important to remark that some uncertainties associated with the tsunami source (e.g., heterogeneous material properties, variability of fault dimensions, timedependence, non-hydrostaticity), or with the tsunami propagation and inundation modeling (e.g., limits of the NLSW model, uniform Manning coefficient, wave dispersion, details of coastal morphology) are not quantitatively addressed in the present study. Neither is addressed the impact associated with potentially inaccurate bathymetry and topography models.

RESULTS
The main goal of our analysis was to compare the current strategy to define inundation maps, used to draw the evacuation areas for the Italian coasts and based on the regional hazard NEAMTHM18, with the hazard results obtained from local SPTHA models, hence evaluating the whole chain of adopted assumptions and approximations described in Design of Italian Inundation Maps section. This analysis was carried out in the two target areas of Catania and Siracusa.
We recall that the GIS-based inundation maps and those from the local SPTHA based on numerical simulations, which are compared to each other in this section, were obtained using the same DTM. In Inundation Maps for the Case-Studies in South-Eastern Sicily section, we have also qualitatively addressed the impact of using a different DTM and different treatments for the rivers onto the design of the inundation maps; we have observed that the latter is more influential than the differences between the DTMs (Figures 5, 6).
As described in Design of Italian Inundation Maps section, for national planning the R p d design value was selected by searching the maximum tsunami intensity in an area containing several POIs (at distances up to 40 km), leading to a value greater than 5 m for both locations (Figure 2). Given the rough discretization of the inundation lines in the GIS-based method, such a value is eventually rounded up to the 10 m and the corresponding contour line is used. On the one hand, it could be argued that the value of R p d , corresponding to the specific POI located in front of each site of interest, could be a more calibrated solution to estimate the local inundation. On the other hand, the use of a regional assessment can be inappropriate for specific land use actions, since the tsunami behaviour in the near-field depends strongly on non-linear effects and detailed terrain morphology not explicitly treated at regional level. Moreover, the local sources may not be sufficiently finely discretized nor investigated closely enough. As a matter of fact, the closest POIs in front of Catania and Siracusa ( Figure 8B) provide the two values R p d,POI (as defined in Inundation Maps for the Case-Studies in South-Eastern Sicily section) equal to 3.6 m and 6.9 m, respectively. The first one is smaller than the value used for planning, which confirms that local effects are relevant, and not homogeneous from site to site. The inundation lines which correspond to these local values are shown in Figure 7 for both Catania and Siracusa, and compared to the line corresponding to R p d 10 m used for planning. A first and direct comparison between the inundation maps from the local SPTHA at 2,500 years ARP for the 84th percentile (Inundation Maps From Local High-Resolution NLSW-Based SPTHA for Catania and Siracusa section) and the GIS-based inundation maps (Inundation Maps for the Case-Studies in South-Eastern Sicily section; the first version, with the same design values as for the national planning) has been performed ( Figure 10). The choice of R p d 10 m for the national planning appears largely conservative, since the extension of the watch/red area exceeds greatly the inundation line calculated with local SPTHA for the same ARP and percentile. On a second glance, by comparing the inundation maps in Figure 10 to those of Figure 7, it appears immediately clear that the local SPTHA predicts an inundation even much more limited than that associated with the POIs at both locations.
Selecting different ARPs and/or different percentiles, we can quantify the extent of this discrepancy. Firstly, we compare the impact of the epistemic uncertainty modeled on the local SPTHA by using different percentiles while keeping the ARP of 2,500 years fixed (Figures 11A,B). As expected, the inundation zone for the 98th percentile extends farther than those for the 50th and the 84th percentiles. However, they remain well within the watch/red inundation area defined with R p d 10 m for both Catania and Siracusa. This can either mean that the discrepancy is greater than it could be explained by the epistemic uncertainty on the 2,500 years ARP, hence the GIS-based method leads to a large overestimation, or that there is a fraction of unmodelled epistemic uncertainty. Secondly, when longer ARPs are considered, the tsunamis from larger but rarer earthquakes progressively gain more importance. Figures 11C-F displays the inundation maps produced for ARPs 25,000 and 100,000 years, considering the 50th, 84th, and 98th percentiles of the epistemic uncertainty. In this case, the watch/red inundation zone still contains the inundation lines in Siracusa, but it is overtaken at some locations in the Catania plain. This could be due to the very different morphology of the coastal area, where the large plain in the south of the harbour could facilitate the water invasion beyond the forecast capability of the empirical rules used for the GIS-based inundation maps. As a result of this comparison, according to the local SPTHA inundation maps, in Catania, the GIS-based inundation maps would correspond to a design ARP between 2,500 and 100,000 years, while in Siracusa to a design ARP >100,000 years.
It is worthy to note that the effect of the epistemic uncertainty for the 100,000 years ARP is smaller than for shorter ARPs (Figure 11). The hazard curves for a single target point (Figure 1) show, as expected when dealing with relatively rare events, that the longer the ARP (smaller PoE), the larger the difference between the tsunami intensity obtained with different percentiles of the epistemic uncertainty. However, when we look at the inundation areas, we observe the opposite behaviour, with the uncertainty on the hazard intensity unexpectedly decreasing for longer ARPs. We then argue that this can be the result of an underestimation of the epistemic uncertainty of the local SPTHA, and that the discrepancy between GIS-based and SPTHA-based maps cannot be entirely due to an overestimation of the former. Indeed, we stress once again that the regional SPTHA includes the epistemic uncertainty associated with the source model (Basili et al., 2021), which is the same as in the GISbased approach. This uncertainty relates to the earthquake mean annual rates, the difference between empirical scaling relations, the slip heterogeneity and depth-dependent features of subduction zones (Scala et al., 2020). Conversely, the local SPTHA does not include the further epistemic uncertainty treatment applied in NEAMTHM18 using amplification factors (Glimsdal et al., 2019), stemming from numerical modeling of the co-seismic displacement, tsunami generation, and tsunami evolution up to inundation. In our local SPTHA, this uncertainty is in principle reduced by the use of NLSW simulations in place of the stochastic inundation treatment included by the local amplification factors, provided that we are not introducing a bias. Several sources of epistemic uncertainty currently exist also when using numerical inundation simulations, Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 628061 as discussed at the end of Inundation Maps From Local High-Resolution NLSW-Based SPTHA for Catania and Siracusa section. Nevertheless, the development of a robust inclusion of epistemic uncertainty into NLSW-based SPTHA is beyond the scope of this paper. As an alternative, one could also forbear from making detailed simulations, and perform a finer-scale version of the NEAMTHM18 PTHA, with more closely-spaced POIs, to better represent the MIH uncertainty through the higher-resolution amplification factors.
Here, as anticipated at the end of Inundation Maps for the Case-Studies in South-Eastern Sicily section, and similarly to Scala et al.
(2020), we roughly and probably much more conservatively quantify these uncertainties following Davies et al. (2017), such as: where Φ is the cumulative distribution function of a log-normal distribution with mean [I(s i ) + β] and standard deviation σ, I(s i ) is the tsunami intensity modeled by the NLSW simulation for the scenario s i , and β is a logarithmic bias. Based on the tables available in Davies et al. (2017) for four large tsunamis (1960Chile, 1964Alaska, 2004Sumatra, and 2011 Tohoku earthquakes), we chose conservative values for parameters β 0.2 and σ 1. This uncertainty has been originally set by comparing modeled intensity and run-up along large portions of coasts, thus it is possibly not ideal for the application to our local case-study; it may also include the contribution from landslide sources for the Alaska 1964 event. We speculate that, if locally calibrated, at least via numerical experiments considering uncertainty in the bathymetry and topography, and with variable bottom friction, the uncertainty on local NLSW simulations is expected to be smaller than that implied by Eq. 2. That being said, and in the absence of any other more specific method, we use this approach in the attempt to quantitatively set an upper bound to the model's uncertainty, hence to the extension of the inundation zone. To apply the uncertainty model, we replace the identity functions of Eq. 1 with the expression in Eq. 2. We note that this correction does not affect the areas where each tsunami scenario did not propagate, being the log-normal distribution applicable only to strictly positive values, but it corrects the MIH estimations in the flooded area, for each scenario (Figure 12). In other words, adding the log-normal correction to the deterministic simulations cannot extend the inundated area. Thus, to evaluate the consequent inundation zones, we apply an approximated procedure mimicking the one described in Design of Italian Inundation Maps and Inundation Maps for the Case-Studies in South-Eastern Sicily sections for the determination of the maximum inundation distance, whose input is the design run-up R p d , corresponding to the maximum expected runup in the area. We estimate the maximum calculated inundation heights for the 2,500 years ARP and 84th percentile over all the simulation domain, finding R p d,LOC

m and R p d,LOC
15.1 m in Catania and Siracusa, respectively. Such values occur relatively close to the coast, before coastal dissipation may become important, reflecting the fact that the maximum inundation height above the sea level is greater or equal than the maximum run-up. These numbers are much larger than the ones obtained from the regional SPTHA, that is 3.6 and 6.9 m, respectively, confirming that the overall effect of the uncertainty model is to increase the hazard, as well known in seismic hazard (e.g. Bommer and Abrahamson, 2006), and as already noted by Gibbons et al. (2020). The last version of the inundation maps is obtained from these values, using the simplified attenuation method described in Design of Italian Inundation Maps and Inundation Maps for the Case-Studies in South-Eastern Sicily sections (Figure 13). In a sense, this is a hybrid version, using both the results of the local SPTHA, with the addition of the log-normal uncertainty on each scenario, and the GIS-based attenuation method for the determination of the maximum inundation distance. We observe that the GIS-based inundation line used for planning (red zone in Figure 13) is now generally exceeded by the one we obtain with the newly proposed approach (orange zone in Figure 13). We argue that the hybrid version provides an upper bound to the uncertainty on the inundation distance, due the conservative parameters we have chosen here for the log-normal. The lower bound for the uncertainty is constituted by the maps obtained with the deterministic NLSW inundation scenarios.
We performed one last experiment, to check how inundation maps based on a worst-case scenario approach would compare with the GIS-based and the hybrid ones we have just introduced by adding uncertainty to the local SPTHA-based maps. We define the worst-case inundation scenario as a composite one from the envelope of all the wet grid points (with wave height >0.01 m), from all the scenarios modeled for the SPTHA, irrespective of the single scenario probability. Even though we have not traced back the scenarios contributing to the maximum inundation distance, we may argue that they are a limited number of relatively high magnitude earthquakes. Actually, such a worst-case corresponds to the SPTHA inundation map for an ARP tending to infinity, which includes all points where a non-identically-zero hazard curve exists (that is exceeding with non-zero probability the threshold of 0.01 m). This composite inundation distance from the envelope of all maximum inundation distances is shown in Figure 14, along with the other inundation maps. In the Catania plain, the worst-case inundation distance exceeds all the others with exception of limited zones along the rivers. In Siracusa, the comparison is more complex: in the northern area the worst-case inundation is similar to the extension of the R p d 10 m inundation map (in red), whereas, in the middle and southern areas, its extension reaches and, occasionally, exceeds the inundation map corresponding to R p

DISCUSSION
The method for designing tsunami evacuation maps for the Italian coasts has been analyzed to understand its limits and advantages, and for addressing to what extent it can be updated through local SPTHA studies. Our findings highlight several important points that should be taken into account for both local high-resolution SPTHA and tsunami coastal and evacuation planning.
The GIS-based method is severely affected by how many secondary riverbeds (depending on their sizes or flow rates) are assumed to be relevant, also when using the same DTM ( Figures 5, 6). The impact of less accurate DTM models had, at least in our case-study, a relatively smaller impact on the definition of inundation zones. However, the situation can be different in other cases depending on the extent of the discrepancies between the DTMs, so we cannot generalize this conclusion. It is noted that DTM uncertainty can be significant,  Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 628061 15 and the treatment can be imposed based on the setting. In PTHA for landslides in Norway for instance, evacuation zones occur on very steep slopes, and here 10 m horizontal distance is added as an additional buffer to the computed hazard zones (e.g., Løvholt et al., 2020). This equals the resolution of the local DTM, and due to the steep slopes, this can involve additional safety factors of the order of 1 m or more.
Regional-SPTHA-informed, GIS-based, inundation maps (defined in Inundation Maps for the Case-Studies in South-Eastern Sicily section) were compared with those derived from a local SPTHA using high resolution inundation simulation, in two specific target areas. Two SPTHA models have been used: one without considering epistemic uncertainty specifically associated with numerical inundation simulations (Inundation Maps From Local High-Resolution NLSW-Based SPTHA for Catania and Siracusa section), and one, dubbed as "hybrid", which includes a possible maximum estimate of this uncertainty (Results section). In a sense, we have defined an "error bar" for the inundation distance related to a given ARP and percentile from local SPTHA. In this way, we are preliminarily exploring the combination of deterministic tsunami modeling with the epistemic uncertainty in tsunami generation, propagation and inundation modeling. We recall that the superposition of uncertainty with deterministic simulations was already considered by some authors, yet in the context of a worst-case scenario approach (Tonini et al., 2011).
Our most important outcome is possibly that the inundation maps resulting from the GIS-based method used for national planning lie within the error bar bound by the NLWS-based model with and without uncertainty treatment. The GIS-based maps overestimated the extension of the inundation zone with respect to the maps obtained from local SPTHA with no uncertainty in the tsunami propagation model. Conversely, they underestimated the extension of the inundation zone when an upper bound of this uncertainty is considered in the hybrid approach.
Another significant outcome is that it results crucial to further characterize these uncertainties to reduce as much as possible this "error bar" or, in more classical terms, to reduce the epistemic uncertainty. Noteworthily, we find that ignoring the epistemic uncertainty of tsunami modeling may have a major impact on local hazard results, potentially producing severe underestimations. Thus, even if local SPTHA studies, in principle, reduce the uncertainty on the definition of local run-up maxima, disregarding completely the related uncertainty may be dangerous. We remark that the parametrization we used is a maximization of the epistemic uncertainty, derived from a global analysis. This highlights two main findings: 1) the important role played by the (usually neglected) uncertainty in tsunami modeling in hazard assessment and, consequently, in the definition of evacuation planning strategies and 2) the need of detailed local analyses based on local data to better calibrate the parametrization of the uncertainties. Nevertheless, it seems that in the absence of more detailed studies, the currently adopted choices for the national planning, which were the outcome of extensive discussions among the decision-makers and different scientific actors at the national level, remain fully justified.
Further considerations on such choices can be certainly made. The use of a regional SPTHA, as main input to determine the value of R p d , may be not ideal for local hazard quantifications. Indeed, results of the regional hazard are projected from the 50 m isobath to inland with mainly 1D assumptions and models; thus, any lateral effect (i.e., parallel to the coast) is difficult to be captured. This source of uncertainty was accounted for by using a multiplicative factor set equal to 3, as explained in Design of Italian Inundation Maps section. This choice maximizes the level of safety since it roughly represents the 98th percentile of the uncertainty estimated for the amplification factor, with the objective of dealing with local "extreme" run-up values, that is with locally maximum inundation height and maximum runup variability around the average inundation height. The maximum inundation distance depends also on the local topography and the hydrodynamics, which are not always simply related to expected/design height. This is for example emphasized by the effect of the possible coseismic displacement altering the coastal elevation for near-source coasts, which is not modeled in the regional SPTHA (Volpe et al., 2019). A predominant uplift was in fact observed at least for the Catania use-case (Gibbons et al., 2020), which to some extent limits the inundation probability. Conversely, in case of subsidence in the target area, inundation would be facilitated, worsening the impact inland. Moreover, different slope steepness, for example ranging from cliffs to plains, or different terrain types or different land use in the coastal areas may determine the extent of nonlinearity and dissipation during the inundation. Here, we have only considered the main river streams to differentiate the attenuation rate of tsunami waves inland, and overlooked minor rivers and waterways. The two considered target areas are very different from one another: the target area in the south of Catania is a flat, straight and long plain mainly characterized by beaches, uncovered terrain and few buildings, which is very different from the steep and rocky coast featured with dense urban infrastructures of Siracusa. Better results could be derived by introducing different coastal dissipation factors as a function of the land use and by acquiring high resolution LIDAR data that completely cover the whole potential inundation area. However, the adopted empirical dissipation rule, in combination with the other choices, seems to be quite robust, according to the comparisons we have presented.
A final remark would be that, as a further experiment, a worst-case inundation scenario was produced and compared to the other approaches. This worst-case is based on maximizing the estimated impact, rather than using the biggest earthquake magnitude available. This is similar to the definition of the 2,500-years Tsunami Design Zone in the United States (Wei et al., 2017), based on the PTHA offshore tsunami amplitude at 100 m depth, in which a "superscenario" is defined according to the maxima of the maxima of all the simulated tsunamis in front of the coast and, then, the corresponding inundation is numerically calculated. As a result of the huge range of scenarios included in our seismic model, the variability for this SPTHA-based worst-case is more extensively explored than usually done for typical worst-case studies, based on a few earthquakes only. We argue that using such an approach may make sense for some critical infrastructures, yet complementing the worst-case effects by using some information from SPTHA about their ARPs. However, as for local-SPTHA, we note that the uncertainty associated with NLSW simulations is not considered also in the definition of this worst-case. On the other hand, if inundation uncertainty was introduced, the definition of worst-case would become even more difficult (e.g., log-normal are unbounded for large values), which challenges by itself the meaning of a worst-case approach not characterized by probability or ARP ranges. Without this, the meaning of the differences with an SPTHAbased approach may remain unclear.

CONCLUSION AND FUTURE DIRECTIONS
Tsunami evacuation maps designed for the Italian coasts have been evaluated and discussed using the results of a local SPTHA performed in the two target areas of Catania and Siracusa, Sicily, for which it was possible to run high resolution tsunami inundation simulations.
We demonstrated that, in lack of a local and detailed hazard analysis, the proposed procedure to define evacuation areas has a good degree of reliability, since the resulting maps are bound by the "error bar" defined by the outputs obtained from a local SPTHA performed with and without the estimation of its epistemic uncertainty.
The analysis underlines the importance to perform local studies using the proper level of data resolution, whenever available, to best serve specific risk mitigation actions, which may require a thorough site-specific uncertainty quantification. Nevertheless, strategies to improve the input regional assessment, such as using a higher density of POIs and/or a refinement of source parameters, would not only improve the regional analysis itself but also enhance the subsequent local analyses.
Awaiting for future site-specific PTHA studies, our results highlight the importance of quantifying the uncertainties that can affect numerical inundation modeling and are instead usually neglected. Not modeling these uncertainties can introduce important biases in the hazard evaluations. However, we have only scratched the surface here, and a specific definition of an uncertainty framework for tsunami inundation models and the calibration of site-specific SPTHA with local data is still required. A detailed local analysis, with massive simulations to explore the impact of the different assumptions and uncertainties is needed. The ongoing PRACE Project "Local Probabilistic Tsunami Hazard Assessment for HPC -TsuHazAP" (see https://prace-ri.eu/hpcaccess/project-access/project-access-awarded-projects/projectsawarded-under-prace-project-access-call-21/) is providing us a framework to start dealing with the above identified issues, and in particular with the sensitivity to local seismic source mechanism discretization, the effects of uncertainty on bathymetry and topography, and treatment of onshore friction.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
RT, PM, SL, JS, and MV conceived the experiment and prepared the first draft of the manuscript. All the authors contributed to the implementation of the analysis, to revising the article and approved the submitted version.