A new MIS 5 to MIS 2 speleothem record from Sandkraal Cave on the South African Cape south coast

Abstract

We present new stable oxygen and carbon isotope composite records (δ¹⁸O, δ¹³C) of speleothems from Sandkraal Cave 1 (SK1) on the South African south coast for the time interval between 104 and 18 ka (with a hiatus between 48 and 41 ka). Statistical comparisons using kernel-based correlation analyses and semblance analyses based on continuous wavelet transforms inform the relationships of the new speleothem records to other proxies and their changes through time. Between 105 and ~70 ka, changes of speleothem δ¹⁸O values at SK1 are likely related to rainfall seasonality. Variations of δ¹³C values are associated with changes of vegetation density, prior carbonate precipitation (PCP), CO₂ degassing in the cave, and possibly variations of the abundance of C₃ and C₄ grasses in the vegetation. The relationships of δ¹⁸O with other proxies shift between ~70 and 48 ka (Marine Isotope Stages 4–3) so that both stable isotope records now reflect CO₂ degassing, evaporation, and PCP. Similar relationships also continue after the hiatus for the deposition phase between 42 and 18 ka. Our findings support modeling results suggesting drier conditions in the study area when the Southern Hemisphere westerlies are shifted north and the paleo–Agulhas Plain is exposed.

Introduction

The present-day Cape south coast is influenced by dynamic and complex climate systems including westerlies, subtropical anticyclones, tropical low-pressure systems, and interactions among them (Tyson, 1986). These climatic systems in combination with variable landscape topography lead to steep present-day gradients in rainfall amounts and seasonality across the region over time (Fig. 1). Glacial–interglacial shifts of the climatic systems are likely to have caused considerable changes. Additionally, the shelf off the Cape south coast is up to ~250 km wide (Dingle et al., 1983). This magnified the effects of Cenozoic sea-level fluctuations, repeatedly exposing and submerging vast tracts of land. The large area that would have been exposed during low sea levels is now referred to as the paleo-Agulhas Plain (PAP; Marean et al., 2014), and its exposure would have increased the inland continental character of present-day coastal sites (van Andel 1989; Fisher et al., 2010; Göktürk et al., 2023).

