Presence and implications of petrogenic organic carbon in High Himalayan Crystalline lake sediment

Abstract

Twelve lacustrine sediment samples from a relict lake in the Kalla Glacier valley were co-dated using AMS radiocarbon (¹⁴C) and infrared stimulated luminescence (IRSL) dating methods. In general, the radiocarbon ages of bulk organic matter were older by a minimum of 1500 years compared to (age depth) modeled luminescence ages after fading corrections. This is observed for the first time in the lake sediments of High Himalayan Crystalline zone. A combination of lipid n-alkane data, Raman spectra and geochemical proxies suggested that this was due to ancient organic carbon (OC_ancient) that is a mixture of pre-aged (OC_pre-aged) and petrogenic (OC_petro) organic carbon within older glacial moraine debris that served as sediment source to the lake. Raman spectra suggest the presence of moderate to highly graphitized OC_petro in all the profile samples. The OC_petro contributed 0.064 ± 0.032% to the sediment and the lake stored 2.5 ± 0.7 Gg OC_petro at variable rates during the last 16 kyr, with the mean burial flux 160 kg OC_petro yr⁻¹. This study implies (1) employing another independent dating method in addition to radiocarbon method using bulk sediment organic matter, if the carbon content is low, to observe any discrepancy, and (2) a need to investigate on the fate of OC_petro as many such small lakes become relict in this region.

Introduction

The radiocarbon chronology of lacustrine sediment in the High Himalayan Crystalline (HHC) zone at ∼4000 m msl (above the treeline) has been, mostly, based on the bulk sediment organic matter (bSOM) (Beukema et al. 2011; Bhushan et al. 2018; Juyal et al. 2004; Kumar et al. 2020; Kusumgar et al. 2016; Srivastava et al. 2017). If any “old” carbon either of rock derived petrogenic or biospheric pre-aged is mixed with the lake sediment, the radiocarbon composition in the bSOM gets diluted and thus overestimate the radiocarbon age. The petrogenic organic carbon (OC_petro) that stored in source rocks of the catchment is liberated during physical erosion and transported by the melt water to the lake. The biospheric pre-aged organic carbon (OC_pre-aged) represents all other particles of organic carbon which are neither of contemporary vegetation (both terrestrial and aquatic), OC_bio nor of petrogenic (Galy et al. 2008). If OC_petro is defined as infinitely old organic carbon i.e., > 60 kyr or F_mod = 0 (Galy et al. 2008) and OC_bio as contemporaneity between the sediment depositional event and the vegetation in and around the lake, the age of OC_pre-aged is between that of OC_bio and OC_petro.

The Higher Himalayan range which is south to the South Tibetan Detachment System (STDS) is predominantly composed of metamorphic rocks whereas sedimentary rocks are at the north of STDS i.e., Tethyan Himalaya. There are two reports on the overestimation of radiocarbon ages of lacustrine sediment attributing to hard water effect from calcareous rock terrain of Tethyan Himalaya (Beukema et al. 2011; Juyal et al. 2004), but there is no report of such radiocarbon age overestimation from the crystalline HHC zone. The Himalayan lithologies contain some organic carbon that has been graphitized during the metamorphic process. This graphitization led to chemical and structural transformation of the organic carbon to crystallized carbonaceous mater (Galy et al. 2011). There are reports on the contribution of the petrogenic organic carbon (OC_petro) to Himalayan rivers (0.03–0.05% of OC_petro) and marine sediments (0.02–0.03% of OC_petro) sourced from HHC zone (Galy et al. 2011) but no reports on the lake sediments. Along with OC_petro, organic carbon older than the contemporary vegetation (OC_pre-aged) with an average age of over 15 kyr (after accounted for the OC_petro) that represents up to 20% of total biospheric OC is stored in Ganges-Brahmaputra basin (Galy and Eglinton 2011).

When transported from land to sea by rivers, OC is vulnerable to oxidization and subsequent release of CO₂ to the atmosphere (Copard et al. 2022; Luo et al. 2022). The vulnerability of oxidization of OC depends on whether it is labile or refractory OC (Petsch 2014). Even graphitized OC_petro, likely poorly graphitized due to low- to high-grade metamorphism (Beyssac et al. 2003), can be oxidized by 70% of it during fluvial transport from the Himalayan region to Bengal fan (Galy et al. 2008). The state of pre-aged OC (OC_pre-aged) is worse than this, because it is less refractory and hence is more vulnerable for oxidation during fluvial transport to the Bengal fan (Galy and Eglinton 2011). Hence, understanding the contribution of burial flux and characteristics of OC_petro in the high altitude lake sediments, and the role of these lakes itself, become important to study the exchange of CO₂ between lithosphere, hydrosphere and atmosphere. So, besides overestimation of radiocarbon age, the presence of OC_petro and OC_pre-aged in lacustrine sediments can be used to understand the continental storage of ancient carbon.

To test our hypothesis that the presence of OC_petro and OC_pre-aged in lake sediments would not only imply overestimation of radiocarbon age, but also on the continental storage of ancient carbon, sediments from a relict lake in HHC zone were dated using radiocarbon and luminescence methods with the later providing a control for the age (Figure 1).

Figure 1. (a) Study area within the High Himalayan Crystalline zone is shown by yellow rectangle on Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) of Uttarakhand (Source: NASA earth data, https://www.earthdata.nasa.gov/). (b) Orographic and precipitation trends from south to north along the study area (horizontal gray rectangle). Precipitation data are from TRMM (Tropical Rainfall Measuring Mission). (c) Location of the study site (red hollow rectangle), traced on Google Earth Pro Imagery (https://earth.google.com/web/) and the catchment area along with the glaciers. (d) Field photograph showing the study site and the moraines. (e) Field photograph of the pit profile used in the present study. The inset shows close up of the sediments used for radiocarbon and luminescence dating (Geological map of the study area is provided as supplementary Figure S1.).

Study area

A relict morainal lake in the Kalla Glacier valley in the Kunti Banar River basin (30°38ʼ13ʼʼN, 79°54ʼ08ʼʼE) (Ali et al. 2020) was revisited for a high-resolution past-hydroclimate reconstruction. It is located in the HHC, between the Main Central Thrust (MCT; in the south) and STDS; in the north) (Ali et al. 2020). Geologically, the catchment area comprises of calc-silicate, quartzite, pelitic gneiss, biotite schist, psammitic gneiss, migmatite, and sillimanite gneiss, and belongs to the Vaikrita group of the HHC (Sinha 1989). The study site lies in the alpine meadow zone, above the contemporary tree line (∼4000 m asl) and is situated on a moraine deposit. The average temperature and precipitation are –9.5°C and 647 mm, respectively (Ali et al. 2020) as per CRU-TS-4.02 (Harris et al. 2017). The majority of the precipitation (∼80%) is contributed by the Indian Summer Monsoon (ISM) with a smaller contribution by winter western disturbances (Figure 1). This lake was fed by the melt water and surface runoff from the catchment (Ali et al. 2022).

