Subsurface CO₂ dynamics in a temperate karst system reveal complex seasonal and spatial variations

Understanding the carbon cycle of the terrestrial critical zone, extending from the tree canopy to the aquifer, is crucial for accurate quantification of its total carbon storage and for modelling terrestrial carbon stock responses to climate change. Caves and their catchments offer a natural framework to sample and analyse carbon in unsaturated zone reservoirs across various spatial and temporal scales. In this study, we analyse the concentration, stable carbon isotopic ratio (δ¹³C), and radiocarbon (¹⁴C) compositions of CO₂ from the atmosphere, boreholes (0.5 to 5 m depth), and cave sampled every 2 months over 2 years at Milandre cave in northern Switzerland. High concentrations of up to 35 000 ppmV CO₂ are measured in the boreholes. The δ¹³C values of CO₂ in the boreholes reflect the δ¹³C of C₃ plants (∼ −26 ‰), which dominate the catchment ecosystem. Shallow meadow boreholes host older CO₂ in winter and modern CO₂ in summer, while forest ecosystems consistently export modern CO₂ (F¹⁴C = ∼ 1) to the unsaturated zone. Cave CO₂ concentrations exceed atmospheric levels and are diluted by temperature-driven seasonal ventilation. Keeling plot intercepts indicate that the cave CO₂, which mixes with atmospheric CO₂, is younger in summer (F¹⁴C = 0.94) and older in winter (F¹⁴C = 0.88), with a δ¹³C consistent with the C₃-plant-dominated catchment. Mixing models utilizing drip water dissolved inorganic carbon ¹⁴C suggest that varying carbonate dissolution and degassing dynamics do not explain the F¹⁴C variation and concurrent δ¹³C stability of the mixing endmember. Rather, contributions from aged carbon reservoirs in the deeper unsaturated zone are likely. This study provides valuable insights into CO₂ source dynamics and cycling within the karstic critical zone, highlighting the impact of seasonal variations and ecological factors on downward carbon export from terrestrial ecosystems.

1 Introduction

Within the context of current and future climate change, investigations into the response of critical zone carbon pools to rising temperatures and changing hydroclimatic conditions are crucial (Brantley et al., 2007). Recent studies have specifically highlighted the unsaturated zone as a potentially important but poorly constrained reservoir of gaseous CO₂ (Mattey et al., 2016; Noronha et al., 2015; Keller, 2019; Stewart et al., 2022). Estimates suggest that between 2 and 53 PgC could be present in the form of CO₂ in the unsaturated zone globally (Baldini et al., 2018). This carbon is particularly vulnerable to changes in the water table level, whereby rises may easily result in the rapid release of CO₂ into the atmosphere (Baldini et al., 2018). Despite its importance, comprehensive assessments integrating spatial and temporal variability of shallow subsurface CO₂ remain scarce.

Understanding unsaturated zone CO₂ dynamics is complicated by the different sources of carbon contributing to the subsurface reservoir. The CO₂ present in the unsaturated zone is often referred to as ground air and was first defined as CO₂ produced by microbial oxidation of organic material which was transported from the surface (Atkinson, 1977). However, ground air can also refer more generally to high CO₂ concentrations in the subsurface, without linking its presence to a particular source. Early evidence for the presence of a ground air reservoir was the observation of high CO₂ concentrations found in deep cave passages with poor connection to the catchment surface, suggesting an endogenous source of CO₂ (McDonough et al., 2016). The age profile of this high-concentration CO₂ reservoir, however, varies by site, suggesting significant variability in its source and processing. CO₂ reservoirs of a lower F¹⁴C, corresponding to an older age, have been reported (Breecker et al., 2012; Noronha et al., 2015; Mattey et al., 2016; Bergel et al., 2017), as well as modern subsurface CO₂ pools, likely derived from ecosystem respiration of recently fixed carbon (Campeau et al., 2019; Tune et al., 2020). Consequently, the primary sources of subsurface CO₂ generally considered are (1) the outside atmosphere (Kukuljan et al., 2021), (2) catchment soil and vegetation respiration (Breecker et al., 2012; Li et al., 2024), (3) microbial respiration in the unsaturated zone (Mattey et al., 2016; Ravn et al., 2020), (4) carbon dissolved from the carbonate bedrock (Milanolo and Gabrovšek, 2015), and (5) volcanic and metamorphic hydrothermal input (Chiodini et al., 2008; Girault et al., 2018).

Carbon isotope analysis of CO₂ can be used to differentiate the contributing sources to the subsurface gas mixture. Specifically, the δ¹³C of CO₂ in the unsaturated zone can provide information about the influence of biological fractionation by overlying catchment vegetation and microbial respiration and the influence of carbonate dissolution and subsequent degassing (Breecker, 2017; McDonough et al., 2016). In addition, due to the contrast between the ¹⁴C content of biospheric carbon, which ranges from bomb-peak-enriched to > 10 000 years old, and carbonate rock, which is typically radiocarbon-dead, ¹⁴C measurements of subsurface CO₂ can provide information for both source apportionment and carbon turnover rates in the terrestrial subsurface (Noronha et al., 2015).

Though there has been extensive work investigating unsaturated zone CO₂ reservoirs, questions remain about the concentration, composition, and sources of ground air on temporal and spatial scales. This study seeks to address this gap by presenting a detailed assessment of CO₂ concentrations and isotopic compositions over 2 years within the Milandre cave karst system in Switzerland. We aim to (1) determine the relationships between pCO₂, δ¹³CO₂, and ¹⁴CO₂ in Milandre cave and its catchment; (2) assess the effect of seasonality on CO₂ concentration and isotopic characteristics; and (3) gain insights into the provenance of CO₂ and how these evolve over time.

Figure 1 Map of the main passages of the Milandre cave network. The locations of the atmospheric sampling (yellow), the soil (green) and unsaturated zone boreholes (blue), and the cross-trip sampling (red) are annotated. Cave survey modified from Gigon and Wenger (1986). Base map from ©swisstopo (Swisstopo, 2024).

2 Site overview

Milandre cave (47.4852° N, 7.0161° E, 373 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) is located in the municipality of Boncourt in the Jura canton, NW Switzerland (Fig. 1). The Milandre karst structure formed within the Jura Mountains, a sub-alpine mountain range which lies laterally in a northwest–southeast direction. The cave is located within the external plateau unit of the Jura Mountains, which consists of thin subhorizontal Mesozoic limestone units that have experienced deformation which produced imbricates and strike-slip faults (Sommaruga, 2011). More specifically, the Milandre karst system formed within the Late Jurassic (Middle Oxfordian) St-Ursanne Formation, which overlies Early to Middle Oxfordian marls (Jeannin, 1998). The unsaturated zone ranges between 40 and 80 m depth. The active karst level developed in the epiphreatic zone, with the underlying phreatic zone extending for ∼ 20 m to reach the low-permeability Lower Oxfordian marls (Liesberg Mb) (Perrin et al., 2003). The cave is formed over 2.7 km, with 10.5 km of galleries forming along N–S-oriented fractures. The Milandre river flows for 4.5 km in a northerly direction in the lower part of the cave and exits via the perennial Saivu spring where it joins the Allaine River and the Bame temporary spring. Clays and silts, deposited during hydraulic backflooding, are particularly present in the downstream part of Milandre cave where they contribute to the sediment load during mild flood events. In contrast, during medium to intense flow, allochthonous sediment transport is often associated with higher microbial concentrations (Vuilleumier et al., 2021). Quantifying the contribution of microbial respiration to the cave carbon cycle was beyond the scope of this study and should be addressed by future research. The Milandre river drains a catchment area of approximately 13 km² (Perrin et al., 2003). There is a series of dolines in the catchment area which form depressions ∼ 10 to 20 m deep in the landscape. Notably, one doline is located in the meadow partially covering the upstream section of Milandre cave (Fig. 1).

Over the past 40 years, land use in the catchment area has been dominated by farmland (37 %), forest (36 %), and meadows (12 %) as estimated by aerial photo analysis (Jeannin et al., 2016). The crops grown primarily include maize and tobacco, two C₄ plants. Long-term land use analysis shows that farmland has decreased slightly in this area and special infrastructures have increased due to the construction of a motorway that overlies part of the cave. The soils of the area are influenced by the ecosystem that is overlying it. Forest soils are shallow Leptosols (< 10 cm deep) that are rich in fragments of the carbonate bedrock. The shallow forest soils transition into deeper (up to ∼ 80 cm) organic-rich Histosols in the meadows, particularly in the dolines. The roots from the forest vegetation likely extend into the epikarst and unsaturated zone; however, no roots are visible within the cave itself.

Table 1 Summary of the boreholes sampled in Milandre cave catchment noting sample depth, number of sampling lines, classification in soil/unsaturated zone, and ecosystem cover type. The IDs are the names given to the boreholes for this study and are concurrent with the supplementary data provided. SISKA (Swiss Institute for Speleology and Karst Studies) IDs refer to the original names of the boreholes at the time of installation.