Figure 1. (a) Map of southern Africa showing annual rainfall amounts (blue shading) and rainfall seasonality (dense dots: more than 60% winter rain; looser dots: 40–60% winter rain). The location of areas west of this receive more winter rainfall and regions to the east are dominated by summer rain. Sandkraal Cave (SK1) is indicated by an orange star. Also included are locations of paleoclimate proxies from the region that were included in correlation analyses and/or mentioned in the text: diamonds, pollen records included in CREST reconstructions (all diamonds were used for Tann reconstruction; brown, records used for PWetQ north; blue, records used for PWetQ center/east; orange, not included in precipitation reconstruction; Chevalier and Chase, 2015); triangles, marine sediment cores (light blue, rainfall proxies [MD08-3167: Collins et al., 2014; CD154-10-06P: Simon et al., 2015; MD96-2094: Stuut et al., 2002; CD154-17-17K: Ziegler et al., 2013]; brown, proxies of temperature and wind [GeoB-1710-3, GeoB1711-4, GeoB1712-4: Kirst et al., 1999]; orange, temperature reconstructions [MD96-2048: Caley et al., 2011; MD02-2594 and MD96-2080: Marino et al., 2013, Dyez et al., 2014].

Over the last decade, our team published speleothem records from five caves in southern South Africa (coastal sites: Crevice Cave, PP29, and Staircase Cave at Pinnacle Point near Mossel Bay, Herolds Bay Cave south of George; interior sites: Efflux Cave in the Little Karoo; Bar-Matthews et al., 2010; Braun et al., 2019, 2020). The records cover the latter part of the middle Pleistocene and the late Pleistocene. Before this, there was just one speleothem record in southern South Africa from Cango Cave in the Little Karoo that spanned the last 50,000 years with a break between ~14–6 ka (Talma and Vogel, 1992). This record was recently amended by the publication of stable isotope records of a second sample from Cango Cave (Chase et al., 2021). The interpretations by our group have mostly attributed changes of δ¹⁸O in speleothems to changes in rainfall seasonality (with higher values reflecting more summer rainfall; Bar-Matthews et al., 2010; Braun et al., 2019, 2020) based on modern records of rainfall isotope seasonality (Braun et al., 2017). The records from Cango Cave were driven by changes of drip water δ¹⁸O and temperature (Talma and Vogel 1992) or thought to be influenced by too many factors to parse (Chase et al., 2021). Speleothem δ¹³C records from the region have been generally interpreted as reflecting changes in the proportions of C3 and C4 plants in the vegetation (Talma and Vogel, 1992; Bar-Matthews et al., 2010; Braun et al., 2019, 2020; Chase et al., 2021). A short Marine Isotope Stage (MIS) 3 record from Maccali et al. (2023) was the first from this region to attribute changes of speleothem δ¹⁸O and δ¹³C to prior carbonate precipitation (PCP) and by extension changes in precipitation amounts based on a comparison with trace element records.

Although considerable advances have been made in reconstructing the response of the complex climatic systems in the region to large-scale glacial–interglacial climate variations to date, much remains to be learned about past climate dynamics on the Cape south coast. Understanding these past climate dynamics in the region is especially important in relation to the long archaeological record found here (Marean et al., 2014; Wadley, 2015) as well as the hyper-diverse extratropical flora of the Greater Cape Floristic Region (GCFR; Braun et al., 2023). Here we present new speleothem stable isotope records of seven samples taken from Sandkraal Cave 1 (SK1) on the Cape south coast near the town of George (Fig. 1). Together, these samples cover the time interval between ~104 and 18 ka (with a break between ~48 and 41 ka) and, for the first time on the south coast, include a record of the Last Glacial Maximum (LGM). We used Gaussian kernel-based cross-correlation analyses and continuous wavelet transform–based semblance analyses to investigate the statistical relationship between these new records and existing speleothem and other proxy records.

Regional setting

Geological background

The geology of the Cape in general is defined by the ~3000-m-high Paleozoic Cape Fold Belt, flanked by a dissected and relatively narrow coastal plain adjacent to the present-day coast. Sandkraal Cave 1 (SK1; 34.021°S, 22.529°E, 9 m above sea level [m asl]) is a sea cave southeast of the town of George on the Cape south coast and is also known as Ballots Baai Cave (Hitchcock, 2007). The cave is located in the Kaaimans Inlier, a geologic window that exposes formations predating the Paleozoic rocks of the Table Mountain Group, which dominate much of the Cape Fold Belt. The host rocks of SK1 mainly belong to the Victoria Bay Formation of the Kaaimans Group and consist of feldspathic quartzites and calcareous layers with boudinaged lenses of calc-silicate rocks (Krynauw and Gresse, 1980). The Rooiklip granites, which intruded the Victoria Bay Formation about 527.2 Ma (Chemale et al., 2011), now make up much of the steep coastal cliffs in the Kaaimans Inlier (Fig. 2c). The most likely source of carbonate for the formation of speleothems are calcareous quartzites in the Victoria Bay Formation and calc-silicate skarns at the contact between the Victoria Bay Formation and the Rooiklip granites, which contain up to 80% calcite (Krynauw and Gresse, 1980; Prent, 2013).

Figure 2. (a) Image taken outside Sandkraal Cave (SK1) showing the metasediment and mafic dyke host rock, speleothem formations that are now outside the cave, and remnants of cemented colluvium under the false flowstone floor. (b) Annotated sketch of the location shown in a. The current entrance to the cave under the false flowstone floor is shown in gray; a blue line indicates the sampling location of the flowstone section of samples 428464, 428465, 428466, and 428467. (c) Image of the slopes surrounding the cave entrance with annotated geology.

As with many caves and rock shelters along the Cape south coast, SK1 was cut into the coastal cliffs along a fault in the host rock during a sea-level high stand (Karkanas et al., 2021). The cave floor near the entrance is ~9 m asl, but the initial wave-cut notch at the back of the cave has its floor at ~13 m asl and the cave roof is level at about 16 m asl across the whole length of the cave (Supplementary Fig. 1). This cave morphology suggests that the main high stand that cut the cave was between 13 and 16 m above the present sea level, which corresponds to the MIS 11 sea level in this region (Roberts et al., 2012) and is at a similar elevation as several caves at Pinnacle Point (Fisher et al., 2010; Karkanas et al., 2021).

After the sea retreated, the cave was sealed by dunes and colluvium from the ~100-m-high surrounding cliffs. A thick flowstone layer (up to 1 m) formed on top of the debris cone that closed the opening, and substantive speleothem formations (some a few meters in diameter) developed on top of this flowstone (Supplementary Panoramas 1 and 2). After these features were laid down, the sediment cone closing the cave was eroded by a renewed sea-level high stand. The flowstone and speleothem formations, however, remained partly intact, still blocking much of the initial opening in the host rock (Fig. 2, Supplementary Panorama 1). The outlines of colluvium pebbles are still visible in the underside of the flowstone and small remnants of the dunes and colluvium are also preserved cemented to the sidewalls of the cave entrance below the flowstone (Fig. 2, Supplementary Fig. 2). The cave can now be entered by crawling through an ~0.7-m-tall and 7-m-wide opening under the flowstone layer. Behind this low entrance, the cave opens into an ~7-m-tall and 7-m-wide chamber (Supplementary Panorama 2). The cave is 52 m long and narrows gradually from its large size near the entrance to a small crawl space at the back (Hitchcock, 2007; Supplementary Fig. 1). Speleothem formations are common throughout the cave, in some places still actively forming (Supplementary Panorama 2). Tufa formation is also common on top of the flowstone surface outside the cave (Supplementary Panorama 1).

Climate and influences on the δ¹⁸O of rainfall

The central Cape south coast in the vicinity of the town of George is characterized by a mild oceanic climate with year-round rainfall that distinguishes the region from much of the rest of South Africa, where rainfall is highly seasonal (Fig. 1). A northern position of the austral westerlies in winter (June–August) brings the southwestern and western coasts of South Africa under the influence of rain-bearing cold fronts, while much of the east and interior of the country experiences dry conditions. When the westerlies shift south during summer (November–February), the west and southwest Cape regions receive less rainfall and a heat low develops over the interior that draws moisture, mostly from the Indian Ocean, toward the east and interior. The central south coast is influenced by complex interactions between the pressure systems that lead to summer and winter rainfall across the country. The most important rain-causing conditions are high-pressure ridges extending from the South Atlantic Anticyclone toward the ocean south of South Africa (46% of annual rain; usually after the passage of a westerly wave) and tropical temperate troughs that form when the continental summer heat low links up with a westerly low-pressure wave to the south (28% of annual rain; Engelbrecht and Landman, 2016). Ridges of the South Atlantic Anticyclone can occur year-round and are associated with stratiform rainfall (Engelbrecht and Landman, 2016). Overall, stratiform rainfall dominates during the winter months along the Cape south coast from the combined influence of ridging anticyclones with other synoptic types such as low-pressure troughs to the southeast and southwest of the region (Engelbrecht et al., 2015; Engelbrecht and Landman, 2016). Although ridging anticyclones also commonly occur in the summertime, tropical temperate troughs add a considerable contribution of convective rainfall during this season (Engelbrecht and Landman, 2016). This increased contribution of convective rainfall leads to amount-weighted average rainfall δ¹⁸O values that are 5‰ higher for summer rainfall than for winter rainfall in Mossel Bay (sampling of individual rain events in 2006–2007 and 2009–2012; Braun et al., 2017). Modeling of LGM climate for southern Africa suggests that northward-shifted westerlies lead to an expansion of winter rainfall in western and interior southern Africa, whereas rain-shadow effects along the coastal mountains cause a decrease in winter rainfall on the present-day Cape south coast (Engelbrecht et al., 2019).

The average annual temperature at SK1 is 16.2°C and mean annual rainfall is 890 mm (Schulze and Lynch, 2007; Schulze and Maharaj, 2007). There is only moderate rainfall seasonality, with each of the 12 months of the year contributing more than 5% to the annual total and only minor peaks in the transitional seasons (April and October).

Environment and vegetation

Southwestern South Africa hosts the diversity center of the Greater Cape Floristic Region (GCFR), one of the most diverse extratropical regions in the world (Bergh et al., 2014; Grobler and Cowling, 2022). The highest diversity within the region is observed in the western winter rainfall–dominated parts of the GCFR, but its biomes extend along the south coast, where rain falls year-round (Procheş et al., 2005; Rebelo et al., 2006). Being a Mediterranean-type shrubland, most of the GCFR plant communities are dominated by C3 plants, but C4 grasses can be common in some Fynbos and Renosterveld communities, especially in the eastern GCFR, where summer rain becomes more prevalent (Vogel et al., 1978; Cowling, 1983; Cramer et al., 2014). Succulent plants using crassulacean acid metabolism (CAM) occur mainly in succulent karoo, strandveld, and subtropical thicket communities (Cowling and Campbell, 1983; Cowling et al., 1994; Rundel et al., 1999).

SK1 is located near the east–west center of the GCFR. The vegetation on top of the cliffs has been heavily altered by agriculture and pasture farming. Without the human activity, Garden Route Granite Fynbos, which is a Sand Fynbos type with abundant C4 grasses in the understory, would be expected to grow here (Rebelo et al., 2006; Williams, 2015). The shrubby Grootbrak Dune Strandveld, which is poor in grasses but has a considerable component of succulent CAM plants, grows along the steep ~100-m-high cliffs at SK1 (Rebelo et al., 2006; Williams, 2015). Deep ravines in this area provide habitats that are more protected from fire than the surface on top of the cliffs and harbor the westernmost patches of the C3-dominated southern Afrotemperate forests on the Cape south coast (Rebelo et al., 2006; Williams, 2015).

The landscape outside the cave would have been considerably different for much of the Quaternary, as lower sea levels would have exposed the PAP on the now-submerged continental shelf to variable degrees. At some of the lowest sea levels, this plain would have been about 90 km wide near SK1 (Fisher et al., 2010). The current coastal lowlands, including the area near SK1, are characterized by rolling hills dissected by incised river valleys, whereas the PAP would have been largely flat with wide meandering rivers (Cawthra et al., 2020). Sand Fynbos most likely remained the dominant vegetation type on top of the cliffs at SK1 during the LGM, and the PAP just offshore would have hosted shale grasslands, floodplain woodlands and fynbos-thicket mosaics (strandveld; Cowling et al., 2020). The presence of grasslands on the continental shelf has long been suggested based on the high abundance of large-bodied grazers in Pleistocene fossil assemblages from sites along the Cape south coast (Klein, 1983; Marean et al., 2014). Fossil trackways of giraffe also suggest a flora more typical of savannah to the east and north (Helm et al., 2018). Recent modeling results suggest that decreased atmospheric CO₂ concentrations may have increased the abundance of C4 grasses in vegetation types with considerable grass components in the understory on the central and eastern south coast and PAP during the LGM (Grobler et al., 2023). It is currently assumed that the vegetation on the exposed PAP would have been similar during successive glacial phases, but this remains to be tested with further simulations and studies of paleo-archives. In summary, during lowered sea levels, SK1, as well as most now-coastal locations, sat at a stark topographic and vegetation ecotone.

Materials and methods

Speleothem material and sampling

A total of seven speleothem samples from SK1 (34.021°S, 22.529°E, ~9 m asl) were analyzed in this study. One stalactite, 162807, and one composite of a stalactite and a drapery, 357387, were found broken off outside the cave. Cross sections through flowstones were sampled in two places, one now outside (428464, 428465, 428466, and 428467) and one inside the cave (460362-1 and 460362-2). The flowstones were cut using a reciprocating saw and a chain saw equipped with diamond-tipped blades. More detailed descriptions and images of the samples can be found in Supplementary Appendix 1 and Supplementary Figure 3, respectively.

Because sample 357387 is a composite, three isotope transects were analyzed on it, one on a section through the drapery and two on opposite sides of the stalactite (Supplementary Fig. 3b and c). Samples 428464, 428465, 428466, and 428467 represent a continuous cross section through a flowstone, and thus all ages analyzed on these samples were used to construct a single age model. The isotope records of samples 428466 and 428467 were excluded from the composite, as their large variability and high detrital content suggests that they may have formed when the cave was still open. Tufa overgrowth was also present in samples 162807, 357387-B2, and 460362-1; stable isotope values of these sections are included in Supplementary Table 1 but were excluded from the normalization and the composite record.

Uranium-series dating

Aliquots for uranium-series disequilibrium dating were drilled across the samples following individual growth laminae using a hand drill equipped with 0.8- to 4-mm-diameter drill bits. Analyses were carried out at the Geological Survey of Israel and at the Isotope Geochemistry Laboratory at the University of Minnesota. Measurements at the Geological Survey of Israel (GSI) followed procedures described in Vaks et al. (2006, 2010) and Grant et al. (2012). Powdered aliquots weighing 0.2–0.3 g were dissolved in 7 N HNO₃. Detrital material was separated in a centrifuge (10 min at 4000 rpm) and treated with HNO₃ and HF. All samples were spiked with a ²²⁹Th + ²³⁶U tracer and recombined with the dissolved and evaporated detrital fractions. Uranium and thorium were separated by conventional anion exchange column chemistry using Bio-Rad AG 1-X8 resin with a 200–400 dry mesh. After the matrix was stripped using 7 N HNO₃, thorium and uranium were eluted using 6 N HCl and 1 N HBr, respectively. Both solutions were evaporated to dryness and dissolved in 5 ml (thorium) and 2 ml (uranium) of 0.1 N HNO₃. Measurements at the GSI were carried out using a Nu Instruments (UK) Multicollector Inductively Coupled Plasma Mass Spectrometer (MC-ICP-MS) equipped with 12 Faraday cups and 3 ion counters. The uranium and thorium solutions were injected into the MC-ICP-MS through an Aridus micro-concentric desolvating nebulizer sample introduction system. Instrumental mass bias was corrected for (using an exponential equation) by measuring the ²³⁵U/²³⁸U ratio using the natural value of 137.88. Ion counters were calibrated relative to Faraday cups for each individual analysis using several cycles of measurement with different collector configurations. The accurate determination of ²³⁴U and ²³⁰Th concentrations was possible by isotope dilution analysis using the ²²⁹Th + ²³⁶U spike. Sample 428464 was also dated at the University of Minnesota. Extraction and purification methods were as described by Edwards et al. (1987) and carried out on a Thermo-Finnigan Neptune MC-ICP-MS following the procedures of Shen et al. (2012) and Cheng et al. (2013).

Previously published decay constants for ²³⁴U and ²³⁰Th (Cheng et al., 2013) allowed the calculation of the ages. All ages were corrected for detrital thorium using the value of material in secular equilibrium with the bulk earth (²³²Th/²³⁸U = 3.8), assuming the initial ²³⁰Th/²³²Th atomic ratio of 4.4 ± 2.2 × 10⁻⁶.

Stable isotope analyses of oxygen and carbon

Samples 162807 and 357387 (A, B1, and B2) were analyzed at the GSI following Bar-Matthews et al. (1997, 2003). Aliquots for stable isotope analyses were drilled at ~0.5 mm intervals along the growth axis of the specimen using a handheld rotary tool equipped with a 0.8-mm-diameter drill bit, trying to achieve a small overlap between the holes to ensure continuous coverage of the stable isotope records. About 600–700 μg of the aliquots were weighed and put into glass vials at room temperature. Vials were dried uncapped in an oven at 60°C overnight. The vials were balanced horizontally, with the carbonate powder at the bottom and a drop of dry phosphoric acid (100%) at the top; capped; and then flushed with pure helium gas for 10 min. After the needles used for flushing were removed, the vials were turned to a vertical position to allow for the reaction of the carbonate with the phosphoric acid. Samples were left at room temperature overnight before analyses. The isotopic composition of CO₂ generated by the reaction was measured in a Finnigan Gas Bench II extraction system attached to a ThermoFinnigan DELTA plus XD mass spectrometer with a CTC Analytics automatic sampler at the GSI. A standard Carrara marble (δ¹³C = 4.74‰; δ¹⁸O = −3.85‰) was measured every eight samples. All δ¹⁸O and δ¹³C values obtained were calibrated against the internal lab standard, which was calibrated against the international standards NBS-19 (δ¹⁸O = −2.20‰; δ¹³C = +1.95‰) and NBS-18 (δ¹⁸O = −23.01‰; δ¹³C = −5.01‰) and are reported in per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB). Analytical reproducibility is better than 0.1‰ for both δ¹⁸O and δ¹³C.

