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, Reference Tyson1986). 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., Reference Dingle, Siesser and Newton1983). 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., Reference Marean, Cawthra, Cowling, Esler, Fisher, Milewski, Potts, Singles, De Vynck, Allsopp, Colville and Verboom2014), and its exposure would have increased the inland continental character of present-day coastal sites (van Andel Reference van Andel1989; Fisher et al., Reference Fisher, Bar-Matthews, Jerardino and Marean2010; Göktürk et al., Reference Göktürk, Sobolowski, Simon, Zhang and Jansen2023).
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, Reference Chevalier and Chase2015); triangles, marine sediment cores (light blue, rainfall proxies [MD08-3167: Collins et al., Reference Collins, Schefuß, Govin, Mulitza and Tiedemann2014; CD154-10-06P: Simon et al., Reference Simon, Ziegler, Bosmans, Barker, Reason and Hall2015; MD96-2094: Stuut et al., Reference Stuut, Prins, Schneider, Weltje, Jansen and Postma2002; 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., Reference Caley, Kim, Malaizé, Giraudeau, Laepple, Caillon and Charlier2011; MD02-2594 and MD96-2080: Marino et al., Reference Marino, Zahn, Ziegler, Purcell, Knorr, Hall, Ziveri and Elderfield2013, 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., Reference Bar-Matthews, Marean, Jacobs, Karkanas, Fisher, Herries and Brown2010; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Cowling, Karkanas, Fisher, Dyez, Zilberman and Marean2019, 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, Reference Talma and Vogel1992). This record was recently amended by the publication of stable isotope records of a second sample from Cango Cave (Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). The interpretations by our group have mostly attributed changes of δ18O in speleothems to changes in rainfall seasonality (with higher values reflecting more summer rainfall; Bar-Matthews et al., Reference Bar-Matthews, Marean, Jacobs, Karkanas, Fisher, Herries and Brown2010; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Cowling, Karkanas, Fisher, Dyez, Zilberman and Marean2019, 2020) based on modern records of rainfall isotope seasonality (Braun et al., Reference Braun, Bar-Matthews, Ayalon, Zilberman and Matthews2017). The records from Cango Cave were driven by changes of drip water δ18O and temperature (Talma and Vogel Reference Talma and Vogel1992) or thought to be influenced by too many factors to parse (Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). Speleothem δ13C 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, Reference Talma and Vogel1992; Bar-Matthews et al., Reference Bar-Matthews, Marean, Jacobs, Karkanas, Fisher, Herries and Brown2010; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Cowling, Karkanas, Fisher, Dyez, Zilberman and Marean2019, 2020; Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). A short Marine Isotope Stage (MIS) 3 record from Maccali et al. (Reference Maccali, Meckler, Lauritzen, Brekken, Rokkan, Fernandez, Krüger, Adigun, Affolter and Leuenberger2023) was the first from this region to attribute changes of speleothem δ18O and δ13C 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., Reference Marean, Cawthra, Cowling, Esler, Fisher, Milewski, Potts, Singles, De Vynck, Allsopp, Colville and Verboom2014; Wadley, Reference Wadley2015) as well as the hyper-diverse extratropical flora of the Greater Cape Floristic Region (GCFR; Braun et al., Reference Braun, Cowling, Bar-Matthews, Matthews, Ayalon, Zilberman, Difford, Edwards, Li and Marean2023). 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, Reference Hitchcock2007). 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, Reference Krynauw and Gresse1980). The Rooiklip granites, which intruded the Victoria Bay Formation about 527.2 Ma (Chemale et al., Reference Chemale, Scheepers, Gresse and Van Schmus2011), 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, Reference Krynauw and Gresse1980; Prent, Reference Prent2013).
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., Reference Karkanas, Marean, Bar-Matthews, Jacobs, Fisher and Braun2021). 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., Reference Roberts, Karkanas, Jacobs, Marean and Roberts2012) and is at a similar elevation as several caves at Pinnacle Point (Fisher et al., Reference Fisher, Bar-Matthews, Jerardino and Marean2010; Karkanas et al., Reference Karkanas, Marean, Bar-Matthews, Jacobs, Fisher and Braun2021).
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, Reference Hitchcock2007; 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 δ18O 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, Reference Engelbrecht and Landman2016). Ridges of the South Atlantic Anticyclone can occur year-round and are associated with stratiform rainfall (Engelbrecht and Landman, Reference Engelbrecht and Landman2016). 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., Reference Engelbrecht, Landman, Engelbrecht and Malherbe2015; Engelbrecht and Landman, Reference Engelbrecht and Landman2016). 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, Reference Engelbrecht and Landman2016). This increased contribution of convective rainfall leads to amount-weighted average rainfall δ18O 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., Reference Braun, Bar-Matthews, Ayalon, Zilberman and Matthews2017). 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., Reference Engelbrecht, Marean, Cowling, Engelbrecht, Neumann, Scott and Nkoana2019).
The average annual temperature at SK1 is 16.2°C and mean annual rainfall is 890 mm (Schulze and Lynch, Reference Schulze and Lynch2007; Schulze and Maharaj, Reference Schulze and Maharaj2007). 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., Reference Bergh, Verboom, Rouget, Cowling, Allsop, Colville and Verboom2014; Grobler and Cowling, Reference Grobler and Cowling2022). 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., Reference Procheş, Cowling and du Preez2005; Rebelo et al., Reference Rebelo, Boucher, Helme, Mucina, Rutherford, Mucina and Rutherford2006). 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., Reference Vogel, Fuls and Ellis1978; Cowling, Reference Cowling1983; Cramer et al., Reference Cramer, West, Power, Skelton, Stock, Allsopp, Colville and Verboom2014). Succulent plants using crassulacean acid metabolism (CAM) occur mainly in succulent karoo, strandveld, and subtropical thicket communities (Cowling and Campbell, Reference Cowling and Campbell1983; Cowling et al., Reference Cowling, Esler, Midgley and Honig1994; Rundel et al., Reference Rundel, Esler and Cowling1999).
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., Reference Rebelo, Boucher, Helme, Mucina, Rutherford, Mucina and Rutherford2006; Williams, Reference Williams2015). 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., Reference Rebelo, Boucher, Helme, Mucina, Rutherford, Mucina and Rutherford2006; Williams, Reference Williams2015). 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., Reference Rebelo, Boucher, Helme, Mucina, Rutherford, Mucina and Rutherford2006; Williams, Reference Williams2015).
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., Reference Fisher, Bar-Matthews, Jerardino and Marean2010). 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., Reference Cawthra, Cowling, Andò and Marean2020). 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., Reference Cowling, Potts, Franklin, Midgley, Engelbrecht and Marean2020). 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, Reference Klein, Deacon, Hendey and Lambrechts1983; Marean et al., Reference Marean, Cawthra, Cowling, Esler, Fisher, Milewski, Potts, Singles, De Vynck, Allsopp, Colville and Verboom2014). Fossil trackways of giraffe also suggest a flora more typical of savannah to the east and north (Helm et al., Reference Helm, McCrea, Lockley, Cawthra, Thesen and Mwankunda2018). Recent modeling results suggest that decreased atmospheric CO2 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., Reference Grobler, Franklin, Marean, Gravel-Miguel and Cowling2023). 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. (Reference Vaks, Bar-Matthews, Ayalon, Matthews, Frumkin, Dayan, Halicz, Almogi-Labin and Schilman2006, Reference Vaks, Bar-Matthews, Matthews, Ayalon and Frumkin2010) and Grant et al. (Reference Grant, Rohling, Bar-Matthews, Ayalon, Medina-Elizalde, Bronk Ramsey, Satow and Roberts2012). Powdered aliquots weighing 0.2–0.3 g were dissolved in 7 N HNO3. Detrital material was separated in a centrifuge (10 min at 4000 rpm) and treated with HNO3 and HF. All samples were spiked with a 229Th + 236U 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 HNO3, 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 HNO3. 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 235U/238U 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 234U and 230Th concentrations was possible by isotope dilution analysis using the 229Th + 236U spike. Sample 428464 was also dated at the University of Minnesota. Extraction and purification methods were as described by Edwards et al. (Reference Edwards, Chen and Wasserburg1987) and carried out on a Thermo-Finnigan Neptune MC-ICP-MS following the procedures of Shen et al. (Reference Shen, Wu, Cheng, Edwards, Hsieh, Gallet and Chang2012) and Cheng et al. (Reference Cheng, Edwards, Shen, Polyak, Asmerom, Woodhead and Hellstrom2013).
Previously published decay constants for 234U and 230Th (Cheng et al., Reference Cheng, Edwards, Shen, Polyak, Asmerom, Woodhead and Hellstrom2013) 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 (232Th/238U = 3.8), assuming the initial 230Th/232Th atomic ratio of 4.4 ± 2.2 × 10−6.
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. (Reference Bar-Matthews, Ayalon and Kaufman1997, Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003). 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 CO2 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 (δ13C = 4.74‰; δ18O = −3.85‰) was measured every eight samples. All δ18O and δ13C values obtained were calibrated against the internal lab standard, which was calibrated against the international standards NBS-19 (δ18O = −2.20‰; δ13C = +1.95‰) and NBS-18 (δ18O = −23.01‰; δ13C = −5.01‰) and are reported in per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB). Analytical reproducibility is better than 0.1‰ for both δ18O and δ13C.
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. CO2 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 (δ18O = –2.01‰; δ13C = 2.1‰) were used for sample standardization to VPDB; both standards were analyzed every 12 samples. Beijing marble (δ18O = –10.9‰; δ13C = 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 δ18O and δ13C.
Age modeling and composite record construction
Age models were constructed based on Bayesian statistics using the software Bacon (Blaauw and Christen, Reference Blaauw and Christen2011) 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. (Reference Williams, King, Zhao and Collerson2005), 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 δ13C 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., Reference Rehfeld, Marwan, Heitzig and Kurths2011). 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., Reference Rehfeld, Marwan, Heitzig and Kurths2011). 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, Reference Cooper and Cowan2008). 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. (Reference Rehfeld, Marwan, Heitzig and Kurths2011) 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., Reference Rehfeld, Marwan, Heitzig and Kurths2011). The kernel-based cross correlation for two data sets x and y is given by:
where $b( {t_j^y -t_i^x } )$
is the kernel, determining how much weight to give to the product of the two observations xi and yj, based on the time offset between them. The kernel is
where d is the temporal distance between the observations and σ is the standard deviation of the kernel distribution. Following Rehfeld et al. (Reference Rehfeld, Marwan, Heitzig and Kurths2011), we use σ= Δtxy/4.
Before cross-correlation analyses, we detrended all time series using a Gaussian kernel smoother (Rehfeld et al., Reference Rehfeld, Marwan, Heitzig and Kurths2011) 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. (Reference Roberts, Curran, Poynter, Moy, van Ommen, Vance and Tozer2017) 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 (Reference Cooper and Cowan2008) 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., Reference Partridge, DeMenocal, Lorentz, Paiker and Vogel1997; Chase, Reference Chase2021), 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, Reference Blaauw and Christen2011). 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 δ18O 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 δ13C 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 δ18O and δ13C values analyzed in tufa and brown speleothem, the values in tufa are generally offset toward higher values (δ18O: −4.67 to 1.73‰; δ13C: −9.97 to 6.76‰), probably because faster degassing and evaporation in the open cave preferentially removed lighter isotopes from the solution.
Figure 4. δ18O (a) and δ13C (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 δ18O (a) and δ13C (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 δ18O and δ13C 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, δ18O 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. δ18O z-scores drop slightly to below average just before the transition from MIS 5c into MIS 5b. δ13C 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 δ18O and δ13C 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. δ18O decreases from 87.9 ka to 83.4 ka and maintains mostly below-average values for MIS 5a. Decreasing z-scores occur in the δ13C record between 85.8 and 77.5 ka. Below-average values of δ13C are observed later in MIS 5a and well into MIS 4 (80.3–69.0 ka). The z-scores of δ18O and δ13C 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), δ18O values vary around the average, whereas δ13C 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) δ18O values continue to vary near the average. δ13C 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 δ18O z-scores show some variability in this interval, but most values are above average. δ13C 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 δ18O and δ13C 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 δ18O and δ13C 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., Reference Mickler, Stern and Banner2006; Dorale and Liu, Reference Dorale and Liu2009; Fohlmeister et al., Reference Fohlmeister, Schröder-Ritzrau, Scholz, Spötl, Riechelmann, Mudelsee and Wackerbarth2012; Daëron et al., Reference Daëron, Drysdale, Peral, Huyghe, Blamart, Coplen, Lartaud and Zanchetta2019). In cases like ours, in which multiple samples are available from the same site, the repetition of trends of δ18O and δ13C 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, Reference Dorale and Liu2009). 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 δ18O and around 2‰ for δ13C 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., Reference Hellstrom, McCulloch and Stone1998; Dorale and Liu, Reference Dorale and Liu2009; Koltai et al., Reference Koltai, Spötl, Luetscher, Cheng, Barrett and Müller2017), 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., Reference Bar-Matthews, Marean, Jacobs, Karkanas, Fisher, Herries and Brown2010; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Cowling, Karkanas, Fisher, Dyez, Zilberman and Marean2019, Reference Braun, Bar-Matthews, Matthews, Ayalon, Zilberman, Cowling, Fisher, Herries, Brink and Marean2020), 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 Reference Talma and Vogel1992; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Zilberman, Cowling, Fisher, Herries, Brink and Marean2020; Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). Comparisons of the normalized records (Fig. 6) show that overall trends overlap, albeit with some temporal offsets. Most of the δ18O records decrease from MIS 5b to MIS 5a, and the lowest δ18O values are often recorded in MIS 5a or early MIS 4 (Fig. 6). Speleothem δ18O 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 δ18O values through MIS 2 (Fig. 6). Overall similar patterns are also seen in δ13C 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 δ13C 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 δ13C record and the Cape Fold Composite by Chase et al. (Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). This might be related to the way that the Cape Fold Composite was constructed to include δ13C records from multiple cave sites and different proxies (speleothems from Talma and Vogel Reference Talma and Vogel1992; Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Zilberman, Cowling, Fisher, Herries, Brink and Marean2020; and hyrax middens from Chase et al., Reference Chase, Chevalier, Boom and Carr2017), which required a reduced major axis transformation to align the various records (Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021).
Figure 6. Comparison of the Sandkraal Cave (SK1) δ18O and δ13C 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 δ18O and δ13C at SK1
To assess the main influences on our new δ18O and δ13C 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., Reference Jouzel, Masson-Delmotte, Cattani, Dreyfus, Falourd, Hoffmann and Minster2007) and the δ18O of the NGRIP ice core from Greenland (NGRIP; NGRIP, 2004; Bazin et al., Reference Bazin, Capron and Landais2013) as well as atmospheric CO2 concentrations (CO2; Bereiter et al., Reference Bereiter, Lüthi, Siegrist, Schüpbach, Stocker and Fischer2012) and relative sea-level records (rsl; Rohling et al., Reference Rohling, Grant, Bolshaw, Roberts, Siddall, Hemleben and Kucera2009, Reference Rohling, Braun, Grant, Kucera, Roberts, Siddall and Trommer2010). 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., Reference Marino, Zahn, Ziegler, Purcell, Knorr, Hall, Ziveri and Elderfield2013; Dyez et al., Reference Dyez, Zahn and Hall2014) 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., Reference Kirst, Schneider, Müller, von Storch and Wefer1999). 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 TexH86, Mozambique Uk’37; Caley et al., Reference Caley, Kim, Malaizé, Giraudeau, Laepple, Caillon and Charlier2011). 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, Reference Chevalier and Chase2015). 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, Reference Chevalier and Chase2015). 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., Reference Ziegler, Simon, Hall, Barker, Stringer and Zahn2013; Simon et al., Reference Simon, Ziegler, Bosmans, Barker, Reason and Hall2015), δ18O of Chinese speleothems as a proxy for East Asian monsoon strength (China spels; Wang et al., Reference Wang, Cheng, Edwards, Kong, Shao, Chen, Wu, Jiang, Wang and An2008), and rainfall reconstructions from Tswaing crater based on sediment composition and deuterium isotopes in plant waxes (Tswaing rainfall, Tswaing δD; Partridge et al., Reference Partridge, DeMenocal, Lorentz, Paiker and Vogel1997; Schmidt et al., Reference Schmidt, Oberhänsli and Wilkes2014). 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., Reference Fischer, Fundel, Ruth, Twarloh, Wegner, Udisti and Becagli2007). 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., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). 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., Reference Stuut, Prins, Schneider, Weltje, Jansen and Postma2002; Collins et al., Reference Collins, Schefuß, Govin, Mulitza and Tiedemann2014). 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., Reference Kirst, Schneider, Müller, von Storch and Wefer1999). 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., Reference Laskar, Robutel, Joutel, Gastineau, Correia and Levrard2004).
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 δ18O and δ13C records from SK1. The advantage of kernel-based cross-correlation analyses (Rehfeld et al., Reference Rehfeld, Marwan, Heitzig and Kurths2011) 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, Reference Cooper and Cowan2008) 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 δ18O and δ13C in Figures 8 and 9, respectively.
Figure 7. Pearson's correlation coefficients between Sandkraal Cave (SK1) speleothem δ18O (a) and δ13C (b) records and other paleoclimate proxies using Gaussian kernel-based cross-correlation analyses (Rehfeld et al., Reference Rehfeld, Marwan, Heitzig and Kurths2011). The 95% confidence intervals were calculated using 25 repetitions of a bootstrapping process with 2000 resampling steps for each repetition (Roberts et al., Reference Roberts, Curran, Poynter, Moy, van Ommen, Vance and Tozer2017). 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) δ18O 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) δ13C 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 δ18O is affected by temperature in the cave and δ18O of drip water, which is in turn related to rainwater δ18O. Influences on rainwater δ18O 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 δ18O.
Our δ18O 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 δ18O 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 δ18O and summer rainfall through time and will be discussed later in the article. The correlation of SK1 δ18O 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 δ18O values at SK1 that coincide with the boundary from MIS 5 to MIS 4 (~70 ka ago). Chase et al. (Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021) observed a similar shift at this time in the Cape Fold Composite record from the Little Karoo. Before ~70 ka, the SK1 δ18O 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 δ18O 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 δ18O is affected by evaporation and recharge biases in areas where mean annual temperature (MAT) is above 10°C (Baker et al., Reference Baker, Hartmann, Duan, Hankin, Comas-Bru, Cuthbert and Treble2019). 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., Reference Heaton, Talma and Vogel1986; Engelbrecht et al., Reference Engelbrecht, Marean, Cowling, Engelbrecht, Neumann, Scott and Nkoana2019; Göktürk et al., Reference Göktürk, Sobolowski, Simon, Zhang and Jansen2023). 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 δ18O 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 δ18O 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., Reference Bar-Matthews, Marean, Jacobs, Karkanas, Fisher, Herries and Brown2010; Braun et al., Reference Braun, Bar-Matthews, Matthews, Ayalon, Cowling, Karkanas, Fisher, Dyez, Zilberman and Marean2019, Reference Braun, Bar-Matthews, Matthews, Ayalon, Zilberman, Cowling, Fisher, Herries, Brink and Marean2020). 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 δ18O. 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 δ18O and δ18O 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 δ18O there reflects less rainfall) and South African summer rainfall (higher SK1 δ18O = 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 δ18O 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 δ18O 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, Reference Chase2021; Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). 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, Reference Partridge1997; Fischer et al., Reference Fischer, Fundel, Ruth, Twarloh, Wegner, Udisti and Becagli2007; Jouzel et al., Reference Jouzel, Masson-Delmotte, Cattani, Dreyfus, Falourd, Hoffmann and Minster2007). Higher δ18O 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 δ18O at this time (Baker et al., Reference Baker, Hartmann, Duan, Hankin, Comas-Bru, Cuthbert and Treble2019). An effect of hydrology, mainly through PCP, on speleothem δ18O 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., Reference Maccali, Meckler, Lauritzen, Brekken, Rokkan, Fernandez, Krüger, Adigun, Affolter and Leuenberger2023). The shift of semblance relationships observed at ~70 ka indicates that the main influence on δ18O 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 δ18O 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 δ18O values are associated with globally cooler conditions, but this relationship changes in MIS 2, when globally cooler conditions are associated with higher δ18O values (Fig. 8a). Semblance relationships between δ18O 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 δ18O 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 δ18O, 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 δ18O 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., Reference Engelbrecht, Marean, Cowling, Engelbrecht, Neumann, Scott and Nkoana2019). 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., Reference Engelbrecht, Marean, Cowling, Engelbrecht, Neumann, Scott and Nkoana2019). This could lead to evaporative enrichment of cave water δ18O 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 δ18O 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 δ18O 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 δ13C of speleothems are the δ13C values of soil CO2 and the carbonate source material, as well as in-cave processes during CO2 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 δ13C value of these source materials has not been analyzed, a good overlap of the speleothem δ13C 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 δ13C to the other locations (−5–1‰). Changes of δ13C in our new speleothem records thus are affected more by changing δ13C of soil CO2 than the contribution from the host rock. Other factors influencing changes of speleothem δ13C 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 CO2 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 δ13C record with our selection of global and regional proxies result in different temporal patterns than the ones seen for δ18O. Notably, the considerable change of the factors driving δ18O variations at the MIS 5–4 transition is not found for δ13C (Fig. 9a). This means that changes of δ13C in MIS 5 through MIS 4-3 were probably driven by consistent influences. The largest changes in the relationships of δ13C 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 δ13C 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 δ13C values are most commonly interpreted as proxies for effective moisture availability, with higher values generally reflecting drier conditions with increased PCP, enhanced CO2 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 δ13C 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 δ13C and rsl in MIS 5–3 (Figs. 7b and 9a), in light of the modeling results of Göktürk et al. (Reference Göktürk, Sobolowski, Simon, Zhang and Jansen2023) showing increased climatic continentality when the PAP is exposed, also supports the interpretation of SK1 δ13C 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 CO2 concentrations were decreased (Grobler et al., Reference Grobler, Franklin, Marean, Gravel-Miguel and Cowling2023). Consistently negative kernel-based correlations of our δ13C record with atmospheric CO2 concentrations also support this additional interaction between effective moisture, vegetation composition, and δ13C 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 δ13C with inland summer rainfall records, whereas coastal records seem to not have a clear relationship with δ13C (semblance values near 0; Fig. 9c). Kernel-based correlations and semblance comparisons of SK1 δ13C 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 δ13C 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 δ13C (Figs. 7b and 9c) may indicate that effective moisture availability still plays a role in driving δ13C 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 δ13C 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 δ13C. Overall, the δ13C 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 δ18O and δ13C records from SK1 with representative proxy records are shown in Figures 10 and 11, respectively.
Figure 10. Comparison of the Sandkraal Cave (SK1) δ18O isotopic profile to other key climate proxies. (a) Temperature reconstruction for Antarctic Dome C (relative to present day; Jouzel et al., Reference Jouzel, Masson-Delmotte, Cattani, Dreyfus, Falourd, Hoffmann and Minster2007); (b) sea-surface temperature (SST) reconstruction for marine sediment core GeoB-1711-4 from intermediate water depth on the Namibian shelf (Kirst et al., Reference Kirst, Schneider, Müller, von Storch and Wefer1999); (c) Fe/K ratios of sediment core CD154 17-17K from the east coast of South Africa (Ziegler et al., Reference Ziegler, Simon, Hall, Barker, Stringer and Zahn2013); (d) δ13C 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., Reference Kirst, Schneider, Müller, von Storch and Wefer1999); (f) precession parameter (Laskar et al., Reference Laskar, Robutel, Joutel, Gastineau, Correia and Levrard2004).
Figure 11. Comparison of the Sandkraal Cave (SK1) δ13C isotopic profile to other key climate proxies. (a) Atmospheric CO2 concentrations (Bereiter et al., Reference Bereiter, Eggleston, Schmitt, Nehrbass-Ahles, Stocker, Fischer, Kipfstuhl and Chappellaz2015); (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., Reference Marino, Zahn, Ziegler, Purcell, Knorr, Hall, Ziveri and Elderfield2013; Dyez et al., Reference Dyez, Zahn and Hall2014); (c) Fe/K ratios of sediment core CD154 17-17K from the east coast of South Africa (Ziegler et al., Reference Ziegler, Simon, Hall, Barker, Stringer and Zahn2013); (d) δ18O 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., Reference Kirst, Schneider, Müller, von Storch and Wefer1999); (f) precession parameter (Laskar et al., Reference Laskar, Robutel, Joutel, Gastineau, Correia and Levrard2004).
MIS 5
MIS 5 shows some of the largest variability in the SK1 δ18O and δ13C 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 δ18O and δ13C values of the SK1 composite records occur in MIS 5c and MIS 5a. We interpret the close positive relationship between δ18O 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 δ18O) also have some of the lowest δ13C values, probably supporting dense C3 vegetation and reducing PCP and fast CO2 degassing from fast dripping water, whereas during phases with less summer rainfall, slower dripping in the cave enhanced CO2 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., Reference Kirst, Schneider, Müller, von Storch and Wefer1999; Jouzel et al., Reference Jouzel, Masson-Delmotte, Cattani, Dreyfus, Falourd, Hoffmann and Minster2007), atmospheric CO2 concentrations and sea levels decrease (Bereiter et al., Reference Bereiter, Lüthi, Siegrist, Schüpbach, Stocker and Fischer2012; Rohling et al., Reference Rohling, Grant, Bolshaw, Roberts, Siddall, Hemleben and Kucera2009, Reference Rohling, Braun, Grant, Kucera, Roberts, Siddall and Trommer2010), and rainfall in the South African summer rainfall region decreases and shifts away from a precession pacing (Partridge, Reference Partridge1997). This change has also been associated with a shift from low- to high-latitude climate drivers in southern South Africa (Chase et al., Reference Chase, Harris, Wit, Kramers, Doel and Stankiewicz2021). The relationship of SK1 δ18O with temperature records after 70 ka shifts to positive, and the relationship with precession/insolation changes to a negative one (Fig. 7).
The composite δ13C record shows an overall increase from early MIS 4 through early MIS 3, which overlaps with a decrease of atmospheric CO2 concentrations, an increase of regional temperature, and a decrease of summer rainfall (Fig. 11a–c), possibly reflecting the impact of changing CO2 on vegetation (overall increase of δ13C in C3 plants, shift to more C4 plants) and decreased effective moisture on PCP and CO2 degassing in the cave. A peak of δ18O 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 δ13C values at the MIS 4–3 transition also mean that PCP is a possible influence on δ18O. 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., Reference Maccali, Meckler, Lauritzen, Brekken, Rokkan, Fernandez, Krüger, Adigun, Affolter and Leuenberger2023).
MIS 2
The stable isotope records from MIS 2 at SK1 are dominated by a steep increase of δ13C and, to a lesser degree, δ18O. 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 δ18O and δ13C 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., Reference Engelbrecht, Marean, Cowling, Engelbrecht, Neumann, Scott and Nkoana2019). Dry conditions at SK1 could also be related to low sea levels and continentality on the exposed PAP (Göktürk et al., Reference Göktürk, Sobolowski, Simon, Zhang and Jansen2023).
Conclusions
We present new records of speleothem composite δ18O and δ13C 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 δ18O, but additional variation is added by changes in rainfall seasonality and rainfall δ18O. Correlations of the δ13C 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 δ13C, whereas coastal records show variable correlation. Correlations between δ13C 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 δ18O and δ13C variations through time. Between the beginning of the record at 105 and ~70 ka, variation of δ18O 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 δ13C 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, CO2 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 δ18O and other proxies change, and now higher δ18O is associated with higher temperatures and less rainfall. Semblance values for the comparisons of δ13C to other proxies do not see a similar change as those for δ18O. 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 δ18O. Higher δ18O and δ13C 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.
SUPPLEMENTARY MATERIAL
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2024.3
Acknowledgments
We thank Nick Scroxton and two anonymous reviewers for constructive criticism that improved this article. We thank Christo Cloete and his family for granting us access to the cave and giving permission for sampling and for their support with our sampling efforts. We thank Nicol van Niekerk for assisting us with accessing the cave and Jacob Harris, Christopher Shelton, and Elmore Becker for assistance with speleothem sampling. We also thank the Institute of Human Origins (IHO) and School of Human Evolution and Social Change staff at Arizona State University (ASU) and the South African Coast Paleoclimate, Paleoenvironment, Paleoecology and Paleoanthropology Project (SACP4) and Mossel Bay Archaeology Project: Cultural Resources Management PTY. (MAPCRM) staff for their assistance; the Dias Museum for field facilities; and the Geological Survey of Israel for research support. We acknowledge funding from the European Commission Seventh Framework Marie Curie People program FP7/2007–2013 through funding of the Initial Training Network “GATEWAYS” (http://www.gateways itn.eu) under grant no. 238512. Additional funding was received from the U.S. National Science Foundation (grants BCS-0524087, BCS-1138073, and BCS-1460376 to CWM; grant 2002486 to KB; and grant 2002474 to RLE), the Hyde Family Foundations, the Institute of Human Origins (IHO) at Arizona State University, and the John Templeton Foundation to IHO at ASU (grant ID 48952). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of any of these funding organizations.
Conflict of interest statement
The authors declare no conflicts of interest.
Open access