The Jura region experiences a marine west coast, warm summer climate (Cfb classification) (Kottek et al., 2006). Daily temperature measurements from the Fahy MeteoSwiss weather station located ∼ 9.5 km SW from Milandre cave show a mean temperature of 9.4 °C (1991 to 2020). Temperature seasonality is strong, with monthly mean temperatures fluctuating between a minimum of −1.2 °C in January and a maximum of 18.2 °C in July from 1991 to 2020 (MeteoSwiss, 2024). In contrast, temperatures within the cave remain almost constant year-round and vary between 10.3 and 11.0 °C (Affolter et al., 2020). The temperature difference between the outside and in-cave temperature drives dynamic ventilation on seasonal scales (Garagnon et al., 2022). Regional meteoric precipitation shows a mean of 1046 mm yr⁻¹. Monthly precipitation data indicate that precipitation in this area is well dispersed throughout the year (MeteoSwiss, 2024). Monitoring of tritium (³H) in stalagmite drip water was used to estimate the residence time of seepage water at ∼ 5.5 to 6.6 years (Affolter et al., 2020).

3 Materials and methods

3.1 Gas sampling setup

All gas samples were taken every 2 months from December 2021 to January 2024. The samples were collected in 5 L sampling bags (Cali-5-Bond, Calibrated Instruments, USA) using a handheld pump. The gas was dried through a glass tube filled with granular magnesium perchlorate (∼ 83 % purity, Supelco, Germany) to reduce the effects of humidity, which may affect the isotopic composition of the sample. The sampling setup and procedure were identical for all gas samples. The CO₂ concentration was monitored during line flushing using a nondispersive infrared (NDIR) CO₂ sensor (SCD30, Sensirion, Switzerland) to ensure accurate sampling. The magnesium perchlorate was exchanged before each sampling day. Prior to analysis, the samples were stored away from direct sunlight in cool indoor temperatures for a maximum of 6 weeks. The same bags were reused for each sampling location and were flushed with N₂ gas for at least 48 h before sampling to reduce cross-contamination risks.

Atmospheric samples were collected at a defined sampling site above the cave (Fig. 1). Before sampling, the line was manually flushed for 1 min. To reduce the risk of breath contamination, the sampling setup was attached to 3 m long Teflon tubes, and samples were always taken against the airflow direction.

The unsaturated zone air was sampled from six boreholes of varying depths between 0.5 and 5 m (Table 1). Spatially, the boreholes cover a large portion of the cave's upstream hydrological catchment and are overlain by contrasting vegetation cover (mixed deciduous forest or a grass meadow) (Fig. 1, Table 1). The gas was sampled in two lines from depths of 0.5 to 0.85 m (Shallow 1), three lines from 0.6 to 1.5 m (Shallow 2), and single lines from a depth of 5 m (Deep 1, Deep 2, Deep 3, Deep 4). Due to the nature of the installation, it was not possible to assign a specific depth to the multi-line boreholes. All boreholes were drilled in 2013 and are equipped with aluminium tubes compacted by layers of gravel, bentonite, and sand. The sampling line was flushed for 1 min once attached to the borehole, and then samples were taken into the Cali-5-Bond bags. Due to the nature of the borehole installation, we assume that the samples were not taken in steady-state conditions, where the production rate of CO₂ would equal the sampling rate.

Cave air was sampled in two locations within Milandre cave proximal to both entrances, Downstream (topographically lower) and Upstream (topographically higher), with the same procedure as for the atmospheric samples (Fig. 1). To gain insight into possible spatial variability within the cave, a trip through the cave across the main passage was carried out in September 2021, and 14 samples were collected (Fig. 1). The trip was taken against the direction of airflow to reduce the risk of breath contamination, from the Downstream to Upstream entrance.

3.2 Continuous O₂ and CO₂ concentration monitoring

Soil gas O₂ and CO₂ concentrations were measured continuously near Shallow 2 at 40 cm depth and recorded at 10 min intervals between April 2023 and June 2024 using a SCD-41 CO₂ sensor (± 50 ppmV + 5 % of reading) (Sensirion, Switzerland) and SGX-40X O₂ sensor (Amphenol SGX, Sensortech, Switzerland). Calibration was cross-checked with an independent, handheld multi-gas monitor (XAM 5600, Dräger, Germany) before installation.

3.3 CO₂ concentration and δ¹³C analysis

The concentration of CO₂ and the stable carbon isotopic composition (δ¹³C) in all gas samples was measured using a cavity ring-down spectrometer (CRDS) G2131-I isotopic CO₂ instrument (Picarro, USA) at ETH Zurich (CO₂ concentration = 0.2 ppm, δ¹³C = < 0.1 ‰ precision). Standard gases with known CO₂ concentrations of 399.6 and 2000 ppmV CO₂ in synthetic air were measured in addition to standard gases of known δ¹³CO₂ values of −27.8 ‰ and −2.8 ‰ VPDB for offline calibration. As there was no equivalent standard gas available for high CO₂ concentrations as in some samples, linearity of the concentration data produced by the CRDS is assumed. All standards and samples were measured by directly attaching the sample bags to the inlet of the CRDS that was previously fitted with a magnesium perchlorate dryer and measuring gas concentration and isotopic composition for 2 min after reaching steady state. We used the mean and standard deviation of this measurement interval for evaluation. δ¹³C refers to the ratio between the two stable carbon isotopes with respect to the Vienna Pee Dee Belemnite (VPDB) standard.

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left(\mathrm{‰}\right)=\left({\displaystyle \frac{{\left(\frac{{}^{\mathrm{13}}\mathrm{C}}{{}^{\mathrm{12}}\mathrm{C}}\right)}_{\mathrm{sample}}}{{\left(\frac{{}^{\mathrm{13}}\mathrm{C}}{{}^{\mathrm{12}}\mathrm{C}}\right)}_{\mathrm{standard}}}}-\mathrm{1}\right)\times \mathrm{1000}\end{array}$$

3.4 Radiocarbon analysis

The CO₂ from each gas sample was converted into approximately 1 mg of graphite using automatic graphitization equipment (AGE, IonPlus, Switzerland) (Wacker et al., 2010a) coupled with a custom-made carbon inlet system. In this process, each sampling bag was successively attached to the inlet system, which dried the sample through a fine-grained magnesium perchlorate water trap via a vacuum line. The volume of air required to sample was calculated from the CO₂ concentration of the sample, allowing the appropriate amount of CO₂ to be trapped within a stainless-steel coil that was cooled by submersion in a liquid N₂ bath. After trapping, the CO₂ was sublimated by submerging the coil in a room-temperature water bath and then adsorbed onto a zeolite trap using helium as a carrier gas between the coil and the trap. The CO₂ adsorbed onto the trap was then thermally desorbed and filled into the AGE reactor, where it was reduced to graphite with hydrogen over an iron catalyst (Wacker et al., 2010a). A modern standard of Oxalic Acid II gas (F¹⁴C = 1.3407, Oxa II, NIST SRM 4990C, HEKAL AMS Lab, Hungary) and a radiocarbon fossil reference CO₂ gas (F¹⁴C = 0, CO₂ ≥ 99.7 % Vol.abs, Carbagas, Switzerland) were also graphitized and used to calibrate the measurements. The resulting graphite was pressed into targets for radiocarbon analysis through AMS analysis with a MIni CArbon DAting System (MICADAS, Ionplus, Switzerland) (Synal et al., 2007; Szidat, 2020). The resultant radiocarbon data were corrected using the BATS (4.3) software (Wacker et al., 2010b). In this study the ¹⁴C content of samples will be discussed in F¹⁴C notation, which represents fraction modern (defined as the 1950 atmosphere) according to Reimer et al. (2004).

$$\begin{array}{}\text{(2)}& {\mathrm{F}}^{\mathrm{14}}\mathrm{C}={\displaystyle \frac{{\frac{{}^{\mathrm{14}}\mathrm{C}}{{}^{\mathrm{12}}\mathrm{C}}}_{\mathrm{sample}}}{{\frac{{}^{\mathrm{14}}\mathrm{C}}{{}^{\mathrm{12}}\mathrm{C}}}_{\mathrm{standard}}}}\end{array}$$

3.5 Dissolved inorganic carbon analysis

Drip water dissolved inorganic carbon (DIC) ¹⁴C was measured to constrain the DIC degassing endmember. Samples were measured over 1 year from December 2021 and December 2022 at nine drip sites focused on the actively dripping galleries near Upstream and Downstream (Fig. 1). Using a 5 mL syringe, 1 mL of drip water was collected directly from the soda straw stalactites on the cave roof to reduce fractionation effects associated with degassing. The water was injected directly into pre-cleaned Exetainer^® (Round Bottom, Borosilicate, 938W, Labco Limited, UK) vials which had been flushed with helium gas and pre-spiked with 150 µL of 85 % H₃PO₄ (Suprapur^®, 85 %, Merck KGaA, Germany). The ¹⁴C content of the DIC was measured by sampling the vial headspace using the Carbonate Handling System (CHS, Ionplus, Switzerland) and transferring the gaseous sample to the AMS using a gas handling system. Here, the CO₂ was introduced to the gas ion source through the CHS at a flow of 50 mL min⁻¹. The ¹⁴C content for each sample was measured for ∼ 60 cycles. Standard materials with a carbonate matrix of similar composition to the DIC samples were used as standards (IAEA C1 (F¹⁴C = 0) and C2 (F¹⁴C = 0.411), sodium bicarbonate (NaHCO₃, Sigma Aldrich, USA) and potassium bicarbonate (KHCO₃, Sigma Aldrich, USA) (both F¹⁴C = 0)).