Materials and methods

Twelve sediment samples each for radiocarbon dating and luminescence dating were collected from a 4-m trench from the middle of the relict lake as shown in Figure 1. Samples for luminescence dating were collected in galvanized iron (GI) pipes hammered horizontally into the strata at different depths and the sediment inside the pipes were tightly packed so that no mixing of samples happened within the pipes. Based on the field observations, samples containing relatively larger portion of organic matter were collected for radiocarbon dating in zip-locked polythene bags. The samples for both dating methods were collected from nearly the same horizon/strata of the pit profile of 4 m deep. In addition, two modern surface samples near the pit were collected for radiocarbon dating.

In the following sections, we summarize of all the methods that were employed for this study. The sequence of methods are radiocarbon dating, luminescence dating, bulk sediment stable carbon isotopic analysis, C/N ratio measurement, palynofacies analysis, loss on ignition (LOI) measurement, n-alkane data, Raman micro-spectrometry, major oxides using XRF, mineralogical composition using XRD, particle size analysis and statistical analysis.

Radiocarbon dating

The samples were treated using acid–base–acid (ABA) method. Oven-dried samples were treated with 0.5 M concentrated HCl at 60°C for 10 hr to dissolve the carbonates. This same acid treatment was repeated with lesser time (2 hr) after the samples were treated with 0.1 N NaOH at 60°C for 3 hr to remove humic acid. Then, the samples were combusted to liberate CO₂ which was then transformed to graphite in the presence of hydrogen at 600°C using an AGE3 graphitization system with iron as a catalyst. Being bulk sediment organic matter, the base treatment (0.1 N NaOH) is redundant, however, to avoid any overestimation in the radiocarbon ages it was included (Blyakharchuk et al. 2017). Radiocarbon concentration was measured using AURiS (Accelerator Unit for Radioisotope Studies) at the Physical Research Laboratory (PRL), Ahmedabad, that is equipped with a compact 1 MeV Accelerator Mass Spectrometer (AMS) of HVEE, the Netherlands. International standards such as IAEA-C1 and IAEA-C2 (carbonates) and FIRI-E and VIRI-U (humic acid) were used to calibrate/standardize the instrument (Bhushan et al. 2019). The measured radiocarbon ages were calibrated using Bchron package and INCAL20 in R environment. Samples were prepared in Birbal Sahni Institute of Palaeosciences, Lucknow and the radiocarbon samples were measured at the PRL, Ahmedabad.

Luminescence dating

The luminescence ages (LAs) were calculated using 2 parameters, namely the equivalent dose (D_e in Gy) and the dose rate (Gy.kyr⁻¹). The D_e was measured using the fine grain (4–11 μm) polymineral fraction from all the 12 samples. The fine grain fraction was chosen because they would have been well bleached (reset by solar exposure) as they spent relatively longer time in suspension compared to its coarser counterpart. The fine grain polymineral fraction was extracted from the unexposed samples from central portion of sampling pipes and deposited in cleaned Al discs of 9.6 mm diameter following standard protocol (Morthekai and Ali 2014). The single aliquot regeneration (SAR) procedure was followed to estimate D_e and fading rate using infra-red (870 ± 40 nm) stimulated luminescence signals (IRSL) measured at 50°C after preheating to 250°C at 2°C/s with a holding time of 60 s. The emitted photons were detected in violet-blue wavelength (320–460 nm) range. Beta particle irradiations were carried out using an on-plate ⁹⁰Sr/⁹⁰Y beta source and it delivered a dose rate of 0.06 Gy.s⁻¹. Fading rate was measured (Auclair et al. 2003) and fading correction was done using Huntley and Kars method (HKM) (Huntley 2006; Kars et al. 2008) as implemented in Luminescence package (Kreutzer et al. 2018; Kreutzer et al. 2012). The concentration of U, Th and K were estimated using high pure Ge gamma spectrometer by comparing the measured gamma spectra against that of NUSSY standard (Preusser and Kaspar, 2001). Radioactive disequilibrium of U-series was examined by comparing the gamma radioactivity of ²²⁶Ra, ²¹⁴Pb and ²¹⁴Bi. Measured concentrations of U, Th and K were converted into dose rate values after accounting for the measured water content in the sample matrix and cosmic ray dose rate using online Dose Rate and Age Calculator, DRAC (Durcan et al. 2015). Gamma detector efficiency corrected count rate (cts.m⁻¹.mg⁻¹) of ²²⁶Ra (energy = 186 keV; efficiency = 0.11), ²¹⁴Pb (295.2 keV; 0.19) and ²¹⁴Bi (609.3 keV; 0.21, 1120.3 keV; 0.05, and 1700 keV; 0.03) was compared to check for disequilibrium in U-series. All the luminescence measurements were carried out in the Luminescence Dating Laboratory, Birbal Sahni Institute of Palaeosciences, Lucknow.

Stable carbon isotope measurement and analysis

The composition of stable carbon isotope (δ¹³C, ‰) was measured in isotope ratio mass spectrometer with Elemental Analyzer (Flash EA 2000 HT) through an auto sampler. Samples (12 pit samples, and 29 modern plant samples) were combusted and thus produced CO₂ was introduced into the Continuous Flow Isotope Ratio Mass Spectrometer (CFIRMS, MAT 253) coupled with a Con-Flow IV interface for isotopic analysis. The instrument has been calibrated using IAEA CH₃ and CH₆ with the accuracy of ±0.1‰ (1σ) for the CO₂ measurement, and carbon isotopic data have been reported against Vienna Pee-Dee Belemnite, VPDB. Total organic carbon (TOC) was calculated from the peak area obtained from the sum of the integrated m/z of 44, 45, and 46 signals measured in the CFIRMS24. All samples were analyzed in the Stable Isotope Laboratory, Birbal Sahni Institute of Palaeosciences, Lucknow.

