This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.
Introduction
Speleothems are important terrestrial archives used for reconstructing paleoenvironmental conditions. The reliability of speleothem-derived data is a direct result of the close and well-known relationship between physico-chemical parameters (e.g. the isotopic compositions of C and O) and the conditions of speleothem crystallisation (e.g. McDermott, 2004; Lachniet, 2009; Fairchild and Baker, 2012). Additionally, it is possible to reliably date speleothems via U-series (e.g. Scholz and Hoffmann, 2008; Baker et al., 2011; McDermott et al., 2011) or paleomagnetic/magnetostratigraphic (e.g., Morinaga et al., 1994; Lascu and Feinberg, 2011) dating methods. However, neither method is perfect. Uranium-series dating, like all physicochemical methods, is influenced by several limitations; in this case, the U concentration must be at a minimum level of few parts per million (Hellstrom, 2006; Fairchild and Baker, 2012), the admixture of detrital impurities must be as low as possible, and the maximum age range of the method is ~700 ka (Cheng et al., 2016).
Palaeomagnetism can provide only correlative ages and relies on the reconstruction of polarity changes in Earth's magnetic field; if these records are not calibrated with other methods of stratigraphic research, matching with the Geomagnetic Polarity Time Scale (GPTS; Cande and Kent, 1995) is difficult. The erosion or discontinuous deposition of cave sediments renders such correlations even more problematic as hiatal periods can hide substantial amounts of time (Bosák et al., 2000, 2002, 2003; Zupan Hajna et al., 2008, 2010, 2020). Therefore, the polarity changes located in hiatuses may not record the age as expected from the GPTS (Lascu and Feinberg, 2011). So oxygen isotope stratigraphy (OIS) was applied here to acquire a more detailed chronology of the speleothems examined. The analysis of the stable isotopic composition of oxygen is a standard method for tracking paleoclimatic changes in speleothems. The resultant oxygen curves can be used for correlation with a reference record that exhibits global or regional characteristics (Lisiecki and Raymo, 2005; Wang et al., 2010).
Note that many speleothem studies are based on relatively young samples (e.g., Lauritzen and Lundberg, 1999; Vaks et al., 2006; Affek et al., 2014) and samples older than Pleistocene-age are rare (Vaks et al., 2013). Cave sediments in the classical karstic Račiška Pečina Cave in SW Slovenia are a part of this rare category, having a total age range from 3.2 Ma to ∼2 ka (Zupan Hajna et al., 2021). They also contain well-developed and time-limited normal-polarised Olduvai Subchron ages (1.95–1.77 Ma) (Zupan Hajna et al., 2008, 2010, 2020). The Pleistocene and Pliocene ages of the Račiška Pečina profile are interesting because of the climatic and paleogeographic changes that occurred on Earth's surface during these periods (Zacwijn, 1974; Van Couvering, 1996; Lourens et al., 2004). The climate dried, steppisation occurred, and glacial development began, initially in the southern hemisphere and later in the north (Van Couvering, 1996; Gradstein et al., 2012). Characteristic cycles of warm (interglacial) and cold (glacial) periods also began in the Pleistocene (Stanley, 2005).
Carbon and oxygen stable isotopes (expressed as δ13C and δ18O values,
{\rm{\delta }} = \left( {{{R\,\,sample} \over {R\,\,{standard}}} - 1} \right)
, where R is the 13C/12C and 18O/16O ratio in the sample and the standard) represent the most frequently applied proxies for the reconstruction of paleoenvironmental conditions from speleothems (McDermott, 2004; Lachniet, 2009; Fairchild and Baker, 2012; Ünal-İmer et al., 2016; Demeny et al., 2017). The C and O isotope compositions in speleothem calcite are primarily controlled by the isotopic composition of the percolating water and the δ18O of dripping water depends on the source of moisture, its transportation, the intensity of atmospheric precipitation and ambient air temperature. The oxygen isotope compositions of precipitation are altered by local conditions, such as temperature and the amount of precipitation at a given site (Lachniet, 2006; Lambert and Aharon, 2011; Fairchild and Baker, 2012). During the slow degassing of CO2 and without evaporation, speleothem calcium carbonate crystallises near or in isotopic equilibrium with dripping water (Lauritzen and Lundberg, 1999) and as a function of the air temperature of the cave (Hendy, 1971).
Many studies have shown (e.g., Ayalon et al., 2002; Levin et al., 2009; Demeny et al., 2017) that the δ18O values of cave drip waters are constant over a long period of time and can nearly equal the oxygen isotope composition of the mean annual precipitation outside of the cave (e.g., Sundqvist, 2007). The temperature in a cave is also very close to the mean annual (surface) temperature in the region above the cave (Sundqvist, 2007; Fairchild and Baker, 2012). Carbon isotope compositions can provide information on climatic and vegetation conditions (McDermott, 2004; McDermott et al., 2006). The main factor controlling the δ13C (13C/12C ratio) of speleothem calcites are the proportions of CO2 from organic sources and dissolved carbonate rocks (McDermott, 2004; Fohlmeister et al., 2020). The C isotope composition of soil CO2 is affected by the type of vegetation existing above a cave, which depends on photosynthetic pathways (C3 or C4 plants) and soil respiration rates (McDermott, 2004; Spötl et al., 2016; Fohlmeister et al., 2020). For C isotopes, the fractionation effect is temperature-dependent and is initiated by the dissolution of soil CO2 in infiltrating meteoric water (Hendy, 1971). The main objectives of this study were to improve estimated speleothem chronologies based on OIS, which was tested as an alternative dating method, and to develop a regional paleoenvironmental interpretation of the obtained records (including the δ13C and δ18O records).
Geological Setting
The Račiška pečina Cave (RP) (45°30′N; 14°09′E; 609 m a.s.l.) (Fig. 1A) is located in the south-eastern part of the Matarsko Podolje, a part of the classical karstic terrain in SW Slovenia (Zupan Hajna et al., 2008, 2010). The RP represents a paragenetic or epiphreatic relict cave, which is 304 m long, 10 m wide at maximum (Pruner et al., 2003) and a simple south-dipping gallery (Bosák et al., 2004). It is part of an old cave system, which was partly opened to the surface by denudation (Zupan Hajna et al., 2008). The cave formed in thickly bedded limestones, dolomites and breccias of Lower Cretaceous period (Šikič et al., 1972). Allogenic siliciclastic sediments, primarily clays and sands of unknown thicknesses that fill the lower part of the cave, are covered by various layered speleothems (Pruner et al., 2003). A 13-m-long profile through clearly stratified speleothems intercalated by muds is situated in the southern part of the cave, approximately 200 m from the present entrance (Fig. 1B). The primary thickness of the studied section reaches 223 cm (Horáček et al., 2007; Zupan Hajna et al., 2008, 2010) (Fig. 1C) and the profile is attached to the cave wall at the passage junction with a slightly corroded side fissure. The section is horizontally divided into two parts — a longer left segment and shorter right one. The shorter segment shows a different evolutionary history in the upper part; it is partly detached from the longer segment along the crack and has slightly subsided. The boundaries of both segments are expressed by columnar stalagmites at different vertical positions (to the right of the position of the RP profile in Fig. 1C) and three principal parts may be distinguished in a vertical profile (Pruner et al., 2003; Zupan Hajna et al., 2008) (Fig. 1C). The magnetostratigraphic log obtained was sufficiently correlated with the GPTS of Cande and Kent (1995) biostratigraphically (Horáček et al., 2007; Moldovan et al., 2011). The past application of multiple different dating methods (see also Zupan Hajna et al., 2008, 2010, 2020, submitted) means that the RP profile contains the best-dated cave sediments in Slovenia. The cave evolution was related to local contact karst characteristics (Mihevc, 2001, 2007); the location of the recent ponors of allogenic streams flowing from the Eocene Flysch towards springs at the Adriatic Sea was confirmed by tracer test (Krivic et al., 1989). The cave was located deep in karst as a part of a larger cave system and the cave functioned at 3.2 Ma ago (Zupan Hajna et al., 2008) as a conduit for underground water flows from the impermeable flysch catchment area to the springs. Before the speleothem dome began to grow, the cave was disconnected from its hydrological function due to the regional tectonic uplift (Mihevc, 2007), but remained there without the direct ventilation. The cave was later opened to the surface by surface chemical denudation and slope retreat; during this period, cave bears started to enter the cave. The recent phreatic conduits are positioned about 200 m below Račiška pečina and due to an average karst denudation rate of 20–60 m/Ma (Mihevc, 2001) in the area, the cave was opened to the surface in the last 0.5 Ma.
Fig 1
(A) Location of the study site in Slovenia; (B) cave map and profile location (marked with a red circle) (modified from Zupan Hajna et al., 2008); (C) studied profile scheme: grey shading — flowstones; brown shading — clays; diagonal shading (light blue) — collapsed limestone blocks; yellow crosshatching — younger speleothem cover; red line — principal disconformity; F1 — faunal position; USp — position of Ursus speleus sp (cave bear) bone fragments; RP — sampling profile for C and O isotope analyses; and (D) magnetostratigraphic log according Pruner et al., (2009) correlated with the GPTS (Cande and Kent, 1995); black — normal polarity, white — reverse polarity, grey — transitional polarity, wavy line — disconformity. GPTS, Geomagnetic polarity time scale.
The RP profile can be divided into three lithostrati-graphic units namely lower, middle and upper. The lower unit represents three growth stages of a massive vaulted speleothem dome built by brown to reddish-brown massive, highly porous speleothem layers (up to 25 cm thick) with interbedded red muds (1–2 cm thick). The internal speleothem texture is laminated (sub-millimetre laminae to millimetre-sized bands) and mostly composed of columnar and radial-columnar coarse-grained calcite crystals. There are also several very thin (centimetres in thickness), light-coloured bands of fine-grained (micritic) speleothem beds with a porcelaneous macroscopic appearance. The bottoms of the beds above the muddy intercalations often contain clasts of underlying muds. Fenestrate (vuggy) porosity occurs in places, and part of the porosity is surely intercrystalline among large calcite crystals, while some pores represent the effects of later corrosion. Two disconformities separate individual growth stadia. Broken stalagmites are preserved at the top of the upper (top) two sequences. With respect to the GPTS-matched magnetostratigraphic log, the highest possible age limits of the basal section are found to be from >3.2 Ma (bottom) to ∼2.58 Ma (top, i.e., the bottom of the middle profile; Horáček et al., 2007) (Fig. 1C, D).
The middle unit of the RP profile predominantly consists of sub-horizontal, laminated to banded, sometimes wavy bedded, partly porous flowstones that are reddish-brown to brownish red in colour, with rimstone dams filled with red and red-brown muds. The internal walls of the rimstones are covered by crystalline corallites. The fill of one rimstone (to the right of sample RP66, Fig. 1C) contain oribatid mites, which are forest taxa dated to the Early Pleistocene, specifically to the termination of the Olduvai Subchron (Moldovan et al., 2011). Individual thinner beige beds (centimetres in thickness) are fine-grained and there are several red mud intercalations (1-mm to up to 21-cm-thick). Red muds found to contain mainly quartz along with some muscovite/illite and chlorite, hematite and magnetite, or kaolinite and feldspar, based on X-ray diffractometry. Some muds also contain small iron-rich spherules. Thin calcified siltstones with flowstone fragments and iron-rich spherules occur in places. Near the base of this unit, collapsed blocks from the ceiling cover thicker layers of red muds with faunal remains (F1 in Fig. 1C) dated as Early to Middle Pleistocene MN17 (ca. 1.8–2.4 Ma; Horáček et al., 2007). Fauna-containing sediments are underlain by several flowstone layers that unconformably overlay the principal unconformity that terminates the lower section. The bottom of the middle unit is dated to ~2.58 Ma according to magnetostratigraphic and biostratigraphic correlations (Horáček et al., 2007) (Fig. 1C, D).
The upper unit of the RP is mostly represented by light-coloured (honey, ochre etc.), thickly bedded, and densely laminated flowstones. Sub-millimetre to millimetre-scale laminae are sub-horizontal, sometimes slightly wavy, or uneven. The internal texture consists of fine columnar to radial-columnar facies with distinct fine laminae, sometimes with distinct porosity (both primary intracrystal-line and partly corroded). Basal beds show greater crystal sizes and higher porosity, and a reddish colour begins to appear. Light-coloured, porcelaneous, thin beds also occur in places. The uppermost flowstone layers have laminae to thin bands enriched in organic carbon (soot), which can be attributed to repeated early human settlements in the cave. Two lenses of brown-yellowish loams interbedded in the right part of the profile contain bone fragments of Ursus speleus sp. (cave bear) and clearly indicate a Pleistocene age (Fig. 1C). The left segment contains only one such layer, most likely correlating with the upper layer in the right segment. Several old, massive, conical stalagmites, which grow from the top of the middle unit, represent the top of the section. The bottom of the upper part lies at the Brunhes-Matuyama boundary (0.78 Ma; Horáček et al., 2007) (Fig. 1C, D). Muds present in between the flowstone layers in all parts of the RP profile resulted from the successive flooding deposition of well-sieved allochthonous sediments from a suspension (Horáček et al., 2007) passing through the fissure system in the epikarst and vadose zones and emptying at the intersection of the narrow side passage with the wide main cave corridor.
Regional Climate
Slovenia is influenced by several climate types, reflecting its geographic position. The climate is a mixture of (1) continental climate, which influences the major part of the country; (2) alpine climate, which prevails in the mountains of NW Slovenia; (3) coastal Mediterranean climate, which influences the south-western part of the country. Four major types of air masses influence the weather in Slovenia and they are (1) maritime polar air masses, which originate in the North Atlantic and North Sea; (2) maritime tropical air masses, which predominantly originate in the Azores area; (3) continental tropical air masses originating in northern Africa and Anatolia; and (4) continental polar air masses that originate in Scandinavia, Finland, Russia and on the Pannonian (Carpathian) Plain (Mezga et al., 2014).
The RP has a transitional Mediterranean climate, with a mean annual temperature of ~10.4°C according to long-term observations by the Ilirska Bistrica Meteorological Station at 415 m a.s.l., 14 km NE from the cave (http://meteo.arso.gov.si/). The warmest month is July (mean temperature: ~16.5°C), and the coldest month is January (mean temperature: ~5.4°C). The mean annual precipitation is approximately 1.356 cm, with a maximum in October (~164 mm) and a minimum in February (~80 mm). Observations from the Kozina Weather Station, 20 km northwest of the RP, show that the present mean δ18O value of local meteoric water changes from −13.8 ‰ in January to −2.2 ‰ in March (Vreča et al., 2005; Mezga et al., 2014).
METHODS
Sampling
Sampling for paleomagnetic and stable isotope studies was conducted with a very dense spacing (Bosák et al., 2002; Zupan Hajna et al., 2008). The samples of unconsolidated sediments (interbedded muds) used for palaeomagnetic analyses usually had centroid distances of 2–7 cm and were collected using non-magnetic plastic boxes (Natsuhara Giken Co., Ltd., Japan) with a size of 20 × 20 × 20 mm and approximate volume of 6.7 cm3. Speleothem samples were taken from the section by cutting a narrow trench with a powered circular saw (Bosák et al., 2002; Zupan Hajna et al., 2008; Pruner et al., 2009). The continuous slices or columns obtained were then cut into 20 × 20 × 20-mm cubes in the laboratory. All hand specimens were collected in situ. Dip directions were measured with a geological compass, and the north direction was drawn on the samples. All data were recorded in situ in field drawings and notes.
The samples used for stable isotopic analyses had two sources: (1) plates cut from the original palaeomagnetic samples (columns, slices) in the lower part of the section (i.e. in the domed speleothem below the principal unconformity) and at the top of the section; (2) the originally sampled RP profile above paleomagnetic sample No. RP 71 up to nearly the top of the section (RP line in Fig. 1C). The RP profile was situated at the centre of the sedimentary section which is in its longer left segment. Its position was selected with respect to (1) high-resolution stratification, including the preserved interlamination and interbedding of muds and flowstones pinching out toward both the right and left, and (2) previous magnetostratigraphic results used to cover the entire N-polarised magnetozone interpreted as the Olduvai Subchron (e.g., Horáček et al., 2007).
Flowstone laminae and beds were sampled relative to the orientation of the in situ hand specimens. Specimens were collected using a geological hammer and/or chisel, and sometimes taken only by hand. A total of 70 calcite speleothem samples with thicknesses ranging from 0.5–21 cm were obtained; some contained muddy intercalations. Solid samples were separated by laminae to beds of mud that ranged from sub-millimetre up to 15 cm in thickness.
Stable Carbon and Oxygen Isotope Analyses
Samples collected for δ13C and δ18O analyses were drilled using a drill (Dremel, USA) bit of diameter of 0.7 mm. All samples were taken along the RP profile with a mean resolution of 5 mm, producing 375 samples. Analyses of the C and O isotope compositions (δ13C and δ18O, respectively) in carbonates were performed using a KIEL IV Carbonate Device coupled to a Delta Plus isotope ratio mass spectrometer (IRMS) using a dual inlet system (Thermo Fisher Scientific, USA). Sample preparation was carried out automatically.
Approximately 100 mg of calcite was reacted with anhydrous orthophosphoric acid at 70°C for 5 min. The isotopic composition of the released CO2 was then measured using a Finnigan Delta Plus IRMS (Thermo Fisher Scientific, USA). The results were normalised to three international standards: NBS 19, NBS 18 and IAEA CO 8, and were reported relatively to the Vienna Pee Dee Belemnite (VPDB) international standard. The analytical precision (1σ) was yielded 0.03 ‰ and 0.08 ‰ for δ13C and δ18O, respectively. Reproducibility was assessed by measuring two internal standards after every 12 samples (for δ13C: ±0.03 ‰; for δ18O: ±0.08 ‰). Analyses were performed in the Stable Isotope Laboratory at the Institute of Geological Sciences, Polish Academy of Sciences in Warsaw.
U-Series Dating
A series of 0.5 to 1.0 g calcite samples were drilled from samples from the upper part of the RP flowstone (marked as RP66 and RP 1 in Fig. 1C). Only the parts with no macroscopically visible traces of porosity or secondary alteration were selected for dating. Following the thermal decomposition of organic matter, the samples were spiked with a 233U–236U–229Th mixture prior to further chemical treatment. The calcite samples were dissolved in nitric acid. Uranium and thorium were separated from the carbonate matrix via chromatography with TRU Resin (Hellstrom, 2006). Standard and blank samples were prepared simultaneously for every series. The mass abundances of 233U, 234U, 235U, 236U, 229Th, 230Th and 232Th were measured using a double-focusing sector field inductively couple plasma–mass spectrometer (ICP–MS) (Element 2, Thermo Finnigan MAT; Thermo Fisher Scientific, USA) at an ICP–MS laboratory in the Institute of Geology, Czech Academy of Sciences in Prague, Czech Republic. Measured results were corrected to account for background and chemical blanks and the final results were reported as isotopic ratios. Uranium-series ages were iteratively calculated from the 230Th/234U and 234U/238U activity ratios using the following decay constants (year−1): λ238 = (1.55125 ± 0.0017) · 10−10 (Jaffey et al., 1971), λ234 = (2.826 ± 0.0056) · 10−6 (Cheng et al., 2000), λ232 = (4.95 ± 0.035) · 10−11 (Holden, 1990), and λ230 = (9.1577 ± 0.028) · 10−6 (Cheng et al., 2000, 2013).
Reconstructing Speleothem Chronology
The studied RP profile contained several age benchmarks and (Fig. 1C) was not continuous (i.e. calcite layers are bedded by clastic layers). Therefore, it was not possible to estimate one continuous chronology based only on the set of available benchmarks. Many calcite segments do not contain any age benchmarks. The OIS method was used to attempt reconstructing a continuous chronology for the profile. This method is based on correlating the studied δ18O record with a reference δ18O curve (Imbrie et al., 1984). The reference curve should meet two basic requirements. First, its age scale must be well defined and second, the reaction of δ18O values to climatic changes must be known, and it should be possible to compare it with the studied δ18O record.
Here, two records were used as potential reference curves. The first one, the global curve (LR04), is a stack record based on 57 cores with recorded carbonate isotopic variations from foraminiferal tests (Lisiecki and Raymo, 2005). The second curve is the regional Mediterranean (Medi) curve, which is based on stable isotopic records from planktonic foraminifera in the MD84641 core (Wang et al., 2010). Oxygen isotope stratigraphic correlations were performed using Gen-Corr v. 1 software (Pawlak and Hercman, 2016). Three types of benchmarks were implemented in Gen-Corr: (1) ‘fixed’ benchmarks, such as the known age at a specific point in the profile; (2) ‘normal’ benchmarks, such as the age information at a specific point described by a probability distribution, and (3) ‘border’ benchmarks, including age information described by age limits (e.g. younger than 3.2 Ma, etc.). The ages assigned for fixed benchmarks remained constant throughout the correlation procedure, while the ages of non-fixed benchmarks could change according to their probability distributions.
Age estimates for the points between benchmarks were required to meet two conditions: (1) the results must agree with the point stratigraphic position and (2) the ages must be in accordance with the ages of the benchmarks (Pawlak and Hercman, 2016). Here, three paleomagnetic benchmarks were used: (i) the Brunhes–Matuyama boundary, with an age of 780 ka at a depth of 280 mm; (ii) the penultimate benchmark, based on the age of the principal disconformity, which was magnetostratigraphically dated to 2.58 Ma; and (iii) magnetostratigraphic results, assuming that the entire outcropping profile could not be older than ~3.2 Ma (Horáček et al., 2007; Zupan Hajna et al., 2008; Pruner et al., 2009), and one reliable U-series age was used as the age benchmark (sample RP 1). Additionally, there is another sampled core in this profile, RP66 (Pruner et al., 2009; Hercman et al., submitted), covering Brunes-Mathuyama reversal and for which a high-resolution age-depth model exists based on paleomagnetic benchmarks and OIS correlations. The location of the RP66 sample in the RP profile can be easily defined via stratigraphic correlation (Fig. 1). The calculated borders of the RP66 profile (Hercman et al., submitted) were used as the two additional benchmarks for reconstructing the chronology of the studied RP profile (Table 1).
Ages obtained from borders of the RP66 flowstone (on Fig. 1) made by Hercman et al., submitted.
Marine carbonate and speleothem δ18O records respond in a similar way to climatic changes (see companion paper for detailed discussion, e.g. Almogi-Labin et al., 2009; Affek et al., 2014). Over long timescales, the global volume of glacial ice impacts the δ18O value of ocean surface waters (Dansgaard, 1964). In low-temperature conditions, the ocean surface water becomes depleted in 16O since evaporation process prefers this isotope. The development of ice sheets cover during glacial periods leads to tapping enriched in 16O water. Therefore, over long timescales, the global volume of glacial ice impacts the δ18O value of ocean surface waters (Dansgaard, 1964). The changes of speleothem carbonate δ18O value are more complex, they depend on δ18O of dripping water, which is driven by several factors such as source effect, the temperature effect at evaporation and atmospheric precipitation sites. Additionally, the local factors like evaporation at precipitation site or during the water percolation play an important role. All these factors can strengthen each other or eliminate their effects. At the coastal site, the source effect remains strong. Farther in the continental interior the source effect is more modified by continental effects and the temperature effect. Evolving from glacial to interglacial conditions the source effect would shift the speleothem δ18O value to the more negative direction, while the ambient air temperature effect would result in higher speleothem δ18O value. At Mediterranean coastal sites, the δ18O values of speleothem are determined by the temperature during the speleothem formation, the source effect and amount of precipitation water. In summary, during warm periods the deposited speleothem calcite can be enriched in δ18O. Therefore, marine data are presented on the inverted scale for correlation purposes. However, the amplitude of changes and mean values of δ18O for speleothems and marine records differ. Therefore, the normalisation of δ18O data must precede correlation. The normalisation eliminates the effect of amplitudes and shifts in absolute value. In this study, both reference and studied records were normalised by scaling to 0–1 range. After that, the records were correlated using modified Gen-Corr software (Pawlak and Hercman, 2016), which seeks to identify the highest similarly (optimal correlation) between two records using a genetic algorithm. Additionally, Gen-Corr allows the use of a priori age information during the correlation procedure, which allows the use of reliable benchmarks. The Gen-Corr seeks the optimal position between correlated records using similarity indicator. The used similar indicator is the mean distance between the correlated records.
Results
Stable Carbon and Oxygen Isotope Compositions
Analyses of 375 samples resulted in the construction of δ13C and δ18O records (Fig. 2). These records revealed that isotopic compositions changed significantly during the deposition of the RP profile. The δ13C values varied from −3.4‰ to −11.2‰, while the δ18O values varied between −4.2‰ and −7.1‰. The most conspicuous changes were observed near the main disconformity (red solid line in Fig. 2), before which the δ13C values had amplitudes of ~3 ‰ and after were ~7 ‰. For the δ18O values, before the main disconformity, the amplitudes were ~0.5 ‰ and after, they were ~2 ‰. Additionally, at a depth of ~900 mm, there was a visible change in the frequencies of the δ13C and δ18O records. The minimum C isotope values were found at depths of 340 −200 mm; there was also a decline in the O isotope values over the same depth range.
Fig 2
Stable oxygen (A) and carbon (B) isotope compositions in the RP profile. Vertical red lines — principal disconformity; vertical, dotted, and brown lines — main clay layers.
U-Series Results
The age of the lower part of the studied flowstone was beyond the U-series limit, unlike the upper part of the RP profile. Despite very high contamination by clays, iron oxides and hydroxides, and a relatively low U content (~0.1 ppm), it was possible to date only two samples collected in the upper part of the RP flowstone (Table 2; Fig. 1C). The RP 1 sample corresponded to the profile sample at a depth of 65 mm, for which the age was 120 ± 20 ka. The age uncertainty of RP66 is so large that it prevents from dating, so this sample cannot be dated by U-Th series dating.
Calculations use the decay constant of Jaffey et al., 1971 (238U), Cheng et al., 2013 (234U and 230Th) and Holden, 1990 (232Th). Ages do not include uncertainties associated with decay constants. Reported errors are equal to 2s.
AR, Activity ratio.
Age adjusted for detrital contamination, indicated by the presence of 232Th, using the typical silicate activity ratio 230Th/232Th = 0.83 ± 0.42, derived from the following activity ratios: 232Th/238U = 1.21 ± 0.6, 230Th/238U = 1.0 ± 0.1, and 234U/238U = 1.0 ± 0.1 (e.g., Cruz et al., 2005).
Chronology Estimates Based on OIS
Further, to establish the OIS chronology, the RP profile (Fig. 1C) was divided into two segments, separated by the principal disconformity (red line in Fig. 1C). Initially, both parts of the RP profile were correlated with the Medi and LR04 curves. The Gen-Core seeks the optimal position between correlated records using similarity indicator. The similar indicator used (as above) is the mean distance between the correlated records. Both sections showed different similarity indicators with specific reference curves (Fig. 3). The OIS results showed that the δ18O record from the lower segment (below the unconformity) of the profile matched the δ18O reference curve (Medi curve) from ~3.20 to 2.58 Ma, which is much better than that of LR04. Meanwhile, the upper segment (above the unconformity) fit the δ18O global reference curve (LR04) from ~2.48 Ma to ∼80 ka.
Fig 3
Comparison of similarity indicators in the Medi and LR04 records to the RP oxygen isotope composition (calculation in Gen-Corr by Pawlak and Hercman, 2016); (A) boxplots for the upper part of the RP profile and (B) boxplots for the lower part of the RP profile; the vertical axis represents the 2s uncertainty of the similarity indicator.
Age–depth models for RP profile segments were constructed independently from 80 correlations for each pair of curves (Fig. 4). During correlation, the number of depositional breaks was identified (grey sections in Fig. 4). The lower segment was divided into five parts, and the upper segment was divided into 15 parts.
Fig 4
Correlation results of the RP oxygen isotope record with the LR04 and Medi records. All data presented were normalised and the marine d18O data were inverted. (A) RP curve (blue line) for the upper segment of the profile correlated with the LR04 curve (black line); (B) RP curve (blue line) for the lower segment of the profile correlated with the Medi curve (black line); (C) age-depth model for the upper segment of the profile; and (D) age-depth model for the lower segment of the profile; red lines — confidence intervals; grey boxes — hiatuses; red boxes —principal disconformity.
The hiatuses designed by the correlation process agreed with the locations of clay intercalations observed macroscopically in both parts of the profile. The longest hiatus, at the time of the principal disconformity, was ~100 ka long (from ~2.58 to ~2.48 Ma); the longest hiatus in the lower part of the profile was ~60 ka long (from ~2.89 to ~2.95 Ma) and the shortest hiatus was ~30 ka long (from ~2.75 to ∼2.78 Ma). In the upper part of the RP profile, the longest break lasted ~100 ka (from ~650 to ~750 ka) and the shortest break was ∼25 ka long (from ~310 to ∼335 ka).
Discussion
Reliability of the Chronology
An OIS age–depth model was constructed for the RP profile based on 80 iterations. Due to many numbers of correlations, it was possible to establish that the age range of the profile was from ~3.2 Ma to 0.08 Ma, which includes marine isotope stages (MIS's) Km3 to 5 (Fig. 5). The uncertainty of the model was 2s bands. The similarity indicator for the lower segment was 0.082 and for the upper segment, it was 0.19 (Fig. 3). Therefore, the reliability of the lower part appears to be higher. However, the similarities of the segments to the reference curves were not uniform over time. In the case of the older profile segment, the segments were most similar in the following age ranges: from 3.1 Ma to 3.40 Ma and from 3.0 Ma to 2.95 Ma. Meanwhile, in the younger part, the sections in the following age ranges showed the greatest similarities to the reference curve (LR04): from 2.18 Ma to 1.50 Ma, from 1.25 Ma to 1.125 Ma and from 875 ka to 750 ka. The reliability of the correlation was higher in segments when the similarity indicator was lower (Pawlak and Herman, 2016).
Fig 5
The final set of stable carbon and oxygen isotope records of the RP profile; (A) MIS scale (warm — red; cold — blue) (Lisiecki and Raymo, 2005); (B) stable oxygen isotope record (smoothed data — blue; raw data — grey); and (C) stable carbon isotope record (smoothed data — black; raw data — grey); grey boxes — hiatuses; red box — principal disconformity; green boxes — episodes where the R2 value is significantly greater, R2 = 0.4538, R2 = 0.4944, and R2 = 0.6065. MIS, Marine isotope stage.
Factors Shaping the Reference Curves
The LR04 curve is based on 57 deep-sea cores, which record isotopic variations in the carbonate tests of foraminifera (Lisiecki and Raymo, 2005). The δ18O of foraminiferal tests is a function of the temperature and δ18O of the water in which they form, and the δ18O of seawater is a function of global ice volume and ocean salinity (Lisiecki and Raymo, 2005). The water in the deep ocean is more uniform in temperature and salinity than are surface waters. Therefore, the tests of benthic foraminifera reflect global climate change in a better manner (Lisiecki and Raymo, 2005). In contrast, the Medi curve is represented by the stable isotopic records of the planktonic foraminifer Globigerinoides ruber from core MD84641 (Wang et al., 2010). This species lives near the sea surface, and so they are more sensitive to changes in temperature and salinity. Therefore, they better reflect more frequent and smaller changes in the temperature, volume and salinity of seawater (Wang et al., 2010). Moreover, the Mediterranean Sea is well separated from the Atlantic Ocean, and therefore, its isotopic composition is more dependent on local factors than on global ones.
Speleothem δ18O values usually reflect the δ18O of moisture, precipitation, and the temperature difference between precipitation site and source region. The isotopic compositions of moisture and precipitation depend on the isotopic composition of their source (s) and the way they are transported (Darling, 2004; Lachniet, 2006). Therefore, speleothem records can be correlated with marine foraminifera records if the source effect is strong enough or if the speleothem record reflects the global changes of temperature.
Local Factors Influencing the RP Record
Speleothems from the Middle East record δ18O changes like the same patterns as those of the Mediterranean cores (i.e. source effect) (Bar-Matthews et al., 1999; Orland et al., 2012).
At the site of precipitation, isotopic compositions are modified by local conditions, such as the ambient air temperature and amount of rainfall (Różański et al., 1993). The thermal gradient for calcite crystallisation (δ18Oct/dt) varies from −0.18 ‰ to −0.23 ‰, while that for the δ18O of precipitated water (δ18Op/dt) ranges from +0.17 ‰ to +0.9 ‰ (Różański et al., 1993; Lachniet, 2006). Speleothem δ18O responses to temperature are strongly site-dependent (Williams et al., 1999; Mangini et al., 2005). Observations from the weather station in Kozina, Slovenia, 20 km northwest of the cave, show that the present δ18O value of local meteoric water ranges from −13.8 ‰ in January to −2.2 ‰ in March (Vreča et al., 2005; Mezga et al., 2014). However, March is not the warmest month in a region. Therefore, the fact that the March rainwater has the maximum δ18O value shows that not only the temperature but also other factors such as the source of moisture and its transportation that play important roles in shaping the δ18O value of rainwater in the studied region. Moreover, this seasonal difference shows how strongly the mean annual isotopic composition depends on the amount of rainfall in particular seasons.
In winter and autumn, the influence of Atlantic air masses prevails over the study area and is associated with an increased amount of rainfall during this time. In contrast, in spring and summer, influences from the Mediterranean Sea dominate when rainfall is lower throughout the year (Milošević et al., 2017). The source of rainfall also varies by season and is related to seasonal changes in the circulation of air masses. Past studies have shown that, in January, sources from the Atlantic Ocean prevail, and for the rest of the year, rainfall from the Mediterranean region dominates (Różański et al., 1993; Vreča et al., 2005; Mezga et al., 2014). Presently, the main source of moisture for the region is the southern-most Mediterranean, and a small amount of moisture is transported from the Atlantic (Vreča et al., 2005). The water from the Mediterranean source is enriched in δ18O in comparison to the Atlantic-derived moisture due to stronger evaporation effect in the Mediterranean basin. The processes that occur in karst environments can shape the δ18O of speleothems. Generally, in paleoclimatic interpretations, equilibrium isotopic fractionation is preferred because such fractionation only depends on temperature (Mangini et al., 2005). Progressive enrichments of δ13C and δ18O, along with the profile and/or a positive correlation of the two isotope compositions, indicate that fractionation occurred under kinetic conditions (Mickler et al., 2015). Additionally, another way to determine whether a speleothem was formed at isotopic equilibrium involves samples taken along the growth axis lacking a positive correlation in δ13C and δ18O values (Goede et al., 1986; Goede, 1994; Desmarchelier et al., 2000; Mickler et al., 2015). In the case of the whole RP flowstone, the compilation of δ13C and δ18O values showed no positive correlation, with R2 = 0.1235 (Fig. 6). However, there were three episodes in the profile where the R2 value (p < 0.05) was significantly greater: 632–954 ka, R2 = 0.4538; 1.246–1.456 Ma; R2 = 0.4944; 1.525–1.676 Ma, R2 = 0.6065 (Fig. 5B, C). The correlation of δ13C with δ18O may be caused by the effect of kinetic fractionation, but in some cases, it may be caused by climatic conditions (Dorale and Liu, 2009). The δ13C proxy reflects the level of soil development and the residence time of water in the karstic landscape. Both factors depend on the volume of precipitation (e.g., Lachniet, 2009).
Fig 6
Correlation of d13C and d18O values from the entire RP profile (p < 0.05).
Similarities of the RP Record to the Medi and LR04 Reference Curves
The marine data used by the authors are presented on an opposite scale to the speleothem data (see Section 4.4.). The RP profile is divided by the principal disconformity. The older part of the record can be better correlated with the Medi reference curve, while the younger part is more similar to the LR04 reference curve. This main difference between these correlations may reflect the long-term changes in the main source of vapour and atmospheric circulation in the region. However, the similarity of the segments to the chosen reference curves is not uniform (see Section 6.1). The affinity for the Medi reference curve to the RP record may be caused by the source effect. Rainfall over the cave area inherits the isotope signature of the northern Mediterranean Sea or the eastern Atlantic Ocean. The resulting δ18O recorded in speleothems is strongly influenced by marine reservoirs contributing vapour to rain formation (Kolodny et al., 2005; Burstyn et al., 2019). This also means that, over the age ranges in which the RP record is most similar to the reference curves, these sources were of the greatest importance in shaping the isotopic composition of the speleothems. Moreover, one should also consider why the rest of the record from the RP does not match the reference curves as well as the aforementioned segments of the profile. This discrepancy may be due to local factors affecting the isotopic composition, such as differences in local water circulation, mean annual air temperature and the amount of precipitation.
In addition, the youngest part of the RP record has been compared with speleothem records from other caves in the Mediterranean region: the Soreq cave and Peqiin cave in Israel (Bar-Matthews et al., 2003), Antro del Corchia cave located in Alpi Apuane karst in Italy (Drysdale et al., 2005) and Pentadactylos cave on Cyprus (Nehme et al., 2020). The noticeable similarities between the shapes of the RP δ18O and δ13C records to the Soreq, Peqiin, Antro del Corchia and Pentadactylos records during MIS 5e, which shows the climatic changes across the entire Mediterranean during the last interglacial period. In the δ18O curve, there is a noticeable difference in absolute values (~0.2‰), which may be attributable to the differences in the locations of the caves (i.e. the continental effect). Unfortunately, it is not possible to obtain a detailed correlation due to large differences in their resolutions, since one point represents the average time resolution of ~0.1 ka in the two Israeli speleothem records while representing ∼1 ka in the earliest part of the RP record.
Conclusions
The RP profile represents a unique data archive, especially for the length of record (approximately the last 3.2 Ma), as well as the uniqueness of the Mediterranean region. The results of this study supplement existing paleoclimate knowledge, as these types of profiles enable us to obtain high-resolution isotope records. Considering their long growth period, especially at the Pliocene–Quaternary boundary, it is truly an important source of climatological information.
At the time of correlation, 19 hiatuses were detected, which corresponded to clay layers present throughout the profile. The RP profile was deposited from ~3.2 Ma to ~0.08 Ma (MIS Km3 to MIS 5) and was divided into two segments separated by the principal disconformity. The lower segment correlated better with the regional curve, which suggests a greater influence of regional and local factors at the time of deposition in the cave area. The upper segment correlated better with the global LR04 curve, which may indicate a stronger influence of global changes at that time. Therefore, Slovenian caves may serve as vital sources of climatological information for much older periods and provide excellent materials for more detailed paleoclimate studies in the future.
Nevertheless, for a more accurate reconstruction of paleoclimate, a multi-proxy approach is necessary. However, the application of OIS, combined with magnetostratigraphy and U-Th data, enabled us to establish a more detailed chronology for the RP profile.