Stable isotopic compositions of samples 428464, 428465, 428466, 428467, 460362-1, and 460362-2 were analyzed at the Metals, Environmental and Terrestrial Analytical Laboratory (METAL) at the School of Earth and Space Exploration, Arizona State University. Samples were drilled using a computer numerical control–operated Taig tools micromill equipped with an 0.5-mm-diameter drill bit. Samples were milled in a continuous trench at increments of 0.2–0.3 mm. Between 0.45 and 0.55 mg of sample powder were weighed into 12 ml borosilicate vials and dried uncapped but loosely covered with aluminum foil at 70°C overnight. After drying, the vials were capped and flushed with pure helium at 300 ml/min for 60 s. About 0.1 ml of pure phosphoric acid was added through the septum of the cap using a syringe with a 20 gauge needle. Samples were left to react at room temperature overnight before analysis. CO₂ was extracted from the vials using a ThermoScientific GasBench II and analyzed in a Thermo MAT 253 magnetic sector mass spectrometer. International standards NBS-18 and Carrara marble (δ¹⁸O = –2.01‰; δ¹³C = 2.1‰) were used for sample standardization to VPDB; both standards were analyzed every 12 samples. Beijing marble (δ¹⁸O = –10.9‰; δ¹³C = 1.6‰) was analyzed in between standards every 12 samples to test the accuracy and stability of the run. Analytical reproducibility was better than 0.15‰ for δ¹⁸O and δ¹³C.

Age modeling and composite record construction

Age models were constructed based on Bayesian statistics using the software Bacon (Blaauw and Christen, 2011) in its implementation in the R (R Core Team, 2022) package rbacon. This method uses millions of Markov chain Monte Carlo simulations to estimate the accumulation rate between the dated horizons, but a user-defined value (acc.mean) is also considered for extrapolation of ages along the edges of the sample and in hiatuses. We defined acc.mean based on the average deposition rate calculated for individual samples between measured ages and rounded to the nearest multiple of 100 (excluding sections with hiatuses). The accumulation rate used for sample 162807 was 1700 yr/mm. The age model for 357387-A has a hiatus of almost 50 ka length and had to be run in two sections. The accumulation rate was 200 yr/mm for both sections. Accumulation rates for the two sides of stalagmite 357387-B were 800 yr/mm (B1) and 600 yr/mm (B2). The samples 428464, 428465, 428466, and 428467 were combined into one age model with an average accumulation rate of 400 yr/mm, because they reflect a continuous cross section of a flowstone. The stable isotope records of samples 428466 and 428467 were not included in the composite record, because they probably formed before the cave was completely closed. Accumulation rates for samples 460362-1 and 460362-2 were 200 yr/mm and 100 yr/mm, respectively.

Bacon was chosen as the method of age modeling because it allows for the calculation of age models with error ranges even for small samples with limited numbers of dating points. For the larger samples included in this study that had more than three dating analyses run on individual sections between hiatuses (428464–428467, 460362-1, and 460362-2) StalAge age models were also constructed and show very good overlap with the Bacon age models (Supplementary Fig. 4).

The composite record was constructed by converting the stable isotope records of individual samples from SK1 into z-scores and aligning them by age using the mean age of the Bacon age model. The combined time series was then smoothed by a 5-point running mean to reduce the local noise. This approach generally follows the method of Williams et al. (2005), with the only change being that z-scores were used instead of aligning stable isotope records to match the average of the longest record among them. A gap in the data between ~48 and 41 ka is considered a hiatus, and we thus included a break in the composite record here. We also excluded the isotope records of samples 428466 and 428467 from the composite record, as these two samples were rich in sand grains and the δ¹³C values were slightly increased and more variable than in the other samples, which may indicate that the carbonate formed when the cave was still open. We estimated confidence intervals around the 5-point running mean of the stable isotopic compositions by bootstrapping with replacement over 1000 iterations.

