A Forty-Year Karstic Critical Zone Survey (Baget Catchment, Pyrenees-France): Lithologic and Hydroclimatic Controls on Seasonal and Inter-Annual Variations of Stream Water Chemical Composition, pCO₂, and Carbonate Equilibrium

Abstract

The long-term trends and seasonal patterns of stream water chemical composition in a small remote forested karst catchment, were investigated from 1978 to 2018. Calcium, magnesium, and bicarbonates, the dominant ions, increased over the period together with temperature, while sulfates decreased. Carbonate and sulfate mineral dissolution was the main source of these elements. These trends and the seasonal opposite patterns of discharge vs. temperature, calcite saturation index vs. pCO₂ and bicarbonate vs. sulfates, suggested the influence of discharge, of reduced long-range atmospheric pollution, and of increasing air temperature on biological activity and carbonate dissolution. Furthermore, the hydrological regime controlled the seasonal stream water chemical composition and fluxes by: (i) a dilution during the high discharge period, (ii) a change in the contribution rate of the waters draining different lithological areas in the catchment, e.g., the increased sulfates to bicarbonates ratio during summer low flows, with a maximum alkalinity decrease of 24%, and (iii) a “piston” and a “flushing” effects of dissolved elements stored in soils and epikarst with the first autumn heavy rains. Long-term stream water hydrochemical surveys of karst system have proved to be powerful indicators of biogeochemical processes, water sources and pathways under variable natural and anthropogenic environmental pressure conditions.

1. Introduction

The increasing trend of CO₂ has a feedback effect on climate change and noticeably on the frequency of extreme hydro-climatic events [1,2,3]. The atmospheric CO₂ sink and the transfer of dissolved inorganic carbon (DIC) from the atmosphere to the critical zone and to the oceans, via the rivers is greatly controlled by carbonate dissolution [4,5]. Although carbonates constitute between 10–15% of the continental surfaces [6,7,8,9], their chemical weathering has an important role: (i) in the annual production of dissolved elements exported (between 45% and 60%) annually by rivers to the oceans [5,10] and (ii) in the consumption of atmospheric/soil CO₂ (0.38 Gt·yr⁻¹) at the global scale [7]. Even though carbonate dissolution is the main part of the karstification process, it is not considered on a global scale as a carbon sink over geological timescales [11] because carbonate precipitation in the oceans (source of CO₂) compensates carbonate dissolution on the continents (sink of CO₂). However, its role on regulating global atmospheric CO₂ becomes relevant over short timescales. This is due to the quick kinetic dissolution [12] getting to the control of DIC concentrations in almost all carbonated catchment areas [13,14,15,16,17].

In addition, karst aquifers provide almost 25% of the water used for drinking water supply [8] and about half of the world’s karst area lies between latitudes 20° and 40° north [7]. Furthermore, mountainous areas in southern Europe appear to be particularly vulnerable and sensitive to climate change [18,19,20,21,22] and also to land cover changes [23,24]. This is the case of the Pyrenees Mountains (South of France), which are expected to be severely affected by environmental perturbations during the 21st century, mainly with the decrease in water availability [25,26].

The karst system is mainly recharged by rainwater, where water seepage is enabled by the fractured and karstic nature of the calcareous rock [27]. The karst system includes both the infiltration zone (including epikarst) and the saturated zone, and by extension, the hydrological system that flows in the fractures and in the conduits formed mainly by the dissolution of carbonate rock [28]. Epikarst is the quasi-permanent shallow and discontinuous saturated layer, where the infiltrated water dissolves the carbonate mineral and changes the geochemical features of the rock [29]. Indeed, the main agent of the dissolution of carbonates is the biogenic CO₂, supplied by the organic matter mineralization and stored in the soil, a few meters above the epikarst [30]. However, the CO₂–H₂O–Ca_(1−x)Mg_(x)CO₃ system is sensitive to environmental changes due to karstic characteristics such as high transmissivity (lower residence time), high dissolution kinetics, and low filtering role of the infiltration zone [31].

Among the parameters that play a role in the dissolution of carbonates, we can highlight hydrology [29,32], climate [17,33], land-use/vegetation [34,35], lithology [7,36,37], and human activities [38,39].

Several studies have shown the strong control that runoff has over carbonate dissolution flow [7,40]. Moreover, during the high flow period, the huge and rapid water flow through the karstic network could reduce the dissolution due to less water-rock interaction [41] but on the opposite, enhance the dissolution by renewal of the solutions in contact with the mineral [42]. Then, it may also explain a gradual increase in element concentrations as runoff increases.

Temperature is another parameter which strongly influences carbonate dissolution. Thermodynamically, an increase in temperature leads to a decrease in carbonate solubility [43,44,45]. However, increasing temperatures also enhance biological activity and decomposition of soil organic matter, causing an increasing pCO₂, carbonic acid, and organic acid contents in the critical zone, particularly in epikarsts and caves, accelerating carbonate dissolution [34,46,47,48]. In addition, based on a global database of pure carbonate lithologies, Probst [49], Amiotte-Suchet et al. [7], and Gaillardet et al. [50] evidenced a control of temperature over carbonate weathering intensity, with a Ca²⁺ + Mg²⁺ concentration peak in river water under temperate conditions.

The CO₂ mainly originates from atmospheric and soil sources. Vegetation plays a key role in sequestering and storing atmospheric CO₂ to be transformed into organic matter [51]. Subsequently, the decomposition of organic matter and root respiration lead to the production of CO₂ in soils [52]. Between these two CO₂ sources, the direct input of atmospheric CO₂ to the soils is usually low due to high pressure of soil CO₂ and low pH values (lower than 5.6) of rainwater [53]. The dissolution of CO₂, mainly resulting from the oxidation of organic matter, is illustrated by the following chemical reactions:

$${\mathrm{CH}}_2\mathrm{O}+{\mathrm{O}}_2\leftrightarrow {\mathrm{CO}}_2+{\mathrm{H}}_2\mathrm{O},$$

(1)

$${\mathrm{CO}}_2+{\mathrm{H}}_2\mathrm{O}\leftrightarrow {\mathrm{H}}_2{\mathrm{CO}}_3,$$

(2)

$${\mathrm{H}}_2{\mathrm{CO}}_3\leftrightarrow {\mathrm{H}}^{+}+{\mathrm{HCO}}_3^{-},$$

(3)

Furthermore, in recent years, increases of pCO₂ in soils and caves [33,54], total dissolved solids (TDS) and alkalinity in rivers have been observed [55,56,57]. In the meantime, the reduction of annual streamflows, particularly in Mediterranean rivers was evidenced [58,59,60]. These changes suggest the intervention of external driving factors such as climate change [61,62,63] and/or land use change [24,38,64,65].

Carbonate dissolution can also be carried out by chemical agents other than carbonic acid. Anthropic disturbances such as agricultural practices using N-fertilizers [66,67,68], atmospheric deposition from long-distance sources of air pollution [69,70,71,72] and acid mine drainage [73,74], are potential sources of sulfuric and nitric acids that can add to or substitute carbonic acid in carbonate dissolution process. Likewise, there are also natural sources of sulfuric acid from volcanoes and hot springs as well as sulfur minerals in soils and sediments [75,76]. Additionally, by enhancing the capture of dry deposition, coniferous forests can contribute to the increase in sulfuric acid inputs to soils by atmospheric deposition (dry and wet deposition), particularly from throughfalls [77,78] at the catchment scale. Unlike carbonate dissolution by carbonic acid that generates two moles of bicarbonates consuming one mole of CO₂, sulfuric and nitric acids only release one mole of bicarbonates, which could escape from the water in the form of gas (CO₂ gas), or could enter seepage water and dissolves carbonate again (aqueous CO₂) [39].

In this context, the Baget catchment in the Pyrenees Mountains, was selected as a research area to analyze the fluctuations of physico-chemical composition of stream waters over a 40-years period. This catchment has been already extensively studied using time series data for a hydrological perspective [79,80,81] and for some chemical parameters until 2005 [82,83]. However, this is the first time that all the hydrogeochemical data available between 1978 and 2018 were investigated and the pCO₂ in stream water was reconstituted at both interannual and seasonal scales. Then, long-term trends as well as interannual seasonal variations for all the physico-chemical parameters allow to better understand the karst biogeochemical functioning.