Figure 2 δ¹³C and F¹⁴C of all gas samples over the monitoring period. The theoretical carbon reservoirs present in the system are shown with the atmosphere (blue), C₃ plants (yellow), C₄ plants (green), soils of various ages and compositions (pink arrows), and carbonate rock (grey).

3.6 Data analysis and statistics

All data and statistical analyses and all graphs were generated using the Python 3 programming language (Van Rossum and Drake, 2009).

Relationships between CO₂ concentration, δ¹³C, F¹⁴C for each sample type, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather stations were explored using Spearman's rank correlation coefficient analysis (Kokoska and Zwillinger, 2000). A significantly correlated relationship between two variables is defined by p < 0.05, with the correlation coefficient denoted by “rho” and the number of samples involved in the analysis as n.

4 Results

4.1 Atmospheric CO₂

The atmospheric samples had a CO₂ concentration ranging from ∼ 380 ppmV (August 2023) to ∼ 485 ppmV (December 2022) with a mean of ∼ 440 ppmV. The δ¹³C ranged from −12.5 ‰ (February 2023) to −7.6 ‰ (August 2023) and had a mean of −10.1 ‰ (n = 17). The CO₂ sampled was typically modern, though it had some fluctuation from F¹⁴C 0.98 (December 2022) to 1.02 (June 2022) and a mean of F¹⁴C 1.0 (Fig. 2).

Figure 3 Maximum and minimum (a) CO₂ concentration, (b) δ¹³C, and (c) F¹⁴C of the three soil borehole samples from Shallow 1 in the forest (green) and the two sampling lines in Shallow 2 in the meadow (brown). Summer (yellow) and winter (blue) seasons are highlighted by the coloured background. Concentrations are shown as ranges between maxima and minima as they aggregate over three sampling lines in the soil per location (Table 1).

4.2 Soil zone CO₂

Soil boreholes have higher CO₂ concentrations and lower δ¹³C compared to the atmosphere (Fig. 2). Concentrations in boreholes that sample the soil zone (Shallow 1 and Shallow 2) (Fig. 3a) vary significantly depending on sampling depth, cover type, and season. The Shallow 1 boreholes (forest) show seasonal high pCO₂ in the summer months from June to October, which steadily declines during winter until it begins to rise again after April. In these boreholes CO₂ concentrations vary from a maximum of ∼ 18 000 ppmV in August 2022 to a minimum in January 2024 of ∼ 1200 ppmV, with a mean of ∼ 10 700 ppmV (n = 21). The CO₂ concentrations in Shallow 2 (meadow) boreholes range from ∼ 27 000 ppmV in June 2023 to ∼ 4800 ppmV in February 2023 with a mean of ∼ 9900 ppmV (n = 30).

The δ¹³C varies by 7.0 ‰ in shallow forest boreholes over the sampling period, ranging from a minimum of −26.7 ‰ in August 2023 to a maximum of −19.5 ‰ in January 2024 (Fig. 3b). There is ∼ 4.0 ‰ isotopic variation in the shallow meadow boreholes with a minimum of −29.0 ‰ in October 2022 and a maximum of −24.4 ‰ in August 2022. The pCO₂ and δ¹³C are negatively correlated in both shallow forest and shallow meadow boreholes (forest: rho = −0.47, p = 0.02, n = 25. meadow: rho = −0.49, p = 0.01, n = 30), and δ¹³C is on average lower in shallow meadow than shallow forest boreholes.

F¹⁴C in Shallow 1 boreholes (forest) shows little annual variability and fluctuates slightly around a modern value (max = 1.03 in April 2022, min = 1.0 in June 2023) (Fig. 3c). Conversely, the Shallow 2 (meadow) boreholes show distinct annual variation in F¹⁴C, which steadily decreases from modern (1.0) in August 2022 and 2023 to a low in April 2023 and a minimum in January 2024 (0.87). The samples from the shallow forest boreholes have a positive correlation between F¹⁴C and δ¹³C (rho = 0.48, p = 0.02, n = 25). In the shallow meadow borehole samples, we observe a positive correlation between F¹⁴C and MMT (rho = 0.67, p = 0.001, n = 30). All other tested parameters show no significant correlations (see Appendix A).

Figure 4 Continuously measured CO₂ (pink) and O₂ (blue) concentrations from the soil zone near the Shallow 2 boreholes. CO₂ concentrations from Shallow 2 are also plotted (brown).

The continuously measured O₂ and CO₂ concentrations near Shallow 2 in the meadow show an inverse relationship (Fig. 4). The concentrations are relatively stable throughout most of the year with higher O₂ and lower CO₂ concentrations. This pattern is disturbed by a sharp decrease in O₂ and an increase in CO₂ observed around May/June in both 2023 and 2024. These different conditions last around 1 month before returning to the typical concentrations. The O₂ concentration ranges from a maximum of 20.2 % in August 2023 to a minimum of 17.8 % during a brief interval in May 2023. Conversely, the CO₂ concentration peaks at ∼ 17 000 ppmV in May 2023 and is the lowest at ∼ 2500 ppmV in August 2023. The Shallow 2 boreholes also show peaks in CO₂ concentration in May to June 2023 (Fig. 4).

Figure 5 (a) CO₂ concentration, (b) δ¹³C and (c) F¹⁴C of unsaturated zone borehole samples from forest boreholes Deep 1 and Deep 2 (green), and meadow boreholes Deep 3 and Deep 4 (brown). Summer (yellow) and winter (blue) seasons are highlighted by the coloured background.

Figure 6 (a) CO₂ concentration as well as (b) δ¹³C and (c) F¹⁴C of cave air samples from the Downstream (red) and Upstream (orange) sites. Summer (yellow) and winter (blue) seasons are highlighted by the coloured background. The approximate timing of the temperature-driven ventilation direction changes is shown by the annotations.

4.3 Unsaturated zone CO₂

Unsaturated zone CO₂ concentrations varied widely during the sampling period but are generally very high (> 10 000 ppmV). On average, the forest boreholes Deep 1 and Deep 2 have higher concentrations than the meadow boreholes Deep 3 and Deep 4. Concentrations are generally highest in the Deep 1 borehole (forest) (Fig. 5a) with a maximum of ∼ 37 000 ppmV in June 2023 and a similarly high concentration peak in June 2022. The lowest values (∼ 21 000 ppmV) at this borehole are recorded in the winter months. Deep 2 (forest) shows an opposite trend with its highest CO₂ concentrations in autumn, winter, and spring (max = ∼ 29 000 ppmV in December 2022). The lowest CO₂ concentrations were measured during the summer months (min = ∼ 7900 ppmV in June 2023). There is comparatively less variability in the Deep 3 and Deep 4 boreholes (both in the meadow). In Deep 3, concentrations fluctuate from ∼ 6100 to ∼ 19 000 ppmV with no obvious seasonality. Deep 4 has higher concentrations over the winter period and lower during summer, peaking at ∼ 14 000 ppmV in October 2022 and dropping to an overall minimum of ∼ 3000 ppmV in June 2023.

The δ¹³C of the majority of unsaturated zone CO₂ samples is very stable around ∼ −26 ‰ (Fig. 5b). Similar to the shallow boreholes, samples from deep meadow boreholes (Deep 3 and 4) have on average slightly lower δ¹³C values. In the meadow boreholes Deep 3 and Deep 4, the δ¹³C and pCO₂ are negatively correlated (rho = −0.66, p = 0.001, n = 21) (see Appendix A).

The F¹⁴C of unsaturated zone CO₂ samples is relatively stable throughout the sampling period, ranging from 1.01 (Deep 1, June 2022) to 0.98 (Deep 4, October 2022) (Fig. 5c). A possible seasonal pattern can be observed in the meadow boreholes (Deep 3 and Deep 4) with lower F¹⁴C values in the winter and spring months and higher F¹⁴C during summer, similar to the observation from the shallow meadow boreholes. No clear seasonal pattern is observed for Deep 1 and Deep 2. The F¹⁴C values of unsaturated zone CO₂ samples are negatively correlated with their δ¹³C value (rho = −0.54, p = 0.01, n = 21) and positively correlated with MMT (rho = 0.53, p = 0.01, n = 21). All other tested parameters show no significant correlations.

4.4 Cave CO₂

