Tsunami coastal hazard along the US East Coast from coseismic sources in the Açores convergence zone and the Caribbean arc areas

Tsunami coastal hazard is modeled along the US East Coast (USEC), at a coarse regional (450 m) resolution, from coseismic sources located in the Açores Convergence Zone (ACZ) and the Puerto Rico Trench (PRT)/Caribbean Arc areas. While earlier work only considered probable maximum tsunamis, here we parameterize and simulate 18 coseismic sources, with magnitude M8-9 and return periods ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}70–2000 year, using seismo-tectonic and historical data. The largest sources in the ACZ are repeats of the 1755 M8.6-9 Lisbon earthquake and tsunami; other sources are hypothetical. In the ACZ, due to the limited data on faults, each source is parameterized with a single fault plane, while in the PRT, coseismic sources are parameterized based on fault segmentation established during a 2019 USGS workshop of experts, using 10–26 SIFT subfault planes (Gica et al. in NOAA Tech. Memo., OAR PMEL-139, 2008). Tsunamis are simulated for each source using the fully nonlinear and dispersive model FUNWAVE-TVD, in two levels of nested grids. At the considered scales, dispersion is shown to affect tsunami propagation. Coastal hazard is quantified by four metrics computed at many save points (∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}20–30 thousand) defined along the 5-m isobath (due to the coarse resolution), i.e., maximum (1) surface elevation, (2) current, (3) momentum force; and (4) travel time, representing flooding, navigation, structural, and evacuation hazards. Overall, the first three metrics are larger, the larger the source magnitude, and their alongshore variation shows similar patterns of higher/lower values, due to the shelf bathymetric control (refraction). The fourth metric mostly differs between sources from each area, but less so among sources from the same area; its inverse quantifies evacuation hazard. A 1–5 score is given to results for each metric, based on five intensity classes representing low, medium low, medium, high, and highest tsunami hazard. A novel tsunami intensity index is computed as a weighted average of these scores, allowing both a comparison among sources and a quantification of tsunami hazard as a function of their estimated return periods. In the most impacted areas of the USEC, the highest tsunami hazard in the 250–500-year return period range is commensurate with that posed by 100-year category 3–5 tropical cyclones, taking into account the larger current velocities and forces caused by tsunami waves. Results of this work could serve as a basis for a future regional Probabilistic Tsunami Hazard Analysis for the USEC, considering additional source types such as underwater landslides, volcanic flank collapse, and meteotsunamis, that were studied elsewhere.


