Abiotic processes control carbon dioxide dynamics in temperate karst lakes

Study site and sampling

The Ruidera Lake district is located in central Spain (Autonomous Community of Castilla-La Mancha; 39°0′-38°52′N and 2°52′–2°48′E) in La Mancha Húmeda Biosphere Reserve (Fig. 1). Ruidera Natural Park extends over 37.72 km², comprising a chain of 15 interconnected tufa lakes developed along the watercourse of the Upper Guadiana River, with deeply karstified massive Jurassic dolostones, and dolomites, overlying Triassic sediments (Ordóñez et al., 2005). Lakes are mainly groundwater-fed from the Campo de Montiel aquifer, which supplies around 90% of water, while the remaining 10% comes from rainfall-runoff events (Montero, 1994; Martínez-Cortina, 2001). The aquifer is enriched with nitrate (>50 mg L⁻¹) from excessive agricultural fertilization used through intensive irrigation farming (Lima, 2008). This implies excess nitrate inputs into the lakes, mainly in those located upstream as La Conceja Lake.

The climate is subhumid temperate with dry summers (Csa, Köppen-Geiger classification), having two well-distinguished seasons: a cold/rainy winter season (from October to February) and a typical Mediterranean dry/warm summer (June to September). The average annual temperature is 15.8 °C, and the mean annual precipitation is 481 mm (AEMET Tomelloso; 39°10′9″N, 3°0′27″E, 662 m a.s.l.). Calcisols and fluvisols are the predominant soils in the area, and the major land uses in the catchment include broad-leaved and mixed forests, meadows and pastures, and farming.

The study was undertaken in three lakes (Table S1), which were chosen as being representative of the lacustrine heterogeneity (i.e., hydrological and ecological characteristics as well as anthropogenic impact—increased urbanization pressure downstream; Table S1): Conceja Lake (LC) as the most pristine groundwater-fed system located upstream of the lake chain; Colgada Lake (CG), located in the middle section of the lake systems and fed in addition by surface water and treated wastewaters from the nearby town Ossa de Montiel; and, finally, Cueva Morenilla Lake (CM), located downstream which mainly receives surface waters coming from the other upstream lakes and treated domestic wastewaters, which represent the cumulative impact of the entire lake system.

Physical, chemical, and biological characterization

Sampling campaigns were carried out in summer (June 2021) and winter (February 2022). Monitoring consisted of three complete samplings per day for each lake in summer (∼09:00 h, 15:00 h, and 21:00 h), and two samplings in winter (∼09:00 h and 15:00 h), because of the shorter daylight time. Sampling collected in situ profiles of temperature (T, °C), dissolved oxygen (DO, mg L⁻¹), electrical conductivity (K₂₅, µS cm⁻¹), pH (pH units), and photosynthetic active radiation (PAR, µE m⁻²s⁻¹) measured with a CTD (Seabird SBE 19Plus) and a YSI 6600 probe (at 1 m depth) in the maximum depth zone of the lakes. The mixing layer (Z_MIX) was estimated based on the vertical profiles of T and DO. The euphotic zone (Z_EU) was estimated from the surface down to the depth at which PAR reaches 0.1% of PAR at the surface (SPAR) (Palmer et al., 2013).

Water samples were taken at three depths along the vertical profile (half a meter below the surface, mid-water, and one meter above the lake bed) with a standard Niskin water sampler bottle (5 L). Water samples for nutrient analyses were placed in polyethylene bottles, stored in the dark, and frozen (−4 °C) until analysis. The nitrogen as ammonia (N-NH₄), nitrite (N-NO₂), and nitrate (N-NO₃), as well phosphorus as soluble reactive phosphorus (SRP), were analyzed in a segmented-flow autoanalyzer (SEAL Analytical AutoAnalyzer 3) following standardized methods (APHA, 2005). Water samples (60 ml) for dissolved organic carbon (DOC) were filtered through Whatman GF/F filters (0.7 µm pore size) acidified with H₃PO₄ (40%). The samples for DOC concentration were stored in amber high-density polystyrene bottles and then analyzed using Shimadzu TOC-V CSH equipment coupled with a TNM-1. The gravimetric method determined the total concentration of suspended solids (TSS, seston; Jellison & Melack, 2001). Chlorophyll-a (Chl-a) concentration was determined following the procedure by Marker, Crowther & Gunn (1979).

