Tropical North Atlantic sea surface temperatures and low-latitude rainfall covary, but prehistoric subtropical rainfall records are often misaligned. Here, a submarine groundwater discharge record from the northern Bahamas archives regional water balance in the northeastern Atlantic Warm Pool. We compare the reconstruction to the tropical North Atlantic seasonal temperature gradient, which can help inform how northeastern Caribbean rainy seasons are influenced by the Atlantic Warm Pool. A positive water balance in the northern Bahamas aligned with a ~0.9 °C seasonal temperature gradient from 0 to 950 CE, with both covarying on multi-decadal timescales. Aridity began at ~950 CE when a ~2.2 °C seasonal temperature gradient increase likely shortened the wet season. From 1450 to 1850 CE, frequent hurricanes offset aridity in the northeastern Caribbean by elevating rainfall. This record archives time-transgressive changes in hydroclimate forcing, and suggests that projected changes to rainfall seasonality must be considered when assessing tropical water security risk.
Mitigating future water scarcity and economic risk requires establishing the natural drivers of hydro-meteorological disasters (drought, floods, and hurricanes). However, there is spatiotemporal complexity in hydroclimate records from the lower latitudes in the Atlantic Ocean, where seasonal changes to the North Atlantic Subtropical High (NASH), Intertropical Convergence Zone, and the Atlantic Warm Pool (AWP) drive distinct wet and dry seasons. Whereas hydroclimate events like the 800-1050 Common Era (CE) drought are regionally documented, many subtropical North Atlantic paleo hydroclimate records during the CE exhibit divergent variability. At ~1000 CE, for example, the Gulf of Mexico (GoM) and tropical North Atlantic (TNA) started a ~1 °C cooling trend that culminated in ~3 °C boreal winter cooling during the Little Ice Age (LIA: 1450-1850 CE) . While observed sea surface temperature (SST)-hydroclimate relationships suggest that the subtropical North Atlantic should have become more arid, as was recorded by Cuban lagoon sediment, a Cuban speleothem and Bahamian organic geochemical proxies document wetter conditions starting at 1000 CE. The subtropical Atlantic is a key locality to understand prehistoric NASH-AWP dynamics, which is critical because future NASH expansion is projected to increase regional evaporative demand.
Hurricanes also augment regional moisture supply, with distinct forcing mechanisms as compared to other regional hydroclimate variability. While hurricanes catastrophically impact local socioeconomics, they may mitigate future Caribbean water scarcity problems through their contributions to local rainfall. Only rare paleo-hydroclimate records can distinguish between synoptic and hurricane-mediated rainfall signals. On paleo timescales, this is further complicated by recent findings that historically unprecedented hurricane activity has been spatially asynchronous in the Atlantic during the CE. Thus, large uncertainty persists in disentangling the imprints of hurricane activity on paleo hydroclimate records from this water-sensitive region.
Here, Bahamian water balance is inferred from a high-resolution record of benthic foraminiferal paleoecology from a coastal blue hole. On carbonate landscapes, underground meteoric lenses form when a positive water balance favors aquifer recharge (i.e., precipitation > evapotranspiration) as rain infiltrates into the landscape. Aquifers subsequently discharge at coastal zones as submarine groundwater discharge (SGD) , which can promote hydrographic stratification in oceanic blue holes acting as springs. SGD is a key part of the hydrologic cycle that influences coastal biogeochemistry and ecosystems. This work demonstrates that SGD has been an overlooked opportunity to explore landscape-integrated signals of excess moisture. In The Bahamas, time-transgressive changes in hydroclimate forcing means that meteorological data from the instrumental record alone is insufficient to inform on hydroclimate conditions during intervals unlike the last ~150 years.
Hydroclimate in The Bahamas
The northeastern Caribbean hydroclimate zone including Abaco Island (Fig. 1A, see supplementary information) experienced ~1297 mm/yr mean annual precipitation (MAP) from 1855 to 2017 CE (range: 603-2521 mm/yr), in near balance with 1198 mm/yr potential evapotranspiration [PET, 1950-2016 CE]. This highly sensitive water balance leaves little difference between when meteoric lenses are recharged by precipitation (P) or depleted by evapotranspiration (ET) . Whitaker and Smart estimate that 25% of MAP is available for effective aquifer recharge because most Bahamian islands have some meteoric lenses. In contrast, diminished P and aquifer recharge can shrink meteoric lenses within months. MAP at Abaco Island is partly driven by hurricane-mediated rainfall, for example, Hurricane Dorian in 2019 (a category 5 event) delivered 580 mm (45% of MAP) in just 3 days.
In this hydroclimate zone, 75% of MAP is delivered during boreal summer and mediated by AWP-NASH dynamics (Figs. S1-S3). Winter rainfall is partially linked to the El Niño/Southern Oscillation and Pacific influences (mean: 315 mm/yr, range: 206-420 mm/yr), but as this is the dry season, Atlantic processes primarily control MAP in The Bahamas (see supplemental text). As the AWP expands in the early summer and reaches the Bahamian archipelago, both NASH western boundary convergence and the wet season begin. In the late summer, the Caribbean wet season ends first in The Bahamas as cooler SSTs arrive from AWP retreat and the NASH contracts eastward. The Bahamas is an ideal location to observe this phenomenon, as its geography and wet season sensitivity to the NASH and AWP allow it to record spatial variability in prehistoric AWPs. MAP over The Bahamas is increased by larger AWP size and the timing of its arrival. For the latter, a weaker seasonal temperature gradient may indicate a longer wet season, as warmer winter and spring temperatures promote earlier AWP expansion toward The Bahamas and thus increased persistence. Indeed, from 1891 to 1970 CE, increased MAP in The Bahamas is associated with a weaker seasonal temperature gradient in the tropical Atlantic Cariaco region (Fig. 1B).
Freshwater River Blue Hole
Freshwater River Blue Hole (FRSH, Fig. 1, 26.426°N, -77.175°W) is on the leeward carbonate tidal flats of Abaco (~700 km), which have experienced 1.4 ± 0.5 m of sea-level rise during the CE. The sediment-water interface in the blue hole center is 5 m below sea level (mbsl), but caves around the periphery deepen to 50 mbsl. Seasonal water column salinity (18.2 to 36.2 psu), pH (6.7 to 8.1), and dissolved oxygen (<0.3 to 6.2 mg/l) reveal a transient meteoric cap whose vertical extent allows a local pycnocline to persist near the sediment-water interface in the blue hole middle (Fig. S4). P and SGD form the meteoric cap, and it decays by E, seawater tidal exchange, and SGD cessation. Hydrographic monitoring indicates that a meteoric cap can develop in FRSH during any season from increased P or decreased ET, and it may decay in less than a year (Figs. S4, S5). Benthic foraminifera adapted to brackish, organic-rich, and oxygen-depleted conditions dominate the modern blue hole population (e.g., Ammonia, Cribroelphidium), and not the larger miliolids commonly found on the tidal flats (e.g., Peneroplis).
Blue hole sediment cores (FRSH-C2, FRSH-C3) reveal a strikingly laminated and bedded ~230 cm carbonate mud unit preserved by benthic hypoxia during the CE. A robust radiocarbon age model developed with mangrove leaves indicates high sedimentation rates: 0.5 to 1.0 mm/yr from 0 to ~1020 CE gradually increased to 2 mm/yr from ~1020 to 2016 CE (Figs. S5-S7, Table S1). Thicker organic-rich horizons from 0 to 850 CE in FRSH-C3 (mean 35% bulk organic matter) transition to millimeter-scale laminations during the Medieval Climate Anomaly (MCA: 950 to 1450 CE), and thicker organic-rich beds reappear from 1450 CE to present. The organic-rich horizons are dominated by fine-grained organic particles, likely degraded products of aquatic primary productivity. Hurricane Dorian (category 5, Saffir-Simpson Scale) had negligible impact on regional carbonate flat geomorphology (Figs. S8, S9), which is like observations elsewhere.
The water balance index (WBI)
Subfossil benthic foraminifera can be organized into a WBI, which is the arithmetic sum of the proportion of 6 foraminifera taxa with a similar ecological behavior. Their changes in relative abundance and diversity reflects a well-established ecological response to environmental conditions at the sediment-water interface, and these conditions are regulated by SGD (Fig. 2B). First, cluster analysis grouped benthic foraminiferal samples into communities with abundant smaller miliolids (groups 1-2). The WBI includes samples with an increasing proportion of 6 taxonomic units (groups 3-8): Ammonia beccarii, Cribroelphidium poeyanum, Cribroelphidium gunteri, Pseudoeponides anderseni, Psuedoeponides davescottensis, and Trichohyalus aguayoi (Figs. S10, S11). These 6 taxa tolerate lower dissolved oxygen, lower salinity, and high organic matter (see supplementary text). For example, A. beccarii dominates group 8 (mean 79.4%) relative to miliolids (mean 16.7%), whereas groups 1-2 contain mostly small Quinqueloculina and Triloculina (mean 81.6%). Miliolid taxa would favor a more negative water balance (reduced P or greater ET) from reduced stratification, increased benthic ventilation, and microbial organic matter remineralization. Winter cooling and convection during more negative water balance conditions would further favor benthic ventilation and miliolid populations. Since ET is highly correlated to atmospheric temperature, regional warming would make FRSH more favorable to miliolids by increasing evaporative demand.
Ordination and isotopic analysis introduce further confidence that the WBI is reflecting SGD and water column stratification. Detrended correspondence analysis (DCA, Fig. S10) Axis 1 explained 46.8% of the data variance, and samples have a higher DCA Axis 1 value if they have higher proportions of the 6 taxa grouped into the WBI (r = 0.95, p < 0.001). Bulk organic matter was also correlated to DCA Axis 1 (r = 0.53, p < 0.001), and the WBI was independently correlated to bulk organic matter (r = 0.56, p < 0.001). An increasing (decreasing) WBI, increasing (decreasing) organic matter preservation, and decreasing (increasing) miliolids are significantly linked. An alternate hypothesis would state that during periods of less SGD, more evaporation and local aridity would vertically contract the meteoric lens, and in turn increase benthic hypoxia, organic matter preservation, and bathing of the sediment-water interface with more saline groundwater. However, species-specific stable oxygen (δO) and carbon (δC) isotopic analysis indicates the opposite: WBI taxa secreted their tests in seawater more diluted by meteoric water from SGD (Fig. S12, see supplementary text). In a last millennium time series analysis, the difference in δO values on miliolids vs. WBI taxa (ΔO) indicates more frequent SGD of meteoric water during the LIA than MCA. In short, higher (lower) WBI values indicate when a more (less) positive water balance on the Abaco landscape enhanced (depressed) SGD and promoted (hampered) hydrographic stratification in FRSH.
The WBI is limited by not having a direct calibration to SGD, but only indirectly through an assessment of local changes in rainfall and temperature during the instrumental and historical period. These comparisons indicate that the WBI is not a simple proxy for precipitation, temperature, or evaporation alone, but the WBI is most likely reflecting the SGD of meteoric water in response to a positive water balance on the landscape. For example, a notable WBI decrease from ~0.7 to 0.2 at ~1960 CE (Fig. 2B) is within radiocarbon uncertainty of when Bahamian MAT increased from 24.3 ± 0.4 °C (1855 to 1970 CE) to 25.4 ± 0.4 °C (1971 to 2014 CE, Fig. 1C), while MAP stayed relatively stable. Increasing air temperatures (both summer and winter, Fig. S1) and local ET (Fig. 1C) thus likely suppressed SGD to create conditions more favorable to miliolids and resulted in a lower WBI value in years younger than ~1960 CE. In contrast, precipitation alone does not drive foraminiferal diversity. NASH expansion from 1910 to 1940 CE suppressed the wet season and caused a MAP decrease to ~800 mm/yr (Figs. 1C, and S1). While sedimentation rates and radiocarbon uncertainties are sufficient to resolve this event, the WBI does not respond to this multi-decadal MAP change. This could be because of compensating impacts from relatively low MAT during this period helped limit ET and maintain the water balance. A direct calibration of the water balance index is further limited because local hurricane activity during the last 150 years is not representative of the full potential of hurricane-mediated rainfall on the landscape based on multi-centennial reconstructions. Direct monitoring of SGD through FRSH is an important opportunity for creating more quantified relationships between local water balance, SGD and benthic biodiversity and geochemical changes. Nevertheless, benthic foraminifera are highly sensitive environmental indicators that would respond to environmental changes at the sediment-water interface 5 msbl, which is most likely caused by changes in the local water balance altering SGD.
Each WBI datapoint reflects a multi-year average (5-20 years) of environmental conditions at each core site. This is because (i) each foraminiferal sample averages multiple years (0 to 1020 CE: 10-20 years/sample, 1020 CE to present: ~5 years/sample), (ii) the vertical water column structure in FRSH is seasonally unstable, and (iii) benthic foraminifera can exhibit boom-bust population changes in <1 year. A single environmental parameter cannot explain foraminiferal population changes because WBI taxa all tolerate low salinity, hypoxia, and abundant organic matter. When conservatively used the WBI provides a detailed perspective of excess moisture supply and SGD (Fig. 2B, see supplementary text).
Time-transgressive hydroclimate forcing
Instrumental data, paleo studies, and models routinely link circum-Caribbean precipitation variability to tropical Atlantic SSTs, but nowhere is moisture supply more sensitive to the AWP than The Bahamian archipelago in the northeastern Caribbean. During the instrumental record, the TNA seasonal temperature gradient is associated with MAP in the northeastern Caribbean, mostly likely through the timing of AWP arrival, its persistence, and the associated AWP impact on wet season length (Fig. 1B). The only available record to derive a tropical seasonal temperature gradient proxy is the seasonally-resolved SST record from the Cariaco Basin (Fig. 2D). Whereas the Cariaco Basin is ~2100 km southeast of The Bahamas, its SST is highly correlated to broader TNA (6° to 22° N, -75° to -15° W) SST anomalies (r = 0.94, p < 0.001) and to AWP extent (r = 0.89, p < 0.001). Like others, we also infer Mg/Ca ratios in the planktic Globigerina bulloides as a proxy for winter/spring SSTs versus G. ruber to infer summer/fall SSTs (Fig. 2E, F). Indeed, G. bulloides can re-bloom in the late summer and it has a slightly deeper habitat (~30 m depth). However, differencing these signals provides the only annually resolved proxy for the tropical Atlantic seasonal temperature gradient (Fig. 2D). A higher (lower) seasonal temperature gradient in the Cariaco Basin suggests a shorter (longer) wet season in the northern Bahamas.
During the first millennium of the CE, the seasonal temperature gradient averages 0.9 ± 0.5 °C, and the WBI is qualitatively more positive than during the second millennium. However, there is evidence for less SGD at 125 CE, 425 CE, 625, and 800 CE synchronous with multi-decadal increases (>1 °C) in the seasonal temperature gradient (Fig. 2D). Within radiocarbon uncertainty, colder conditions were also detected in lower-resolution SST reconstructions at 450, 600, and 850 CE in the southeastern Caribbean and northern GoM (Fig. S14). In contrast, periods with the most positive water balances (centered at 50 CE, 275 CE, 525 CE, and 725 CE) coincide with reduced (-1 to 1 °C) seasonal temperature gradients, which suggests longer wet seasons in the northeastern Caribbean. The close visual covariance between Abaco hydroclimate and the Cariaco Basin seasonal temperature gradient introduces further confidence in the Cariaco Basin as a robust reconstruction of seasonal SSTs (Fig. 2E, F).
From ~950-1450 CE, lower WBI values indicate less SGD and water column stratification. Sedimentation rate and adult-size miliolid foraminifera increased (e.g., Quinqueloculina tenagos), which is consistent with aridity decreasing SGD and promoting abiotic carbonate deposition. Evidence for AWP changes at ~1000 CE was previously documented in the northern GoM based on (i) decreasing Globigerinoides sacculifer abundance, a planktic foraminiferal proxy for Loop Current penetration, and (ii) cooler SSTs based on Mg/Ca-δO measurements of planktic foraminifera. This suggests that a change in ocean circulation also occurred at 1000 CE, which may have contributed to seasonal temperature gradient changes. Nearby Cuban lagoons document a synchronous aridity trend that deepens through the LIA, along with Caribbean-wide droughts from 800 to 1200 CE. The MCA experienced ~1 °C atmospheric warming in the TNA, which would have amplified local evaporative demand similar to the period from 1960 CE to present (Fig. 1C). While summer/fall SSTs are greatest during the MCA (~28.75 °C, Fig. 2E), the sharp winter/spring SSTs cooling at ~1000 CE markedly increased the seasonal temperature gradient to 2.2 ± 0.4 °C during the second millennium of the CE (Fig. 2D). If correct, this should have caused a concomitant change in the AWP, as supported by lower WBI values.
Counterintuitively, WBI values indicate greater moisture supply from 1425 to 1625 CE (Fig. 2B) when the southeastern Caribbean and northern GoM are colder, and there is a greater seasonal temperature gradient (Figs. 2D, and S14). We hypothesize that the greatest increase in hurricane activity observed over the last 1500 years in the northern Bahamas compensated local water balance (Fig. 2A). Specifically, Abaco Island transitioned from ~6 to ~13 hurricanes per century by 1500 CE as observed elsewhere at Thatchpoint Blue Hole (Fig. S11). On the other hand, the WBI decreased from ~1650 to 1725 CE along with local hurricane activity. From 1750 to 1850 CE, the seasonal temperature gradient exhibits a largely secular decreasing trend (Fig. 2D), and Bahamian hurricane activity increased (Fig. 2A, C), which together likely favored increased local moisture supply (Fig. 2B). While there were fewer local hurricanes during the 1960s-1980s CE, the onset of modern Caribbean warming and aridity (Fig. 1C) likely caused the recent WBI decrease at ~1960 CE, which is strikingly like the WBI observed during the MCA.
Assessing future hydroclimate in the context of CE paleoclimate
Coarsely resolved state-of-the-art climate models cannot realistically simulate highly regional hydroclimate processes, so projections of current and future northern Caribbean hydroclimate are highly uncertain. Nevertheless, current and future seasonal temperature gradient changes are important to quantify given (i) long-term linkages between the Cariaco region seasonal temperature gradient and northern Caribbean hydroclimate demonstrated herein, and (ii) the large variability in the seasonal temperature gradient over the CE (Fig. 2D). Curiously, the only two simulations spanning the CE do not exhibit large seasonal temperature gradient variability like the Cariaco record (Fig. 3A). Nevertheless, both observations (-0.53) and historical simulations (ensemble mean of -0.25) from the Coupled Model Intercomparison Project phase 6 (CMIP6) suggest that the seasonal temperature gradient has already increased from anthropogenic forcing (88% of simulations show an increase: Fig. 3B). The seasonal temperature gradient is projected to further increase under both a lower (SSP2-4.5) and higher (SSP5-8.5) emissions scenario, with an ensemble mean change of -0.39 (66% of simulations show an increase) and -0.46 (63% of simulations show an increase) by the end of the 21st century, respectively (Fig. 3B). In two of the SSP5-8.5 simulations the projected seasonal temperature gradient change is larger than the difference between the first and second millennium of the CE (Fig. 3A -- fifteen simulations exceed the 25th percentile of the differences between all 30-year periods in the first and second millennium of the CE). While considerable winter SST cooling increased the seasonal temperature gradient in the second millennium of the CE, summer SSTs are warming faster than their winter counterparts (not shown) in both model simulations and recent observations. Despite these differences, the large projected changes to the seasonal temperature gradient suggest the potential for even starker drying in the northern Caribbean by the end of the 21st century, perhaps exceeding the reconstructed aridity observed (i) during the MCA and (ii) since 1960 CE by the WBI Index (Fig. 2B). A future research priority should contextualize projected changes to the seasonal temperature gradient in terms of wet season duration and AWP-NASH dynamics. It remains a challenge to quantify first to second millennium climate changes because most CE simulations begin in 850 CE. Understanding the MCA onset and using it as a target for climate model validation may require past2k simulations like those used herein.