Statistical analyses

Proxy records of past climates and environments are almost always sampled at irregular time intervals due to changes in sedimentation or growth rates through time. To be statistically compared, such irregularly sampled time series usually need to be interpolated and resampled, which introduces artificial and uncertain data points, especially for highly irregular time series (Rehfeld et al., 2011). Commonly applied regression methods then compare these resampled records and provide a single correlation coefficient representative of the relationship between the data set for the whole period of overlap. To account for the complexity of dynamic climatic and environmental systems, we employ two methods to assess the relationships between our new stable isotope records and a range of proxies representative of global and regional changes in temperature, rainfall amounts, and insolation. The first method is a Gaussian kernel-based cross-correlation method (Rehfeld et al., 2011). The advantage of this method is that it does not require interpolation and resampling of the original time series and thus avoids the introduction of uncertainties associated with resampling. Like many traditional methods of regression, however, this method only provides a single correlation value between the two time series. It thus does not account for the inherent dynamics in the climate system and the possibility of changes in the driving forces of the stable isotopes through time. The second method of comparison we apply is semblance analysis (Cooper and Cowan, 2008). This method may introduce some uncertainty into the comparisons of the time series, because it does require the time series to be compared to be sampled to the same time step and resolution. By using a continuous wavelet transform, the semblance comparison of two time series allows the comparison of local phase relationships and frequency contents of two data sets. This means that it provides measures of correlation between time series through time and at different frequencies. Combining these two methods allows us to assess larger-scale relationships of our stable isotope records as a whole against other proxies and to estimate possible changes in these relationships through time.

Kernel-based cross correlation

The Gaussian kernel-based method of cross-correlation analysis (gXCF) of Rehfeld et al. (2011) allows the analyses of unevenly spaced time series without interpolation and resampling. Instead, the method calculates Pearson's cross correlation between measured data points using a function that weighs them depending on their temporal offset. It has been demonstrated that this approach reduces bias in the correlation analysis compared with interpolated and resampled data sets, especially for data sets with very irregular sampling (Rehfeld et al., 2011). The kernel-based cross correlation for two data sets x and y is given by:

$$\hat{\rho }_{x, y} = \displaystyle{{\mathop \sum \nolimits_i^{Nx} \mathop \sum \nolimits_j^{Ny} x_iy_ib( {t_j^y -t_i^x } ) } \over {\mathop \sum \nolimits_i^{Nx} \mathop \sum \nolimits_j^{Ny} b( {t_j^y -t_i^x } ) }}$$

where $b( {t_j^y -t_i^x } )$ is the kernel, determining how much weight to give to the product of the two observations x_i and y_j, based on the time offset between them. The kernel is

$$b = \displaystyle{1 \over {\sqrt {2\pi \sigma } }}e^{-{\vert d \vert }^2/2\sigma ^2}$$

where d is the temporal distance between the observations and σ is the standard deviation of the kernel distribution. Following Rehfeld et al. (2011), we use σ= Δt^xy/4.

Before cross-correlation analyses, we detrended all time series using a Gaussian kernel smoother (Rehfeld et al., 2011) with a bin width of one-third of the respective data set length. Confidence intervals of the correlation coefficients were calculated using 25 repetitions of the bootstrapping protocol of Roberts et al. (2017) with 2000 resampling steps. This protocol accounts for persistence in the original data sets. In cases where the median confidence interval of the 25 repetitions does not span 0, we consider the correlation significant at the 95% confidence level. The composite stable isotope record has a hiatus, and therefore cross-correlation analyses were done separately for the two sections.

Semblance

We used the methods of Cooper and Cowan (2008) for semblance analyses of our records compared with other proxies. The method uses a continuous wavelet transform to allow the comparison of local phase relationships and frequency contents of two data sets. Data sets to be compared by semblance analyses had to be resampled to the same time step and resolution. To avoid oversampling low-resolved data sets, we first calculated the average resolution of all data sets. We then resampled the SK1 stable isotope records and the data set for comparison to the lower average resolution between the two in R using the function approxfun (R Core Team, 2022), rounding to the next lower 0.1 ka time step. For example, if the lower average resolution of the two data sets was 0.12 ka, we resampled both at 0.2 ka time steps. Semblance analyses were performed on the resampled data sets using the MATLAB implementation of the method separately for the two sections of the isotope records before and after the hiatus. Because precession has been previously named as an important orbital parameter affecting climate in this region (Partridge et al., 1997; Chase, 2021), we focus on the precession band in our discussion of the results of semblance analyses. For this, we averaged semblance results for the periods between 15 and 25 ka.

Results

Dating results

Measured ages used for age model construction in this study range between 120.5 ± 1.3 ka (428467-A) and 13.1 ± 0.3 ka (162807-B; Fig. 3, Supplementary Table 2). Of the 59 age analyses, 7 had errors above 5% of the age value (Supplementary Table 2). Ages are plotted against depths, and age models used for stable isotope analyses are shown in Figure 3. Images of the samples and drilling locations of dating analyses and stable isotopic records are shown in Supplementary Figure 3.

Figure 3. Measured ages for samples from Sandkraal Cave (SK1) with age models and 95% probability intervals constructed using Bacon (Blaauw and Christen, 2011). Outlier ages are shown in gray. (a) 162807, (b) 357387-A, (c) 357387-B1, (d) 357387-B2, (e) 428464-67, (f) 460362-1, and (g) 460362-2. The locations of hiatuses in samples 162807 and 357387-A are indicated by horizontal lines.

Petrography

Thin sections of two of the samples (162807, 357387-A) showed elongated columnar fabric with crystals of a few centimeters in length (Supplementary Fig. 5). According to Frisia (2015), such a fabric in speleothems forms from constant drip with a comparatively low supersaturation and intermediate Mg/Ca ratios.

Stable isotope changes through time

Results of stable isotope analyses for individual samples are given in Supplementary Table 1 and plotted in Figure 4. Composite records constructed from them with bootstrapped confidence intervals are plotted in Figure 5, and data are given in Supplementary Table 3. Overall values of δ¹⁸O in speleothems from SK1 range between −4.5 and −1.8‰, mean values of individual samples are between −3.5 and −2.5‰, with standard deviations ranging between 0.27 and 0.41‰. Values of δ¹³C range between −12.5 and −5.7‰, individual sample means are between −9.5 and −8.1‰, with standard deviations between 1.28 and 0.55‰. Some of the samples from SK1 have layers of lighter-colored, laminated, and porous tufa in addition to brown speleothem (Supplementary Fig. 3). Tufa layers can usually be distinguished from speleothem by color/appearance, and a hiatus often separates the deposition from the darker brown layers. Because tufa layers usually form when the cave is open, they were excluded from the composite record and z-score normalization. While there is considerable overlap between the range of δ¹⁸O and δ¹³C values analyzed in tufa and brown speleothem, the values in tufa are generally offset toward higher values (δ¹⁸O: −4.67 to 1.73‰; δ¹³C: −9.97 to 6.76‰), probably because faster degassing and evaporation in the open cave preferentially removed lighter isotopes from the solution.

Figure 4. δ¹⁸O (a) and δ¹³C (b) and measured ages with 2σ uncertainty (c) of individual samples from Sandkraal Cave (SK1). Numbers and letter at the top of the panel denote Marine Isotope Stages (MIS) and substages.

Figure 5. Five-point running mean composite of z-scored δ¹⁸O (a) and δ¹³C (b) records from Sandkraal Cave (SK1, black lines). Individual analyses are shown in gray; bootstrapped 5th and 95th percentiles of the composite record are shown in blue. Numbers and letters at the top of the panel denote Marine Isotope Stages (MIS) and substages.

The composite stable isotopic record from SK1 covers the time interval between 104 and 18 ka with a break in deposition between 48 and 41 ka. Values of the z-scored smoothed δ¹⁸O and δ¹³C values range between 2 SD below the mean and almost 3 SD above the mean values (Fig. 5, Supplementary Table 3).

