Global simulations at 70 and 82 ka were conducted using the fast version of the Norwegian Earth System Model (NorESM1-F). A horizontal resolution of ∼2° and 26 vertical levels are present in the atmosphere component. The ocean-sea ice component is on a tripolar grid with a nominal 1° resolution and 53 vertical layers in the ocean. Greenhouse gas levels (Martínez-Botí et al., 2015), orbital parameters (Berger and Loutre, 1991), and global mean sea-level (GMSL) drops (Rohling et al., 2014) at 70 and 82 ka can be found in Supplementary Table S1. The changes in the land-sea distribution arising from GMSL variations are also considered in the model setup. For instance at 70 ka, as a result of the 75-m drop of the GMSL with respect to today, the Bering Strait becomes closed and much of the Sahul, Sunda and East Siberian continental shelves are exposed. The simulations at each epoch were initialised from the Levitus temperature and salinity (Levitus and Boyer, 1994). The model provides output with a temporal resolution of 6 h for the atmosphere, 24 h for land surface, and 1 month for the ocean. The reader is referred to Guo et al. (2019) and Göktürk et al. (2023) for a more detailed description of the global model, including its technical specifications.

We used the Weather Research and Forecasting Model (WRF) as the RCM. WRF is a limited area atmospheric model that provides geographically detailed information, and is widely used in climate research (Powers et al., 2017). A single domain that covers the entire African continent south of the equator (0°–50°S, 0°–55°E; Figure 1B) was used to run WRF version 4.5.1 (Skamarock et al., 2019). The horizontal resolution is 12 km, while there are 35 vertical levels between the surface and ∼50 hPa level. The data generated by the global NorESM1-F simulations (described in the previous subsection) were downscaled using WRF, in order to determine climatic changes from 82 ka to 70 ka. There are two 30-year simulations at each epoch, with present-day and reconstructed coastlines (see Table 1 and the next subsection for details). This enables us to disentangle the effects of large scale changes from those of local coastline shifts. Two continuous WRF runs of 10 and 20 years were performed to form one complete downscaling simulation of 30 years. The global model output to be downscaled by WRF was chosen such that the 30-year means of temperature and precipitation at the 20 NorESM grid points around coastal southern Africa have minimal deviation from their long-term averages in the entire NorESM output available. The aim here was to force WRF with boundary conditions that are representative of the NorESM time slice, and to minimize the influence of internal variability. In other words, our model setup helps avoid results being impacted by anomalously warm, cold, wet or dry years, or decades. As model spin-up time an extra year was added to the beginning of each WRF simulation period, to be discarded from the analysis later. Solar radiation calculations within the WRF source code were adjusted using orbital parameters at 82 and 70 ka, as done in Yoo et al. (2016). To obtain high resolution sea surface temperature (SST) input, the anomaly method (Haywood et al., 2010; Zhang and Yan, 2012) was used. This approach makes use of the anomalies between the present-day SSTs and the global model’s pre-industrial SSTs. The method is superior to using the raw SST output from the global model directly, as it ensures an improved representation of the SSTs in the regional model. The resulting SST input fields for the WRF simulations are on the ERA5 reanalysis (Hersbach et al., 2020) grid.

To simulate the climatic effects of coastal land extent, the WRF land-sea mask was modified for the exposed land (ExpLand) simulations at each epoch (Figure 1A). This was done by assuming sea levels of −42 and −80 m at 82 and 70 ka, respectively, with respect to present-day sea level, which is used for the control (CTRL) simulations. These values, with ranges given in Table 1, are the maximum likelihood estimates in the global mean sea level reconstruction by Grant et al. (2012). This record matches well with the sea level index points determined by Cawthra et al. (2018; 2020) for the southern continental shelf of South Africa. It also has higher temporal resolution and lower age uncertainties compared to other available sea level reconstructions (Spratt and Lisiecki, 2016). On the other hand, we acknowledge that paleo sea level estimates are associated with uncertainties and may not account for local variations in specific regions. For instance, for MIS 5a (∼82 ka), Pico et al. (2017) and Dorale et al. (2010) report sea levels that are around present-day values, from the United States mid-Atlantic region and Mallorca, respectively. Yet, we believe the values we use in our regional model are realistic, as sea level for southern Africa tracks global mean sea level estimates, as revealed from the works by Cawthra et al. (2018, 2020) in the area.

