Revising evidence of hurricane strikes on Abaco Island (The Bahamas) over the last 700 years

The northern Bahamas have experienced more frequent intense-hurricane impacts than almost anywhere else in the Atlantic since 1850 CE. In 2019, category 5 (Saffir-Simpson scale) Hurricane Dorian demonstrated the destructive potential of these natural hazards. Problematically, determining whether high hurricane activity levels remained constant through time is difficult given the short observational record (< 170 years). We present a 700-year long, near-annually resolved stratigraphic record of hurricane passage near Thatchpoint Blue Hole (TPBH) on Abaco Island, The Bahamas. Using longer sediment cores (888 cm) and more reliable age-control, this study revises and temporally expands a previous study from TPBH that underestimated the sedimentation rate. TPBH records at least 13 ≥ category 2 hurricanes per century between 1500 to 1670 CE, which exceeds the 9 ≥ category 2 hurricanes per century within 50 km of TPBH since 1850 CE. The eastern United States also experienced frequent hurricanes from 1500 to 1670 CE, but frequency was depressed elsewhere in the Atlantic Ocean. This suggests that spatial heterogeneity in Atlantic hurricane activity since 1850 CE could have persisted throughout the last millennium. This heterogeneity is impacted by climatic and stochastic forcing, but additional high-resolution paleo-hurricane reconstructions are required to assess the mechanisms that impact regional variability.


SUPPLEMENTAL 1. Hurricane track density distribution map methods
The map of hurricane track density distribution was based on International Best Track Archive for Climate Stewardship (IBTrACS) v04 hurricane track segments from 1851 to 2019 CE 1,2 . Hurricane intensity (Saffir-Simpson scale) at each segment is from the IBTrACS v4 subset USA_Agency_SSHS, which is derived from HURDAT_ATL maximum wind speed data. The USA_Agency subset is a hierarchical composite of multiple meteorological organizations measured maximum wind speed data, but is primarily populated by HURDAT2 data. For more information on how data from the USA_Agency subset was compiled, see Knapp et al. 1,2 . To create this map, a grid of 50x50 km cells was overlaid across the North Atlantic Basin (including the GOM and Caribbean) in ArcMap v10.7.1 software. IBTrACS v04 storm tracks were then filtered to exclude any storm track segments that were not ≥cat. 2 (or 3 in Supplementary Fig. S1) hurricanes, meaning that segments of storms with maximum wind speeds <154 km hr -1 were not counted at any point, even if the storm reached ≥cat. 2 intensity at other points along its track. Once the data was filtered, the storm-tracks were "dissolved" (i.e., merged) into a single track for each qualifying storm using the Storm Identification number (SID) as a qualifier. This left 566 individual storms segments that qualified as ≥cat. 2 hurricanes overall. These 566 hurricane tracks were then spatially joined to the 50x50 km grid cells so that the number of segments that fell within each could be summed. These grid cells were then spatially joined to points located at the center of each cell so that the point could be assigned the z-value of the sum of ≥cat. 2 or 3 storm tracks within the grid cell.
In order to spatially characterize the gradients in storm counts throughout the basin, we utilized Inverse Distance Weighted (IDW) interpolation with the ArcMap Geostatiscal Analyst toolbox. This interpolation method allows the user to specify the degree to which surrounding point values influence one another. We utilized a power factor of 8, meaning that the IDW interpolation weighted proximal value points far more heavily that distal points. This is an appropriate method given that hurricanes strikes are relatively localized events and each grid cell, and the point with event count values were richly and evenly disturbed across the north Atlantic. Given that storm tracks are dynamic features that are geographically stochastic (albeit, not necessarily climatologically), IDW interpolation using a higher power factor of 8 will better preserve extremes and spatial count gradients in data values that are evenly spaced and relatively abundant.
The resulting interpolated output was visualized using six in event frequency bins. We acknowledge increased levels of uncertainty in the track and measured intensity of storms that occurred prior to the satellite era of hurricane observation in the 1960s, assuming that observational error was consistent from region to region over time, we can be confident in the results of this spatial analysis as a metric of relative geographic variability in hurricane track density.

