The Potential of Tufa as a Tool for Paleoenvironmental Research—A Study of Tufa from the Zrmanja River Canyon, Croatia

Abstract

Tufa is a fresh-water surface calcium carbonate deposit precipitated at or near ambient temperature, and commonly contains the remains of macro- and microphytes. Many Holocene tufas are found along the Zrmanja River, Dalmatian karst, Croatia. In this work we present radiocarbon dating results of older tufa that was found for the first time at the Zrmanja River near the Village of Sanaderi. Tufa outcrops were observed at different levels, between the river bed and up to 26 m above its present level. Radiocarbon dating of the carbonate fraction revealed ages from modern, at the river bed, up to 40 kBP ~20 m above its present level. These ages fit well with the hypothesis that the Zrmanja River had a previous surface connection with the Krka River, and changed its flow direction toward the Novigrad Sea approximately 40 kBP (Marine Isotope Stage 3). Radiocarbon AMS dating of tufa organic residue yielded a maximum conventional age of 17 kBP for the highest outcrop position indicating probable penetration of younger organic material to hollow tufa structures, as confirmed by radiocarbon analyses of humin extracted from the samples. Stable carbon isotope composition (δ¹³C) of the carbonate fraction of (−10.4 ± 0.6)‰ and (−9.7 ± 0.8)‰ for the Holocene and the older samples, respectively, indicate the autochthonous origin of the carbonate. The δ¹³C values of (−30.5 ± 0.3)‰ and (−29.6 ± 0.6)‰ for organic residue, having ages <500 BP and >5000 BP, respectively, suggest a unique carbon source for photosynthesis, mainly atmospheric CO₂, with an indication of the Suess effect in δ¹³C during last centuries. The oxygen isotopic composition (δ¹⁸O) agrees well with deposition of tufa samples in two stages, the Holocene (−8.02 ± 0.72‰) and “old” (mainly MIS 3 and the beginning of MIS 2) (−6.89 ± 0.34‰), suggesting a ~4 °C lower temperature in MIS 3 compared to the current one.

1. Introduction