Net ecosystem production (NEP) was measured at 0.5 m depth using the light-dark bottle method (Wetzel & Likens, 2013). Incubation lasted 4–8 h, and changes in DO concentration were measured with a HACH portable oximeter (HQ40d). Finally, carbonate precipitation was assessed through the calcite saturation index (SIc; APHA, 2005), an indicator of the degree of calcium carbonate saturation in water. SIc was calculated for each lake based on T, K₂₅, pH values, and HCO₃⁻ and Ca⁺ concentrations. SIc values above 1 indicate calcite precipitation.

Dissolved inorganic carbon and CO₂ flux measurements

Triplicated water samples for dissolved inorganic carbon (DIC) were collected at the three mentioned water depths along the vertical profile of the lake, ensuring no bubbles formed. An amount of 60 ml of the water samples was filtered through Whatman GF/F filters, acidified with H₃PO₄ (85%; pH ≈1), and kept cool and in the dark for 24 h, forcing the carbonate system into equilibrium (all DIC converted to CO₂). The headspace equilibration method (Borges, Abril & Bouillon, 2018) using He as carrier gas was used to extract the CO₂ which was then stored in 12 mL Exetainer vials (Labco Ltd., Lampeter, UK) until further analysis. Also, the headspace method, avoiding water sample acidification, was used to measure in triplicate the dissolved CO₂ concentration in the lake water (C_CO2; Goldenfum, 2010). All samples were analyzed by gas chromatography (Agilent model 6890N) equipped with a single-stage dual-packed column, where CO₂ was detected using a thermal conductivity detector (TCD).