The carbon isotopic composition of individual n-alkanes was measured using Trace GC Ultra (Thermo Fisher Scientific) coupled with a MAT-253 IRMS linked via a GC Isolink (combustion interface) and Thermo Fisher Scientific Conflo IV interface. Instrument performance was routinely checked using international standards A5 (C₁₆–C₃₀) and Fluka alkane mixture with known δ¹³C values. The reproducibility in the international standard mixtures A5 and Fluka alkane mixture was found to be ±1.0‰ and ± 0.4‰ for carbon isotope. The samples were analyzed in duplicate, and the mean carbon isotope ratios are reported with respect to VPDB. The extraction and compound-specific isotope analysis of n-alkanes were conducted in the Stable Isotope Laboratory of the Indian Institute of Science Education and Research Kolkata.

Extraction, identification, and quantification of n-alkanes

The detailed methodology for the extraction and isotope analysis of n-alkanes from 4 modern plants, 2 surface sediments (Figure S5 a), and pit (soil) samples (KBPR 01-12; Figure S5 b) is outlined in Sarangi et al. (Sarangi et al. 2019; Sarangi et al. 2022). Briefly, dried, crushed, and homogenized plant and pit samples were extracted with an accelerated solvent extractor (Dionex ASE-350, Thermo Fisher Scientific) using a mixture of dichloromethane and methanol (93:7) to obtain the total lipid extract (TLE). Consequently, the non-polar hydrocarbon fraction (n-alkane) was separated by short column silica gel chromatography. To check the recovery during the TLE extraction process, a few of the samples were spiked with an internal standard 5α-androstane. The n-alkanes were analyzed by gas chromatography (7890A GC system; Agilent Technologies) equipped with a flame ionization detector (FID). Individual n-alkanes were identified based on the characteristic retention time (RT) obtained from the calibration standards SUPELCO C₈–C₄₀ alkane and Fluka alkane mixture (C₁₀–C₄₀). The concentration of the 5α-androstane showed a recovery of >85%, and the replicate measurements of the calibration standards showed an uncertainty of <4% in the relative concentrations of n-alkanes.

The n-alkanes were characterized using ratio between the concentration of higher (≥C₂₅) and lower-chain (≤C₂₄) n-alkanes (HC/LC), carbon preference index (CPI) of long-chain n-alkanes (C₂₅-C₃₇; CPI_long), average chain length (ACL) of long-chain n-alkanes (C₂₃-C₃₇; ACL_long), and proxy ratios between terrestrial/aquatic plant groups (P_aq) values. The relative concentrations of the n-alkane homologues were used to calculate the CPI_long, ACL_long, and P_aqvalues. The CPI_long(Marzi et al. 1993), ACL_long(Eglinton and Hamilton 1967), and P_aq values (Wang et al. 2014) were calculated using Eqns. 1–4, which are as follows:

(1) $${\rm{HC}}/{\rm{LC}} = {\rm{ }}\Sigma {{\rm{C}}_{{\rm{25}} - {\rm{37}}}}/\Sigma {{\rm{C}}_{{\rm{13}} - {\rm{24}}}}$$

(2) $${\rm{CP}}{{\rm{I}}_{{\rm{long}}}} = {\rm{ }}0.{\rm{5 }} \times {\rm{ }}\left[ {\left( {\Sigma {{\rm{C}}_{{\rm{25}} - {\rm{37}}}}/\Sigma {{\rm{C}}_{{\rm{24}} - {\rm{36}}}}} \right){\rm{ }} + {\rm{ }}\left( {\Sigma {{\rm{C}}_{{\rm{25}} - {\rm{37}}}}/\Sigma {{\rm{C}}_{{\rm{26}} - {\rm{38}}}}} \right)} \right]$$

(3) $${\eqalign{{\rm ACL}_{\rm long} = & {\ }(23 \times {\rm C}_{\rm 23} + 25 \times {\rm C}_{\rm 25} + 27 \times {\rm C}_{\rm 27} + 29 \times {\rm C}_{\rm 29} + 31 \times {\rm C}_{\rm 31} + 33 \times {\rm C}_{\rm 33} \cr & + 35 \times {\rm C}_{\rm 35} + 37 \times {\rm C}_{\rm 37})/({\rm C}_{\rm 23} + {\rm C}_{\rm 25} + {\rm C}_{\rm 27} + {\rm C}_{\rm 29} + {\rm C}_{\rm 31} + {\rm C}_{\rm 33} + {\rm C}_{\rm 35} + {\rm C}_{\rm 37})}}$$

(4) $${{\rm{P}}_{{\rm{aq}}}} = {\rm{ }}\left( {{{\rm{C}}_{{\rm{23}}}} + {\rm{ }}{{\rm{C}}_{{\rm{25}}}}} \right)/\left( {{{\rm{C}}_{{\rm{23}}}} + {\rm{ }}{{\rm{C}}_{{\rm{25}}}} + {\rm{ }}{{\rm{C}}_{{\rm{29}}}} + {\rm{ }}{{\rm{C}}_{{\rm{31}}}}} \right)$$

Raman micro-spectrometric measurements and analysis

Raman spectra were measured from all the pit samples (n=12) that were prepared for radiocarbon dating using AMS. The excitation wavelength was 780 nm with a beam width of 3 nm and < 5 mW power. Wherever the black opaque materials were (at ×50 magnification), the excitation beam was focused on those spots as we were looking for graphite materials. The grating was having 1200 grooves/cm and the resolution was 1.02 cm^–1. R₁ ratio was calculated using the peak intensity ratio of a physico-chemical defect band, D₁ (1320 cm^–1) and graphite band, G (1580 cm^–1). Two artificial graphite were also used as reference. All samples were measured in the Raman Spectrometer, Birbal Sahni Institute of Palaeosciences, Lucknow.

C/N ratio measurements