The oldest section of the records in MIS 5c (~104–92 ka) shows considerable variability in both isotopes. Overall, δ¹⁸O z-scores tend to be below average between the beginning of the record at 103.6 ka and 96.8 ka followed by above-average z-scores until about 93.7 ka. δ¹⁸O z-scores drop slightly to below average just before the transition from MIS 5c into MIS 5b. δ¹³C z-scores overall follow a decreasing trend across MIS 5c—with mostly above-average scores between the beginning of the record at 103.6 ka and 97.5 ka followed by more common scores below average. Both δ¹⁸O and δ¹³C show increasing z-scores at the transition between MIS 5c and MIS 5b, and values are above average for much of MIS 5b. The two isotopes also show a decrease in z-scores at the transition from MIS 5b into MIS 5a, but the timing differs. δ¹⁸O decreases from 87.9 ka to 83.4 ka and maintains mostly below-average values for MIS 5a. Decreasing z-scores occur in the δ¹³C record between 85.8 and 77.5 ka. Below-average values of δ¹³C are observed later in MIS 5a and well into MIS 4 (80.3–69.0 ka). The z-scores of δ¹⁸O and δ¹³C tend to be lower than average early in MIS 4 and reach some of the highest values just before the transition into MIS 3 (60.0 ka). During MIS 3, the variability is lower for both isotopes than during MIS 4 and 5. In early MIS 3 (~59–48 ka), δ¹⁸O values vary around the average, whereas δ¹³C values are generally above average. The hiatus in the composite record lasts between ~47.5 and 40.6 ka. In late MIS 3 (~41–29 ka) δ¹⁸O values continue to vary near the average. δ¹³C values indicate a short low at ~39.5 ka followed by an increasing trend. The SK1 records cover the section between 29 and 17.5 ka in MIS 2. The δ¹⁸O z-scores show some variability in this interval, but most values are above average. δ¹³C z-scores are near average at the MIS 3–2 transition and increase to some of the highest values near the youngest part of the records.

Discussion

Isotopic equilibrium

The use of speleothem δ¹⁸O and δ¹³C in paleoclimatology is based on the premise that the fractionation of isotopes between the water/dissolved inorganic carbon and the secondary carbonate is defined by a temperature-dependent fractionation factor, a condition usually referred to as “stable isotopic equilibrium.” Correlation of δ¹⁸O and δ¹³C along the speleothem time series, as is seen in our samples (Supplementary Fig. 6), is often thought to be caused by disequilibrium processes. However the deposition of modern cave calcites commonly occurs out of isotopic equilibrium, and previous studies have shown that this does not exclude a paleoenvironmental interpretation of the records (Mickler et al., 2006; Dorale and Liu, 2009; Fohlmeister et al., 2012; Daëron et al., 2019). In cases like ours, in which multiple samples are available from the same site, the repetition of trends of δ¹⁸O and δ¹³C changes through time in the different samples is a useful criterion to assess whether or not speleothem records represent larger scale climatic and environmental factors as opposed to drip site–specific kinetic effects (Dorale and Liu, 2009). Figure 4 shows a comparison of the raw stable isotope data of the samples included in this study. Overall, there is good overlap between the direction of trends in the samples. Differences between samples of up to 1‰ for δ¹⁸O and around 2‰ for δ¹³C are observed in some sections of the records, especially for speleothems with low sampling resolutions. This degree of in-cave differences of coeval speleothems is a common feature for speleothem studies (e.g., Hellstrom et al., 1998; Dorale and Liu, 2009; Koltai et al., 2017), and our method of composite record construction is designed to even out possible drip site–specific effects and highlight overlapping trends in the records.

Comparison to other speleothem records from the region

Speleothem records that cover the same time period as the SK1 record have previously been published from three sites along the Cape south coast (Crevice Cave, PP29, and Herolds Bay Cave; Bar-Matthews et al., 2010; Braun et al., 2019, 2020), and three records have been published from two inland locations (Cango Cave, Efflux Cave, and Cape Fold Composite from both these sites; Talma and Vogel 1992; Braun et al., 2020; Chase et al., 2021). Comparisons of the normalized records (Fig. 6) show that overall trends overlap, albeit with some temporal offsets. Most of the δ¹⁸O records decrease from MIS 5b to MIS 5a, and the lowest δ¹⁸O values are often recorded in MIS 5a or early MIS 4 (Fig. 6). Speleothem δ¹⁸O values tend to be above average in late MIS 4 and then often drop to below average in early MIS 3 (Fig. 6). Another common pattern is a gradual increase from this early MIS 3 low of δ¹⁸O values through MIS 2 (Fig. 6). Overall similar patterns are also seen in δ¹³C records from the region: above-average values tend to occur in MIS 5b and MIS 4, whereas MIS 5a has below-average values (Fig. 6). The increase of δ¹³C at the beginning of MIS 4 appears to happen later at SK1 than at most of the other locations. Major differences are observed between the SK1 δ¹³C record and the Cape Fold Composite by Chase et al. (2021). This might be related to the way that the Cape Fold Composite was constructed to include δ¹³C records from multiple cave sites and different proxies (speleothems from Talma and Vogel 1992; Chase et al., 2021; Braun et al., 2020; and hyrax middens from Chase et al., 2017), which required a reduced major axis transformation to align the various records (Chase et al., 2021).

Figure 6. Comparison of the Sandkraal Cave (SK1) δ¹⁸O and δ¹³C to previously published records from the region. Records were normalized as z-scores (if they were not published as such). The records from Crevice Cave and Efflux Cave were smoothed with a 3-point running mean to make the trends more easily comparable. Numbers and letters at the top of the panels denote Marine Isotope Stages (MIS) and substages.

Correlations with other records from the region and factors influencing δ¹⁸O and δ¹³C at SK1

To assess the main influences on our new δ¹⁸O and δ¹³C records, we used two statistical methods to analyze relationships between our records and proxies for global glacial–interglacial change, regional temperatures, summer rainfall in southern Africa, and runoff and upwelling along the southern African west coast. We consider this approach to be more reliable for revealing larger-scale relationships than selecting single records for comparisons. The data sets, references, and acronyms used for comparison are listed in Table 1. We include the following categories and specific records: for global proxies we used a temperature record from Antarctic ice core EPICA DomeC (DomeC T; Jouzel et al., 2007) and the δ¹⁸O of the NGRIP ice core from Greenland (NGRIP; NGRIP, 2004; Bazin et al., 2013) as well as atmospheric CO₂ concentrations (CO₂; Bereiter et al., 2012) and relative sea-level records (rsl; Rohling et al., 2009, 2010). For regional temperature records, we included sea-surface temperature reconstructions from the coasts surrounding southern Africa, including the eastern Cape Basin off the southwestern tip of Africa (Cape Basin SST; Marino et al., 2013; Dyez et al., 2014) and three records forming a cross section from shallow to deep water across the Namibian shelf (Namibia SST off, Namibia SST inter, Namibia SST near; Kirst et al., 1999). Namibia SST near was constructed closest to the coast, and this temperature record is affected by changes in upwelling bringing water masses of different temperatures into this area. Additional regional temperature records include reconstructions from the Mozambique Channel (Mozambique Tex^H₈₆, Mozambique U^k’₃₇; Caley et al., 2011). These records have a low resolution, and cross-correlation analyses for the younger section of our composite record failed to give reliable results. Another temperature record is the pollen-based reconstructions of land temperatures for the last ~45 ka using the CREST method (CREST Tann; Chevalier and Chase, 2015). Proxy reconstructions for precipitation in the present-day summer rainfall region were also included to investigate influences of changing rainfall seasonality at SK1. The records include pollen-based reconstructions using the CREST method for separate northern and central eastern data sets (CREST PWetQ n, CREST PWetQ ce; Chevalier and Chase, 2015). Because the CREST records only cover the last ~45 ka, we were only able to compare them with the younger section of our composite records. Other summer rainfall indicators include proxies for the strength of chemical weathering in watersheds near the east coast (Fe/K 10-06P, Fe/K 17-17K; Ziegler et al., 2013; Simon et al., 2015), δ¹⁸O of Chinese speleothems as a proxy for East Asian monsoon strength (China spels; Wang et al., 2008), and rainfall reconstructions from Tswaing crater based on sediment composition and deuterium isotopes in plant waxes (Tswaing rainfall, Tswaing δD; Partridge et al., 1997; Schmidt et al., 2014). Well-resolved records of winter rainfall in western South Africa that cover a large part of the time period of the SK1 speleothems are harder to find. We therefore use a combination of proxies with some caveats concerning their interpretation. Fluxes of Na⁺ ions to Antarctic EPICA Dronning Maud Land (EDML) ice core are a proxy for the extent of sea ice surrounding Antarctica (EDML Na flux; Fischer et al., 2007). In southern Africa, this has been interpreted to affect the position of the austral westerlies and the amount of winter rainfall, with more sea ice (higher Na flux) shifting the westerlies north and causing higher winter rainfall (Chase et al., 2021). While we do not reject this interpretation, the EDML Na flux record mainly appears to be a reverse record of the global temperature proxies mentioned earlier, so correlation analyses generally showed opposite values/trends through time compared with the global proxies. We include two records of sediment fractions most likely reflecting river runoff (Terr MD08-3167, Terr MD96-2094; Stuut et al., 2002; Collins et al., 2014). While these have been interpreted to reflect more winter rainfall when more fluvial sediment is detected, they might record a more complex signal if sediments from the Orange and Cunene Rivers reach the coring locations. We also include records of total organic carbon on the Namibian shelf that record the strength of upwelling and are related to the intensity of the trade winds (Namibia TOC near, Namibia TOC inter, Namibia TOC off; Kirst et al., 1999). Finally, we also included comparisons to the precession parameter and summer (December) insolation at 30°S in our analyses (Precession, Insol Dec 30S; Laskar et al., 2004).

