Three‐dimensional stratigraphic architecture, secondary pore system development and the Middle Pleistocene transition, Sandy Point area, San Salvador Island, Bahamas

The Bahamian Archipelago has been the subject of extensive studies of stratigraphy, carbonate island morphology and diagenetic overprinting for many decades. A recently‐acquired dataset comprised of high‐resolution light detection and ranging, tightly spaced cored boreholes with image logs, thin sections, porosity and permeability measurements, and isotopic geochemical analyses provides a unique opportunity to perform detailed sequence stratigraphic analysis. This multi‐faceted study first establishes a new high‐resolution sequence stratigraphic framework for the Sandy Point area of San Salvador Island, encompassing six unconformity bound sequences (exposure surfaces) representing varying lengths of missing time. The stratigraphically lowest sequence is Late Pliocene in age and is dominated by reef facies capped by an extensive laminar caliche that represents 1.5 million years of exposure. The overlying Early Pleistocene is split into two sequences, EP1 and EP2, which are both dominated by low‐energy subtidal facies. The upper three sequences are tied to the Marine Isotope Stage 11, 9 and 5e interglacials, and are distinctly different in facies composition and architecture, being dominated by higher energy facies distributed in a more complex mosaic. The new stratigraphic framework highlights significant changes in depositional style and facies architecture that can be linked to increasing sea‐level amplitude oscillations and greater climate gradient at the Mid Pleistocene Transition. The well‐constrained framework is then used to provide key context for observed complex diagenetic evolution through repeated sea‐level inundation and subaerial exposure. These results suggest that well‐developed subaerial exposure surfaces drive increases in non‐matrix features in the subsequent stratigraphic package following exposure as the well‐cemented exposure surface acts to concentrate fluid flow and dissolution. Integration of non‐matrix features within this framework highlights the importance of dissolution/cementation patterns on formation of these modified unconformity surfaces for focusing later meteoric diagenesis and creating enhanced fluid flow pathways for continued multicyclic diagenetic events.