Nearly ∼1 gm of 12 pit samples were taken after coning and quartering to measure C/N ratio. Samples were dried at 60°C and ground to fine powder using a mortar and pestle in order to be homogenized. A subsample was used directly in tin cups for combustion and determination of total nitrogen (TN). For total organic carbon (TOC), a subsample was oven dried, crushed and decarbonated using 0.5 M concentrated HCl at 60°C for 10 hr and taken in silver cups for the combustion. Measurements were made in an in-line isotope ratio mass spectrometer (IRMS) for measuring stable carbon and nitrogen isotopes (δ¹³C, δ¹⁵N) with EA-IRMS-AGE combo. The instrument has been calibrated with various standards of carbon (IAEA 601, ACA, Ox-II and USGS 24) and nitrogen (IAEA N1, IAEA N2, ACA and USGS 26) where the measured and reported values of respective elements were consistent with R² = 0.999. Carbon and nitrogen contents (TOC wt. % and TN wt. %) of samples were measured using an elemental analyzer (EA1112; Thermo®) (Agnihotri et al. 2020) of Birbal Sahni Institute of Palaeosciences, Lucknow.

Palynofacies analysis

Palynofacies analysis was carried on the 12 samples that have been used for the radiocarbon dating. Approximately five gram of each sample was treated with 30% HCl and 35% HF to remove the carbonates and silicates, respectively. The sieved samples (25 μm) were used to make strewn slides using polyvinyl alcohol and mounted using Canada balsam. The slides were scanned under the Olympus BH-2 microscope and DP-25 camera was used to capture the palynofacies. All samples were chemically treated and scanned in the Birbal Sahni Institute of Palaeosciences, Lucknow.

Loss on ignition (LOI) measurement

Five grams of the left over luminescence samples from the pipes (n = 12), were used to measure LOI. The moisture content (110°C for 12 hr), organic matter (550°C for 2 hr) and the inorganic content (after heating to 950°C for 2 hr, CaCO₃, %) present in the sediments were removed by the sequential heating in a muffle furnace. Then, the LOI was calculated

(5) $$LO{I_{110}}\left( {\rm{\% }} \right) = {{{{n_0} - {n_1}}}\over{{{n_0}}}} \times 100$$

(6) $$LO{I_{550}}\left( {\rm{\% }} \right) = {{{{n_0} - {n_2}}}\over{{{n_0}}}} \times 100$$

(7) $$LO{I_{950}}\left( {\rm{\% }} \right) = {{{{n_0} - {n_3}}}\over{{{n_0}}}} \times 100$$

(8) $$LOI\left( {\rm{\% }} \right) = {\rm{average\;of\;LO}}{{\rm{I}}_{110}},{\rm{\;LO}}{{\rm{I}}_{550}},{\rm{\;and\;LO}}{{\rm{I}}_{950}}$$

(9) $$CaC{O_3}\left( \% \right) = \,1.36\times {\rm{LO}}{{\rm{I}}_{950}}$$

where n₀, n₁, n₂ and n₃ are the weight of samples (no heat), after heating to 110°C, 550°C and 950°C, respectively (Dean 1974). The factor 1.36 is the weight ratio of carbon dioxide (CO₂) and carbonate (CO₃ ²⁻).

XRF measurements and calculation geochemical indices

XRF measurements were made on the pressed powder pellets of KBPL samples (n = 12; luminescence dating samples) using wavelength dispersive X-Ray Fluorescence (WD-XRF) Model: Axios Max; 4 KW, PANalytical make at BirbalSahni Institute of Palaeosciences, Lucknow. The instrument has been calibrated with the international standard, QLO1, with the accuracy of measurement is better than 2–5% and precision < 2%. Chemical alteration index (CIA, %) and index of compositional variability (ICV) were calculated (Armstrong-Altrin 2015; Baiyegunhi et al. 2017; Cox et al. 1995; McLennan 1993; Nesbitt and Young 1989) as follows:

(10) $${\rm{CIA}}\left( {\rm{\% }} \right) = {{{{\rm{A}}{{\rm{l}}_2}{{\rm{O}}_3}}}\over{{{\rm{A}}{{\rm{l}}_2}{{\rm{O}}_3}} + {\rm{CaO}} + {\rm{N}}{{\rm{a}}_2}{\rm{O}} + {{\rm{K}}_2}{\rm{O}}}} \times 100$$

(11) $${\rm{ICV}} = {{{{\rm{F}}{{\rm{e}}_2}{{\rm{O}}_3} + {\rm{K}}2{\rm{O}} + {\rm{N}}{{\rm{a}}_2}{\rm{O}} + {\rm{CaO}} + {\rm{MgO}} + {\rm{MnO}}}}\over{{{\rm{A}}{{\rm{l}}_2}{{\rm{O}}_3}}}}$$

XRD measurement and analysis

XRD measurements were made in PANalytical Xpert’3 Powder irradiated with Cu target (Kα₁ = 1.540598 A° and Kα₂ = 1.544426 A°) with a continuous scan type. The diffraction pattern was recorded at room temperature in the 2θ range from 5° to 90° (Chaddha et al. 2021). These measurements were done on 3 sets of 12 samples each. The first set was the grains that were predominantly of feldspar as they had been magnetic separated (Frantz model) at higher current, 1.5 A (Porat 2006). The second set was the samples that were prepared for radiocarbon measurements in AMS. The third set was bulk samples. The background subtracted XRD spectra were analyzed by comparing against rockjock mixtures dataset using powdR package (Butler and Hillier 2021). XRD measurements were done to quantify the plagioclase, illite and graphite in the samples.

Particle size distribution and end member analysis (EMA)

Particle size distribution analyses were done on the same sediment samples for which the LAs and RAs were determined. The samples were treated with 10% HCl and 3% H₂O₂ to remove carbonate and organic carbon, respectively (Andreola et al. 2004; Battarbee et al. 2001; Schumacher 2002; Vaasma 2008). To avoid grain flocculation during pretreatments, 1% solution of sodium hexametaphosphate was also used (Andreola et al. 2004; Murray 2002). We measured particle size distribution using a Laser Diffraction Particle Size Analyser (Beckman Coulter LS™ 13 320). End member analysis (EMA) was performed with AnalySize v1.2.0 (Paterson and Heslop 2015).

Statistical analysis

Welch two sided t-test was done to check whether the intermediate samples, between depth 180 cm to 235 cm, are significantly different than the rest of the profile in terms of the measured/derived variables’ means. This was performed using t.test() method in R environment (R Core Team 2022).

Results and discussion

The calibrated radiocarbon ages, RA_cal and other associated parameters regarding radiocarbon ages are given in Table 1. The range of RA_cal with the confidence interval (CI) are also given. The average RA_cal was calculated from the expected value of the probability density function computed in between the range of RA_cal, <RA_cal>. The uncertainty was calculated from the differences between the higher range of RA_cal and <RA_cal>, and between <RA_cal> and the lowest range of RA_cal. Among these 2 differences, the maximum value was used for upper error, ue, and the minimum value was for lower error, le (Data S1). The RA_cal will be called RA hereafter