Table 1. Summary of proxy records used in statistical comparisons with the Sandkraal Cave 1 (SK1), South Africa, stable isotope records.

We have used two methods to compare the proxy records mentioned above to the new δ¹⁸O and δ¹³C records from SK1. The advantage of kernel-based cross-correlation analyses (Rehfeld et al., 2011) is that irregularly sampled time series can be compared directly. However, the method generally does not account for possible changes in the relationships between climate proxy records through time. We separated our composite record at its hiatus and therefore do have independent cross-correlation values for the two growth periods of the speleothems that are plotted in Figure 7. Wavelet-based semblance analyses (Cooper and Cowan, 2008) reveal changes in the relationship between proxies through time and on various periodicities of variation. For this method, we focus mainly on periods related to precession (15–25 ka period). The proxies that are compared do need to be resampled for semblance analysis, which introduces uncertainty. The results are plotted separately for δ¹⁸O and δ¹³C in Figures 8 and 9, respectively.

Figure 7. Pearson's correlation coefficients between Sandkraal Cave (SK1) speleothem δ¹⁸O (a) and δ¹³C (b) records and other paleoclimate proxies using Gaussian kernel-based cross-correlation analyses (Rehfeld et al., 2011). The 95% confidence intervals were calculated using 25 repetitions of a bootstrapping process with 2000 resampling steps for each repetition (Roberts et al., 2017). Correlation coefficients with confidence intervals that span 0 are not considered significant. MIS, Marine Isotope Stage.

Figure 8. Average semblance values for the periods between 15 and 25 ka for comparisons of the Sandkraal Cave (SK1) δ¹⁸O record to global proxies (a), regional temperature proxies (b), precipitation proxies from the summer rainfall region with precession and December insolation at 30°S (c), and proxies of rainfall and upwelling from the west coast (d).

Figure 9. Average semblance values for the periods between 15 and 25 ka for comparisons of the Sandkraal Cave (SK1) δ¹³C record to global proxies (a), regional temperature proxies (b), precipitation proxies from the summer rainfall region with precession and December insolation at 30°S (c), and proxies of rainfall and upwelling from the west coast (d).

To a first order, speleothem δ¹⁸O is affected by temperature in the cave and δ¹⁸O of drip water, which is in turn related to rainwater δ¹⁸O. Influences on rainwater δ¹⁸O can be as varied as continentality, moisture source, rain amount, transport distance, and more. Evaporation from the soil and seasonal biases in infiltration can also affect speleothem δ¹⁸O.

Our δ¹⁸O record shows negative or no cross-correlation for both deposition periods with most proxies that are related to global and regional temperature change included in our comparison (although many of these correlations are not considered statistically significant; Fig. 7a). Kernel-based correlation analyses suggest a tendency toward positive correlation between SK1 δ¹⁸O and summer rainfall region precipitation for MIS 5–3 (Fig. 7a). The section of the composite record covering MIS 3–2 shows more variable correlations with a tendency to negative values (Fig. 7a). This suggests a change in the relationship between δ¹⁸O and summer rainfall through time and will be discussed later in the article. The correlation of SK1 δ¹⁸O with precession, summer (December) insolation at 30°S proxies from the southern African west coast, and EDML Na flux are more often positive than negative, but mostly not statistically significant (Fig. 7a).

Based on the results of the semblance analyses plotted in Figure 8, there is a change in the driving forces that cause variations of δ¹⁸O values at SK1 that coincide with the boundary from MIS 5 to MIS 4 (~70 ka ago). Chase et al. (2021) observed a similar shift at this time in the Cape Fold Composite record from the Little Karoo. Before ~70 ka, the SK1 δ¹⁸O record has negative relationships with proxies for global change and regional temperatures along with positive relationships with hydrological proxies from the summer rainfall region, sea ice extent (position of the westerlies), and summer insolation/precession (Fig. 8a–c). Relationships with runoff on the west coast and upwelling seem to be variable between 103 and 85 ka and negative between 85 and 65 ka (Fig. 8d).

The negative correlation between δ¹⁸O and temperature in MIS 5 (Figs. 7a and 8a and b) reflects the impact of temperature variations on the fractionation of oxygen isotopes between cave drip water and carbonate. Recent analyses of global data from cave monitoring have shown that drip water δ¹⁸O is affected by evaporation and recharge biases in areas where mean annual temperature (MAT) is above 10°C (Baker et al., 2019). Temperature modeling and proxy reconstructions suggest that MAT at SK1 was above this threshold during the LGM and thus probably throughout the deposition periods of the speleothems included here (Heaton et al., 1986; Engelbrecht et al., 2019; Göktürk et al., 2023). The negative correlation between SK1 and temperature records in MIS 5 (Figs. 7a and 8a and b) makes a recharge bias more likely than an evaporation bias at the site, meaning that SK1 δ¹⁸O values mainly reflect the main rainfall season or strong rainfall events. The positive relationship with proxies of precipitation amounts in the summer rainfall region in MIS 5 (Figs. 7a and 8c) suggests that higher SK1 δ¹⁸O values are associated with increases in summer rainfall on the Cape south coast that are in phase with the summer rainfall area, as has been previously suggested (Bar-Matthews et al., 2010; Braun et al., 2019, 2020). Considering the interglacial boundary conditions of MIS 5, the simultaneous positive relationship with the extent of Antarctic sea ice, summer insolation at 30°S, and proxies of summer rainfall (Fig. 8 c) indicate that complex atmospheric interactions govern changes in SK1 δ¹⁸O. Under present-day conditions, tropical temperate troughs are major sources of summer rainfall, especially in the South African interior. These atmospheric systems are formed when a continental summer low-pressure system links up with the low-pressure waves of the westerlies, and such connections are more likely to occur when the westerlies are in a northern position (higher EDML Na flux) and summer insolation is high. Semblance analyses indicate a strong positive relationship between SK1 δ¹⁸O and δ¹⁸O of Chinese speleothems before 67 ka (not shown). The positive relationship suggests an anticorrelation between rainfall amounts in the Chinese East Asian monsoon region (higher δ¹⁸O there reflects less rainfall) and South African summer rainfall (higher SK1 δ¹⁸O = more summer rainfall) during MIS 5. There is an ongoing debate over the relationship between South African monsoonal rainfall and that of the Northern Hemisphere and how it has changed through time. Our record suggests that the speleothem δ¹⁸O on the South African south coast for much of MIS 5 was influenced by biases toward the main rainfall season and tropical–temperate connections led to increases in rainfall in phase with summer rainfall region precipitation that were out of phase with Northern Hemisphere monsoon strength.

