impact on the US East Coast

Abstract. We perform numerical simulations of the coastal impact of large co-seismic tsunamis, initiated in the Puerto Rican trench, both in far-field areas along the upper US East coast (and other Caribbean islands), and in more detail in the near-field, along the Puerto Rico North Shore (PRNS). We first model a magnitude 9.1 extreme co-seismic source and then a smaller 8.7 magnitude source, which approximately correspond to 600 and 200 year return periods, respectively. In both cases, tsunami generation and propagation (both near- and far-field) are first performed in a coarse 2′ basin scale grid, with ETOPO2 bathymetry, using a fully nonlinear and dispersive long wave tsunami model (FUNWAVE). Coastal runup and inundation are then simulated for two selected areas, using finer coastal nested grids. Thus, a 15″ (450 m) grid is used to calculate detailed far-field impact along the US East Coast, from New Jersey to Maine, and a 3″ (90 m) grid (for the finest resolution), encompassing the entire PRNS, is used to compute detailed near-field impact along the PRNS (runup and inundation). To perform coastal simulations in nested grids, accurate bathymetry/topography databases are constructed by combining ETOPO2 2′ data (in deep water) and USGS' or NOAA's 15″ or 3″ (in shallow water) data. In the far-field, runup caused by the extreme 9.1 source would be severe (over 10 m) for some nearby Caribbean islands, but would only reach up to 3 m along the selected section of the East coast. A sensitivity analysis to the bathymetric resolution (for a constant 3″ model grid) of runup along the PRNS, confirms the convergence of runup results for a topographic resolution 24″ or better, and thus stresses the importance of using sufficiently resolved bathymetric data, in order to accurately predict extreme runup values, particularly when bathymetric focusing is significant. Runup (10–22 m) and inundation are found to be very large at most locations for the extreme 9.1 source. Both simulated spatial inundation snapshots and time series indicate, the inundation would be particularly severe near and around the low-lying city of San Juan. For the 8.7 source, runup along the PRNS would be much less severe (3–6 m), but still significant, while inundation would only be significant near and around San Juan. This first-order tsunami hazard analysis stresses the importance of conducting more detailed and comprehensive studies, particularly of tsunami hazard along the PRNS, for a more complete and realistic selection of sources; such work is ongoing as part of a US funded (NTHMP) tsunami inundation mapping effort in Puerto Rico.

selected section of the East coast.A sensitivity analysis to the bathymetric resolution (for a constant 3 model grid) of runup along the PRNS, confirms the convergence of runup results for a topographic resolution 24 or better, and thus stresses the importance of using sufficiently resolved bathymetric data, in order to accurately predict extreme runup values, particularly when bathymetric focusing is significant.Runup (10-22 m) and inundation are found to be very large at most locations for the extreme 9.1 source.Both simulated spatial inundation snapshots and time series indicate, the inundation would be particularly severe near and around the low-lying city of San Juan.For the 8.7 source, runup along the PRNS would be much less severe (3-6 m), but still significant, while inundation would only be significant near and around San Juan.This first-order tsunami hazard analysis stresses the importance of conducting more detailed and comprehensive studies, particularly of tsunami hazard along the PRNS, for a more complete and realistic selection of sources; such work is ongoing as part of a US funded (NTHMP) tsunami inundation mapping effort in Puerto Rico.