The CO₂ concentrations at the Downstream and Upstream sampling points in the cave are inversely correlated (Fig. 6a). Downstream shows a higher pCO₂ from June to October in 2022 and 2023, with a maximum of ∼ 14 000 ppmV in August 2023. The lowest concentrations are seen between December and April, reaching nearly atmospheric concentrations with a minimum of ∼ 420 ppmV in December 2021. Conversely, the highest CO₂ concentrations at the Upstream site were measured during the winter months, with a maximum of ∼ 30 000 ppmV in December 2021. The lowest CO₂ concentrations at Upstream were observed during the summer months with a minimum of ∼ 2300 ppmV in August 2022. Overall, concentrations at the Upstream site reach higher maxima and do not approach atmospheric values at their minimum, as at the Downstream site. The isotopic composition of CO₂ is inversely related to its concentration, with higher CO₂ concentrations coinciding with lower F¹⁴C and δ¹³C values. Seasonal trends differ between sites; at the Downstream site δ¹³C and F¹⁴C values are higher in winter and lower in summer, whereas at the Upstream site, δ¹³C and F¹⁴C values show the opposite pattern (Fig. 6b and c). The δ¹³C in Downstream ranged between −25.9 ‰ in June 2023 and −13.6 ‰ in December 2021, and the F¹⁴C ranged between 0.91 in October 2023 and 0.98 in February 2023. At Upstream, the opposite relation with CO₂ concentration is well expressed for F¹⁴C but not as clear for δ¹³C, with minimal variability over the entire study period. The δ¹³C in Upstream fluctuates slightly around a mean of −24.8 ‰.

The F¹⁴C of CO₂ in Upstream is overall lower than that of Downstream, increasing over late spring 2022 and into summer, peaking during August 2022 at 0.95. Decreasing F¹⁴C is observed during winter and into spring with lows of ∼ 0.88 in February 2023. At the Downstream site, the pCO₂ and δ¹³C (rho = −0.9, p = 0.0001, n = 14), pCO₂ and F¹⁴C (rho = −0.83, p = 0, n = 14), and δ¹³C and MMT (rho = −0.79, p = 0.0001, n = 14) are negatively correlated. Moreover, the pCO₂ and MMT (rho = 0.53, p = 0.05, n = 14) and the δ¹³C and F¹⁴C (rho = 0.74, p = 0, n = 14) are positively correlated. At the Upstream site, pCO₂ and MMP are positively correlated (rho = 0.55, p = 0.04, n = 14), and pCO₂ and MMT are negatively correlated (rho = −0.54, p = 0.05, n = 14) (see Appendix A). All other tested parameters show no significant correlations.

Figure 7 Keeling plots for δ¹³C (a) and F¹⁴C (b) of cave (Upstream and Downstream sites) and atmospheric samples from all monitored months. The dashed regression lines denote the months with the maximum and minimum y intercepts for δ¹³C and F¹⁴C. (c) Keeling plot y intercepts for δ¹³C (blue) and F¹⁴C (orange) over time. Summer (yellow) and winter (blue) seasons are highlighted by the coloured background.

To constrain the isotopic composition of the endmembers contributing CO₂ to the cave air mixture, we use the Keeling plot approach (Fig. 7a and b). This approach assumes that cave air is a mixture of two main sources, the atmosphere (with known concentration and isotopic values) and a second source of a priori unknown composition. The y intercepts of the Keeling plots represent the isotopic composition of the contributing endmember which mixes with atmospheric air inside the cave (Keeling, 1961; Pataki et al., 2003) (Fig. 7a and b). We find that over time, the isotopic composition of the endmember varies for both δ¹³C and F¹⁴C (Fig. 7c). The endmember δ¹³C value varies by ∼ 2.0 ‰ from −26.8 ‰ in June 2023 to −25.0 ‰ in February 2022. The variation in the F¹⁴C is larger, ranging from 0.88 in February 2023 to 0.94 in June 2022, August 2022, October 2022, and June 2023 (ca. ∼ 0.06 F¹⁴C). Maxima in the F¹⁴C value derived for the endmember through the Keeling plot generally correspond to decreases in δ¹³C, with two F¹⁴C maxima occurring in June to October 2022 and in June to August 2023.

Figure 8 (a) CO₂ concentration, (b) δ¹³C, and (c) F¹⁴C of cave air samples from the different sites sampled during the trip through the cave. Sampling locations are shown in Fig. 1. The blue crosses and blue bar notate where samples were taken standing in the river. The direction of cave air ventilation is shown by the lower arrow. The approximate locations of the Upstream and Downstream sampling sites are annotated.

Figure 9 (a) Mixing model between f_soil and f_DIC for the measured range of F¹⁴C DIC values and F¹⁴C cave values of 0.94 (blue), 0.9 (orange), and 0.88 (green). The grey bar shows physically possible mixing ratios between 0 and 1. (b) Mixing model between f_soil and f_DIC using the extended DIC F¹⁴C range and cave air F¹⁴C values. (c) Cave air δ¹³C for varying DIC F¹⁴C and cave air F¹⁴C values with equilibrium fractionation at −9.54 ‰ (circles) and kinetic fractionation at −30 ‰ (stars). The grey bar shows the cave air δ¹³C range which fits the measured values between −25 ‰ and −27 ‰.

Samples taken along the main passage of the cave during the trip through the cave show increasing pCO₂ from the lower entrance (location 1, corresponding to the Downstream monitoring site, ∼ 6900 ppmV) to a plateau of ∼ 14 000 ppmV from sample point 4 to 9 towards the centre of the cave (Fig. 8a). The pCO₂ decreases sharply after sample point 9 until point 11 (∼ 6600 ppmV). The pCO₂ increases briefly at sample point 12 and then decreases again upon approach to the cave exit (Upstream). The δ¹³C values remain essentially constant at −25.5 ‰ through the initial segments of the cave, beginning to increase slightly past sample point 9 (Fig. 8b). The δ¹³C value is −25.8 ‰ at sample point 9 and progressively increases as pCO₂ levels drop between sampling points 9 and 11, followed by a strong increase toward the Upstream exit (−11 ‰). The F¹⁴C values are constant at ∼ 0.93 between sample points 1 and 5, followed by more variability and generally lower F¹⁴C values in the middle of the transect, coinciding with the river passage of the cave (minimum 0.89 at sampling point 11; Fig. 9c). Towards the Upstream exit, F¹⁴C begins to increase strongly (maximum 0.98 at sampling point 15; Fig. 8c). Sample 13 is excluded due to a ruptured bag during sampling.

5 Discussion

5.1 Sources and variability of shallow-depth CO₂

The observed fluctuations in CO₂ concentrations in the Shallow borehole gas samples vary with depth, vegetation cover, and season (Fig. 3), suggesting complex carbon cycling dynamics taking place in the soil zone and shallow fractured epikarst. Overall, samples from both forest and meadow locations have similar mean CO₂ concentrations. In the Shallow 1 forest boreholes, higher CO₂ concentrations occur during summer months (Fig. 3a), aligning with increased autotrophic (occurring in plant roots, leaves, and stems) and heterotrophic (microbial) respiration rates during warmer seasons. This results in higher CO₂ production compared to losses from soil efflux to the atmosphere and downward transport in gaseous or dissolved form. Due to our discrete sampling approach, we cannot discern whether higher production or decreased losses are responsible for higher CO₂ concentrations. The decline in CO₂ concentrations during winter suggests that reduced autotrophic and soil microbial respiration in response to lower temperatures is likely an important driver of the seasonal cyclicity. The lower CO₂ concentrations during winter in the forest boreholes (Shallow 1) compared to the meadow boreholes (Shallow 2) could be an indication of frost effects. Soil frost induced by colder temperatures can reduce heterotrophic respiration in soils. Experiments inducing deep frost in forested soils by removing snow cover caused a greater decrease in heterotrophic respiration during winter months compared to snow-covered control soils (Muhr et al., 2009). Furthermore, frost heave can physically damage the extraradical hyphae on plant roots (Lekberg and Koide, 2008). Due to the shielding of the tree canopy, the forested areas in the Milandre catchment likely experience less snow cover compared to the open meadow where snow can accumulate freely. The meadow soils may experience reduced frost effects, meaning that heterotrophic respiration is not impacted to the same extent as the forest soils and can accumulate CO₂ year-round, as seen in borehole Shallow 2. However, during the monitored period, the Milandre catchment experienced little snow cover, suggesting that frost effects were likely minimal overall.

Though CO₂ production is reduced in winter, concentrations of several thousand ppmV are measured at our sites, perhaps due to (i) slow, but persisting metabolic microbial activity at low temperatures (possible even below freezing temperatures as low as −5 to −7.5 °C; Kähkönen et al., 2001), (ii) reduced, but persisting autotrophic respiration of shallow roots, or (iii) reduced transportation of CO₂ due to higher soil water content (Hashimoto and Komatsu, 2006). Similar seasonal trends in CO₂ concentrations have been reported in other studies (Billings et al.; 1998; Pumpanen et al., 2003; Zhang et al., 2023).