Analyses of the general sedimentary matrix of the MSA deposits in South African cave sites show that they contain a significant amount of organic material derived from terrestrial plants and trees (Goldberg et al., 2009; Collins et al., 2017; Haaland et al., 2020; Wadley et al., 2020). This organic material contains well-preserved lipid compounds derived from leaf epicuticular waxes, the wax coating covering the outer surface of the leaf cuticle in land plants (Eglinton and Hamilton, 1967). Despite the fact that plant waxes are omnipresent and abundant in many soils today, they are not preserved in all depositional settings and are susceptible to biodegradation. Hence it is not trivial to successfully establish a plant wax biomarker sequence from cave deposits. In 2019, 12 sediment samples were collected from the south wall vertical section at BBC for long-chain (C25-C33) n-alkane analyses (Table 2; Figure 2).

Organic geochemical and stable isotope analyses were performed at the organic geochemistry laboratory of Lamont Doherty Earth Observatory (LDEO) and the FARLAB stable isotope facility at University of Bergen (UiB). For this study, we analysed 12 samples that were extracted at semi-regular intervals from the section of the BBC sediment stratigraphy (Table 2; Figure 2).

Each sample was freeze-dried and homogenised using a mortar and pestle. Lipids were then extracted from each 13-41 gr sample with a 9:1 (v/v) solution of dichloromethane (DCM) and methanol using a Dionex 350 Accelerated Solvent Extractor (ASE). The resulting Total Lipid Extract (TLE) was dried under a stream of nitrogen gas in a Biotage TurboVap.

Compounds were next separated using activated silica gel flash column chromatography. Samples were sequentially eluted with four bed volumes of hexane, DCM, and methanol to isolate the aliphatic, ketone, and polar fractions respectively. n-alkanes (in the aliphatic fraction) were identified and quantified using an Agilent 7,890 gas chromatograph equipped with a mass selective detector (MSD), flame ionisation detector (FID). Gas chromatography was performed with a 30m DB-5 capillary column. The GC oven was held at 60 °C for 1.5 min before the temperature ramp: first to 150 °C at a rate of 15 °C per minute, and then to 320 °C at 4.5 °C per second. We identified and quantified n-alkanes by comparing sample chromatograms to a commercially available standard containing C10-C40 homologues. Hydrogen isotope measurements were performed on a Thermo Trace GC coupled to a Thermo Delta V Isotope Ratio Mass Spectrometer (IRMS) via an Isolink Conflo IV interface. To convert compound-specific δ²H values to the VSMOW scale, we measured the authenticated A7 n-alkane mixture provided by Arndt Schimmelmann from the University of Indiana. The A7 was measured after 6 samples, at least. Each run was monitored for instrument drift with the A7 before calibration. Sample replicates were measured in separate runs and thus had separate calibrations for each run. Each sample was measured at least twice, and up to five times. The precision of replicate measurements was always better than 2.6‰ (1σ). An n-alkane mixture (even n-alkanes only) was measured as an unknown control throughout the 5 runs, and δ²H precision was measured to 1.95‰ and 1.83‰ (1σ, N=10) for C30 and C32 respectively.

To assess the level of degradation of the plant waxes, we calculated the carbon preference index, CPI (Bray and Evans, 1961, (Supplementary Figure S2)). The CPI reflects the molecular distribution of odd-to-even n-alkanes within a specific carbon number range (here n-C21 to n-C33), (Herrera-Herrera et al., 2020). High CPI values (>5) indicate a significant contribution of waxes from terrestrial plants and low levels of degradation while low CPI values indicate low input of plant waxes from terrestrial plants or high degradation (Eglinton and Hamilton, 1967). The average chain length (ACL_21-35) was calculated following Poynter et al. (1989) to further help assess changes in vegetation (Supplementary Figure S2). Ratios between specific n-alkanes can be calculated as they are sensitive to environmental influences (e.g., Schefuß et al., 2003; Carr et al., 2014). The ratio of the n-alkane homologues C29 and C31 was used to calculate the Norm31 index (Supplementary Figure S2); (e.g., Carr et al., 2014). Using this index together with the general n-alkane distribution patterns can reveal the differences between plant types, and their affinity to distinct biomes. In South Africa, soils with the highest ACL and Norm31 values have been found in the arid-adapted Succulent Karoo biomes, followed by Nama Karoo. Soils with lower Norm31 index were attributed to more humid-adapted biomes such as Grassland, Savanna and Fynbos (Carr et al., 2014; Herrmann et al., 2016).

