Modeled Flooding by Tsunamis and a Storm Versus Observed Extent of Coral Erratics on Anegada, British Virgin Islands—Further Evidence for a Great Caribbean Earthquake Six Centuries Ago

Models of near‐field tsunamis and an extreme hurricane provide further evidence for a great precolonial earthquake along the Puerto Rico Trench. The models are benchmarked to brain‐coral boulders and cobbles on Anegada, 125 km south of the trench. The models are screened by their success in flooding the mapped sites of these erratics, which were emplaced some six centuries ago. Among 25 tsunami scenarios, 19 have megathrust sources and the rest posit normal faulting on the outer rise. The modeled storm, the most extreme of 15 hurricanes of category 5, produces tsunami‐like bores from surf beat. In the tsunami scenarios, simulated flow depth is 1 m or more at all the clast sites, and 2 m or more at nearly all, given either a megathrust rupture 255 km long with 7.5 m of dip slip and M8.45, or an outer‐rise rupture 130 km long with 11.4 m of dip slip and M8.17. By contrast, many coral clasts lie beyond the reach of simulated flooding from the extreme hurricane. The tsunami screening may underestimate earthquake size by neglecting trees and shrubs that likely impeded both the simulated flows and the observed clasts; and it may overestimate earthquake size by leaving coastal sand barriers intact. The screening results broadly agree with those from previously published tsunami simulations. In either successful scenario, the average recurrence interval spans thousands of years, and flooding on the nearest Caribbean shores begins within a half‐hour.