The δ¹³C values of CO₂ in shallow boreholes located in both the meadow and forest correspond to the typical isotopic signature of C₃ plants which dominate the catchment area, supporting our interpretation of a dominant biogenic source of carbon in the soil gas (Figs. 2 and 3b). Plants using the C₃ metabolic pathway (e.g. most temperate vegetation) produce carbon with a δ¹³C of ∼ −23 to −34 ‰ (Staddon, 2004). Furthermore, due to the nature of the sampling, the soil was not sampled in steady-state conditions, potentially disturbing any existing CO₂ trends with depth. Another consideration is the symbiotic relationship between arbuscular mycorrhizal fungi (AMF) and plant roots, which is ubiquitous globally, occurring in 80 % of terrestrial plants (Fitter et al., 2000). AMF grows on plant roots and forms arbuscules within root cells and hyphae which extend into the soil (Zhu and Miller, 2003). Plants benefit from this relationship because of the increased uptake of nutrients such as phosphorus and nitrogen and improved resistance to environmental stress through more stable soil structure and enhanced water absorption (Kiers et al., 2011; Boyno et al., 2024). The fungi benefit through access to plant-derived carbohydrates for energy and a stable environment for growth and reproduction (Kiers et al., 2011). The respiration of AMF can contribute to the subsurface sequestration of carbon (Treseder and Allen, 2008) and may indeed contribute to the subsurface CO₂ which we measured at Milandre cave. There appears to be a temperature dependance of roots with AMF colonization, where respiration rates increase at higher temperatures (Atkin et al., 2008). We observe an increase in CO₂ concentration during the higher-temperature summer months (Figs. 3a and 5a), which may be caused in part by a contribution of CO₂ from enhanced respiration of AMF. The isotopic composition of the CO₂ produced by AMF respiration is not currently well constrained. Studies analysing the δ¹³C of individual AMF spores showed that the δ¹³C broadly corresponded to the C₃ or C₄ plant that it was associated with (Nakano et al., 1999). Therefore, it can be assumed that CO₂ produced by AMF respiration is similar to that of the C₃ host-plant-dominated catchment of Milandre cave. As carbon is fixed from modern atmospheric CO₂ during host plant photosynthesis, the F¹⁴C signature should be modern (∼ 1). Hence, the respiration of AMF may contribute to the modern CO₂ measured in the boreholes in both meadow and forest sites (Figs. 3 and 5).

Though both meadow and forest boreholes have C₃ plant δ¹³C signatures, the δ¹³C values of CO₂ in meadow soils are on average lower than those in the forest soils. Previous analysis of soil CO₂ from various karstic sites in Slovenia showed a similar trend, whereby the δ¹³C values in meadow sites were ∼ 1 ‰ lower compared to forested sites (Stickler et al., 1997). The meadow and forested areas at our site have developed contrasting soil compositions, with the meadow soils deeper, more compacted, and with a higher organic content than the forest soils, which are shallow, unconsolidated, and with more roots. Comparisons of root-free and root-containing soils have suggested that root respiration contributes lower δ¹³CO₂ to the soil gas (Diao et al., 2022). This is the opposite trend to what our findings would suggest. Efforts to disentangle the δ¹³C of autotrophic and heterotrophic contributions to soil CO₂ have shown a range of results, where autotrophic respiration can produce the same (Wu et al., 2017), higher (Moyes et al., 2010), or lower (Risk et al., 2012) δ¹³C than heterotrophic respiration. However, different soil compositions can affect soil moisture, which can cause increases (Unger et al., 2010), decreases (Powers et al., 2010), and no change (Diao et al., 2022) in soil δ¹³C. It is therefore difficult in the context of this study to determine the exact reason for the difference in δ¹³C between surface covers. Interestingly, we also do not find evidence for contributions to subsurface gas from C₄ plants, possibly due to the absence of sampling locations in the fields in the catchment where C₄ crops are grown, suggesting limited lateral gas transport.

The CO₂ F¹⁴C of meadow boreholes in Shallow 2 shows statistically significant seasonal behaviour (rho = 0.67, p = 0.001, n = 30), with more aged CO₂ dominating the winter and spring months (F¹⁴C ∼ 0.88) and modern CO₂ in summer (F¹⁴C ∼ 1.00) (Fig. 3c). These seasonal shifts suggest a potential influence of temperature-sensitive processes on soil carbon dynamics, with higher autotrophic respiration rates and heterotrophic respiration of very recently fixed photosynthates dominating in the warm months (Campeau et al., 2019). Though typically the productive growing season for meadow grasses occurs in spring (Wingler and Hennessy, 2016), there is a lag in the F¹⁴C response, which peaks between June and October. This may be due to the time it takes for enough modern CO₂ to accumulate in the soil and dominate the average gas age. Similar impacts of seasonality on the contributions of young carbon pools to the soil CO₂ budget have been observed in soils of various environments whereby autotrophic respiration and heterotrophic decomposition of younger carbon dominate in warmer months (Trumbore, 2000; Chiti et al., 2011; Vaughn and Torn, 2018). A contribution of fossil carbon from the dissolution of carbonate bedrock fragments in the soil zone could also be considered to explain the low soil CO₂ F¹⁴C values in winter. However, such a bedrock contribution would result in a concomitant increase in δ¹³C, as marine carbonates, which typically constitute the karstified host rock, are characterized by δ¹³C of values of ∼ 0 ‰ (Planavsky et al., 2015) and up to +2 ‰ in the region of Milandre cave (Weissert and Mohr, 1996) (Fig. 2). This would contribute a gaseous δ¹³C ∼ 9.54 ‰ lighter when at equilibrium with DIC from carbonate dissolution at 10.5 °C (Mook et al., 1974). Since we do not observe increases in δ¹³C of soil gas during the periods of lower F¹⁴C, this is an unlikely scenario. A further explanation for the seasonal variability could be due to the seasonal ventilation dynamics in the karst system, which promote upward movement of cave air during the winter and downwards movement of soil air in summer (see Sect. 5.4).

The modern F¹⁴C CO₂ signature in the Shallow 1 forest borehole (Fig. 3c) suggests a year-round dominant source of CO₂ from autotrophic root respiration and the decomposition of very recently fixed soil organic matter. As the root systems are deeper and more developed in the forest soils, there is a higher input of modern carbon compared to the meadows. The soils in the forested area are shallow (in some areas < 10 cm deep) and are dominated by large pebbles and fragments of the carbonate bedrock and shallow roots. Thus, shallow-depth carbon sequestration may be suppressed in the forest soil, resulting in a CO₂ profile dominated by deeper root respiration (Hasenmueller et al., 2017; Tune et al., 2020).

5.2 Sources and variability of CO₂ at 5 m depth

CO₂ concentrations at 5 m depth vary considerably over the sampling period and between individual sites (Fig. 5a). CO₂ concentrations are generally higher in the forest boreholes (Deep 1, Deep 2) than in the meadow boreholes (Deep 3, Deep 4), likely due to the high volumes of CO₂ produced by autotrophic respiration of the mature woodland root system, which penetrates deep into the unsaturated zone. The CO₂ δ¹³C values across all boreholes reflect the C₃-dominated ecosystem of the catchment and do not indicate any input from C₄ plants in the nearby fields. Notably, the ecosystem composition of the overlying vegetation (meadow vs. forest) influences the isotopic composition of CO₂ at 5 m in the same way as the shallow boreholes, with slightly lower δ¹³C values in the meadow than the forest. The F¹⁴C values of unsaturated zone CO₂ vary across all deep boreholes by ± 0.02 around a mean value of 1.00 (Fig. 5c), suggesting that across the catchment, autotrophic tree root respiration may be the dominant source of CO₂ at these depths (Breecker et al., 2012; Tune et al., 2020). Furthermore, the year-round modern CO₂ F¹⁴C values in the deeper meadow boreholes Deep 3 and Deep 4 compared to the seasonally pronounced signal in the shallower boreholes of Shallow 2 highlight the potential influence of the meadow doline morphology. Boreholes Deep 3 and Deep 4 were drilled on the edge of the doline beneath shallow Leptosols like those found in the forest, while Shallow 2 was drilled fully within the doline, where the soil is thicker and much richer in organic matter. Unfortunately, there is no unsaturated zone depth borehole installed within the lower basin of the meadow to verify whether the seasonal signal measured at Shallow 2 translates to the deeper subsurface.