Table 1. The AMS measured radiocarbon concentration in terms of fraction modern carbon (F_mod), uncalibrated radiocarbon age (RA, year), expectation value of calibrated radiocarbon age <RA_cal> with lower error (le) and upper error (ue), and the range of calibrated (IntCal20) radiocarbon ages are tabulated. <RA> was calculated from the higher contributing peak of the probability density function that was outputted by Bchron package after calibration using BchronCalibrate() method. The δ¹³C (‰) values that were used to correct for the radiocarbon ages, and TOC (%) are also given here. For the calculations and plotting of RA_cal the larger error among ue and le are used

Luminescence ages were measured using IRSL arising from fine (4–11 μm) polymineral grains. The measured D_e (Gy) and dose rate (Gy.kyr^–1) are given with the depth of the profile (Figure 2a). Combined dose response curves (standardized dose response curve) and sensitivity corrected IRSL signals (L_n/T_n’s) of all samples are given in Figure 2b. Measured, un-faded and fading correction dose response curve of KBPL 01 are shown in Figure 2c. The measured luminescence ages were then fading corrected, and the fading corrected age (LA_measured; plotted as filled black spheres in Figure 3b) are given in Table 2 with other required parameters ( Data S2). The RA (RA_cal) are also plotted in Figure 3b (red filled spheres) and the litholog is given in Figure 3a (for a detailed lithology description, please see Figure S2). With the exception of sample at 55 cm, 180 cm and 235 cm, other nine RAs are higher than LA_measured (Figure 3b). Both the RA and LA_measured of three samples at depth from 180 cm to 230 cm (intermediate depth/samples) were older than the downward samples.

Figure 2. (a) The variation of equivalent dose (D_e, Gy) and dose rate (Gy.kyr^–1) with the profile depth (cm). (b) Combined standardized growth curves (SGC) of all samples (n = 12) with their respective natural L/T’s give an overall picture about the luminescence characteristics of the samples dated. L_n/T_n’s of KBPL 01 and KBPL 02 and the three intermediate samples (KBPL 05, KBPL 06 and KBPL 07) are meeting the SGC after the characteristic dose (D₀ = 198 Gy). (c) Fading correction was carried out using Huntley’s method (Huntley 2006) as implemented by Kars et al. (Kars et al. 2008). The growth curves (GC) of measured, unfaded (corrected for the laboratory irradiation time) and faded (calculated with natural dose rate) with its natural L/T’s of KBPL 01 are shown as representative of all samples.

Figure 3. (a) Litholog of the lake sediment profile. A detailed description of the lithology is given in Figure S2. (b) The measured luminescence ages (LA_measured, black filled spheres), and calibrated radiocarbon ages (RA, red filled spheres) are plotted with depth. The vertical shaded region indicates the moraine ages (23 ± 4 kyr) on which the lake is situated.

Table 2. The concentration of radioactive nuclides (U, Th and K), water content (WC), alpha efficiency (a-value), dose rate, equivalent dose (D_e), observed over-dispersion (OD, %), number of aliquots used (n), fading rate parameters (g_tc=2days-value and ρ′) and the fading corrected age (LA_measured) are tabulated. All the other parameters that are used to estimate luminescence ages are given in Data S2

Because of the fact that the studied relict lake is situated on a ∼24 kyr old moraine (Ali et al. 2022), those 3 ages of the intermediate samples were reasonably assumed to be overestimated. The observations of the presence of EM 4 and EM 5 (Figure 4b; orange bar at 180 and 235 cm and red bar at 205 cm; EMs are characterized as shown in Figure 4a) and unexpectedly too large—almost saturated—accumulated luminescence/dose (Figure 2a,b) in the intermediate samples support our assumption that these intermediate samples would have poorly bleached and hence overestimated LA_measured. The large equivalent dose is most probably the unbleached dose than the burial dose because the dose rate for the intermediate samples is similar to the other samples (Figure 2a). In addition to that, as the overestimated ages were equivalent to moraine ages (∼24 kyr) where that lake is situated, it might be possible that the sediment were eroded from the moraine and deposited quickly hence did not get enough time to bleach. Testing whether a fine grain polymineral sample was poorly bleached or not is extremely difficult because each aliquot has 1000s of grains in it and thus gross averaged. So, we have provided circumstantial evidence that suggest those 3 intermediate samples, and the sample at 55 cm depth, would have been poorly bleached and hence luminescence ages were overestimated.

Figure 4. (a) End member analysis (EMA) decomposed five characteristic EMs from the grain size distribution data. (b) EM1 and EM2 have their modes developed in silt domain (mean EM1: 14 µm and EM2: 27 µm respectively) whereas rest three have their modes in sand domain (mean size EM3: 64 µm, EM4: 91 µm, and EM5: 151 µm). Out of five EMs, EM4 is least sorted (high Standard Deviation) and is characterized by broad and multi-modal distribution pattern. EM5 is fine-skewed and have a narrow and well resolved peak in medium sand domain. EM1, EM2 and EM3 are mostly unimodal having relatively broader peaks. The effectiveness of sediment sorting can be ascribed by the shape of individual EM distribution. Broader and poly-modal peaks are associated with poor sorting than the narrow and unimodal one.

For a further comparison of both age estimates, LAs were calculated using Bayesian age-depth modeling (by excluding the above mentioned 4 overestimated LA_measured as shown in Figure S4c) for the depths of RA. RA_offset were calculated by subtracting modeled luminescence ages (LA_modeled = LA) from the respective RA (red filled spheres, Figures 3b and 5a) for all the samples. Using the model, the radiocarbon ages of all the samples were overestimated and RA_offset (RA–LA) vary from 1.5 kyr to 15.6 kyr (Figure 5b; Table 3).

Figure 5. (a) LAs were calculated using Bayesian age-depth modeling (by excluding the above mentioned 4 LA_measured) for the depths of RA. RA_offset were calculated by subtracting modeled luminescence ages (LA_modeled) from the respective RA (red filled spheres) for all the samples. LA_modeled (LA) is also plotted (blue filled spheres and line). (b) RA_offset (RA–LA) with depth is given.