Static floating chambers (5 L polyethylene), equipped with 12V electronic fans to ensure gas mixing inside the chamber, were used to measure CO₂ emissions and to compute total evasion fluxes (F_CO2) from the water surface (see below; DelSontro et al., 2016). Chamber deployment is recognized to likely enhance gas transfer through disturbance of the surface boundary layer, but the effect is mainly noted in extremely low-wind environments (Matthews, St.Louis & Hesslein, 2003) as well as using anchored chambers in running waters (Lorke et al., 2015). In Ruidera Lakes, low wind speeds are uncommon, although they can appear between 0:00 and 02:00 h (65% of cases; Sánchez-Carrillo, unpublished data). Moreover, any artificial turbulence created during gas exchange is minor for drifting (floating) chambers (Lorke et al., 2015). Therefore, floating chambers are considered ideal for the environment in which this study took place. Eight chambers were randomly placed in the center of each lake for 30–45 min to monitor the changes in CO₂ concentration inside the chamber. CO₂ concentration measured in the headspace of the chamber displayed a linear trend, as indicative of rate constant emission from diffusive gas flux. However, on some occasions, we observed signs of overpressure inside the chamber, when at the end of the measurements, the gas concentration described a parabola. It only occurred in five cases, when the incubations lasted for 45 min, and the procedure used to avoid these errors was to discard those final gas concentration measurements (2–3) to avoid including these data in the computation of the linear function that defined the gas emission rate during the incubations with the floating chambers. In summer, a syringe was used to extract the gas samples via a butyl rubber stopper in the upper chambers. The samples were stored in Exetainer vials (12 mL) until analysis by gas chromatography. In winter, a Gasera ONE PULSE based on photoacoustic spectroscopy through NDIR-PAS technology (mechanically chopped broadband IR source with optical bandpass filters) was used to measure in situ variations in concentration within the gas chamber. Gas measurements using the Gasera device were taken automatically at 8–9 min intervals through 2 m Teflon tubing, recirculating the measured gas back into the chambers to avoid negative pressure. Prior intercalibration of both analytical methods (chromatography and photoacoustic spectroscopy) yielded comparable control values [R² = 0.987; [Spectros (in ppm CO₂) = 0.99 × chromat (in ppm CO₂)].

The slope of the change in chamber CO₂ concentration (ppm) over the sampling time (t) was used to calculate the F_CO2. The area (A_C) and the volume (V_C) of the chamber were also considered (Eq. 1). The concentration was adjusted for pressure and temperature according to the ideal gas law. (1)${\displaystyle {\mathbf{F}}_{\mathbf{CO}2}=\frac{\Delta \mathbf{C}}{\Delta \mathbf{t}}\left(\frac{{\mathbf{V}}_{\mathbf{c}}}{{\mathbf{A}}_{\mathbf{c}}}\right).}$

Stable isotope analyses and calculations

In each season, triplicate exetainers containing DIC (after acidification with H₃PO₄) and CO₂ (without acidification), collected as explained before, were shipped to the stable isotope lab for δ¹³C analyses (δ¹³C-DIC and δ¹³CO₂; see below). The isotopic ratios were determined by direct injection in a Thermo Fisher Scientific Delta V Plus mass spectrometer (Thermo Fisher Scientific, Waltham, USA) connected to a GasBench system GC combustion interface (CO₂). These analyses were carried out in the Stable Isotope Laboratory at the University of California-Davis (SIF-UC Davis). The gas concentration in water was calculated based on the ratio between the volume of headspace and the volume of water, according to the efficiency of gas extraction (the mean extraction efficiency of DIC was 99% as determined against standard solutions of NaCO₃). For most gases, the heavier isotope presents a higher solubility (Jancsó, 2002).

Since gases with higher solubility reach equilibrium between phases quicker (Garcia-Tigreros et al., 2016), the rate of gas exchange for the heavier ¹³CO₂ isotope should be faster. Lake sediments were also analyzed to determine the δ¹³C of organic matter (δ¹³C-OM) at the Environmental Isotope Laboratory of the University of Arizona, using a Finnigan Delta PlusXL (Thermo Fisher Scientific, Waltham, MA, USA) coupled to a Costech EA (Costech Analytical Tech Inc, Valencia, CA, USA). The isotopic data are described as: (2)${\displaystyle \delta \left(‰\right)=\left(\frac{R_{sample}}{R_{std}-1}\right)\times 1000}$ where R represents the isotopic ratio of the heavy isotope with the lighter isotope (¹³C/¹²C), and R_std refers to the international standard V-PDB (fossil calcite of the Peedee formation belemnite, South Carolina), using as laboratory reference material the NIST 8545 (lithium carbonate; National Institute of Standards and Technology, Gaithersburg, MD, USA).

Carbon isotope fractionation during the precipitation of calcium carbonate was estimated, considering the fractionation of the carbon isotopes ¹²C and ¹³C in the equilibrium system CO₂ (gas)-HCO₃⁻: (3)${\displaystyle CaC{O}_3\left(s\right)+H_{2}C{O}_3\left(aq\right)\rightleftharpoons C{a}^{2+}\left(aq\right)+2HC{O}_3^{-}\left(aq\right)}$       ⇓ (4)${\displaystyle H_{2}O\left(l\right)+C{O}_2\left(g\right)}$ where (s) represents the solid-calcite, (aq) the aqueous solutions, and (l) and (g) the liquid and gas states of H₂O and CO₂, respectively. Both Eqs. (3) and (4) relate to each other, in such a way that when CO₂ is removed from the system (e.g., into the atmosphere), H₂CO₃ breaks apart, and the reaction shifts toward the left to replace lost bicarbonate and then calcite precipitates.

Water temperature dependence of the fractionation for ¹³C between the DIC (aq) and CO₂ (g) was calculated according to the Mook, Bommerson & Staverman (1974) equation: (5)${\displaystyle {\delta }^{13}C-C{O}_{2-fractionation}\left(\text{\%}\right)=\left({\delta }^{13}C-DIC\right)+23.644-\left[\frac{9701.5}{\left(T+273.15\right)}\right]}$ where T is the water temperature expressed in (° C), and the δ¹³C-DIC is the measured value of DIC.

The Miller-Tans Miller & Tans (2003) and Keeling plots Keeling (1958) were used to explore observed δ¹³C-CO₂ values and identify potential sources and fate of our lakes. Both graphical techniques are based on the principle of conservation of mass and assume mixing of isotopically distinct C sources: (6)${\displaystyle {\delta }^{13}Cobs\times Cobs=\left({\delta }^{13}CS\times CS\right)+\left({\delta }^{13}CB\times CB\right)}$ where δ¹³Cobs and Cobs are the observed C isotopic composition and concentration, respectively, in each sample, resulting from a mixture of the C isotopic composition and concentration of the source ”S,” for example, biogenic or geogenic DIC source (δ¹³CS × CS) and background “B” (δ¹³CB × CB), represented in the atmospheric CO₂. This principle can be expressed as a linear relationship as in Eq. (7) in terms of the Keeling plot regression and as in Eq. (8) in the Miller-Tans plot regression: (7)${\displaystyle {\delta }^{13}Cobs=CB\left({\delta }^{13}CB-CS\right)\times \left(\frac{1}{Cobs}\right)+{\delta }^{13}CS}$ (8)${\displaystyle {\delta }^{13}Cobs\times Cobs=\left({\delta }^{13}CS\times Cobs\right)-\left[{\delta }^{13}CB\left({\delta }^{13}CB-{\delta }^{13}CS\right)\right].}$ δ¹³CS is found in the intercept of the linear regression in Eq. (4) or the slope of the linear regression in Eq. (6). The differences between the two regression models imply that the lake δ¹³C-DIC values in equilibrium with atmospheric CO₂ (δ¹³CB × CB) must remain fixed across observations in (Eq. 5). However, this requirement can be disregarded in Eq. (6) since δ¹³CB × CB is found in the residual variation of the regression line. The Miller-Tans plot technique is particularly suitable for approximating δ¹³Cs when including observations from multiple lakes that have undergone various degrees of CO₂ evasion. Both models assume linearity without further fractionation processes (Pataki et al., 2003). These assumptions can be violated in inland waters since kinetic fractionation processes occur when CO₂ is outgassed from the lake water, and other in-lake biogeochemical processes may also be involved. Accordingly, care should be taken when interpreting δ¹³CS derived from these mixing equations.

Statistical analyses

The descriptive statistics (mean, median, and coefficient of variation) were calculated for each lake to get an overall dataset representation. The normality and homogeneity of variance of all measured variables (T, DO, pH, K₂₅, Z_EU, SRP, N-NO₃, N-NO₂, N-NH₄, TSS, Chl-a, NEP, DIC, DOC, C_CO2, and F_CO2) were examined using the Shapiro–Wilk test (Shapiro & Wilk, 1965) by the shapiro.test() function in the “stats” package (R Core Team, 2021). Due to the non-normal distribution of the variables, we used Mann–Whitney and Kruskal-Wallis tests to determine significant differences between the seasons (summer and winter), lakes (LC, CG, and CM), times of day (morning, evening, and night), and water sample depths (surface, middle, and bottom). The tests were performed by the functions wilcox.test() and kruskal.test() in the “stats” package (R Core Team, 2021). A significance level of p < 0.05 was used to determine the differences between the data sets. The variations among factor levels were tested using Dunn’s post-hoc test by the function dunnTest() in the “FSA” package (Ogle et al., 2023). We then performed a principal component analysis (PCA) for each season using the “FactoMineR” package (Lê, Josse & Husson, 2008) to remove spurious variables not representative in the majority of DIC and C_CO2-associated processes. A redundancy analysis (RDA) using the rda() function from the “vegan” package (Oksanen et al., 2022) was also undertaken to specifically confirm the impact of environmental variables on the variability of C_CO2for each season, but findings were statistically inconclusive and were not included in results. Finally, we conducted multiple stepwise linear regression analyses, considering T, DO, pH, K₂₅, TSS, Chl-a, nutrients, DIC, and DOC as independent variables and C_CO2 as dependent. To perform the multiple regression models (MRMs) analysis, we used the “olsrr” package’s ols_step_forward_p() function (Hebbali, 2020). If necessary, the variables were transformed using a logarithmic function to meet the requirements of normality and heteroscedasticity. Finally, ordinary least squares linear regression models were performed to determine the correlations between C_CO2 and F_CO2, and in the Keeling and Miller-Tans plots. All statistical analyses and graphs were performed using RStudio (2021.09.2; RStudio Team, 2021).

Physical, chemical, and biological characterization

Table S2 summarizes the main physical and chemical variables of the lakes in a seasonal comparison. Although T ranged from 9.1 to 25°C and was about 12 °C higher in summer than in winter, there was no thermal stratification in any lake during the summer season. Z_MIX and Z_EU in all lakes were equivalent to the entire water column. Nitrogen fractions were significantly higher in lakes during summer whereas SRP showed the opposite pattern (Table 1). By far, the N-NO₃ was the most abundant inorganic nitrogen fraction for both seasons.

Seasonal significant differences were found in DOC, Chl-a and NEP rates, but the post-hoc tests revealed similar values recorded during summer and winter in the lakes (Table 2). SIc values were similar in both seasons with values ranging from 10.0–11.3 (Table 1), this being indicative of calcium carbonate (CaCO₃) supersaturation and probably calcite precipitation. Only T and the SRP concentration showed temporal variation along the diel cycle (KW test, p < 0.05).

Among the studied lakes, LC showed the highest concentration of N-NO₃ (Table 1), but CM and CG had the highest N-NO₂ and N-NH₄ concentrations, respectively (Table 1). Lakes exhibited similar SRP, SIc, Chl-a, TSS, NEP, and DOC values (Tables 1 and 2). Chl-a and TSS concentrations showed vertical variation along the water column (KW test, p < 0.05).

DIC and dissolved CO₂ concentration

The studied lakes were consistently CO₂ supersaturated relative to the atmospheric equilibrium concentration (∼10.9 µmol L⁻¹) during both seasons, with higher average concentrations in summer than in winter (Fig. 2). There was no significant variation observed in C_CO2 and DIC concentrations throughout the day in any season (KW test, p = 0.6 for DIC and p = 0.2 for C_CO2). C _CO2 and DIC exhibited a weak but significant correlation in summer (Spearman correlation, R² = 0.16, p < 0.05). LC showed the highest concentration of DIC, whilst higher C_CO2 was observed in LC and CG, (Fig. 2). No significant vertical variation of DIC or C_CO2 was observed along the water column in the lakes (KW test, p = 0.06 for DIC and p = 0.8 for C_CO2).

PCAs revealed that C_CO2 was positively related to DO as opposed to photosynthesis-respiratory cycles (Fig. 3). In summer, a PCA explaining 63% of the observed variance showed DIC in the first axis together with N-NO₂, T, N-NO₃, and DOC. In winter, DIC also appeared in the first axis of the PCA, but related to K₂₅, and Chl-a (Fig. 3). The explained variance of this PCA was 60% in winter. We further identified predictors of C_CO2 by using MRM (Table 3). During the summer, C_CO2 in the water column was negatively related to T, pH, DO, and N-NO₃ which accounted for 80% of the C_CO2variability. T was the most significant variable in the summer model (Table 3). During winter, the MRM for C_CO2 included N-NO₂, pH, and log₁₀(TSS) inversely as independent variables, which together explained 70% of the C_CO2 variability. N-NO₂ was the most significant variable in the winter model (Table 3).

CO₂ evasion fluxes (F_CO2)

CO₂ was consistently emitted from the Ruidera Lakes into the atmosphere during the day in both seasons, with an estimated mean annual emission rate of 3.64 ±2.59 g CO₂ m ⁻² d⁻¹. Higher average seasonal F_CO2 was observed in summer, and the post-hoc test confirmed this pattern at LC and CM lakes (Fig. 4). During both seasons, F_CO2 values were higher in the evening and at night compared to the morning in LC and CM in summer, and LC and CG in winter (Fig. 4). Highest F_CO2 values were measured in LC and CG during summer and winter (Fig. 4).

Using all data, there was a positive and moderate correlation between the mean C_CO2 and F_CO2 (Pearson correlation: F_CO2 = 8.2 log₁₀(C _CO2)–9.8, p < 0.05, R² = 0.43; Fig. 5). However, these relationships were not significant when testing seasonally (summer: R² = 0.32, p > 0.05, and winter: R² = 0.42, p > 0.05).

δ¹³ C-CO₂ values and sources

There were no seasonal differences in δ¹³C-DIC, δ¹³C-CO₂, and δ¹³C-OM (Table 4). δ¹³CO₂₋fract did exhibit significant seasonal differences supported by the contrasted values recorded in LC-winter and CM-summer (Table 4). δ¹³C-DIC values significantly differed between LC and CM during summer (Table 4). δ¹³C-DIC and δ¹³CO₂became enriched downstream during both seasons (Table 4). Finally, δ¹³C-OM was similar seasonally (Table 4).

δ¹³C-CO₂ values were uncorrelated with the δ¹³C-OM signatures (Fig. 6A). Contrarily, δ¹³C-CO₂ did depend on δ¹³C-DIC, most strongly during summer (Fig. 6B). Computed signatures considering the fractionation of the carbon isotopes during calcite precipitation (δ¹³C-CO₂-fract) also displayed a strong correlation with observed values, mainly in summer (Fig. 6B).

The keeling plot for DIC showed a significant relationship in winter: the increase in DIC concentration led to δ¹³C-DIC enriched values; in summer, the increase in DIC concentration did not affect δ¹³C-DIC signatures (Fig. 7A). Contrarily, the increase in C_CO2 resulted in more negative δ¹³C-CO₂ values in both seasons, although less so in summer; however, both estimated regressions have very similar slopes for the two seasons (Fig. 7A). There were significant relationships in the Miller-Tans plot for both seasons and the whole annual data set (Fig. 7B). The linear regression models for DIC were found to be seasonally distinguishable, while those for CO₂were not. The δ¹³C-DIC source values, approximated from the slope of the Miller-Tans regression models, showed strong seasonal differences, a more positive δ¹³C-DIC influence in winter (−5.9‰) than in summer (−14.5‰; Fig. 7B). The δ¹³C-CO₂ source values resulted very similar between seasons, consistent with the average δ¹³C-CO₂ calculated from the fractionation of calcite precipitation (–18.96‰  vs. −18.47‰, respectively; Fig. 7B). Finally, the isotope effect observed on δ¹³C-DIC values as CO₂ degassing downstream in the lake chain was well represented in the model depicted in Fig. 7C for the summer season; however, the effect in winter was not conclusive.