Between 70 and 48 ka (MIS 4–3) many of the relationships mentioned earlier are reversed. Higher values of SK1 δ¹⁸O are now associated with higher values in global change proxies, higher temperatures, less summer rainfall, a southward shift of the westerlies, weaker trade winds, and lower summer insolation (Fig. 7a). Similar changes at this time have been attributed to a change in the climate drivers affecting southern Africa between tropical (before 70 ka) and high-latitude (after 70 ka) forces (Chase, 2021; Chase et al., 2021). This shift is accompanied by an overall decrease of temperatures, increase in ice volume, and a decrease of rainfall amounts in the summer rainfall area at the MIS 5–4 transition (Partridge, 1997; Fischer et al., 2007; Jouzel et al., 2007). Higher δ¹⁸O values in MIS 4–3 are related to overall drier conditions (higher temperatures, less summer rainfall, southward position of the westerlies [less winter rainfall and fewer tropical temperate troughs]; Fig. 8), indicating that evaporation biases between rainwater and cave water are likely to affect speleothem δ¹⁸O at this time (Baker et al., 2019). An effect of hydrology, mainly through PCP, on speleothem δ¹⁸O has recently also been inferred for a speleothem record from Bloukrantz Cave (~170 km west of SK1) on the South African south coast during MIS 3 (Maccali et al., 2023). The shift of semblance relationships observed at ~70 ka indicates that the main influence on δ¹⁸O switches from infiltration bias in MIS 5 to evaporation bias in MIS 4–3. This means that rather than mainly reflecting the main rainy season, higher δ¹⁸O values are now associated with the dry season.

After 41 ka (MIS 3–2), semblance analysis results become less uniform (Fig. 8). At the end of the hiatus, lower δ¹⁸O values are associated with globally cooler conditions, but this relationship changes in MIS 2, when globally cooler conditions are associated with higher δ¹⁸O values (Fig. 8a). Semblance relationships between δ¹⁸O and regional temperatures are very variable and inconsistent (Fig. 8b). Comparisons with the summer rainfall region proxies across the MIS 3–2 growth period show mostly lower δ¹⁸O associated with higher rainfall for inland proxy records (CREST PWetQn, CREST PWetQce, Tswaing rainfall, Tswaing δD), whereas relationships with coastal summer rainfall records (Fe/K 10-06P, Fe/K 17-17K), Chinese speleothem δ¹⁸O, precession, and Insol Dec 30S are the opposite (Figs. 7a and 8c). Semblance comparisons with proxies from the west coast are also mostly positive in MIS 3–2 (Fig. 8d). While there is a lot of uncertainty due to the variable relationships of δ¹⁸O with other records, some of these results may support climate simulations of the LGM that suggest an overall increase of winter rainfall in large parts of the west and center of southern Africa (Engelbrecht et al., 2019). Along the south coast, the same model indicates a decrease of winter rainfall related to downwind effects along the coastal mountain ranges and the increased occurrence of dry Berg winds (Engelbrecht et al., 2019). This could lead to evaporative enrichment of cave water δ¹⁸O and the positive correlation with proxies of global change and west coast runoff in MIS 2 (Fig. 8a and c). Prior calcite precipitation is also a possible driver of higher δ¹⁸O during drier phases in this time interval. The highly variable relationship with temperature and summer rainfall region records (Fig. 8b and c), however, possibly contradict this, making the interpretation of this section of the record complex.

Overall the suggested influences on speleothem δ¹⁸O at SK1 change from higher values associated with higher rainfall amounts of isotopically heavy summer rain in MIS 5 to higher values due to evaporation in dry periods of MIS 4 and probably also MIS 3–2.

The main influences on δ¹³C of speleothems are the δ¹³C values of soil CO₂ and the carbonate source material, as well as in-cave processes during CO₂ degassing and carbonate precipitation. At SK1, the carbonate sources are the calc-silicate skarn at the contact of the Rooiklip granite-gneiss with the metasediments of the Victoria Bay Formation and calcareous quartzites within the Victoria Bay Formation. Although the δ¹³C value of these source materials has not been analyzed, a good overlap of the speleothem δ¹³C values measured at SK1 with other sites in the region (Fig. 6) suggests that speleothem deposition at SK1 happened under at least partly open system conditions and/or that the source material had a similar δ¹³C to the other locations (−5–1‰). Changes of δ¹³C in our new speleothem records thus are affected more by changing δ¹³C of soil CO₂ than the contribution from the host rock. Other factors influencing changes of speleothem δ¹³C values include: the type of vegetation above the cave (C3, C4) and its density (which is driven by precipitation and temperature), as well as processes of CO₂ degassing and PCP. We evaluate the potential impact of these factors based on the correlation patterns with other proxy records from the region and beyond.

The patterns of correlation of the SK1 δ¹³C record with our selection of global and regional proxies result in different temporal patterns than the ones seen for δ¹⁸O. Notably, the considerable change of the factors driving δ¹⁸O variations at the MIS 5–4 transition is not found for δ¹³C (Fig. 9a). This means that changes of δ¹³C in MIS 5 through MIS 4-3 were probably driven by consistent influences. The largest changes in the relationships of δ¹³C with other records are seen during the transition from the first (MIS 5–3) to the second deposition phase (MIS 3–2; Fig. 9).

In MIS 5–3, kernel-based correlations with global temperature-related proxies and regional temperature reconstructions indicate variable relationships (Fig. 7b). Semblance comparisons on the precession band, however, suggest mostly negative relationships (Fig. 9a and b). Kernel-based correlation and semblance comparisons agree that relationships of SK1 δ¹³C with summer rainfall proxies, precession, and summer insolation in MIS 5–3 are mostly negative, whereas relationships with west coast proxies tend to be positive (Figs. 7b and 9c and d). Speleothem δ¹³C values are most commonly interpreted as proxies for effective moisture availability, with higher values generally reflecting drier conditions with increased PCP, enhanced CO₂ degassing from slower dripping of water in the cave, decreased soil bioproductivity and vegetation density, and, in regions where C3 and C4 plants overlap, a possible increase of C4 grasses. The observed negative relationship with proxies of summer rainfall in MIS 5–3 (Fig. 9 c) supports such an interpretation for SK1 δ¹³C values as well, provided that low rainfall in the summer rainfall region is also associated with drier conditions on the south coast. The negative relationships between SK1 δ¹³C and rsl in MIS 5–3 (Figs. 7b and 9a), in light of the modeling results of Göktürk et al. (2023) showing increased climatic continentality when the PAP is exposed, also supports the interpretation of SK1 δ¹³C in terms of water availability. The GCFR hosts a mix of C3 and C4 grasses, especially in sections with substantial contributions of warm-season rainfall along the Cape coastal lowlands. Drier conditions in these sections of the GCFR and on the exposed PAP may also be associated with an increase of C4 grass proportions in the understory, particularly during times when atmospheric CO₂ concentrations were decreased (Grobler et al., 2023). Consistently negative kernel-based correlations of our δ¹³C record with atmospheric CO₂ concentrations also support this additional interaction between effective moisture, vegetation composition, and δ¹³C values in speleothems.

Some of the relationships to other proxies change for the section of the record covering MIS 3–2: kernel-based correlation with global temperature–related proxies become consistently negative (although not all statistically significant), whereas semblance comparisons indicate a positive relationship on precessional timescales (Figs. 7b and 9a). The relationships with regional temperature records indicate some variability in MIS 3–2, but semblance values tend to be positive (Figs. 7b and 9b). Kernel-based correlation analyses tend to be negative for inland summer rainfall records and positive for coastal records (Fig. 7b). Semblance comparisons for MIS 3–2 also show negative relationships of SK1 δ¹³C with inland summer rainfall records, whereas coastal records seem to not have a clear relationship with δ¹³C (semblance values near 0; Fig. 9c). Kernel-based correlations and semblance comparisons of SK1 δ¹³C with precession, insolation, Chinese speleothems, and Antarctic sea ice extent all indicate positive relationships in MIS 3–2 (Figs. 7b and 9c). While kernel-based correlation analyses of the SK1 δ¹³C against records from the west coast show variable results for this time period, semblance comparisons suggest that a positive relationship is maintained on precessional timescales (Figs. 7b and 9d). The continued tendency toward negative correlations between inland summer rainfall records and SK1 δ¹³C (Figs. 7b and 9c) may indicate that effective moisture availability still plays a role in driving δ¹³C values in MIS 3–2. Positive semblance relationships with global and regional temperature proxies, precession, and insolation (Fig. 9 a–c) might reflect the influence of higher temperatures on effective rainfall and SK1 δ¹³C on the precessional band, but the variable kernel-based correlation (Fig. 7 b) values may mean that for other periods of variability, in-cave and disequilibrium processes also affect δ¹³C. Overall, the δ¹³C values of the SK1 record most likely reflect changes in effective moisture availability, its effects on the soils and vegetation above the cave, and especially during the MIS 3–2 deposition period, in-cave processes.