Tufa is a secondary carbonate wide spread within the karst areas worldwide. It is often the only available secondary carbonate within an area of important interest for paleoenvironmental investigation. Tufa precipitates from fresh waters (super)saturated in CaCO₃ at ambient temperature and usually includes the remains of micro- and macrophytes, invertebrates, and bacteria [1,2,3]. Dating of tufa and its use as paleoenvironmental tool has always been connected with difficulties. Researchers have attempted to obtain reliable dates from tufa samples using various approaches, such as establishment of the direction of tufa growth [4,5,6,7,8], radiocarbon dating after correction for the effects of “dead” carbon, or using U-Th dating. Additionally, tufa may erode, incorporate some allochthonous terrestrial material and “dead” carbonate, and experience recrystallization [9,10,11]. Reliable dating of tufa enables interpretation of the physico-chemical proxies and can provide many answers on Holocene and Pleistocene palaeoclimate [1,4,6,7,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. It is essential that the results obtained by different analyses can be put into a correct time frame.

During tufa formation, a certain amount of “dead” carbon is incorporated (quantified as the “dead-carbon proportion”, or DCP). This dead carbon originates from carbonate in the source bedrock (limestones and dolomites) which does not contain ¹⁴C. Bedrock dissolution from H₂CO₃ occurs due to the reaction of meteoric water with CO₂ produced from plant degradation in soils. The HCO₃^- ions formed by this dissolution contains both ¹⁴C-dead bedrock carbon and modern or recent ¹⁴C from plants, lowering the initial radiocarbon activity (a₀) of secondary carbonates precipitated from this water. Such waters emerging at springs are rich in HCO₃^- ions, and are characteristic of karst [26,27]. Measured radiocarbon dates appear older and must be corrected for DCP to obtain an accurate ¹⁴C age. The DCP can be determined since a₀ can be obtained from radiocarbon activity of freshly sampled: (1) dissolved inorganic carbon (DIC); (2) carbonate precipitated on an inert material; or (3) actively precipitating tufa. The first method is the most reliable for reasons stated below. It is important to emphasize, however, that it is not known if the DCP established in this manner is constant through time, as is the case with the marine reservoir effect (MRE [28,29], and references therein, [30]), although the MRE and DCP have different origins. The DCP is often used in discussion of speleothem ¹⁴C ages [31,32,33]. However, in secondary carbonates, such are tufa and lake sediments of the Croatian karst, DCP in precipitated carbonate may be higher due to the incorporation of weathered limestone fragments into the precipitate [34]. For these reasons, herein we will use terms “hard-water effect” or “fresh-water effect” rather than the DCP.

U-Th dating is an alternative to radiocarbon dating although it is not free of inherent limitations. It can be applied if a tufa sample is compact, hard, and/or crystalline [11,32]. If tufa is not compact, dating can be obtained by tufa leachate isochron method [12,35,36,37], but this technique does not necessarily solve the dating problem [38]. Recent work suggests many U-Th tufa dates are rejected and dating of both the carbonate fraction and organic residue from tufa by the radiocarbon method is the most reliable method [6,7,8], in spite of the difficulties with this latter approach [39]. It has been suggested that the only possibility of accurately dating secondary so-called “dirty” carbonates, such as tufas, is radiocarbon dating of the organic residue of the sample, since this part should not be influenced by erosion and recrystallization as is the case with the carbonate fraction [8]. Humin is generally regarded as the most resistant organic component and the best representative of true age, and is often used in geochronology [40,41,42,43,44].

Consequently, in this study radiocarbon dating of both carbonate and organic fractions has been implemented, as well as the stable isotopic composition of carbon (carbonate and organic fractions, extracted humin) and oxygen (carbonate fraction), in order to investigate the chronology of Zrmanja River Canyon development (Croatia, Dalmatian karst). Chronology will be discussed on the basis of 12 tufa samples found over a range of elevations from the river bed to 26 m above the stream level.

2. Materials and Methods

2.1. Site Description

The Zrmanja River is a relatively small river in Dinaric karst, Dalmatia, Croatia (Figure 1) with a total length of 72.5 km [45]. It receives water from many smaller springs and creeks on its route. The biggest contribution comes from the Krupa River, whose confluence with the Zrmanja River occurs 47 km from the spring of Zrmanja, downstream of Sanaderi. Similar to many rivers in the Croatian karst (e.g., Cetina [46]) according to its geomorphological properties the Zrmanja has a composite river valley characterized by an interchange of carbonate canyons and zones of lateral valley widening. The canyons are of ~20 m depth with variable widths. Generally, canyon slopes are rocky with poor vegetation (typically Mediterranean type), and at the Sanaderi location the tufa is partially concealed by grasses and trees.

The Zrmanja River changes its flow direction between the Palanka and the Sanaderi villages (Kravlja Draga) and turns from N-S to an almost E-W direction (Figure 1c). This is where the Zrmanja River may have had a surface connection with the Krka River (approximately 40 kyears ago [47]). The fossil river bed can be discerned in places while the groundwater connection with the Krka River has been shown from dye-tracing experiments performed during 1980s [48,49].

2.2. Sampling

Four tufa samples were collected at the Zrmanja river bed and eight fossil tufa samples were collected at the higher levels of Zrmanja River Canyon near Sanaderi (Figure 1) during sampling in 2015. Fluvial tufa samples were collected along the course of the river (samples #1–4), up to 2 m above water level (Table 1, Figure 1). Fossil tufa samples (#5–12) were collected at different levels above the river-bed (17–26 m) (Table 1, Figure 1) found on Cretaceous limestones and dolomites. Approximately 1 kg per sample was taken at each site, placed into a labelled plastic bag and transported to the laboratory.

Tufa samples from the investigated locations differ in color and morphology (Table 1, Figure S1). The majority can be classified as encrusted mossy deposits, except for samples #5–8, which can be classified as algal laminated crusts [51].

2.3. Measurement Methods

Specific radiocarbon activity a¹⁴C in carbonate samples/fractions was determined by the liquid scintillation technique (LSC). Between 50 and 70 g of sample was dissolved in HCl (18%) in an inert N₂ atmosphere and the obtained CO₂ was converted to C₆H₆ in a vacuum synthesis line. Benzene samples were measured on an LSC Quantulus 1220 [52,53].

The bulk residue after the carbonate dissolution was neutralized, dried, and prepared as graphite for a¹⁴C AMS measurements; herein it is referred to as the organic fraction. About 10–15 mg of each residue sample was combusted in pre-evacuated sealed quartz tube (850 °C). The humin samples were obtained from the bulk residue of the organic fraction by acid-alkali-acid method [54,55,56] and combusted in sealed quartz tubes. The obtained CO₂ was converted to graphite by Zn reduction [54]. Graphite samples were pressed into aluminum holders and sent to the accelerator mass spectrometry (AMS) facility at the Centre for Applied Isotope Studies (CAIS), University of Georgia, USA. An aliquot of CO₂ obtained from organic tufa residue combustion was vacuum sealed in a glass tube and sent to the same laboratory for δ¹³C determination. Graphite samples were measured by 0.5 MeV accelerator mass spectrometer (AMS) and δ¹³C was measured by isotope ratio mass spectrometer (IRMS) [57].

For both LSC and AMS techniques, Oxalic Acid II (SRM NIST4990C) was used as the primary reference material, and anthracite and Carrara marble as the background materials. As a control sample (a sample of well-defined a¹⁴C prepared as an unknown sample), humic acid T (a¹⁴C = 65.82 pMC) from the VIRI intercomparison [54] exercise and ANU sucrose (a¹⁴C = 150.6 pMC) were used for the AMS and LSC techniques, respectively.

Results of the radiocarbon measurements, normalized by δ¹³C to −25‰, are expressed as relative specific activity a¹⁴C (in units pMC, percent Modern Carbon) and as a conventional radiocarbon age (expressed in year Before Present, BP, where 0 BP = 1950 AD). Calibrated ages (cal BP) were obtained using OxCal software [58,59] and the IntCal20 calibration curve with included respective initial activities (hard-water corrections) for carbonate and organic fractions [60]. Correction for the hard-water effect is reported here either as initial activity a₀ or/and reservoir age defined as:

R = −8033 × lna₀.(1)

Stable oxygen and carbon isotope values (δ¹⁸O and δ¹³C) of carbonates were measured at the University of Melbourne, Australia using an AP2003 continuous flow mass spectrometer. Approximately 0.8 mg of carbonate sample was analyzed [61], with analytical uncertainties better than 0.10‰ for δ¹⁸O and 0.05‰ for δ¹³C.

3. Results

Results of radiocarbon dating and stable isotope compositions of tufa carbonate fraction of 12 samples from the Sanaderi location are shown in Table 2. Samples adjacent to modern river level (0–2 m) gave Holocene ages, while samples from tufa outcrops from the upper parts (elevations 17–26 m above the water level) gave radiocarbon conventional ages above 22,000 BP. δ¹³C values of the carbonate component range from −10.4‰ to −8.5‰, while δ¹⁸O values range from −7.4‰ to −6.4‰.

Results of radiocarbon dating of tufa organic fraction of samples from the Sanaderi location (except sample #9), shown in Table 3, yielded the Holocene dates for seven out of eleven samples. The organic residue of the sample #6 (12,060 ± 30 BP) was significantly older than sample #8 (7365 ± 30 BP) and corresponds to the end of the last glaciation, although it was sampled from approximately the same layer of the same tufa outcrop. Two samples from higher elevations, samples #11 and #12, were dated to 17,000 ± 40 BP and 15,980 ± 35 BP, respectively, i.e., within the glacial period. Four samples were selected for humin dating. The humin was extracted from samples #3 (Holocene), #6, #7, and #11 (older than the Holocene period, later referred as “old”). Two results from the humin component, compared to organic residue, gave slightly lower a¹⁴C values: the sample #3 and #7 for approximately 1.2 pMC and 2.2 pMC, respectively (Table 3). The humin extracted from samples #6 and #11, compared to organic residue, yielded higher a¹⁴C values for 19.6 pMC and 11.2 pMC, respectively (Table 3).

4. Discussion

4.1. Isotope Analyses of Carbonate and Tufa Ages

The obtained conventional ages of secondary carbonates in karst areas usually appear older than true due to the hard-water effect, i.e., reservoir ¹⁴C age. In such cases, the initial radiocarbon activity, a₀, should be determined [26]. In this work a₀ was determined from the a¹⁴C of dissolved inorganic carbon (DIC) in water from which tufa precipitates, giving the most reliable information among the possibilities described earlier. For the Sanaderi location, a¹⁴C of DIC was 84.5 ± 0.3 pMC resulting in reservoir age of 1350 ± 20 years [62].

The mean values of δ¹³C of tufa carbonates (−9.6 ± 0.8‰, Table 2) and δ¹³C of DIC (−11.4‰ [62]) indicate authigenic precipitation of carbonates from DIC. Both values are similar to the mean values of DIC and authigenic carbonate from another Dinaric karst system, Plitvice Lakes (−10.9 ± 1.3‰ and −8.9 ± 1.1‰, respectively) [34,63]. At the same time there is practically no difference in δ¹³C values between the Holocene (mean value −10.04 ± 0.57‰, n = 18) and older samples (“old”) (mean value −9.73 ± 0.78‰, n = 6) (Figure 2a). Here, the Holocene δ¹³C value data for the Zrmanja River tufa samples were taken from [64]. These values indicate that during the tufa deposition carbon source did not change. There is also no indication of a change in the hard-water effect in the past and we may assume that the a¹⁴C of DIC (84.5 ± 0.3 pMC) presents the initial activity for all tufa samples (Table 4).

The dating results of samples from higher positions at Sanaderi, samples #5–12 (Table 1 and Table 2), imply that tufa grew during Marine Isotope Stage 3 (MIS 3) and the beginning of MIS 2, a finding not previously found in Croatia. This would also mean that the Zrmanja River bed was approximately 20 m higher some 40 kyears ago compared to its current position. Furthermore, this is in accordance with the assumption that the Zrmanja River changed its watercourse approximately 40 kyears ago [47].

4.2. Possibility of Tufa Isotopic Composition of Oxygen as Paleoclimate Proxy

The oxygen isotopic composition of tufa samples analyzed within this work ranges from −7.4‰ to −6.4‰ (Table 2). Similar δ¹⁸O values of tufa from interglacial have been found in other regions of the Dinaric karst [2,3,4]. This indicates that all tufa samples from this work precipitated during warm periods. Considering the δ¹⁸O values, ¹⁴C dates of carbonate fractions are assumed correct, even for samples of ages close to the radiocarbon detection limit. Other indications of the relevant dates of the samples can be obtained from the difference of δ¹⁸O values between the Holocene samples (younger than 10,000 BP (11,700 cal BP)) and the “old” samples (older than 20,000 BP), Figure 2b. The Holocene samples, which have a mean value of δ¹⁸O = −8.02 ± 0.72‰ (n = 18, present data and data from [64]), precipitated from water of approximate temperature of 16.4 °C calculated from paleotemperature Equation (2) [65]:

T (°C) = 15.7 − 4.36 × (δ_c − δ_w) + 0.12 × (δ_c − δ_w)²,(2)

where δ_c is δ¹⁸O mean value vs. VPDB for either the Holocene or old carbonate samples, and δ_w is δ¹⁸O mean value vs. VSMOW of water from Miljacka Spring, −7.8‰ (Figure 1c).

Water temperature and δ¹⁸O data were obtained from one of the Krka River Springs, Miljacka Spring [66,67], which is completely fed by Zrmanja water [68,69]. The temperature range of Miljacka water obtained by 2-year seasonal measurement [66] was 8.8–14.5 °C so we may assume that Equation (2) [65] is applicable here. The difference between calculated and measured temperature values is expected since tufa precipitates more intensively during spring and summer [4] and due to cooling of Zrmanja water through the underground flow. The mean value of δ¹⁸O = −6.89 ± 0.34‰ (n = 6) for the “old” samples indicates their precipitation at approximate water temperature of 12.5 °C. The temperature difference of ≈ 4 °C between the Holocene and “old” tufa calculated here is similar to the temperature difference between the Holocene and MIS 3 presented in [70,71,72]. This further means that “old” tufa could have precipitated during MIS 3 and beginning of MIS 2.

The Holocene δ¹⁸O values obtained from the Zrmanja area (present results and data from [64]) are higher (between −9.5‰ and −7.0‰) than the Holocene authigenic carbonates from a well investigated area of the Plitvice Lakes system, where δ¹⁸O ranges between −11.9‰ and −8.7‰ [2,3,4,34,64,73]. This can be explained by environmental settings, and continental and altitude effects [4,64,74,75,76,77,78]. The Plitvice Lakes are at higher elevations (between 400 and 700 m asl) and are characterized by the moderate mountain climate (Cfb) compared to the Zrmanja locations with Csa climate type [79,80,81] and at lower elevations (<100 m asl). The Holocene tufa from the Krka River have δ¹⁸O values in the same range as tufa from Zrmanja (from −8.5‰ to −7.8‰) since both locations are situated within the same geographic area [4,66,67].

4.3. Tufa Age Implications from Isotope Analyses of Organic Residues and Extracted Humin

¹⁴C analyses of tufa organic residues gave higher a¹⁴C values compared to that of the carbonate fractions (Table 2, Table 3 and Table 4). The difference is expected due to a different intensity of the hard-water effect in organic and carbonate fractions [6,39,62,63]. In order to accept the obtained dates as reliable, the difference in a¹⁴C values between the organic residue and the carbonate should not be larger than 15 pMC [39,82]. However, recent authigenic and consolidated sediments from the Plitvice Lakes showed that the a¹⁴C differences between carbonate and organic fractions could be up to 20 pMC [63,83] and we consider the value of 20 pMC as a criterion for acceptable results. Six out of 11 samples have differences of a¹⁴C between organic residue (humin is not regarded here) and carbonate fraction of less than 20 pMC (Table 4).

The hard-water effect is also observed in organic matter since it is a mixture of different organics that photosynthesize from DIC and atmospheric CO₂ [26,62,84,85,86]. Age correction due to a hard-water effect for organic matter is usually lower than that in the carbonate fraction [26]. The δ¹³C values of organic residue in analyzed tufa samples are within the narrow range between −30.85‰ and −28.6‰ (Table 3) which indicates mainly constant origin of carbon in organic matter, namely that of atmospheric CO₂.

The hard-water effect of the organic matter was determined from a¹⁴C of the fresh moss sample, Eucladium verticillatum. The δ¹³C value of −34.2‰ and a¹⁴C value of 97.0 ± 0.9 pMC (Z-5948) of this moss indicate that atmospheric CO₂ is the main source of carbon with a possibility of incorporating DIC if flooded [62]. Additionally, the Suess effect [87], i.e., lower δ¹³C values of atmospheric CO₂ today than in the past, may also contribute to lower δ¹³C of moss. The Suess effect is confirmed by the clear difference between δ¹³C values for recent samples (<500 years: −30.5 ± 0.3‰, n = 3) and older samples (−29.6 ± 0.6‰, n = 8), as shown in Figure 3. There is also no difference between the Holocene (5–10 kyears: δ¹³C = −29.4 ± 0.4‰, n = 4) and the “old” samples (>10 kyears: −29.8 ± 0.6‰, n = 4). These results imply that the carbon source of organic residue has been constant in the past. Since the a¹⁴C of atmospheric CO₂ today is almost equal to 100 pMC [88,89,90] a₀ for organic residue is set as 97 pMC (reservoir age of 240 ± 70 years) for age calculation for all organic residue samples (Table 4). Because of the small value of the reservoir age for organic residue, the ages of dated samples did not change significantly (Table 4). As a consequence, the difference between corrected and calibrated carbonate and organic fraction ages diminished for the Holocene tufa (Table 4) and a good correlation between calibrated dates is obtained (Pearson’s r = 0.99, Figure 4a). On the contrary, there is no correlation between the ages of the carbonate and organic fraction in older samples (r = 0.05, Figure 4b). The absence of a correlation for the “old” samples (Figure 4b) and the large difference in a¹⁴C values of organic and carbonate fractions (Table 4) indicate that there should be an additional source of these effects. The ages of the organic fractions are especially suspicious since they suggest that tufa was mainly deposited during the MIS 2 cold period, [91,92,93,94] which is generally regarded as unfavorable for tufa precipitation [4,95].

A plausible explanation can be proposed after considering the humin analyses. The humin fractions extracted from samples #3 and #7 have almost the same values of a¹⁴C as the bulk organic residue (1.2 and 2.2 pMC difference, respectively) and their δ¹³C values were equal to those of the bulk organic measurements (Table 3). Humin extracted from samples #6 and #11 yielded higher a¹⁴C values for 19.6 and 11.2 pMC, respectively (Table 3). The differences in radiocarbon ages between the bulk organic residue and the humin component are probably due to the penetration of roots and plants through porous tufa after its deposition [38,39]. The sample #6 has the biggest difference between its humin and organic fraction (19.6 pMC), and it is the only sample with a significant difference in δ¹³C values of 0.9‰ (Table 3). Such a difference in δ¹³C can be expected if the Suess effect is considered (Figure 3) and this would indicate that very recent plants penetrated an old tufa deposit, shifting the organic a¹⁴C in the direction of a younger age for the tufa.

4.4. The Zrmanja River Canyon Development

Radiocarbon ages obtained on the different fractions (carbonate, organic) give two possible scenarios of tufa deposition and, consequently, Zrmanja River Canyon development. All radiocarbon data of tufa found just above the water level (samples #1–4) confirm tufa deposited during the Holocene (Figure 5). Data obtained from samples #5–12, found at elevations between 17 and 26 m above the present water level, differ significantly. Carbonate fraction dates suggest deposition of higher tufa samples during MIS 3 and the beginning of MIS 2 and formation of the canyon principally after MIS 3 (Figure 5). Organic residue results offer an alternative of canyon development in MIS 2 and the Holocene (Figure 5). In both cases, two groups are clearly distinguished: recent samples positioned just above the river bed, and older samples positioned at levels approximately 20 m above the river bed.

In the External Dinarides, the Last Glacial Maximum (LGM) was characterized by a cold climate, as indicated by numerous glacial morphological features, e.g., at the nearby Velebit Mt. found above 800–900 m asl [96,97,98,99] where average annual temperature was estimated to be 10 °C lower than today. Local conditions enable the formation of large snow packs and regional glaciations with estimated ice thickness of around 200 m [97]. Tufa formation is not expected in such cold-climate conditions. However, along the eastern Adriatic coast, growth of speleothems is recorded during the LGM [100,101,102]. During the LGM, conventionally 21 kBP, sea level was at its lowest at −121 ± 5 m below present [103]. Thus, a significant part of the Adriatic Sea area was dry land [100,104]. However, since that period the relative sea level (RSL) started to rise, so a change in fluvial system is expected, from incision toward deposition. Consequently, the incision of ~20 m between 8 kyears and 6 kyears as suggested using ¹⁴C dating of organic residue and humin (Figure 5) is unlikely.

During MIS 3, the RSL generally dropped [100,105,106]. Due to the erosional base level lowering, river valleys were generally getting longer, inducing the dominant trend of vertical erosion in the valleys. However, periods of RSL stability in the second part of MIS 3 [107,108] could have favored tufa formation. The further, rapid sea level lowering during early MIS 2 would have supported active fluvial incision processes.

Furthermore, the accumulation of tufa deposits in Mediterranean regions seems to start generally at the end of MIS 2 or at the MIS 2—MIS 1 transition, i.e., after 14 kyears in Spain ([109] and reference therein) and after 13 kyears on different locations in France ([110] and reference therein). Tufa deposition episodes during 32 ± 10 kyears were also defined on several locations in Spain [110]. Accordingly, we strongly support the period of tufa deposition at Sanaderi location principally during MIS 3.

Consequently, the two episodes of tufa formation during MIS 3 to beginning of MIS 2 and during MIS 1 (Middle and Late Holocene) reflect fluvial system response driven predominantly by climate change and RSL change. So, we assumed, as generally accepted, that tufa deposits formed during temperate/warm interglacial periods and interstadials (e.g., [109]), while downcutting occurred under rather colder conditions, i.e., at Sanaderi location, roughly between 25,000 and 6000 cal BP with the mean incision rate of around 1.1 mm/yr. However, incision could have also started earlier. All obtained results are close to the conclusions of earlier work of Fritz [47], who defined the main downcutting period in a single phase that lasted from 40,000 to 8500 BP.

5. Conclusions

Isotope analyses of tufa from the Zrmanja River near Sanaderi, Dalmatia, Croatia, are presented. For the first time, tufa from a period pre-dating the Holocene, was found in the Zrmanja River area and this is the first time that tufa from MIS 3 and the beginning of MIS 2 has been found in Croatia. The carbon stable isotope composition of the carbonate fraction indicates that tufa formation in all periods was authigenic, without incorporation of old detrital limestone and therefore suitable for paleoenvironmental research. Radiocarbon dating of tufa samples found at the river-bed level revealed similar ages of both organic and carbonate fractions (the Holocene). Radiocarbon ages of tufa samples found at higher positions (~20 m above the river bed) show significant discrepancies between carbonate and organic fractions. The organic residue dates (20–8 cal kBP) implies tufa deposition mainly during the MIS 2 period and partially in the early Holocene. Radiocarbon dates obtained from the carbonate fraction (42–25 cal kBP) indicate mostly MIS 3 deposition. The carbonate ages are argued to be more reliable since tufa deposition principally during MIS 3 is supported by the oxygen isotopic composition of the carbonate fraction. Organic-fraction dates are rejected since tufa deposition during the deep glacial of MIS 2 is highly unlikely. Additionally, humin dates indicate contamination of the organic fraction of tufa with younger plant-sourced ¹⁴C, as dates obtained from humin residues were younger than bulk organic residues, especially in the case of older tufa deposits. Finally, radiocarbon dates of the carbonate fraction of 20–40 kyr confirm previous estimates that the palaeo-Zrmanja River course differed from its current position (led probably into the Krka River) and that the new part of the Zrmanja river bed developed after 40 kBP, with a major incision period between 25,000 and 6000 cal BP, and an estimated incision rate of 1.1 mm/yr.

This research opens numerous possibilities for further investigation along the Zrmanja River and contributes to the understanding of regional paleoenvironmental changes.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/geosciences11090376/s1, Figure S1: Characteristic tufa samples from the Zrmanja River Canyon.

Author Contributions

Conceptualization, J.B., S.F. and I.K.B.; methodology, J.B., A.S., R.N.D. and D.B.; formal analysis, J.B.; investigation, J.B. and S.F.; resources, I.K.B.; data curation, J.B., A.S., I.L.M. and D.B.; writing—original draft preparation, J.B.; writing—review and editing, J.B., I.K.B., S.F., R.N.D. and A.S.; visualization, J.B., I.K.B. and I.L.M.; project administration, I.K.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Croatian Science Foundation grant number HRZZ-IP-2013-11-1623 (Reconstruction of the Quaternary environment in Croatia using isotope methods—REQUENCRIM).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All relevant data are presented in this paper.

Acknowledgments

The authors thank Janislav Kapelj, Maša Surić, Robert Lončarić, and Dražen Perica for help in sampling. The authors thank Nada Horvatinčić for continuous interest and help in the research of karst phenomena.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Geographical position of Croatia within Europe; (b) position of the investigated area within Croatia; (c) geological map of investigated area [50]; for full understanding of this study, Kravlja Draga, where the Zrmanja River changes its course, is also presented as well as the Krka River and Miljacka Spring; (d) enlarged area of the Sanaderi location, black dots ● mark sampling locations, and the numbers represent tufa sample numbers, see Table 1 for details.

Figure 2. Box-plot of (a) δ13C and (b) δ18O values for Holocene and “old” samples. Data for the Holocene tufa samples taken from Table 2 and [64].

Figure 3. δ13C values of tufa organic residues for different groups of ages.

Figure 4. Calibrated ages of the organic residue vs. calibrated ages of the carbonate component of the tufa samples from Sanaderi. (a) the Holocene samples; (b) the “old” samples.

Figure 5. Positions of tufa above river bed vs. calibrated ages of Sanaderi samples: carbonate ages are corrected for hard-water effect of 1350 years and organic and humin for 240 years.

Table 1

Geographical coordinates, approximate elevations, and height above the water level of samples at the Sanaderi location (Figure 1). Comments regarding microlocations and sample description are also given.

Sample No.	Geographical Coordinates	Elevation (m asl)	Height above Water Level (m)	Comment	Description
1	N 44°08′57″E 15°53′13″	52	0	Approximately at 1 m horizontal distance from river bed	Porous, moss ramains and soil incorporated, very hard, white—light yellow, not laminated, not clear if autochthonous or not
2		54	2	4 m horizontal distance from river bed	Porous, soil incorporated with lots of moss, very soft, reddish, not laminated
3		52	0	1 m horizontal distance from river bed	Porous, soil incorporated with moss, very soft, reddish, not laminated
4		53	1	Exactly below sample 2	Porous, hard, slightly laminated with moss overgrowth, grey—yellowish
5	N 44°09′02″E 15°53′12″	72	20	Approximately 50 m horizontal distance from Zrmanja R.	Nicely laminated, non-porous, brown and yellow laminae, moderately hard, compact
6	N 44°09′06″E 15°53′08″	73	21	Outcrop along the road, approximately 50 m horizontal distance from Zrmanja R.	Laminated, medium porosity, white and yellowish laminae, soft
7		72	20	1 m below sample #6, The same outcrop as #6
8		73	21	1 m horizontal distance from sample #6, from the same outcrop	Nicely laminated, non-porous, brown and yellow laminae, moderately hard, compact
9		72	20	1 m horizontal distance from sample #7 at the same height (approx. 20 m above water level) and 1 m below the sample #8 of the same outcrop along the road as #6	Laminated, medium porosity, white and yellowish laminae, soft
10	N 44°09′06″E 15°53′10″	69	17	Approximately 25 m horizontal distance from Zrmanja R. lBelow outcrop with samples #6–9	Non-porous, no visible remains of plants and soil, white-yellowish, laminated
11	N 44°09′06″E 15°53′07″	77	25	Outcrop approximately 100 m horizontal distance from Zrmanja R., bottom of outcrop	Porous, no visible remains of plants and soil, white-yellowish, soft
12		78	26	Top of the same outcrop as #11	Non-porous, no visible remains of plants and soil, white-yellowish, laminated, soft

Table 2

Results of isotopic analyses of tufa samples, a¹⁴C, δ¹³C and δ¹⁸O in carbonate fraction. Sample numbers defined in Table 1 and Figure 1.

Sample No.	Height above Water Level (m)	Lab. No. Z-xxxx	a14C (pMC)	Conventional Age (BP)	δ13C (‰)	δ18O (‰)
1	0	5951	40.35 ± 0.23	7290 ± 45	−8.96	−7.40
2	2	5952	89.5 ± 1.0	890 ± 95	-	-
3	0	5953	83.11 ± 0.55	1485 ± 55	−9.13	−7.42
4	1	5954	77.61 ± 0.60	2035 ± 60	-	-
5	20	5955	2.69 ± 0.12	29,040 ± 360	−9.85	−7.01
6	21	5956	0.78 ± 0.12	39,000 ± 1250	−10.33	−7.19
7	20	5957	5.29 ± 0.15	23,600 ± 230	−9.94	−7.25
8	21	5958	1.05 ± 0.13	36,600 ± 1000	−10.5	-
9	20	5959	2.03 ± 0.12	31,300 ± 480	-	-
10	17	5960	4.06 ± 0.14	25,740 ± 280	−8.47	−7.04
11	25	6295	0.81 ± 0.13	38,700 ± 1300	−10.39	−6.35
12	26	6296	6.46 ± 0.16	22,000 ± 200	−8.64	−6.49

Table 3

a¹⁴C, radiocarbon dates and δ¹³C values of the organic residues and the humin of tufa samples from the Sanaderi location.

Sample No.	Lab. No. Z-xxxx	a14C (pMC)	Conventional Age (BP)	δ13C (‰)
1	6458	51.8 ± 0.2	5285 ± 25	−29.96
2	6630	99.7 ± 0.3	25 ± 20	−30.85
3	6454	95.5 ± 0.3	370 ± 20	−30.27
3	7608, humin	94.3 ± 0.3	470 ± 20	−30.26
4	6631	95.7 ± 0.3	355 ± 20	−30.39
5	6455	26.7 ± 0.1	10,600 ± 30	−28.92
6	6459	22.3 ± 0.1	12,060 ± 30	−29.52
6	7610, humin	41.9 ± 0.2	6995 ± 45	−28.60
7	6460	29.3 ± 0.1	9870 ± 25	−28.80
7	7611, humin	27.1 ± 0.1	10,500 ± 30	−28.80
8	6461	40.0 ± 0.1	7365 ± 30	−29.38
10	6462	34.2 ± 0.1	8630 ± 30	−29.59
11	6456	12.1 ± 0.1	17,000 ± 40	−29.84
11	7609, humin	23.3 ± 0.1	11,700 ± 35	−29.93
12	6457	13.7 ± 0.1	15,980 ± 35	−30.70

Table 4

The a¹⁴C difference between organic/humin and carbonate fractions, ¹⁴C reservoir corrected calibrated ages of tufa carbonate fractions from the Sanaderi location (a₀ = 84.5 ± 0.3 pMC, R = 1350 ± 20 years), and calibrated ages of organic residues (a₀ = 97.0 ± 0.9 pMC, R = 240 ± 70 years).

Sample No.	a14C Difference (pMC)	Carbonate Fraction		Organic Residue and Humin
		Calibrated Age Span (cal BP) (68.3%)	Calibrated Age Median (cal BP)	Calibrated Age Span (cal BP) (68.3%)	Calibrated Age Median (cal BP)
1	11.45	6988–6874	6931	5869–5765	5811
2	10.2	1956–1957 cal AD (14.6%)2006–now cal AD (53.6%)	2007 cal AD	1956–1956 cal AD	1956 cal AD
3	12.39	275–160	218	196–104	148
3 humin	11.2			283–209	245
4	18.09	842–696	770	180–84	131
5	25.01	32,534–31,579	32,060	12,411–12,291	12,344
6	21.52	42,417–40,780	41,684	13,770–13,656	13,714
6 humin	41.1			7631–7526	7578
7	24.01	26,821–26,349	26,582	11,131–11,052	11,089
7 humin	22.3			12,277–12,102	12,188
8	38.95	40,791–39,291	40,042	7972–7884	7928
9	-	34,958–33,994	34,491	-	-
10	30.14	28,976–28,400	28,689	9449–9358	9397
11	11.3	42,277–40,558	41,501	20,366–20,225	20,294
11 humin	22.5			13,373–13,270	13,320
12	7.24	25,220–24,786	25,003	19,100–18,958	19,031