The objectives were: (i) to investigate the variation of discharge, temperature, CO₂, SIc (Calcite Saturation Index), and major element concentrations in the river water using long-term surveys at various time scales (instantaneous, monthly, and annual values); (ii) to depict the impact of global changes on carbonate dissolution and hydrochemistry of stream waters; (iii) to elucidate the key factors controlling carbonate dissolution (such as temperature, discharge, and vegetation), and the respective contributions of water sources, particularly from karst and epikarst; and iv) to understand the dynamics of CO₂, carbonate dissolution, and stream water chemistry.

According to described uplines, several hypotheses were set up to explain the observed trends. At the year scale, the discharge would cause an increase in carbonate dissolution, but with a dilution effect during the highest flow conditions. The increase in air temperature should finally enhance indirectly carbonate dissolution due to the rise in organic matter mineralization and production of carbonic and organic acids, even if carbonate solubility would decrease with temperature increase [82]. Therefore, an increase in global air temperature and change in land use, i.e., from crops or grasslands to forest, should increase carbonate dissolution on the long-term. However, other external anthropic pressures might add other drivers to the karstic system and change the dynamic of carbonate dissolution.

2. Materials and Methods

2.1. Site Description

The small Baget catchment (BC, 13.25 km²; 42°57′42″ N 0°58′30″ E; Figure 1) is located in the lower part of the Pyrenees Mountains, in the Ariege department, in the Southwest of France. Its altitude ranges from 498 to 1417 m.s.l., with a water flow orientation from West to East [84]. Slopes are very steep, particularly in the middle part of the basin, where they reach 60% [79].

BC drains a karst area, mainly in the calcareous formation of the Jurassic and Lower Cretaceous, relatively homogeneous at the outcrop and with a mineralogical composition usually greater than 98% of calcite. However, they are associated with pyrite, which is sometimes abundant in the area of contact with the impermeable formation of black flysch (Figure 1b) [83]. Some schists (with small pyrite impurities), minor gypsum rocks and dolomitic rocks were also observed in the BC [85].

The catchment is essentially forested. The land cover is dominated by fir-beech forests and it is particularly dense on the north-facing slope, whereas, grassland occupies partly the upper part and the lower part close to the outlet of the south-facing slope. Although there is an extensive vegetation, the soils are generally very thin. The limestones are covered by discontinuous brown soils, sometimes with rendzina soils, generally occupying the limestone pavement slits. Podzols are developed on detrital formations. The grasslands are located on podzolic soils well exposed to the south near to the Lachein stream water and also, in the bottom of the valley, downstream, on the alluvium [83]. BC is weakly exposed to local anthropogenic pollution with only few remaining mountain stockbreeding on the upper part of the south facing slope.

BC is under the influence of the Atlantic oceanic climate characterized by a clear mountain tendency with an average annual air temperature of 12.3 °C, an average annual rainfall close to 1700 mm and an estimated annual evapotranspiration of 540 mm [79,87,88]. Those values are close to the daily average air temperature of 12.0 ± 6.3 °C and higher than the annual average rainfall of 975.5 mm, registered at the closest Meteo-France meteorological station of St. Girons (43°00′19″ N 1°06′25″ E, elevation 414 m), located 8.3 km downstream from the outlet of BC. The precipitation gradient between the two stations is of +138 mm/100 m, which is close to that reported for the Jura Mountains [34], but much greater than that mentioned by Joly [89] for the Pyrenees Mountains (>25 mm/100 m). Indeed, this gradient depends on factors such as wind direction, relief, distance to the ocean and latitude, highly variable in mountainous areas [90].

The runoff which drains calcareous lithology can be lost through seepage or in losses of the Lachein stream water (e.g., losses of “La Hille” and “La Peyrère”, this last one performing as a spring during flood period, Figure 1a), thus feeding the Baget groundwater karstic network. Downstream of the Baget valley, the Lachein stream water (usually dried up in summer) connects with the main perennial spring of the Baget karst system called “Las Hountas” (Figure 1a) [79]. The high-water flow period (from November to April) displays a bimodal rainfall regime with usually a first peak in December and a second one in February [79]. Considering the mean altitude, the storage of water in the form of snow is low and snowmelt does not result in significant increases in discharge. The average annual discharge is 0.48 m³·s⁻¹ with a specific discharge of 36 L·s⁻¹·km⁻² [79,91], as measured at the outlet station (Baget B1, Figure 1a).

This experimental catchment belongs to the French Karst Network (SNO Karst, Jourde et al., 2018) [92] and is one of the observatories of the French Research Infrastructure called OZCAR (Observatoire de la Zone Critique: Application et Recherche) [93], which is part of the European Research Infrastructure eLTER (Long Term Ecosystem Research in Europe) [94,95].

2.2. Sampling and Analytical Methods

The hydrochemical monitoring station is located at the outlet of the catchment (Baget B1, Figure 1a). The water level was measured continuously using a mechanical limnigraph (OTT 20 1/5) and a float-type water level sensor (OTT Thalimedes) [96]. Discharge is estimated from the rating curve calibrated for this gauging station [79]. For water quality purpose, we used, the samples collected weekly from 1978 to 2006 by the CNRS “Laboratoire Souterrain de Moulis” [82] and every six months from 2007 to 2014 by the Adour-Garonne Water Agency (accessible via the ADES database, BRGM) [96]. By 2014, water sampling and hydrochemical parameters analyses were performed by EcoLab (Laboratoire Ecologie fonctionnelle et Environnement in Toulouse), fortnightly and more frequently during flood events using an automatic sampler. It should be noted that in 1990, a long-term pumping test dried the Baget spring (“Las Hountas”), so no data was available for part of this year.

From 1978 to 2006, pH, total carbonate alkalinity (in mmol·L⁻¹; 0.2 TAC = HCO₃⁻) and total water hardness (in mmol·L⁻¹; 0.1 TWH = Ca²⁺ + Mg²⁺, R² = 0.92, N = 541) were measured using a titrimetric method. Major ions were analyzed by ion chromatography.

Since April 2014, pH, pressure (water level), conductivity, water temperature, turbidity, dissolved O₂, N-NO₃⁻, and Cl⁻ concentrations were registered at a high frequency (every 10 mn) in the mid-depth of the water column (at Baget B1) using a multiparameter probe (YSI 6920V2-01, YSI Incorporated, Yellow Springs, OH, USA) equipped with different sensors and an atmospheric pressure correction. Moreover, from 2017 to present-day, a portable multi-parameter (WTW probe, Xylem Analytics Germany Sales GmbH & Co., Weilheim, Germany) was used to measure physicochemical parameters of water in the field, including water temperature (WT), pH, and electrical conductivity (EC) at the time of sampling.

Samples were collected in polypropylene plastic bottles of 1 L, then filtered immediately back to field in the laboratory using a Millipore 0.22 µm cellulose nitrate membrane. Since 2014, the protocol for sample treatment in the laboratory is as follows: firstly, an aliquot of 125 mL was taken to measure cation concentrations. This sub-sample was acidified with 3 drops of 16N HNO₃ (to prevent complexation and precipitation) and stored until the analyze by inductively coupled plasma optical emission spectrometer (ICP-OES; Iris Intrepid II XLD, Thermo Electron, Thermo Fisher Scientific, Waltham, MA, USA); secondly, an aliquot of 2 mL was used to analyze major anion concentrations (SO₄²⁻, NO₃⁻, and Cl⁻) by ion chromatography (Dionex apparatus ICS 5000⁺, Thermo Fisher Scientific, Waltham, MA, USA); lastly, an aliquot of 20 mL was taken to measure alkalinity concentration using standard acid titration method with HCl 0.02N and a Metrohm titrant (716 DMS Titrino, Metrohm, Riverview, Forida, USA) within 12 h after sampling. DOC was analyzed on a Shimadzu TOC 5000 analyzer (Shimadzu Corporation, Kyoto, Japan) using the 680 °C combustion catalytic oxidation method and CO₂ detection using an infrared gas analyzer (NDIR). All samples were stored in darkness at 5 °C before analysis.

Concerning the quality assurance (QA) and the quality control (QC), the data set is too long to have a QA/QC procedure over the whole period from 1978 until 2018. Anyway, concerning the QA, the protocols for sampling, filtration, pre-treatment and storage of the stream water before analysis were carefully controlled during the whole period. Concerning the QC, all the methods of analysis are normalized following French (NF), European (EN), and international (ISO) standards. The accuracy of analyses was controlled using blanks and two Certified Reference Materials (ION-96.4 and ION SUPER 05) depending on the range of ion concentration. The detection limits for major ions range between 0.2 and 65 μmol·L⁻¹ according to the elements.