Climate variations at SK1 through time

Comparisons of the composite δ¹⁸O and δ¹³C records from SK1 with representative proxy records are shown in Figures 10 and 11, respectively.

Figure 10. Comparison of the Sandkraal Cave (SK1) δ¹⁸O isotopic profile to other key climate proxies. (a) Temperature reconstruction for Antarctic Dome C (relative to present day; Jouzel et al., 2007); (b) sea-surface temperature (SST) reconstruction for marine sediment core GeoB-1711-4 from intermediate water depth on the Namibian shelf (Kirst et al., 1999); (c) Fe/K ratios of sediment core CD154 17-17K from the east coast of South Africa (Ziegler et al., 2013); (d) δ¹³C of the SK1 speleothems (this study); (e) total organic carbon content of nearshore sediment core GeoB 1712-4 from the Namibian shelf (Kirst et al., 1999); (f) precession parameter (Laskar et al., 2004).

Figure 11. Comparison of the Sandkraal Cave (SK1) δ¹³C isotopic profile to other key climate proxies. (a) Atmospheric CO₂ concentrations (Bereiter et al., 2015); (b) sea-surface temperature (SST) reconstruction for marine sediment cores MD02-2594 and MD96-2080 from the Cape Basin using Globigerinoides ruber (Marino et al., 2013; Dyez et al., 2014); (c) Fe/K ratios of sediment core CD154 17-17K from the east coast of South Africa (Ziegler et al., 2013); (d) δ¹⁸O of the SK1 speleothems (this study); (e) total organic carbon content of nearshore sediment core GeoB 1712-4 from the Namibian shelf (Kirst et al., 1999); (f) precession parameter (Laskar et al., 2004).

MIS 5

MIS 5 shows some of the largest variability in the SK1 δ¹⁸O and δ¹³C records. The highest values of both isotopes are recorded near the transition from MIS 5d into MIS 5c at the beginning of the record and during MIS 5b (Figs. 10 and 11). The lowest δ¹⁸O and δ¹³C values of the SK1 composite records occur in MIS 5c and MIS 5a. We interpret the close positive relationship between δ¹⁸O and precipitation in the summer rainfall region to mean that summer rainfall was also an important source of moisture on the south coast at this time. Periods with the highest summer rainfall (highest δ¹⁸O) also have some of the lowest δ¹³C values, probably supporting dense C3 vegetation and reducing PCP and fast CO₂ degassing from fast dripping water, whereas during phases with less summer rainfall, slower dripping in the cave enhanced CO₂ degassing, PCP was increased, vegetation was sparser, and C4 grass was more common in the understory.

MIS 4–3

The boundary between MIS 5 and MIS 4 marks the transition from interglacial to glacial conditions, with boundary conditions changing considerably—overall temperatures are lower (e.g., Kirst et al., 1999; Jouzel et al., 2007), atmospheric CO₂ concentrations and sea levels decrease (Bereiter et al., 2012; Rohling et al., 2009, 2010), and rainfall in the South African summer rainfall region decreases and shifts away from a precession pacing (Partridge, 1997). This change has also been associated with a shift from low- to high-latitude climate drivers in southern South Africa (Chase et al., 2021). The relationship of SK1 δ¹⁸O with temperature records after 70 ka shifts to positive, and the relationship with precession/insolation changes to a negative one (Fig. 7).

The composite δ¹³C record shows an overall increase from early MIS 4 through early MIS 3, which overlaps with a decrease of atmospheric CO₂ concentrations, an increase of regional temperature, and a decrease of summer rainfall (Fig. 11a–c), possibly reflecting the impact of changing CO₂ on vegetation (overall increase of δ¹³C in C3 plants, shift to more C4 plants) and decreased effective moisture on PCP and CO₂ degassing in the cave. A peak of δ¹⁸O at the MIS 4–3 transition overlaps with comparatively high temperatures within the glacial stage, a peak in trade wind strength, and low summer rainfall (Fig. 10a, c, and e), suggesting an effect of evaporation. Relatively high δ¹³C values at the MIS 4–3 transition also mean that PCP is a possible influence on δ¹⁸O. Lower deposition rates for the speleothems in SK1 during MIS 3 than in the previous intervals and the hiatus in deposition may also support drier conditions. A similar interpretation in terms of changes of PCP has also been proposed for a nearby speleothem from Bloukrantz Cave (Maccali et al., 2023).

MIS 2

The stable isotope records from MIS 2 at SK1 are dominated by a steep increase of δ¹³C and, to a lesser degree, δ¹⁸O. This probably still reflects dry conditions and PCP/evaporation and supports proxies suggesting low summer rainfall amounts, low temperatures, and high upwelling/trade wind strength (Figs. 10 and 11). Both δ¹⁸O and δ¹³C also show a positive relationship with Antarctic sea ice extent/the latitudinal position of the westerlies (Figs. 8 and 9), which supports LGM paleoclimate modeling indicating decreased winter rainfall near SK1 when the westerlies are shifted north (Engelbrecht et al., 2019). Dry conditions at SK1 could also be related to low sea levels and continentality on the exposed PAP (Göktürk et al., 2023).

Conclusions

We present new records of speleothem composite δ¹⁸O and δ¹³C values from Sandkraal Cave (SK1) on the Cape south coast, covering the time interval between 104 and 18 ka, with a hiatus from 48 to 41 ka. Kernel-based correlation analyses and semblance comparisons based on continuous wavelet transforms inform the relationships between our records and other proxies for the two deposition periods as a whole and their changes through time, respectively. Correlation analyses show that there is a relationship between the temperature at the cave site and speleothem δ¹⁸O, but additional variation is added by changes in rainfall seasonality and rainfall δ¹⁸O. Correlations of the δ¹³C record with other proxies are more complex, with no clear relationship with temperature records. Inland summer rainfall region records tend to be negatively correlated with δ¹³C, whereas coastal records show variable correlation. Correlations between δ¹³C and precession/summer insolation change from negative in MIS 5–3 (deposition before the hiatus) to positive in MIS 3–2 (after the hiatus)

Semblance analyses further show that there are changes in the factors driving δ¹⁸O and δ¹³C variations through time. Between the beginning of the record at 105 and ~70 ka, variation of δ¹⁸O at SK1 is most likely related to rainfall seasonality, with higher values generally associated with more summer rainfall probably from increased interactions between tropical and temperate circulation systems and higher frequencies of tropical temperate troughs. Values of δ¹³C are generally lower when summer rain is high, suggesting that more summer rainfall was associated with denser vegetation and less summer rainfall and drier conditions increased PCP, CO₂ degassing in the cave, and possibly C4 grass abundance in the vegetation and decreased vegetation density and soil bioproductivity. Between ~70 and 48 ka, the semblance values between SK1 δ¹⁸O and other proxies change, and now higher δ¹⁸O is associated with higher temperatures and less rainfall. Semblance values for the comparisons of δ¹³C to other proxies do not see a similar change as those for δ¹⁸O. In general, higher values of both isotopes probably are associated with drier conditions in this interval, reflecting enhanced degassing, evaporation, and PCP. After 42 ka, the semblance relationships change again for δ¹⁸O. Higher δ¹⁸O and δ¹³C values now at times with northward extended westerlies may support findings of paleoclimate modeling, suggesting a decrease of winter rainfall on the South African south coast due to downwind effects along the coastal mountains and increased continentality on the exposed PAP.