Introduction
Disastrous Asian tsunamis in 2004 and 2011 have spurred reassessment of earthquake and tsunami potential at subduction zones worldwide.Subduction rate and subducting-plate age have been deemed unreliable as controls on maximum earthquake size (Stein & Okal, 2005); perhaps any subduction zone long enough to make a giant earthquake should be assumed capable of doing so (McCaffrey, 2008).The reassessments commonly include extending earthquake and tsunami history thousands of years into the geologic past.In some instances, geologic traces of past tsunamis have aided in evaluating numerical models of tsunami sources and inundation (MacInnes et al., 2010;Sawai, 2020).
Such assessment of megathrust hazards can be clouded by normal faults and tropical cyclones.Normal faulting on the outer trench wall and the adjoining outer rise, produced by bending of the incoming oceanic plate, not only produces higher stress drop and a higher ratio of radiated energy to seismic moment than in megathrust earthquakes (Ammon et al., 2008), but also generates tsunamis efficiently because the normal faults dip steeply and rupture to the surface beneath deep water (Craig et al., 2014;Satake & Tanioka, 1999).Examples include twentieth-century tsunamis in northeast Japan (Kanamori, 1971(Kanamori, , 1972;;Okal et al., 2016), Indonesia (Gusman et al., 2009;Lynnes & Lay, 1988), Mariana Islands (Satake et al., 1992), Kuril Islands (Ammon et al., 2008;Fujii & Satake, 2008;Lay et al., 2009;Steblov et al., 2008), and the northern Tonga Trench (Beavan et al., 2010;Fritz et al., 2011;Lay et al., 2010).In further examples, outer-rise faults were activated by the giant 2011 Tohoku earthquake (Lay et al., 2011) and have been identified as likely sources of tsunami hazard in Central America (Alvarez-Gomez et al., 2012).Outer-rise earthquakes are typically located on the outer trench wall trenchward of the bathymetric outer rise, because of the concentration of bending stresses and in-plane stresses there (e.g., Craig et al., 2014).As for tropical cyclones, 2013 Typhoon Haiyan showed that storm waves, as they pile up against a coral reef, may spawn tsunami-like bores.Such bores came ashore in the Philippines during 2013 Typhoon Haiyan (Roeber & Bricker, 2015;Soria et al., 2018).
Geologic effects of outer-rise faulting and tropical cyclones may be mistaken for those of megathrust faulting along the Puerto Rico Trench.There, Mesozoic lithosphere of the North America plate converges obliquely, at 2 cm/yr, with the Caribbean plate.The resulting megathrust has been posited as a source of tsunami hazards not just in the Caribbean near-field but also on the Atlantic seaboard of the United States and Canada (ten Brink & Lin, 2004;Grilli et al., 2022).But no Puerto Rico Trench earthquake, back to about 1500 C.E., has a confirmed size of M8.0 or larger (ten Brink et al., 2011); the plate convergence is highly oblique; and this part of the megathrust is weakly coupled (DeMets et al., 2010;ten Brink & López-Venegas, 2012;Symithe et al., 2015).On the outer wall of the trench, scarps and grabens conspicuously attest to normal faulting of unknown recency (Andrews et al., 2013;ten Brink, 2005)-potential earthquake and tsunami sources that some Caribbean hazard assessments have neglected (Mueller et al., 2010;Zimmerman et al., 2022) and tsunamis (IOC, 2019).Geologists have struggled, in some instances, to distinguish between tsunamis and storms by means of their Caribbean deposits (Engel et al., 2016).The Caribbean history of catastrophic storms continued in 2017 with hurricanes Irma and Maria (Spiske et al., 2021).This paper aims to clarify earthquake and tsunami potential along the Puerto Rico Trench by using numerical models to simulate hypothetical causes of a catastrophic flood from the sea that is evidenced geologically and which took place during the last centuries before Columbus.The flood was ascribed to either a tsunami or an unusual storm on the basis of deposits and scours on Anegada (Atwater et al., 2017), a low carbonate island perched beside the trench 125 km to its south (Figures 1 and 2).The same event, interpreted as a tsunami, has been invoked to explain deposits on St. Thomas (Fuentes et al., 2017) and Anguilla (Biguenet et al., 2021).Preliminary model results for Anegada favored Puerto Rico Trench faulting as the cause for severe flooding over an extreme hurricane or a transatlantic tsunami and pointed to an outer-rise source (Buckley et al., 2012).In further tsunami simulations, either megathrust or normal faulting in the Puerto Rico Trench was found compatible with additional geologic evidence from Anegada (Lodhi, 2020) and from St. Thomas and Anguilla as well (Cordrie et al., 2022).The Anegada evidence has been ascribed to flooding of three generalized ranks in all (Atwater et al., 2017).Our simulations, which are tied to an updated compilation (Atwater, 2023), seek to explain evidence for the greatest of this flooding-a catastrophe from about six centuries ago.We screen modeled inundation against observed coral boulders and cobbles that were deposited hundreds of meters southward from the north shore, away from the Puerto Rico Trench (Figure 2).
The paper builds on previous work in three main ways.First, we try (and fail) to explain Anegada's coral erratics by means of Haiyan-like bores from a particularly impactful hurricane of category 5. Second, we refine rough estimates of earthquake size by screening a total of 25 Puerto Rico Trench tsunami sources-19 megathrust, 6 outer rise-both by the coral clasts and by scaling relationships among rupture dimensions and seismic slip.Finally, as measures of hazard, we estimate earthquake recurrence intervals and highlight quick tsunami arrival times.

Historical Hurricanes
A catalog of tropical cyclones identifies 98 that have passed with 100 km of the British Virgin Islands since 1851 (https://coast.noaa.gov/hurricanes/).Thirty-nine of these attained hurricane force, including cyclones of category 5 and eight of category 4.
Among the hurricanes of category 5, two skirted the British Virgin Islands but made landfalls in southeast Puerto Rico.The third, 2017 Hurricane Irma (Figure 1), passed through the British Virgin Islands themselves; it crossed Tortola with a maximum wind speed of ∼290 km/hr and a minimum atmospheric pressure of 915 hPa.Traces of Irma inundation on the north shore of Anegada were surveyed in 2018 at elevations as high as 3.8 m (Spiske et al., 2021).
At category 4, 1960 Hurricane Donna crossed Anegada itself.An eyewitness interviewed in 2008 recalled wading in storm surge that was driven into The Settlement (location, Figure 2b) by winds of the trailing left quadrant (Atwater et al., 2012).Less flooding in that same area resulted from 2010 Hurricane Earl, also of category 4 as it passed Anegada 30 km to the island's north (Atwater et al., 2014).Storm surge in The Settlement reached heights of 2.5 m (Donna) and 1.5 m (Earl) above sea level.Neither 1960 Hurricane Donna (Atwater et al., 2012), nor 2010 Hurricane Earl (Atwater et al., 2014), nor even 2017 Hurricane Irma of category 5 (Spiske et al., 2021), built overwash fans on Anegada that extended farther inland than about 50 m.Irma's most dramatic geologic effect on Anegada was to shift inland, by 10 m, the low crest of a pre-existing coral rubble ridge on the north shore at Soldier Wash.

Historical Tsunamis
Faulting along the Puerto Rico Trench, whether at the plate boundary or in the descending North American Plate, is not known to have produced a tsunami in Caribbean written history that goes back a little more than 500 years.The historical catalog tells instead of transatlantic waves of the 1755 Lisbon tsunami, and of tsunamis from sources within the Caribbean Plate.Accounts of the Lisbon tsunami, which has a poorly understood source offshore Iberia (Text S2,Figures S2.27,S2.28 in Supporting Information S1), are available for six Caribbean islands (Figure 1; Barkan et al., 2009;Lander et al., 2002;McCann et al., 2011).
Normal faults on the outer wall of the Puerto Rico Trench are evidenced less by historical earthquakes than by deep-sea bathymetry.Scarps up to 1,500 m high corrugate the outer rise between 62°W and 66°W (Figure 1).Some bound half-grabens containing as much as 1 km of sedimentary fill (ten Brink, 2005).A moderate outer-rise earthquake-a 1939 shock centered north of the trench at the longitude of eastern Puerto Rico-apparently resulted from normal slip and attained M6.4 (Doser et al., 2005).Though absent from most hazard assessments, outer-rise faults north of Puerto Rico were included as conjectural tsunami sources in preliminary inundation mapping of the island (Mercado Irizarry & Sepulveda, 2003).
Among historical tsunamis from intraplate sources in the northeast Caribbean, perhaps the largest was generated during the large Virgin Islands earthquake of 1867 (Reid & Taber, 1920).It scarcely affected north-and eastfacing shores fronting on the Atlantic Ocean because the tsunami source was a fault or faults in the Anegada Passage, south of St. Thomas (Reid & Taber, 1920;Barkan & ten Brink, 2010).A 1918 tsunami in northwest Puerto Rico (Reid & Taber, 1919) resulted from a large earthquake in nearby Mona Passage (Doser et al., 2005;Mercado Irizarry & McCann, 1998;ten Brink et al., 2023).Faulting in or near the Lesser Antilles arc likely set off a 1690 tsunami that reached St. Thomas and Nevis (O'Loughlin & Lander, 2003).
We know of no reports of tsunamis in the past five centuries at San Juan, Puerto Rico.This negative evidence is telling because San Juan faces the Puerto Rico Trench, stands on low ground, and has a rich written history of earthquake damage (McCann et al., 2011).Major earthquakes north of Puerto Rico and the Virgin Islands in 1785 and 1787 did not produce documented tsunamis (ten Brink et al., 2011;and references therein), nor did a large 1974 earthquake from deep normal faulting (McCann et al., 1982).
A submarine slide in Mona Passage was recently discounted as a source for the 1918 tsunami in western Puerto Rico (ten Brink et al., 2023).However, both the Puerto Rico Trench and Anegada Passage have steep submarine slopes where tsunamis may be triggered by submarine landslides.Slides on the trench wall nearest Puerto Rico have approximate areas of 45-160 km 2 and estimated recurrence intervals of 70,000-253,000 years (ten Brink et al., 2006).No turbidites younger than 18-19 ka have been identified on the trench floor north of Puerto Rico (Chaytor & ten Brink, 2014).
The Anegada evidence has been calibrated in the field to traces of modern hurricanes and of a sea flood between 1650 calibrated years (cal yr) C.E. and 1800 C.E. Shore-hugging sand and coral rubble resulted from storms of the past century or more (Atwater et al., 2014;Spiske, 2016;Spiske & Halley, 2014;Xu et al., 2015).These modern storm deposits extend onshore several tens of meters at most.Not even 2017 Hurricane Irma moved coarse clasts farther inland on Anegada (Spiske et al., 2021).As observed elsewhere (Ferrario et al., 2014), storm waves can dissipate on coral reefs and on shallow flats between reef and shore.The sea flood of 1650 cal yr C.E. and 1800 C. E., probably the 1755 Lisbon tsunami, is evidenced most visibly by breaches in low beach ridges, and stratigraphically by an inland sheet of sand and shells that was derived mostly onshore (Atwater et al., 2012;Pilarczyk & Reinhardt, 2012;Reinhardt et al., 2012).
Appraised in this paper is the far greater precolonial flooding that is evidenced on Anegada most conspicuously by scattered coral boulders and cobbles and by inland fields of limestone boulders and cobbles.Most of the coral clasts are brain coral, Pseudodiploria (formerly Diploria) strigosa.Their distribution, along with that of anomalous limestone boulders and cobbles, is shown in Figure 2. Radiocarbon ages of inland coral clasts date their southward emplacement solely to the time window 1200-1480 cal yr C.E. (Atwater et al., 2017, p. 325), even though relative sea level in the northeast Caribbean has been nearly stable in the past 3,000-4,000 years (Khan et al., 2017, their Figure 6).On that basis, just one southward sea flood in recent millennia managed to overwhelm much of Anegada.
Anomalously coarse-grained deposits in that same time window have been identified beneath salt ponds of St. Thomas and Anguilla.Fuentes et al. (2017) found most layers of overwash sand beneath coastal ponds on St. Thomas and nearby cays.They attributed these layers to as many as five tsunamis, one of which Tuttle and Fuentes (2021) dated to 1280-1390 cal yr C.E.On Scrub Island, Anguilla, the thickest and coarsest salt-pond deposit identified by Biguenet et al. (2021) has been dated to 1364-1469 cal yr C.E.These St. Thomas and Anguilla deposits have yet to be compared with overwash deposits on Barbuda that Biguenet et al. (2023) ascribed to 2017 Hurricane Irma.

Previous Models
Limestone boulders on Anegada provided the basis for screening by Buckley et al. (2012) of potential causes of the late precolonial flooding.They used the largest of the boulders that Watt et al. (2012) had surveyed south of Windlass Bight (fields labeled "w" in Figure 2a).Simplified equations for boulder transport gave estimates for the minimum flow velocities required to move the largest of the clasts.The thresholds were then compared with simulated flow depths and flow speeds from a one-dimensional Boussinesq-type numerical model.Models were made for a category-5 hurricane, a Lisbon tsunami, and two tsunamis from faulting in the Puerto Rico Trench.For one of the near-field tsunami simulations, the hypothesized source has 8 m of thrust displacement on a rupture 640 km by 90 km.Its rupture extends southeastward from the longitude of Anegada into the Lesser Antilles Subduction Zone, and it corresponds to a thrust earthquake of M8.7.The other near-field source is centered due north of Anegada.It has 7 m of slip on a normal-fault rupture 130 km by 45 km, on the outer wall of the Puerto Rico Trench, for an earthquake of M8.0.Buckley et al. (2012) found that the outer-rise scenario best explains the transport of the largest of the limestone boulders, and that the thrust scenario may also account for them.
Limestone boulders on Anegada were also used in tsunami-source screening by Lodhi (2020).She focused on boulders beside southeast Red Pond and in the East End (locations, Figure 2b).A revised equation for initiation of boulder transport (Lodhi et al., 2020) gave estimates of flow speeds that she applied to limestone knolls from which the modeled boulders trail.Comparing these estimates with flow speeds in simulations of 13 tsunami scenarios, Lodhi (2020) found agreement with earthquakes no smaller than M8.4 on the subduction thrust and a minimum of M8.1 on the outer rise.
Simulations by Cordrie (2021) and Cordrie et al. (2022) evaluated tsunamis from a total of 37 hypothetical sources, both intraplate and interplate, in the northeast Caribbean.They compared model results with sedimentary deposits on Anegada, St. Thomas, and Anguilla.From Anegada they made use of brain-coral erratics and bioclastic carbonate sand.Among the tsunami sources considered, their models explain evidence from the three islands most completely with a subduction-thrust rupture (their scenario T2) having 100 km down-dip width, 200 km along-strike length, and 20 m average dip-slip displacement.They found the deposits less consistent with a normal-fault rupture (their OR3) 150 km by 25 km with 10 m slip.These scenarios correspond to earthquakes of M8.7 (megathrust) and M8.0 (outer rise).

Roadmap
This paper relates brain-coral erratics on Anegada to simulations of an extreme hurricane and of 25 hypothetical tsunamis.For the hurricane modeling (Section 3.2), we first modeled the paths, wave heights, surge levels, and wave periods of 15 hurricanes of category 5 (Figure 3).We tested these storm models against observations from 2017 Irma (Text S3 in Supporting Information S1).Then, for the most extreme of the 15 storms, we fed its offshore waves into high-resolution simulations that produced surf beat.Finally, we compared the resulting flooding with the mapped distribution of coral boulders on a part of the island where an absence of tall, sandy beach ridges probably gives all the modeled storm waves full opportunity to move coral clasts inland, without first eroding such ridges (Figure 4).
For the tsunami modeling (Section 3.3), we outlined eight rupture segments on the megathrust and one on the outer rise (Figure S3.1 in Supporting Information S1 and Table 3).We then devised 25 tsunami scenarios in all by varying megathrust-segment combinations and amounts of slip.For each of these scenarios we computed inundation on Anegada to obtain flow depths at brain-coral sites island-wide (Figure 5).We also screened the parent earthquakes in these 25 scenarios by means of scaling relationships (Figure 6).As due diligence we modeled, in addition, a generous source for the 1755 Lisbon tsunami (Text S2, in Supporting Information S1).The storm models and the tsunami models all employ gridded bathymetry and topography (Section 3.5).Simplifying assumptions contribute to uncertainty in the bottom-line estimates of megathrust and outer-rise earthquake size (Section 3.6).
The geologic screening approach in Figure 5

Modeling Code
We modeled storm wave dynamics in case some combination of storm-surge height and wave conditions would produce the Anegada catastrophe.Our storm-wave modeling follows the method of Cheung et al. (2003) and Li et al. (2014) that integrates circulation, spectral wave, and Boussinesq models to simulate the coastal inundation from storm tides to individual waves.
First, a circulation model based on the non-hydrostatic Evolution of Ocean Wave (NEOWAVE) (Yamazaki et al., 2009(Yamazaki et al., , 2011) ) computes the phase-averaging storm tide water level driven by the atmospheric forcing and the tidal surface elevation and velocities.Next, the spectral wave model SWAN (Booij et al., 1999), taking into account wave breaking and dissipation over fringing reefs in tropical environment (Filipot & Cheung, 2012), simulates energy growth, propagation, and dissipation from the open ocean to nearshore under wind forcing in a phase-averaging process.
Finally, Boussinesq approximations of Roeber et al. (2010) and a spectral wave-maker are used to resolve the surf and swash-zone dynamics wave-by-wave to realistically simulate the coastal inundation induced by hurricane surge and waves.The spectral wave-maker acquires the SWAN-computed wave spectrum at the 50-m depth offshore north of Anegada.This same modeling approach was previously used to simulate tsunami-like surf beats, or infra-gravity waves, that were observed during 2013 Typhoon Haiyan.These waves, which built up behind a coral reef, and which resembled tsunami bores, washed away houses 200 m inland at Hernani, Philippines (Roeber & Bricker, 2015), while also moving reef clasts and transporting fine sediments 1 km inland (Soria et al., 2018).

Validation With 2017 Hurricane Irma
We validated the storm models in Figure 3 by simulating effects of 2017 Hurricane Irma where they had been observed at deep-sea buoys, at tide gauges on Puerto Rico and the U.S. Virgin Islands, and at a coral-rubble berm at Soldier Wash.For these comparisons we used the best track of 2017 Hurricane Irma from historical hurricane tracks of the National Oceanic & Atmospheric Administration (NOAA) (https://coast.noaa.gov/hurricanes/).The results are shown in Figures S3.1-S3.3 in Supporting Information S1.In the deep sea, the simulation mostly reproduces the significant wave heights, peak periods, and mean directions observed at three offshore buoys 41043, 41044, and 42060 (Figure S3.1 in Supporting Information S1).At the tide gauges, the computed water levels slightly underestimate the observations (Figure S3.2 in Supporting Information S1)-a discrepancy is probably produced by a regional grid resolution (1 arc min, ∼1,850 m) too coarse to resolve the nonlinear wave dynamics, such as wave setup, in the nearshore area.
The wave conditions at 50 m depth offshore north of Anegada predicted by SWAN show a significant wave height of ∼10 m, a peak wave period of 12 s, storm tide water level of 0.6 m (the sum of storm surge and tidal level).These conditions provide boundary conditions for the wave-by-wave inundation simulation with the Boussinesq model.It shows that the hurricane-spawned waves, cresting 4 m above sea level, barely overtop the coral-rubble berm at Soldier Wash (Figure S3.3 in Supporting Information S1)-an estimate indistinguishable from an inundation level of 3.8 m above sea level that was surveyed in 2018 by Spiske et al. (2021).

Posited Storms
We set up 15 extreme scenario storms by using methods in Section 3.2.1.Each is a hurricane of category 5 that attains the lowest recorded Atlantic atmospheric pressure (882 hPa, during 2005 Hurricane Wilma) and the highest recorded tropical-cyclone wind speed (315 km/hr, in 2013 Typhoon Haiyan).
The storms follow five tracks in each of three swaths (Table 1, Figure 3a-3i).One of the three swaths crosses Anegada from southeast to northwest (SE-NW), as did 2010 Earl (Figure 3a).Another runs from northeast to southwest (NE-SW), as did 2007 Noel (Figure 3e).Along the third swath the storms proceed from southwest to northeast (SW-NE), as did 1999 Lenny (Figure 3i).The central track in each swath crosses Anegada.Flanking each central track are two to the left and two to the right.The flanking tracks are at distances of 20 and 40 km to improve the odds that Anegada is intersected by a 25-mile radius of maximum wind speed.

Modeling Code
We computed tsunami generation, propagation, and inundation using the Method of Splitting Tsunami (MOST) model.MOST is a suite of numerical simulation codes using a finite-difference approximation of the characteristic form of the depth-averaged nonlinear shallow water equations (Titov & Synolakis, 1995, 1998).For propagation, MOST solves the shallow-water wave equations in spherical coordinates.The physical process of frequency dispersion can be approximated by numerical dispersion in MOST (Burwell et al., 2007;Zhou et al., 2012).For inundation, MOST employs unidirectional coupling with computational grids to telescope down to the high-resolution area of interest for inundation computation.The outer grid provides the inner grid with computed wave amplitudes and flow velocities at all four boundaries by linear interpolation from coarser resolution, while the inner grid has no effect on the outer grid.MOST has been validated exhaustively (Synolakis et al., 2008;Tang et al., 2009Tang et al., , 2012;;Titov, 2009;Titov et al., 2005;Uslu et al., 2011;Wei et al., 2008Wei et al., , 2013;;Zhou et al., 2012).
The initial condition in MOST replicates, at the sea surface, the vertical coseismic deformation of the ocean floor.This deformation is computed with the elastic earthquake-faulting theory of Okada (1985).The earthquake rise time is not considered; coseismic uplift and subsidence are instantaneously transferred to the ocean surface.
MOST uses the Manning's coefficient to account for frictional resistance to flow across the bed.We used a uniform Manning's coefficient 0.03 m 1/3 /s throughout all grids-a simplifying assumption discussed in Section 3.6.1.

Megathrust Sources
The subduction thrust provides the most extensive of the potential tsunami sources along the Puerto Rico Trench.Changes in the orientation of the frontal thrust in the multibeam bathymetry data (Andrews et al., 2013) provide the basis for the eight megathrust rupture segments in Figure 5 and Table 2.The fault plane in each case is assumed to dip 15°.Each scenario name identifies segments by letter.For each of six segment combinations we devised three scenarios (four for segment combination CDE) that differ in fault slip, giving a total of 19 megathrust scenarios (Table 3).The fault slip, simulated as uniform, is in the range 3.5-7.5 m, stepped mainly at increments of 1.1 m.Each scenario name includes a corresponding earthquake size, computed with a bulk rigidity of 3 × 1010 Pa.The six segment combinations, from west to east: ABCDE spans most of the length of the Puerto Rico Trench.The rupture extends about 570 km between the longitudes of the eastern Dominican Republic and Anguilla.We considered slips of 5.2-7.5 m that correspond to M8.57-M8.68.

Outer-Rise Sources
Segment N provides a hypothetical normal-fault rupture 130 km long that is centered north of Anegada among scarps and grabens on the outer wall of the Puerto Rico Trench (Figures 1 and 5, Figure S3.1 in Supporting Information S1).The fault plane dips 60°south.Uniform slip of 3.5-14.1 m corresponds to earthquake magnitudes of M7.8-M8.24(Table 3), given a conservative bulk rigidity of 3 × 10 10 Pa.

Iberian Source
Text S2 in Supporting Information S1 illustrates the Iberian source that we used to simulate the 1755 Lisbon tsunami.In brief, as shown in Figures S2.27 and S2.28 in Supporting Information S1, we modified the Barkan et al. (2009) source by increasing the earthquake magnitude from 8.7 to 9 following Muir-Wood and Mignan ( 2009) and by increasing rupture length along strike to 600 km to compensate for using a lower rigidity (3 × 10 10 Pa) than did Barkan et al. (2009).This simple, generous source violates geologic constraints by running across faults mapped by Zitellini et al. (2009) (Figure S2.27 in Supporting Information S1).

Local Geology
We screened all 25 near-field tsunami scenarios by now completely and deeply they drown the observed Anegada sites of 171 coral boulders and cobbles to minimum flow depths of 1 and 2 m.This screening, with caveats in Section 3.6, provides a rough guide to scenario success and failure.The smallest slip that drowns all the brain corals to a simulated depth of 1 m or more provides an uncertain estimate of the smallest earthquake required.Substituting a minimum flow depth of 2 m or more increases those magnitude estimates.
All of those 171 boulders and cobbles are clasts of brain coral.We use brain coral for two reasons.First, it has smoothly curved margins, similarly curved growth bands, and cross-cutting radial structure that combine to make this kind of coral erratic easy to spot.Second, unlike anomalous deposits of limestone or sand, brain-coral clasts can be dated directly, and all 26 of the radiocarbon ages obtained on them are consistent with emplacement in a single extraordinary sea flood in 1,200-1,480 cal yr C.E. (Atwater, 2023;Atwater et al., 2017).Some of the boulders consist of entire coral colonies.Measured lengths of 23 of the boulders, compiled in Atwater (2023), are in the range 1.0-2.7 m.Given a density of 2040 kg/m 3 for water-saturated brain-coral (Spiske et al., 2021, their Table S1 in Supporting Information S1), the largest boulders when freshly deposited weighed two tonnes or more.Sixteen of the 171 clast sites are walls or piles to which the clasts were moved, probably by less than 50 m, during eighteenthor nineteenth-century land clearing.All 171 sites are plotted in Figure 2a. Figure 5 summarizes the percentages flooded to minimum depths of 1 m and 2 m-measures of scenario success displayed in map view in Supporting Information S1 Text S2, Figures S2.2-S.2.26.

Global Scaling Relationships
Figure 6 screens hypothetical tsunami sources by means of scaling relationships of earthquake source parameters.Empirically, rupture length and width scale broadly with earthquake magnitude (Blaser et al., 2010;Strasser et al., 2010;Thingbaijam et al., 2017;Wells & Coppersmith, 1994).The scaling relationship from Strasser et al. (2010) does not distinguish between reverse and normal faults, but all the others do.The Puerto Rico Trench scenarios plotted in Figure 6 include all 25 in Figure 5, along with the megathrust and outer-rise scenarios in Buckley et al. (2012) and Cordrie et al. (2022).

Bathymetric and Topographic Grids
The model grids used for both tsunami and storm surge modeling are derived from digital elevation models constructed from bathymetric and topographic data identified in Supporting Information S1 Text S1 (Figure S1, Table S1).Telescoped grids for four computational domains were used in MOST to simulate the tsunami wave dynamics throughout the generation, propagation from deep water to nearshore, and the inundation onshore.The spatial grid resolutions are 30 arc sec (∼900 m), 6 arc sec (∼180 m), 1 arc sec (∼30 m) and 1/3 arc sec (∼10 m) covering computational domains for the entire Caribbean Sea, the northeastern Caribbean, the British Virgin Islands, and Anegada, respectively.
For hurricane modeling, the circulation model and the spectral model are both simulated in two nested grids at the spatial resolution of 1 arc min (∼1,800 m) for deep ocean in the northeastern Caribbean and the neighboring Atlantic waters (Figure S3.1a in Supporting Information S1), and at the spatial resolution of 6 s (∼180 m) for the nearshore area of Anegada (Figure S3.1b in Supporting Information S1).In the phase-resolving Boussinesq simulation, we use a high grid resolution of 5 m at the study site, Soldier Wash on the north shore of Anegada, that is needed to resolve individual swell waves (Figure S3.3 in Supporting Information S1).

Limitations
Notably excluded from the methods above is comparison with threshold velocities for boulder transport.Such comparisons go beyond the scope of this report, for several reasons.First, in tsunami models uncertainties in simulated flow velocities are commonly larger, by 20%-200%, than are water-surface elevation and flow depth, if benchmarked to physical experimental and field measurements (Lynett et al., 2017).Reducing these uncertainties in simulated flow velocities requires local field calibrations.
Even with our simplifying use of simulated flow depth as a screening tool, Anegada geology gives rise to uncertainties that have opposing effects on estimated earthquake size.On the one hand, earthquake magnitude may be underestimated where roughness and vegetation prevented real flows from reaching the depths of the simulated ones.Magnitude may be also underestimated through neglect of beach-ridge erosion-and perhaps also of wave dispersion.These complications particularly affect the coral clasts on which we and Cordrie et al. (2022) rely.In two previous Anegada reports, simulated flow speeds were compared with threshold velocities for transporting limestone boulders (Buckley et al., 2012;Lodhi, 2020).For those boulders, the probable sources of the modeled limestone boulders were observed nearby, onshore, where they had been reached by flows that crossed or breached low beach ridges and continued, without much further obstruction by woody plants, into interior ponds and salt flats.By contrast, sediment-transport modeling of the coral clasts would need to contend with their entrainment offshore, their transport over tall beach ridges along the Atlantic shore of southeast Anegada, and their continued transport on wooded, karstic limestone uplands.

Vegetation and Roughness
Woody plants on Anegada probably both impeded the precolonial flows and shortened landward distances of precolonial brain-coral transport.An early sketch of Anegada shows woody vegetation on its skyline (Jefferys, 1775).Trees and shrubs probably also abounded in precolonial times, when human use of the island may have focused on nearshore fishing of the queen conch (Davis & Oldfield, 2003).Coral erratics in transport were likely slowed, and perhaps sieved or blocked as well, by trunks and branches that withstood the flooding or, having been broken off, built sieve-like barriers against obstacles farther inland.The observed distribution of coral boulders in Figure 2 may then understate the inland extent of flows that were otherwise capable of moving them.
Our use of a Manning's coefficient of 0.03 m 1/3 /s throughout all grids (Section 3.1.1)would be appropriate for floodplains with light grass or few crops (Chow, 1959).But this roughness value may permit greater simulated inundation than would higher Manning's coefficients suitable for landscapes that are shrubby or forested (Figure 2b), corrugated by beach ridges (Figure 2a), and (or) pitted by limestone karst (sinks in Atwater, 2023).We allowed roughness to vary across the Anegada landscape in a run of N M8.1 (Figure S4.1 in Supporting Information S1).We followed the guidance of Arcement and Schneider (1989) to assign suitable Manning's coefficients by taking into account surface irregularity, shape and size of the cross section, obstruction, and vegetation.We assigned values of 0.02-0.03 to areas now covered by near-sea-level salt ponds, and values of 0.04-0.05 to wooded uplands on karstic limestone and to tall beach ridges beside the north shore.The results show that more realistic values for roughness locally reduce maximum flow depth (Figure S4.2b in Supporting Information S1, circled areas) while scarcely affecting maximum flow speed (Figure S4.2d in Supporting Information S1).

Eroded Landforms
The sea flood roughly six centuries ago variously beveled and incised sandy beach ridges along the north shore of Anegada (Atwater et al., 2017, p. 313).This erosion, if accomplished during an initial wave, may have enabled later waves to run farther inland, and to move entrained clasts there more readily.Such dynamic topography, if not accounted for, can allow a tsunami model to underestimate inundation and its effects Tehranirad et al. (2017).In our tsunami simulations, that outcome could result in overestimating earthquake size, by not just by understating the inland extent of the precolonial sea flood, but also by permitting the mapped boulders to have traveled farther inland than they otherwise would have.
We tried to limit the effects of dynamic morphology in our test of simulated flooding by the extreme hurricane.We did so by evaluating detailed storm-inundation modeling at Soldier Wash and Deep Bay (locations, Figure 2b).In both these areas, flow from the Atlantic can reach inland limestone without first encountering sandy beach ridges that rise 6-8 m above sea level-beach-ridge heights typical along much of the island's Atlantic shore (Figure 2a).The shore at Soldier Wash is instead flanked by a modern coral-rubble berm that crests 4 m above sea level (Spiske, 2016;Spiske & Halley, 2014)-no higher than a broad limestone ridge behind, across which precolonial coral boulders and cobbles were scattered (Atwater et al., 2017, their Figure B4).At Deep Bay, a modern sandy beach ridge similarly rises 4 m from sea level, and a limestone plain behind is littered with coral clasts (Atwater et al., 2017, their Figure B5).

Wave Dispersion
By using depth-averaged flow velocity, MOST does not directly take wave dispersion term into account.Although the numerical scheme can mimic physical dispersion to some degree by tuning the relationship between spatial and temporal steps (Burwell et al., 2007), the model as we used it does not fully solve wave dispersion.However, because wave dispersion increases with propagation distance (Glimsdal et al., 2013), it is unlikely to have much effect on our simulations of near-field waves from Puerto Rico Trench sources.

Simulated Storms
The highest offshore waves are generated by the extreme hurricane that follows track SE-NW 3 (Figure 3).This track parallels that of 2010 Earl, but instead of passing to the north the hypothesized storm makes landfall on the island.At the 50-m depth offshore north of Anegada, the SE-NW 3 category-5 hurricane scenario produces waves with 13 m significant height (Figure 3b) and 15.6 s period (Figure 3d).The sum of storm surge level (1.33 m, Figure 3c) and an assumed high ambient tide level (0.2 m) gives a 1.53-m storm tide level relative to Mean Sea Level (MSL).This storm tide level serves as the vertical datum in Figure 4b; for gauge g2, near Soldier Wash, the peak amplitudes accordingly correspond to an elevation close to 4.5 m above MSL, or about 0.5 m above the crest height (in lidar topography) of the coral-rubble berm of Spiske and Halley (2014) and Spiske (2016).
The spectral wave-maker and the Boussinesq approximations succeed in producing surf beat, much as did these methods in simulating tidal bores in Philippines from 2013 Typhoon Haiyan (Roeber & Bricker, 2015).The extreme hurricane waves of our SE-NW 3 scenario produced surf beat by piling up, in sets, on shoals of Windlass Bight.
Not even this extreme scenario, however, floods the sites of a majority of the dozens of coral erratics mapped inland of the low coastal berms at Soldier Wash and Deep Bay (Figure 4). Figure 4b plots three modeled time series of the water surface elevation: one in deep water (g1); another in shallow water (g2) offshore Soldier Wash, where the barrier reef grazes the shore; and in shoals of Windlass Bight (g3), which separate that reef from the island's north shore.The time series show little surf beat at Soldier Wash.By contrast, infra-gravity waves of Windlass Bight have distinct periods of 2-3 min.
The Windlass Bight time series confirms that our high-resolution, phase-resolving modeling can reproduce surf beat akin to that observed and modeled for Haiyan.But the paucity of modeled surf beat at Soldier Wash explains why, even with the advantage of 1.5-m storm tide, the simulated storm inundation reaches few of the coral clasts inland from there (Figures 4c and 4d).Likewise at Deep Bay, the modeled inundation for the extreme SE-NW 3 scenario falls short of most of the precolonial brain-coral boulders.Neither of these areas is sheltered by beach ridges tall enough to block the simulated storm waves (Section 3.6.2).

Simulated Tsunamis
The geologic screening in Figure 5 casts doubt on any of the 25 near-field tsunami scenarios that correspond with earthquakes smaller than M8.1 (outer rise) or M8.4 (megathrust).The smallest tsunamis that flood all the coral-clast sites to depths of 1 m or more correspond to earthquakes of M8.17 (11.4 m slip on the outer-rise normal fault) and M8.45 (7.5 m slip on the combined thrust segments C, D, and E).The further geologic screening in Figure S2.29 in Supporting Information S1 disallows the tested M9 Lisbon source.

Outer-Rise Scenarios
Tsunami scenario N M8.17 (Figures 7 and 8) achieves flow depths of 1 m or more at 100% of the Anegada braincoral sites and depths of 2 m or more at 91% of them (Figure 5b).Taken at face value, this scenario suffices to explain the Anegada distribution of brain-coral boulders on Anegada (Figure 8).Uncertainties noted in Section 3.2.3probably permit, as alternatives, scenarios N M8.10 (8.8 m slip) and N M8.24 (14.1 m slip).In M8.24, for instance, simulated flow depths are no less than 2 m at 99% of the coral-clast sites.None of the outer-rise scenarios violates scaling relationships for earthquakes on normal faults (Figures 6c and 6d).
Smaller outer-rise earthquakes might produce successful scenarios under different assumptions about the fault ruptures.Our outer-scenarios presuppose uniform slip on fault ruptures with an upper limit is 5 km below the ocean floor (Table 2).Non-uniform slip, greater near the surface than at depth, and a shallower upper limit would both have the effect of generating larger tsunami waves of shorter period.

Megathrust Scenarios
Tsunami scenario CDE M8.45 (Figures 9 and 10) floods all the brain-coral sites to depths of 1 m or more, and it drowns 95% of them to minimum depths of 2 m (Figure 5c).The scenario has 7.5 m of slip on a rupture 255 km long, as a great subduction-thrust earthquake along the eastern Puerto Rico Trench centered north of Anegada.A lesser alternative, CDE M8.40, with 6.3 m of slip, comes close at 98% (1 m) and 72% (2 m), while drowning 83% (1 m) and 44% (2 m) may cast doubt on CDE 8.34, with 5.2 m slip.
Among the full set of 19 megathrust scenarios, simulated flooding on Anegada is insensitive to rupture far along strike, either westward to the Dominican Republic or southeastward to Guadeloupe.Adding western segments A and B to the CDE rupture, in scenario ABCDE M8.68, while retaining 7.5 m of slip, yields flooding to minimum depths of 1 and 2 m at 100% and 95% of the brain-coral sites.Similarly, with CDE extended southeastward as CDEFGH M8.76, again with 7.5 m of slip, those figures are 100% and 90%.Conversely, segment C, which extends 130 km westward from the longitude of Anegada, is shown to be important by experiments that include or exclude it.At the same slips, CDE scenarios succeed far more than do the DE scenarios.CDEFGH similarly outperforms DEFGH.Segment E is also important in flooding of Anegada, but less so, as judged by comparing histograms for CD with those for CDE.
As for scaling relationships, the CD and CDE scenarios fit comfortably in regressions of rupture length and rupture width against thrust-earthquake magnitude (Figures 6a and 6b, points no.4-6 and 7-10).The poor fit with rupture length in the CDEFGH scenarios (no.11-13) reflects the modest slips assumed and the consequent upper bound of M8.76.None of the megathrust scenarios in Figure 5 violate scaling of earthquake size with rupture width in Figure 6b.

Outer-Rise Versus Megathrust
Screening in Figures 5 and 6 provides little if any basis for favoring outer-rise N M8.17 over megathrust CDE M8.45, or vice versa.In Figures 8 and 10, the main difference is greater flooding of south-central Anegada under the megathrust scenario.
Coseismic land-level change, if detected geologically, would provide an observational basis for favoring the megathrust over the outer rise.The expected subsidence at Anegada, while near zero in N M8.17 (Figure 7b), is close to three-quarters of a meter in CDE M8.45 (Figure 9b).The simulated flooding takes account of this difference in modeled coseismic subsidence on Anegada.But whether the island subsided suddenly some six centuries ago is an open geologic question.

Comparisons With Previous Findings
Our screening results broadly accord with previous field observations and numerical simulations as evidence for late-precolonial faulting along the Puerto Rico Trench.This overall agreement is obtained despite differences in methods.
Our results break new ground by discounting even Haiyan-like surf beat as a cause of the Anegada catastrophe.In the modeling of Buckley et al. (2012), the simulated hurricane was of category 5, but it did not produce surf beat.Consequently, in light of tsunami-like bores from surf beat during 2013 Typhoon Haiyan, it was difficult to rule out an unusual storm as the cause of the extraordinary precolonial sea flood on Anegada (Atwater et al., 2017).But now, with coral clasts beyond the reach of the surf beat simulated here as well as the overwash during 2017 Hurricane Irma (Spiske et al., 2021), we recommend putting the storm alternative to rest.
Our minimally successful tsunami scenarios of M8.1 (outer rise) and M8.4 (megathrust) closely resemble those of Lodhi (2020).Her fault segments were derived from ours, but her ground truthing was different: she compared flow speeds simulated in GeoClaw tsunami models against flow speeds computed for the onshore initiation of limestone-boulder transport.
As for the Iberian alternative, nearly all the Anegada brain-coral erratics lie beyond the inland limit in our Lisbon M9.0 simulation (Figure S2.30 in Supporting Information S1).This finding agrees most directly with radiocarbon ages on these erratics, none of which are young enough to overlap with 1755 C.E. (Atwater et al., 2017).Also in agreement are two previous models that compare Iberian scenarios with Anegada limestone boulders.Onshore flow in one of these simulations, in Buckley et al. (2012), was found insufficient to transport the largest limestone boulders that Watt et al. (2012) measured south of Windlass Bight.The simulated Lisbon flooding in Lodhi (2020, p. 84) reaches but half of the 633 anomalous limestone clasts mapped in Figure 2a.

Precolonial Inundation on Other Islands
The simulations for outer-rise scenario N M8.17 and megathrust scenario CDE M8.45 subject the British Virgin Islands to extreme waves (Figures 7e and 7f, 9e-9f, and 11).In both scenarios, a high-energy beam trends southeastward toward the passage between Tortola and Virgin Gorda.The simulated wave heights for N M8.17 are mostly lower than those in CDE M8.45.In both simulations the tsunami waves arrive at all islands with a leading depression (Figure 11)-an arrival that can presage larger runup than in a tsunami with a leading elevation (Tadepalli & Synolakis, 1994).
Two supplementary figures provide a three-way comparison among maximum coastal flow depths achieved by near-field scenarios CDE M8.45 and N M8.17, and by far-field scenario Lisbon M9.0-onSt. Thomas (Figure S2.29 in Supporting Information S1), and in the vicinity of San Juan, Puerto Rico (Figure S2.30 in Supporting Information S1).The St. Thomas flow depths for the near-field scenarios are probably consistent with precolonial overwash that Fuentes et al. (2017) inferred from salt-pond stratigraphy.None of the three scenarios cause much flooding near San Juan.The simulated flooding at San Juan is greatest overall for megathrust scenario CDE M8.45, probably because the deep water of Puerto Rico Trench impedes and diverts tsunami waves cross it (Mofjeld et al., 2000).Among the three scenarios, only CDE M8.45, generated south of the trench, escapes this effect.The simulated flooding in the San Juan area would surely be much greater if sources comparable to CDE M8.45 and N M8.17 were placed due north of the city.

Tsunami Arrival Times
As near-field tsunami scenarios, both N M8.17 and CDE M8.45 leave little time between parent earthquake and tsunami arrival.Just offshore northern Anegada, the simulated CDE M8.45 tsunami begins with a leading depression 7 min after the earthquake, followed at 15 min by a positive wave 6 m high (Figure 11, virtual gauge 1).Arrival times under an hour prevail among almost all coastlines in the northern Caribbean, from Puerto Rico to Montserrat.First arrivals reach many of the islands sooner under the megathrust scenario than with the outer-rise scenario, particularly where the megathrust rupture is closer (as in Figure 9) than is the outer-rise rupture (Figure 7).The time difference at Anegada amounts to 6 min on the north shore (virtual gauge 1) and 25 min on the south shore (virtual gauge 2).The megathrust tsunami also precedes the outer-rise tsunami at most other virtual gauges that face the Puerto Rico Trench (gauges 4-9).The waves in both scenarios reach the south shores of many islands after proceeding through the passage between Anegada and Anguilla (Figures 7e and 7f, 9e and 9f).

Earthquake Recurrence
The tsunami on Anegada about six centuries ago was an unusual event, whatever its source.As noted above, radiocarbon ages of stray coral clasts imply that it far exceeded all others on the island in recent millennia   5a.Longitudes and latitudes are defined at the center point of the down-dip side of the subducting fault.Minimum depth refers to the distance between the sea floor and the up-dip limit of fault rupture.
(Section 2.3.1).Furthermore, as shown here, the earthquakes inferred from the successful scenarios have recurrence intervals thousands of years long, whether on the megathrust or the posited normal fault.
Megathrust recurrence intervals of 1,000 years or more can be deduced from the posited slip of 7.5 m in megathrust scenario CDE M8.45.The intervals range from 1,250 years if that slip expends an orthogonal component of plate convergence that accumulated at 6 mm/yr.The interval is three times longer, 3,750 years, if the trench-normal seismic-slip rate is 2 mm/yr.These estimates begin with the full convergence rate of ∼20 mm/y (DeMets et al., 2010), and with the convergence direction north of Anegada is 16º-21°from the trench axis.That obliquity gives an orthogonal component of interplate convergence close to 6 mm/yr.Grilli et al. (2022), using orthogonal convergence of 6-7 mm/yr, estimated return periods of 250, 550, and 1,000 years for Puerto Rico Trench earthquakes of M8.3 (slip 1.5 m), M8.7 (slip 3.7 m) and M9.0 (slip 7.1 m), respectively.But the fault loading depends on interplate coupling, which has been shown geodetically to be weak (ten Brink & Lin, 2004;ten Brink, 2005;ten Brink & López-Venegas, 2012;Symithe et al., 2015), A combinatorial optimization method gave a seismic slip rate along the Puerto Rico Trench of 2 ± 1 mm/yr, along with a mean return time greater than 1,000 years for earthquakes of M > 8.1 (Geist & ten Brink, 2021).
A longer recurrence interval is likely for the 11.4 m of normal slip under scenario N M8.17.In that case, the rate of fault loading can be taken as the rate of elongation along the surface of the flexed plate as the plate develops tensile bending stresses and normal faulting in its upper part before entering the trench.We follow Craig et al. (2014, their equation 2) in calculating the horizontal elongation rate, Xʹ, at the plate's upper surface (z = 0) Where V x is the plate convergence rate perpendicular to fault N in Figure S3.1 in Supporting Information S1, Z n is the depth to the neutral plane of the flexed plate (the depth where tensile bending stresses above change to compressive stresses below), and R c is the approximate radius of plate curvature in the faulted region.The plate convergence rate is 20 ± 0.5 mm/yr at azimuth of 254° (DeMets et al., 2010), and its perpendicular component to the modeled normal fault is 9.17 mm/yr.There are not sufficient earthquakes in the outer rise to determine the neutral plane from earthquake focal mechanisms, but the depth to the neutral plain under the outer rise of other subduction zones with similar crustal ages (Kuril and northern Tonga, Müller et al., 2008) was determined empirically to be ∼25 km (Craig et al., 2014).Plate curvature was calculated from three points along the flexed plate and gave an approximate radius of 233 km.The calculated maximum horizontal surface elongation rate is therefore close to 1 mm/yr.It will take approximately 5,800 years to accumulate 11.4 m of slip on a 60°inclined fault, necessary to generate the N M8.17 tsunami scenario.This is a rough estimate because of the above assumptions and because we neglected the possibilities of aseismic slip and of rupture on multiple faults.

Conclusions
Extensive geologic evidence for catastrophic prehistoric flooding on Anegada, British Virgin Islands was previously attributed to either a tsunami or an unusual storm (Atwater et al., 2017).Correlative overwash was inferred from salt-pond stratigraphy on St. Thomas, US.Virgin Islands (Fuentes et al., 2017) and in Anguilla (Biguenet et al., 2021).This and other extreme-wave geology provided ground truth for simulating tsunamis from hypothetical fault ruptures along the Puerto Rico Trench (Buckley et al., 2012;Cordrie et al., 2022;Lodhi, 2020).Note.Tests with coral clasts on Anegada are summarized in Figure 5.The relationship between rupture length, L, and magnitude follows Wells and Coppersmith (1994)  Here we have offered additional simulations to evaluate overwash on Anegada by extreme tropical cyclones and by thrust faulting and normal faulting in the trench.The most extreme of the hypothesized storms emulates 2013 Typhoon Haiyan by producing tsunami-like bores from surf beat.Nineteen of the tsunami scenarios posit faulting on the megathrust that coincides with the Puerto Rico Trench, and another six tsunami scenarios vary slip on a simplified high-angle normal fault on the outer rise north of the trench.Each tsunami scenario was tested by comparing simulations of flow depths at the observed sites of 171 brain-coral boulders and cobbles that were scattered inland some six centuries ago.We searched for the least coseismic slip that generates a tsunami large enough to flood the resting places of at least 95% of these mapped coral clasts to depths of 1 m or more.We also checked the proposed ruptures for consistency with fault scaling relationships.The comparisons of simulated tsunami flooding with observed locations of coral clasts neglects potentially large uncertainties from roughness, vegetation, and tsunami erosion.The main results are as follows: 1.The extreme storm scenario fails to account for a majority of the coral clasts in areas scarcely sheltered by modern beach ridges.This negative result is consistent with post-hurricane surveys on Anegada that documented meager erosion and deposition by 2010 Hurricane Earl (Atwater et al., 2014) and, more importantly, 2017 Hurricane Irma (Spiske et al., 2021).Taken together, these negative findings leave a tsunami as the only likely cause of the catastrophic precolonial flooding at Anegada. 2. In a leading tsunami scenario, simulated tsunami flooding covers all the mapped coral clasts to a depth of 1 m or more given 7.5 m of slip on a megathrust rupture about 255 km long centered due north of Anegada.The corresponding earthquake is of M8.45.This outcome scarcely changes, given 7.5 m of slip, if the megathrust rupture is extended eastward or westward from the limits in Figure 5a, even though earthquake magnitude increases with fault length.As for normal faulting on the outer rise, slip of 11.4 m is the smallest that enables tsunami flooding, to 1 m or more, of all the brain-coral clasts.The corresponding earthquake is of M8.17. 3. The successful megathrust and outer-rise scenarios each produce rapid flooding in the Virgin Islands.The shortest simulated arrival times, on Anegada, are 7 min for the megathrust and 13 min for the outer rise.4. For the M8.45 megathrust scenario, its 7.5 m of slip corresponds an average recurrence interval between 1,250 and 3,750 years, depending on the seismic coupling assumed.For the M8.17 outer-rise scenario, a simplified estimate gives a minimum recurrence interval of about 6,000 years.

Figure 1 .
Figure 1.Tectonic structure and historical earthquake/tsunami events in the northern Caribbean.The numbers indicate the year of the earthquake plotted in source area, queried where location of source is very uncertain.Blue circles are locations where records of the 1755 Lisbon tsunami heights are available.The light blue dash line is the 2018 Irma Hurricane track.

Figure 2 .
Figure 2. Brain-coral and limestone clasts deposited on Anegada during a sea flood during the last centuries before Columbus(Atwater et al., 2017; data in  Atwater, 2023).
plays out in Figures 7-10 for two tsunami scenarios deemed successful, one each for the outer rise and the megathrust.These examples show posited coseismic deformation (Figures 7a-7d and 9a-9d), simulated maximum water level throughout the Virgin Islands (Figures 7e and 7f and 9e, 9f), simulated flow depths on Anegada (Figures 8 and 10), and tsunami wave trains computed at synthetic tide gauges (Figure 11).

Figure 3 .
Figure 3. Posited hypotherical hurricane scenarios for simulations of hurricane-induced coastal inundation on Anegada: (a) five southeast-northwest (SE-NW) scenarios mimicking the 2010 Earl track; (b) Modeled significant wave height of the SE-NW tracks; (c) Modeled storm surge height of the SE-NW tracks; (d) Modeled peak wave period of the SE-NW tracks; (e) Five northeast-southwest (NE-SW) scenarios mimicking the 2007 Noel Track; (f) Modeled significant wave height of the NE-SW tracks; (g) Modeled storm surge height of the NE-SW tracks; (h) Modeled peak wave period of the NE-SW tracks; (i) Five southwest-northeast (SW-NE) scenarios mimicking the 1999 Lenny track.(j) Modeled significant wave height of the SW-NE tracks; (k) Modeled storm surge height of the SW-NE tracks; (l) Modeled peak wave period of the SW-NE tracks; In panels (a, e, i), the black solid line represents NOAA's best track of the base hurricane path, the red solid line represents the hypothetic track making landfall on north shore of Anegada, and the dashed lines represent the variants of the hypothetic landfall track.

Figure 4 .
Figure 4. Simulation results along the north shore of Anegada, between Windless Bight and Deep Bay, obtained from the phase-resolving Boussinesq modeling of the hypothetic hurricane scenario SE-NW 3. (a) Modeled area, where the black box represents the model domain shown in panels (c, d); (b) Offshore and nearshore time series of wave amplitudes that clearly indicates the formation of infra-gravity waves (2-3 min wave period) in the nearshore area as a result of the aggregation of shortperiod swell waves offshore; (c) Simulated hurricane-induced storm water surface elevation along the north shore of Anegada.Spiske et al. (2021) gives the most recent update on the modern storm berm.(d) Simulated hurricane-induced coastal inundation depth along the north shore of Anegada.The black and purple circles in panels (c, d) indicate locations of coral and limestone clasts, respectively.Crowds of purple circles south of Windlass Bight represent the two limestone boulder fields of Watt et al. (2012) and Buckley et al. (2012).

Figure 5 .
Figure 5. Maximum tsunami flow depths on Anegada resulting from model simulations of all near-field source scenarios.(a) Source subfaults used to configure nearfield thrust-and normal-faulting source listed in Table 2. (b) Percentages of brain-coral clasts where the maximum simulated flow depth equals or exceeds 1 and 2 m for six outer rise normal-faulting source scenarios.(c) Percentages of brain-coral clasts where the maximum simulated flow depth equals or exceeds 1 and 2 m for 19 thrustfault source scenarios.

Figure 6 .
Figure 6.(a) Thrust rupture lengths, (b) thrust rupture downdip widths, (c) outer rise normal fault lengths, and (d) outer rise downdip widths of the postulated sources evaluated using multiple scaling relationships of the earthquake source parameter.
CD, 180 km long, extends from north of Culebra to north of Anegada.The modeled slips of 5.2-7.5 m represent earthquakes of M8.24-M8.35.CDE, centered due north of Anegada, extends 255 km between the longitudes of Culebra and Anguilla.The slip and magnitude ranges are 4.2-7.5 m and M8.28-M8.45.DE tests the contribution of segment C to modeled inundation on Anegada.Segment C, 130 km long, is located north-northwest of Anegada, and its eastern limit is due north of the island.DE is the shortest of the combinations, with a rupture length of about 125 km.Modeled slips of 5.2-7.5 m yield M8.13-M8.24.CDEFGH curls into the Lesser Antilles Subduction Zone, extending nearly 770 km eastward and southward from north of Culebra to the east of Guadeloupe.Modeled slips of 5.2-7.5 m accordingly correspond to M8.66-M8.76.DEFGH, like DE, evaluates absence of rupture in segment C. DEFGH extends 640 km and parallels the eastern Puerto-Rico Trench east of Anegada.The fault slips are 5.2-7.5 m and the earthquake sizes, M8.60-M8.71.

Figure 7 .
Figure 7. Co-seismic deformation resulting from the outer rise normal fault M8.17 earthquake scenario.(a) shows the co-seismic deformation, where the solid and dashed white lines represent the contour lines (at 0.5 m interval) of co-seismic uplift and subsidence, respectively; the black solid line is the profile used to demonstrate deformation in panel (b).(c) shows the bathymetry and topography, vertically exaggerated 10 times, across the British Virgin Islands, the Puerto Rico Trench, and the outer rise region over the bending portion of the North American Plate before it subducts underneath the Caribbean Plate; (d) indicates, without vertical exaggeration, the extent of normal fault on outer rise and relative slip directions; (e) shows the maximum tsunami wave amplitude computed at the regional scale of northeast Caribbean (computed at the grid resolution 30 arc sec (∼900 m)), where the white box indicates the model domain shown in panel (f); (f) shows the maximum tsunami wave amplitude at the regional level of US Virgin Islands and British Virgin Islands (computed at the grid resolution of 6 arc sec (∼180 m)).

Figure 8 .
Figure 8. Simulated tsunami flow depths on Anegada resulting from the outer rise normal fault N M8.17 scenario.

Figure 10 .
Figure 10.Simulated tsunami flow depths on Anegada resulting from the pure thrust CDE M8.45 scenario.

Figure 9 .
Figure 9. Co-seismic deformation resulting from the pure thrust CDE M8.45 earthquake scenario.(a) shows the co-seismic deformation, where the red and blue solid lines represent the contour of co-seismic uplift and subsidence at 0.5 m interval, respectively; the black solid line is the profile used to demonstrate deformation in panel (b).(c) shows the bathymetry and topography, vertically exaggerated 10 times, across the British Virgin Islands, the Puerto Rico Trench, and the outer rise region over the bending portion of the North American Plate before it subducts underneath the Caribbean Plate; (d) indicates, without vertical exaggeration, the extent of subduction thrust fault and relative slip directions; (e) shows the maximum tsunami wave amplitude computed at the regional scale of northeast Caribbean (computed at the grid resolution 30 arc sec (∼900 m)), where the white box indicates the model domain shown in panel (f); (f) shows the maximum tsunami wave amplitude at the regional level of US.Virgin Islands and British Virgin Islands (computed at the grid resolution of 6 arc sec (∼180 m)).

Figure 11 .
Figure 11.Tsunami arrival time and waveforms.(a) Contours of the tsunami arrival time in mins from the outer rise M8.17 event.The contours are set at an interval of every 5 min.(b) Contours of the tsunami arrival time in mins from the thrust M8.45 event.The contours are set at an interval of every 5 min.(c) Time series of simulated waveforms at virtual gauges on islands in the northern Caribbean.Datum on ordinate is storm-surge level about 1.5 m above mean sea level.In panels (a, b), the small white circles indicate the locations of the virtual gauges, and the circled numbers correspond to ordering of the virtual gauges in panel (c) where the simulated waveforms are presented together with the arrival times in minutes.

Table 1
Extreme Category-5 Hurricane Scenarios Posited in Figure3All tracks are simulated with the atmospheric pressure 882 hPa and the wind speed 315 km/hr.

Table 2
Subduction Thrust (A-H) and Outer-Rise (N) Fault Segments Simulated, by Combinations Listed in Table3 Note.Segments A-H and N are mapped in Figures S3.1 in Supporting Information S1 and Figure

Table 3
Scenarios Simulated for Ruptures of Subduction Thrust (A-H) and Outer-Rise (N, NW, NE) Segments Defined in Table2