Table 3. The calibrated RAs (mid-point of the range of calibrated radiocarbon age) and the modeled luminescence ages (LA_modelled/LAs) for the same depth as of RAs are given. The difference between RA and LA_modelled is radiocarbon age offset (RA_offset). The usage of LA_modelled, not LA_measured (as shown in Table 2), to calculate RA_offset is discussed in the caption to Figure 5

Were luminescence ages underestimated?

The important part of this work is to ascertain whether LAs are underestimated or RAs are overestimated. Luminescence ages can be underestimated if the equivalent dose (D_e) is underestimated or dose rate overestimated. Equivalent dose will never be underestimated but only be overestimated if (1) resetting/bleaching of luminescence signal was heterogeneous before deposition, or/and (2) anomalous fading rate (g-value or ρ′) is large (due to its logarithmic relationship larger the fading rate, larger the corrected age. So, a small change in the fading rate will produce not only ages with large uncertainty but also spuriously older ages). Dose rate can be overestimated if (1) the radioactive nuclides (U, Th, and K) are overestimated, or/and (2) alpha efficiency (a-value) is overestimated.

Regarding underestimation of D_e, it is highly unlikely a poor bleaching (presence of unbleached geological dose) would underestimate, because it can only overestimate. The range of fading rate parameters (g_{tc = 2days}-values: 1.1–3.5%/decade; ρ′: 0.76 × 10^–6 – 2.34 × 10^–6) is observed as normal. Apart from Huntley and Kars method (HKM) (Huntley 2006; Kars et al. 2008), the fading correction was performed using Huntley and Lamothe (Huntley and Lamothe 2001) (HLM) method as well (Singh 2023). As L_n/T_n’s are not in the linear range for older samples (KBPL 07 - 01), fading corrected ages by HLM were younger (Morthekai et al. 2011; Morthekai et al. 2007) than that of HKM (Singh et. al. 2023). Due to this reason, we have used only HKM correction procedure for further calculation. Considering the above, the final luminescence ages (LA) might not have underestimated because D_e will only overestimate. However, KBPL 01, KBPL 02 and the three intermediate samples might be older than what we presented (Table 2) because their L_n/T_n’s were meeting the dose response curve/growth curve after characteristic dose (D₀) as shown in Figure 2b. The dose recovery test (measured/given dose = 59.6 ± 0.5 Gy/60 Gy) performed on KBPL 03 after laboratory bleaching, indicates that the error in measured D_e is < 1%.

Regarding overestimation of dose rate, the NUSSY standard that was used to calibrate our high pure Ge gamma spectrometer was measured 29 times and the uncertainty in the reproducibility was within 2% for the concentration of K (0.96 ± 0.02%) and Th (7.3 ± 0.4 ppm) whereas U was underestimated by 20 ± 5% (of 2.7 ± 0.1 ppm). The values within parentheses of the earlier sentence were the concentration of K, Th and U in NUSSY (Preusser and Kasper 2001). Considering the fact that a significant contribution of U to the dose rate, the dose rate might have been underestimated only, but not otherwise. The range of measured a-values (0.09–0.11) is well within the reported values (Kreutzer et al. 2014). Hence, the LAs might have been overestimated but never underestimated.

Radioactive disequilibrium can also alter estimates of the dose rate. Due to the geochemical nature and the relatively shorter half-life of the radioactive daughters, the Th series is little prone for disequilibrium compared to that of U series (Degering and Degering 2020; Olley et al. 1996; Preusser et al. 2023). By comparing the radioactivity of ²²⁶Ra, ²¹⁴Pb and ²¹⁴Bi, the radioactive disequilibrium in U series was tested. A similar count rate (cts.m^–1.mg^–1) of ²²⁶Ra (186 keV), ²¹⁴Pb (295.2 keV) and ²¹⁴Bi (609.3 keV, 1120.3 keV, and 1700 keV) suggest that there is no disequilibrium in the U-series (Figure S3a). Another observation is the lesser than unity ratio of Th/U in the lake profile (Figure S3b) which is not expected in a crystalline rock terrain such as HHC zone. This may suggest either a decrease in Th or an increase in U concentration during burial which will result in erroneous luminescence ages. However, the observation of normal values of Th/U ratio in the moraines of this region (Ali et al. 2022) suggest the Th/U ratio was not modified during burial time and the Th/U ratio remained same since deposition. The increase in U concentration may be attributed to the prevailing oxidizing conditions, except during the intermediate samples’ period as inferred below (section Provenance of OC_petro).

This is to comment on the worst situation that is regarding accidental exposure of samples to the light either in the field or in the laboratory. The samples were processed in two sets namely, Set 1: KBPL 1, KBPL 2, KBPL 5, KBPL 6, KBPL 11, KBPL 12 and Set 2: KBPL 3, KBPL 4, KBPL 7, KBPL 8, KBPL 9, KBPL 10). Hence, it is least likely that a constant amount of white light exposure in the laboratory for the two sets of samples which would reduce the luminescence to a systematic lower value (lead to systematic underestimation of LAs).

The sample at 340 cm underestimate due, probably, to a higher authigenic carbonate (CaCO₃) as suggested by correlation between loss on ignition vs uranium concentration being r = 0.53 (p-value = 0.08, n = 12). A three thousand years offset in RA by linear extrapolation of calibrated RAs to modern surface sediment ( ± 1200 years) (Figure S4a) and linear regression of RAs and LA_measured ( ± 1000 years) (Figure S4b) was also observed (An et al. 2018; Hou et al. 2012; Jena et al. 2022). These analyses suggest the LA_measured are reliable and hence the observed discrepancy is real. Thus, a systematic overestimation of RAs compared to LAs is unambiguously observed.

Is OC _petro the cause for overestimated radiocarbon ages?

Considering the size (0.6 km²) and depth (∼4 m) of this relict lake, which remain frozen for not more than two months in the winter, the freshwater reservoir effect (FRE) may be ruled out as an explanation for the overestimation of RAs because CO₂ in lake waters might be in equilibrium with atmosphere (Carrizo et al. 2019; Doran et al. 1999; Hou et al. 2012; Lockot et al. 2016). Another potential cause for RA overestimation is the mixing of ancient organic carbon (OC_ancient) with lake sediments that dilute the radiocarbon concentration of the bSOM. The OC_ancient represents organic carbon particles that are older than the contemporary vegetation in and around the lake (biospheric OC; OC_bio) which includes OC_pre-aged and OC_petro. The age of OC_petro is > 60 kyr and that of OC_pre-aged is less than 60 kyr and older than contemporary vegetation.