Whilst many studies point to the degradation of exported aged organic matter being the main source of karstic ground air (Breecker et al., 2012; Mattey et al., 2016; Bergel et al., 2017), unsaturated zone borehole CO₂ in the Milandre cave catchment has predominantly modern F¹⁴C values and does not show seasonal isotopic variability. This suggests that the majority of the CO₂ at 5 m depth is likely supplied by contemporaneous C₃ tree root respiration. Similarly modern (106.4 pMC (percent modern carbon, corresponding to 1.064 in F¹⁴C)) biologically sourced CO₂ was also detected in boreholes from Nerja cave, Spain (Vadillo et al., 2016). Modern but not substantially above modern F¹⁴C also excludes the contribution of substantial amounts of decadal-aged soil material containing bomb spike carbon due the enrichment in ¹⁴C by thermonuclear weapons testing during 1950s and 1960s (Trumbore, 2000; Shi et al., 2020). This result supports the body of literature stating that, in addition to the export and respiration of older carbon, large volumes of modern carbon are also transported into the unsaturated zone, likely in large part through deep root respiration (Breecker et al., 2012; Campeau et al., 2019; Tune et al., 2020). This suggests that deep roots, as found in mature forest ecosystems, might be more important than previously thought in contributing to the ground air budget.

5.3 Ventilation-driven isotopic variation in cave air

The CO₂ dynamics in the cave show distinct variations in concentration and isotopic composition between the two sampling points Upstream and Downstream (Fig. 6a–c). Cave ventilation dynamics modulate the mixing ratio between the atmosphere and the other contributing pools and are usually the dominant source of variability (Kukuljan et al., 2021; Buzjak et al., 2024). Recent modelling of airflow dynamics in Milandre cave found that outside temperature controls 95 % of flow variability (Garagnon et al., 2022). During colder months when the outside temperature is ≤ 8 °C, the air in the cave flows topographically upwards from the northern downstream entrance close to the Downstream site to the higher entrance closer to Upstream. When the temperature increases above 8 °C, the ventilation regime reverses. This seasonality is reflected in the alternating concentrations and isotopic composition of the CO₂ at the cave sampling sites as they experience varying amounts of dilution with atmospheric air during the year. CO₂ degassing from the cave river at Upstream reduces the dilution effect in the upper passage. The effects of ventilation can also be observed in the higher-spatial-resolution sampling during the trip through the cave (Fig. 8a–c) whereby the summer regime ventilation resulted in lower concentrations, more positive δ¹³C, and increased F¹⁴C due to atmospheric mixing closer to the upper entrance (Upstream site). Spatially, a marked decrease in CO₂ concentration and F¹⁴C occurs at sampling locations 10 and 11 (∼ 750 to 1000 m from the cave entrance, respectively) though little change is observed in the δ¹³C. This air input could be associated with an older organic matter pool which, however, contributes only little to the entire CO₂ mass flux.

5.4 Sources and variability of CO₂ in cave air

The Keeling plot y intercepts show that the composition of the gas pool (karst endmember) that mixes with atmospheric air in the cave varies seasonally (Fig. 7c). It is important to note that, while the isotopic composition of cave CO₂ is influenced by ventilation-driven dynamics, the Keeling plot y intercept represents the composition of the subsurface gas that mixes with the atmospheric air during ventilation and is not affected by it. The δ¹³C of the karst endmember is mostly stable, varying ± 2 ‰ around a mean value of −26 ‰, a value typically associated with the CO₂ produced by the C₃ plants which dominate the catchment ecosystem. This value is also very similar to the δ¹³C CO₂ values of soil and unsaturated zone gas and reinforces the notion of the karst endmember being dominated by ecosystem respiration. Previous work at Milandre cave investigating aquifer dynamics using ²²²Rn suggested that the majority of the cave CO₂ comes from the overlying soils (Savoy et al., 2011). The F¹⁴C of the Keeling plot y intercept, on the other hand, shows seasonal variability with distinctly more “modern” (closer to F¹⁴C = 1) during the summer and “aged” (F¹⁴C < 1) CO₂ during the winter. The transition from an F¹⁴C of 0.94 to 0.88 represents a change in the mean apparent age of CO₂ from ∼ 100 to 1000 years, well beyond the assumed residence time of water in the unsaturated zone of ∼ 5.5 to 6.6 years (Affolter et al., 2020).

Several mechanisms may be responsible for the seasonal F¹⁴C fluctuations of the karst endmember. A contribution of carbon from the oxidation of methane, either from methanotrophs or through the evolution of pore fluid waters, must be considered. Methane can be produced by both thermogenic (originating from ancient organic-rich sediments) and biogenic (particularly through methanotrophic bacteria) processes (Suess and Whiticar, 1989). The δ¹³C CO₂ range of published values for the CO₂ resulting from methane oxidation from both sources is broad, with biogenic spanning from ∼ −20 ‰ to −90 ‰ and thermogenic from ∼ −30 ‰ to −65 ‰ (Suess and Whiticar, 1989). The δ¹³C CO₂ produced by thermogenic sources can also evolve to less negative values of ∼ −25 ‰ to −34 ‰ in settings where there is nearly total conversion of CH₄ to CO₂ (Schoell, 1980; Frieling et al., 2016). Methane that undergoes biogenic oxidation is typically sourced from the atmosphere, meaning the resulting CO₂ will have a modern F¹⁴C signature (F¹⁴C ∼ 1). In contrast, CO₂ produced by thermogenic methane oxidation will have a fossil F¹⁴C signature because the methane originates from ancient organic-rich sediments that contain no measurable radiocarbon (F¹⁴C ∼ 0). While measuring methane concentrations quantitatively was not possible with our analytical setup, qualitative data show that methane concentrations in Milandre cave air and catchment borehole samples were exceedingly low (∼ 90 to 970 ppb), suggesting methane oxidation may occur in the Milandre subsurface environment. Methanotrophic bacteria have been identified in 98 % of soils sampled from 21 North American caves at a mean bacterial abundance of 0.88 % (Webster et al., 2022). Methanotrophic bacteria in caves significantly contribute to methane consumption, with a study in Pindal Cave (NW Spain) showing how methanotrophic bacteria associated with moonmilk deposits remove 65 % to 90 % of atmospheric methane entering the cave (Martin-Pozas et al., 2022). As our cave and borehole CO₂ samples have a δ¹³C of ∼ −26 ‰ and a variable F¹⁴C in the cave air and shallow meadow samples, we cannot rule out a contribution of methane oxidation from biogenic processes actively occurring in Milandre. However, this would not explain the depleted radiocarbon signature that we measured episodically in the cave air. Methane evolved abiogenically from pore fluids in Oxfordian marls may have an isotopic ratio similar to that measured in our samples if thermogenically derived with a dominant conversion of CH₄ to CO₂. However, it is unlikely that this process occurred in the Milandre cave sediments given the lack of organic-rich sedimentary rocks in the stratigraphy. Thus, while we cannot quantify the importance of methane oxidation for the subsurface carbon budget at our site, it is unlikely that methane oxidation processes are responsible for the variation in radiocarbon ages in cave air CO₂.

Despite the apparent limiting factors for microbial growth in caves, it is known that different organisms are able to thrive within this niche including bacteria, archaea, algae, and fungi. In terms of bacteria, Proteobacteria is typically reported as the dominant group in caves (Zhou et al., 2007; Tomczyk-Żak and Zielenkiewicz, 2016). However, in some cases Actinobacteria is the largest group of microbes (Wiseschart et al., 2018). In two Chinese caves with high CO₂ concentrations (> 3000 ppm), Proteobacteria, Actinobacteria, Bacteroidetes, and Thaumarchaeota were the dominant phyla (Chen et al., 2023). Hundreds of fungi taxa have been described in caves globally, with Ascomycota being the most commonly reported (Gherman et al., 2014). Though not yet thoroughly investigated, diverse communities of protists have also been identified within caves (Gogoleva et al., 2024), as well as archaea (Biagioli et al., 2023) and algae (Kosznik-Kwaśnicka et al., 2022). As biological processes generally favour the incorporation of lighter ¹²C compared to ¹³C, the respiration of these organisms would contribute an organic signature (lower δ¹³C and a F¹⁴C corresponding to the age of the carbon source being consumed, which can range from modern to more aged F¹⁴C) to the CO₂ composition of cave air at Milandre cave. With the isotopic analysis used in this study, and because we did not analyse the microbial community diversity present in Milandre cave and its catchment, the CO₂ contributions from different types of microbial processes cannot be distinguished from one another. We are only able to distinguish between “organic” processes (in-cave microbial, soil, root, and fungi respiration), which produce a CO₂ with a lower δ¹³C signature and a range of F¹⁴C depending on the age of the organic material consumed, and “inorganic” processes (carbonate degassing and atmospheric ventilation), which contribute CO₂ with a higher δ¹³C and a range of F¹⁴C.

Varying contributions of host-rock-derived carbon to cave air through shifts in the host rock dissolution regime can affect the isotopic composition of cave air. Dissolution of the host rock carbonate can occur under a wide range of conditions between two extreme cases: completely closed and completely open systems. In a completely closed system, the aqueous solution becomes isolated from the soil CO₂ reservoir after passing into the unsaturated zone, resulting in a theoretical 50:50 contribution ratio of carbon ions in solution from the carbonate host rock and from soil air and overall decreasing F¹⁴C to as low as 0.5 (Fohlmeister et al., 2011; Milanolo and Gabrovšek, 2015). On the other hand, in a completely open system the aqueous solution can continuously exchange with an unlimited soil CO₂ reservoir, resulting in most carbon atoms being sourced from the soil CO₂ reservoir and the rock contribution being minimal (Fohlmeister et al., 2011). Most natural systems exist in an intermediate state between completely open and completely closed settings.