Introduction
Since 2010, under the auspice of the US National Tsunami Hazard Mitigation Program (NTHMP; http:// nthmp. tsuna mi. gov/ index. html), the authors and colleagues have performed tsunami modeling work to develop high-resolution (10 m) tsunami inundation maps for the US East Coast (USEC), starting with the most critical or vulnerable areas, but with the goal to eventually cover the entire coast. These "first-generation inundation maps" were constructed as envelopes of maximum inundation caused by the most extreme nearand far-field tsunami sources in the Atlantic Ocean basin, a.k.a., Probable Maximum Tsunamis (PMTs), parameterized based on historical or hypothetical events (Schambach et al. 2020b). In this first-generation work, no return periods were considered for each PMT, and a probabilistic tsunami hazard analysis (PTHA) was left out for future work.
The present work is part of the initial preparatory steps necessary to perform a complete PTHA for the USEC and create the next generation of probabilistic NTHMP inundation maps. In PTHA, instead of only considering PMTs, the hazard maps combine results of simulations for a collection of tsunami sources of different magnitudes and return periods, for each type of sources. Here, we only consider far-field coseismic sources in the North Atlantic Ocean Basin (NAOB), which are the only seismic sources that can be sufficiently tsunamigenic to cause a significant tsunami hazard along the USEC, where near-field seismicity is quite moderate (ten Brink et al. 2008Brink et al. , 2014. Additionally, coseismic sources, as a group, are believed to have lower return periods than the SMF and volcanic collapse sources in the NAOB (ten Brink et al. 2008Brink et al. , 2014; meteotsunamis are more frequent but are typically associated with a lower level of hazard on the USEC (e.g., Geist et al. 2014). Hence, this work also provides, to our knowledge, the first comprehensive assessment of regional scale coseismic tsunami hazard for the entire USEC, for a range of return periods including some low enough to overlap with those considered in hazard analyses of coastal/ ocean structures and infrastructures, from other natural hazards such as tropical cyclones (e.g., 100 to 500 years).
More specifically, we consider and model coseismic tsunamis generated by a collection of sources located near and around ( Fig. 1): (i) the Açores Convergence Zone (ACZ; Fig. 1 and see Fig. 18 in ten Brink et al. (2014)), which includes the estimated location of the 1755 Lisbon (LSB) earthquake source; and (ii) the Puerto Rico Trench (PRT) and Caribbean arc area (Lander 1997;ten Brink et al. 2011;Harbitz et al. 2012;Grilli et al. 2016). Several large tsunamigenic earthquakes have occurred in these areas. In the ACZ/LSB area, which is located along the western segment of the Eurasia-Nubia Plate Boundary, between the Azores archipelago and the Strait of Gibraltar (Baptista 2020), some of these historical earthquakes have caused large transoceanic tsunamis, with the Lisbon 1755 earthquake, of estimated magnitude M8.6-9 (Muir-Wood and Mignan 2009; Barkan et al. 2009), triggering the largest known historical tsunami. In this catastrophic event, in the near-field, 5-15-m-high waves impacted the coast of Portugal near and around Lisbon, causing extensive destruction and tens of thousands of fatalities (in combination with ground shaking and fires). In the far-field, large waves reached the coasts of Morocco, England, Newfoundland, Brazil, and the Antilles. In particular, after a long transoceanic propagation in the NAOB, the tsunami still caused a few meters of inundation in the eastern Lesser Antilles (e.g., Baptista et al. 1998aBaptista et al. , b, 2003Baptista et al. , 2009Zahibo and Pelinovsky 2001;Barkan et al. 2009;ten Brink et al. 2014;Baptista 2020). Large coseismic tsunamis have also occurred in the PRT and Caribbean arc area (e.g., Lander 1997; ten Brink et al. 2008;Grilli et al. 2010;Mercado and McCann 1998), which is a highly seismic area running parallel to the north shore of Hispaniola, Puerto Rico, and the northeastern lesser Antilles (e.g., Calais et al. 1995;Zahibo and Pelinovsky 2001;ten Brink et al. 2008ten Brink et al. , 2011ten Brink et al. , 2014Grilli et al. 2010;Harbitz et al. 2012). (See Fig. 3 for the observed distribution of earthquakes, in magnitude and depth in this area.) The PRT is the only subduction zone in the NAOB, in which the North American Plate subducts under the Caribbean Plate, with a nearly E-W relative motion (i.e., largely left lateral strike slip; see large black arrows in Fig. 3) and only a small component of perpendicular convergence (3-6 mm/yr) (Calais et al. 1995;Benford et al. 2012). Despite this highly oblique convergence, in: (i) 1842, a M8 earthquake occurred  Table 4). Location are marked for the two areas of historical/hypothetical tsunami coseismic sources (red oval) considered, near the Açores Convergence Zone (ACZ), including Lisbon 1755 (LSB), and near and around the Puerto Rico Trench (PRT). The Madeira Torre Rise (MTR; see also Fig. 2) is the shallower ridge located on the north of the ACZ circled area and the Horseshoe Plain (HSP) is to the east of the ACZ and MTR. Yellow/red symbols within the regional grids mark locations of numerical wave gauge stations where time series of surface elevation are calculated in simulations for validating the one-way grid coupling (Table 2) in the western segments of the Septentrional fault (SF; see Fig. 3), which runs nearshore parallel to the north shore of Hispaniola, triggering a large tsunami that caused significant destruction in Haiti (Calais et al. 2010;Gailler et al. 2015;Grilli et al. 2016;(ii) 1787, a M8.1 occurred in the PRT, triggering a moderate tsunami; and (iii) 1918, a M7.3 earthquakes occurred in the Mona Passage, 15 km off the northwest coast or Puerto Rico (see Fig. 3), triggering a tsunami that caused 116 fatalities and up to 6-m runup (Mercado and McCann 1998;López-Venegas et al. 2015).
Since the Indian Ocean M9.2 earthquake and tsunami (e.g., Grilli et al. 2007;Ioualalen et al. 2007), we know that highly oblique subduction zones can generate devastating tsunamis if the earthquake rupture has a large thrust component. Hence, although this is still somewhat controversial, large tsunamigenic M8.7-9 earthquakes have been proposed and modeled in the PRT (Knight 2006;Grilli et al. 2010), showing that a devastating tsunami could occur in the near-field, which, in the far-field, particularly along the upper USEC, could still cause a 2-3-m inundation. A recent meeting of experts at the USGS Powell Center (May 2019) reached the conclusion that such large PRT events were possible; however, the upper bound magnitude would likely require that fault segments in parts of Hispaniola on the west and the Caribbean arc (up to Guadeloupe; see Fig. 3) on the east be also involved (this will be further detailed later).
In this work, based on seismo-tectonic and historical data, and on the conclusions from the Powell center meeting for the PRT area, we parameterize and site a collection of M8-9 hypothetical coseismic sources in the two selected areas; for each source, we then simulate tsunami generation and propagation to the USEC. In the ACZ/LSB area, we consider ten (10) sources, with two different locations and orientations (strike) selected to maximize tsunami hazard on the upper and lower USEC, respectively (Grilli and Grilli 2013a); four of these sources are the M9 PMT sources already considered in earlier NTHMP inundation mapping work (Grilli and Grilli 2013a;Schambach et al. 2020b). In the PRT area, following recommendations on fault plane segmentation made during the 2019 Powell Center meeting, and we consider eight (8) sources with, for comparison, one of these being the M9 PMT source already used and simulated in earlier NTHMP inundation mapping work (Grilli et al. 2010;Grilli and Grilli 2013b;Schambach et al. 2020b).
For these 18 coseismic sources, tsunami simulations are performed using the fully nonlinear and dispersive Boussinesq model FUNWAVE-TVD (Wei et al. 1995;Shi et al. 2012;Kirby et al. 2013), by one-way coupling in a series of nested spherical or Cartesian coordinate grids. FUNWAVE has been extensively used to simulate historical (e.g., Watts et al. 2003;Grilli et al. 2007Grilli et al. , 2013Grilli et al. , 2016Grilli et al. , 2019Ioualalen et al. 2007;Tappin et al. 2008Tappin et al. , 2014Kirby et al. 2013;Schambach et al. 2020aSchambach et al. , 2021 and hypothetical (e.g., Grilli et al. 2010Grilli et al. , 2017aAbadie et al. 2012;Tehranirad et al. 2015;Schnyder et al. 2016;Schambach et al. 2019) tsunami case studies. This model was also validated against a collection of tsunami benchmarks (Tehranirad et al. 2011;Horrillo et al. 2015). As we aim at assessing the regional coastal tsunami hazard caused by the selected sources along the entire USEC, here, we only use two levels of nested grids ( Fig. 1): (i) two large scale spherical grids over the North Atlantic Ocean (or half of it), with a 1 arc-min resolution; and (ii) three smaller regional shore-parallel Cartesian grids overlapping along the USEC, with a 450-m resolution. This moderate resolution is not fine enough to accurately resolve tsunami coastal hazard for complex coastal morphologies and/or highly developed areas. However, it should be adequate to assess the regional tsunami hazard intensity along the USEC, particularly for the long barrier-beaches and barrier-islands that skirt the USEC from Florida to Massachusetts. In fact, Ioualalen et al. (2007) used a 450-m nested grid with FUNWAVE to simulate the impact of the 2004 Indian Ocean tsunami along the fairly straight coasts of southwest Thailand, and showed a good agreement between the predicted and observed runups at 58 locations. In view of the moderate resolution used, however, to reduce coarse grid effects, in this work tsunami hazard will be assessed for the inundated areas based on a few tsunami metrics computed for each source, as a proxy, at a short distance away from shore along the 5-m isobath with respect to mean water level (MWL).
High-resolution inundation and runup simulations will be left out for future work. To do so, similar to earlier work done for the collection of PMTs simulated by Schambach et al. (2020b), additional levels of nested grids would be used (e.g., with a 112.5-m, 30-m and 10-m resolution; see, for instance, ), starting with coastal areas having the largest regional hazard or deemed to be critical. High-resolution inundation maps would then be modeled for each of the new sources, based on incident tsunami waves such as computed here in our finest nested grids. Examples of first-generation maps covering a small area (typically a US coastal county) at a 10-m resolution can be found on our NTHMP project webpage: https:// www1. udel. edu/ kirby/ nthmp. html.
To more easily identify coastal areas with high tsunami hazard, the coastal tsunami hazard simulated for each source will be quantified using four hazard metrics computed based on FUNWAVE results along the 5-m isobath: (1) maximum surface elevation above mean water level (MWL); (2) maximum current velocity magnitude; (3) maximum momentum force; and (4) tsunami travel time. Overall, the first three factors are larger, the larger the source magnitude, and the fourth factor is expected to differ mostly between sources from each area (i.e., ACZ and PRT), but less so among sources from the same area. To provide a single indicator of tsunami hazard at any coastal location (here, at many save points defined along the 5-m isobath), the four hazard metrics will be combined into a novel Tsunami Hazard Intensity Index (TII). To do so, a score is first attached to selected ranges, or classes, of each hazard metric (five classes are selected for each metric); the TII is then computed as a weighed average of these scores, thus providing a single tsunami hazard intensity value at each location. Similar indices have been proposed in earlier work, but using a smaller number of metrics, and shown to be useful to discriminate between low, medium, and high coastal tsunami hazard areas (e.g., Nemati et al. 2019;Boschetti et al. 2020).
In the following, we first detail recent studies of seismic sources in the NAOB and, on this basis, propose our tsunami source selection and parameterization for the two considered areas (ACZ and PRT). We then present the modeling methodology and its results in terms of the various tsunami hazard metrics computed along the USEC. Finally, based on these, we compute and discuss the TII values obtained for each source.

Overview
In this section, we define fault/subfault parameters and compute the coseismic deformation for a collection of ACZ and PRT sources, based on earlier published work and, for the PRT and Caribbean arc area, on the fault segmentation/parameterization established during the May 2019 workshop of expert at the USGS Powell Center. In the ACZ, since little information is available on the detailed nature of active faults and their parameters, coseismic sources are simply represented by a single fault plane. In contrast, in the PRT 1 3 and Caribbean arc area, which is part of subduction zones, so-called SIFT subfaults and their parameters were defined by Gica et al. (2008); hence, these are used here to define the coseismic sources as a series of subfaults in each identified segment of the considered area.
While actual events would have a non-uniform slip distribution over each source area, here, in the absence of specific information, we simply use a uniform slip S for each source. Using a non-uniform slip would cause differences in tsunami generation that could significantly affect the near-field tsunami hazard. In the far-field, however, as we shall see, the sensitivity of tsunami waves (and their coastal hazard) to detail the source, gradually decreases, the larger the distance from the source. Hence, since our main goal is to assess the far-field tsunami hazard along the USEC, caused by each selected source, we simply use a uniform slip S for each source. In the context of a future PTHA, however, considering multiple sources with randomly defined slip distributions would be important to quantify the resulting uncertainty on tsunami generation and far-field hazard.
As is standard in the parameterization of coseismic sources, fault slip S is related to the event seismic moment M o as (e.g., Grilli et al. 2010), where (W n , L n ) are the width and length of each of N subfault planes defining the source, is the medium material constant (i.e., Coulomb/shear modulus; an average value of 40 GPa is used here in all cases; Gica et al. (2008)), and the source magnitude M is defined as, Additionally, for each subfault plane in the source ( n = 1, … , N ), we specify the centroid location ( x 0,n , y 0,n ) (Lon.-Lat.) and depth d 0,n , and the three angles orientating the fault plane: dip n , rake n , and strike n ; with these definitions, d o,n = d n + (W n ∕2) sin n , where d n is the depth at the highest point on the subfault plane.
Based on these fault/subfault plane parameters, Okada (1985)'s method is used to compute the resulting coseismic seafloor deformation (uplift/subsidence), that is used as surface elevation to initialize the tsunami propagation model (see details later). This method solves a set of linear three-dimensional elasticity problems, for a series of homogeneous half-spaces corresponding to the subfaults, each with a dislocation specified over an oblique plane. The problems are expressed as a set of boundary integral equations, which are solved over a specific Cartesian grid encompassing the subfaults (here with a 1-km resolution); due to the linearity of the equations, the total seafloor deformation is computed as the sum of that caused by each subfault.

ACZ/LSB area
The tectonic setting of the ACZ/LSB area is complex, with several active faults (e.g., Figs. 3 and 4 in Baptista 2020), and has been the object of many studies, particularly with the aim of identifying the source and parameters of the 1755 event. The reader is referred to those works for more information and figures (e.g., Baptista et al. 1998aBaptista et al. , b, 2003Barkan et al. 2009;Baptista et al. 2016;Baptista 2020). Udias et al. (2020) analyzed historical information on the occurrence of earthquakes and tsunamis in the region of Cape Saint Vincent-Gulf of Cádiz and identified large tsunamigenic earthquakes that occurred SW of Iberia, in the ACZ/LSB area ( Fig. 1), before the great Lisbon 1755 earthquake. Separating events that occurred before and after 500 A.D, the authors concluded that there had been other large events during this period in this area, i.e., in 241/216 B.C., 881, 1356 and 1531. Hence, LSB 1755 is merely the most recent in a series of recurring large events in the ACZ/LSB area and there is thus a high likelihood for a similarly large event to occur in the future in this area. The exact location and parameters of the 1755 Lisbon earthquake are still unknown and subject to debate. Various studies have placed its magnitude in the M8.6 to 9 range and its most likely location in the Horseshoe Plain thrust fault (HPTF) area, to the East of the ACZ and Madeira Torre Rise (MTR; Fig. 1; Muir-Wood and Mignan 2009; Barkan et al. 2009), an underwater ridge that likely caused westward propagating tsunami waves to be somewhat diverted from directly aiming at the USEC; hence, the MTR may have offered some level of protection to the coast. By performing backward ray tracing from locations of known arrival of the Lisbon 1755 tsunami, Baptista et al. (1998aBaptista et al. ( , 1998b) located its source just north of the HPTF, in the area between the Gorringe Bank Fault (GBF) and the southwestern end of the Portuguese coast. Based on field surveys and tsunami observations, Baptista et al. (2003) later proposed a composite source in the Marquês de Pombal fault (MPF) and Guadalquivir Bank faults (GuBF), which are both located in the same area north of the HPTF. Consistent with this earlier work, Barkan et al. (2009) considered three potential sources for the Lisbon 1755 event, respectively, located in the: (i) GBF; (ii) MPF; and (iii) Gulf of Cadiz Fault (GCF) (southeast of the HPTF). Based on a set of 16 potential tsunami source simulations performed in the area of these three major faults, they inferred that the most likely source would have been located in the HPTF, which has a NW/SE strike orientation; however, this fault may just be a paleo plate boundary (Baptista 2020). Omira et al. (2009) simulated tsunamis from sources located in (i)-(iii), plus the Portimao Bank, and showed that a source in the GBF would have radiated most of its energy towards the NE and Morocco, with a minor impact along the Gulf of Cadiz. They concluded that the most likely location was in the HPTF (for actual locations of these faults see figures in Barkan et al. 2009;Omira et al. 2009;Baptista 2020).
There have been a few recent tsunamigenic events in the ACZ area, which, although they did not generate large tsunamis, confirmed that large future earthquakes could potentially occur on both sides of the ACZ/MTR area. Hence, in comprehensive tsunami hazard assessment studies, both of these locations should be considered for siting similar future events. Specifically, in 1941, a M8.3 strike-slip earthquake occurred in the Gloria fault, northwest of the ACZ. A small tsunami was registered at the tide stations of Cascais, Lagos, Portugal, Morocco, Madeira, Azores (and in the UK), with a maximum observed height of 0.45 m (peak to peak) in Casablanca in Morocco. In 1969, a M7.9 earthquake occurred in the Horseshoe Plain, which created a small tsunami with a maximum height of 0.06 m measured on the Portuguese coast in Lagos. And, in 1975, a M7.9 earthquake occurred south of the Gloria Fault, within the ACZ, southeast of the MTR ( Fig. 1; Baptista et al. 2016;Baptista 2020), which generated a tsunami recorded at a set of coastal tide gauges, with a amplitude of up to 0.3 m measured in Lagos.
Based on the conclusions of Barkan et al. (2009)'s study, as part of work done for NTHMP, Grilli and Grilli (2013a) designed and modeled a set of extreme PMTs originated in the ACZ area, as repeats of the Lisbon 1755 event. To cover the range of uncertainty in source location and parameters, they modeled twelve M9 sources, each with a single fault plane, and sited them at various locations in the ACZ. These source parameters were selected based on earlier published work, such as Barkan et al. (2009)'s study. As in the latter work, to maximize tsunami generation, each source's fault plane dip and rake angles were set to = 40 ° and = 90 °, respectively, and depth to d = 5 km. Other parameters, such as fault plane dimensions (W, L), slip S, and strike angle , were varied for each source. Strike angle, in particular, which strongly affects the tsunami initial directionality, was widely varied to identify which angles caused maximum hazard on various sections of the USEC. Specifically, values of = 15 . 30, 195, 345 and 360° were specified, and it was found that = 15 and 345° caused maximum tsunami hazard on the upper and lower USEC, respectively, with other angles being redundant. Moreover, as expected from earlier work, among all the sources, Grilli and Grilli (2013a) found that those sited in the MTR caused a larger tsunami hazard along the USEC, and the largest maximum flow depth at the coast, ∼1-2 m, occurred for the sources that had the smaller fault plane area among those suggested by Barkan et al. (2009) for the M9 sources, i.e., W = 126 km and L = 317 km. Finally, in all cases, the tsunami travel time to the USEC was between 8 and 12 h, from north to south.
According to this preliminary work, in recent modeling of coastal tsunami hazard caused by PMTs along the USEC as part of work done for NTHMP, Schambach et al. (2020b) considered four (4) so-called LSB M9 PMTs, sited at locations both west of the MTR (M9.0-MTR1 and M9.0-MTR2) and east of the MTR (M9.0-HSP1 and M9.0-HSP2) in the Horseshoe Plain. In the absence of more detailed data on fault plane parameters, each source was made of a single fault plane ( N = 1 ), with parameters based on Barkan et al. (2009) and Grilli and Grilli (2013a), i.e., = 15 or 345°, W = 126 km, L = 317 km, d = 5 km, = 40 °, and = 90 °, which using Eqs. (1-2) yielded S = 20 m (Table 1). In combination, the PMTs triggered by these sources were expected to cause maximum coastal hazard along the USEC, from tsunamis originating in the ACZ/LSB area.
Here, in addition to these extreme sources, and to prepare for future PTHA work along the USEC, we parameterize and model six (6) additional tsunamigenic coseismic sources Table 1 Parameters of coseismic source collection modeled in ACZ/LSB area ( Fig. 1) Consistent with earlier work (Barkan et al. 2009;Grilli and Grilli 2013a), all sources are defined as a single-shallow fault plane ( N = 1 ), with depth d = 5 km at the mid-fault highest point, dip angle = 40 °, and rake angle = 90 °; other parameters are given in table, with slip S computed from magnitude M based on a Coulomb modulus = 40 GPa using Eqs. (1-2). Estimates of return period T r for each event are provided for information only and have a large uncertainty. Figure 2 shows the initial elevation computed for each source. The M9 sources are Lisbon 1755 proxy sources, while the smaller magnitude M8, 8.3 and 8.7 are modeled to prepare for future PTHA along the USEC. Except for M9-HSP1 and M9-HSP2, which are sited in the Horseshoe Plain at the likeliest location for the 1755 event, all sources are centered at the same location west of the MTR (Fig. 1) in the ACZ/LSB area, using 3 smaller magnitudes, M8, 8.3, and 8.7, and 2 strike angles = 15 or 345°. These smaller magnitude sources also have a single fault plane, with dimensions proportionally reduced, compared to the M9 sources, derived based on the various references listed before (e.g., Barkan et al. 2009); using Eqs. (1-2) this yields a slip S = 2.87 , 5.95, and 18.14 m, for each source, respectively (see Table 1). To maximize tsunami generation, as for the M9 sources, these additional sources have a shallow fault plane depth d = 5 km, and dip and rake angles = 40 ° and = 90 °, respectively. Finally, to maximize the tsunami hazard on the USEC, they are all sited at the same location west of the MTR as two of the M9 Lisbon proxy sources.
As indicated before and acknowledged in the literature, there is little information and hence large uncertainty on the ACZ/LSB sources' fault plane parameters and locations. However, as the ACZ/LSB area is quite distant from the USEC, tsunami propagation and refraction over the entire NAOB width, and more importantly nearshore, will make the effects of the selected fault plane parameter values relatively less important for tsunami hazard in any particular area, than the controlling effect of the wide USEC shelf (e.g., Tehranirad et al. 2015). Nevertheless, should a proper PTHA be carried out in future work, it could include more sources, with random variations of the fault parameters around values selected in this study.
The fault plane footprints of the ten ACZ sources and their corresponding coseismic deformation are plotted in Fig. 2.

PRT and Caribbean arc area
When our NTHMP USEC work was initiated in 2010, there had been many studies regarding the seismo-tectonic context and past near-field tsunamis that had affected the PRT and Caribbean arc area (e.g., Calais et al. 1995;Lander 1997;Zahibo and Pelinovsky 2001;ten Brink et al. 2008;Harbitz et al. 2012;Mercado and McCann 1998) (see Fig. 3). However, no comprehensive study had considered the largest tsunamigenic coseismic sources that could affect the USEC from this area, except perhaps in some work performed for the US nuclear regulatory commission (ten Brink et al. 2008Brink et al. , 2011. In early work done to estimate tsunami hazard on the USEC, in the wake of the 2004 Indian Ocean Tsunami and in the absence of specific guidance, Knight (2006) had modeled as an extreme PMT for the USEC a M9.1 source parameterized as a single fault plane in the PRT (Segments #2 and #3 in Fig. 3), which is a 600-km-long deep trench that runs parallels to the Puerto Rico North shore. In other work also predating our NTHMP study, Grilli et al. (2010) modeled this and another similar source with a smaller M8.7 magnitude, to assess both the nearfield tsunami hazard on Puerto Rico and the far-field hazard on the USEC. Based on the work of DeMets (1993), Mercado and McCann (1998), and Zahibo and Pelinovsky (2001), they assumed a predominantly (lateral) strike-slip motion of the Caribbean plate at 20 mm per year with respect to the North American Plate, in an ENE direction (a 10-20 degree angle with respect to the PRT axis). Analyzing the historical earthquakes and the tsunamigenic events among those, they noted that 12 earthquakes of at least a M7 magnitude had occurred in and near the PRT area in the past 500 years (Dawicki 2005), with two of these having a M8.1 magnitude, and three having generated a tsunami with a 5-7-m runup on Puerto Rico. Combining these observations with the plate convergence rate, they estimated that a M7.5-8.1 event in the PRT would have an 80 year return period, a M8.7 event at least a 200 year return period, and a M9 event at least a 600 year return period. In view of more 1 3 recent work on the segmentation and fault locking in the area (see below), it is likely that these return periods were under-estimated. Simulating Knight's extreme M9 event, Grilli et al. (2010) showed that the generated tsunami caused up to a ∼2-3-m flow depth along the USEC shore and tsunami travel times were between 2.5 and 6 h, depending on the location from south to north.
In our subsequent NTHMP first-generation work, in light of new data provided on subfault plane parameters in the PRT/Caribbean arc area as part of NOAA's SIFT dataset (Gica et al. 2008), Grilli and Grilli (2013b) designed and modeled a different M9 PMT source, which approximately covered the same area of the PRT as Knight (2006)'s source, i.e., 600-km long by 150-km wide, but was parameterized as 6 by 2 (i.e., N = 12 ) individual SIFT subfaults, each 100-km long by 50 km wide (in the oblique fault plane), which better represented the convex geometry of the PRT. Although this was an extreme scenario, particularly for Puerto Rico, our work on the 2004 Indian Ocean M9.2 earthquake and devastating tsunami Ioualalen et al. 2007) had convinced us, and others (ten Brink et al. 2008Brink et al. , 2011) that a large megathrust event could occur in the PRT because of the similarities between the Puerto Rico and Sumatra-Andaman trench geometry (both trenches are arched) and plate dynamics. This was also supported by our work on the Tohoku 2011 M9 earthquake and tsunami that occurred in the similar size Japan Trench (Grilli et al. 2013;Kirby et al. 2013;Tappin et al. 2014). In their recent, higher-resolution simulations of tsunami hazard caused by PMTs on the USEC, Schambach et al. (2020b) used as one of their sources the PRT PMT of Grilli and Grilli (2013b), and confirmed the range of runup and travel time values predicted on the USEC by Grilli et al. (2010) and Grilli and Grilli (2013b) in their earlier work.
In the present work, besides the earlier PMT, we design and model a collection of hypothetical, but realistic, coseismic tsunami sources of various magnitude in the PRT/Caribbean arc area that would significantly impact the USEC. As indicated before, a meeting of experts was organized in May 2019 at the USGS Powell Center, attended by the lead author, whose agenda was devoted in large part to establishing a logic tree of coseismic sources in the PRT and Caribbean arc area, in preparation for future PTHA work on the USEC. Such a tree defines the fault plane parameter ranges for each source, while attaching a probability to each branch of the tree. The workshop approach was based on the Delphi method, which is a process used to arrive at a group opinion by surveying the experts attending the venue. Through several rounds of questionnaires, the responses provided by the experts on various source parameter values were transformed into tables of probability for classes of parameter ranges.
As a first step, considering the existing seismo-tectonic knowledge in the area, the experts attending the meeting established a realistic segmentation of the entire Caribbean arc and estimated each segment's parameter ranges. The defined segments were those most likely to fail either individually or together in clusters, during seismic events of various magnitude in the considered area. Coseismic tsunami sources of various magnitude could then be designed by combining various segments. Then, the experts Note, the M9.0 sources represent Lisbon 1755 proxies, located in either the ACZ, west of the MTR ("MTR"), or Horseshoe Plain areas (Fig. 1). Color scale is surface elevation (in meter), with same scale used in all plots for comparison; black contours denote depth (in meter) established realistic ranges for the level of plate locking and magnitude of the faultnormal convergence rate for each segment that would most contribute to causing large tsunamigenic earthquakes (Benford et al. 2012), and noted in particular that the level of plate locking had a large uncertainty in the PRT area. Figure 3 shows the likeliest segmentation established during the meeting, which is composed of five segments that include various faults, encompassing subduction zones in the Hispaniola trench on the west (Segment #1), the PRT in the middle (Segments #2,3), and in the northern Lesser Antilles arc (Segments #4,5). Note, additional information regarding the Hispaniola sector (Segment #1) can be found in Rodríguez-Zurrunero et al. (2020). Each of these segments approximately overlaps with subfaults defined in the SIFT dataset (Gica et al. 2008) (see insert in Fig. 3). For each SIFT subfault, of dimension L = 100 by W = 50 km in the trench-parallel and trench-normal directions, respectively, the dataset provides the strike and dip angles, with the rake angle assumed to be = 90 ° in all cases (to maximize the seafloor deformation), and the lat-lon coordinates of the fault centroid, and depth d of the fault highest point. Table 2 provides these parameters for all the SIFT subfaults overlapping with Segments #1 to 5 in Fig. 3, with the subfault locations marked in the figure inset.
The second step considered by the experts was to define the likeliest grouping of individual segments that would fail together in a single event, thus creating sources from smaller to larger magnitude, as well as the associated magnitude for each such event and their likelihood (i.e., probability). Although the final proceedings of the workshop are still pending, based on draft material from the workshop, we defined  (Gica et al. 2008) in the considered area (see Table 2 for parameter values). Large black arrows show the nearly E-W relative plate motion of the North America Plate subducting under the Caribbean Plate. Labels locate: A: Septentrional Fault (SF); B: Mona Passage; C: Guadeloupe (see Fig. 1 for additional location information). Note, the basis for part of this figure is a public domain USGS figure (https:// www. usgs. gov/ media/ images/ locat ionearth quakes-north easte rn-carib bean), that was modified by U. tenBrink for presentation in the workshop, and then to mark the established segmentation; the authors further modified this figure for use in this work.
[See Rodríguez-Zurrunero et al. (2020) for additional information on the Hispaniola sector segmentation.] seven (7) new sources in the PRT/Caribbean arc area, with magnitudes M8.3, M8.7, and M9.0, regrouping two to four segments from Fig. 3 (Table 3); an additional eighth source, M9.0-PRT3, is the extreme PMT source used in earlier work to date (Grilli and Grilli 2013b;Schambach et al. 2020b). Based on the SIFT subfault dimensions and the number N = 10 − 26 of subfaults used for each source listed in Table 3, an associated average slip S was computed using Eqs. 1 and 2, assuming = 40 GPa as recommended in the SIFT study.

3
The fault plane footprints of the 8 PRT/Caribbean Arc area sources and their corresponding coseismic deformation are plotted in Fig. 4.

Estimate of return periods of the ACZ and PRT/Caribbean arc sources
As discussed above, only very limited historical data is available regarding tsunamigenic coseismic sources in our two selected areas, which prevents performing a rigorous  Table 2. Color scale is surface elevation (in meter), with same scale used in plots a-g, for comparison; black contours denote depth (in meter). Note, M9.0-PRT3 is the extreme PMT source modeled in earlier work statistical analysis. Accordingly, in the following, we provide first-order estimates of the return periods associated with the ACZ and PRT/Caribbean Arc coseismic tsunami sources (see Tables 1 and 3 for the sources and their parameters), based on the limited statistical data available and on a characteristic magnitude-frequency recurrence model (with a plate coupling of 1); hence, the estimated return periods, which are also listed in the tables, have large uncertainty. More rigorous analyses and probabilistic considerations will be addressed in future PTHA work.
The ACZ/LSB area is a convergence zone featuring multiple faults, where only a few historical tsunamigenic events were reported. Hence, accurate return period estimates cannot be made for each of the sources in Table 1. For the largest LSB 1755 type events, however, it appears that the average recurrence is on the order of 500 years (Udias et al. 2020) (4-5 medium-large events in 2000 years). Hence, prorating the slip S and the square root of fault area A = L W of the smaller events to those of the M9 event, one would infer that M8, M8.3, and M8.7 events would have a return periods of 35, 80, and 285 years. However, accounting for historical observations of the two smaller magnitude events and the uncertainty of this estimate, one would propose the range of return periods listed in Table 1 for each event.
In the PRT/Caribbean Arc, if the maximum convergence rate of ∼ 20 mm/y (i.e., 1 m of slip in 50 years) estimated for the strike-slip motion of the Caribbean plate with respect to the North American Plate (DeMets 1993; Mercado and McCann 1998;Zahibo and Pelinovsky 2001) contributed to the slip values listed in Table 3, the return period for each source would be ∼75, 185, and 355 years, respectively. However, the relative plate subduction is highly oblique (10-20°) on the western side of the considered area and less so in the lesser Antilles (Fig. 3), and the accumulated fault slip must be multiplied by the sine of the convergence angle. Hence, assuming an average of 20° for the entire area yields a factor of Table 3 Parameters of coseismic source collection modeled in PRT/Caribbean Arc area (Fig. 1) Sources are defined in terms of segments selected during the USGS 2019 Powell Center workshop and corresponding SIFT subfault planes (Fig. 3); their seafloor deformation is modeled with Okada (1985)'s method to define a tsunami source. Slip S is computed from magnitude M based on a Coulomb modulus = 40 GPa using Eqs. (1-2); other SIFT source parameters used in each segment are given in Table 2. Estimates of return period T r for each events are provided for information only and have a large uncertainty. Figure 4 shows the initial surface elevations computed for these sources. The M9.0-PRT3 source is the extreme PMT source used in earlier work (Grilli et al. 2010;Grilli and Grilli 2013b;Schambach et al. 2020b one-third or so and approximate return periods on the order of 250, 550, and 1,000 years for each source magnitude, respectively, which is consistent with estimates made during the 2019 Powell workshop; applying the same considerations to the M9.0-PRT3 source, we find an approximate return period of 2,000 years (i.e., about 3 times the earlier estimate of Grilli et al. (2010) that did not consider the oblique subduction).

Overview
For each of the selected coseismic sources, tsunami propagation to the USEC is modeled using the nonlinear and dispersive two-dimensional (2D) Boussinesq long wave model (BM) FUNWAVE (Wei et al. 1995), in a series of nested grids of increasing resolution, by a one-way coupling method. We use FUNWAVE-TVD V.3, the newer implementation of FUNWAVE, which is fully nonlinear in Cartesian grids (Shi et al. 2012) and weakly nonlinear in spherical grids (Kirby et al. 2013). The model was efficiently parallelized for use on a shared memory cluster, with a reported 90% scalability (Shi et al. 2012), which allows using large grids; the present simulations were run using 100-200 processors on NSF-XSEDE HPC resources, in a reasonable computational time of a few up to 8-10 hours each. Note, FUNWAVE-TVD now also has an efficient multi-GPU implementation that allows using such hardware that typically feature more than 4,000 computing cores each (Yuan et al. 2020). FUNWAVE and then FUNWAVE-TVD are open-source codes available on GitHub that have been widely used to simulate tsunami case studies (e.g., Grilli et al. 2007Grilli et al. , 2010Grilli et al. , 2013Grilli et al. , 2016Grilli et al. , 2017aGrilli et al. , 2019Ioualalen et al. 2007;Tappin et al. 2008Tappin et al. , 2014Abadie et al. 2012;Kirby et al. 2013;Tehranirad et al. 2015;Schambach et al. 2020aSchambach et al. , 2021. As discussed in the introduction, since 2010, the authors have used FUN-WAVE and related methodology as part of NTHMP work to simulate coastal hazard from PMTs and compute tsunami inundation maps along the USEC (see, http:// china cat. coast al. udel. edu/ nthmp. html, e.g., Abadie et al. 2012;Grilli and Grilli 2013a, b;Grilli et al. , 2017aTehranirad et al. 2015;Schambach et al. 2019). The same modeling approach was also used to perform several other tsunami hazard assessment studies of coastal nuclear power plants in the U.S. Both spherical and Cartesian versions of FUNWAVE-TVD were validated through benchmarking and, through this process, officially approved by NOAA for NTHMP work (Tehranirad et al. 2011;Horrillo et al. 2015;Lynett et al. 2017).

Frequency dispersion
As mentioned above, both FUNWAVE (Wei et al. 1995) and its more recent implementation FUNWAVE-TVD (Shi et al. 2012;Kirby et al. 2013) used here are fully nonlinear Boussinesq long wave models that feature extended physics compared to the more commonly used tsunami models based on the non-dispersive nonlinear shallow water (NSW) equations. While FUNWAVE's nonlinear properties are the same as those of NSW models, it features extended frequency dispersion effects (with respect to the original weakly nonlinear Boussinesq equations, e.g., Peregrine 1967), which allows accurately simulating linear dispersion properties up to nearly the deep water wave limit (i.e., where the local water depth h is equal to half the wavelength ); for long tsunami waves, this adequately covers most oceanic water depths. To achieve such extended dispersion properties, based on earlier work by Nwogu (1993), the horizontal velocity components (u, v) computed as a function of time in each cell of the model grid, together with the surface elevation , are defined at 0.531 times the local depth, while the horizontal velocity profile is quadratic over depth (instead of uniform in the NSW models). Using a dispersive long wave model is necessary to accurately simulate landslide tsunamis (e.g., Glimsdal et al. 2013;Schambach et al. 2019), even in the near-field, as these are typically made of shorter and hence more dispersive waves than coseismic tsunamis (Watts et al. 2003). However, dispersion is also important for modeling the propagation over long distances of coseismic tsunamis (Kirby et al. 2013), or for that matter of any tsunami. Indeed, as shown by Glimsdal et al. (2013), the cumulative effects of dispersion gradually become larger, the larger the propagation distance L. Assuming periodic waves propagating over constant depth, they quantified dispersive effects by way of a non-dimensional parameter = 6h 2 L∕ 3 and showed that these become significant when > 0.1 . For instance, considering an average ocean depth h = 3 km and a distance of propagation L = 3000 km, this implies that any wave of wavelength < 117 km will feature significant dispersive effects (this also corresponds to a long wave period T = ∕ √ gh ≤ 11.4 min). These equations also imply that the wavelength or period threshold for dispersive effects to become significant is ∝ L 1∕3 . Hence, for the same average depth, doubling or halving the propagation distance L increases or reduces those thresholds by 25% or 21%, respectively. In the present study, the USEC is in the far-field of both coseismic tsunamis triggered in the ACZ or the PRT/Caribbean Arc areas (Fig. 1), with distances of propagation from the source of L = [4600 − 6300] km and [1300-2700] km, for each source, respectively (calculated using Google Earth TM 's ruler tool); for an average depth of 3 km, this yields ranges of wave periods thresholds below which dispersive effects are important of T ≃ [13.2 − 14.6] min and [8.6-11.0] min, respectively. As we shall see in simulation results, time series of incident tsunami wave trains on the USEC shelf will feature many waves that based on these thresholds are short enough to have been affected by dispersion. As shorter waves have different phase speeds in a dispersive model, they can combine through constructive or Table 4 Parameters of model grids used in FUNWAVE to compute the propagation of far-field coseismic sources from the ACZ/LSB and PRT areas (Figs. 1, 5) "Res." refers to the resolution of spherical (S)-or Cartesian (C)-type grids, and N x and N y indicate the number of grid cells in each direction. Letters in parenthesis indicate whether coordinates are for: (SW) southwest corner, (NW) corner, or (CT) center of grid. Angle refers to the orientation (trigonometrically from east) of the long side of the rotated grids G1, G2 and G3  (Tappin et al. 2008), and similarly for the 2018 Palu landslide/coseismic tsunami (Schambach et al. 2021). For coseismic tsunamis, although a large part of the tsunami modeling community is still using non-dispersive NSW models, we confirmed long ago and more recently an increasingly large number of tsunami case studies showed, for actual events, that dispersive effects, as would be expected from the above theoretical considerations, can cause significant differences in the elevation of successive waves in a tsunami wave train, particularly over long distances of propagation. For instance, Ioualalen et al. (2007) compared dispersive versus non-dispersive simulations with FUNWAVE for the impact of the 2004 Indian Ocean tsunami in Thailand, and showed up to 25% differences in maximum surface elevations. A similar comparison was made by Kirby et al. (2013) for the 2011 Tohoku tsunami, using the spherical version of FUNWAVE-TVD that showed up to 60% differences in maximum tsunami elevation in the far-field between dispersive and non-dispersive simulations; in fact, dispersive effects were much larger than Coriolis effects, which only caused up to a 5% difference on surface elevations. In both cases, there were areas were tsunami elevations decreased and others where they increased, without any predictable pattern, with the largest wave in the tsunami train often becoming the second or third wave, rather than the leading wave as is typically predicted using NSW models. A recent study in the Japan Trench even confirmed the importance of simulating dispersion for the near-field impact of coseismic tsunamis caused by outer-rise normal faults (Baba et al. 2021).
A final and important aspect of dispersion for coastal tsunami hazard assessment, which is not illustrated in the present study as it only uses fairly coarse grids, is the modeling of undular bores that can be generated nearshore near the crest of incident long waves (Madsen et al. 2008;Schambach et al. 2019). For instance, using many levels of nested coastal grids with FUNWAVE (down to a 7.5 m resolution), Schambach et al. (2019) showed that undular bores increased tsunami surface elevations, current magnitude, and momentum forces at the shore by up to 25-50%. This increase is in addition to dispersive effects that may have occurred during tsunami propagation. 1 3

Shoreline algorithm and nearshore dissipation
Along the shore, FUNWAVE has an accurate moving shoreline algorithm that identifies wet and dry grid cells, allowing to dynamically compute tsunami inundation and runup. Although in view of the moderate resolution of the coastal grids, we do not consider these tsunami hazard metrics here, having an accurate simulation of these processes in the model is nevertheless important to prevent excessive reflection or numerical dissipation at the coast, which would affect the tsunami metrics computed along the 5-m isobath. Regarding actual dissipation simulated in FUNWAVE, bottom friction is modeled in all grid cells by a quadratic law, but of course only becomes significant when tsunami waves propagate over the shallow USEC shelf. Although the bottom friction coefficient can and has been adjusted in the model, as a function of geology and land use, in higher-resolution simulations (e.g., Grilli et al. 2017a), here, considering the coarser grids and the absence of site specific data, we use the standard constant value C d = 0.0025 , which corresponds to coarse sand. Tehranirad et al. (2015) have shown that over the wide USEC shelf, bottom friction can significantly reduce the incident tsunami elevations; hence, the moderately low  Table 1), at t = a 1, b 2, c 5, d 6.5, e 8, f 9.5, g 10.5, and h 11 h. Some isobaths are plotted for reference, but without labels to simplify the figures. Higher-resolution results are used wherever available value used here is conservative as far as coastal tsunami hazard. Dissipation from breaking waves is simulated in FUNWAVE using a front tracking (TVD) algorithm and switching to NSW equations in the model grid cells where breaking is detected (based on a standard breaking criterion). Earlier work has shown that numerical dissipation in NSW models closely approximates the physical dissipation in breaking waves (Shi et al. 2012). An accurate simulation of the nearshore dissipation due to breaking waves is also important for the accurate modeling of both incident tsunami waves and their reflection at the coast, which both affect the tsunami hazard metrics computed along the 5 m-isobath.  Table 1 and initial surface elevations shown in Fig. 2i, j). Color scale is surface elevation in meter

Nested model grids
Simulations of tsunami propagation and coastal impact are performed in nested grids ( Figs. 1 and 5), by one-way coupling. In this method, time series of surface elevations and depth-averaged currents are computed at a large number of stations/numerical wave  (Table 3; Fig. 4h), at t = a 0.5, b 2, c 3, d 4, e 4.5, and f 5 h. Some isobaths are plotted for reference, but without labels to simplify the figures. Higher-resolution results are used wherever available gages, defined in a coarser grid, along the boundary of the finer grid used in the next level of nesting. Computations are performed in the coarser grid and then are restarted in the finer grid, using the station time series as boundary conditions. As the latter include both incident and reflected waves computed in the coarser grid, this method closely approximates open boundary conditions. It was found that a nesting ratio with a factor 3-4 reduction in mesh size allowed achieving good accuracy in tsunami simulations. Finally, to prevent reflection in the first grid level, sponge layers are used along all the offshore boundaries (see, e.g., Tappin et al. 2014;Schambach et al. 2019Schambach et al. , 2020aSchambach et al. , 2021, for a few recent examples).
For each coseismic source, simulations are first run for 24 h in a 1 arc-min (about 1800 m) resolution grid, in which the tsunami source elevation is specified as an initial condition at t = 0 , with a zero velocity as is standard in such simulations. For the ACZ/ LSB sources, this is the ocean basin scale spherical grid "Large G0", which covers the footprint of Fig. 1, and for the PRT/Caribbean arc sources, this is the smaller "Local G0" (Fig. 1; Table 4). In both grids, 200-km wide sponge layers are specified along the west, south and north offshore boundaries, and 400-km wide on the east boundary, to eliminate reflection. To compute the tsunami hazard along the USEC, simulations are restarted in three regional, angled, Cartesian nested grids with a 450-m resolution, G1, G2, and G3 (Figs. 1, 5; Table 4), also up to the full 24h duration, but for more efficiency starting from the time of first arrival of the tsunami along each nested grid boundary.
For each grid level, bathymetry and topography are interpolated from the most accurate data available, commensurate with each grid resolution. In deeper water, we use NOAA's 1 arc-min ETOPO-1 data (Fig. 1); in shallower water and on the continental shelf, we use NOAA's NGDC 3" (about 90 m) Coastal Relief Model data (Fig. 5). Save  (Table 4; Fig. 1) for the M9-PRT3 coseismic source (with parameters listed in Table 3 and initial surface elevations shown in Fig. 4h). Color scale is surface elevation in meter 1 3 points are defined in each grid (P-G1, P-G2, and P-G3 in Table 5; Fig. 1), where time series of surface elevations and currents are computed in order to verify that results in nested grids (G1, G2, G3) are consistent with those in the coarser parent grid (G0). Additional numerical wave gauges (P1 to P6 in Table 5; Fig. 1) are specified to visualize  (Table 1), at the 9 save points defined in Table 5 (shown in Fig. 1): (blue lines) in grid Large G0; (red lines) in nested grids G1-G3. For point P4, yellow indicates surface elevation computed in G2, red indicates surface elevation computed in G3 1 3 details of incident tsunami wave trains (i.e., number of waves, wave amplitudes and periods), with P1 to P5 being located along the 200-m isobath, at the continental shelf edge, and P6 being in 1,000-m depth (southernmost point north of the Great Bahama Banks off of Florida).  (Table 3), at the 9 save points defined in Table 5 (shown in Fig. 1): (blue lines) in grid Local G0; (red lines) in nested grids G1-G3. For point P4, black indicates surface elevation computed in G2, red indicates surface elevation computed in G3 Note that transformations from (Lon.-Lat.) to Cartesian (x, y) coordinates and back are performed using MATLAB's mapping toolbox, using a transverse Mercator projection (UTM), with an origin at each grid's centroid.

Initial tsunami elevations
Following the standard approach for large coseismic tsunamis, the initial tsunami surface elevation (with zero velocity) of each coseismic source is assumed to be equal to the maximum seafloor deformation caused by the earthquake; this is acceptable since water is nearly incompressible and raise times are usually small. Tables 1 and 2, 3 list the parameters of the 18 coseismic sources modeled in the ACZ/LSB and the PRT/Caribbean arc areas, respectively. Based on these, Fig. 2 shows the initial surface elevations for the ACZ/ LSB sources, computed as discussed before using Okada (1985)'s method. The four largest M9 sources in Fig. 2g-j cause a nearly 12 m maximum elevation and a 4-m maximum depression. The other, smaller magnitude, sources plotted in Fig. 2a-f have similar patterns of uplift and subsidence, but smaller surface elevations distributed over gradually smaller areas. Figure 4 shows the similarly computed initial surface elevations for the PRT/ Caribbean arc sources. The new M9 sources in Fig. 4f, g (PRT1 and PRT2) have an over 4-m maximum elevation and 2-m maximum depression, about half the maximum value of the M9.0-PRT3 PMT source used in earlier work (Fig. 4h). The other, smaller magnitude, sources plotted in Fig. 4a-e have similar patterns of uplift and subsidence, but gradually smaller surface elevations as their magnitude decreases.

Tsunami propagation to the USEC
For the ACZ/LSB area, Fig. 6 shows an example of surface elevations computed for the Lisbon 1755 proxy source M9-HSP1 ( Fig. 2i; Table 1) at times t = 1 to 11 h, with the figures at later times focusing on the impact computed in grid G3 off of Florida. The snapshots illustrate how waves radiate away from the source, diffracting around islands and refracting as a function of the bathymetry. As a result of frequency dispersion, nearshore tsunami wave trains are made of at least three major wave crests and troughs (e.g., Fig. 6h). Here, despite the 15° strike angle of the M9-HSP1 source, the largest tsunami waves are aimed at Florida in the far-field. This strong wave-guiding effect of the deep water bathymetry is further shown in Fig. 7, which shows the envelopes of maximum surface elevations above MWL computed for the M9-HSP1 and M9-HSP2 sources, which have the same magnitude, but strike angles = 15 and 345°, respectively. While generating quite different wave envelopes in the near-field of the source, in the far-field, both sources are directing their major wave energy southwestward, toward Florida. In addition to this main directionality, as we shall see more clearly in detailed coastal hazard results, due to their strike orientation, the first source causes a relatively larger tsunami hazard on the upper USEC, while the second source causes a larger hazard on the lower USEC and the Caribbean Islands. Finally and importantly, Fig. 6e-h illustrate how the continental slope and shelf cause intense refraction of incident tsunami waves, which increasingly bend to become nearly parallel to the isobaths as they gradually slow down. This phenomenon affects all the ACZ/LSB sources in the same manner and causes them to focus or defocus their energy on identical coastal areas, depending on whether there is a convex/concave isobath geometry, respectively, similar to effects of optical lenses for light rays. This alongshore modulation pattern can be seen in the envelopes of Fig. 7 and will also be clear in detailed coastal hazard results. A movies showing animated results of the M9-HSP1 source simulation is provided as supplementary material (HSP1.mp4).
For the PRT area, Fig. 8 shows an example of surface elevations computed computed for the M9.0-PRT3 source ( Fig. 4h; Table 3), at times t = 0.5 to 5 h, with the figures at later times focusing on grid G1, off of the upper USEC (centered in New England) where the most impacted areas are located, due to the northward directionality of tsunami energy for this source. This main northward directionality of the PRT/Caribbean Arc sources is clearly confirmed in the maximum envelope of surface elevation shown in Fig. 9 for the M9-PRT2 source. Similar to the tsunamis caused by the ACZ/LSB sources (e.g., Figs. 6), 8c-e show that an intense refraction of incident tsunami waves occurs at the shelf break and over the continental shelf, which makes wave crests become increasingly parallel to the local isobaths. As a result of frequency dispersion, the nearshore tsunami wave trains are again made of at least three major wave crests and troughs (e.g., Fig. 8e). Comparing these results to Fig. 6f-h, we see that wave focusing-defocusing caused by refraction over the wide and shallow USEC shelf, particularly around ridges and canyons in the bottom bathymetry, yields a similar along-shore modulation of tsunami coastal impact that is independent from the directionality of the incident tsunami in deeper water. Detailed results of coastal hazard, discussed later, will confirm this phenomenon, showing that besides an effect of the initial tsunami directionality, coastal hazard metrics will have similar alongshore modulation patterns for the ACZ/LSB and PRT sources, which are located eastward and southward of the USEC, respectively. Hence, where there is a wide shallow shelf, the same areas of the coast will always face relatively more or less tsunami hazard, whatever the tsunami source origin. This property was already noted in earlier work, for other tsunami sources affecting the USEC, by Grilli et al. (2017a, b); Tehranirad et al. (2015); Schambach et al. (2019). A movies showing animated results of the M9-PRT3 source simulation is provided as supplementary material (PRT3.mp4).
Finally, Figs. 10 and 11 show typical time series of surface elevation computed for the M9-HSP1 and M9.0-PRT1 sources, respectively, at the 9 save points (i.e., numerical wave gauges) defined earlier near and on the continental shelf (Fig. 1, Table 5), in the coarser grids Large G0/Local G0 and the finer nested grids G1, G2 and G3. These both confirm the relevance of the one-way coupling method and illustrate the effects of frequency dispersion on the far-field tsunami wave trains. Indeed, at most locations, there is a good agreement between the coarser and finer grid simulations for the first few hours of each tsunami wave train. Once reflected waves from the coast propagate back Table 6 Class limits [C i k,min − C i k,max [ (with k = 1, … , 5 ) of hazard intensity (low, low-medium, medium, high, to highest hazard), for four hazard metrics M i ( i = 1, … , 4 ) used to plot each hazard and compute the Tsunami Intensity Index (TII) (with Eq. (3)), along the USEC [0.5-1.5[ [1.5-2.5[ [2.5-4[ -5[ [5-10[ [10-25[ [25-50[ > 50 M 4 [0.1-0.17[ [0.17-0.25[ [0.25-0.50[ > 0.50 to the save points, however, differences become larger as reflection is more accurately computed in the finer grids, as would be expected from their better resolution of shoreline processes. Differences are largest at point P-G1, which is on the shelf in shallower water and closer to shore, east of Cape Cod. Overall, the agreement of time series computed in different grids confirms the relevance of the one-way coupling methodology in nested grids used here. Regarding dispersion, most time series of incident wave trains from both areas show many smaller period waves riding over longer period waves, with the shortest periods being 4-5 min. in both cases, i.e., clearly below the lower bound of the period threshold (13 and 9 min.) estimated above for dispersive effects to become important for tsunamis propagating to the USEC from the ACZ/LSB and PRT areas, respectively. Time series of surface elevations at save points look qualitatively similar for the other ACZ/LSB and PRT sources and will not be repeated here.

Definition of hazard metrics and classes
For each of the ten ACZ/LSB (Table 1; Fig. 2) and eight PRT/Caribbean arc (Table 3; Fig. 4) coseismic sources simulated with FUNWAVE, coastal hazard metrics are computed along the 5-m isobath (with respect to MWL) that parallels the USEC, by interpolating results computed in the three Cartesian nested grids (G1, G2, G3) at a large number of save points defined by their latitude-longitude coordinates along the isobath. The four computed tsunami hazard metrics, referred to as M i ( i = 1, … , 4) , are: (1) maximum tsunami surface elevation above MWL M 1 = max , (2) current velocity magnitude M 2 = U max , and  (Table 1; Fig. 2i): a envelope of maximum surface elevations above MWL max over grids G1, G2, G3 (color scale in meter; colored dots along the coast correspond to hazard intensity C 1 k from b); b-d color-coded hazard intensity classes C i k ( k = 1, … , 5 ; yellow, green, blue, magenta and red; Table 6) with corresponding values of hazard metrics M 1 = max , M 2 = U max , and M 3 = F max , respectively ( J = 28, 401 save points are defined on the 5 m isobath, including large bays) (3) momentum force M 3 = F max , as well as (4) M 4 = 1∕t a , which is function of tsunami travel time t a (here a large travel time corresponds to a low hazard).
These metrics quantify tsunami inundation hazard, navigational hazard nearshore and in harbors, hazard from tsunami impact forces on structures, and hazard resulting from low warning time, respectively. Note that the maximum momentum force, is the force measure used in coastal engineering standards (e.g., ASCE 2016) for tsunami related structural design. It is the force, per unit length of coast, resulting from the tsunami impact from bottom to free surface, on a 1-m wide vertical wall normal to the velocity vector. Physically, it results from the blocking of tsunami momentum by a structure, hence causing a horizontal force corresponding to the incident horizontal momentum flux (Newton's second law in a control volume approach). FUNWAVE outputs directly provide the surface elevation and the horizontal current components (u, v) (at 0.531 times the local depth) at each grid point as a function of time, based on which max , U max = max{ √ u 2 + v 2 } , and F max = max{ (h + ) U 2 } (h the undisturbed water depth with respect to MWL) are computed over the entire duration of simulations at all grid points. (Note that in shallow water, the horizontal velocity is assumed to be uniform over depth.) Tsunami travel time is calculated along the same isobath, as the time when a positive or negative surface elevation first occurs over a threshold, i.e., t a (in hours) is the minimum time such that | (t a , x j , y j ) |≥ ; where (x j , y j ) ( j = 1, … , J ) denotes the save points along the isobath and here = 0.01 m. Since tsunamis are very small amplitude waves relative to depth in most of their propagation, their celerity is well approximated as a function of depth by the linear long wave celerity, c = √ gh (with g the gravitational acceleration), which is not amplitude dependent; hence different tsunamis originated from the same area propagate similarly along the same "wave rays." This similarity of propagation to shore is further reinforced by refraction that takes place in large depth for long tsunami waves and causes each tsunami to propagate similarly over the wide USEC shelf, whatever its origin. Consequently, tsunamis of different magnitude originated from the same area, LSB/ACZ or PRT/Caribbean arc, will be found to have very similar travel times along the USEC. One caveat is, for the weakest LSB/ACZ M8 and M8.3 sources that approach areas of the USEC featuring bays and more complex shoreline geometries, with a small amplitude, and hence only reach the threshold later in their propagation, due to shoaling and reflection of the tsunami wave trains. This tends to increase the travel time for these sources in some areas. Details will be shown later. Figures 7 and 9 show examples of maximum surface elevations above MWL computed over the entire computational domain for some of the tsunami sources. These and similar results for the other metrics are then interpolated at the save points along the 5-m isobath to be used for assessing tsunami hazard at the coast. With its moving shoreline algorithm, FUNWAVE computes in real time the wet and dry cells and associated water levels (and currents/forces) in each wet cell, which in flooded areas end up being located on what was originally dry land. Here, rather than considering the largest values of each tsunami hazard metric within the inundated areas (e.g., runup), which could be affected by the low grid resolution, these are computed for the incident tsunami, as a proxy, along the 5-m isobath (based on MWL). This isobath, hence, is always located well within the wet cells, allowing for an accurate interpolation of model results to be done as a post-processing; thus: (i) we first use the MATLAB function contour to compute the 5-m isobath in the bathymetry defined with respect to MWL; then, (ii) for the large number of (lat-lon) "save" points defining this contour, we apply the MATLAB function interp2 (bi-directional cubic interpolation) over the wet area of the computational grid to interpolate tsunami metrics at each of the isobath's point. Based on results, there appears to always be enough wet points on the wide shallow shelf of the USEC to provide enough interpolated results along the isobath. With this method, the water depth and location of the selected contour never varies, independent of tsunami impact and, for instance, the maximum sea surface elevation, used as a proxy for the maximum inundation (flow) depth on the coast, is always that measured above MWL at the location of the 5-m contour line.
Hazard metrics interpolated at save points along the 5-m isobath are plotted for each source as a function of the increasing latitude along the USEC. However, because depending on the coastal bathymetry, the save points coordinates are not necessarily monotonously increasing along the isobath, in some areas with bays and inlets, multiple values may occur for the same latitude. For this reason, and to ensure clarity, in most of the results and figures, we use save points that exclude large bays (i.e., Cheasapeake and Delaware bays, and Long Island Sound; in this case there are J = 18, 201 save points); when large bays are considered, additional save points are specified within the bays (in this case, the total number of save points is J = 28, 411).
Once the four hazard metrics computed, to more easily identify areas facing lesser or larger hazard, similar to classes defined by Boschetti et al. (2020) for the first two metrics in their tsunami intensity scale, we define five intensity classes, referred to as C i k ;k = 1, … , 5 for each of the four metrics ( i = 1, … , 4 ), with ranges of values [C i k,min − C i k,max [ from low to increasing hazard intensity (or severity; see Table 6). These classes of hazard for each metric can also be referred to as: low, low-medium, medium, high, and highest hazard. As discussed in Boschetti et al. (2020), who reviewed other relevant work to date, the selected values for maximum inundation max (here estimated by its proxy value computed along the 5-m isobath), correspond, for adult pedestrians, to being up to knee-tight deep for class C 1 1 , up to chest/head deep for class C 1 2 , up to head to overhead deep for C 1 3 while losing the ability to feel the terrain, and then very deep; in classes C 1 3 and higher, people would be forced to swim or have to find high ground or vertically evacuate to be safe. Likewise, for maximum currents U max (see also, Lynett et al. 2017), adult pedestrians would only be able to fight the current in classes C 2 1 and C 2 2 , while navigation would start being gradually and then severely impended nearshore and in harbors for classes C 2 2 and above. Maximum momentum force F max classes range from forces for which most structures are resisting the tsunami impact, up to forces for which most structures suffer significant damage or Table 7 Statistics of simulation results computed at J = 18, 201 save points along the 5-m isobath (excluding large bays) for LSB/ACZ sources (Table 1; Fig. 14): mean, standard deviation, root-mean-square and average of the top 33, 10, 1 and 0.1 percentiles of maximum elevation (m), flow velocity U (m/s), and momentum force F (kN/m) (m) destruction, except for the strongest concrete or steel-built (an/or elevated) structures (e.g., ASCE 2016). The F max class thresholds were then selected to roughly represent percentile classes based on a statistical analysis of model results at save points, assuming a lognormal distribution. Finally, travel time classes were selected based on results of pedestrian evacuation models for developed coastal areas (e.g., Wood et al. 2016), and discussions with emergency managers. They reflect having a large enough time ( t a = 10 h) to evacuate most of the population from high hazard areas to not having enough time ( t a = 2 h) for warning and evacuating a meaningful fraction of the population at risk.

Overall hazard metrics and classes along the USEC
For each source considered here, values of the four hazard metrics M i ( i = 1, … , 4 ) were computed along the 5-m isobath and sorted out into hazard classes C i k (Table 6). Each Table 8 Statistics of simulation results computed at J = 18, 201 save points along the 5-m isobath (excluding large bays) for PRT sources (Table 3; Fig. 15): mean, standard deviation, root-mean-square and average of the top 33, 10, 1 and 0. 1 3 metric was then plotted along the isobath, with the data points color coded as a function of their corresponding hazard class, i.e., in yellow, green, blue, magenta, and red, for k = 1, … , 5 . Figures 12 and 13 show examples of such results for the M9-HSP1 ACZ/LSB (Table 1; Fig. 2i) and the M9.0-PRT3 PRT/Caribbean arc (Table 3; Fig. 4h) sources, for the three, physical, metrics computed in grids G1,G2,G3 along the 5-m isobath (here, at J = 28, 401 save points including large bays). In each figure, panel (a) shows the envelope of maximum surface elevation above MWL, and panels (b,c,d) show the color-coded maximum elevation, current velocity magnitude, and momentum force, respectively (note, the class color code of surface elevation is also marked in (a) for points along the coast). While both of these M9 sources cause similar ranges of metric values along the coast, the absolute magnitude of each metric is affected by the initial tsunami directionality resulting from the source location and orientation. Specifically, as seen in earlier large scale results, the eastward ACZ/LSB source with a 15 o strike (M9.0-HSP1) affects more the lower USEC and the Caribbean Islands in Fig. 12, and the southward PRT source (M9.0-PRT3) affects more the upper USEC in Fig. 13. Additionally, as observed for the envelopes of maximum surface elevations (Figs. 7 and 9; see also Figs. 12a and 13a), for both sources, shelf refraction causes repetitive patterns of tsunami hazard modulation along the coast. Because these two sources are among the largest magnitude ones modeled here, many locations have hazard metrics in the third and fourth classes and a few in the fifth, highest hazard, class. Results would look qualitatively similar for the other, smaller magnitude, ACZ/LSB and PRT sources and will not be detailed individually. However, their overall impact on the USEC is compared with each other next. Figures 14 and 15 compare the four tsunami hazard metrics computed for the ten ACZ/ LSB sources (Table 1) and eight PRT/Caribbean arc sources (Table 3), respectively (at J = 18, 201 save points, here excluding large bays for more clarity). Overall, due to source location and directionality combined with deep water refraction, as noted before, the physical hazard metrics of the ACZ/LSB sources are larger on the southern USEC (south of 35 deg. N in Fig. 14), whereas it is the opposite for the PRT/Caribbean arc sources (Fig. 15). Additionally, as expected from the shelf refractions effects discussed before, the physical hazard metrics of the sources originating from the same area all have identical fine scale patterns of focusing/defocusing along the USEC, with overall each metric being larger, the larger the source magnitude. Finally, as expected, the tsunami travel times are very similar along the coast for all the sources originating from each area, except in the ACZ area for the two sources located further east in the HSP that have longer travel times by about 30 min. (Fig. 14d). Travel times for the ACZ/LSB sources in Fig. 14d range between 7.75 and 12.5 h, except for a few larger spurious values for the weakest ACZ sources, and are mostly in hazard class 2, with a smaller fraction of values in hazard class 1. For the PRT sources, Fig. 15d shows that travel times are shorter than for ACZ sources (2.5-7.25 h), but follow a similar pattern along the coast, due to the similar refraction over the wide shelf. The PRT travel times are nearly the same for all sources (within a few minutes of each other) and fall mostly within the medium and high hazard classes (3 and 4). Tables 7 and 8 provide statistics computed based on results plotted in Figs. 14 and 15, for the: mean, standard deviation, root-mean-square and average of the top 33, 10, 1 and 0.1 percentiles, of results for the first three physical metrics; these definitions are consistent with standard statistical analyses of ocean wave heights and period. Results show the expected trend, with the intensity of each metric statistics increasing with the source magnitude, except for the 0.1 percentile, which only averages a small sample of 18 individual results and hence is more sensitive to a few outliers. One exception is the two M9.0-HSP1 and -HSP2 sources in the ACZ area, which due to the effect of the MTR, cause a slightly lower overall impact on the USEC than the two M9.0-MTR1 and -MTR2 sources which are sited west of the MTR.
The hazard metric results for the 18 coseismic sources affecting the USEC are further discussed in the next section.

Detailed results of hazard metrics and classes along the USEC
For the ACZ/LSB sources, results of the three physical metrics in Fig. 14a-c and in Table 7 show that the two M8.0 sources, which have 32 times less energy than the M9.0 sources (see Eq. 1), only cause a ∼ 0.03 m average maximum elevation along the USEC, with a 0.03 or 0.013 m standard deviation, and the two M8.3 sources, which have 11 times less energy than the M9.0 sources, a ∼ 0.09 m average maximum elevation with a 0.07 standard deviation. Average maximum currents velocity magnitude and forces are commensurate with these values and small for the most part, with currents velocities being less than two knots (1 m/s) at all but a few locations (in 0.1 percentile). Note that larger current velocity (and corresponding force) values can result here from the coarse grid resolution and insufficient energy dissipation by bottom friction nearshore. Overall, statistics show that the physical metrics for the M8 and M8.3 sources cause a low to medium-low hazard. For elevation, only 0.1% of elevation values are in a hazard class higher than 1; 0.1% of current velocity values are in a hazard class higher than 2, and 99.9% of the force values are in hazard class 1. For the larger ACZ sources, results in Figs. 14a-c and in Table 7 show that, as expected, the largest hazard for the first three metrics is caused by the four M9.0 sources, with many values falling in the high/highest hazard classes 4/5 for each metric (Table 6). Among those, at many locations, slightly larger values of the three metrics are caused by the sources located west of the Madeira Torre Rise (MTR); but there are also locations where the sources located in the Horseshoe Plain (HSP) cause larger values. Looking at the detailed impact in the figures, south of 35.5 Lat. N. (approximately Cape Hatteras,NC), there are only small differences for the three metrics, between the sources located west of the MTR (M9-MTR1, M9-MTR2) and those east of it (M9-HSP1, M9-HSP2), with even some locations in the south where current velocities and forces are slightly larger for the latter two sources, indicating that here the Madeira Torre Rise (MTR) either does not affect tsunami propagation or reinforces it slightly. In contrast, north of 35 Lat. N., the tsunami hazard for the three metrics is larger for the M9 sources located west of the MTR than for the HSP sources, indicating that the shallow MTR ridge offers some protection from coseismic sources located in the area of the Lisbon 1755 event. Here, while maximum surface elevations above MWL reach similar maximum values as in the south, mostly in the high hazard category, current velocities and forces reach values in the highest hazard category at a larger number of locations.
Regarding the medium-high magnitude M8.7 LSB/ACZ sources, results in Fig. 14a-c and in Table 7 for the three metrics show that overall, these sources, which are about 3 times less energetic than the M9 sources (see Eq. 1), cause less impact at most locations than the M9 sources, while following the same pattern of highs and lows along the coast as for the M9 sources. However, at some isolated locations in the south, maximum elevation reaches up to 1.9 m for the M8.7 sources and in the north, particularly above 39 Lat. N (about Atlantic City, NJ), at a few locations current velocities and forces reach values nearly as large as those of the M9-HSP1 and M9-HSP2 sources. Clearly, this results from the complex tsunami propagation from the source, including wave interactions with the MTR near the source, dispersive effects during their transoceanic propagation, and interaction with the shelf slope and shelf bathymetry and coastal features along the USEC. Resonances (seiching) can also be triggered nearshore due to different frequency content in the incident wave train, that may increase current velocities even if tsunami elevations are not significantly affected (see, e.g., Nemati et al. 2019, for an analysis of such effects in their particular location).
For the PRT/Caribbean Arc sources, Fig. 15a-c show that as a result of refraction and a nearly northward directivity for each source, despite their location in different fault segments, the physical hazard metrics all exhibit the same patterns of focusing/defocusing onto the same coastal areas, with just a modulation of the intensity commensurate with the source magnitude. Overall, for these northward sources, the intensity of the three physical metrics is 2-3 times smaller south of 34 Lat. N. (approximately south of Cape Hatteras, NC) than north of it. So clearly, the PRT/Caribbean arc sources mostly impact the part of the USEC less impacted by the ACZ/LSB sources.
Statistics of the three physical metrics for the PRT/Caribbean arc sources (Fig. 15a-c; Table 8) show that the three M8.3 sources, which have about 21 times less energy than the M9.0 sources (see Eq. 1), only cause a ∼ 0.05 − 0.15 m average maximum elevation, with a 0.03-0.011 m standard deviation, and 0.4-0.7 m for the 0.1 percentile values. The two M8.7 sources, which have 11 times less energy than the M9.0 sources, have a ∼ 0.22 − 0.38 m average maximum elevation with a 0.15-0.23 m standard deviation, and 0.8-0.9 m for the 0.1 percentile values. Average maximum current velocity magnitudes and forces for these sources are commensurate with these inundation values and small for the most part, with current velocity being less than 1.25 m/s (2.5 knots) at all but a few locations (in 0.1 percentile). Note that larger current velocity (and corresponding force) values can again result from the coarse grid resolution and insufficient energy dissipation by bottom friction nearshore. Overall, statistics show that the physical metrics for the M8.3 and M8.7 sources cause a low to medium-low hazard. Only the 0.1% of maximum current velocity values fall in hazard class 3. As expected, tsunami hazard is most severe for the three larger M9 PRT sources. Despite resulting from quite different fault areas, results in Fig. 15a-c and in Table 8 show fairly similar overall impacts for the three M9 sources, with moderate  (Table 1; Fig. 2), based on results of the four hazard metrics in Fig. 14: a, [0.40, 0.30, 0.15, 0.15], and color-coded hazard intensity classes are: (yellow) low; (green) medium low; (blue) medium (magenta) medium high, (red) high, and (black) extreme. Note, for the last 2 sources, the TII is also provided within the large bays 1 3 average values of the three metrics. However, there is large alongshore variation of the metrics, as confirmed by values in the two highest percentile (1 and 0.1 %), that all fall in hazard classes 4/5, i.e., high to highest. Among the three sources, M9.0-PRT3, which was used as an extreme PMT source in earlier work, causes the maximum inundation and forces, but M9.0-PRT2 causes slightly larger maximum current velocities for the 10, 1 and 0.1 percentiles, which are all in the highest hazard category.

Tsunami intensity index (TII) values along the USEC
To provide an overall indicator of the intensity of tsunami hazard along the USEC, resulting from each coseismic source, values of the four coastal hazard metrics sorted into five classes are combined into a novel Tsunami Intensity Index (TII). The TII is calculated as a weighted average of the class index values k I (x j , y j ) , interpolated onto a 0 to 5 scale (from low hazard to highest hazard) for the four metrics M i (i = 1, … , 4) ( Table 6), at each save point location (x j , y j ) ( j = 1, … , J ), i.e., where w i are weights associated with each metric, with ∑ i w i = 1 , and the k i s are computed by first finding the hazard class within which each metric falls and then applying Eq. (3). The TII, hence, also has a 0 to 5 range and represents the overall tsunami hazard intensity, from low to highest, at any given save point (x j , y j ) . If all the metrics are given equal weight then we simply have w i = 0.25 ; however, as we shall see, the weights allow giving more importance to some metrics than others.
In the following, we arbitrarily give more weight to the metrics that directly affect the ability of the pedestrian population to evacuate, which makes the TII reflect more the hazard related to injuries/fatalities than structural damage, i.e., w i = [0.4, 0.3, 0.15, 0.15] for each metric, respectively. Specifically, in this example, maximum inundation height above MWL and current velocity magnitude metrics have a larger weight and, hence, forces on structures and tsunami travel time have relatively less weight. The rationale for giving less weight to travel time here is that these are all far-field coseismic tsunamis for which advanced warning would be given within minutes (at least in the US by the NOAA-operated tsunami warning centers) upon the triggering of an earthquake with parameters recognized to be tsunamigenic (e.g., magnitude, depth,...); hence, this warning would somewhat mitigate this aspect of the hazard. The hazard that cannot really be avoided is that caused by the tsunami flow itself in the inundated area, in terms of surface elevations and currents, as both of these would still affect the evacuation and survival of people who would not have evacuated soon enough, or would have procrastinated. Additionally, flow depth and current would cause structural damage and, when occurring in combination, create large momentum forces that would intensify this damage. Since the momentum force is somewhat correlated with maximum elevation and current, we gave it a lower weight as well to reduce redundancy. However, for the same simulation results, any other combination of metric weights could be applied and results for the TII easily generated. Figures 16 and 17 show the alongshore distributions of TII values computed with these weights using Eq. (3), based on the four hazard metric results computed for the ten LSB/ ACZ (Fig. 14) and eight PRT/Caribbean arc (Fig. 15) sources, respectively. We see that the TII class values (which are color coded here for different levels of hazard using the same color scheme as that of the hazard intensity classes) well summarize the overall tsunami hazard and confirm the detailed observations made above for the individual hazard metrics results. The TII, however, as it can be more easily understood than the raw hazard metrics, provides a unique and simple indicator of the hazard intensity expected from each source at the many save points along the coast that can help rapidly convey to emergency managers whether they need to worry about this event and hazard or not. Emergency managers can also more easily find the TII values for the locations of their critical facilities, rather than looking at individual metrics and making judgment calls. Although these figures are plotted at large scale, all results produced can be easily zoomed into and/or imported into a GIS, to evaluate the hazard faced by any particular area or specific community.

Estimating the effect of return period
Although a complete PTHA of tsunami hazard along the USEC is not the object of this study and will be addressed in future work, the estimated return periods provided in Tables 1 and 3, for each of the 18 coseismic sources whose hazard was modeled along the USEC, allow drawing some preliminary conclusions as of the expected maximum intensity of tsunami hazard on the USEC as a whole, for different return periods. Specifically, the statistics of the three physical hazard metrics computed for each source (Tables 7 and 8)) can provide a measure of the intensity of the tsunami hazard as a function of the return period T r estimated for each event (Tables 1 and 3). Altogether, the sources simulated here cover an estimated range of return periods, T r ≃ [70 − 2000] years, with the low end being the M8 ACZ sources and the high end the extreme M9.0-PRT3 source. While detailed results of each tsunami metric can be extracted for any save point along the coast for each of the sources, to be conservative we use the largest values of the overall statistics of the three physical tsunami hazard metrics computed for the entire USEC, independent of location, i.e., the 0.1 percentiles of maximum elevation above MWL, current velocity, and momentum force to estimate the maximum values one would expect to occur on the USEC as a function of T r . These are listed in Table 9 and provide a single order-of-magnitude estimate of the maximum value expected for each metric as a function of T r .

Conclusions
In this study, we simulated for the first time the regional tsunami hazard caused on the USEC by a collection of (18) realistic coseismic sources in the NAOB, covering a range of magnitude M8-M9 and estimated return periods of 70-2000 years, from the two most intense seismic areas in the far-field of the USEC, the ACZ/LSB and PRT/Caribbean Arc areas. This work builds on 10 years of development of high-resolution tsunami inundation maps, in the context of the US-NTHMP, for the areas of the USEC facing the highest tsunami hazard from PMTs (from coseismic, landslide and volcanic collapse sources), without consideration of return periods, and was meant to establish a preliminary baseline for a comprehensive PTHA that would be performed in future work.
Tsunami coastal hazard was represented by four metrics, which due to the coarse, regional, resolution used here were computed from model results interpolated along the 5 m isobath: (1) maximum surface elevation above Mean Water Level (MWL); (2) maximum current velocity magnitude; (3) maximum momentum force; and (4) tsunami travel time. Overall, the first three, physical, metrics, which quantify flooding, navigation, and structural hazards, are larger the larger the source magnitude, and larger in the main direction of incident tsunami propagation, which depends on source location/orientation and deep water wave guiding effects (i.e., southwestward for the ACZ/LSB sources and northward for the PRT sources). Additionally, for all the sources, hazard intensity also has a strong  (Table 3; Fig. 4), based on results of the four hazard metrics in Fig. 15: a, [0.40, 0.30, 0.15, 0.15], and color-coded hazard intensity classes are: (yellow) low; (green) medium low; (blue) medium (magenta) medium high, (red) high, and (black) extreme. Note, for the last source, the TII is also provided within the large bays alongshore variation, with similar patterns of higher and lower values, due to bathymetric control causing similar wave refraction patterns of focusing/defocusing over the wide USEC shelf, for any incident tsunami. This partial "memory loss" effect due to refraction makes tsunami coastal hazard somewhat less dependent on the tsunami source location and parameters (except for magnitude) and, hence, mitigates the effects of the large uncertainty and many assumptions made, in particular when defining the source parameters, on the intensity of tsunami coastal hazard. The fourth metric (4), which quantifies evacuation hazard, differs mostly between sources from each area (ACZ and PRT), but less so among sources from the same area; its inverse is used as a measure of increased hazard associated with short warning/evacuation times.
These results were used to compute a novel Tsunami Intensity Index (TII), calculated as a weighted average of the four metrics classified into five hazard classes (low, medium low, medium, high, and highest hazard) that can be used to quickly assess the tsunami hazard intensity resulting from each event, at any given location along the USEC. The TII values would allow emergency managers to easily and quickly identify coastal areas and critical infrastructures facing greater hazard from coseismic tsunamis and to more easily compare the tsunami hazard caused by sources of various magnitude and location. This information could then be used for establishing future coastal development and zoning plans.
Pending a comprehensive PTHA, the results presented here (e.g., in Tables 1,3 , 7, 8, 9) are also meant to provide some preliminary answers to the emergency management community of the 14 states covering the USEC, from Florida to Maine, as of the intensity of tsunami hazard that would be expected from coseismic sources of different return period (and probability) in the NAOB, particularly in comparison with other natural hazards affecting the different areas. To this effect, results from this regional hazard study for some of the sources presented here were already included in a GIS by an emergency management agency in charge of the USEC, in which our modeled hazard metrics were mapped against various layers of societal development and critical facilities. The overall estimates made here of tsunami hazard intensity for several physical variables, as a function of the event return period, can help put in perspective the tsunami hazard faced by the USEC with that resulting from other more prominent natural disasters such as tropical cyclones. Here, for the highest hazard areas (0.1 percentile), the tsunami hazard in the 250-500 year return period range (see Table 9) might be commensurate with that posed by the 100-year category 3-5 tropical cyclones that strike the USEC with an increasing frequency, taking into account the larger current velocities and forces caused by tsunami waves. Indeed, while the storm surge caused by such hurricanes might reach 4-6-m along the USEC (e.g., Grilli et al. 2017;Torres et al. 2019), the resulting coastal inundation is typically not associated with flow velocity and momentum forces as large as those caused by long, fast moving, tsunami waves which, hence, could be more damaging even when associated with a smaller coastal inundation. These results show that while not being yet a complete PTHA, the present study constitutes a meaningful step in this direction.
Based on results of the present work, obtained for eighteen coseismic sources of various magnitude (and hence return period), and similar results for other types of sources to be modeled in future work (i.e., for nearshore submarine landslides, volcanic flank collapse, and meteotsunamis), a PTHA analysis could be conducted in future work, first at the regional coarse grid scale (as done in the present work). Then for areas identified in this regional PTHA as having the highest tsunami hazard for a given return period (e.g., 100, 500 y, depending on how critical the considered facilities are), higher-resolution probabilistic inundation maps, similar to the current first generation NTHMP maps based on PMTs, could be developed. For the nested grid simulations required to develop such maps, archived results such as those of the present work could provide initial and boundary conditions for the coarser higher-resolution coastal grids.
Clearly, the simulations required to develop probabilistic high-resolution maps would be very computationally demanding since a comprehensive PTHA might require considering hundreds to thousands of individual sources of the various types (i.e., coseismic, landslide, volcanic, meteotsunami), and the question of whether to use a fully dispersive model in all simulations might need to be considered to reduce computational time. For instance, dispersive effects might only be considered for the grids simulating tsunami propagation over long (i.e., transoceanic) distances, where these are significant (e.g., Glimsdal et al. 2013), and for the highest resolution nearshore grids (i.e., order 5-10-m resolution) where undular bores would need to be resolved to accurately simulate the coastal tsunami hazard (e.g., Madsen et al. 2008;Schambach et al. 2019). As indicated before, typical computational times on NSF-XSEDE HPC resources, using 100-200 processors, were ∼ 2-10 h for each simulation with FUNWAVE presented here. So for a regional PTHA study on the USEC at a coarse grid level such as in this paper, if 1000 sources were considered, this could require 1000 s to 10,000 h of computer cluster time, using 100-200 processors, so up to 2 million CPU hours. If a multi-GPU cluster was used, with each GPU featuring at least 4000 processors, and using say 20 GPUs as is now available (or 80000 processors), a complete regional PTHA could be theoretically simulated for the USEC in a day or so, on current computer hardware.
at http:// github. com/ fengy anshi/ FUNWA VE-TVD/. The authors gratefully acknowledge the two anonymous reviewers for their in-depth and constructive reviews of our initial manuscript.
Author contributions SG led the work, acquired funding, defined methodology, processed results and wrote most of the paper. MM performed most model simulations, did background literature work, post-processed some results, made some figures and proofread paper. LS created model grids and performed model simulations for a few extreme sources, did background literature work, developed post-processing routines, and did some post processing of results. AG co-led the work, helped to acquire funding, defined some of the methodology, performed early tsunami simulations and developed initial post-processing routines, and proofread paper.
Funding This work was supported in the United States by NOAA's National Weather Service (NWS), through grants administered by the US National Tsunami Hazards Mitigation Program (NTHMP): NA-13-NWS4670014, NA-14-NWS4670041, NA-15-NWS4670029, NA-16-NWS4670034, NA-17-NWS4670010, NA-18-NWS4670073 and NA-19-NWS4670074. Numerical simulations reported in this work used HPC resources, as part of the Extreme Science and Engineering Discovery Environment (XSEDE) (project BCS-170015), which is supported by the National Science Foundation (NSF) grant number ACI-1548562.
Data availability Two Google Drive repositories have been created, with the input data that was used to set up and run the tsunami models and all the model result data that was used to create figures and tables in the paper for the ACZ/LSB and PRT sources. These can be accessed at: https:// drive. google. com/ drive/ folde rs/ 1NTTh BauiLr_ wnBdz w787U Ual6C 3pTexu? usp= shari ng and https:// drive. google. com/ drive/ folde rs/ 1q9XV rZfYs QgtEh nOs-zfyrm nbgzQ hOco? usp= shari ng, respectively.