Introduction
Tsunami hazard assessment is critical for coastal communities, emergency services, industry, and to develop regional contingency plans in response to catastrophic events.Large numbers of fatalities and widespread destruction were Published by Copernicus Publications on behalf of the European Geosciences Union.S. T. Grilli et al.: First-order hazard from co-seismic tsunami sources in the Puerto Rico Trench caused by recent catastrophic tsunamis that struck unprepared coastal populations without any warning: the 1998 Papua New Guinea landslide tsunami was responsible for over 2000 deaths (e.g., Tappin et al., 2008), and the 2004 Indian Ocean tsunami, for over 230 000 deaths and many missing in 9 countries (e.g., Grilli et al., 2007;Ioualalen et al., 2007).In the US, to improve tsunami preparedness and awareness, the "National Tsunami Hazard Mitigation Program" (NTHMP) has undertaken the development of tsunami inundation maps, both in high risk areas (such as Hawaii and the West coast; see http://nthmp.tsunami.gov/index.html) as well as other areas not traditionally thought to be exposed to high tsunami hazard (e.g., the East Coast).Such maps, which are based on modeling the "envelope" impact of all relevant tsunami sources on a specific coastal area, have been developed and released for many of the higher risk US coastlines, but are still in development or in preparation for lesser risk areas, such as the East Coast.
Although not usually identified as a high tsunami-risk area, historical analyses of tsunami events in the Caribbean Sea have catalogued 27 likely candidates (Lander et al., 2002;O'Loughlin et al., 2003;Caribbean Tsunami Hazard, 2006).One of the most deadly event among those, the 1918 Puerto Rico tsunami, caused by a 7.3 magnitude earthquake in the Mona Passage, caused major damage on the West Coast of Puerto Rico (up to 6 m runup) and 116 fatalities (Mercado et al., 1998).Besides such historical analyses, we were recently reminded by the catastrophic 7.0 magnitude earthquake that hit Haiti (on the island of Hispaniola just West of Puerto Rico, see Fig. 1) on 12 January 2010, heavily damaging Port-au-Prince and killing over 217 000 people in the process, that the entire area overlying the small Caribbean plate, which includes Puerto Rico and pushes its way eastward (at less than 25 mm a year) against the much larger (subducting) North American and South American plates, is prone to large, dangerous and potentially tsunamigenic earthquakes (see, Zahibo et al., 2001, their Fig. 1, for the geodynamic context of faults in the area; Jansma, 2008).While the 12 January 2010 earthquake, which was mostly land-based, only generated a small tsunami, a large ocean-based earthquake in the Puerto Rican Trench (PRT), as we shall see, could generate a significant tsunami that might have catastrophic effects in the near-field on the lower lying coastal areas of the Puerto Rico North Shore (PRNS; e.g., San Juan), as well as induce significant far-field effects on some distant shores, including the upper US East Coast.
In a recent independent investigation of tsunami hazard along the northeastern United States coastline, among other tsunami sources, we analyzed the impact of a large coseismic tsunami initiated in the PRT/subduction zone (Fig. 1; Pérignon, 2006), due to a 9.1 magnitude earthquake occurring over 600 km of the trench extension, which is nearly the entire East-West length of the trench (Knight, 2006; details are provided later).This source, which was aimed at representing the maximum (long return period) event that could be expected from the PRT region, was overall geologically sound, but perhaps not fully realistic in its details.In this paper, we first briefly revisit far-field effects on the upper US East Coast due to a similarly large co-seismic tsunami initiated in the PRT, and then perform a detailed study of near-field tsunami hazard, in terms of runup and inundation, along the directly exposed PRNS.To do so, we first conduct a historical analysis of earthquake and tsunami events in the PRT area, to establish a basis for the likely magnitude and return period of co-seismic tsunami sources in the PRT.Based on this, we define and perform simulations of tsunami propagation and impact on the PRNS for two extreme coseismic tsunami sources: (i) the catastrophic 9.1 magnitude event used in our initial work (Pérignon, 2006); and (ii) a more realistic 8.7 magnitude event; return periods for these will be estimated at 600 and 200 years, respectively.
A variety of nested model grids, with 2 to 3 cell size, will be used to simulate both far-and near-field tsunami propagation.For the latter, we will create an accurate nearshore bathymetry by merging data sets available in the deep ocean (2 ETOPO-2 grid; Fig. 1) and nearshore (3 NGDC data; Divins and Metzger, 2008).Simulations of tsunami propagation and coastal impact will be performed using the fully nonlinear and dispersive Boussinesq model referred to as FUN-WAVE.This model, whose initial implementation was aimed at coastal waves (Wei et al., 1995), has later been parameterized for and successfully applied to a variety of tsunami case studies (both co-seismic and landslide; e.g., Watts et al., 2003;Days et al., 2005;Grilli et al., 2007;Ioualalen et al., 2007;Tappin et al., 2008;Karlsson et al., 2009; see Appendix A for additional details).FUNWAVE is a long wave model with full nonlinearity, similar to that of Nonlinear Shallow Water (NSW) equations models more traditionally used for studying tsunamis.However, FUNWAVE has more complete physics, since it also includes dispersive effects inherent to intermediate and deep water wave processes, which may occur in shorter waves that are part of tsunami trains.Such shorter waves may affect tsunami propagation and runup, through wave-wave interactions.When included in the equations of a long wave model, dispersive effects only manifest themselves if the application (e.g., landslide tsunamis) or the physics (e.g., shorter wave trains) call for it.If these do not, dispersive effects will remain negligible in the numerical results and FUNWAVE will essentially solve NSW equations.Finally, while Boussinesq models are more computationally intensive than NSW models, the efficient parallel implementation of FUNWAVE used here makes it possible simulating in a very reasonable time on a small computer cluster, nearshore tsunami impact in very large and finely resolved grids.We shall see, in time series of computed nearshore surface elevation on the PRNS, that nearshore tsunamis appear as a series of long waves, with shorter and more dispersive waves riding on top of those, similar to the Boussinesq simulations recently reported by Madsen et al. (2008).
FUNWAVE does not have an automatic mesh refinement or nesting implementation.Hence, to refine grid size nearshore, where (as we shall see) it becomes essential to properly resolve salient bathymetric features in order to correctly model runup (particularly in areas where bathymetric focusing is important), we use a manual nesting approach, similar to that successfully used by Grilli et al. (2007) and Ioualalen et al. (2007), to model observed runup in Thailand from the 2004 Indian Ocean tsunami.This nesting method consists in using overlapping coarser and finer grid domains, with the latter being large enough to allow for the meaningful part of the incoming tsunami wave train to enter it, while its front does not significantly interact yet with important coastal bathymetric features.This will be detailed later in the context of both far-field and near-field simulations in finer nested grids, for tsunamis initiated in the Puerto Rico Trench.
To assess the numerical accuracy and convergence of runup predictions, we will perform a sensitivity analysis of near-field tsunami impact to the bathymetric resolution, while maintaining the model grid constant.While tsunami hazard will be estimated and discussed along the entire PRNS, we will particularly focus on the lower lying and populated city of San Juan, which would be particularly vulnerable to tsunami impact.
Because we will only consider 2 extreme potential tsunami sources in the PRT, and these are not based on a detailed geological analysis, we refer to our work as a first-order analysis of tsunami hazard on the PRNS.Both more geologically detailed and diverse tsunami sources should be used in the future, to further refine these first-order hazard assessments.Additionally, while both our digital elevation maps (DEMs) and finer model grid have reasonably accurate 3 cells (or about 90 m), the inundation figures presented at the end of this paper only represent so-called Category 2 or secondgeneration inundation maps (according to NTHMP mapping and modeling guidelines that are about to be released).More accurate maps should be developed in the future, for both better estimating inundation hazard and drawing evacuation plans, which should be Category 3, i.e., based on 1/3 or so (or 10 m) accurate DEMs and inundation model grids.It should be noted that, as part of NTHMP supported work, both a more comprehensive analysis of tsunami sources affecting Puerto Rico, and the development of more detailed inundation maps are ongoing (A.Mercado, personal communication, 2010; see http://poseidon.uprm.edu/).tonic setting for Puerto Rico, with the island lying within the East-West trending plate boundary zone, between the WSW moving North American Plate (to the North and right on sub-figure) and the ENE moving Caribbean Plate (left on the sub-figure; Mercado and McCann, 1998; see also Zahibo and Pelinovsky, 2001 for the general geodynamic context of the region).Overall, the North American Plate subducts under the Caribbean plate, at a rate that has been estimated for the general Caribbean plate area from about 20 mm per year (DeMets, 1993) to about 37 mm per year (Sykes et al., 1982).
As in other recent studies near the PRT area (Zahibo and Pelinovsky, 2001;USGS, 2001;tenBrink and Lian, 2004;tenBrink, 2005;Jansma, 2008), we assume in the following a predominantly (lateral) strike-slip motion of the Caribbean plate at 20 mm per year with respect to the North American Plate, in the ENE direction, at a 10-20 degree angle with respect to the trench axis.
Frequent earthquakes occur in the region, as a result of the large component of relative strike-slip motion of the plates (see, Zahibo et al., 2003, their Fig. 1, andtenBrink, 2005, for historical maps of seismicity in the larger Caribbean Sea area).tenBrink (2005) mapped the depth and intensity of the large number of observed earthquakes of magnitude greater than 2.5 in the PRT region; earthquake locations are clearly aligned with the boundary of the subducting plates.By contrast, the same analysis only identifies 6 large www.nat-hazards-earth-syst-sci.net/10/2109/2010/Nat.Hazards Earth Syst.Sci., 10, 2109-2125, 2010 historical events of magnitude M w =7 or greater (a typical threshold for potentially large tsunamigenesis), for the past 220 years in or near the PRT; these are listed in Table 1.Among these, two events occurred with (estimated) magnitude greater than M w =8, and four of these events reportedly generated a tsunami, with three of those causing 5-7 m runup on Puerto Rico (cases 3-5 in Table 1; USGS, 2001;Zahibo et al., 2001Zahibo et al., , 2003;;Lander et al., 2002).For completeness, it was also reported by Dawicki (2005) that twelve earthquakes of magnitude M w =7 or greater occurred near Puerto Rico in the past 500 years, but no additional tsunami records, other than those presented in Table 1, were given.
Although a rigorous analysis of earthquake and tsunami return periods would be difficult to perform, due to the paucity of observations of large seismic events and tsunamis in the region (only one historical event is specifically located in the PRT), it appears from data in Table 1 that there were 3 large tsunamigenic events affecting Puerto Rico during an 80 year period and 5 large earthquakes during a 160 year period in the Puerto Rico area, two of those with magnitude greater than 8. Hence, as a first approximation, one can associate a magnitude M w =7.5-8.1 seismic event in the area around Puerto Rico with a 30 to 80 year return period.Similarly, Dawicki's (2005) data would yield an average 42 year return period, for events of magnitude M w =7 or greater, which is consistent with the latter data.Longer period events have not been observed, but based on estimated plates' subduction rate and approximate maximum length and width of the PRT area that could move during a future large scale event, one could try and estimate the magnitude of extreme seismic events in the trench, as a function of their return period.This is discussed below.
To prepare for future major tsunamis in the Caribbean region, Knight (2006) developed a first-order estimate of the most extreme earthquake that could occur in the PRT, and assessed the resulting potential tsunami hazard (mostly for the Caribbean islands).Specifically, Knight assumed a simple homogenous source model (i.e., without considering effects of shores, small islands and archipelagos), covering a 600 km by 150 km area of the Puerto Rico trench (i.e., nearly the full E-W extension of the PRT by three times its actual width), with a fault plane orientation based on the PRT geology (angles are given in a following section).This extreme source corresponded to a magnitude 9.1 earthquake, which using Okada's (1985) method yields an average slip of ¯ =11.9 m (and corresponding maximum slip of =19 m; see details of our slightly modified version of the method in Appendix A).Based on the estimated 20 mm per year subduction rate in the region, this average slip would yield a long return period earthquake of 600 year or so.Now, if the two largest historical events of magnitude 8.1 listed in Table 1 had affected the same (entire) area of the PRT, one could estimate their return period by prorating average slip to the released energy, as compared to the 9.1 event (note, under Okada's, 1985, method assumptions, total energy released by an earthquake is proportional to the assumed surface area and average slip; see Appendix A for details).This can be done using Hanks and Kanamori's relationship between energy M o [J] released by an earthquake and its magnitude: M w =(log M o /1.51) -6 (where 1.51=log 32.36), which yields an average slip: ¯ 8.1 = ¯ 9.1 /32.36 (9.1−8.1)=0.37 m.Based on the estimated subduction rate in the PRT, this would only represent a 20 year or so return event, while historical data points to a longer 30-80 year return period.This implies, as could be expected, that in such smaller but still significant events only a fraction of the length of the PRT was likely mobilized by the earthquake.For instance, for the same M w =8.1 event, due to the proportionality of released energy to slip and surface area, reducing the affected length of PRT to 300 km would increase average slip to 0.74 m and the estimated return period to 40 years or so.Note, this would also have the effect of proportionally increasing the initial tsunami source and further concentrating its effects on lands and islands closest to the earthquake area, around the PRT.However, since no or few observations were made for these historical events, it is nearly impossible to further constrain the tsunami source based on hydrodynamic observations (as, e.g., was successfully done for the widely observed 2004 Indian Ocean Tsunami; e.g., Grilli et al., 2007;Ioualalen et al., 2007).
While earthquakes cannot be predicted, tenBrink and Lian's (2004) recent survey of the PRT uncovered evidence of current seismic activity and internal stress build-up in the subduction (slip) zone near the PRT, which supports the "impending" occurrence of an earthquake in the trench, and justifies the urgency for estimating tsunami hazard in the region.This will first be done in the present work for an extreme, long return period events, by using the "600 year" 9.1 source developed by Knight (2006) (and already used in our preliminary work to estimate extreme tsunami hazard on the upper US East coast; Pérignon, 2006).Additionally, as it is desirable to estimate the likely maximum event that could occur in the PRT in the near future, we develop another source, based on accumulated potential slip in the trench since the last known event that significantly affected it in 1787, i.e., 223 years ago; this represents a ¯ =4.46 m average slip.Again, assuming first that the same 600×150 km entire area of the PRT would be affected by this earthquake (as for the 9.1 source) and applying the same proportionality methodology as above, yields the potential for a very significant M w =8.81 event occurring at present time ( ¯ 8.81 = ¯ 9.1 /32.36 (9.1−8.81)= 4.46 m).However, because it is unlikely that, for this smaller earthquake, all of the PRT would move at once (or in the same manner), and in the absence of detailed geological information, we arbitrarily reduce the total energy released by this event to about twothird of it, which corresponds to the lower bound correction of the return period range estimated for the 5 large historical events in the area (i.e., 20/30); this yields a M w =8.7 magnitude earthquake (8.7≈8.81 + log(2/3)).
Finally, to reduce the earthquake energy from an 8.81 to a 8.7 magnitude, for lack of geological information and to better allow for a direct comparison of tsunami simulations with the 9.1 event, rather than decreasing the affected area to 2/3 (e.g., 400×150 km) using the same value of average slip, we further assume that this smaller event would still affect the same entire PRT area (i.e., 600×150 km), but with a proportionally reduced slip (which conserves the total released energy).This yields a M w =8.7 magnitude earthquake with average slip ¯ 8.7 = ¯ 9.1 /32.36 (9.1−8.7)= (2/3) • 4.46=2.96m, which we still assume to have a return period of 200+ years or so (similarly, maximum slip for this event is found as 8.7 =4.72 m).In summary, in this work, we assess first-order coastal tsunami hazard both in the far field along the upper US East coast and in the near field on the PRNS, by simulating effects of large tsunamigenic earthquakes in the PRT.For lack of better information (both historical and geological), we first perform tsunami propagation and coastal impact simulations, both in the far field and along the PRNS, using the extreme co-seismic tsunami source of magnitude M w =9.1 developed by Knight (2006) (with average slip ¯ 9.1 =11.9 m and return period about 600 years).Then, we perform simulations for the "more reasonable" but still very significant M w =8.7 magnitude earthquake (with average slip ¯ 8.7 =2.96 m and return period 200+ years), affecting the entire PRT, representing a potential event in the near future.
3 Co-seismic tsunami simulations for sources in the PRT

Computation of co-seismic tsunami sources in the PRT
Following the standard practice for coseismic-tsunami sources, we specify the initial surface elevation as a "hot start" in the propagation model (i.e., maximum static elevation with no flow velocity; e.g., see, Grilli et al., 2007;Ioualalen et al., 2007).Thus, based on specified geometrical, material, geological and seismological parameters defining each seismic event, we first calculate the deformation of the ocean floor according to Okada's (1985) method.The latter represents the deformation of an elastic half space, resulting from a dislocation occurring along an oblique plane, representing the slip plane in between two subducting plates.Assuming water is incompressible and maximum seafloor deformation occurs over a few seconds (during which tsunami propagation is assumed to be negligible), the initial elevation of the water surface is set equal to the seafloor deformation in the earthquake source area.Parameters of Okada's method are three angles orienting the "slip plane" (strike ϕ, dip δ, rake λ), the length l and width w of the rupture area A = lw, and material (Lamé) parameters (µ, λ).[These parameters are equal in an isotropic elastic medium with Poisson ratio, ν = µ/(µ + λ)=0.5.]Using Okada's method, we calculate both the maximum and mean fault slips from the specified earthquake moment magnitude M w and values of other parameters listed above, and then calculate the vertical displacement of the ocean floor.Note, some changes were made to the original Okada (1985) method (see details in Appendix A).In particular, to more realistically represent the unknown, but certainly non-constant, natural slip distribution within the dislocation plane, a Gaussian-like distribution of slip is assumed, with mean value identical to that obtained from the moment magnitude (i.e., the same total energy release).
As detailed in the previous section, based on historical records, geodynamics, and other analyses, we elected to simulate 2 co-seismic tsunami sources in the PRT.Both of these sources, of magnitudes 9.1 and 8.7, respectively, correspond to en masse slip motion of the trench, over an area A, l=600 km long and w=150 km wide.Geological parameter values for the selected M w ≈ 9.1 event, with epicenter located at 19.5 • N lat.and 66 • W long., in local water depth d=7 km, with maximum slip 9.1 = 19 m, average slip ¯ 9.1 =11.9 m, and equivalent scalar Moment M o =4.5E+22 J, were obtained from Knight (2006) as: ϕ = 92 • , δ = 15 • , λ = 50 • with a shear modulus µ = 4.2×10 10 Pa. Figure 2a shows the resulting initial co-seismic tsunami elevation; we see, maximum seafloor uplift/elevation reaches nearly 6 m, while the maximum subsidence/trough near the source SE corner reaches over -2 m.For lack of better information, the 8.7 source is developed using the same geometric and geological parameters as for the 9.1 source, by prorating slip to the released energy as discussed above, i.e., 8.7 =4.72 m and ¯ 8.7 =2.96 m.The initial surface elevation for this smaller source is plotted in Fig. 2b, where we see that maximum elevation reaches nearly 1.5 m, which is about the same factor four reduction as for slip, as compared to the 9.1 source maximum elevation, and thus is consistent with Okada's formulation in Appendix A.

Far-field tsunami propagation and coastal impact
The 9.1 magnitude source of Fig. 2a is first used to initialize FUNWAVE and perform simulations of far-field tsunami propagation in a coarse 2 ×2 grid (about 3.5 km size mesh), covering the area and bathymetry shown in Fig. 1.Bathymetric data is obtained from the 2 ETOPO2 (gravitational anomaly) database.The computational grid has 1501×1351 cells, and time steps of about 5 s are used, based on a constant mesh Courant number and deep water tsunami celerity.Considering the E-W orientation of the source, parallel to the PRT, it can be anticipated that larger far-field waves will be generated northward, as well as in the NW direction, towards the Gulf of Mexico.This is confirmed by the envelope of maximum computed surface elevations shown in Fig. 3.In the W and NW directions, simulations predict significant far-field tsunami impact (up to 10 m runup) on the coasts of the Dominican Republic (Fig. 3b) and some other Caribbean islands (the Bahamas; Fig. 3a).No significant waves reach as far as the Gulf of Mexico, although up to 2 m runup occurs in Georgia (GA) and South Carolina (SC; Fig. 3a).In the North, Bermuda and further away Nova Scotia, which are directly north of the PRT, would suffer the greatest impact, with up to 8 m and 6 m runup or so, respectively (Fig. 3a, c).In Novia Scotia, large runups would be due in great part to bathymetric focusing in shallower water areas.Southwest of this region, large elevations shown in Fig. 3a correspond to the shallows of Georges Bank, which protect Maine (ME) and northern Massachussets (MA) by inducing strong breaking and dissipation of incoming tsunami waves, before they can reach the coast.Along the rest of the upper US East coast (New Jersey, NJ; Long Island, NY; Rhode Island and Cape Cod, MA), potential tsunami impact from this source would be smaller, although significant, with up to 3 m maximum runup in this coarse grid (more detailed simulations of tsunami impact in this area, which is marked by a yellow box in Fig. 1, are presented below).Runup is of course very severe on the PRNS, directly S of the source, with up to 10-15 m runup in this coarse grid (Fig. 3b), but this is the object of more detailed simulations discussed in the next section (this region is also marked by a yellow box in Fig. 1).Finally, directly East of Puerto Rico, some of the Antilles, mostly the Virgin Islands, Anguilla and Saint Kitts, would suffer from up to 8-10 m runup (Fig. 3a, b).
It should be stressed that the above runup predictions, which are made in a very coarse 2 grids, are only orders of magnitudes, that must be confirmed and refined based on more detailed simulations performed in finer nested regional grids, in areas identified to be at higher risk in the coarse grid simulations.
Since this work focuses on the near-field tsunami hazard along the PRNS, such accurate simulations of far-field coastal tsunami impact are only presented here for the upper East Coast, where the impact is the largest in the continental US.To do so, a 15 ×15 (about 450 m) regional grid is designed to encompass the coastlines of New Jersey (NJ), Long Island (NY), Connecticut (CT) Rhode Island (RI), and Cape Cod (MA) (yellow box in Figs. 1 and 4a).For this grid, depth is obtained by interpolating from a 100 m resolution USGS database of coastal bathymetry and topography.To initialize computations, the incoming tsunami wave trains simulated in the coarse 2' grid are re-interpolated in the finer grid (both surface elevation and horizontal flow velocity) at time t=200 min from the start of the event, which corresponds to the arrival of the first 4 major waves in the fine grid.Time step is adjusted according to the new mesh size (to about 0.25 s), to proceed with further simulations.
Figure 4a shows the envelope of maximum computed surface elevations in the fine grid, which confirms that coastal runup would be less than 3 m in most places, as a result of the 9.1 PRT source.Only southern Nantucket and Martha's Vineyard would experience up to 5 m runup in localized areas.Also note, bathymetric focusing of tsunami waves occurs in a number of locations, such as the left and right of the Hudson river Canyon (leading to large runup in NJ near Atlantic City and in Western Long Island), and a variety of other smaller canyons and ridges.
Figure 4b shows time series of tsunami elevation computed at 4 numerical gages (marked 1-4 in Fig. 4a), with depth 14.4, 3.2, 18.1, and 14.8 m, respectively.We see, the tsunami arrives as a series of 3-5 major waves, the first one having about a 15 period and the second one a 10 period at all locations.In parts of the records, much smaller waves with periods on the order of 1 can also be seen to ride on the longer waves.At point (1), we see waves (initial rises of water) reach the area of Atlantic City (NJ), 4 h 2 min after the start of the tsunami event.At points (2,3), waves reach Western and Eastern Long Island (NJ), after 4 h 15 min and 3 h 52 min, respectively.Finally, at point (4), waves reach Martha's Vineyard (MA), after 3 h 55 min.As a result of bathymetric focusing, incident waves are largest at points (2) and ( 4), with about a 2.6 and 2.2 m maximum height, respectively.A number of stripes of increased elevation can be seen in Fig. 4a.A closer inspection of animations of model results would show that leading incident waves are reflected off the coast, particularly in Long Island and NJ, which combine with other incident waves to cause some stripes of increased maximum elevation at various offshore locations seen in Fig. 4a.

Near-field tsunami propagation and impact on PRNS
Despite their coarse grid resolution and poor coastal bathymetry, results of basin-scale simulations in Fig. 3b showed, as could have been expected from the proximity to the 9.1 PRT source, that near field tsunami impact would be severe along the PRNS, with wave elevations reaching more than 10 m.Notably, two areas of increased coastal impact were predicted for the western and eastern extremities of the island's North shore.In the following, we confirm these predictions by performing more accurate and detailed simulations of near-field tsunami impact along the PRNS, for both the 9.1 and 8.7 co-seismic PRT sources.As was done for the upper East coast, to be able to have a high resolution mesh over a large enough area, while keeping the grid size (and computational time) reasonable, the finest resolution nearshore simulations are performed in a smaller nested regional grid encompassing the PRNS (specifically, for the 9.1 source using a 3 resolution grid (or 90 m) covering the area marked by a yellow box in Fig. 1).Various mesh resolutions will be tested in this grid, with the coastal bathymetry obtained in each case from a 3 accurate DEM database (see details below).
As before, computations in the nested grid are initialized using results of the coarse 2 basin-scale grid, at a time when most of the tsunami wave train has entered the grid area.However, because of the proximity of the PRNS to the PRT tsunami source (unlike the East Coast nested grid), whose effect extends all the way to shore as a small negative surface elevation (see Fig. 2; according to Okada's method), the selection of both the proper size and initialization time of the nested grid is not trivial.In fact, this was the object of many trials and errors, in view of the following requirements: (1) the nested grid must extend far enough North to allow for: (i) the main tsunami wave train to enter it at initialization time; and (ii) the front of the train not too significantly interact with the shoreline prior to this time.(2) Yet, because of its fine resolution, the nested grid domain must not be so large as to yield a prohibitive computational cost.This grid selection process is further detailed in a following section.
As we shall see, due to its fine resolution, simulations in the nested coastal grid will require using a much larger mesh than for the basin scale simulations and thus be more computationally demanding.Hence, these simulations are performed using the recently developed parallel version of FUNWAVE (Pophet, 2008;Pophet et al., 2010).Thanks to the Message Passing Interface (MPI) language, this version of the model achieves computational speed-ups nearly linear with the number of processors.

Bathymetry/topography
While adequate in the deep ocean, ETOPO-2 data is too coarse and not accurate enough in coastal areas to perform simulations of coastal tsunami propagation over fine bathymetric grids.To this effect, a high-resolution coastal bathymetric database is designed by merging the 2 resolution ETOPO-2 data shown in Fig. 2, North of 19 • N Lat., with NGDC's 3 resolution topography and bathymetry (Coastal Relief Model; Divins and Metzger, 2008), using a krigging interpolation.[Computerized digital images and associated databases are available from the National Geophysical Data Center, NOAA, US Department of Commerce, http://www.ngdc.noaa.gov/.]A relevant part of the resulting bathymetry and topography is shown in Fig. 5. Note, the over 7000 m deep PRT to the North as well as numerous coastal canyons and ridges, visible as a series of "undulations", for instance, in the 100 and 500 m depth contours.These are expected to induce focusing and defocusing of incoming tsunami waves, which must be correctly reproduced in the propagation simulations.A closer inspection of the contour lines indicates that wave focusing effects should be particularly significant on both the western and eastern extremities of the PRNS, where the coarse 2 simulations predicted maximum coastal impact.

Initial conditions and numerical parameters for near-field simulations
Near-field simulations were first performed in the nested coastal grid for the M w =9.1 source.Various grid resolutions were tested, the finest being 3 (about 90 m), and results (runup) sensitivity to both grid parameters and bathymetric resolution was analyzed.Based on these results (shown later), simulations were then performed for the M w =8.7 source using a 15 resolution, which was deemed adequate and made it possible using a single grid of reasonable size, encompassing both the PRT source area and the PRNS, without requiring nesting.
To initialize simulations in the nested coastal grid for the M w =9.1 source (shown in Fig. 2a), FUNWAVE is first run in the 2 coarse grid, which as before has 1501×1351 cells and covers the entire basin area shown in Fig. 1.Here, however, for consistency with near-field simulations, the more accurate bathymetry of Fig. 5 is used in the part of the coarse grid closer to the PRNS.Selection of the nested grid parameters is discussed below.For each simulation in the nested grid, both initial surface elevation and velocity are interpolated from coarse grid results over the finer grid, using bi-directional cubic splines; simulations are then Nat.Hazards Earth Syst.Sci., 10, 2109-2125, 2010 www.nat-hazards-earth-syst-sci.net/10/2109/2010/ S. T. Grilli et al.: First-order hazard from co-seismic tsunami sources in the Puerto Rico Trench 2117 re-started using a proportionally smaller time step, and performed at least up to the time of maximum runup on the PRNS.
A preliminary sensitivity analysis was performed to select both the best size and initialization time for the nested coastal grid, considering a number of constraints.To speedup simulations, this was done using a 15 resolution mesh (or about 450 m); once all parameters were selected, final computations in the nested grid were performed using a 3 resolution mesh (or about 90 m).The size of the nested grid domain was selected as follows, while trying to limit the number of meshes to a reasonable value.To allow for the deepest part of the PRT to be within the grid, the northern boundary was located at 19.50 • N Lat.Similarly, to include salient topographic features, both West and East of the PRNS, the southern boundary of the grid was located at 18.33 • N Lat.and its Western and Eastern boundaries at −67.4 and −65.4 • W Long., respectively.This corresponds to the area shown in Fig. 5 and also as a yellow box insert in Fig. 1.
Various stages of the incident tsunami train calculated in the coarse grid were successively used to initialize nested grid simulations, corresponding to computational times of 200 to 600 s (by steps of 50 s) after the start of the earthquake.In each case, it was first verified whether the initial tsunami seamlessly propagated from the coarse grid to the fine grid, by analyzing surface elevations computed at the first few time steps of simulations in the fine grid.In all cases, small initial free surface oscillations and adjustments only occurred at or near the incident northern boundary, due to truncating the tail of the tsunami train, and these were smaller, the larger the initialization time (corresponding to more of the main tsunami train entering the fine grid).For initialization times of 450 s or less, no significant initial changes were observed for the remainder (and most significant part) of the incident tsunami, which at this stage is made of very long small amplitude waves, most of these being located in very deep water where small scale bathymetric features do not matter.For longer initialization times, some adjustments occurred near the PRNS (see below).
A second verification was done in each case, at the front end of the tsunami train close to the PRNS, to check that incoming tsunami waves had not yet significantly interacted with fine scale bathymetric features or the shore in the coarse grid, before being moved to the finer grid.As a result of the topology of the co-seismic tsunami source (Fig. 2a), the initial tsunami specified in the nested grid typically exhibits a long and shallow negative elevation wave at the front, which is moving towards the PRNS; this is followed by both taller and somewhat shorter elevation and depression waves.For initialization times up to 450 s, both the negative elevation of the tsunami at the shore was found small, at the onset of simulations in the nested grid (order −1 to −2 m maximum; this will be seen later in results at numerical gages along the PRNS), and did not significantly change over the next few time steps.For longer initialization times, more significant interactions with the coast and shore had started occurring in the coarse grid, which precluded using those values.
Finally, a third verification was done regarding runup values.Thus, simulations were performed in the nested grid for all cases, up to obtaining maximum runup on the PRNS.Although all runup distributions exhibited the same overall variation along the PRNS, maximum runup increased with initialization time up to 350 s and then stabilized for longer initialization times.This is likely due to truncating too much of the main part of the incoming tsunami waves for the smaller initialization times, whereas for 350 s, most of these waves have entered the fine grid domain and hence can contribute to maximum runup.
Overall, it results from this sensitivity analysis and the various verifications conducted, that the initialization time of 350 s provides an adequate compromise between the various constraints at the back and front end of the tsunami train.This time was thus selected for further simulations, using an even finer resolution of 3 , which yields a 2401×1401 cell nested grid.[Note that all verifications listed above were also successfully conducted for this resolution as well.]Numerical parameters and results are given in the following for this finally selected nested grid simulation, using the M w =9.1 source.Simulations are first performed in the coarse grid with a: (i) time step of 2 s; (ii) total number of time steps of 5000; (iii) computational time of 9.65 h CPU using 6 processors of an 8 processor mini-cluster (speedup with 6 processors is 5.99 or a 99% efficiency of the parallel solution).Initialization from coarse grid results at 350 s yields the surface elevation shown in Fig. 6a, for the nested grid area closer to the PRNS.As expected from the above sensitivity analysis, at this time, the main tsunami waves are only approaching the shallower parts of the PRNS, without significantly interacting with it; initial elevations show small negative or positive values near the coast at most place.Model parameters for the nested grid simulations are a: (i) time step of 0.25 s; (ii) total number of time steps of 14400; (iii) computational time of 42.21 h CPU, using 16 processors of a 72 processor cluster (speedup with 16 processors is 11.56 or a 72% efficiency of the parallel solution).

Tsunami surface elevation and runup on PRNS, for M=9.1 source
Figure 6 shows a time sequence of tsunami elevation in the fine 3 grid (using the 3 base bathymetry of Fig. 5), for 24 min of simulations after the initialization time.These snapshots show the arrival (and withdrawal) of multiple large elevation and depression waves, with reinforced impact on each extremity of the PRNS.Significant inundation is also seen to occur for the eastern half of the PRNS, with particular impact in the area of San Juan (located in Fig. 7).
It should be noted that no reflection from the model lateral The sequence (a-j) corresponds to 0-24 min after the arrival of the tsunami onshore (350 s after the start of the event, at initialization of elevations in the small regional grid, from 2 basin results of Figs. 1  and 3).[Small linear waves seen incoming in (h-j) are spurious waves resulting from small grid initialization, which do not influence computed runups and inundations.]Fig. 7. Same simulation as in Fig. 6.Envelope of maximum surface elevation (in meters) calculated using FUNWAVE in a 3 model and bathymetric/topographic grid (Fig. 5).The thin black lines are isobaths marked in meters.The thick black line marks the mean coastline location (zero elevation) prior to tsunami arrival.boundaries, over which numerical sponge/absorbing layers are specified in FUNWAVE, can be seen on the picture sequence (as well as in longer animations of model results not shown here).The only artifact of the small grid initialization (due to wave train truncation) is a train of smaller and shorter spurious waves, which show up to the North in Fig. 6h-j late in the process, but do not affect runup or inundation results on the PRNS.This confirms the relevance of the model grid parameters and initialization selected for this nested simulation of coastal tsunami impact.
Figure 7 shows the envelope of maximum predicted surface elevations for these simulations.Maximum elevations (i.e., runup) are over 15 m, inland of the two westward and eastward locations previously identified as highest impact areas, in the 2 grid simulations and in the snapshots of Fig. 6.Around and east of San Juan, the maximum inundation significantly penetrates inland and reaches up to 10 m.
Figure 8 shows the runup variation calculated alongshore in the fine 3 grid, as compared to that of the 2 basin scale grid (using the base bathymetric map of Fig. 5).[Runup is defined here as the maximum water surface elevation occurring at the original mean sea level water line.]In the fine grid, Nat.Hazards Earth Syst.Sci., 10, 1-17, 2010 www.nat-hazards-earth-syst-sci.net/10/1/2010/The sequence (a-j) corresponds to 0-24 min after the arrival of the tsunami onshore (350 s after the start of the event, at initialization of elevations in the small regional grid, from 2 basin results of Figs. 1  and 3).[Small linear waves seen incoming in (h-j) are spurious waves resulting from small grid initialization, which do not influence computed runups and inundations.]Fig. 7. Same simulation as in Fig. 6.Envelope of maximum surface elevation (in meters) calculated using FUNWAVE in a 3 model and bathymetric/topographic grid (Fig. 5).The thin black lines are isobaths marked in meters.The thick black line marks the mean coastline location (zero elevation) prior to tsunami arrival.boundaries, over which numerical sponge/absorbing layers are specified in FUNWAVE, can be seen on the picture sequence (as well as in longer animations of model results not shown here).The only artifact of the small grid initialization (due to wave train truncation) is a train of smaller and shorter spurious waves, which show up to the North in Fig. 6h-j late in the process, but do not affect runup or inundation results on the PRNS.This confirms the relevance of the model grid parameters and initialization selected for this nested simulation of coastal tsunami impact.
Figure 7 shows the envelope of maximum predicted surface elevations for these simulations.Maximum elevations (i.e., runup) are over 15 m, inland of the two westward and eastward locations previously identified as highest impact areas, in the 2 grid simulations and in the snapshots of Fig. 6.Around and east of San Juan, the maximum inundation significantly penetrates inland and reaches up to 10 m.
Figure 8 shows the runup variation calculated alongshore in the fine 3 grid, as compared to that of the 2 basin scale grid (using the base bathymetric map of Fig. 5).[Runup is defined here as the maximum water surface elevation occurring at the original mean sea level water line.]In the fine grid, runup reaches 22 m and 18 m at the westward and eastward extremities of the PRNS, respectively, as a result of bathymetric focusing over nearshore bottom features discussed before.To the west of the island, the shape of bathymetric contours shown in Fig. 5 from 1000 to 100 m depth, is clearly such as to focus tsunami "wave rays" on the island extremity.To the east, the isobaths shown in Fig. 5, except for a few undulations, are fairly parallel and oriented W-E from 1000 to 100 m depth, while the coast is oriented ESE; hence, most of the tsunami wave refraction occurs in shallower water.In fact, animations of model results (not shown) clearly show that strong shallow water focusing occurs on the Eastern side of the island, as seen for instance in the snapshots of Fig. 6f-i.In the central part of the PRNS, by contrast, such strong focusing does not occur since, except for small undulations, isobaths are all more or less parallel to shore (Fig. 5); runup thus varies between 4-12 m, around an average of 7 m or so, mostly as a result of wave focusing and defocusing on the small scale undulations of the bathymetric contours.While runup in the middle part of the island is reasonably well predicted in the coarse 2 grid, the strong increase in runup towards the PRNS extremities is totally non-existent and runup is significantly under-predicted.This is likely because bathymetry is very under-resolved in the 2 grid and small-scale bottom features that act like lenses for tsunami waves have been smoothed out.This will be more apparent in the systematic sensitivity analysis to bathymetric resolution performed in the next section.

Sensitivity analysis of maximum surface elevation and runup on PRNS, to bathymetric resolution, for the M=9.1 source
To assess the convergence of computed runup and coastal inundation with bathymetric resolution, a sensitivity analysis is performed for the M w =9.1 source, by varying bathymetric resolution independently from that of the computational grid, which is kept at 3 to ensure an accurate numerical solution of FUNWAVE model equations.To do so, the 3 base map of Fig. 5 is downgraded by krigging re-interpolation, to a 6 , 12 , 24 or 48 bathymetric resolution; the resulting bathymetric contours are shown in Fig. 9.While the deep water bathymetry is little affected by this downgrading, it is clear that small scale features of the shallow water bathymetry and coastal topography gradually disappear (e.g., undulations in depth contours, shoreline tortuosity).As we shall see, this significantly affects tsunami propagation and coastal inundation and runup.The same initial condition as before, obtained from coarse grid simulations, is used in these four cases.Tsunami simulations are then performed in the nested grid, using FUNWAVE as detailed before.
The envelope of maximum surface elevation and runup are computed for each case and plotted in Figs. 9 and 10, respectively.In Fig. 9, we see that both maximum elevation and inundation (strength and extent of penetration) are signifi-S.T. Grilli et al.: please complete 11 runup reaches 22 m and 18 m at the westward and eastward extremities of the PRNS, respectively, as a result of bathymetric focusing over nearshore bottom features discussed before.To the west of the island, the shape of bathymetric contours shown in Fig. 5 from 1000 to 100 m depth, is clearly such as to focus tsunami "wave rays" on the island extremity.To the east, the isobaths shown in Fig. 5, except for a few undulations, are fairly parallel and oriented W-E from 1000 to 100 m depth, while the coast is oriented ESE; hence, most of the tsunami wave refraction occurs in shallower water.In fact, animations of model results (not shown) clearly show that strong shallow water focusing occurs on the Eastern side of the island, as seen for instance in the snapshots of Fig. 6f-i.In the central part of the PRNS, by contrast, such strong focusing does not occur since, except for small undulations, isobaths are all more or less parallel to shore (Fig. 5); runup thus varies between 4-12 m, around an average of 7 m or so, mostly as a result of wave focusing and defocusing on the small scale undulations of the bathymetric contours.
While runup in the middle part of the island is reasonably well predicted in the coarse 2 grid, the strong increase in runup towards the PRNS extremities is totally non-existent and runup is strongly under-predicted.This is likely because bathymetry is very under-resolved in the 2 grid and smallscale bottom features that act like lenses for tsunami waves have been smoothed out.This will be more apparent in the systematic sensitivity analysis to bathymetric resolution performed in the next section.

Sensitivity analysis of maximum surface elevation and runup on PRNS, to bathymetric resolution, for the M=9.1 source
To assess the convergence of computed runup and coastal inundation with bathymetric resolution, a sensitivity analysis is performed for the M w =9.1 source, by varying bathymetric resolution independently from that of the computational grid, which is kept at 3 to ensure an accurate numerical solution of FUNWAVE model equations.To do so, the 3 base map of Fig. 5 is downgraded by krigging re-interpolation, to a 6 , 12 , 24 or 48 bathymetric resolution; the resulting bathymetric contours are shown in Fig. 9.While the deep water bathymetry is little affected by this downgrading, it is clear that small scale features of the shallow water bathymetry and coastal topography gradually disappear (e.g., undulations in depth contours, shoreline tortuosity).As we shall see, this significantly affects tsunami propagation and coastal inundation and runup.The same initial condition as before, obtained from coarse grid simulations, is used in these four cases.Tsunami simulations are then performed in the nested grid, using FUNWAVE as detailed before.
The envelope of maximum surface elevation and runup are computed for each case and plotted in Figs. 9 and 10, respectively.In Fig. 9, we see that both maximum elevation and inundation (strength and extent of penetration) are signifi-  cantly affected by the resolution, the coarser the bathymetry (particularly coarser than 24 ), and more so at each extremity of the PRNS, where bathymetric focusing is most intense, as www.nat-hazards-earth-syst-sci.net/10/1/2010/Nat.Hazards Earth Syst.Sci., 10, 1-17, 2010  compared to results for the 3 resolved bathymetry (Fig. 7).This is even clearer in Fig. 10, which compares runup for each case.At both westward and eastward extremities, where the largest runups occur, results stay fairly consistent down to a 24 bathymetric resolution, but significantly degrade for a 48 resolution.For the middle section of the PRNS, however, from −65.8 to −67 • W Long., grid resolution effects seem less important at most locations.As indicated before, this is likely because of the fairly straight and shore-parallel bottom contours in this region.By contrast, the marked ridges and canyons at both PRNS extremities cause strong tsunami focusing and hence must be sufficiently resolved in order to accurately predict runup and inundation in model simulations.

Coastal inundation on PRNS for M=9.1 source
Coastal inundation for the 9.1 magnitude source was visible on Figs. 6 and 7.In Fig. 7, we see that the tsunami penetrates significantly inland at many locations, particularly along the flatter eastern half of the PRNS, where water elevation reaches up to 10 m or more.
Three computed time series of tsunami arrival on the PRNS at locations (numerical gages) marked (1-3) in Fig. 5, are shown in Fig. 11.These locations are respectively on the West, central (near San Juan) and East side of the PRNS and have depths of 51.2, 10.2, and 22.7 m, respectively.As discussed before, the elevation at each gage at initialization time in the fine grid (t=0 and 350 s from the start of the event) is small and negative (order of −1 to −2 m), representing the extent of interactions of the incoming tsunami with the shoreline prior to this time, i.e., essentially a small initial withdrawal, which is well captured in coarse grid simulations.
Similar to the time series in Fig. 4b for the upper East Coast, the tsunami arrives on the PRNS as a series of 4 to 5 major long waves seen in Fig. 11, the first having a 11-12 min period.Here as well, many shorter and steeper waves of period order 1 min or less ride on the longer tsunami oscillations, and the 3 grid allows such waves to be much better resolved than in the 15 regional East Coast grid (these are also most developed and intense for the shallowest gage 2).Such shorter wave trains result from the expression of dispersive effects (not included in standard NSW model equations, but inherent to FUNWAVE's), and are fully consistent with the recent work of Madsen et al. (2008) and their (and many other) reported observations of the arrival of the large 2004 Indian Ocean tsunami on many shorelines, as long waves causing inundations, with a series of overriding much shorter waves, more akin to typical gravity swells, which end-up arriving on the shore as a series of strongly breaking bores.In the time series of Fig. 11, the lowest and highest elevations are (−12, +8), (−5, +8) and (−14, +10) m for points 1-3, respectively, which is consistent with runup values reported in Fig. 10.
The large city of San Juan (around 66 • 09 W Long.; Fig. 12a), which is a fairly flat area in the region of maximum penetration of the tsunami inundation (Fig. 7), is selected for further analysis because of its dense population.For comparison, a similar analysis is performed in the maximum runup region on the west side of the PRNS, near the city and airport of Isabela, which are located on top a steep cliff (around 67 • 05 W Long.; Fig. 12b).In and near San Juan, both the simulated inundation and runup reach up to 10 m, and the tsunami penetrates a large distance inland (Figs. 7,10,and 11b).In the Isabela area, where runup reaches over 20 m, by contrast, due to the steep cliff-like shoreline, the inundation only penetrates a short distance inland (Figs. 7 and 10).
Figure 13 shows a more detailed series of inundation snapshots based on results computed in the fine grid, in and around San Juan, at 2-3 min intervals over a 20 min period following initialization, for the same case as in Figs. 6, 7, 10, and 11. Figure 13 displays simulated tsunami elevations above local ground level onshore, and below sea level offshore, as a color scale.Inundation simulations are laid over an aerial picture of the city, similar to that depicted in Fig. 12a, which helps visualizing the inundation penetration Nat.Hazards Earth Syst.Sci., 10, 1-17, 2010 www.nat-hazards-earth-syst-sci.net/10/1/2010/ compared to results for the 3 bathymetry (Fig. 7).This is even clearer in Fig. 10, which compares runup for each case.At both westward and eastward extremities, where the largest runups occur, results stay fairly consistent down to a 24 bathymetric resolution, but significantly degrade for a 48 resolution.For the middle section of the PRNS, however, from −65.8 to −67 • W Long., grid resolution effects seem less important at most locations.As indicated before, this is likely because of the fairly straight and shore-parallel bottom contours in this region.By contrast, the marked ridges and canyons at both PRNS extremities cause strong tsunami focusing and hence must be sufficiently resolved in order to accurately predict runup and inundation in model simulations.

Coastal inundation on PRNS for M=9.1 source
Coastal inundation for the 9.1 magnitude source was visible on Figs. 6 and 7.In Fig. 7, we see that the tsunami penetrates significantly inland at many locations, particularly along the flatter eastern half of the PRNS, where water elevation reaches up to 10 m or more.
Three computed time series of tsunami arrival on the PRNS at locations (numerical gages) marked (1-3) in Fig. 5, are shown in Fig. 11.These locations are respectively on the West, central (near San Juan) and East side of the PRNS and have depths of 51.2, 10.2, and 22.7 m, respectively.As discussed before, the elevation at each gage at initialization time in the fine grid (t=0 and 350 s from the start of the event) is small and negative (order of −1 to −2 m), representing the extent of interactions of the incoming tsunami with the shoreline prior to this time, i.e., essentially a small initial withdrawal, which is well captured in coarse grid simulations.
Similar to the time series in Fig. 4b for the upper East Coast, the tsunami arrives on the PRNS as a series of 4 to 5 major long waves seen in Fig. 11, the first having a 11-12 min period.Here as well, many shorter and steeper waves of period order 1 min or less ride on the longer tsunami oscillations, and the 3 grid allows such waves to be much better resolved than in the 15 regional East Coast grid (these are also most developed and intense for the shallowest gage 2).Such shorter wave trains result from the expression of dispersive effects (not included in standard NSW model equations, but inherent to FUNWAVE's), and are fully consistent with the recent work of Madsen et al. (2008) and their (and many other) reported observations of the arrival of the large 2004 Indian Ocean tsunami on many shorelines, as long waves causing inundations, with a series of overriding much shorter waves, more akin to typical gravity swells, which end-up arriving on the shore as a series of strongly breaking bores.In the time series of Fig. 11, the lowest and highest elevations are (−12, +8), (−5, +8) and (−14, +10) m for points 1-3, respectively, which is consistent with runup values reported in Fig. 10.
The large city of San Juan (around 66 • 09 W Long.; Fig. 12a), which is a fairly flat area in the region of maximum penetration of the tsunami inundation (Fig. 7), is selected for further analysis because of its dense population.For comparison, a similar analysis is performed in the maximum runup region on the west side of the PRNS, near the city and airport of Isabela, which are located on top a steep cliff (around 67 • 05 W Long.; Fig. 12b).In and near San Juan, both the simulated inundation and runup reach up to 10 m, and the tsunami penetrates a large distance inland (Figs. 7, 10, and 11b).In the Isabela area, where runup reaches over 20 m, by contrast, due to the steep cliff-like shoreline, the inundation only penetrates a short distance inland (Figs. 7 and 10).
Figure 13 shows a more detailed series of inundation snapshots based on results computed in the fine grid, in and around San Juan, at 2-3 min intervals over a 20 min period following initialization, for the same case as in Figs. 6, 7, 10, and 11. Figure 13 displays simulated tsunami elevations above local ground level onshore, and below sea level offshore, as a color scale.Inundation simulations are laid over an aerial picture of the city, similar to that depicted in Fig. 12a, which helps visualizing the inundation penetration in the city and identifying the most affected areas (within the 90 m computational grid accuracy).Specifically, Fig. 13a-b shows the arrival of a small leading depression wave (visible in Fig. 11b), followed in Fig. 13c-f by the main elevation wave; maximum inundation penetration and impact are depicted in Fig. 13e for the West side of the city and Fig. 13f for the East side.These correspond to 9 and 11 min, respectively, after the first arrival onshore of the tsunami (t=0 in Fig. 11b), itself occurring 350 s (about 6 min) after the beginning of the event and earthquake in the PRT. Figure 13 g shows the arrival of a strong depression wave (visible for t=14 min in Fig. 11b), with the simultaneous quick withdrawal of the inundation.The latter continues in Fig. 13h for the middle part of the city, while a new elevation wave arrives (visible for t=16 min in Fig. 11b) on both the West and East sides of the city, causing further inundation.The last two pictures show the arrival of another smaller depression wave (also seen in Fig. 11b), with the eventual withdrawal of the inundation.
Figure 14 shows a similar sequence results for the area around Isabela (Fig. 12b) where, as expected, for the same succession of depression and elevation waves as in Fig. 13 (albeit with a marked phase shift, as seen in Fig. 11a), little or no inland penetration of the large inundation occurs, due to the steep cliff-like features facing the ocean to the north.

Maximum elevation and coastal inundation on
PRNS for the M=8.7 source Similar computations are performed for the less catastrophic, but nevertheless still large, M w =8.7 magnitude PRT source, corresponding to a 200 year or so earthquake (Fig. 2b).As indicated before, in view of the smaller initial tsunami and based on the above sensitivity analyses to both the nested grid parameters and bathymetry for the 9.1 source, this simulation is performed a single larger 15 regional grid (about 450 m), using the same 3 bathymetric base map developed earlier, and encompassing both the PRNS and the source area in the PRT.Hence no grid nesting was required here.
As shown in Fig. 2, since only maximum slip is changed from the 9.1 to the 8.7 magnitude source, initial elevations are geometrically similar, but significantly reduced for the 8.7 source (by a factor of 4 or so).Accordingly, the simulated envelope of maximum elevation/coastal inundation, and runup, shown in Figs. 15 and 16, respectively, are also similarly reduced for the latter source.Only the low elevation regions West of San Juan would still be significantly flooded.Runup is reduced by a factor of 2-4, depending on location, reaching up to 5-6 m at each extremity of the PRNS and 2-3 m at other in between locations (hence bathymetric focusing still occurs at the extremities).Clearly, while still being quite damaging for the nearshore coastal area, such runup and related inundation would not reach the disaster level that a M w =9.1 source would cause, except perhaps for San Juan.

Conclusions
We performed a first-order tsunami hazard analysis of farfield effects in the NW Atlantic basin (in greater detail for the upper US East coast) and near-field effects on the Puerto Rico North shore, of large or extreme co-seismic tsunamis generated by earthquakes triggered in the Puerto Rico trench.We first used a magnitude 9.1 catastrophic earthquake source, which would correspond to a 600 year or so return period.Then, based on this source, we designed and used a smaller 8.7 magnitude source, which would correspond to a 200 year or so return period.In both cases, coastal runup and inundation were computed and discussed both in general terms along the PRNS and more specifically for two selected key agglomerations near regions of intense impact (San Juan and Isabela).
In the far field, tsunami waves and impact caused by the catastrophic 9.1 source would be most severe on nearby Caribbean islands West of Puerto Rico, the Antilles to the East, and the Bahamas and Bermudas.On the upper US East Coast, although significant, the tsunami would cause up to Fig. 15.Bathymetry (black isobaths in meter) and maximum surface elevation and (coastal) inundation (color scale in meters) along PRNS, simulated for the M=8.7 co-seismic tsunami source (using a 15 model grid, with the 3 topography and bathymetry depicted in Fig. 5).The thick black line marks the mean coastline location prior to tsunami arrival (zero elevation).
2-3 m runup at most places, similar to the storm surge caused by a large tropical cyclone, and tsunami travel time would be on the order of 4 h, hence allowing for sufficient warning.
In the near-field, along the PRNS, we identified for both sources two areas of elevated runup, on the western and eastern sides, as compared to runup in the middle part of the PRNS.This represents up to 18-22 m and 5-6 m runup, for each source, respectively.These elevated runups appear to be due to wave energy focusing, through refraction over nearshore bathymetric features.For other locations along the PRNS, where isobaths are fairly straight and parallel to shore, runup only reaches up to 10 m and 3 m, respectively.For the 9.1 source, a sensitivity analysis of runup to the bathymetric resolution (3 to 48 ), while maintaining the model grid resolution at 3 (or about 90 m), indicated result convergence provided the bathymetric resolution is at least 24 or better.This stresses the importance of having both accurate and sufficiently resolved bathymetric data for nearshore tsunami transformation simulations.
Maximum inundation elevation was calculated along the PRNS, for both sources.For the 9.1 source, the inundation significantly penetrates onshore for the eastern half of the coast, while the western side is somewhat protected from flooding by steep or cliff-like coastal features (despite the intensity of tsunami impact and runup).Inundation was found to be most intense in the lower-lying region surrounding the city of San Juan.Detailed inundation simulations near San Juan and Isabela, superimposed on aerial pictures, showed the time evolution of the water line as a function of the arrival of successive tsunami waves (also documented in time series at numerical gages).Such results provide a means of assessing tsunami impact on specific small scale coastal features or structures (for the latter, however, model grid resolution would need to be further increased, perhaps to 10-30 m or so, in an even finer nested inundation grid).For the much less energetic 8.7 source, inundation was only significant in the region West of San Juan.It should be stressed, as indicated before, that this work is only a first-order analysis based on fairly crude tsunami sources, which do not account for detailed geology in or near the PRT as well as the presence of nearby Caribbean islands.Also, considering the PRT is more or less parallel to the PRNS (Figs. 1 and 2) and has not ruptured for 223 years, one could design other potential co-seismic sources, in which the same earthquake energy (e.g., corresponding to a M=8.7 magnitude) would be concentrated over a shorter length of the trench (rather than the entire 600 km E-W extension of the PRT used here).This would locally yield proportionally larger slip and energy density, with the potential for stronger tsunami generation over a smaller length of the PRT facing a shorter section of the PRNS, and hence potentially larger runup in some locations of the PRNS.However, only detailed analyses of the geology and recent seismicity in the PRT, which are beyond the scope of this work, could help in selecting such more realistic earthquake scenarios and mechanisms.In any case, one can certainly conclude that coastal impact for the extreme 9.1 source provides an upper bound for the first-order tsunami hazard that could reasonably be expected along the PRNS, while the impact for the 8.7 source represents a lower bound for a large event that could be expected in the near future.
These preliminary first-order results stress the importance of conducting more detailed and comprehensive studies, particularly of tsunami hazard along the PRNS; such work is ongoing as part of a US funded (NTHMP) tsunami inundation mapping effort in Puerto Rico.

Figure 1
Figure 1 (sub-figure) shows the topography and bathymetry around the deep PRT, north of the island.The trench is approximately 770 km long and 50 km wide, with a depth reaching over 7000 m (up to 8340 m at one location).The northeastern portion of the Caribbean Plate is the general tec-

Fig. 1 .
Fig. 1.ETOPO-2 bathymetry (2 accuracy) for the NW Atlantic Ocean basin (axes are in degrees of N Lat.and W Long.).Both the island of Puerto Rico (bottom) and the upper US East Coast region (top), are marked by yellow rectangular boxes.The lower box is approximately 204 by 126 km.The sub-figure shows enhanced topography and bathymetry in and around Puerto Rico; the 770 km long, 50 km wide, and 7 km deep Puerto Rico Trench (PRT; pink color) is to the north of the island (USGS, 2001).

Fig. 2 .
Fig. 2. ETOPO-2 bathymetry (black isobaths are marked in meters) and tsunami source elevation (color scale in meters) for M w = (a) 9.1 and (b) 8.7, magnitude co-seismic sources, calculated using Okada's method as initial condition for FUNWAVE, on a 2 ×2 grid (axes are in degrees of N Lat.and W Long.).

Fig. 3 .
Fig. 3. FUNWAVE simulations over a 2 ×2 grid, using ETOPO-2 bathymetry, for the 9.1 co-seismic tsunami source of Fig. 2a (axes are in degrees of N Lat.and W Long.).(a-c) Maximum surface elevations (color scale in meter) at anytime during computations.(b) Zoom near and around Puerto Rico; (c) zoom near Nova Scotia.
Fig. 4. (a) Bathymetry (black isobaths marked in meters; USGS 100 m resolution) for New Jersey, Long Island, Rhode Island and Cape Cod areas, and nearshore tsunami simulations for 9.1 PRT coseismic source of Fig. 2a, in nested 15 regional grid (450 m; initialized with Fig. 3 results): (a) maximum surface elevation at any time during computations (color scale in meter; axes are in degrees of N Lat.and W Long.); (b) time series of tsunami arrival at locations 1-4 (numerical gages) marked as symbols (•) in (a), with depth 14.4, 3.2, 18.1, and 14.8 m, respectively , for points (1) to (4); times correspond to the start of the tsunami event.
S. T.Grilli et al.:  First-order hazard from co-seismic tsunami sources in the Puerto Rico Trench

Fig. 6 .
Fig.6.Snapshots of surface elevations (color bar scale in meter; dark blue below MWL; light blue above MWL; light yellow: emerged), simulated in a 3 model grid, for the M=9.1 PRT coseismic tsunami source (Fig.2a), using the 3 bathymetry of Fig.5.The sequence (a-j) corresponds to 0-24 min after the arrival of the tsunami onshore (350 s after the start of the event, at initialization of elevations in the small regional grid, from 2 basin results of Figs.1 and 3).[Small linear waves seen incoming in (h-j) are spurious waves resulting from small grid initialization, which do not influence computed runups and inundations.]

Fig. 6 .
Fig.6.Snapshots of surface elevations (color bar scale in meter; dark blue below MWL; light blue above MWL; light yellow: emerged), simulated in a 3 model grid, for the M=9.1 PRT coseismic tsunami source (Fig.2a), using the 3 bathymetry of Fig.5.The sequence (a-j) corresponds to 0-24 min after the arrival of the tsunami onshore (350 s after the start of the event, at initialization of elevations in the small regional grid, from 2 basin results of Figs.1 and 3).[Small linear waves seen incoming in (h-j) are spurious waves resulting from small grid initialization, which do not influence computed runups and inundations.]

Fig. 9 .
Fig. 9. Bathymetry (black isobaths marked in meters) and maximum tsunami surface elevation (color scale in meters), for the same source as in Figs. 6 and 7, using a 3 model grid with downgraded bathymetry obtained by re-interpolating the 3 s.data of Fig. 5 by krigging, on 6, 12, 24, and 48 s grids (top to bottom).

Fig. 10 .
Fig. 10.Maximum runup (in meter) along PRNS, for cases in Figs.7 and 9.For clarity, all solid lines denote smoothed curves based on raw computational data (not shown).

Fig. 9 .
Fig. 9. Bathymetry (black isobaths marked in meters) and maximum tsunami surface elevation (color scale in meters), for the same source as in Figs. 6 and 7, using a 3 model grid with downgraded bathymetry obtained by re-interpolating the 3 s.data of Fig. 5 by krigging, on 6 , 12 , 24 , and 48 grids (top to bottom).

Fig. 11 .
Fig. 11.Same simulations as in Figs. 6 7 (fine grid).Time series of surface elevation at locations 1-3 marked in Fig. 5, along PRNS.Depth at gages is: (a) 51.2;(b) 10.2; and (c) 22.7 m from MWL. [Time is with respect to fine grid simulations (i.e., starting at 350 s).] Fig. 12. Aerial photographs of PRNS areas selected for detailed inundation studies of Figs. 13 and 14, near the cities of: (a) San Juan; and (b) Isabela (Google Map, 2010).Note the many submarine canyons and ridges visible in figure (a), off of the middle section of the PRNS near San Juan, as seen in isobath undulations in Fig. 5 and discussed in the text.

Fig. 16 .
Fig. 16.Comparison of coastal runup along PRNS, for M=9.1 (.--.) and 8.7 (o---o) co-seismic tsunami sources in PRT.Computations are performed in 3 and 15 model grids, respectively, using the 3 resolution topography and bathymetry depicted in Fig. 5. Solid lines are smoothed representations of the symbols indicating model data points.