If the seasonal variability in CO₂ F¹⁴C in Milandre cave reflected seasonal shifts in the “openness” of the system, summers (with higher F¹⁴C) would be characterized by a more open system with higher rates of exchange between the aqueous solution and the soil CO₂ reservoir. The observed shift to lower F¹⁴C in winter would reflect a more closed system with a higher contribution of the F¹⁴C free carbonate bedrock. Shifts in the open–closed continuum would likely result in changes in the δ¹³C of the endmember (Fairchild and Baker, 2012), which is not seen in the cave air, with an expected lower δ¹³C in a more open system and an increase in a more closed system because of the higher carbonate host rock carbon contribution (δ¹³C ≈ 0 ‰). However, varying amounts of isotopic fractionation can occur when CO₂ is degassed from the drip water DIC pool into the cave atmosphere, potentially influencing the isotopic composition of the measured cave air CO₂ (Mickler et al., 2019).

We constructed a mixing model which allows us to evaluate how much variability in the Keeling plot y intercept (i.e. the second source of carbon to cave air besides the atmosphere) can be explained by variation in DIC contribution from dissolution regime changes and degassing fractionation. This model assumes that the Keeling plot y intercept is itself a mixture of two endmembers that contribute CO₂ to cave air, the modern soil (F¹⁴C_soil = 1) and the DIC (F¹⁴C_DIC), where the F¹⁴C_DIC represents a mixture of dissolved soil CO₂ and DIC from carbonate mineral dissolution, which we can constrain with our measured DIC values. The F¹⁴C_soil represents direct input of CO₂ soil gas into the cave.

The mixing ratio between F¹⁴C_soil and F¹⁴C_DIC results in the F¹⁴C of CO₂ of the karst endmember (F¹⁴C_cave). Thus, the cave air CO₂ F¹⁴C can be expressed as follows (where f is fraction, defined as between 0 and 1).

$$\begin{array}{}\text{(3)}& {\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{cave}}=f_{\text{soil}}{\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{soil}}+f_{\text{DIC}}{\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{DIC}}\end{array}$$

Since $f_{\text{DIC}}=\mathrm{1}-f_{\text{soil}}$, the equation can be written as follows.

$$\begin{array}{}\text{(4)}& {\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{cave}}=f_{\text{soil}}{\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{soil}}+(\mathrm{1}-f_{\text{soil}}){\mathrm{F}}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{DIC}}\end{array}$$

For F¹⁴C, fractionation is generally negligible during carbonate dissolution and degassing compared to δ¹³C due to unit conversion, i.e. % vs. ‰ (Fohlmeister et al., 2011). We use the mixing calculation (Eq. 4) to estimate the contributions of f_soil and f_DIC using the highest (0.94), an intermediate (0.90), and the lowest (0.88) measured cave CO₂ F¹⁴C values (F¹⁴C_cave) and F¹⁴C DIC (F¹⁴C_DIC) based on the highest (0.98), intermediate (0.88), and lowest (0.78) measured values of drip water DIC. We find that when F¹⁴C_DIC = 0.98, the mixing model results in impossible scenarios with f_soil > 1 and f_DIC < 0 because this value is higher than the highest F¹⁴C_cave. When F¹⁴C_DIC = 0.88 it results in a highly unlikely scenario based on our observations and the literature where f_soil=0 and f_DIC=1, as this corresponds to the lowest F¹⁴C_cave value (Fig. 9a). Only the scenario with the lowest F¹⁴C_DIC = 0.78 produced realistic mixing fractions. Overall, this implies that there is a discrepancy between the measured F¹⁴C_DIC and what is derived by a simplified mixing model based on the F¹⁴C of cave CO₂, with F¹⁴C_DIC values too high to result in the F¹⁴C_cave composition that we observe. Thus, carbonate dissolution and subsequent degassing of CO₂ from drip water DIC to cave air are unlikely to be the dominant processes explaining the shifts in cave air F¹⁴C over time.

As our mixing ratios derived from measured DIC values were mostly not realistic according to the two-endmember mixing scenario, we considered how a wider theoretical range of F¹⁴C_DIC values (from closed to mainly open system dissolution conditions, 0.50 to 0.90 F¹⁴C) would affect the contributing fractions of f_soil and f_DIC. Variations in F¹⁴C_DIC from 0.50 to 0.85 result in viable mixing fractions across the range of F¹⁴C_cave, with the f_soil (20 % to 76 %) and f_DIC (24 % to 80 %) differing widely (Fig. 9b).

As the mean value of cave CO₂ δ¹³C is −26 ‰ (we take a range from −25 ‰ to −27 ‰), we created a similar two-endmember mixing calculation for δ¹³C values to test whether the mixing ratios derived from F¹⁴C translate to observed δ¹³C values. For this we used the mean δ¹³C of soil CO₂ (δ¹³C_soil) = −26 ‰ and the δ¹³C of the DIC (δ¹³C_DIC) = −15 ‰ (estimated from a mean of the δ¹³C produced during the F¹⁴C DIC AMS measurement; see Appendix B).

Fractionation effects during degassing must be considered when modelling δ¹³C. If equilibrium fractionation between DIC and CO_2(g) occurs, the fractionation factor (${\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{eq}}$) at 10.5 °C between HCO₃− and CO₂ is 9.54 ‰, with the gas being 9.54 ‰ depleted in ¹³C relative to the fluid (Mook et al., 1974). Kinetic fractionation can occur during rapid degassing of CO₂ as well as during rapid precipitation of carbonate minerals (Mickler et al., 2019). If kinetic fractionation occurs during degassing, this results in a greater depletion of the dissolved CO₂ in ¹³C compared to equilibrium fractionation (i.e. a more positive δ¹³C_DIC). The possible ranges reported for the composition of CO₂ generated from kinetic fractionation during degassing in caves vary widely, with some reporting ${\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{kin}}$ close to that resulting from equilibrium fractionation (Dulinski and Rozanski, 1990) and others a range of 10 ‰ to 65 ‰ (Mickler et al., 2019). The extent of kinetic fractionation is expected to depend on the amount of degassing that occurs, which is a function of the disequilibrium between the fluid and the surrounding air (Frisia et al., 2011). In our system with limited geochemical data regarding fluid compositions, it is difficult to determine what extent of kinetic fractionation may have occurred. Hence, we first assume that the fractionation is < −9.54 ‰ (i.e. greater compared to equilibrium conditions) and take a moderate but arbitrary value ballpark of ${\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{kin}}$ = 30 ‰.

A mixing calculation is then applied to assess if the measured range of the karst endmember δ¹³CO₂ given by the Keeling plot can be consistent with the extended range of F¹⁴C_DIC (0.50 to 0.85 F¹⁴C) and measured F¹⁴C_cave under conditions with equilibrium degassing and kinetic degassing.

Equilibrium degassing:

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{cave}}=f_{\text{soil}}{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soil}}+f_{\text{DIC}}\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}-{\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{eq}}\right)\end{array}$$

Kinetic degassing:

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{cave}}=f_{\text{soil}}{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soil}}+f_{\text{DIC}}\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}-{\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{kin}}\right)\end{array}$$

For the maximum (0.94), intermediate (0.90), and minimum (0.88) values of the endmember F¹⁴C, only equilibrium fractionation during degassing can explain the observed δ¹³C_cave value of approximately −26 ‰ (Fig. 9c). Kinetic fractionation scenarios result in very negative δ¹³C_cave values which do not fit those measured, except when using very small kinetic fractionation factors similar to that of the equilibrium fractionation (${\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{HCO}}_{\mathrm{3}}-{\mathrm{CO}}_{\mathrm{2}},\mathrm{kin}}$ of −10 ‰, see Appendix C). It is possible that kinetic fractionation occurs between Milandre cave waters and the cave air due to atmospheric ventilation which occurs year-round. The atmospheric air which enters from the upper entrance in the summer and the lower entrance in the winter creates a concentration gradient between the CO₂ present within cave waters and the concentration in the cave atmosphere. However, the rapidly increasing pCO₂ with distance from the cave entrance (Fig. 8a) implies that the concentration gradients required for kinetic fractionation do not occur throughout most of the cave passage and that kinetic influences are probably minimal.

Overall, the F¹⁴C_DIC measured beyond F¹⁴C = 0.90 is too high to reproduce our measured cave air CO₂ F¹⁴C and δ¹³C. The F¹⁴C_DIC is based on the analysis of nine drip water sites of varying drip rates and ecosystem coverage over the course of 1 year. It is possible that the sampling setup did not capture the true range of F¹⁴C_DIC in the drip water, though this hypothesis is highly unlikely due to the diversity in drip rate and hydrological response of drips measured, and we believe our DIC data to be representative of the system.