The absolute concentration of C₁₀-C₄₀ for the plant, modern soil and pit sediments are shown in Figures S5 a, b. While all the samples have contributions from OC_ancient, the intermediate samples received relatively larger proportion of OC_ancient as characterized collectively by lower values of n-alkane derived parameters CPI_long and H/L chain ratios (Figures 6a and 7a, Data S4). Characterized by CPI_long, H/L, OEP and to some extent AOM, there are apparently two clusters among 12 samples (Figure 7a, Figure S6). The top 5 samples and the bottom-most sample are placed in one cluster, and their RAs are overestimated by 3.0 ± 0.6 kyr (RA_offset), and the remaining samples (excluding the intermediate samples) are in the second cluster with an offset of 4.6 ± 1.1 kyr. Relatively higher δ¹³C _n-alk values (by 0.7–2‰) of the C₂₉ homologue, and contrasting values (varying from 0.1 to 0.3‰) in modern plants, together, support our conclusion of higher contribution of relatively more refractory OC_ancient to the intermediate samples (Diefendorf and Freimuth 2017; Jambrina-Enríquez et al. 2018) (Figure 6b). Also, the δ¹³C_bulk values (–24 to 27.5‰) of the intermediate samples (Figure 7a) are similar to the crystalline bedrock of this region (–25.9 ± 1.7‰) (Menges et al. 2020). The range of R1 ratio (intensity at 1350 cm^–1 to that of 1580 cm^–1) from 0.024 – 0.99 calculated from the Raman spectra suggest that all the samples have moderate to highly graphitized OC_ancient (Beyssac et al. 2004; Beyssac et al. 2003) (Figure 6c). The presence of OC_ancient in ordered graphite form (French 1964) (2θ = 26.3°–26.6°) was observed in the XRD spectra in two intermediate samples (KBPR 05 and KBPR 07) and KBPR 10 with contribution of 6.2%, 13.9% and 12.9% of respectively (Figure S7a,b, Data S3).This suggests that some of OC_petro is moderate to highly graphitized. The available data is insufficient to conclude that all the OC_petro is in graphite form.

Figure 6. (a) Variation in carbon preference index (CPI_long) calculated from the relative concentration of long chain n-alkanes, (b) compound specific stable carbon isotope of n-alkanes (for C₂₉ and C₃₁ homologues) with the profile depth. Compound specific stable carbon isotope values measured from the modern plant samples (n = 4) and surface soil (n = 2) are also given for comparison. (c) The Raman spectra measured using 780 nm excitation on the selective opaque spots on the samples. The organization of organic matter index (R1 ratio) were calculated using graphite band (G) and D1 band. Two artificially made graphite samples were used for comparison. It is to be noted that the larger Raman peak doesn’t mean that particular sample has more concentration of graphite, it only says the measured graphite was of larger size.

Figure 7. Depth variation of selected (a) biotic variables and (b) abiotic variables. All the measured biotic and abiotic variables are given in Data S5. Palynofacies component that show significant variation in the profile is the amorphous organic matter (AOM; 4, 43, 25.8 small, large and average value). Variation in isotopic value for whole pit (δ¹³C; –26.8‰, –22.8‰, –24.6‰) and TOC (0.09 %, 1.1 %, 0.33%).

Provenance of OC _petro

The marked decrease (by ∼81.5%) of long chain n-alkanes (> C₂₄) concentration, relatively larger presence of unresolved complex mixture (UCM; C₁₆-C₃₀), and the higher amorphous organic matter in the samples after 6 kyr collectively indicate post-depositional microbial alteration of bSOM (Feakins et al. 2007; Freimuth et al. 2019). Lower δ¹³C_bulk values of the intermediate samples also support the degraded nature of OM (Krull and Retallack 2000; Sarangi et al. 2022). However, the bottom most sample (KBPL 01) is characterized by dominance of long chain n-alkanes with CPI_long higher than that of the modern plants and soil (Figure 6a), and suggests on the contrary, a non-degraded OM. Hence, our analysis suggests that OM was transported to the lake as already-degraded material, as opposed to having been degraded in-situ.

On the distinctly different nature of the intermediate samples compared to the rest (Figure S8), two inferences can be drawn: 1) enhanced hydrological conditions in the catchment either in terms of rain water or melt water (via lower C/N ratio, more negative δ¹³C_bulk, lower U, higher P_aq values, lower TAR), and 2) rapid transportation of more weathered but less recycled sediments (via higher CIA, higher ICV, lower plagioclase and higher illite) in the lake (Figure 7b, Figure S9, Data S5, cols. 23, 24, 21, 22). The anti-correlation between plagioclase and illite (r = –0.64, p-value = 0.03) indicates illite formation within plagioclase during weathering (CIA vs illite; r = 0.59, p-value = 0.04) specific to semi-arid regions (Bétard et al. 2009) (Figure 7b, Figure S9). Generally, a larger CIA implies a lower ICV i.e., more recycled and mature sediment must be highly weathered, but a positive correlation is observed in this study (r = 0.84, p-value = 0.0006). This suggests that the weathering happened before these materials were brought into the lake. The immature and less recycled sediments are also reflected in the lower sensitivity(Stokes et al. 2001) of feldspar grains during this period (1.2 ± 2.0 < 8.5 ± 15.6, × 10⁶ cts. Gy^–1.mg^–1).

We contend that the sediments that are highly weathered (and less recycled) are older moraine materials that were left behind by both the retreating glaciers after the Younger Dryas (YD, 11.7 kyr, MOHIT 1K) and two local advances of glacier at around 8 kyr (MOHIT 1H) and 6 kyr (MOHIT 1F) (Ali et al. 2022; Murari et al. 2014). During the advancing of glaciers, the OC_ancient was unearthed and made available for erosion when melt water was available (Zhou et al. 2016) (Figure 8). These materials were then transported in higher energy conditions (EM 4 & EM 5) (Figure 4b) at a higher erosion rate (higher Ti/K; Data S5, col. 35) to the lake during 10–6 kyr when Indian summer monsoon (ISM) was observed to be moderately intensified (Ali et al. 2020). As per the litholog, this period comprises a sand-gritty layer with coarser pebbles implying deposition under more energetic runoff conditions than before or after (Figure 3a). The concomitant overestimation in the LA_measured (average of 20 ± 5 kyr compared to actual modeled ages) can be explained by poor bleaching which has already been inferred from the presence of EM4 (Figure 4) and unexpectedly large luminescence/dose in the intermediate samples (Figure 2a).