For each sample analysis, the Net Inorganic Charge Balance (NICB % = 100 × (TZ⁺ − TZ⁻)/(TZ⁺ + TZ⁻)) between the sum of cations (TZ⁺) and the sum of anions (TZ⁻) (considered both in meq·L⁻¹) was calculated and was inferior to 5% for almost all samples, indicating that there is a good equilibrium between cation and anion charges, and consequently that the contribution of organic anions to the anionic charge remain negligible. This is confirmed by the low DOC concentrations (0.82 ± 0.79 mg·L⁻¹, N = 39, from November 2016 to September 2018).

Furthermore, since several decades, the physico-chemical analysis laboratory of EcoLab participates each year to the ICP-Water Intercomparison with the Norwegian Institute for Water Research (NIVA), which follows a QA/QC procedure.

2.3. Data Treatment

2.3.1. pCO₂ and Calcite Saturation Index (SIc)

The partial pressure of carbon dioxide (pCO₂) and calcite saturation index (SIc) were calculated using the mass equation. The calculations were computed with PhreeqC Interactive 3.4 software and the database “phreeqc.dat” [44], using the analyzed major element concentrations in stream water as well as the physicochemical parameters measured in the field during the water sampling.

In addition, stream water temperature was necessary to calculate pCO₂. The air temperature available for the full period and the water temperature available from 2014 to 2018 (N = 73) confirmed that the water temperature in stream varied very slightly during the year (from 9.5 to 11 °C). It could thus be estimated using the air temperature in centigrade degrees (T_water = 0.03 T_air + 9.93, R = 0.6, N = 59, p-value < 0.0001).

Due to the absence of concentration data for some major ions (SO₄²⁻; Cl⁻; NO₃⁻, Na⁺, and K⁺) during certain time slots, pCO₂ was also calculated manually, assuming the dissolution of pure calcite and using the following chemical reaction [97]:

$${\mathsf{\alpha }}_1{\mathrm{CaCO}}_{3\left(\mathrm{s}\right)}+{\mathsf{\alpha }}_2{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}+{\mathsf{\alpha }}_2{\mathrm{CO}}_{2\left(\mathrm{g}\right)}\leftrightarrow {\mathsf{\alpha }}_1{\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{2+}+{\mathsf{\alpha }}_1{\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-}+{\mathsf{\alpha }}_2{\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-},$$

(4)

The coefficients α₁ and α₂ had the same value but showed that, when the dissolution of carbonates by carbonic acid, one bicarbonate ion came from the mineral dissolution and the other one came from the consumption of atmospheric/soil CO₂.

$${\mathrm{CaCO}}_{3\left(\mathrm{s}\right)}\leftrightarrow {\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{2+}+{\mathrm{CO}}_{3\left(\mathrm{aq}\right)}^{2-};{\mathrm{K}}_{\mathrm{s}}=\left({\mathrm{Ca}}^{2+}\right)·\left({\mathrm{CO}}_3^{2-}\right),$$

(5)

$${\mathrm{CO}}_{2\left(\mathrm{g}\right)}+{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}={\mathrm{H}}_2{\mathrm{CO}}_{3\left(\mathrm{aq}\right)}^{};{\mathrm{K}}_{\mathrm{H}}=\frac{\left({\mathrm{H}}_2{\mathrm{CO}}_3^{}\right)}{\left({\mathrm{pCO}}_2\right)},$$

(6)

$${\mathrm{H}}_2{\mathrm{CO}}_{3\left(\mathrm{aq}\right)}^{}={\mathrm{H}}_{\left(\mathrm{aq}\right)}^{+}+{\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-};{\mathrm{K}}_1=\frac{\left({\mathrm{H}}^{+}\right)·\left({\mathrm{HCO}}_3^{-}\right)}{\left({\mathrm{H}}_2{\mathrm{CO}}_3^{}\right)},$$

(7)

$${\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-}={\mathrm{H}}_{\left(\mathrm{aq}\right)}^{+}+{\mathrm{CO}}_{3\left(\mathrm{aq}\right)}^{2-};{\mathrm{K}}_2=\frac{\left({\mathrm{H}}^{+}\right)·\left({\mathrm{CO}}_3^{2-}\right)}{\left({\mathrm{HCO}}_3^{-}\right)},$$

(8)

$${\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}={\mathrm{H}}_{\left(\mathrm{aq}\right)}^{+}+{\mathrm{OH}}_{\left(\mathrm{aq}\right)}^{-};{\mathrm{K}}_{\mathrm{w}}=\left({\mathrm{H}}^{+}\right)·\left({\mathrm{OH}}^{-}\right),$$

(9)

where the round parenthesis represents the activities of the respective species, and K, the equilibrium constants calculated using the temperature dependence of thermodynamic constants defined by Plummer and Busenberg [98].

Assuming that the dissolution of carbonates took place under open system conditions and in equilibrium with regard to calcite, the Ca²⁺ or pCO₂ could be estimated using the following equation:

$$({\mathrm{Ca}}^{2+})={\left(\frac{{\mathrm{K}}_1{\mathrm{K}}_{\mathrm{H}}{\mathrm{K}}_{\mathrm{S}}}{4{\mathrm{K}}_2{\mathsf{\gamma }}_1^2{\mathsf{\gamma }}_2}{\mathrm{pCO}}_2\right)}^{\frac{1}{3}},$$

(10)

where γ₁ and γ₂, the ion activity coefficients for single-charge and dual-charge species, respectively. These activity coefficients were calculated using the Davies equation [99] with the PhreeqC Interactive 3.4 software [44].

These pCO₂ values calculated using the PhreeqC Interactive 3.4 software and manually calculations for the periods 1994–2006 and 2016–2018 confirmed that both results are consistent (R² = 1.00).

Likewise, the calcite saturation indexes were calculated using the thermodynamic constant K_s (solubility product of calcite) at given water temperature, as given below by Clark and Fritz [97]:

$$\mathrm{SIc}=\log \left(\left({\mathrm{Ca}}^{2+}\right)·\left({\mathrm{CO}}_3^{2-}\right)\right)/{\mathrm{K}}_{\mathrm{s}},$$

(11)

2.3.2. Statistics

All calculations were computed on a hydrological year basis (October–September). For simplicity, only the second year of each period will be indicated, for example the period from October 1978 to September 1979 will be shown as the year 1979. The mean concentrations were discharge-weighted. Statistical data treatments like correlation matrix and principal component analysis (PCA) were done using RStudio software (CRAN, the Comprehensive R Archive Network, version 3.5.1, RStudio, Inc., Boston, MA, USA). PCA was performed on raw data transformed with log-ratio data using the “rgr” package. The detection of the change-points was calculated using the Buishand test [100,101]. Due to non-normal distribution of the data: (i) trends were determined using the Mann Kendall trend test on linear regression and (ii) the degree of association between two variables was evaluated by the Spearman’s correlation coefficient. Plotting for ternary diagram (Piper diagram, [102]) was performed using the method developed by Hamilton [103].

3. Results

3.1. Stream Water Physico-Chemical Characteristics

3.1.1. Hydrochemical Composition of the Baget Stream Water

The hydrochemical properties are reported in

Table 1. Comparison of the main physico-chemical parameter characteristics of the Baget stream water for the two periods: 1978–2018 (first number) versus 1994–2006 plus 2016–2018 (second number), both periods separated by a slash, respectively. Max., maximum; Min., minimum, Std. Dev., standard deviation; CV, coefficient of variation $\left(\mathsf{\sigma }/\overline{\mathrm{X}}\right)$. The units are m³·s⁻¹ (Q, discharge), °C (T°, air temperature), µS·cm⁻¹ (Cond, conductivity), meq·L⁻¹ (major elements), mg·L⁻¹ (TDS) and atm (pCO₂). The mean concentrations are discharge-weighted. SIc is the calcite saturation index.

Parameter	N	Max.	Min.	Mean $\left(\overline{\mathbf{X}}\right)$	Median	Mode	Std. Dev. (σ)	CV
Q	14609/592	10.10/5.41	0.02/0.04	0.44/0.39	0.22/0.20	0.13/0.15	0.67/0.54	1.52/1.40
Air T°	14604/592	29.5/28.2	−11.1/−5.0	12.0/12.2	11.9/12.0	9.2/8.4	6.33/6.35	0.53/0.52
pH	1506/592	8.4/8.4	6.5/6.9	7.7/7.7	7.7/7.7	7.6/7.7	0.18/0.18	0.02/0.02
Cond	637/381	398/398	223/268	306/317	307/316	319/319	24.90/18.97	0.08/0.06
Ca2+	1504/592	3.65/3.55	2.40/2.57	2.97/3.03	2.96/3.01	3.00/3.15	0.18/0.17	0.06/0.05
Mg2+	1504/592	0.75/0.75	0.23/0.23	0.36/0.37	0.36/0.36	0.34/0.40	0.05/0.06	0.13/0.15
Na+	880/592	0.07/0.06	0.03/0.03	0.05/0.05	0.05/0.05	0.05/0.05	0.01/0.00	0.11/0.10
K+	880/592	0.19/0.04	0.00/0.00	0.01/0.01	0.01/0.01	0.01/0.01	0.01/0.00	0.57/0.28
HCO3−	1501/592	4.40/3.66	2.43/2.57	2.90/2.98	2.88/2.96	2.76/2.76	0.22/0.19	0.07/0.06
SO42−	763/592	0.84/0.78	0.01/0.04	0.35/0.35	0.34/0.34	0.27/0.27	0.13/0.13	0.37/0.36
Cl−	893/592	0.25/0.13	0.00/0.02	0.05/0.05	0.05/0.05	0.04/0.04	0.01/0.01	0.27/0.23
NO3−	705/592	0.17/0.17	0.00/0.00	0.03/0.03	0.03/0.03	0.03/0.03	0.02/0.02	0.57/0.56
TDS	642/592	331.5/327.5	229/229	272.8/273.6	272.0/272.8	264.3/264.3	15.02/14.41	0.06/0.05
pCO2	1498/592	1 × 10−1.64 /1 × 10−1.68	1 × 10−3.23 /1 × 10−3.23	1 × 10−2.45 /1 × 10−2.48	1 × 10−2.50 /1 × 10−2.52	1 × 10−2.44 /1 × 10−2.50	0.002/0.001	0.53/0.46
SIc	1482/592	0.78/0.78	−0.88/0.64	0.03/0.09	0.03/0.07	0.03/0.04	0.18/0.18	6.09/2.12

. All available data from 1978 to 2018 were considered in the Section 3.1.1 and Section 3.1 whereas the others sub-sections only concerned the dataset from the period October 1994 to September 2006, plus October 2016 to September 2018 (N = 592). This allowed to homogenize the periods of time studied since the parameters had different sampling frequencies. The values show averages with a 95% confidence interval (±2σ).

To ensure the consistency of common parameters, the averages were compared for the two studied periods (Table 1). It could be noted that the average and median values for the physico-chemical parameters were similar, with an average variation always less than 5% for almost all the considered parameters. The exception was mainly for discharge (Q) whose mean value decreased from 0.44 to 0.39 m³·s⁻¹ and the calcite saturation index (SIc) which increased from 0.03 to 0.09, comparing the whole period (1978–2018) to the 1994–2006 plus 2016–2018 period, respectively.

As a whole, Q exhibited a high variability (coefficient of variation $\left(\mathsf{\sigma }/\overline{\mathrm{X}}\right)$, CV = 1.52) compared to the variability in air temperature (T°, CV = 0.53), sulfate (SO₄²⁻, CV = 0.37), partial pressure of CO₂ (pCO₂, CV = 0.53), magnesium (Mg²⁺, CV = 0.13), alkalinity (HCO₃⁻, CV = 0.07), and calcium (Ca²⁺, CV = 0.06).

Over the whole period, Q and T° values registered daily in BC were in the range of 0.02 to 10.10 m³·s⁻¹ and −11.1 to 29.5 °C, respectively. Stream water pH averages 7.7 ± 0.4 and the conductivity varied from 223 to 398 µS·cm⁻¹.

The total dissolved solid content (TDS = (Ca²⁺) + (Mg²⁺) + (Na⁺) + (K⁺) + (HCO₃⁻) + (SO₄²⁻) + (Cl⁻) + (NO₃⁻) + (H₄SiO₄)) varies from 229 to 331.5 mg·L⁻¹, with a mean value of 273 ± 30 mg·L⁻¹. The relative abundance of each ion or group of ions (ex: monovalent or divalent cations) can be expressed as the percentage of the total respective sum of cations or anions in equivalents per liter (Figure S1). Whatever the season and the hydrological conditions, the stream water chemistry was dominated by calcium and alkalinity (mainly corresponding to bicarbonate), which are controlled by carbonate dissolution as follows:

$${\mathrm{Ca}}_{1-\mathrm{x}}{\mathrm{Mg}}_{\mathrm{x}}{\mathrm{CO}}_3+{\mathrm{H}}_2{\mathrm{CO}}_3{}_{}\to \left(1-\mathrm{x}\right){\mathrm{Ca}}^{2+}+\left(\mathrm{x}\right){\mathrm{Mg}}^{2+}+{\mathrm{HCO}}_3^{-}+{\mathrm{HCO}}_3^{-},$$

(12)

Calcium was the most abundant cation with a concentration varying from 2.40 to 3.65 meq·L⁻¹, i.e., 80.1–92.0% of the sum of cations. The second most important one was Mg²⁺ with a mean value of 0.36 ± 0.10 meq·L⁻¹.

Then, the mean Ca²⁺/Mg²⁺ equivalent ratio was 8.4 ± 1.8. This ratio in Baget was higher than in the Beipanjiang River (3.6) from the Asian Karts Region [70], in the Herault River (2.5) from southern France [104], in the Alpine region (3.7 ± 2.8) [105], in the Ibrahim River–Lebanon (5.8) [106], and exceeded the average value of the world’s carbonate basins (4.4) [4]. Nevertheless, it was lower than in the springs of the karstic area of the Jura Mountains (22) in France [34] or in the La Sigouste River (54 ± 39) from the French Southern Alps [107]. In any case, it was rather consistent with the mean value proposed by Meybeck [4] for the French monolithologic basins draining carbonate rocks (12).

Likewise, HCO₃⁻ was the most abundant anion (75.5–99.6% of total anions) with a mean concentration of 2.90 ± 0.42 meq·L⁻¹, while SO₄²⁻ varies from 0.01 to 0.84 meq·L⁻¹. In addition, the equivalent ratio of Ca²⁺ + Mg²⁺ to HCO₃⁻ + SO₄²⁻ was close to a constant value of 1 (1.02 ± 0.06 N = 734). The SO₄²⁻/HCO₃⁻ equivalent ratio was more variable and ranged from 0 to 0.29 and exhibits a mean value of 0.12 ± 0.10 (N = 735), indicating during some hydroclimatic conditions that the stream water was slightly enriched in SO₄²⁻.

Concentrations of sodium and potassium were much lower, with mean values of 0.05 ± 0.01 and 0.01 ± 0.00 meq·L⁻¹, respectively. Chloride and nitrate ions exhibited also low concentrations, with mean values of 0.05 ± 0.01 and 0.03 ± 0.02 meq·L⁻¹, respectively. The respective contribution of these strong acid anions did not exceed 6.5% of the total anions.

The pCO₂ in stream water, which was mainly a function of organic matter degradation and soil respiration, varied from 1 × 10^−3.23 to 1 × 10^−1.64 atm with a mean of 1 × 10^−2.45 atm. The highest pCO₂ values are during low water flow period from June to October. Stream water pCO₂ was always higher than the atmospheric pCO₂ (1 × 10^−3.5). SIc varied in the range of –0.88 to 0.78, and 57% of the samples had SIc > 0. Moreover, SIc in BC was lower in the dry season than in the wet season.

3.1.2. Multivariate Statistical Analysis

The principal component analysis (PCA, Figure 2) sets the studied variables according to factors explaining their relative variations [108]. The correlation between the variable and the factor was indicated by a factor loading varying from 1 (strong positive correlation) to −1 (strong negative correlation) [109]. In order to improve the representation of the variables in the PCA, it was decided to consider only 11 variables, the least ones autocorrelated and representative of the hydrochemical processes in the BC. Thus, for example, neither pH nor SIc were considered in the PCA, both strongly correlated with pCO₂ (see Equations (10) and (11)), or the TDS which was significantly correlated with Ca²⁺ and HCO₃⁻, these two ions representing more of the 75% of TDS. The results of PCA are presented in Figure 2, Tables S1 and S2.

The percentage of variance in the data explained by the first four axes was 83.6%. The first axis (Table S1), responsible for 42.6% of the variance, contrasted Ca²⁺ (87% of contribution to this axis), air temperature (79%), HCO₃⁻ (74%), Na⁺ (68%), Mg²⁺ (60%), and SO₄²⁻ (30%) on the upper positive side of the PCA biplot, with NO₃⁻ (40%) and discharge Q (17%) positioned in the lower negative side.

The second axis (Tables S1 and S2), which extracted 18.9% of the variance, contrasted Q (76%) and HCO₃⁻ (14%) in the right-hand side with SO₄²⁻ (53%), NO₃⁻ (37%), and Mg²⁺ (16%) positioned in the left-hand side. This distribution is reinforced by the positive correlation between Q and HCO₃⁻ and the negative one between SO₄²⁻ and Mg²⁺ (Table S2), suggesting that during low-water level period, SO₄²⁻ and Mg²⁺ concentrations increased, while HCO₃⁻ decreased. Also, Ca²⁺ with HCO₃⁻ on one side and Mg²⁺ with SO₄²⁻ on the other side, appeared in opposite position relatively to Axis 2 in the PCA plot (Figure 2a, Table S1).

The third axis (Figure 2b and Table S1) accounted for 11.2% of the variance (Table S2). It opposed pCO₂ (55%) and K⁺ (34%) with NO₃⁻ (16%), meanwhile the fourth axis, contributing 10.9% to the variance, contrasted K⁺ (47%) and Cl⁻ (33%) with pCO₂ (24%). Besides, Cl⁻ exhibited a positive correlation with K⁺ and NO₃⁻ (Table S2).

3.1.3. Alkalinity Decrease (∆_Alk)

When the main lithology of a catchment is composed of carbonate rocks, the natural dissolution of carbonate by carbonic acid (see Equation (4)) led to an equivalent ratio Ca²⁺ + Mg²⁺/HCO₃⁻ of 1 since for two equivalent of calcium, two equivalents of bicarbonate were released, one coming from the soil CO₂ and the other one from the calcite. However, if the carbonate is dissolved by strong acids such as H₂SO₄ (see Equation (13)) or HNO₃ (see Equation (14)), all the bicarbonates released in solution originated from the CO₃ of the minerals as follows:

$$2{\mathrm{CaCO}}_{3\left(\mathrm{s}\right)}+{\mathrm{H}}_2{\mathrm{SO}}_4\leftrightarrow 2{\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{2+}+2{\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-}+{\mathrm{SO}}_{4\left(\mathrm{aq}\right)}^{2-},$$

(13)

$${\mathrm{CaCO}}_{3\left(\mathrm{s}\right)}+{\mathrm{HNO}}_3\leftrightarrow {\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{2+}+{\mathrm{HCO}}_{3\left(\mathrm{aq}\right)}^{-}+{\mathrm{NO}}_{3\left(\mathrm{aq}\right)}^{-},$$

(14)

These strong acids can be supplied by atmospheric acid deposition over the catchment or by pyrite oxidation (Equation (15)) if there are any pyrite minerals in the geological formation:

$${\mathrm{FeS}}_{2\left(\mathrm{s}\right)}+15/4{\mathrm{O}}_{2\left(\mathrm{g}\right)}+7/2{\mathrm{H}}_2\mathrm{O}\leftrightarrow \mathrm{Fe}{\left(\mathrm{OH}\right)}_{3\left(\mathrm{s}\right)}^{}+2{\mathrm{SO}}_{4\left(\mathrm{aq}\right)}^{2-}+4{\mathrm{H}}_{\left(\mathrm{aq}\right)}^{+},$$

(15)

Then, in calcareous environments, the sulfuric acid released can be rapidly neutralized by carbonate dissolution [110], following the Equation (13). In that case, the equivalent ratio Ca²⁺ + Mg²⁺/HCO₃⁻ increases to 2 if the carbonate is only dissolved by the sulfuric acid.

Nevertheless, in the Baget catchment, there is a mixing of carbonate dissolution by carbonic acid and by strong acids like sulfuric acid, releasing sulfates into solution. Finally, the sum Ca²⁺ + Mg²⁺ is generally equilibrated by the sum HCO₃⁻ + SO₄²⁻ and thus the equivalent ratio Ca²⁺ + Mg²⁺/HCO₃⁻ + SO₄²⁻ tends to 1.

Consequently, strong acid inputs also result in a decrease of the stream water alkalinity, mainly composed by HCO₃⁻ [72,111]. This acidification process due to strong acids can be estimated by the decrease of alkalinity. Indeed, the amount of Ca²⁺ + Mg²⁺ released by natural dissolution of carbonates by carbonic acid produces a higher alkalinity than that observed (Equation (4)). The decrease of alkalinity (Δ_Alk, expressed in equivalent percentage), can be thus calculated as follows for carbonated systems [67]:

$${\mathsf{\Delta }}_{\mathrm{Alk}}(\%)=\frac{\left({\mathrm{Ca}}^{2+}+{\mathrm{Mg}}^{2+}\right)-{\mathrm{HCO}}_3^{-}}{\left({\mathrm{Ca}}^{2+}+{\mathrm{Mg}}^{2+}\right)}\times 100$$

(16)

As seen in Figure 3a, there is a significant positive relationship between Ca²⁺ + Mg²⁺ and HCO₃⁻ in the stream waters of BC but the slope (equivalent ratio) was slightly greater than 1 (Figure 3a), and the cluster of points is above the theoretical line of carbonate dissolution by carbonic acid.

Then, the lag observed in Figure 3a could be both attributed to a relative decrease of alkalinity and to a relative increase of Ca²⁺ + Mg²⁺ concentration, from theoretical alkalinity and Ca²⁺ + Mg²⁺ values corresponding to Equation (4) (line 1/1 in Figure 3a). Nevertheless, the lag observed in Figure 3a could be also attributed to gypsum dissolution as follows:

$${\mathrm{CaSO}}_{4\left(\mathrm{s}\right)}+2{\mathrm{H}}_2\mathrm{O}\leftrightarrow {\mathrm{Ca}}_{\left(\mathrm{aq}\right)}^{2+}+{\mathrm{SO}}_{4\left(\mathrm{aq}\right)}^{2-}+2{\mathrm{H}}_2\mathrm{O}$$

(17)

In order to detect a seasonal and hydrological influence, instantaneous values were gathered into three groups: from July to October (predominantly warm and dry conditions), from November to February (winter time), and from Mars to June (mainly cold and wet conditions). The summer-autumn group exhibited the highest ratios with a mean slope of 1.18 for the regression (dashed line on Figure 3a; p ≤ 0.001; n = 204).

The alkalinity decrease due to strong acids and to gypsum dissolution was correlated with the SO₄²⁻ to Ca²⁺ + Mg²⁺ equivalent ratio, with an average ∆_Alk of 15.2 ± 7.5% which corresponds to an average ratio of 0.13 ± 0.07 (Figure 3b). Thus, as this ratio increased, the ∆_Alk also increased, reaching 24% when the ratio was close to 0.2. These highest ∆_Alk values were associated with the highest equivalent ratios between Ca²⁺ + Mg²⁺ and HCO₃⁻ (Figure 3a), also displayed mainly between July and October (red points). On the opposite, the lowest ∆_Alk values were observed during the winter (blue points) and spring (green points) periods.

3.2. Long Term Trends in Stream Water Chemistry: Instantaneous Values and Mean Annual Data

Long-term trends of the dissolved elements at the outlet of BC for both instantaneous (Figure 4) and mean annual values (Figure 5) were considered to elucidate their behavior over time and their control by environmental factors.

Over the whole 40 years, the air temperature, Ca²⁺ + Mg²⁺ and HCO₃⁻ trends slightly increased, meanwhile the discharge and the pCO₂ trends are not significant (Table S3). The moving averages, which allow to smooth the seasonal variations [112], exhibit the same trends. The two statistical tests applied on these trends (Mann Kendall trend test and linear regression of instantaneous values) were only significant (p < 0.01) for temperature, Ca²⁺ + Mg²⁺ and HCO₃⁻ (Table S4).

During the period 1978–2006, similar increasing trends were observed for air temperature, Ca²⁺ + Mg²⁺ and HCO₃⁻. On the opposite, discharge and SO₄²⁻ concentrations decreased (Table S3). Similarly, as for the longest time period, these trends were only significant (p < 0.01) for temperature, discharge, Ca²⁺ + Mg²⁺, and HCO₃⁻ (Table S4).

The Buishand test performed on instantaneous values evidenced the common ruptures in the time series (Table S4). Concerning Ca²⁺ + Mg²⁺, a breaking point was observed in September 1988, just following a breaking point for discharge in July 1988, at the beginning of the driest period of the whole series: 1988–1989 (0.22 m³·s⁻¹, deficit of 49% compared to the inter-annual average) and 1989–1990 (0.24 m³·s⁻¹, deficit of 45% compared to the inter-annual average). Regarding HCO₃⁻, the breaking point was in October 1991 and occurred six months after the breaking point of pCO₂ in the same year. All these breakpoints are in the period of the pumping test carried out in 1990.

Looking at the mean annual values over the study period, as shown by the deviation to the interannual mean (Figure 5), the air temperature showed an increasing trend with the lower values at the beginning of the survey period (1978–1981, 1983–1986, and 1990–1993), while the highest temperatures were registered at the middle and end of the period (1989–1990, 2002–2003, and 2014–2018). Furthermore, the discharge displayed a succession of dry periods (1988–1990, 1997–1998, 2005–2006, and 2016–2017) with a mean discharge of 0.25 m³·s⁻¹, deficit of 43% compared to the inter-annual average) and humid ones (1980–1981, 1993–1994, 2012–2014, and 2017–2018) with a mean discharge of 0.67 m³·s⁻¹, increase of 52% compared to the inter-annual average).

Additionally, the three strongest flood events over the whole period were recorded (in decreasing order of mean daily discharge) in October 1993 (10.10 m³·s⁻¹), February 2015 (8.94 m³·s⁻¹) and May 2002 (8.54 m³·s⁻¹). Regarding the temperature, the highest daily air temperature was observed in August 2003 (29.5 °C). Furthermore, even if there is no clear relationship between pCO₂, discharge and temperature, pCO₂ showed the highest values during humid periods with low temperature (1990–1994), the highest one being in October 1991 (1 × 10^−1.64 atm).

3.3. Seasonal Patterns Based on Long Term Data Series in Stream Water Chemistry

The seasonal variations of discharge, temperature, pCO₂ and the dynamics of ion concentrations throughout the year were assessed in Figure 6 using the mean monthly values calculated over the whole period 1994–2006 plus 2016–2018.

The discharge and air temperature exhibited seasonal patterns which are characteristic of the Pyrenees Mountains [88,113]. The fluctuations of mean monthly discharge deviation from the interannual monthly average showed two high flow level periods: the first and highest one from December to February in winter (0.66 m³·s⁻¹, 5.0 °C) and a second and lower one between March and May in spring (0.54 m³·s⁻¹, 11.3 °C). On the opposite side, the lowest water flow period was from June to October during summer (0.17 m³·s⁻¹, 19.3 °C). The air temperature fluctuations presented an exact inverse pattern to discharge (Figure 6a).

Moreover, pCO₂ increased in summer, with the highest average peak in August (1 × 10^−2.42 atm) and it began to decrease slowly in autumn reaching its minimum average peak in March (1 × 10^−2.64 atm). Inversely, spring waters were in equilibrium or undersaturated with respect to calcite at the end of summer and oversaturated with respect to calcite the rest of the year with the highest average peak in March (Figure 6b), corresponding to the lowest mean pCO₂.

Regarding major elements, K⁺ and Cl⁻ concentrations were higher during the rainy period in winter and decreased as summer approaches (Figure 6d,f, respectively). Conversely, NO₃⁻ and SO₄²⁻ showed a different pattern with low concentrations during both high flow periods and an accentuated increase in summer especially for SO₄²⁻ (Figure 6e,f).

In addition, (Ca²⁺ + Mg²⁺), HCO₃⁻, and Na⁺ varied slightly during the year (Figure 6c–e, respectively). In fact, (Ca²⁺ + Mg²⁺) and HCO₃⁻ exhibited the maximum concentrations at the end of autumn (October and November), meanwhile, the lowest concentrations occurred in spring for (Ca²⁺ + Mg²⁺) and in summer for HCO₃⁻.

4. Discussion

4.1. Natural vs. Anthropogenic Sources of Elements

To elucidate the export dynamics of major ions, we analyzed in detail their origin and each of the factors that could influence the export patterns of the major elements in BC.

The dissolved elements in stream water result of the interaction between the anthropogenic activities, the atmospheric deposition, the biosphere, the soils and the chemical weathering of rocks and sediments [114].

The natural weathering process due to the reaction of the carbonic acid on minerals, provides dissolved elements such as Ca²⁺, Mg²⁺, Na⁺, K⁺, HCO₃⁻, and SO₄²⁻. These elements originate mainly from the dissolution of silicates, carbonates, and sometimes from evaporite dissolution or from the pyrite oxidation [115,116,117,118]. This is evidenced in the BC stream waters by the distribution of these dissolved elements along the Factor 1 in the PCA (Figure 2a). According to the major lithology of BC (Figure 1), the main source of Ca²⁺ and HCO₃⁻ is the calcite dissolution (Equation (2)), these two dominant ions (Figure S1) being strongly positively correlated (Figure 2a, Table S1). Dolomite and magnesian calcite dissolutions are part of the origin of Mg²⁺ [83], however a positive correlation linked Mg²⁺ and SO₄²⁻ (Table S1 and Figure 2). This relationship is related to their higher relative contribution for both elements during the summer period (Figure 6c,e), but only the first one would come from carbonates (Figure 1).

Indeed, sulfate in stream water may have several origins. According to the different lithologies in the basin, two minerals can be sources of sulfate. The first one is gypsum (CaSO₄), which dissolution produces Ca²⁺ and Mg²⁺, respectively, but does not produce HCO₃⁻. A gypsum origin must be taken into account since although the gypsum outcrop is only about 0.2% [83,85], its dissolution remains still possible due to its very high solubility, even after the calcite and dolomite reach saturation [27]. The second mineral is pyrite (FeS₂), whose oxidation produces sulfuric acid (see Equation (15)). The pyrite presence in limestones and schists is confirmed by Bakalowicz [83] and Debroas [85]. Many authors have mentioned that sulfuric acid can substitute carbonic acid in the dissolution of carbonates [70,76,107,119]. This substitution can be proposed as a hypothesis because of the relative decrease in HCO₃⁻ during summer whereas the SO₄²⁻ concentration increases during low water period (Figure 6e) as illustrated also by the PCA factor 2 (Figure 2). However, one must also keep in mind that another hypothesis can be put forward: the higher relative contribution of enriched-sulfur sources at the outlet during the summer period compared to high water conditions. This might also explain this relative balance between HCO₃⁻ and SO₄²⁻ (see Section 4.2).

Indeed, the influence of sulfuric acid on the dissolution of carbonate could be demonstrated by the molar ratio Ca²⁺ + Mg²⁺ over HCO₃⁻ slightly greater than 0.5 (Figure 3a). It indicates that: (1) some Ca²⁺ + Mg²⁺ comes from the gypsum dissolution, but also (2) some carbonate dissolution could be induced by sulfuric acid. Therefore, because both sulfuric and carbonic acids are weathering actors in BC, the new weathering equation would be as follows [39]:

$$\begin{gathered} \left(\mathsf{\alpha }+2\mathsf{\beta }\right){\mathrm{Ca}}_{1-\mathrm{x}}{\mathrm{Mg}}_{\mathrm{x}}{\mathrm{CO}}_3+{\mathsf{\alpha }\mathrm{H}}_2{\mathrm{CO}}_3{}_{}+{\mathsf{\beta }\mathrm{H}}_2{\mathrm{SO}}_4 \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\to \left(1-\mathrm{x}\right)\left(\mathsf{\alpha }+2\mathsf{\beta }\right){\mathrm{Ca}}^{2+}+\left(\mathrm{x}\right)\left(\mathsf{\alpha }+2\mathsf{\beta }\right){\mathrm{Mg}}^{2+}+{\mathsf{\beta }\mathrm{SO}}_4^{2-}+\left(2\mathsf{\alpha }+2\mathsf{\beta }\right){\mathrm{HCO}}_3^{-} \end{gathered}$$

(18)

This equation is under alkaline conditions and it is pH-dependent process [120]. Also, when bicarbonate ions released will be exported downstream, a significant amount could be transformed to dissolved CO₂ if the pH of stream waters decreases downstream [121].

The influence of sulfuric acid on carbonate dissolution and the gypsum dissolution are both supported by the significant positive correlation between the alkalinity decrease and the relative sulfate content in the Baget stream water (Figure 3b). Nevertheless, during the low water period (Figure 7a,b and Figure 8a), the influence of sulfuric acid on carbonate dissolution is consistent with: (i) the increase of sulfate and the slightly lower pH values in the stream water and (ii) the opposite ratio Mg²⁺/Ca²⁺ pattern observed between HCO₃⁻ and SO₄²⁻ indicating an enhanced dolomite dissolution, which is less soluble than calcite. Moreover, the increase of carbonate dissolution by sulfuric acid during summer causes an undersaturation of stream waters regarding calcite (Figure 8b).

Finally, recent results [122] on spatial isotopic analysis of ${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}$ carried out on spring and stream waters during summer 2019 based on the lithology of their sub-catchments, give arguments for the major origin of sulfates in stream waters. The presence of sulfate released by pyrite oxidation was detected mainly in waters draining the black flysch (${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}=$−8.1‰) and the sulfate from gypsum dissolution in waters draining the evaporite clays (${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}=$13.3‰). However, the ${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}$ isotopic signature of the riverine sulfates (${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}=$11.3‰) at the outlet of BC was close to the gypsum signature, showing that sulfate at BC outlet would be mainly supplied by gypsum dissolution [122], even if the contribution of sulfuric acid to carbonate dissolution is higher during the summer period.

Atmospheric deposition is also a potential source of sulfate. The ${\mathsf{\delta }}^{34}{\mathrm{S}}_{{\mathrm{SO}}_4}$ in precipitation is unknown in BC but Puig et al. [123] reported values between +3.3‰ and +0.2‰ for Vallcebre (Pyrenees) with a potential atmospheric deposition regional background of 7.2‰, and Probst A. and Fritz B. (unpublished data, 1994) reported values between +5.2‰ and +1.8‰ for the Strengbach catchment (Vosges). These data do not allow to discriminate the atmospheric contribution of sulfates for the BC stream waters. However, despite the BC remote position, Binet et al. [82] proposed that a non-negligible part of sulfate would be deposited as dry deposition, particularly by the early sixties. Indeed, some French unpolluted catchments also exhibit high sulfate concentrations [36]. These catchments, which are covered by coniferous forests, were probably also affected by enhanced sulfuric acid inputs due to atmospheric pollution, particularly captured as dry deposits, as already shown for example by Probst et al. [78] and Pierret et al. [124] for the coniferous Strengbach catchment in the Vosges Mountains (France).

After a significant increase until the 1980s, like in other regions in the northern hemisphere [125], acid atmospheric deposition was reduced by half from 1980 to 2005, and has probably decreased at the Baget site like in other sites in France located in remote areas [126]. This might explain the general decreasing trend observed in stream water (SO₄²⁻ (meq·L⁻¹) = −6 × 10⁻⁶ (date) + 0.6; n = 763; p < 0.001 from 1982 to 2018, Figure 4 and Table S3). This trend is very consistent (similar slope) with that observed in the Strengbach stream, which is strongly linked to the decreasing atmospheric sulfates input by precipitation [127].

Even if the atmospheric contribution to riverine sulfate in BC might be lower than the lithological contributions, the long-term influence on carbonate dissolution and CO₂ consumption, has been proved to be non-neglectable at least during the last decades (16–25% of the total Ca²⁺ + Mg²⁺ concentration could have been due to an increase of carbonate dissolution by acid deposition [82]). In addition, sulfuric acid atmospheric deposition decreases during these last decades would have led to stream water Ca²⁺ + Mg²⁺ decrease contrary to what happens, indicating that other processes are involved such as indirect effect of temperature increase or land use change (see Section 4.2).

In the BC, chloride mainly originates from rainfall because there is no anthropogenic activities (such as agriculture) which could supply chloride and there is no evaporitic formation containing halite or sylvite (Figure 1b) [36,128]. This is supported in the PCA (Figure 2b), where this conservative anion is opposed to pCO₂ and Cl⁻ concentration increasing when pCO₂ decreases during the rainy periods. The positive correlation between K⁺ and Cl⁻ (Figure 2a, Table S1), and their joined substantial increase during the rainy season in winter, argue for an atmospheric origin of the two elements, and possibly a soil contribution for K⁺ since it is a largely recycled element not used by vegetation during winter time (Figure 6d,f). Furthermore, Cl⁻ and NO₃⁻ increase in autumn with the first flood events at the beginning of the hydrological year. These high concentrations of Cl⁻ and NO₃⁻, with peaks in October and November, would be due to the onset of rainfall and the leaching of the fir–beech forest or epikarst (see Section 4.3, for NO₃⁻), respectively (Figure 6f) [36].

4.2. Role of Discharge, Temperature and Vegetation on Stream Water Composition and Trends

The hydrological cycle plays a major role in controlling the carbon dynamics in the karst critical zone [129]. The stream discharge, which shows the highest variability during the year, is a key factor in relation to mineral dissolution processes (Equation (4)). Its primary control exerting over fluxes of dissolved elements has long been demonstrated [7,34,37,56,130,131].

However, due to dilution processes, discharge also contributes to concentration decrease. This is mainly due to the input of rainwater and surface runoff, which are less concentrated than groundwater flow [132,133,134]. This could be illustrated by the decrease in both pCO₂ and ions loads due to high discharge during winter and spring in BC (Figure 8c–e).

In addition, it was observed that the decrease in discharge would lead to an increase in Ca²⁺ and HCO₃⁻ in BC (Figure 6a,c,e). A similar effect was observed on the Houzhai catchment, southwestern China, where a decrease of rainfall between 2007 and 2013 led to an annual mean increase of dissolved element concentrations in stream water [53,129].

Moreover, in BC, during the low water flow period in summer, some tributaries dry out completely, causing particular areas to increase their relative contribution of dissolved elements to the waterbody. An evidence of this process can be supported by the contrast between discharge and Mg²⁺ and SO₄²⁻ on factor 2 in PCA (Figure 2a) and the significant relative increase of SO₄²⁻ in summer (Figure 6e). As evoked before, the change in the relative contribution of water sources can explain the variation of around 10% in Mg²⁺ (coming from dolomitic lithology) and of around 20% in SO₄²⁻ (coming from pyritic lithology) in the piper diagram (Figure S1).

In addition to discharge, the temperature is another factor of similar importance, but with an opposite tendency (Figure 2 and Figure 6a, Table S1). This antagonism related to the climate characteristics of the Pyrenees Mountains (high water flow with lower temperature in winter and early spring, and low water flow with higher temperature during summer and early autumn). It is also supported by the negative relationship between temperature and both Cl⁻ and K⁺ (Figure 2), which are rainfall or vegetation dependent elements.

Indeed, temperature is an important catalytic agent in the dissolution reactions of minerals. However, carbonate dissolution was shown as an exception. The increase in temperature contributes to a decrease in the dissolution because the solubility of calcite increases when temperature decreases [12,98]. In BC, this is opposite to the observed long-term trend patterns of increasing temperature and Ca²⁺ + Mg²⁺ and HCO₃⁻ (Figure 4), indicating the influence of other key factors. Indeed, air temperature is considered as the main factor related to seasonal changes in soil pCO₂ as demonstrated by the nice anti-clockwise hysteresis observed in Figure 8f. In summer, low soil moisture could enhance the diffusion of biogenic CO₂ to the atmosphere and could also limit the growth rate of plants by decreasing their role in the soil [135]. Despite these processes, Zhao et al. [39] showed that the increase in temperature not only improves pCO₂ derived from root respiration, but also stimulates the decomposition of organic matter in the soil. Therefore, the increase in temperature would rise the generation of biogenic CO₂, which can enhance the dissolution of carbonates. The factor 1 in PCA (Figure 2a) shows this influence, in which the temperature is resembled with SO₄²⁻ but also with Ca²⁺ and HCO₃⁻. Additional evidence is shown with the positive relationship between Ca²⁺ + Mg²⁺ and HCO₃⁻ with pCO₂ (Figure 8c,d), denoting that a higher pCO₂ level in water promotes higher carbonate dissolution rates. Figure 8f also shows a net pCO₂ decrease from end-of-summer to winter, probably due to the increase in water flow level and to the decrease of temperature [136,137]. But outside the dilution effect during the high flow period, the lower content of bicarbonate could be attributed to the consumption of DIC by secondary formation of CaCO₃ and also to losses by CO₂ degassing into the atmosphere (inverse reaction of Equation (4)). This hypothesis is supported by a higher SIc and a lower pCO₂ during high flow conditions (Figure 8b) [138]. The lowest pCO₂ values could be also related to a lower residence time during the spring period. After this period, stream water pCO₂ increases again with temperature increase until August. Previous studies showed similar behavior at stream water BC: pCO₂ was 1 × 10^−1.9 atm in summer and 1 × 10^−2.6 atm in winter [83].

At the annual and inter-annual scale, it is difficult at this stage (more data in Ulloa-Cedamanos PhD in progress) to provide complete evidence that the trends observed for Ca²⁺ and HCO₃⁻ are mostly related to climate change or land-use change (Figure 5). The observed significant increasing trends over time of Ca²⁺ + Mg²⁺ and HCO₃⁻ (Figure 4) could be explained by the temperature increase (p < 0.01) as described above and the discharge decrease (p < 0.01 from 1978 to 2006). However, surprisingly, under usual hydroclimatic conditions without extreme events, there is no relationship between pCO₂ level and Ca²⁺ + Mg²⁺ and bicarbonates (Figure 5). The pCO₂ peaks observed in the early 1990s could be due to a side effect of the long-term pumping test of 1990 (Figure 4 and Figure 5). At this period, the decrease of the discharge corresponds to a pH decrease. This could be interpreted as: (1) a decrease of the relative contribution of waters draining carbonates and (2) a relative increase of the contribution of waters draining lithological formations containing pyrites. This pH decreases of the water would have caused the transformation of bicarbonate ions into dissolved CO₂ (Equation (19)). In fact, Mangin [79] has shown that a very severe low water flow conditions can impact the hydrological functioning of a karst system for a long period (at least 5 years in the case of Baget).

At the seasonal scale, most of stream waters in BC had SIc above 0 (Figure 8b), indicating that Equation (4) would occur in the opposite direction in response to calcite oversaturation, allowing calcite precipitation and CO₂ degassing to the atmosphere. Indeed, it has been shown by several studies [38,53,139] based on a theoretical diffusion model of CO₂ fluxes between river waters and atmosphere, that CO₂ evasion is important in the riverine carbon cycle. Conversely, August is the only month showing an undersaturation with regard to calcite and holding the highest pCO₂. In fact, the pCO₂ increase leads to the increase in calcite dissolution needed to reach the equilibrium.

Additionally, pCO₂ of stream water in winter and spring tends towards low values but remaining always greater than the atmospheric pCO₂ (>1 × 10^−3.5 atm). This shows that the direct atmospheric CO₂ contribution to DIC can be ignored in BC. At the same time, this could be a clue showing that groundwater stored near the surface in the epikarst could be in relation to a gas phase permanently enriched with the soil pCO₂ [28]. Therefore, the role of epikarst on the dynamics of carbonate dissolution is important to consider.

4.3. Respective Control of Karst and Epikarst on Stream Water Chemical Composition and Trends

The epikarst drives rainwater into a network of drainage pathways [140,141]. During the rainy season, a part of rainwater which recharges the karst aquifer is partly stored in the epikarst causing three different processes: (i) evaporation processes of a large fraction (with concentration of dissolved salts), (ii) slow infiltration through the fine cracks and pores, and (iii) delayed infiltration during the heavy rains recharging the epikarst (“piston effect”) [28,142,143]. These processes can be observed directly and indirectly in BC. For example, all dissolved elements, even SO₄²⁻, would be concentrated in the epikarst during summer, a period of intense mineralization in the soils. Later, all these major elements subjected to a delayed infiltration by storage near the surface, are released with the first heavy rainfalls in autumn, in particular NO₃⁻ (Figure 6f and Figure 8c,d). This “piston effect” may be displayed by the positive correlation between dissolved elements coming from the dissolution of carbonates, mainly calcite (Ca²⁺ and HCO₃⁻) and those coming from rainwater or vegetation interface recycling (Cl⁻ and K⁺) (Figure 2, Table S1).

The decrease of sulfate concentrations observed in autumn with the inputs of the first rainfalls, illustrate the contribution of water draining lithological areas less concentrated in sulfate (see November position on Figure 8a). Subsequently, since the pre-stored water in the epikarst has been released in winter, the ions exhibiting the greatest increase are those mainly originating from rainfall (i.e., Cl⁻), or less recycled by the vegetation due to the dormant period (i.e., K⁺) (Figure 6f,d, respectively). Lastly, during the spring, part of rainwater incoming will be used to recharge the epikarst [79] and then, during summer, the epikarst water is progressively mineralized.

The epikarst also appears as a reservoir with higher pCO₂ values than those measured in the surrounding soils and caves, due to an active exchange between the sub-surface flow and deep flow [53], and to a CO₂ accumulation effect [47,144]. Therefore, during the temperature conditions in summer, carbonate dissolution is enhanced by high concentrations of soil CO₂. Then, Ca²⁺ and HCO₃⁻ concentrations increase in epikarst aquifers. Afterwards, large amounts of Ca²⁺ + Mg²⁺ and HCO₃⁻ can be released during the epikarst drainage in autumn (Figure 8b,c, respectively). Also, the higher exchange between epikarst zone and the atmosphere can be evaluated by the pCO₂ decrease and the calcite oversaturation (SIc > 0) during the epikarst drainage (Figure 8b).

Lastly, the rise in ion concentrations observed in stream water during autumn with the first stormflow events occurring after the strong summer evaporation period, is supposed to come from soil horizons and epikarst drainage. This “flushing effect” was evidenced by Walling and Foster [145] and later on by several authors on catchment draining various lithologies and under various climate conditions: on silicate basins such as the Mont-Lozère, the Strengbach and Ringelbach [72,111], the silicate-rich Chari–Logone River catchment in Chad [146], the Vantaa River boreal catchment in southern Finland [147], the Heidingzi agricultural catchment in China [137], and also on calcareous basins such as the Oka River forested catchment in Spain [148], the Middle Bussento Karst System in Italy [149], and the Chenqi agricultural karstic catchment [150].

5. Conclusions

This study investigated the response of a small, remote, and forested karst catchment to the impact of inter-annual and seasonal changes on stream water hydrochemistry using a 40-year period survey, through instantaneous, monthly, and annual data series.

Over the long-term, the significant rising trends in Ca²⁺ + Mg²⁺ and bicarbonate concentrations could be related to the global warming as shown by the air temperature increase, which enhanced the biogenic activity and then the carbonate dissolution. Whereas discharge and pCO₂ trends are not significant, except for discharge during the period 1978–2006 (decreasing trend). In contrast, a decreasing trend was observed for sulfates, which might be related to atmospheric deposition decrease over the last decades. Indeed, in the context of global warming and hydroclimatic fluctuations, the increasing trends of Ca²⁺ + Mg²⁺ and bicarbonates could be a consequence of the buffering effect of carbonate dissolution on the global increasing trends of atmospheric CO₂, and then on short time scales, an important CO₂ sink.

The seasonal patterns allowed to decipher the origin and the export dynamics of major ions in the BC: (i) a major lithological origin for the main base cations (Ca²⁺ and Mg²⁺) and anions (bicarbonates and sulfates) from carbonates and/or sulfate mineral dissolution; (ii) a major atmospheric origin for chloride and potassium, and (iii) a biological origin for nitrates. Besides, calcite undersaturation and carbonate dissolution increase were observed during the warm climatic conditions and low flow period, in relation with the biogenic CO₂ production.

Outside the classical controls of water–rock interactions (depending on water residence time), dilution effect (by precipitation input), and concentrating effect (by evapotranspiration output), the water discharge regime induces a variable importance of water source contribution from the different lithologies. This is illustrated by the increased sulfate to bicarbonate ratio up to two-times higher during low water flow period. Moreover, the role of the epikarst contribution was highlighted, and particularly illustrated by the first heavy rains in autumn, which cause by “piston effect” and “flushing effect”, the release of large amounts of Ca²⁺ + Mg²⁺ and bicarbonates previously stored in the epikarst solutions.

Such long-term studies of the critical zone must be promoted since they have proven to be efficient in identifying: water and element sources and trends, weathering and biogeochemical processes, at various time scales, and the key controlling parameters. In a period of contrasting regional climate patterns, such critical zone survey provides robust validation data for the development of predictive models.