Several significant large-scale changes are observed from 82 to 70 ka. For instance, seasonality over land is more pronounced at 70 ka, with up to 4 °C higher (lower) mean surface air temperatures during austral summer (winter) compared to 82 ka (Figure 3). Widespread increases in rainfall during summer are also evident (Figure 4). Summer warming over land areas occurs despite the cooling effects of increased cloudiness and evaporation (not shown), as well as from the higher overall elevation (caused by the lower sea level) at 70 ka (Göktürk et al., 2023). These epochal differences in temperature and precipitation are mainly associated with an increase (decrease) in summertime (wintertime) insolation from 82 to 70 ka, owing to a half-period advance of the orbital precession cycle (Thompson et al., 2005). The accompanying free troposphere (500 hPa) changes feature anomalous northerly and easterly flows (with respect to 82 ka) south of 20°S (Figure 4), which indicate enhanced monsoon circulation and increased moisture transport into southern Africa from the tropics and the Indian Ocean during summer at 70 ka. The southern boundary of the wetter area during summer (Figure 4) is located roughly where the anomalous northeasterly flow is replaced by a northwesterly flow from the cooler Atlantic Ocean, highlighting the role of the large-scale circulation influence on surface climate. There are precipitation increases at 70 ka over inland southern Africa in winter too, though they occur largely in the winter-dry region, thus are generally small in absolute values.

Over the oceans, difference patterns in mean temperature and precipitation are broadly compatible with SST changes. SSTs are generally lower at 70 ka around the continent, for instance with an Agulhas Current that is locally more than 2 °C cooler than at 82 ka off the coast of southwestern Cape region in winter (Supplementary Figure S1). Weaker westerlies (Figure 4), as well as the aforementioned drop in SSTs at 70 ka lead to reductions in rainfall of up to 40% over the ocean off the southern Cape during winter (Figure 4), compared to 82 ka. In similar fashion, rainfall reductions (up to 40%) are simulated over the Indian Ocean off the south eastern Cape coast at 70 ka in summer, when the Agulhas Current is 0.5°C–1 °C cooler than 82 ka. This is in stark contrast with the increasing summer temperatures and rainfall inland, where higher insolation at 70 ka controls the pattern of change (Figure 4).

Overall, large-scale differences between 70 and 82 ka are dominated by insolation-driven changes over land and SST-driven changes over the ocean.

The BBC samples contained plant wax-derived n-alkanes with chain lengths ranging from C21 to C35 carbon atoms (Figure 12). Their concentrations and elevated CPI values (4.5–16.5, with an average of 9.8) in all samples document moderate input of non-degraded plant-waxes (Supplementary Figure S2). The C31 n-alkane was consistently the most abundant compound in all samples with the C33 n-alkane being the second most abundant homologue. Leaf wax n-alkane concentrations are highly variable; they exhibit a range from 0.1 to 1.87 μg g-1 dw (for the C31 n-alkane) and as such are similar to present day Lowland Fynbos soils (0.4-5.6 μg g-1 dw), (Herrmann et al., 2017). Three sampled layers: CP (M3 unit), CGAC and CGAA/CGAB (M2 lower unit) yielded insufficient n-alkanes for analysis (Table 2).

Leaf-wax ACL_21-35 exhibit stable values with an average of 30.2 (±0.2), (Supplementary Figure S2). Norm31 values are variable and range between 0.64 and 0.76 with an average of 0.68 (±0.4) (Supplementary Figure S2). Throughout the sequence the ACL and CPI values display no overall trend (Supplementary Figure S2). The ACL-CPI cross plot reveals that values of the BBC samples differ significantly from samples from burning experiments of Wiesenberg et al. (2009) as they fall within the range of unburnt straw and soils. As such, we conclude that heating/burning has left no impact on our samples or that the plants were not heated or burned in those layers (Supplementary Figure S2).

Hydrogen isotope composition (δ²H_wax) values of the n-C31 n-alkane range throughout the MSA sequence in BBC from −114‰ ± 7.4 to −157‰ ±0.9 (Figure 12). There is significant inter-layer δ²H_wax variability within depositional units (Figure 12). Samples from unit M2: CGAC and CGAA/CGAB yielded n-alkane abundances too low for reliable δ²H determination. All reported plant-wax data are available in our Supplementary Material.

The distribution of long-chain (C27-C33) compounds indicates that n-alkane homologs predominantly derive from terrestrial vegetation with all samples closely resembling the average of modern Lowland Fynbos vegetation (an evergreen, sclerophyllous shrub) (Figure 12); (Carr et al., 2014). Throughout the sampled sequence (91.8± 6.5–68.3 ± 5.4 ka), no overall trend in Norm31, ACL and CPI indices and n-alkane chain length distribution occurs (Supplementary Figure S2). This might imply that no significant vegetation changes occurred over time in the area surrounding BBC. Given the known climate changes driven by orbital and SST forcing and the strong local changes arising from sea level drop suggested by the regional model, the lack of variation in vegetation type might indicate that the collecting habits of the inhabitants were biased towards collecting the plants resembling the modern Fynbos biome despite changing climate conditions and vegetation around the cave. However, this assumption requires further investigation.