Figure 8. The strength of Indian summer monsoon (ISM) which contribute 80 % of the total precipitation in Central Himalaya, was calculated for the last 19 kyr by synthesizing the inferences drawn by 27 research publications that are available for this region. The weighted average palaeoclimate index (WApCI), thus, indicate the strength of the ISM (line). The data points are the RA, calibrated radiocarbon ages. The overestimation of RAs of intermediate samples are correlating with the intensification of ISM.

A negative correlation was observed between TOC and RA_offset (r = –0.44, p-value = 0.23, n = 9; Figure 9a), excluding the intermediate samples. This suggests a relatively larger proportionate contribution of OC_petro to TOC that control the degree of overestimation of RAs. Assuming a constant proportion, with time, of OC_petro was buried into the lake along with OC_bio and a minor contribution from OC_pre-aged (Supplementary Materials 1, Figure S10a,b), the contribution of OC_petro to the lake sediment was calculated to be 0.064 ± 0.032% which was comparable to the published values in that region. This calculation yielded a similar value of δ¹³C_bio from this calculation (–23.6 ± 0.1 ‰) and measurements (–28.5 ± 1.7‰) using 29 modern plants in the study area after accounting for ∼4‰ enrichment during the transfer from plant to profile samples (Staddon 2004). The percentage of OC_petro to the lake sediment was calculated and found the RA_offset was positively correlated to it (r = 0.74, p-value = 0.22, n = 9; Figure 9b). It is clear that the sediment with lower TOC most probably have relatively larger OC_petro which resulted in overestimation of RAs.

Figure 9. Correlation between RA_offset with (a) TOC (%) and (b) percentage of OC_petro to TOC.

Implications

There are two implications of this study. The first one is on the radiocarbon dating of Quaternary sediments, which has regional significance. The age offset of the Kedarnath (30°47′24″N, 79°03′36″E) sediments is 0.9 kyr (Srivastava et al. 2017), while for the Benital Lake (30°09′38″N, 79°14′45″E), it is –0.9 kyr (Bhushan et al. 2018), using linear extrapolation. A 2.4 m core retrieved from the Chandra Tal (32°29′43″N, 77°36′48″E) also suffer an age offset of 1.2 kyr when the top three RA are used (Kumar et al. 2020). An inconsistent RA_offset among the sites in the HHC region suggest that RA_offset is variable. In the Trans Himalaya the lithology is dominated by carbonate rocks, and RA_offset of 4 ± 3 kyr (during 11–16 kyr) and 7.5 ± 4.0 kyr (during 11–18 kyr) have been estimated from the Burfu and Garbayang Lakes respectively (Beukema et al. 2011; Juyal et al. 2004). Apart from using compound-specific RA (Berg et al. 2020), age offset diagnoses can be done, (1) by employing another dating method as luminescence dating technique is used in this study (Doran et al. 1999; Juyal et al. 2009) and (2) by qualitatively assessing the reliability of RAs by comparing with the OC content (Strunk et al. 2020).

The second implication concerns the presence and the fate of OC_petro in the lake sediments from this region. Based on the conservative calculations (Supplementary Materials 1) using the contribution of OC_petro to the lake sediment (0.064 ± 0.032%), the average burial flux of OC_petro was calculated to be 160 ± 40 kg OC_petro.yr^–1 (Figure S10c) which is 33% of what is observed in High Himalayan upstream segment of Narayani River (Märki et al. 2021). This calculated burial flux of OC_petro to this lake is < 1% of the average global value per lake (31000 kg OC_petro.yr^–1.lake^–1; for OC_petro = 0.07 %) (Mendonça et al. 2017). The burial flux of OC_petro to this lake is less than one millionth of the flux of OC_petro to the headwaters of Himalayan rivers (1.15 × 10⁹ kg OC_petro.yr^–1; calculated from the 70 % of oxidized OC_petro during fluvial transport; Galy et al. 2008). Although the studied lake buried a meager amount of OC_petro in it, considering there are many such small lakes (∼1353 in number) (Pratima et al. 2021) in this region, the climate-sensitive OC_petro (Eglinton et al. 2021; Hein et al. 2020) that is stored in Central Himalayan freshwater lakes/ponds for a certain period (until it is filled) should be quantified.

Present date, there is no space for storage of OC_petro in the studied lake sediments because it is relict now, so the OC_petro is being transported to the Bengal fan via Kunti Banar River, a tributary to Ganges. Thus transported OC_petro is being either stored in marine sediments or oxidized during transportation and entrains the carbon cycle contributing to global warming. Although the OC_petro can withstand multiple cycles of erosion (Sparkes et al. 2020), it is reportedly vulnerable to oxidization during fluvial transport in the Ganges (Galy et al. 2008) as not all the OC_ancient was fully graphitized. There is about 36% of bedrock OC_petro oxidized in the catchment of Narayani River (Märki et al. 2021). Along with OC_petro, the lake also apparently received significant contribution from OC_pre-aged which is more prone for oxidation.

Conclusion

Twelve samples from a 4 m sediment deposit of a relict lake in the High Himalayan Cyrstalline zone were co-dated using radiocarbon and luminescence methods. All the AMS measured radiocarbon ages were observed to be overestimating compared to the fading corrected luminescence ages based on polymineral fine grains. The n-alkane data of the sediment profile, modern sediment and modern plant samples, and their stable carbon isotopic compositions suggested that there was a possible mixing of ancient organic carbon in the lake sediment. Raman spectra and X-ray diffraction spectra suggested the presence of graphite in the sediment. The presence of ancient organic carbon and some in graphite form collectively suggest that diluted radiocarbon concentration in the sediment and hence the overestimation of radiocarbon ages that is negatively correlated with total organic carbon in the sediment. It is also observed that the old sediment with low organic carbon content were excavated during glacial advance and its subsequent transport to the lake during deglacial time with relatively intensified monsoonal rain. This old sediment with relatively large proportion of OC_petro overestimated the radiocarbon ages and luminescence ages by poor bleaching due to rapid transport. This result warns us to be alert on the radiocarbon ages in this region when the organic content in the sediment is low. A more detailed study to constrain the burial flux of OC_ancient to the lakes of this region is also needed.