Supplementary Fig. S1|
Frequency of exposure to ≥category 3 hurricane winds a 50 km radius from 1850 to 2019 CE throughout the North Atlantic. As with Fig. 1 from the primary manuscript and described above, but a density map of higher intensity >category 3 events. Storm track and intensity data is from the IBTrACS v04 1,2 dataset. Saffir-Simpson storm intensity is specifically derived from the IBTrACS v4 subset USA_Agency_SSHS, which is derived from HURDAT_ATL maximum wind speed data. Thatchpoint Blue Hole is indicated with a red star. The regions that we refer to as modern hurricane hotspots are pinkish white in color (11-16 ≥category 2 events within 50 km s). Basemaps were downloaded from DIVA-GIS 3 . Map and assciated spatial calculations were performed in ArcMap 10.7.1 software using North America Albers Equal Conic Area projected coordinate system 4 .  5 ; however, it is now thought that the is a push core collected on SCUBA in August, 2011 was collected from underneath the wall that overhangs to the east of the yellow circle as marked by the yellow arrow. Aerial imagery was collected by Pete van Hengstum in January 2015 using a DJI Phantom 3 with a mounted GoPro Hero3 digital camera. Supplementary Table S1| Radiocarbon results from Thatchpoint Blue Hole. Of the 16 samples from TPBH-C3 that were radiocarbon dated, only 11 of were included in the Bayesian statistical age model as 5 samples had measured ages older than stratigraphically deeper dates. This suggests that these samples were reworked from older sediments before being redeposited in to TPBH and are not representative of stratigraphic-age. Both dates from TPBH-C2 were included in the age-model. Radiocarbon dates were calibrated using IntCal13 or MarineCal13 calibration curves 6 as described in the METHODS in the manuscript. The text "NaN ± NaN" in the "Conventional 14 C age" column of sample 1 denotes that the sample is younger than 1950 CE, and was therefore calibrated using the Northern Hemisphere Zone 2 (NHZ2) post-comb calibration data set in CALIbomb 7 . All possible 1σ and 2σ age calibrations are provided for each submitted date.  Downcore radionuclide activity for 137 Cs and 210 Pb was measured to better constrain the modern chronology of the upper portion of TPBH-C3. 137 Cs is a manmade radionuclide that is useful for identifying the onset of hydrogen thermonuclear weapons testing in 1954 and its moratorium in 1963 13 . To quantify 137 Cs activity, bulk sediment was sampled at 5-15 cm intervals, desiccated, and powdered. Powdered samples were then placed on a Canberra GL2020RS low energy Germanium gamma well detector where the gamma decay of the sample was measured for a 24-hour period.

Index
The peak in 137 Cs activity associated with 1963 CE moratorium on nuclear weapons testing falls at 80-81 cm.
Independent of this constraint, the sedimentary age-model predicts as median age of 1963 CE at 79-80 cm (Supplementary Fig. S4). Further, the onset of detectable levels of 137 Cs around 1954 CE occurs between 115 to 100 cm in the composite TPBH record. This is within the age and depth error range of the Bayesian age-model (Fig. 3)     Event beds E10-11 were deposited in the last half of the 19th century, which corresponds with the greatest uncertainties in the sedimentary age-model during the historical period (2σ age: ±28 years, n = 62, Fig. 4a) and the observational dataset 48-51 (precision/accuracy of storm track coordinates and intensity measurements). As such, we consider the attributions made from 2014 to 1900 CE to be far more reliable for characterizing the sensitivity of the composite record. With this in mind, event bed E10 has a median age of 1883 (2σ age-range: 1909 to 1851 CE) and event bed E11 has a median age of 1873 (2σ age range: 1900 to 1881 CE). 14 storms passed within 115 km of TPBH between 1855 and 1899 CE (Fig. 4e), but only three were likely to deposit event beds. Three of these events passed about 40 km west of the TPBH on south to north trajectories on 26 August 1893, 9 September 1883, and 2 October 1866 as category 3, 2, and 4 storms, respectively. We attribute event bed E10 to the 1893 hurricane given that it was more intense than the 1883 event, and event bed E11 to the category 4 hurricane in 1866 CE. The only other intense storms to make a close passage to TPBH were three category 3 events that all passed ~20 to ~40 km northeast of the site (1896 CE, 1893 CE, and 1887 CE; Fig. 4e). However, peak winds from storms following this trajectory were likely from the north meaning that TPBH was likely shielded by Abaco Island, thus inhibiting event bed deposition.