The samples analysed from the M3 phase (CP, CJ, CIBh2 and CH) consisted of the lowest and oldest parts in the MSA section of BBC and are dated to fall within MIS 5c/b. The δ²H_wax values in M3 range from −137‰ ± 1.4 for layer CJ; the most depleted values in the record are −157‰ ± 0.9 for layer CIBh2 and values of −140‰ ± 1.4 in layer CH, the youngest in the M3. The averaged value (orange star) for the M3 phase (n=3) is −144‰ ± 11, (Figure 2; Figure 11). The two samples from the Upper M2 phase (dated to MIS 5a), layers CFB/CFC and CFA, had stable values of −130‰ ± 2.6 and −133‰ ± 0.6, respectively. The phase’s average is −132‰ ± 2.3. The δ²H_wax values in the M1 phase range from −115‰ ± 7.4 in the CD layers, the lowest (oldest in that phase) to −142‰ ± 0.9 in the CC layer and −138‰ ± 1.2 in the layer CB. The latter sample is positioned on top of the phase stratigraphically sitting at the MIS 5a/4 transition. The phase’s mean average is −131‰ ± 14.9. Positioned at the boundary between CA and the dune layer (MIS 4) the δ²H_wax value for the youngest sample analysed in the sequence is −119‰ ± 1.9.

The range of δ²H_wax values throughout the sequence is in line with the contemporary soil samples from the YRZ (Herrmann et al., 2017; Strobel et al., 2020) suggesting that sedimentary δ²H is representative of the mean δ²H of vegetation from the region surrounding the cave.

Increasing δ²H_wax values throughout the sequence suggest a shift in hydroclimate at BBC from MIS 5c/b to glacial MIS 4 with a coeval drop in sea level (Figure 12). The interpretation is not dependent on ice volume correction of seawater δ²H, which means, applying a correction for global ice volume changes on the proxy record does not change the observed pattern. Correcting for global ice volume (∼8‰ δ²H enrichment of seawater from increased global ice volume during glacial periods (Duplessy et al., 2002)) is common practice in paleowater isotope studies. However, it assumes that the same net change is experienced everywhere on the planet. This view was recently challenged (Windler et al., 2020).

One way to interpret these changes is as a decrease in precipitation amount over time (Sachse et al., 2012) in the landscape surrounding BBC. Reduced overall rainfall amounts would lead to an enrichment process by which the relative abundance of the heavier isotope increases, making δ²H of rainfall more positive. When it is drier, δ²H is more enriched as raindrops evaporate more in a drier atmosphere (Dansgaard, 1964). Evapotranspirative enrichment of leaf water is an additional driver of δ²H_wax. Changes in regional climate which alter the balance between precipitation and evapotranspiration can also influence the δ²H_wax signal (Schefuß et al., 2005), with higher values suggesting increased evapotranspiration associated with long, dry summers, i.e., higher surface air temperatures in southern South Africa and/or lower relative humidity.

Alternatively, the δ²H_wax signal could be interpreted as changing rainfall seasonality contributing precipitation from different moisture sources (Figure 12). Southern Africa is characterized by large differences in the isotopic composition of precipitation. Rainfall brought by the westerlies in winter is more depleted in δ²H, while it is more enriched when brought in via easterlies in summer (Harris et al., 2010; Braun et al., 2017). It should be noted that the modern-day intra-zonal values of the southern-YRZ (Figure 12) are large and overlap with the range of δ²H_wax values for the SRZ. This makes it challenging to determine whether the δ²H_wax data indicate a trend towards proportionally more summer precipitation or more uniform rainfall throughout the year, as time progresses.

We view the possibility that the shift in values may reflect input of more shrub-like vegetation rather than grasses as unlikely. This is because the input of more shrub-like vegetation rather than grasses throughout the sequence, is not consistent with the CPI_25-33, ACL_21-35 and Norm31 indices reconstructions; that is, throughout the sampled sequence (91.8± 6.5–68.3± 5.4 ka) no overall trend in these indices occurred (Supplementary Figure S2). Supporting this is the micro-mammal reconstruction from BBC which indicates that shrubland oscillated, with the greatest abundance associated with the M3 phase and an overall decline towards MIS 5b/a in the sequence (Nel and Henshilwood, 2021).

Due to the lowering of the sea-level over time, the relative position of the BBC with regards to ocean moisture sources changed. Elevation difference can influence rainfall δ²H signatures which become deuterium-depleted with altitude (ca. 10%–15% per 1,000 m; Gonfiantini et al., 2001). However, the change in altitude is too small for this effect to be significant (<100 m, i.e. 1-1.5‰ deuterium depletion), and the trend in the δ²H_wax signal over time is the opposite of what would be expected due to the altitude effect on deuterium.