INTRODUCTION
The Pliocene through to Pleistocene stratigraphic record of the Bahamian archipelago records complex interactions between eustatic fluctuations related to ice-house climatic conditions, sediment production and depositional patterns, and subsidence/accommodation signals (or lack thereof).
There have been extensive studies on sedimentation of Bahamian carbonates (e.g.Shinn et al., 1969;Lynts, 1970;Hine & Neumann, 1977;Mullins & Lynts, 1977;Harris, 1979;Rankey & Reeder, 2010;Harris et al., 2015;Purkis & Harris, 2016) and a few lesser known contributions on the surface exposures (Garrett & Gould, 1984;Hearty & Kindler, 1993, 1995;Carew et al., 1995;Mylroie & Carew, 1995) that include mainly strata of the major interglacials between the Holocene (Marine Isotope Stage 1 -MIS 1) and 410 ka, including MIS 5e, 9 and 11, and possibly MIS 5a.What is lacking from this suite of modern analogue observations is a record of how these sediments are transformed into a fully three-dimensional stratigraphic record of a scale and complexity similar to that observed in the deep-time record.
Key questions still exist regarding: (i) how highly topographically/geomorphically complex interglacial deposits transform into the more 'layer-cake' strata like those observed in ancient records; (ii) how sediment accumulation, erosion and dissolution interact to create the final or 'nearfinal' diagenetically stable systems; and (iii) what is the impact of the alternation between the brief highstand depositional record (10-15 kyr) versus the extended lowstand phase of dissolution and mineralogical stabilization, and does this prolonged period of exposure and dissolution completely overprint the depositional record or simply modify or enhance the depositional units.
A critical role of outcrop analogue studies for guiding subsurface reservoir characterization is to illustrate true scales of lateral and vertical heterogeneity that can commonly only be speculated away from well control and below seismic resolution.Previous studies of Pleistocene-Pliocene subsurface stratigraphy of the Bahamas (Beach & Ginsburg, 1980;Pierson, 1982;Williams, 1985;Eberli & Ginsburg, 1987;Budd & Land, 1990;Whitaker & Smart, 1990, 1994;Aurell et al., 1995;Beach, 1995;Kievman, 1998;Melim, 1996;Whitaker & Smart, 1997a,b, 1998, 2007;Whitaker, 1998;Ginsburg, 2001;Melim et al., 2001;Cunningham et al., 2004;Eberli, 2013) have well spacings of tens of kilometres and were largely targeted at delineating the gross stratigraphic picture and diagenetic trends of the Bahamian carbonate platforms.This study uses tightly spaced (<500 m) cores within a relatively small area (15 sq km), paired with a newly collected suite of airborne light detection and ranging (LiDAR), petrographic thin sections, optical and acoustic image logs, and plug and whole-core porosity-permeability data to better characterize the evolution from Holocene sediments to 'ancient' Pleistocene-Pliocene strata.Age dating constraints were limited, but a combination of Sr 87/86 profiles, rare index fossils, and correlation with the palaeomagnetic work of McNeill et al. (1988) allow the timing of deposition and hiatal intervals across the cored interval to be reasonably interpreted, certainly by the standards of older Phanerozoic records.This study primarily helps to clarify the scale and complexity of primary facies heterogeneity on an icehouse isolated platform, further highlighting the important difference between the Middle Pleistocene and Late Pleistocene icehouse patterns in terms of stratigraphic and diagenetic complexity.These results highlight that carbonate strata that are within the first 2 to 3 Myr of evolution and at burial depths that do not exceed 100 m, undergo a dramatic transformation from unconsolidated, mineralogically diverse sediments to welllithified, mineralogically stable, locally strongly leached and dissolved products that bear close resemblance to strata orders of magnitude older.A critical message from studying the diagenetic and stratigraphic transformations that occur in these sediments is that the vast majority of their alteration occurs very early in their history.The detailed characterization of touching vug/cavernous pore networks within this tightly constrained stratigraphic framework also allows for a reevaluation of the nature and timing of dissolution and pore network evolution, suggesting that Ó 2023 ExxonMobil.Sedimentology published by John Wiley & Sons Ltd on behalf of International Association of Sedimentologists, Sedimentology, 71, 207-240 lithological controls set up by subaerial exposure surfaces are the controlling mechanism, rather than dynamic geochemical interfaces such as suggested by mixing zone and flank margin cave models.These data and results point to the importance of the Pleistocene stratigraphic record as a guide for understanding fundamental processes of the evolution of modern carbonate sedimentary sequences into 'ancient' carbonate strata.
The present-day Archipelago is draped by Holocene sediments, but the island core mostly consists of carbonates deposited during Marine Isotope Stages 11 and younger, representing the high frequency, high amplitude glacioeustatic fluctuations associated with icehouse climatic conditions (Hearty & Kindler, 1995).Deposition occurred during the highstand portion of the sealevel curve (Fig. 2) (Goldhammer et al., 1990;Wright, 1992;Hearty & Kindler, 1995) when the flat-topped platforms were flooded, representing an average duration of 6 to 10 kyr (Fig. 2).Periods of deposition and exposure/erosion are marked by Marine Isotope Stages with greater numerical values corresponding to increasing age.Odd numbers represent sea-level highstands, starting in the Holocene (MIS 1).The Late to Middle Pleistocene are represented by MISs 3 to 25 (Emiliani, 1955;Shackleton, 1995;Lisiecki & Raymo, 2005;Railsback et al., 2015).
Subsidence rates for the Bahamas vary regionally but the Bahamian archipelago has been tectonically quiescent during the Pleistocene, permitting a detailed rock record of the Plio-Pleistocene aged stratigraphy.Subsidence rates have been inferred to be 1 to 2 m/100 kyr for the Three-dimensional stratigraphic architecture 209 last 100 Myr (McNeill et al., 1988;McNeill, 2005).However, these subsidence rates should be used with caution because there is little data to verify subsidence rates for the last 125 kyr, which makes up 0.1% of such record (Hearty, 1998).
There have been numerous studies that focus on various Bahamian islands, breaking out detailed stratigraphic frameworks that most commonly focus on the deposition during the highstand portion of sea-level rises.Multiple nomenclature schemes have been developed for the early to Late Pleistocene deposits common to the island chains.Early work on the subsurface stratigraphic framework of the Bahamian archipelago proposed a broad grouping of Pleistocene strata into the Lucayan Limestone (see Supko, 1970Supko, , 1977;;Beach & Ginsburg, 1980;Williams, 1985;Mullen, 1993;Caracuel et al., 1996), which primarily consisted of oolitic facies.Work on San Salvador Island (Carew and Mylroie, 1985, 1987, 1991) proposed three broad groupings for the Bahamas stratigraphy, a Holocene unit (Rice Bay Formation), a Late Pleistocene unit (Grotto Beach Formation) and a Middle Pleistocene unit (Owl's Hole Formation).This stratigraphic framework was later revised by Hearty & Kindler (1993) to include eight units extending from the Middle Pleistocene to Late Holocene.Most Recently Kerans et al. (2019b) have investigated the Late Pleistocene units providing an updated stratigraphic framework from 129 to 11.7 kyr.The various revisions and divisions of the stratigraphic nomenclature for the Bahamas that have occurred throughout the years are collated into a stratigraphic column for the Island of San Salvador, which will be used within this study (Fig. 3).

San Salvador Island
The island of San Salvador is an isolated carbonate platform, located within the Bahamian Archipelago between 24°11 0 and 23°55 0 North latitude and 74°29 0 and 74°31 0 East longitude (Fig. 1).It is located on the eastern edge of the Bahamas island chains, experiencing easterly trade winds (Neuman et al., 1978;Wanless et al., 1988).San Salvador Island is smaller than other large islands of the Bahamas such as Andros or Great Abaco, covering 209 km 2 .Much of the island is comprised of arcuate dune ridge complexes and intervening lakes.The island is surrounded by a shallow (1 to 15 m), narrow (1 to 6.5 km wide) shelf.Its margins are steep walled and display common failure scarps.Beyond the shelf edge, water depths reach 4500 m on the San Salvador Island Abyssal Plain (Mulder et al., 2019).
The Bahamas has long been studied as type localities for the investigation of early meteoric dissolution and cementation processes (e.g.Matthews, 1968;Gavish & Friedman, 1969;Land, 1970;Steinen & Matthews, 1973;Steinen, 1974;Beach, 1982;James & Choquette, 1984).Meteoric diagenesis is most commonly thought to be dominated by dissolutional processes; however, cementation is also common in this realm (James, 1972;James & Choquette, 1984;Boardman et al., 1995).Subaerial exposure surfaces (commonly referred to as caliche crusts, terra rossa, calcretes or palaeosols) are observed widely across the Bahamian Archipelago, and are attributed to sea-level lowstands associated with glacial maxima mainly of stages 6, 8, 10 and 12, as well as older stages that do not crop out (Vacher et al., 1989(Vacher et al., , 1992;;Boardman et al., 1995), and typically results in a zone of increased cementation from meteoric processes.On San Salvador Island, cementation processes during sealevel lowstand are evidenced in the extensive development of subaerial exposure surfaces (Boardman et al., 1995;Carew & Mylroie, 1995;Hearty, 1998).As pioneered in numerous outcrop and core studies, the presence of these surfaces across the island have been pivotal in helping to bracket stratigraphic units of similar age (Hearty & Kindler, 1993;Carew & Mylroie, 1995).
In addition to the presence of widespread subaerial exposure surfaces on the island, there is extensive evidence of larger scale karstification processes at play.Blue holes, caves, banana holes, pit caves and touching vug porosity are common across the island.San Salvador Island provides an opportunity to better understand the processes behind the formation of such karst related features and the hydrological conditions that drive their formation (Supko, 1970(Supko, , 1977;;Mylroie & Carew, 1990, 1995;Sampson & Guilbeault, 2013).
The outcrop or surface stratigraphy of San Salvador Island has been studied by numerous individuals, including White & Curran (1988), Hearty & Kindler (1993), Carew & Mylroie (1995) and Kerans et al. (2022;Fig. 3).Supko (1970), performed a detailed study of a core drilled on the north side of San Salvador island, to the depth of 168 m; this work focused on the description of the depositional and diagenetic features present in the subsurface, focusing on using geochemistry to identify mineralogy and date the various deposits in the core, with particular focus paid to the presence and mechanism of dolomitization.
Following the work of Supko (1970), McNeill et al. (1988) performed a detailed magnetostratigraphic dating study using 91 m of the same core.Using a total of 136 samples collected from the 91 m, McNeill et al. (1988) identified magnetic polarity consistent with the Late Pliocene through to the Late Pleistocene-Holocene.The sequence of magnetic reversals documented within the Supko core reveal six chronostratigraphic zones.This refinement highlights lithological and faunal changes that are coincident with magnetic reversal, with the disappearance of corals in the Late Pliocene and a shift from a skeletal dominated lithology to a non-skeletal dominated lithology in the Pleistocene.Further work by McNeill et al. (2001) showed the same trends, giving credence to the notion that this is not just a local occurrence.Fig. 3.A stratigraphic column for the Island of San Salvador that integrates both surface and subsurface terminology (Hearty & Kindler, 1993;Carew & Mylroie, 1995;Hearty, 1998;Godefroid, 2011;Kerans et al., 2019b).Column on the furthest left designates age and respective Epoch, followed by the various Marine Isotope Stages, approximate ages ranges, and finally the stratigraphic unit's nomenclature.The time between the deposition of each stratigraphic package is highlighted between the individual rows.The present study builds off of the surface stratigraphy established by previous authors and provides a detailed subsurface stratigraphic framework with laterally spaced cores from 100 to 500 m (on average) to constrain the complex unconformity-bound icehouse sequences and their internal facies architecture.Using the combination of detailed surface observations and core facies descriptions, it has been possible to establish, for the first time on the island of San Salvador, the stratigraphic architecture and palaeo-sea-level position spanning not only the MIS 5e but also the MIS 9 and 11 units, as well as the older Pleistocene intervals of less distinct age position.

FIELD DATA, ANALYTICAL METHODS AND CORE DESCRIPTIONS
As previously mentioned, the subsurface stratal architecture of the Plio-Pleistocene for the Sandy Point area of San Salvador Island (Figs 1 and 4) was reconstructed using a combination airborne LiDAR, of core descriptions from newly drilled boreholes, thin sections, image logs and geochemistry.

Drilling techniques
Twenty-two boreholes, 73 mm in diameter, were drilled over three field seasons using a Hydracore Prospector© diamond core drill (Hydracore Drills, Delta, BC, Canada; www.hydracore.com) and size NTW drill rods.Of the 22 boreholes drilled, CL-13, CL-17 and CL-18 were abandoned in the shallow subsurface due to poor hole conditions and, thus, little core was recovered in these locations (Table 1).Drilling of the other 19 locations produced cores with a diameter of 54 mm to depths between 14.4 m (CL-13) and 53.2 m (CL-16) in 1.5 m intervals.Lost core recovery in the shallow subsurface (ca 3 to 6 m) occurred due to obliteration of poorly cemented bedrock during drilling, while core losses in the deeper subsurface were indicative of caves and touching vugs.Void apertures were estimated in the field from bit drops, lost core recovery, and observations of vuggy and recrystallized core (e.g.Breithaupt et al., 2021;Fig. 5).Sixteen of the wells drilled in Sandy Point experienced bit drops/voids as identified by the driller, and all 22 wells experienced lost returns (see Breithaupt et al., 2021, for more details on core recovery, lost returns and bit drops).In the laboratory, core recovery and gaps were reconstructed using image logs where available to better assess drilling-related losses versus dissolution-created voids (caves).Both image logs and the optical borehole imager were critical for assessment of vuggy and cavernous porosity, which was performed in a similar manner to the detailed work of Cunningham et al. (2004).11).(B) Bare Earth LiDAR images for the same locale as the latter.Well locations are marked as well as the mapped surface geology (Kerans et al., 2019a).

Borehole images
Slimhole image logging probes were used to characterize the open macropore system at the wellbore wall because they provide continuous and oriented 360°views of the wellbore wall.On each wellbore at least one complete run was collected using an Acoustic Borehole Imager (ABI™) or Optical Borehole Imager (OBI™) manufactured by Advanced Logic Technology (ALT SA, Redange, Luxembourg).
The OBI uses a ring of lights to illuminate the borehole and a camera with a conical reflector housed in a transparent cylindrical window to generate a 360°photograph of the wellbore wall.OBI logging speed was 0.5 m/min to ensure quality of the optical image.The OBI run was typically collected before any ABI to maximize fluid clarity that could be reduced by particles put in suspension during logging operations.Light intensity was adjusted between boreholes depending on the colour of the rocks to maximize image quality.OBI images can be collected in air or clear fluid filled wellbores and, thus, were acquired below and above the water table.The quality of the OBI was compromised in certain intervals due to the presence of unflushed drilling mud that produced a coating on the borehole walls.
The ABI uses an ultrasonic pulse-echo configuration, where the amplitude and transit time of the reflected acoustic signal are recorded and used to generate azimuthal amplitude and travel-time images of the wellbore wall.Travel time can also be used to generate a calliper log.ABI logging speeds were less than 1 m/min and data was collected within one year of drilling the well.This logging speed is slower than the 1 to 3 m/min reported by Williams & Johnson (2004) and is slower than most conventional geophysical logging, resulting in higher resolution images collecting a vertical measurement every 0.002 mm.ABI images were collected in water-filled intervals of the boreholes (i.e.below water table) because the ABI does not function in air filled boreholes.Borehole enlargements related to open pores, such as fractures or caves, scatter energy from the acoustic beam, reducing the signal amplitude and producing recognizable features on the images (Paillet et al., 1990).Acoustic impedance contrast between the borehole fluid and the well is indicative of the relative density of the borehole wall.Filled secondary porosity features might be detected even when there is no change in borehole diameter if there is sufficient acoustic contrast.
The primary use of the OBI and ABI tools was to identify large pores, like open caves, that otherwise cannot be cored or might result in a heavily degraded core quality.This study has defined three classes of secondary porosity that can be interpreted from image logs: (i) caves; (ii) voids; and (iii) touching vugs (Fig. 5).These secondary porosity features were classified as nonmatrix features.Where non-matrix features are defined as being non-fabric selective, typically much larger than the matrix (depositional) pores, such as fractures, vugs and caves.This also includes fabric selective (mouldic porosity) that has been dissolutionally enhanced/ enlarged.These non-matrix features are typically distributed randomly, where local mechanical and chemical processes play a role in their Touching vugs refer to intervals with a very dense network of 0.5 to 5.0 cm 2 vugs that in three dimensions are likely connected and often are responsible for losses of wellbore fluids while drilling.Other non-matrix features identifiable in borehole images that were described in core include breccias and speleothems that are interpreted to correspond to cave collapse or infill materials (Fig. 5).Borehole image logs were depth-aligned with core to integrate both datasets and improve interpretation of image logs where there was no core control (Fernández-Ibá ñez et al., 2018).Exposure surfaces are identifiable in optical images due to their low porosity and reddish colours.Also, in acoustic amplitude images, such features appear as high acoustic impedance events representative of their extreme increase in density.

Isotope geochemistry
A single core (CL-16) was sampled for carbon, oxygen and strontium isotope geochemistry.Twenty-four individual points across the length of the CL-16 core were sampled.100 grams of material were collected, via Dremel (hand-held drill) and sent to Rutgers University for analysis (Fig. 6).
The Sr was separated from the matrix in water and rock samples by ion chromatography using 80 μl of Sr-spec resin (Eichrom), 2 M nitric acid to remove major elements, and 7 M to remove Ba in rock samples, and thus recovering Sr in water.Isotopes were run on ThermoScientific Neptune Plus MCICPMS (Thermo Scientific, Waltham, MA, USA).Measurement of the standard NBS 987 was included with all analyses, and average 0.710245 AE 0.000006 (2σ, n = 40) by TIMS or 0.710270 AE 0.000004 (2σ, n = 13) by MC-ICPMS.Analyses by MC-ICPMS were then corrected based on those values of NBS 987 values to match those of the TIMS instrument.
The insights gained from CL-16 isotope geochemical profiles were combined with core and log data to calibrate the stratigraphic framework for the Sandy Point area (Table 2).Local age data obtained from the Sr 87/86 ratios were compared and matched with the global signal which was essential for deriving a relative age.However, there is likely some error present in the deepest portions of geochemical profiles due to the pervasive dolomitization present.Irrespective of the diagenetic overprint from dolomitization, the correlation of the relative age dates in conjunction with the presence of Stylophora affinis, a coral species found in Pliocene, aid in placing our various stratigraphic units within relative time.

Core description and thin section analysis
Detailed facies descriptions were conducted for all cores in Sandy Point.Cores were described at a centimetre-scale using a standard Dunham (1969) scheme, and a combination of acid-etched core slab and thin section analysis helped with grain identification and recognition of diagenetic overprints.Core descriptions were depth corrected via distinguishing features present in image logs.Optical, amplitude and travel time derived image logs were used to interpret and classify sedimentary structures and non-matrix features.Figure 7 details the suite of data used to make detailed core descriptions.When present, the optical image log, wrapped using the calliper log, is particularly useful in identifying regions where there are well bore excursions, indicating possible open non-matrix features and their corresponding depth in the core.In addition to stacking patterns and facies association, particular attention was paid to the presence of subaerial exposure surfaces, non-matrix features and related diagenetic overprinting, evidence of cavernous porosity (for example, speleothems) and brecciation (Fig. 8).Petrographic analyses of more than 1000 Alizarin red S stained thin sections were completed to aid in the detailed core descriptions.Core descriptions were compiled into two regional crosssections to aid in interpreting stratigraphic relationships in this highly heterogeneous locale.

Methodology of cross-section construction
Nine lithofacies were present in core (Fig. 9; Table 3) from the six unconformity-bound stratigraphic units extending to a depth of approximately 30 m below surface and were compiled into two regional cross-sections (Figs 10 and 11).The cross-sections (Figs 10 and 11) are referenced to a zero mean sea-level (msl) datum.Ground surface elevations were derived from RTK-GPS mapping and compared with elevations derived from the island-wide airborne LiDAR map with agreement being within 10 cm from the two approaches.The uppermost surface in both cross-sections is the bare-earth ground surface derived from airborne LiDAR.Below 12 m water depth, multibeam imagery guides the youngest surface expression, that extends the profile to 125 m.Cross-sections are Three-dimensional stratigraphic architecture 215 Fig. 6.The CL-16 core description is provided here alongside vertical profiles of Sr 87/86 data as well as carbon and oxygen isotopic profiles.The stable isotope profile for CL-16 is of too low resolution to define individual exposure surfaces, and the shift to heavier values in the Pliocene unit reflects the dolomitized character of the unit rather than either influence of meteoric fluids or a secular trend in ocean chemistry.In contrast, the strong shift to more radiogenic values across the Pliocene-Pleistocene contact is consistent with well-documented secular trends in ocean chemistry for this time period (see Banner, 2004, for summary).*Supko core strontium values provided for reference (Swart et al., 1987).
purposefully constructed to intersect as many well locations as possible.The cross-sections will be described in more detail in the following section.

SANDY POINT STRATIGRAPHY
The Sandy Point area (Fig. 4) on the far southwest tip of San Salvador Island was selected for detailed stratigraphic and diagenetic modelling for several reasons.First, it provided a relatively small (15 sq km) self-contained region that contained all elements of the surface facies scheme (subtidal shelf, upper shoreface, foreshore and dune ridge).Second, the area is known to have abundant karst features such as flank margin caves (Mylroie & Carew, 1990;Mylroie et al., 2020), blue holes (Sampson & Guilbeault, 2013;Breithaupt et al., 2022) and pit caves such as Owl's Hole that were of direct interest to the project.Third, both at Owl's Hole and Watling's Quarry, exposures of pre-MIS 5e strata had been documented, demonstrating a reasonably complete sampling of known Pleistocene stratigraphy.Finally, a network of roads through the area provided ready access to drill sites across the region allowing for a detailed 3D assessment of the subsurface setting.
Geological mapping of the surface facies of MIS 5e affinity (Fig. 4B) illustrates a 270-degree arrangement of platform interior and shoreline facies.Starting with a relatively low-relief interior lagoon developed behind an aeolian ridge, facies transitions seaward in all directions into foreshore, upper shoreface, and eventually shallow subtidal shelf strata.The position of the MIS 5e reef facies is not well-constrained by our core data in the Sandy Point area.The Gulf-01 core just to the east, as well as coastal outcrops along the southern coast, contain a thin Acropora palmata facies attributed to a reef crest environment.The true shelf edge reef may have been farther outboard, but this cannot be constrained by our data.
With this arrangement of facies/depositional environments, the most effective way to represent lateral facies transitions is through a pair of stratigraphic cross-sections, one oriented eastwest in the northern part of the study area cutting perpendicular to most of the facies trends, and one section oriented NNE-SSW along the axis of the interior lagoon and then cutting directly across the dune ridge at the southernmost position (Figs 10 and 11).As shown by Purkis & Harris (2016), and guided by observation of present day facies arrangements, the authors anticipated a degree of inheritance of facies patterns driven by aeolian ridge topography.However, the degree to which older aeolian ridges were preserved through the subsequent transgression was anticipated to be low (note in Figs 10 and 11 the rare occurrence of MIS 9/11 topography preserved prior to MIS 5e deposition).
The oldest unit encountered within the study area is that of an algal-doloboundstone.This dolomitized unit closely fits those described from the Supko core at the north end of the island, including the presence of the Stylophora affinis coral species.A Pliocene age is supported by this correlation and by the Sr 87/86 data from  Facies display rare to common meniscus cementation and preferential dissolution of ooid cortices.In core, this lithofacies commonly has extensive burrowing and rhizoconcretions, and is weakly consolidated.Three-dimensional stratigraphic architecture 221 helpful for constraining age, but as yet no material that is sufficiently pristine has been located.Continuous core mapping also allowed the authors to recognize that a distinct high-energy oolitic unit with foreshore elevations of À7.5 to 12 m rested below the proposed MIS 9 and could be traced across the entire San Salvador Island platform (unpublished core mapping data).Although there is no radiometric dating control for the oolitic unit, its higher palaeo-shoreline elevation relative to MIS 9 and its distinctly different (oolitic versus skeletal) composition favours the interpretation of deposition associated with MIS 11.Importantly, the next unit above deposition associated with the MIS 11 high energy oolitic complex unit is distinctly different, being almost entirely composed of skeletal grains, mostly foraminifera, coralline algae and mollusc shell material.It is mapped in continuity with the previously described Owl's Hole locality and was postulated deposition to be attributed to be either MIS stage 7, 9 or 11 by previous workers (Hearty & Kindler, 1993;Carew & Mylroie, 1995).With our core control the authors were able to illustrate that foreshore deposits of this next oldest unconformity bound unit occur at 0 to +1 m relative to present day sea level, eliminating the possibility that it is an MIS 7 deposit.The MIS 7 marine facies across the Bahamas is known not to have reached present-day sea-level elevation (Lisiecki & Raymo, 2005;Sherman et al., 2014) thus by a process of elimination making the best age model for this second-oldest unit MIS 9.The youngest facies/units in our cross-section are known to be laterally equivalent to MIS 5e on the basis of mapping to areas where U-Th ages have been determined (Carew & Mylroie, 1995).The identification of this MIS 5e depositional unit is dominated by oolitic foreshore and dune facies, as is consistent with deposition associated with MIS 5e across the island (Hearty & Kindler, 1993;Carew & Mylroie, 1995).Aspects of each of the six unconformity bound units will be summarized in succession, from oldest to youngest, discussing, in order, the lithology, thickness, age depositional facies composition and arrangement, and nature of the upper bounding surface (typically a subaerial exposure surface) in greater detail in the following sections.

Pliocene reef complex depositional unit
The oldest (and deepest) unit penetrated in the Sandy Point area is composed of a single facies, a pervasively dolomitized algal boundstone (Figs 10 and 11).The algal boundstone is a tan, massive dolostone.Both finger corals (Stylophora affinis) and head corals are common (Fig. 9A).Geopetal sediments composed of peloids, coral fragments, red algae and lithoclasts occur in primary shelter cavities.In situ corals and coral fragments display extensive encrusting by coralline red algae.Goniolithon is common to abundant.Peloids, composite grains and lithoclasts are present and are heavily micritized.
The algal boundstone is pervasively dolomitized by microsucrosic fabrics with average crystal size of 10 to 20 μm.All aragonitic allochems have been leached or recrystallized to low magnesium calcite or dolomite.Interparticle pore space is commonly occluded by isopachous dolomite cements.Peloid and skeletal grain boundaries display extensive micritization.Dolomitization results in an overall lower total porosity and permeability for the unit, where very little original depositional porosity remains, a process referred to as overdolomitization by Lucia & Major (1994) where the dolomitization process continues beyond replacement and results in dolomite cementation The average depth of the top of this unit is 19.5 m in the Sandy Point wells, displaying ca 5.5 m of differential surface relief.The deepest cored well, which extends to 45.5 m subsea, does not identify the base of the unit.The dolomitized algal boundstone facies observed in the lower part of all cores that reached sufficient depth is consistent with a similar depositional facies and setting observed at the north end of San Salvador Island and attributed to a Pliocene age based on palaeomagnetic data and other biostratigraphic evidence (Dawans & Swart, 1988;McNeill et al., 1988).The Strontium 87/86 isotope geochemistry performed on CL-16 further suggests that this package is likely of Late Pliocene age (Fig. 6).Every core drilled in the Sandy Point encounters this reefal unit.(Fig. 12).
The Pliocene depositional unit is capped by a well-developed subaerial exposure surface (0-40 cm thick laminar caliche) where pedogenic overprint is common, with alveolar fabric, soil pisoids and rhizoconcretions found across all cores that penetrate this unit.The exposure surface is also commonly marked by breccia fragments of the host rock floating in a poor to wellindurated matrix.

Pleistocene depositional units
Resting above a well-developed caliche surface capping the interpreted Pliocene section is a succession of five unconformity-bound limestone sequences that range from high-energy reef and beach facies on the exterior to low-energy burrowed wackestones and packstones in the interior.The upper three unconformity-bound units are reasonably well-understood in terms of age, representing the MIS 11, MIS 9 (only locally developed) and MIS 5e (surface) interglacials as discussed previously.Below these units and resting on top of the caliche surface that caps the Pliocene, are two less well-defined units that are dominated by burrowed whole-mollusc packstones and wackestones analogous to restricted marine lagoonal deposits of the Holocene on San Salvador Island (Hagey & Mylroie, 1995).The authors postulate that these lower two units, which are separated only locally by a clear caliche crust, represent the two highest eustatic peaks of the Early Pleistocene, being MIS 25 and MIS 31.Greater detail will be required to confirm such an interpretation, but for now they will be referred to as Early Pleistocene (EP) 1 (lower) and 2 (upper).These five Pleistocene successions will be discussed in order from oldest to youngest.Early Pleistocene 1 depositional unit The Early Pleistocene 1 unit covers the entire cored area with relatively constant 3 to 4 m thick burrowed peloid packstones to floatstones (Figs 9B, 10 and 11) and a rare algal boundstone facies penetrated by a single core.The peloid grain-dominated packstone to floatstone is primarily composed of large Peneroplis foraminifera and molluscs floating in a fine-grained (peloidal) matrix.This facies varies from a packstone to floatstone (same composition) as the molluscs vary in size with depth.There is preferential dissolution of the large bivalve fragments.The uppermost portion of this unit is occasionally marked by the presence of faint planar laminations, rare coral fragments and red algae.Blocky calcite cementation and micritization of selected grains is common in this facies.Only minor facies variation is observed, with locally developed grain-dominated packstones representing areas of increased wave activity.Unmistakable facies vertical asymmetry is observed, and the capping caliche that separates this unit from the overlying EP2 displays a patchy distribution.

Early Pleistocene 2 depositional unit
The Early Pleistocene 2 unit is closely comparable in thickness and facies development to the EP1, with the separation of units suggested by the local occurrence of a patchy laminar caliche (for example, CL-22, 27 m; Figs 9C, 10 and 11).The skeletal mud dominated packstone to floatstone displays articulated and whole, disarticulated mollusc shells floating in a peloid-skeletal matrix.This unit is largely devoid of sedimentary structures other than burrowing (mostly Thalassinoides), displaying a low overall diversity facies assemblage.Core displays common burrowing and a touching vug network following dissolution of abundant skeletal fragments.Photomicrographs (Fig. 9C) reveal mouldic porosity following the dissolution of mollusc and gastropod shells.This unit is predominately calcitic mineralogy.
Subaerial caliche crusts locally separate the EP1 and EP2 units from one another, with a welldeveloped caliche marking the boundary from EP2 and the overlying depositional units.

Pleistocene Marine Isotope Stage 11 depositional unit
Within the MIS 11 depositional unit, it is possible to observe a full spectrum of facies (Figs 9D, 9E, 9F, 10 and 11).Facies observed within this unit include peloid grainstone, trough cross-stratified ooid grainstone and peloid graindominated packstone, and a large-ooidgrainstone that displays planar laminated facies with fenestral porosity that commonly caps the unit where a strongly developed capping caliche make it easy to recognize across all cores.
Massive burrowed peloid grainstone facies (Fig. 9D) are present in the platform interior, ranging from medium to coarse grain in size.This lithofacies is commonly weakly cemented and is massive to faintly laminated with abundant burrowing.Photomicrographs (Fig. 9D) reveal an assemblage that is heavily leached, with resulting mouldic porosity.Grain boundaries are commonly micritized.Lithofacies is poorly sorted.Foraminifera (Miliolids and abundant Peneroplis), skeletal fragments, composite grains, peloids and rare red and green algae make up this facies.To a lesser extent a peloid grain-dominated packstone to grainstone lithofacies is present within the MIS 11 unit.This lithofacies displays a well-sorted nature and a high percentage of cementation with common evidence of pedogenic overprinting.In core, this lithofacies is faintly laminated with individual laminations displaying solution enhancement.Clasts of laminated grainstone are entrained within.Etched core slabs and photomicrographs display common foraminifera, Halimeda and red algal grains.Grains show a lesser degree of dissolution and a higher percentage of blocky calcite cementation than the skeletal mud dominated packstone to floatstone lithofacies of the Early Pleistocene.
The planar laminated ooid fenestral grainstone comprises the interpreted foreshore depositional environment present within MIS 11.This lithofacies is characterized by a well-sorted, fine to medium-grained texture with abundant peloids, red algal fragments, large ooids and skeletal fragments.Blocky calcite occludes interparticle pore space with common micritization of grain boundaries.The main difference in the MIS 11 unit when compared to the overlying MIS 9 strata or underlying older Pleistocene units is the widespread development of low-angle seaward-dipping laminae that display common to abundant fenestral (keystone) pore systems and the earliest occurrence of well-developed ooids.Importantly, the CL-21 core shows for the first time the development of a dune ridge that is approximately 10 m in thickness (Fig. 10).This pattern is mimicked by the overlying MIS 9 (Figs 10 and 11).This variable thickness across the unit is likely inherited in part by the original island topography (either depositional or erosional) in combination with sediment depocentres, and wind/wave action acting to focus sedimentation and the development of high-relief dune.

Pleistocene Marine Isotope Stage 9 depositional unit
Strata attributed to MIS 9 are relatively easy to differentiate from the underlying MIS 11 and overlying MIS 5e units.The contacts of this unit with the overlying MIS 5e and underlying MIS 11 are marked by well-developed calichelaminar crust fabrics.In terms of grain composition this unit is also distinct, being devoid of ooids, and instead consisting of two facies, a very weakly cemented skeletal grainstone with abundant planar laminations (cement <1%), and a skeletal grain-dominated packstone (Fig. 9G).The most abundant facies encountered, the skeletal grainstone is primarily composed of a coralline algae, foraminifera and skeletal fragments of highly variable thickness that illustrates both aeolian, as at the type section in Owl's Hole and Watling's Quarry, and shallow marine facies, as observed in core and less commonly in exposures on the south coast of Sandy Point (Hearty & Kindler, 1993;Carew & Mylroie, 1995).The primary difference between the aeolian and shallow marine facies is the aeolian type section (also present in CL-21) displays abundant, multi-directional high angle dipping laminations while the shallow marine facies exhibit low angle, seaward dipping planar laminations and commonly contains fenestral porosity.
Skeletal grain-dominated packstone (Fig. 9G) as observed in core ranges from 0.25 to 1 m in thickness, and in some core sections appears to have lapped out or been removed by erosion (for example, CL-4 in Fig. 11).This lithofacies displays faint laminations and common burrowing.The facies is overall poorly sorted, medium to coarse-grained texture, with abundant coated grains, Peneroplis foraminifera, peloids, skeletal fragments and rare coral fragments are present.This lithofacies is distinct in its lack of cementation and high degree of grain leaching and micritization compared to all other lithofacies encountered.

Pleistocene Marine Isotope Stage 5e depositional package
This unit can be broken down into four lithofacies that are present in core and outcrop within the Sandy Point study area, including high angle stratified ooid grainstone, swash-laminated ooid grainstone, massive to trough-crossstratified ooid grainstone, and massive, burrowed ooid grain-dominated packstone (Fig. 9H).The Gulf-01 core captures reef lithofacies not penetrated by cores in the Sandy Point study area.The lithofacies of this unit often displays well-sorted peloids and ooids that are moderately cemented by meniscus and blocky calcite spar.
The high angle stratified lithofacies often display multi-directional dips, with differential cementation across the core intervals grading into low angle, seaward low angle planar dipping ooid grainstones.This unit displays differential thickness across the study area, influenced from antecedent topography from the underlying units, where the thickest section occurs proximal to modern day sea level.This unit is very commonly friable/rubble when recovered during coring but displays high porosity (often greater than 35%) and multi-darcy permeabilities.
The interior of this unit is dominated by a heavily burrowed occasionally faintly laminated, differentially cemented ooid grain-dominated packstone.This facies is typically between 1 m and 2.5 m in thickness, but represents a relatively flat top datum, filling in the underlying antecedent topography.
This lithofacies displays a fair degree of overprinting, with common rhizomorphs and laminated caliche stringers throughout the unit.A patchy, often discontinuous, caliche can be found on top of this unit, both in outcrop and in the subsurface.

DISTRIBUTION AND ABUNDANCE OF NON-MATRIX PORE SYSTEMS
The surface and subsurface geology and detailed sequence stratigraphy at Sandy Point shed light on the high proportion of non-matrix features and subaerial exposure surfaces encountered throughout the study area, highlighting their impact on diagenetic patterns, fluid flow, erosion and deposition through time.Greater detail will be provided about the significance and importance of these features in the following section.
Seventeen boreholes were used to characterize non-matrix pore networks in the Sandy Point area.Core and overlapping image logs were available in 13 of the 17 boreholes.The remaining four boreholes (CL-1, CL-3, CL-14 and CL-16) Three-dimensional stratigraphic architecture 227 collapsed, preventing acquisition of images, and thus were limited to detailed core description for characterization.Image logs are a critical dataset to account for open cavities (with size ≥ core diameter) that otherwise could be easily overlooked in core due to the absence of record or occurrence of a rubble zone (Fig. 12).On the other hand, core provides a scale of observation that allows for interpretation of non-matrix features that are currently filled with either cement or sediment.To quantify the non-matrix proportion at the wells, the cumulative thickness of each feature type was divided by the total interval of observation (i.e.length of logged boreholes).The reliability of estimates of non-matrix proportion may vary from well to well depending on length of observation interval (wells are of similar length, ca 30 m) and types of datasets used (core versus core and image logs).In the case of coreonly observations, the number of open features might be under-represented as a result of core recovery issues and the inability to properly capture the thickness of open features with the core (i.e. a large vug or cave in a core interval is simply welded together when the bit drops to the lower level).An attempt to quantify open intervals has been made using bit drops and missing core recovery to help with this problem in wells with no image logs (Breithaupt et al., 2022).
Of the non-matrix intersected by the boreholes, 75% correspond to open features, while the remaining 25% are features filled with speleothems, or internal cave sediments ranging in grain size from micritic to breccia fragments up to several decimetres.The proportion of nonmatrix in the boreholes sums to a total of 8.9%, where the remaining 91.1% is deemed matrix.The types of non-matrix are distributed as follows: 4.3% touching vugs, 3.5% caves and 1.1% voids.To evaluate vertical trends in non-matrix distribution, a frequency distribution plot was constructed that aggregates all non-matrix observations in the Sandy Point wells.The mid-point depth of each interpreted non-matrix feature was used, both open and filled, and mean sea level was used as common datum.Results show a bimodal distribution of non-matrix features with two distinct peaks (Fig. 13): a shallow peak at ca 8 m below mean sea level (bmsl) and a deeper peak at ca 16 m bmsl.While the distribution is dominated by open features, the amount of filled features is greater near the peak at 16 m bmsl (Fig. 13B).Both peaks in the frequency distribution show a strong correlation with the location of two major exposure surfaces interpreted in core: the top Pliocene and top Early Pleistocene (EP2).The box plots to the right of the histograms in Fig. 13 show the range of depths observed for the exposure surfaces at the top of the Pliocene, EP2 and MIS11.The exposure surface following EP1 was not plotted due to its patchy, discontinuous nature across the study area.
Interestingly, the majority of the non-matrix features appear to develop just above the mean depth of the exposure surfaces, with a rapid decay in frequency immediately below.In order to take a more in-depth look at the vertical distribution of open features around exposure surfaces, the proportion of open features within a 5 m window centred at the depth of the top Pliocene and top EP2 within each well was evaluated.Results show total 9.9% proportion of non-matrix at the EP2 window, where 5.8% of the non-matrix sits above the exposure surface, while 4.1% occurs below.The dominant types of non-matrix identified within this window are touching vugs with a total 5.6% (1.6% below and 4% above top EP2), followed by 2.6% of caves (1.4% and 1.2% above and below, respectively).Non-matrix features at the Pliocene window account for a total of 3.6%, with 3.5% developed above the top Pliocene.In this case, caves are the dominant feature, with 2.4%, followed by touching vugs, 0.9%.

DISCUSSION
The tightly spaced boreholes and subsequent collection of core allows for a detailed reconstruction of the subsurface stratigraphy within the study area.The combination of geochemistry, stacking patterns, superposition, and welldefined subaerial exposure surfaces allow the authors to build a robust, age constrained stratigraphic framework that is tied to the sea-level curve through the Late Pliocene.This work presents a detailed picture of how the high frequency, high amplitude sea-level fluctuations associated with icehouse conditions affect both the deposition and erosion of an isolated carbonate platform.The development of subaerial exposure surfaces during sea-level lowstand, that are often marked by well-defined and heavily cemented caliche surfaces drive clear differences in porosity and permeability above and below said surfaces.This relationship is clear when the presence of non-matrix features are evaluated in regard to their depth of occurrence.

Age relationships, stratigraphic implications, sea level and subsidence
Six unconformity-bound stratigraphic units comprise the sections in Sandy Point, with the oldest unit most likely late Pliocene and the youngest at ground surface being of Marine Isotope Stage 5e (125 ka).The oldest reefal dominated unit was encountered by every core drilled in the Sandy Point study area suggesting that it is representative of relatively stable, widespread subtidal reef growth, as with other coralcoralline algal facies seen on many Bahamian platforms (Vahrenkamp & Swart, 1981;Beach & Ginsburg, 1988;McNeill, 2005;Figs 10 and 11).Given the high surface relief of this unit and the brecciation across the cores drilled, this surface is interpreted to represent one that has experienced extensive dissolution/reprecipitation and erosion but shows no signs of intra-sequence exposure.The evidence for extensive subaerial exposure coupled with an age date of Late Pliocene from the presence of Stylophora affinis corals and strontium age dating (Sr 87 /Sr 86 values from 0.709020 to 0.709060), agrees with the interpretation from the Sr isotope data suggesting a 2.5 o 3.0 Ma age which is consistent with the ages derived from the Supko core (McNeill et al., 1988).This would imply that the island surface was subaerially exposed for 1.0 to 1.5 Myr, during which mean relative sea level was below current day sea level (Fig. 14) until the Middle Pleistocene .
Previous estimates of mean subsidence rates for Bahamian platforms have been summarized by McNeill (2005), showing San Salvador Island to have a mean subsidence rate of 26 m/Myr.Examining data from the CL-16 core (Fig. 8) an estimate of 6 m/Myr was derived, substantially lower than most modelled subsidence rates.This lower value comes from the general trend in elevation of the top Pliocene, which trends deeper northward, occurring at 15.5 m below present day in CL-16 and 24.7 m in the Supko core.A nearby core gathered for a related study, the Line Hole deep well, intersected the top Pliocene surface at À25 m, with the correlative top-Pliocene dolomite unit at À41 m.It is likely that in fact the 6 m/Myr value given for the CL-16 core is a maximum estimate because it appears that the top Pliocene in CL-16 is not the true top Pliocene but a somewhat eroded (i.e.deeper/older) portion of this unit.
A low-energy shallow subtidal shelf setting is postulated for both EP1 and EP2 on the basis of their mono-specific (often large, thin-walled molluscs), massive facies assemblages that often lack sedimentary structures other than burrowing and lack variability.This suggests a moderate-amplitude eustatic/climatic forcing.The absence of a transition to reefal or shelfedge grainstones within these units within the cored area suggests that either these higherenergy facies are not present (such as a Cay Sal Bank model) or they only exist close to the platform margin outside of the study area.Little evidence for energy asymmetry across this narrow arm of the island complex supports the deeper but well-oxygenated seafloor, analogous to platform interior facies of the Great Bahama Bank which average 5 m water depths (Harris et al., 2015).The low energy environments coupled with often patch caliche crusts are consistent with the shorter-term 40 kyr obliquity forcing and the lower amplitude of these events.
The depositional unit proposed herein as MIS 11 represents a significant shift in depositional style from any previous units, which were all relatively uniform in facies composition and were strictly subtidal whereas MIS 11 is dominated by high energy, oolitic facies.These facies display a Fig. 14.A pseudo wheeler diagram for the A-A 0 cross-section with the seven stratigraphic packages identified, from Pliocene to Holocene in age. Figure depicts deposition laterally and vertically for the various lithofacies.This diagram depicts periods of nondeposition as empty space, highlighting the significance of the periods of subaerial exposure this island has gone through during sea-level lowstands.Ages for each unit are derived from previous surface and subsurface studies of San Salvador Island and the Bahamas (Supko, 1970(Supko, , 1977;;McNeill et al., 1988;Hearty & Kindler, 1993;Carew & Mylroie, 1995;Hearty, 1998;Godefroid, 2011;Kerans et al., 2019b).
high diversity assemblage of fauna, large, welldeveloped ooids, planar faintly seaward dipping laminations with abundant fenestral porosity.This unit is assigned to the MIS 11 sequence by comparison with the Star Town unit from West Caicos that bears a U-Th age of 410 ka (Kerans et al., 2019b) and on the basis of its position directly underlying the Owl's Hole MIS 9 unit (Hearty & Kindler, 1993).Considering subsidence between MIS 11 and MIS 9, this relationship supports a peak interglacial sea level for MIS 11 higher than that for MIS 9.
An important observation for deposition attributed to MIS 9 is the presence of well-developed low-angle seaward-dipping foreshore strata with fenestral porosity just at 0 m (present day sea level).In the coastal exposures along the south margin of Sandy Point massive poorly sorted skeletal grainstone-rudstone with large fragments of coralline algae and molluscs are found, separated from the overlying MIS 5e grainstones by a welldeveloped caliche.This observation allows the authors to infer that this unit is of MIS 9 affinity rather than MIS 7 (the other skeletal-dominated unit of the mid-late Pleistocene).Previous studies (Hearty & Kindler, 1993;Carew & Mylroie, 1995) lacked the subsurface control generated here, and thus only noted aeolian facies for the Owl's Hole unit, making it difficult to differentiate between MIS 7 or MIS 9 interglacials.The 0 m foreshore elevation is significant because it allows confident identification of MIS 9, because MIS 7 deposits that are in some respects similar to MIS 9 did not reach above À10 m relative sea level.The sedimentary structures and fabrics are distinctly shallow marine in origin and not aeolian.Such an association precludes a placement as MIS 7 because the MIS 7 shoreline reached a maximum of approximately 18 m below present-day (Dutton et al., 2009).Other distinctive aspects of the MIS 9 unit are the development locally of aeolian dune ridges, such as in the Watling's Quarry, and in CL-21 core.For MIS 9 it is clear that the transgression was of significantly lesser magnitude than for MIS 11 or MIS 5e, as the unit pinches out in several places on positive topography.It is possible that transgressive reworking by MIS 5e may have removed some of the weakly lithified MIS 9 strata but, in contrast to MIS 5e, there is no low-energy interior lagoon/tidal flat, with thin grainstones blanketing the interior and foreshore deposits occurring on the south and west flanks of the promontory (Fig. 15).
The MIS 5 unit is exposed at surface across the Sandy Point area, and is assigned to the MIS 5e substage as it was mapped both in core samples and in outcrop and can be shown to be coeval and contiguous with strata of MIS 5e referred to as Grotto Beach Formation (Carew & Mylroie 1985, 1995;Hearty & Kindler, 1993, 1995).Welldeveloped dune ridges, >20 m in height were intersected in the CL-19, CL-21 and CL-22 cores, which transition downward into foreshore and in some cases shallow subtidal shoreface deposits or, in the case of the Gulf-01 core, coral reef facies highlighting a full suite of depositional environments represented.During MIS 5e, a welldeveloped dune ridge complex was deposited, creating a low energy, protected lagoon/intertidal environment within the interior of the study area (Fig. 15).The sea-level position for most of the observed foreshores in Sandy Point cores is 2 to 4 m above present-day, with a few occurrences of +6 m elevations (for example, CL-22).An excellent correspondence exists between the LiDARbased geomorphology and the core-based facies scheme, allowing a detailed reconstruction of the MIS 5e time-slice (Fig. 15).
The aforementioned relationships highlight three main styles of carbonate deposition and facies architecture observed in the Pliocene-Pleistocene history of the Sandy Point region: (i) Pliocene dolomitized reef complex; (ii) subtidal lowenergy lime wackestones, packstones and floatstones of the Early Pleistocene (units EP2 and EP1); and (iii) higher energy grain-dominated shoreface-foreshore-dune ridge systems of the Later Pleistocene (MIS 11, MIS 9 and MIS 5e; Fig. 16).All three phases are distinct both in terms of deposition and diagenetic overprint.The uniform subtidal reef cap of the Pliocene is consistent with this well-known warm period of the Neogene.Widespread dolomitization of this interval has been linked to seawater pumping by Kohout convection (Vahrenkamp et al., 1991).
The EP1 and EP2 units are very similar both in terms of their predominantly low-energy welloxygenated deeper subtidal platform-top facies pattern and the absence of well-developed laminar caliche exposure caps marking the boundary between these two units.Although the Early Pleistocene is typically linked by most stratigraphers to classical icehouse conditions, the isotopic and sea-level record of this subdivision of the Pleistocene has distinctly lower amplitude sea-level oscillations than those that characterize the Late Pleistocene (Fig. 2).The shift in amplitude and frequency of the climate signal from the 41 kyr obliquity pacing of the Late Pliocene-Early Pleistocene changes to the 100 kyr high amplitude Three-dimensional stratigraphic architecture 231 Fig. 15.Interpreted aerial extent of each lithofacies and their progression through time, from left to right.Each well location is identified on the map along with the modern day shelf edge and island extent.Lithofacies legend is colour coded and consistent for each time step.Overall, each borehole drilled in the study area penetrated a dolomitized reef unit at 20 m depth, suggesting that the southern part of the island was dominated with reefal communities.This unit is believed to be Pliocene in age.During the earliest Pleistocene (EP1), a small aeolian dune (catenary island) developed on the southern tip of the island, consequently controlling energy, where a broad low energy intertidal/lagoonal platform interior developed.EP2 represents a time when the study are experiences a period of restriction where all boreholes encountered the same low diversity assemblage floatstone with large mouldic mollusc fragments.Restriction here could be caused by a broad platform interior, an energy barrier (barrier reef?) or by the possible underlying erosional antecedent topography.During MIS 11 time the platform has returned to a high energy environment with the development of thin aeolian dune ridges and extensive foreshore deposits.During MIS 9, deposits onlap (drape) the MIS 11 topography, often filling lows with no deposition on the interior of the study area.Foreshore deposits are extensive on the platform during this time.During MIS 5, aeolian dune ridges rim the edges of the island, creating a broad lagoon on the interior with energy from the north-west.Thin foreshore deposits rim the edges of the dune edges, often forming strand plains.The Holocene time period is very reminiscent of MIS 5 time period, except that the south part of the island (dune ridge) has been eroded.Foreshore strand plan deposits are more widespread.Three-dimensional stratigraphic architecture 233 signal of the Late Pleistocene across a transition known as the Middle Pleistocene Transition or MPT (Lisiecki & Raymo, 2005;Miller et al., 2011).Much discussion of the nature and origin of this change in climate forcing exists, but the details of this change are beyond the scope of this presentation.However, the stratigraphic result of this transition is particularly relevant.Other workers (e.g.Yamamoto et al., 2006) have demonstrated a transition from non-coral reef bearing to coral reefbearing sequences in Japan tied to the MPT.In the bank margin record of the western margin of the Great Bahama Bank, detailed study of the Clino and Unda wells (Kievman, 1998;McNeill et al., 2001) also highlight a shift from Early Pleistocene parasequences with subaerial exposure only in the platform interior to Late Pleistocene cycles where exposure surfaces extend across the flat-topped aggradational platform, requiring an increase in high-frequency eustatic amplitude.Herein it is proposed that the distinct change in stratigraphic style from more uniform low relief flooded shelf cycles to those with complex preserved aeolian topography and high lateral facies heterogeneity is a reflection of this increase eustatic amplitude but also likely an increased climate gradient.The evolving climate across the Middle Pleistocene Transition was the likely driver for this shift.Further, it is likely that the increased sea-level amplitude and longer cycle time of the 100 kyr forcing drove greater meteoric dissolution during sea-level lowstands, and thus the increased development of caliche deposits and cave systems (James, 1972).Breithaupt et al. (2021) have presented compelling evidence that these laminar caliches are central to patterns and localization of meteoric dissolution; the authors add to the discussion below.

Stratigraphic controls on non-matrix distribution
The strong correlation between exposure surfaces and maximum occurrence of non-matrix features (fractures, vugs and caves) suggests a genetic control of exposure zones on karst distribution (Fig. 13).The broadly accepted model of freshwater lens diagenesis in island karst would predict non-matrix features to occur below major exposure surfaces (Mylroie & Carew, 1995).The present observations, however, suggest caves and other non-matrix features, preferentially developing atop of exposure zones.The top Pliocene is the best example of such vertical arrangement, with 2.9% of caves and 1.86% of touching vugs (out of a total of 8.9% non-matrix encountered in all units) intersected by boreholes (Fig. 13).The high frequency of nonmatrix features above the Pliocene unit is likely linked to the 1.0 to 1.5 Myr of exposure prior to the continuations of carbonate deposition, which would have developed diagenetically modified fabrics (caliches) of very low permeability (<1 md, averaging 0.5 m in thickness) that effectively behaved as a barrier or baffle to vertical fluid flow in overlying units (Fig. 8B).The time gap at the top of EP2 is likely on the order of 100 to 500 kyr, resulting in an exposure zone of similar permeability and thickness to that at the top Pliocene, and 1.8% of caves and 5.5% of touching vugs forming above it.In contrast, those units exposed for relatively short timeframes (MIS 11, MIS 9 and MIS 5e), have thinner exposure fabrics, and tend to show little to no relationship between the exposure surface and increased non-matrix feature development in close proximity.The lack of non-matrix features clustered around units with relatively short exposure durations before the next episode of deposition suggests that the development of non-matrix features occurs after the exposure has developed, and not penecontemporaneously with the deposition of the bedrock.
Evidence that caves and vugs form as a result of low permeability exposure fabrics limiting vertical fluid flow is consistent with the model proposed by Breithaupt et al. (2021), similar to those proposed by Beach (1995) and Melim (1996).Beach (1995) proposed a three-step process that ultimately effects porosity, permeability and fluid flow in the surface.Beach (1995) suggested that diagenetic stage I was marked by vadose diagenesis and subaerial exposure, and the development of a well-indurated surface.During stage II, shallow burial and inundation of ephemeral freshwater phreatic conditions occurs, where sediments alter to low magnesium calcite and secondary porosity develops.Step III is marked by deeper burial and prolonged episodes of residence freshwater phreatic and mixing zone waters that result in widespread dissolution.
The model proposed by Breithaupt et al. (2021) differs slightly from that of Beach (1995) where meteoric waters enter the subsurface through vadose fast flow routes and create ephemeral aquifers perched on low permeability horizons (an exposure zone for instance).Waters are forced laterally into the bedrock above the perching interface driving dissolution as undersaturated waters flow away from the locus of injection (Breithaupt et al., 2021).Between rainfall events, waters that Three-dimensional stratigraphic architecture 235 have equilibrated with respect to carbonate minerals drain through perforations in the exposure zone, removing saturated waters which primes the system for dissolution during each subsequent recharge event.Over numerous episodes of filling and draining cycles, voids enlarge, with some collapsing and others remaining open in the subsurface as isolated globular chambers.
Applying this model to the Sandy Point region, focusing dissolution by low permeability zones is a plausible explanation that explains preferential concentration of non-matrix features above exposure surfaces (Fig. 17).The length of the exposure dictates the diagenetic maturity and thickness of the exposure zones which subsequently promotes the formation of transient aquifers (filled with groundwater during highstand and drained dry during lowstands).Dissolution within the bedrock overlying these exposure zones can result in hydraulically efficient non-matrix networks.As a result, the permeability architecture is tied directly to the stratigraphy, with high and low permeability end members forming in pairs around major sequence boundaries.
Preservation of open non-matrix features may not always occur because they may collapse, become filled with precipitates, and/or the host rock containing them may be entirely stripped by denudation.As a result, the preservation potential of open features relies, in part, on deposition outpacing denudation to preserve the host rock and provide sufficient structural overburden to prevent roof collapse in the vadose zone (Nolting et al., 2021).Despite sufficient overburden constraints, large features such as caves may still collapse (Breithaupt et al., 2022;Nolting et al., 2021) and all non-matrix features have the potential to become filled by precipitates as porosity is redistributed thorough the system (Budd, 1988).The authors conclude that the stratigraphic and permeability architecture of a carbonate platform sets up the framework for non-matrix development, while the abundance of non-matrix features is tied to both their formation and preservation potential.This finding has strong implications for developing concepts for subsurface models of karst distribution and, potentially, connectivity.

CONCLUSIONS
The massive dataset collected in the Sandy Point area on the south-west tip of San Salvador Island serves as a high-resolution model for carbonate facies, stratigraphic architecture and patterns of meteoric diagenesis and karst development during the last 3 Myr of Earth history.Six unconformitybound sequences have been correlated and mapped throughout the 15 sq km area.Representative detailed cross-sections and palaeoenvironmental maps summarize the high-resolution framework created (Figs 10,11 and 14).The oldest sequence is of Late Pliocene age according to Sr 87/86 isotope data (Fig. 6) as well as correlation with palaeomagnetically constrained strata from the Supko well (McNeill et al., 1988).Coral-coralline algal reef facies dominated the Late Pliocene sequence, containing Stylophora affinis corals which are consistent with a Pliocene age range.An extensive karst and laminar caliche profile developed atop the Late Pliocene unit which is consistent with the longest-extant subaerial hiatus in the described stratigraphic succession.The next two sequences, Early Pleistocene 1 and 2, are composed of facies attributed to a distinctly lower energy environment than the underlying and overlying units.The inference of a MIS 25 to 37 age is proposed on the basis of the relatively high sea-level position for these interglacials (EP2 and EP1).The next unit is best tied to the MIS 11 interglacial and is composed of mixed skeletaloolitic sequence that displays a well-developed laminar caliche cap across the majority of the study area.The highest present-day relative sealevel indictor for MIS 11 is the foreshore to dune transition seen in CL-19.The overlying MIS 9 unit is mapped from the surface type section at Owl's Hole and Watling's Quarry into the subsurface and is a thin draping unit that locally pinches out on the underlying MIS 11 (for example, CL-4, Fig. 10).A full range of high energy deposits are observed including aeolian, foreshore and shoreface facies with a relative sealevel indicator in foreshore facies of 0 m.The MIS 5e sequence is the dominant unit in outcrop in the Sandy Point area, and contains aeolian dune, foreshore, shoreface and lagoonal/intratidal facies.The predominant grain types in MIS 5e are ooids, with lesser peloids and, in some settings, skeletal fragments and coral-coralline algal material.The MIS 5e sequence (Grotto Beach Formation) has been dated via U-Th measurements from coral material at Grotto Beach and the Gulf area, placing it clearly within the last interglacial (Donohue, 1991;Carew & Mylroie, 1995).
The secular trend from warmer, higher sea level Pliocene reefal deposits, low-energy, low-faciesdiversity Early Pleistocene sequences, and the complex high-energy grainstone systems of the Mid-Late Pleistocene can be matched to the climatic forcings recorded in the marine isotopic record for the last 3 Myr.Climatic forcing that drove the transition from the 40 kyr dominant obliquity signal of the Late Pliocene-Early Pleistocene to the 100 kyr high amplitude signal of the Mid-Late Pleistocene known as the Middle Pleistocene Transition was not related to orbital forcing but rather other climate drivers (Clark et al., 2006;Lisiecki et al., 2021).These distinct sequence compositions and architectures show that a lumping of Late Pliocene-Pleistocene successions into 'icehouse' (cf.Read et al., 1990) is not sufficient to capture the variability observed.
San Salvador Island is representative of many Cenozoic Bahamian carbonate platforms in displaying extensive karst and cave development.The Island of San Salvador served as the type example of the flank margin cave model (Mylroie & Carew, 1990) and extensive mapping and measuring of karst features has occurred.The spatial relationship of secondary pore networks and vuggy pore systems is essential for understanding fluid flow in both aquifers and subsurface reservoirs.Using the high-resolution image log dataset, combined with the newly created stratigraphic model, the spatial distribution of vuggy pore systems was mapped and linked to periods of prolonged subaerial exposure.More importantly, this detailed 3D framework confirms the work of Breithaupt et al. (2021), illustrating the importance of differential dissolution and caliche formation at exposure surfaces for controlling lateral focusing of touching vug pore networks (Figs 12,15 and 17).These non-matrix pore systems are concentrated in particular at the top Pliocene unconformity and the top Early Pleistocene (EP2), surfaces that reflect time gaps of 1.5 Myr and 0.5 Myr, respectively.This linkage does not imply that all non-matrix features developed during these specific periods of subaerial exposure, but rather that the cementation and development of low-flow caliche horizons mantling these surfaces generated aquitards within the stratigraphy that served to concentrate lateral flow during younger episodes of subaerial exposure, and resulted in concentrated aggressive groundwater dissolution and abundant vuggy porosity in comparison to the lower parts of the other four units.
Together, the refined petrographic, sequence stratigraphic and chronostratigraphic aspects of the Sandy Point model provide a new appreciation for the impact of changing climatic forcing on facies type and architecture of ancient carbonate platforms.Integration of non-matrix features within this framework highlights the importance of dissolution/cementation patterns on formation of these modified unconformity surfaces for focusing later meteoric diagenesis and creating enhanced fluid flow pathways for continued multicyclic diagenetic events and later reservoir behaviour.

Fig. 1 .
Fig. 1.Location map of the Bahamian Archipelago with a zoom in view of San Salvador Island showing its respective locations within the Archipelago.The study area, Sandy Point, is highlighted in the south-western tip of the island.

Fig. 2 .
Fig. 2. (A) Sea-level curve modified after Lisiecki & Raymo (2005) through the Early Pliocene (5 Ma).(B) A zoomed in portion of the sealevel curve for the last 400 ka, where periods of deposition are highlighted in blue.Marine Isotope Stages (MIS) are labelled for the obvious peak sea levels (highstands).Numerous lowstands have occurred within the Pleistocene and Pliocene Epochs, where carbonate platforms have the potential to be subaerially exposed and undergo extensive early diagenetic processes.

Fig. 4 .
Fig. 4. (A) Aerial image of the Sandy Point study area on San Salvador Island.Location of each borehole that was drilled within the area is designated along with cross-section traverses (see Figs 10 and 11).(B) Bare Earth LiDAR images for the same locale as the latter.Well locations are marked as well as the mapped surface geology (Kerans et al., 2019a).

Fig. 5 .
Fig. 5. Example of non-matrix features observed in borehole images.Interval of interest is highlighted by the grey bar on the far right.From left to right: measured depth (in metres), minimum (min CALI) and maximum (max CALI) acoustic calliper, optical borehole wall image, amplitude and travel time borehole image.A low amplitude and high travel time signature correspond to features open at the wellbore: (A) 30 cm tall open cave, note stalactites captured by optical and acoustic image logs; (B) example of centimetre-scale open voids; (C) dense network of touching vugs; circulation was lost while drilling this interval, resulting in all drilling fluids lost to the formation; (D) breccia (likely related to cave collapse); (E) cavity mostly filled by a speleothem, note crystal growths observed in amplitude image log; (F) caliche associated to exposure surface.

Fig. 7 .
Fig. 7. Example of data suite used for this study.Example details the CL-5 borehole.Data used within this study included detailed core descriptions (inch scale), thin section descriptions on every foot scale, calliper log, optical, acoustic and travel time image logs (when available).The aforementioned data were used to identify facies types and variations, subaerial exposure surfaces, relative age relationships, and secondary porosity (non-matrix) features present in the subsurface.

Fig. 10 .
Fig. 10.A-A 0 cross-section, trending west/north-west to east/south-east (see Fig. 4 for exact location) cross-cutting the study area.The topographic expression is controlled by bare earth LiDAR image that has been merged with a multi-beam.Core descriptions for the individual boreholes are overlaid on the exact location of each well.Blue dots signify depths of thin sections for each borehole, red asterisks correspond to thin sections used in Fig. 9. Dashed stratigraphic horizons mark regions where there is a lack of well control, thus representing regions where facies correlations are solely interpretive based on predictable facies associations from present day.Depositional environments are coloured according to the legend provided on the figure.Inset box details the zoomed in view of Fig. 12.A detailed breakdown of stacking patterns, correlations, thickness, exposure surfaces and ages is provided in the text.

Fig. 11 .
Fig. 11.B-B 0 cross-section, trending south/south-east/north to north/north-west (see Fig. 4 for exact location) cross-cutting the study area.The topographic expression is controlled by bare earth LiDAR image that has been merged with a multi-beam.Core descriptions for the individual boreholes are overlaid on the exact location of each well.Blue dots signify depths of thin sections for each borehole, red asterisks correspond to thin sections used in Fig. 9. Depositional environments are coloured according to the legend provided on the figure.A detailed breakdown of stacking patterns, correlations, thickness, exposure surfaces and ages is provided in the text.

Fig. 12 .
Fig. 12. Zoomed in view of cross section A-A 0 (Fig.10) with lithofacies greyed out.Wrapped optical image log overlaid on the borehole location with interpreted non-matrix feature on the left hand side.CL1 does not have image logs.Figure highlights non-matrix features, including brecciation, touching vug porosity, speleothems, fractures, caverns, and exposure surfaces and their byproducts were flagged and put into spatial context relative to the stratigraphic framework (see inset, Fig.10).

Fig. 13 .
Fig. 13.Depth frequency distribution of non-matrix features in the Sandy Point area.On the left panel, depth distribution of all non-matrix features as interpreted from borehole image logs and core.On the right, non-matrix features are colour coded according to whether they are open or filled.Depths are referred to mean sea level.The box plot to the right of the first panel shows the range of depths of the top Pliocene, Early Pleistocene 2 (EP2) and Marine Isotope Stage (MIS) 11 exposure surface picks from core.Note how the maximum frequency of non-matrix features occurs at and above the location of the exposure surfaces, and it quickly drops below the exposure surfaces.m bsl: metres below sea level.

Fig. 16 .
Fig. 16.Comparison between stratigraphic and facies architecture versus climatic/insolation forcing functions.The Late Pliocene coral-coralline algal reef complex that forms the lowest unit cored at Sandy Point formed during the global warm temperatures and muted 41 ka obliquity forcing of that time interval, showing little evidence of environmental/facies variability internally.The two low-energy Early (?) Pleistocene units mapped across Sandy Point are devoid of high-energy shoreline facies or distinct exposure surfaces, consisting mostly of subtidal molluscan wackestone-floatstone.In contrast to the Pliocene and Early Pleistocene units, the Mid-Late Pleistocene MIS 11-9-5e units are characterized by high-relief dune complexes and high-energy shoreline facies dominated by ooid grainstones, with coral reef facies observed in both MIS 11 and 5e at least.Each unit is bound by welldeveloped exposure surfaces and laminar caliches, with the best developed being those at the top of MIS 11 and MIS 5e.

Fig. 17 .
Fig. 17.Conceptual model of an isolated carbonate platform through time.Lowstand systems tract (LST) and highstand systems tract (HST).(1) Following deposition, sea-level falls, subaerially exposing the top of the carbonate platform, with a freshwater lens developing.The freshwater lens drives the development of secondary porosity (non-matrix) features, while there is increased cementation due to the development of a caliche (brown) on the top of the platform.(2) During sea-level highstand deposition occurs and a freshwater lens develops on top of the previous stratigraphic package, above the subaerial exposure surface (zone of increased cementation and decreased porosity and permeability.(3) Following sea-level lowstand, the subaerial exposure surface, permeability barrier, perches the freshwater lens on top of the oldest stratigraphic package while subsequently exposing the platform, forming a second, younger exposure surface (permeability barrier).

Table 1 .
Sandy Point Well locations with total depth and core recovered.
Ó 2023 ExxonMobil.Sedimentology published by John Wiley & Sons Ltd on behalf of International Association of Sedimentologists, Sedimentology, 71, 207-240Three-dimensional stratigraphic architecture 217 our deep control well at Sandy Point, CL-16 (Fig.6).A minimum of two, and possibly more, unconformity-bound units that are dominated by shallow-shelf low-energy burrowed molluscan packstone occur above the algaldoloboundstones and below the proposed MIS 11.Additional constraints from older coral material within MIS 11 would be particularly

Table 3 .
Summary table of facies encountered in Sandy Point.