The F¹⁴C_DIC of the cave river was not measured; however, previous studies investigating the Milandre cave river have suggested that in non-flooding conditions the river is fed by slow diffuse flow as the unsaturated zone acts as a buffer, storing water and dampening the signal of normal rainfall. During flooding events, the unsaturated zone aquifer is bypassed and the stream is mainly fed through fracture flow (Perrin et al., 2003; Savoy et al., 2011). A diffuse-flow-fed river is likely to have isotopic characteristics similar to that of the drip waters, which have higher F¹⁴C, suggesting more open-system dissolution conditions. In contrast, times of flooding and soil saturation result in less exchange with the soil gas reservoir, leading to dissolution under more closed conditions (Perrin et al., 2003; Savoy et al., 2011). However, the drip water F¹⁴C_DIC is relatively stable over 1 year of monitoring in a variety of hydrological conditions, implying that variations in host rock dissolution and degassing dynamics are not the explanation for the variation in F¹⁴C and δ¹³C of the endmember.

Figure 10 Comparison between the (a) δ¹³C and (b) F¹⁴C of the shallow depth boreholes in the forest (Shallow 1, green) and meadow (Shallow 2, brown), as well as the Keeling plot endmembers δ¹³C (blue) and F¹⁴C (orange) also show in Fig. 7. Summer (yellow) and winter (blue) seasons are highlighted by the coloured background. The direction of cave ventilation is shown by the arrows with downwards ventilation (from Upstream to Downstream) in summer and upwards ventilation (Downstream to Upstream). Periods of decoupling between the F¹⁴C of the meadow borehole and the cave endmember are shown in black boxes.

The meadow soil boreholes in the doline (Shallow 2) have similar patterns in CO₂ composition to the karst endmember, with a δ¹³C of −26 ‰ and variable F¹⁴C (Fig. 10a and b). In 2022, the F¹⁴C in the meadow soil and karst endmember align closely over the seasonal cycle, but then decouple slightly in 2023 during the transitional periods between winter and spring where the karst endmember F¹⁴C changes before the meadow. This suggests that, instead of the soils influencing the cave air due to seasonal changes in soil respiration, the cave ventilation regime may be influencing the soil CO₂ F¹⁴C. The cave ventilates upwards in winter from Downstream to Upstream in a chimney effect (and vice versa in summer). Here, the lower-F¹⁴C karst endmember CO₂ may flow upwards through the highly fractured unsaturated zone of the doline and ventilate the meadow soils, resulting in a mixture of low-F¹⁴C CO₂ from the karst endmember and modern CO₂ from the soils. This possibly reflects a contribution of older CO₂ from within the karst itself, likely from a “ground air” reservoir of older respiring organic material in the unsaturated zone as found by several other studies (Breecker et al., 2012; Noronha et al., 2015; Mattey et al., 2016; Bergel et al., 2017). It is possible that several reservoirs of older organic material exist throughout the downstream catchment of Milandre, which may contribute the low-F¹⁴C CO₂ to the ventilating cave air in winter. In summer the ventilation reverses, and the higher-F¹⁴C CO₂ from the soil is transported downwards into the cave, increasing the F¹⁴C of the cave endmember CO₂. The F¹⁴C of the karst endmember is always lower than the meadow soil F¹⁴C, implying that there is always mixing between the soils and the karst endmember (which additionally also reflects the DIC contribution). Similar upwards transport of CO₂ in karst boreholes has been observed in Nerja Cave, Spain (Benavente et al., 2010; Benavente et al., 2015), and in the Gibraltar karst (Mattey et al., 2016). We do not observe significant effects of the winter upwards ventilation in the F¹⁴C CO₂ in any of the deeper boreholes in the meadow or any of the boreholes in the forest. This could be related to the location of Shallow 2 in the doline, whilst the rest of the boreholes are situated higher on the banks (Deep 3 and Deep 4) or further away from the doline in the forest (Deep 1 and Deep 2). The topography difference between the higher forest and lower meadow likely reflects influences of the secondary porosity in the bedrock, which is likely more highly fractured in the meadow doline formation compared to the forest. As the doline structure likely has a higher secondary porosity and is more highly fractured than the higher bedrock, ventilation effects may be more important here than at the other locations.

The accumulation of the older organic carbon in the unsaturated zone likely occurs over time. Soil layers in karst regions often accumulate on bedrock with deep vertical fractures (Cheng et al., 2023). This decreases the ability of long-term storage of older soil carbon pools due to preferential fracture flow, also known as subsurface leakage, which destabilizes the oldest soil layer and leads to export of older carbon into the unsaturated zone (Sánchez-Cañete et al., 2018; Wan et al., 2018). Several studies have acknowledged this potential source of CO₂ produced autochthonously inside the unsaturated zone by the microbial degradation of this old soil, plant, or root material which was previously washed into the unsaturated zone, contributing older CO₂ to the ground air mix (Noronha et al., 2015; Bergel et al., 2017; Ding et al., 2023).

6 Conclusion

We conducted a comprehensive investigation over 2 years into the dynamics of CO₂ concentrations and δ¹³C and F¹⁴C values across different land covers and soil types within the critical zone at Milandre cave in northern Switzerland. Our analysis found distinct seasonal fluctuations in CO₂ levels in soil boreholes in both forest and meadow areas, with higher concentrations during summer compared to winter. Low δ¹³C values of ∼ −26 ‰ indicate respiration of C₃ plants, which dominate the catchment area. The stable modern F¹⁴C of ∼ 1.00 in shallow forest boreholes indicated a year-round modern CO₂ contribution from tree and plant roots to the subsurface, consistent with a well-ventilated soil. On the other hand, samples from meadow boreholes situated in a doline with a thick, well-developed soil cover displayed seasonality in F¹⁴C from ∼ 0.88 in the winter and spring to ∼ 1.00 in summer.

In the unsaturated zone we observed substantial variability in CO₂ concentrations across the catchment (∼ 3000 to ∼ 37 000 ppmV), with higher concentrations in boreholes with forest cover than those with meadow cover. Despite differences in surface vegetation, the isotopic compositions of δ¹³C and F¹⁴C remained stable in both forest and meadow environments, reflecting the dominant modern C₃ vegetation signature contributing to the ground air. The year-round modern unsaturated zone CO₂ in meadow boreholes Deep 3 and Deep 4 contrasts with the seasonally attenuated signal in the thicker meadow soils.

Large variations in cave air CO₂ concentrations (atmospheric to ∼ 30 000 ppmV), δ¹³C (∼ −9 to −25 ‰), and F¹⁴C (0.88 to 0.98) are controlled by seasonal temperature-driven cave ventilation. The seasonal changes in the isotopic composition of the Keeling-plot-derived karst endmember cannot be explained with the observed DIC F¹⁴C range, suggesting that bedrock dissolution and degassing dynamics are likely not the cause of the isotopic variation in the cave air. The isotopic characteristics of the cave air are comparable to those of the meadow doline soil. The seasonal variability in karst endmember CO₂ F¹⁴C precedes that in the soil and is always older. This suggests that cave ventilation contributes older CO₂ to the doline soils, likely sourced from a reservoir of aged organic material in the unsaturated zone, during the upward winter ventilation regime. In summer, the reversed downwards ventilation contributes younger soil CO₂ to the cave, mixing with the older unsaturated zone reservoir.

Our study provides important insights into how carbon is transported to the subsurface in karstic critical zones. Understanding these processes is crucial for accurate estimation of the size of subsurface CO₂ pools (ground air) and to refine terrestrial CO₂ budgets. Moving forward, it would be beneficial to undertake more detailed investigation of CO₂ transport into the subsurface by implementing ¹⁴C analysis to further constrain the sources of ground air. Furthermore, higher-resolution monitoring over periods of interest could assist us in understanding short-term variations.

Appendix A

Figure A1 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from soil depth borehole Shallow 1. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Figure A2 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from soil depth borehole Shallow 2. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Figure A3 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from unsaturated zone depth boreholes Deep 1 and Deep 2. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Figure A4 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from unsaturated zone depth boreholes Deep 3 and Deep 4. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Figure A5 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from the Downstream cave site. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Figure A6 Spearman's rank correlation matrix displaying the Spearman's rho of CO₂ concentration, δ¹³C and F¹⁴C, mean monthly temperature (MMT), and mean monthly precipitation (MMP) from Fahy weather station from samples from the Upstream cave site. Parameters with a significant correlation (positive or negative) are bolded, and the strength of the correlation is indicated by the blue to red colour bar.

Appendix B

Table B1 Drip water dissolved inorganic carbon F¹⁴C and AMS-derived δ¹³C from the actively dripping galleries near Upstream and Downstream.

Appendix C

Figure C1 Cave air δ¹³C for varying DIC F¹⁴C and cave air F¹⁴C values with equilibrium fractionation at −9.54 ‰ (solid lines) and kinetic fractionation at −10 ‰ (dashed lines).