Glacier loss and vegetation expansion alter organic and inorganic carbon dynamics in high-mountain streams

High-mountain ecosystems are experiencing the acute effects of climate change, most visibly through glacier recession and the greening of the terrestrial environment. The streams draining these landscapes are affected by these shifts, integrating hydrologic, geologic, and biological signals across the catchment. We examined the organic and inorganic carbon dynamics of streams in four Alpine catchments in Switzerland to assess how glacier loss and vegetation expansion are affecting the carbon cycle of these high-mountain ecosystems. We find that the organic carbon concentration and fluorescence properties associated with humic-like compounds increase with vegetation cover within a catchment, demonstrating the increasing importance of allochthonous dissolved organic carbon sources following glacier retreat. Meanwhile, streams transitioned from carbon dioxide sinks to sources with decreasing glacier coverage and increased vegetation coverage, with chemical weathering and soil respiration likely determining the balance. Periods of sink behavior were also observed in non-glaciated streams, possibly indicating that the chemical consumption of carbon dioxide could be more common in high-mountain, minimally vegetated catchments than previously known. Together, these results demonstrate the dramatic shifts in carbon dynamics of high-mountain streams following glacier recession, with significant changes to both the organic and inorganic carbon cycles. The clear link between the terrestrial and aquatic zones further emphasizes the coupled dynamics with which all hydrologic and biogeochemical changes in these ecosystems should be considered, including the carbon sink or source potential of montane ecosystems.

1 Introduction

The effects of climate change on high-mountain areas are dramatic, with temperatures increasing approximately twice as quickly as in lower-elevation areas (IPCC, 2021). With glacial retreat, the streams draining these landscapes are experiencing significant change in the timing, magnitude, and source of flows (Kneib et al., 2020; Mackay et al., 2019). The terrestrial environment is also shifting with the expansion of vegetation spatially (i.e., to higher elevations) and temporally (i.e., longer growing season), thereby impacting the hydrology, biogeochemistry, and ecology of catchment streams (Knight and Harrison, 2014; Brighenti et al., 2019). In the Swiss Alps, recent work has highlighted the rapid “greening” of high-mountain areas and decreasing snow and ice cover (Rumpf et al., 2022). While the implications of climate change for terrestrial ecosystems have been examined broadly (Finstad et al., 2016), the impact that these changes will exert on the streams draining these landscapes is much less explored (Beniston et al., 2018). Given the global extent and integral role of streams in connecting high-mountain areas with downstream ecosystems (Milner et al., 2017), exploring how these landscape alterations will affect the carbon dynamics of streams is critical to contextualize their role in the global cycle (Horgby et al., 2019c).

High-mountain streams are tightly linked to the catchment that they drain (Milner et al., 2009; Brighenti et al., 2019). In particular, the presence of glaciers dominates stream hydrology (Kneib et al., 2020), with significant biogeochemical and ecological implications. For example, as glaciers generally provide the majority of water to their proglacial streams, solute concentration and flux are frequently controlled by the contents and magnitude of glacier meltwater (Bergstrom et al., 2021). In particular, dissolved organic carbon (DOC) ice-locked within glaciers can be the dominant source of DOC to the proglacial stream upon melting (Colombo et al., 2019). The lability of this glacier-derived DOC is often high, serving as a major source of carbon fueling downstream metabolism (Hood et al., 2009). Glaciers are also associated with high rates of geochemical weathering, both underneath the glacier (Anderson et al., 1997) and in the proglacial stream (St. Pierre et al., 2019). The weathering of both carbonate and silicate minerals can consume atmospheric carbon dioxide (CO₂), whereby CO₂ dissolved in water is converted to bicarbonate through these reactions (Donnini et al., 2016). These reactions involve significant transformations of dissolved inorganic carbon (DIC) and potentially consume large amounts of CO₂ in the process (Hodson et al., 2000).

As glaciers shrink, there is generally a concomitant increase in soil development and vegetation cover within catchments (Guelland et al., 2013). Higher vegetation cover and soil development provides a pool of organic carbon for export to the streams (Garcia et al., 2015). From this change, increases in stream DOC concentration are likely. Indeed, increased DOC in aquatic ecosystems globally has been directly linked to the greening of the terrestrial landscape (Finstad et al., 2016). Elevated aquatic DOC has implications for ecosystem respiration, productivity, and water quality (Roulet and Moore, 2006; Hongve et al., 2004). This change in DOC source also implies changes in the quality of stream organic matter (Zhou et al., 2019), which could further alter stream metabolic regimes by promoting heterotrophy (Bernhardt et al., 2018; Duarte and Prairie, 2005; Boix Canadell et al., 2020). In terms of inorganic carbon, soils frequently represent the dominant source of CO₂ to streams, as the products of soil respiration are transported to the stream via groundwater (Hotchkiss et al., 2015). Thus, as soils develop and allow for the expansion of vegetation in mountain catchments, emissions of CO₂ from the aquatic system may be promoted as the products of soil respiration are transported to the stream and emitted.

Given these complex relationships, consideration of both glacial influence and the terrestrial environment at large is key to fully contextualize how climate change may alter carbon flows to and from mountain streams. Moreover, both the organic and inorganic carbon components must be evaluated to complete this cycle, providing perspective on the relative influence of different catchment properties. In this study, we aim to evaluate landscape effects on dissolved organic and inorganic carbon dynamics in high-mountain streams across a glacial, vegetation, and elevation gradient. By comparing dissolved carbon concentration and fluxes across these gradients, we can directly assess the relative impact of glacial retreat and catchment greening. We hypothesized that the presence of glaciers would drive CO₂ consumption (St. Pierre et al., 2019) and that the loss of glacier influence would elevate the role of catchment soils as a source of CO₂ (Crawford et al., 2015). We also expected that these landscape transformations would increase the role of allochthonous dissolved organic carbon in the stream (Fasching et al., 2016), with consequential changes to the quality of organic matter.

2 Methods

Samples of DOC and DIC as well as in situ sensor measurements of dissolved carbon dioxide (pCO₂) were collected in 12 streams in the high-mountain area of the western Swiss Alps over 5 years, from 2016 to 2020. The sampling locations covered a broad range of catchment glacier coverage, vegetation coverage, and elevation, providing for a space-for-time substitution approach in which streams draining lower-elevation, lower-glacier-cover, and higher-vegetation-cover catchments represent potential future conditions of higher-elevation, higher-glacier-cover, and lower-vegetation-cover catchments. With this spatial design, we can to evaluate how these carbon constituents may evolve with ecosystem processes following anthropogenic climate change. Consideration of various other water quality and catchment properties (e.g., dissolved oxygen, inorganic carbon isotopes, and dissolved organic matter fluorescence) provides further insight into changes in the relative contribution of geochemical weathering, in-stream processes, and terrestrial inputs within these streams.

Figure 1 Map of the 12 study sites within four catchments of the Alps in southwestern Switzerland (glacial cover and stream network from swissTLM3D; Geodata © Swisstopo).

Table 1 Catchment characteristics.

2.1 Site description

Our 12 stream sampling locations were equally distributed within Vallon de Nant, Champéry, Valsorey, and Val Ferret, four catchments in the western Swiss Alps (Fig. 1). These sites are part of the METALP (Metabolism in Alpine Streams) project (https://metalp.epfl.ch, last access: 12 December 2022), which has been described extensively in previous studies (Ulseth et al., 2019; Boix Canadell et al., 2019; Horgby et al., 2019c); as part of this project, numerous hydrological and biogeochemical parameters have been monitored in the aforementioned catchments since 2016. The drainage areas vary from 0.31 to 23.2 km², and mean catchment elevation ranges from 1778 to 2892 m (Table 1). Vegetation cover is highest at lower elevations and ranges from approximately 94 % to 21 % coverage of the catchment.

The geology of the Champéry and Vallon de Nant catchments is dominated by the presence of limestones, calcareous shales, and flysch (Burri et al., 1999). Val Ferret is characterized by limestones and sandstones showing pronounced schistosity, whereas Valsorey is underlain primarily by a metamorphic lithology that is comprised of gneisses, crystalline shales, and blue-gray schists. Additionally, in Valsorey, glacial cover accounts for about one-third of the total surface, with blue-gray schists most prominent at high elevations beneath and near the glaciers. Soils in this region are generally young, poorly developed Leptosols and Fluvisols (Egli et al., 2010), with an organic content that increases at lower elevations (Hoffmann et al., 2014). While soil organic content has not been measured in these catchments, measurements in nearby areas of the Swiss Alps have observed a soil organic content ranging from < 0.1 kg C m⁻² in soils at around 2150 m elevation to around 2.0 kg C m⁻² at 1000 m elevation (Hoffmann et al., 2014; Egli et al., 2010).

2.2 Grab sampling and sensor measurements

Grab sampling of various physical and chemical parameters was carried out at all sampling sites at approximately monthly intervals during the snow-free season. These parameters include DOC, dissolved organic matter (DOM) fluorescence, major ions, pCO₂, DIC, and alkalinity. The analysis of these analytes has been described previously (Horgby et al., 2019b; Boix Canadell et al., 2019). Concentrations of major anions (Ca²⁺, Mg²⁺, K⁺, and Na⁺) and cations (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and Cl⁻) were measured on stream water filtered through 0.22 µm filters (mixed cellulose ester; Millipore, Billerica, MA, USA) using ion chromatography (Metrohm 930 Compact IC Flex; Aargau, Switzerland).

Samples for the quantification of DOC and DOM fluorescence are filtered through pre-combusted 0.45 µm GF/F filters (Whatman, Cytiva, Massachusetts, USA) into acid-washed and pre-combusted 40 mL amber glass vials. Samples were kept refrigerated and analyzed for concentration within 24 h of collection. DOC concentration was measured using a Sievers M5310 C TOC analyzer (GE Analytical Instruments, New York, USA) with an accuracy of ± 2 % µg C L⁻¹, a precision of < 1 % relative standard deviation, and a detection limit of 22 µg C L⁻¹. Samples were taken in triplicate, and the mean concentration was used; outliers were removed if the concentration exceeded 3 standard deviations from the mean. Calibration standards ranged from 0.05 to 1 mg C L⁻¹. DOM fluorescence excitation–emission matrices (EEMs) were created by measuring the fluorescence intensity of samples within a 1 cm cuvette across a range of excitation (240–450 nm, 2 nm increment) and emission (211.19–620.23 nm, 2 nm increments) wavelengths using an Aqualog^® optical spectrometer (Horiba, Kyoto, Japan). Validation scans were performed prior to sample analysis to validate instrument performance, such as measuring the water Raman signal-to-noise ratio and emission calibration using a sealed, standard cuvette of Milli-Q water (Type-1 laboratory reagent water) to ensure that the Raman peak position was at 297 nm ± 1 nm. Milli-Q water was used as a blank, which was employed to remove background fluorescence from the spectra. Absorbance was also measured within a 10 cm cuvette with a PerkinElmer LAMBDA 850 spectrophotometer (Massachusetts, USA).

Duplicate samples for dissolved inorganic carbon (DIC) concentration and the relative stable carbon isotopic composition (δ¹³C-DIC) were filtered through 0.2 µm membrane filters into acid-washed 12 mL glass Exetainer vials and stored refrigerated until analysis. A total of 2 mL of sample was injected into synthetic air-filled, septum capped 12 mL Exetainer vials containing 200 µL of 85 % phosphoric acid to convert all DIC to gaseous CO₂. The resulting gas phase was then analyzed on a cavity ring-down spectrometer equipped with a Small Sample Isotope Module 2 (CRDS-SSIM2, Picarro Inc., California, USA) and converted into DIC concentration.

Additionally, each monitoring station was instrumented with sensors measuring physical and chemical parameters of the water or air at a 10 min frequency, including water temperature, dissolved oxygen, dissolved carbon dioxide (pCO₂), and depth. The specifications and the calibration and maintenance procedures of these sensors have been described previously (Horgby et al., 2019c; Boix Canadell et al., 2020). Stream pCO₂ was measured using a CARBOCAP^® GMP252 probe (Vaisala, Vantaa, Finland) within a porous polytetrafluoroethylene (PTFE) semipermeable membrane. The probes were then protected with a fine-grained mesh and a metal casing. Raw data were adjusted, according to the manufacturer's recommendations, for barometric pressure and water temperature. All pCO₂ sensors were tested in the laboratory using certified gas mixtures of CO₂ diluted in synthetic air to final concentrations of 0, 400, and 2000 ppmv prior to deployment.

Discharge was calculated using rating curves relating water depth to discharge (Boix Canadell et al., 2021), where direct measurements of discharge were made using slug injections of sodium chloride (NaCl) as a conservative tracer (Gordon et al., 2004). Additionally, when stream conditions allowed, 10 random measurements of stream depth were collected to provide a measure of average stream morphology for comparison with measurements recorded by the sensor installed on the side of the stream.

2.3 CO₂ saturation and efflux

From the 10 min sensor data, the daily median concentration of CO₂ was found for all sample locations during the monitoring period (Horgby et al., 2019c). The saturation and efflux for these values were then estimated using measurements of stream water temperature as well as estimates of barometric pressure, atmospheric concentration of CO₂, and gas exchange velocity. Barometric pressure was obtained from the MeteoSwiss weather station network (Swiss Federal Office and Meteorology and Climatology). The Col du Grand St Bernard station (elevation 2473 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) was used for the Valsorey and Val Ferret catchments, and the Evionnaz station (482 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and Les Diablerets (2964 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) stations were used for Champéry and Vallon de Nant stations. Barometric pressure at each monitoring station (P_site, mbar) was adjusted for site-specific elevation and temperature as follows:

$$\begin{array}{}\text{(1)}& P_{\text{site}}=P_{\mathrm{0}}\exp \left({\displaystyle \frac{-gM(h-h_{\mathrm{0}})}{RT}}\right),\end{array}$$

where P₀ (mbar) is the barometric pressure measured at the MeteoSwiss station; h₀ and h (m) are the altitude of the meteorological and monitoring stations, respectively; g is gravity acceleration (9.81 m s⁻²); M is the molar mass of air (0.0289644 kg mol^−L); and R is the universal gas constant (8.31432 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). The temperature of air T_air (^∘C) at the METALP stations is estimated through the temperature T₀ (^∘C) measured at the MeteoSwiss station, where the regional temperature gradient $\mathrm{\Delta }T/\mathrm{\Delta }h$ is set to −0.54 ^∘C$/$100 m, obtained from air temperature data collected during the 1990–2020 period by the Evolène-Villa (1427 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and the Montana (1825 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) weather stations (MeteoSwiss; Deluigi et al., 2017).

$$\begin{array}{}\text{(2)}& T_{\text{air}}=T_{\mathrm{0}}-\left((h-h_{\mathrm{0}})\cdot {\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }h}}\right)\end{array}$$

Sensor measurements of pCO_2,raw (ppm) were then adjusted to site-specific temperature and barometric pressure following the ideal gas law:

$$\begin{array}{}\text{(3)}& p{\mathrm{CO}}_{\mathrm{2},\mathrm{corr}}=p{\mathrm{CO}}_{\mathrm{2},\mathrm{raw}}\cdot {\displaystyle \frac{P_{\text{site}}}{\mathrm{1013}}}\cdot {\displaystyle \frac{\mathrm{298}}{T_{\text{water}}}},\end{array}$$

where P _site (mbar) is the barometric pressure at each location and T_water (K) is the measured water temperature. The dissolved CO₂ concentration (CO_2,water, µmol L⁻¹) was then derived by multiplying the corrected pCO_2,corr by Henry's constant K_H ($\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{atm}}^{-\mathrm{1}}$) at each site:

$$\begin{array}{}\text{(4)}& {\mathrm{CO}}_{\mathrm{2},\mathrm{water}}=p{\mathrm{CO}}_{\mathrm{2},\mathrm{corr}}\cdot K_{\mathrm{H}}.\end{array}$$

K_H is a function of the water temperature in kelvin (T_water), where A is 108.3865, B is 0.01985076, C is −6919.53, D is −40.4515, and E is 669365 according to Plummer and Busenberg (1982):

$$\begin{array}{}\text{(5)}& K_{\mathrm{H}}={\mathrm{10}}^{A+B\cdot T_{\text{water}}+\frac{C}{T_{\text{water}}}+D\cdot {\text{log}}_{\mathrm{10}}\left(T_{\text{water}}\right)+\frac{E}{T_{\text{water}}^{\mathrm{2}}}}.\end{array}$$

A corresponding dissolved equilibrium concentration of CO₂ (CO_2,sat, µmol L⁻¹) was calculated for each sensor measurement at each site using an estimate of daily mean atmospheric CO₂ (CO_2,air),

$$\begin{array}{}\text{(6)}& {\mathrm{CO}}_{\mathrm{2},\mathrm{sat}}={\mathrm{CO}}_{\mathrm{2},\mathrm{air}}\cdot K_{\mathrm{H}},\end{array}$$

by adjusting measurements of CO₂ concentration at Jungfraujoch (freely available at http://www.climate.unibe.ch, last access: 12 December 2022) for differences in barometric pressure and temperature,

$$\begin{array}{}\text{(7)}& {\mathrm{CO}}_{\mathrm{2},\mathrm{air}}={\mathrm{CO}}_{\mathrm{2},\mathrm{Jungfraujoch}}\cdot {\displaystyle \frac{P_{\text{site}}}{P_{\text{Jungfraujoch}}}}\cdot {\displaystyle \frac{T_{\text{Jungfraujoch}}}{T_{\text{site}}}}.\end{array}$$

The standard gas transfer velocity (k₆₀₀, m d⁻¹) was calculated using the relationships developed by Ulseth et al. (2019) and extrapolated from the same 12 streams used in this study:

$$\begin{array}{}\text{(8)}& {\displaystyle }\ln \left(k_{\mathrm{600}}\right)\text{ for }{e}_{\mathrm{D}}>\mathrm{0.02}=\mathrm{1.18}\cdot \ln \left(e_{\mathrm{D}}\right)+\mathrm{6.63},\text{(9)}& {\displaystyle }\ln \left(k_{\mathrm{600}}\right)\text{ for }{e}_{\mathrm{D}}<\mathrm{0.02}=\mathrm{0.35}\cdot \ln \left(e_{\mathrm{D}}\right)+\mathrm{3.10}.\end{array}$$

Here, e_D is the stream energy dissipation rate (m² s⁻³), which is obtained by multiplying the gravity acceleration (9.81 m s⁻²) by slope (S, unitless) and streamflow velocity (V, m s⁻¹):

$$\begin{array}{}\text{(10)}& e_{\mathrm{D}}=g\cdot S\cdot V.\end{array}$$

Velocity was calculated according to the hydraulic geometry scaling proposed by Horgby et al. (2019c) for these streams:

$$\begin{array}{}\text{(11)}& V=\mathrm{0.668}\cdot Q^{\mathrm{0.365}},\end{array}$$

where Q is discharge (m³ s⁻¹). To convert k₆₀₀ to $k_{{\mathrm{CO}}_{\mathrm{2}}}$ (Eq. 11), we used the temperature-dependent Schmidt scaling according to Wanninkhof (2014):

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}{\mathit{Sc}}_{{\mathrm{CO}}_{\mathrm{2}}}=& \mathrm{1923.6}-\mathrm{125.06}\cdot T_{\mathrm{Water}}+\mathrm{4.3773}\cdot T_{\mathrm{Water}}^{\mathrm{2}}\\ & -\mathrm{0.85681}\cdot T_{\mathrm{Water}}^{\mathrm{3}}+\mathrm{0.00070284}\cdot T_{\mathrm{Water}}^{\mathrm{4}},\end{array}\text{(13)}& {\displaystyle }{k}_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{k_{\mathrm{600}}}{{\left(\frac{\mathrm{600}}{{\mathit{Sc}}_{{\mathrm{CO}}_{\mathrm{2}}}}\right)}^{-\mathrm{0.5}}}}.\end{array}$$

The CO₂ efflux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$, $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}-\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) was then calculated as follows:

$$\begin{array}{}\text{(14)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=k_{{\mathrm{CO}}_{\mathrm{2}}}\times ({\mathrm{CO}}_{\mathrm{2},\mathrm{water}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{sat}}).\end{array}$$

2.4 Parallel factor analysis (PARAFAC) modeling

PARAFAC modeling of fluorescence EEMs was used to identify and determine the main fluorescence components of DOM present across collected water samples and was conducted using the staRdom (Pucher et al., 2019) and eemR (Massicotte, 2019) R packages. Preprocessing of EEMs was necessary prior to PARAFAC development (Murphy et al., 2013; Stedmon and Bro, 2008). Briefly, spectra were corrected for instrument-specific effects, blank subtraction, and inner-filter effects. First- and second-order Rayleigh scattering was removed, and corrected EEMs were normalized to Raman units – RU (Murphy et al., 2010). A total of 220 samples were included for model development. The final PARAFAC model was validated using split-half analysis. The resulting components were compared to previously published fluorescence components from aquatic ecosystems in the OpenFluor database (Murphy et al., 2014).

2.5 Statistical analyses

All statistical analyses were performed in MATLAB and Statistics and Machine Learning Toolbox (MATLAB, 2021, version 9.10, R2021a, MathWorks, Massachusetts, USA). Differences in concentration or saturation between groups of streams were investigated using Kruskal–Wallis tests. Simple linear regression was used to evaluate relationships between DOC concentration or CO₂ saturation with catchment properties and water quality parameters. The Pearson correlation coefficient (r) and coefficient of determination (r²) were used to determine the strength of correlations, with the Pearson correlation coefficient used to show the direction of interaction.

The highly correlated nature of potential explanatory variables limited interpretability for CO₂ saturation; thus, we used partial least squares (PLS) regression to identify variables important for predicting median CO₂ saturation at each site. PLS is a method that is well designed for datasets with many collinear predictor variables and when the number of observations is small relative to the number of predictor variables (Wold et al., 1984; Carrascal et al., 2009; Nash and Chaloud, 2011). Here, our response variable is the median CO₂ saturation of each stream location, and 39 predictor variables (standardized within the PLS model) are included (Table S1 in the Supplement).

A Monte Carlo cross-validation method assessed the predictive ability of the resulting PLS model, where the model was fitted with a subsample of data. The calibration-to-validation ratio was set to 0.8, following Onderka et al. (2012), and the resulting fitted models were then tested on the validation set. This process was repeated 500 times. The mean cross-validated goodness of prediction (Q²) was then compared to the original model fit (R²Y). Following this, the strength of each predictor variable within the model was analyzed using variable importance in the projection (VIP), where highly important variables had a VIP > 1.0 (Eriksson et al., 2001). Additionally, moderately important (0.8 < VIP < 1.0) or less influential (VIP < 0.8) variables were identified.

Finally, catchment areal fluxes of CO₂, DIC, and DOC were calculated using catchment area and estimates of stream surface area. We focus on the snow-free period, 1 July through 31 October (Deluigi et al., 2017), to exclude snow cover as a confounding factor affecting gas exchange. Concentration and gas exchange rates are considered constant within subcatchments. An estimation of the network stream area was computed as the product of the stream length and width during this snow-free period. Perennial stream length was extracted from the large-scale topographic landscape model of Switzerland (swissTLM3D) and compared to a 2 m resolution digital elevation model stream network (swissALTI3D). Considering the complexity of the network and its remoteness, stream widths were estimated using aerial images with a 25 cm pixel resolution, with a minimum of one width measurement per stream order. An average of 187 width estimates were made per catchment. The calculation of areal flux for CO₂ is particularly uncertain, as stream surface area (Paillex et al., 2020), gas exchange (Ulseth et al., 2019), and pCO₂ (Horgby et al., 2019b) are each highly dynamic in high-mountain river networks. Thus, these estimates remain approximations intended to provide perspective on the relative balance of dissolved carbon constituents in these stream networks rather than robust calculations of flux. We consider CO₂ as a vertical flux, either into or out of the stream, whereas DOC and DIC are downstream fluxes. The downstream DIC flux inherently includes the downstream transport of CO₂.

Table 2 Median concentration of DOC and DIC, percent saturation of CO₂ and O₂, and isotopic composition of DIC for the 12 sites. Concentration and isotopic composition are summarized from grab samples, whereas CO₂ and O₂ saturation are summarized from sensor data.

3 Results

3.1 Dissolved carbon concentrations

The overall median concentration of DOC was 0.22 mg C L⁻¹, with site-specific median concentrations ranging from 0.12 mg C L⁻¹ at the upper Val Ferret site (FEU) to 0.45 mg C L⁻¹ in the tributary stream at Vallon de Nant (RIC) (Fig. 2a, Table 2). All measured DOC concentrations (212 samples) were below 1.00 mg C L⁻¹. From simple linear regression, the median DOC concentration at a site varied most strongly with catchment vegetation cover (r = 0.76), δ¹³C-DIC values (r = −0.75), and catchment glacier cover (r = −0.53).

Figure 2 Box plots of (a) DOC and (b) DIC concentration (mg L⁻¹) from grab samples and of (c) CO₂ saturation (%) derived from sensor measurements.

Concentrations of DIC were generally greater and more varied than DOC, with an overall median DIC concentration of 1.77 mg C L⁻¹ across 191 samples, ranging between 0.79 and 2.65 mg C L⁻¹ (Fig. 2b). DIC concentration was most strongly correlated to decreasing mean catchment elevation (r = −0.67), with the three relatively high-elevation Valsorey locations exhibiting significantly lower median concentrations than the other nine sites (p < 0.01). The median δ¹³C-DIC value across sites was −6.14 ‰ (Table 2). The Champéry locations exhibited the most depleted δ¹³C-DIC values (median = −9.28 ‰), which were significantly lower than the remaining nine streams (p = 0.02).

Across all streams, the median saturation of CO₂ was 95.1 %, with the lowest median saturation of 68.1 % measured at the upstream location at Valsorey (VAU) and the highest median saturation of 137 % measured at the upstream location at Val Ferret (FEU) (Fig. 2c). All sites exhibited periods of oversaturation and undersaturation, except for VAU, where undersaturation was always observed. CO₂ saturation was significantly positively correlated with specific conductivity, alkalinity, DIC, and calcium, and negatively correlated with glacier coverage and specific UV absorbance at 254 nm (SUVA₂₅₄). However, the variance explained by any of these individual variables was low (r² < 0.3). A three-component PLS model was extracted that explained roughly 49 % of the variance in median CO₂ saturation (R²Y = 0.49), with moderate predictive power (Q² = 0.42). Ten variables were deemed highly influential (VIP > 1). These include catchment characteristics of mean catchment elevation, catchment area, glacier cover, and vegetation cover. Additionally, water quality parameters deemed influential were specific conductivity, sulfate and calcium concentration, total suspended solids, and discharge. Additionally, DOC was identified as a moderately influential variable.

Dissolved oxygen saturation was much less variable than CO₂ across sites, with median values between 98 % and 100 % and periods of over- and undersaturation for all sites (Table 2). Similarly, the interquartile range of CO₂ saturation across all sites was large (38.1 %) when compared with that of dissolved oxygen (2.3 %). The major cation across sites was Ca²⁺, and the major anion was ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Table S2 in the Supplement). The log ratios of Mg²⁺ and Ca²⁺ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ are similar across sites (Fig. 3), clustering closest to carbonate end-members (Torres et al., 2017). Within larger catchments, only the tributary site within the Val Ferret catchment (PEU) differs significantly from the main stem stream locations.

Figure 3 Stoichiometry of dissolved ions in the 12 study streams and from a global database of 95 glacier-fed streams (Torres et al., 2017). The range of each lithological end-member is shown by the boxes. The tributary stream in the Val Ferret catchment (PEU) is shown, as it is clearly distinguished from the main stream locations.

Figure 4 Intensity of the four components within the PARAFAC model against DOC concentration from grab samples, with catchment vegetation cover shown by color. (a) Component 1 and (b) component 2 represent humic-like compounds, whereas (c) component 3 and (d) component 4 represent proteinaceous compounds. The coefficient of determination (r²) is shown for each linear regression.

3.2 PARAFAC modeling results

PARAFAC modeling resulted in a four-component model (Fig. S1 in the Supplement). In comparing these components to the OpenFluor database, the first (C1) and second (C2) components are likely of terrestrial humic origin, while the third (C3) and fourth (C4) components are proteinaceous, likely of microbial origin (Kida et al., 2019). The components resemble those reported from other freshwater and glacial environments (e.g., Spencer et al., 2014; Imbeau and Vincent, 2021; Kida et al., 2021). When compared to EEM fluorophore peaks assigned by Coble et al. (1990, 1998), C1 appears to reflect the A and C peaks, which are associated with humic-like compounds from biodegradation of terrestrial plant matter, while C2 contains peak M, which is linked to humic-like compounds related to primary production. Similarly, C3 appears like the T peak and C4 appears like the B peak, both of which are suggested to be proteinaceous compounds of microbial origin. In general, the humic-associated components were found in greater intensity (median = 0.038 and 0.024 RU for C1 and C2, respectively) than the protein-associated components (median = 0.019 and 0.008 RU for C3 and C4, respectively). Both of the humic-associated components were significantly positively correlated with DOC concentration across all sites, C1 (r = 0.86) and C2 (r = 0.69) (Fig. 4). The protein-associated peaks showed little correlation with DOC concentration (r²≤0.2).

Figure 5 Estimated annual fluxes of the dissolved carbon components (CO₂, DOC, and DIC) normalized for catchment area. DOC and DIC are downstream fluxes, whereas CO₂ is a vertical flux. The DIC flux includes downstream transport of CO₂.

3.3 Catchment carbon fluxes

Total areal fluxes of dissolved carbon during the snow-free period ranged from −0.027 to 0.052 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ at the upstream Valsorey and downstream Champéry locations, respectively (Fig. 5). Considering absolute fluxes, the vertical flux of CO₂ was the largest component of the dissolved carbon flux, contributing a median of 67 %. The downstream flux of DIC (which includes dissolved CO₂) contributed 29 % to the total carbon flux, and DOC contributed the least (4 %). Negative net fluxes of carbon represent occasions when the stream is estimated to be a net sink of CO₂, and this sink exceeds the downstream transport of DOC and DIC. This occurred in only a single catchment (Valsorey).

4 Discussion

Comparing the dissolved carbon constituents in stream water within the space-for-time framework provided by these 12 study sites highlights how the changing nature of high-mountain catchments will have dramatic effects on the stream carbon cycle. There is a clear difference in DOC between higher-elevation and lower-elevation sites, likely as allochthonous carbon becomes more important with increasing vegetation cover at lower elevation. The saturation of CO₂ appears related to these DOC inputs, not only as a potential source of carbon for in-stream respiration but also as an indicator of an increasing importance of soil-derived CO₂ to the stream. Geochemical weathering remains a significant sink of CO₂, most strongly in glaciated catchments. However, the relevance of geochemical weathering to the CO₂ budget is not limited to glaciated catchments, as periods of undersaturation were observed in non-glaciated streams. Thus, the dissolved carbon dynamics of montane streams are critically tied to the dissolved carbon dynamics of high-mountain streams.

4.1 Increasing allochthonous DOC in high-mountain streams

The observed relationships between DOC concentration and catchment vegetation cover (and, even more strongly, the humic-like components of the DOM pool) suggest that allochthonous sources drive the increase in DOC concentration across these high-mountain streams. A higher stream DOC concentration has been attributed to greater terrestrial inputs and increasing vegetation cover (Zhou et al., 2019; Pain et al., 2020) as well as decreasing glacier influence (Fellman et al., 2010). The routing of water through catchment soils should thus play an increasingly large role in determining the timing and magnitude of allochthonous carbon export to high-mountain streams generally. For example, as the terrestrial environment becomes a more important source of DOC to streams, so too should hydrologic transport (Gómez-Gener et al., 2021). For example, in our study streams, DOC export has been shown to be strongly related to discharge patterns, with snowmelt mobilizing additional DOC compared with other seasons (Boix Canadell et al., 2019; Boyer et al., 1997). Similarly, rain events should then be related to increased humic-like DOM inputs from terrestrial sources, as transport from hillslope to stream is amplified (Fasching et al., 2016). While the indication of allochthony from DOM optical properties is imperfect (Begum et al., 2023; Guillemette and del Giorgio, 2012), the complementary DOC concentration patterns with vegetation cover and elevation reinforce the interpretation of these data.

Vegetation cover, as used in this study, serves as a broad indicator of soil development within these catchments, where accumulation of soil material allows for vegetation expansion (Hagedorn et al., 2019; Henne et al., 2011). The use of vegetation cover as a proxy for soil development following deglaciation is accurate in early successional stages (Klaar et al., 2014), which is true of the catchments in this work. The development of soil, as indicated by increasing vegetation cover, can increase the pool of organic carbon in glacier forelands (Wietrzyk-Pełka et al., 2020; Dümig et al., 2011; Egli et al., 2010), which can then be a source or organic carbon to the proglacial stream (Zah and Uehlinger, 2001). For example, in glacier-fed streams in Canada, the stream DOC concentration increased with catchment soil development (slope ≈ 0.2 mg C L⁻¹ per percentage of catchment area covered by soil; Lafrenière and Sharp, 2004), similar to our relationship with vegetation (slope = 0.05 mg C L⁻¹ per percentage of catchment area covered by vegetation). Either of these metrics, vegetation or soil, is indicative of significant catchment change with implications for terrestrial–aquatic carbon transfers.

Considering the greening of the terrestrial environment in the Alps (Rumpf et al., 2022), it follows that the streams draining these landscapes may be expected to experience an increase in the DOC concentration of terrestrial origin. Our results support this hypothesis, in which stream DOC concentration and the humic-like components likely of allochthonous origin increase with catchment vegetation cover. These changes have potentially important implications for these streams as well as their downstream ecosystems, such as altering metabolic regimes by promoting heterotrophy (Hall et al., 2016), limiting primary productivity (Kritzberg et al., 2019), and causing higher drinking water production costs (Hongve et al., 2004). Even while relatively low in concentration, the foundational physical, biochemical, and ecological nature of DOC within streams magnifies the impact of these changes in DOC concentration and highlights the substantial consequences of vegetation expansion following glacial retreat.

4.2 Terrestrial biogeochemical processes drive aquatic CO₂ saturation patterns

With regards to CO₂, extensive periods of undersaturation are relatively rare in riverine systems, but they are likely explained by geochemical weathering (St. Pierre et al., 2019). In our study, the isotopic signature of DIC provides the primary evidence of geochemical weathering, where depleted δ¹³C-DIC values (approximately −9 ‰ to −3 ‰) relative to atmospheric equilibrium are indicative of weathering (Skidmore et al., 2004). This agrees well with glacier-fed streams in Alaska (−7 ‰ to 0 ‰; St. Pierre et al., 2019), and mineral sources of DIC have been highlighted in Swiss high-mountain streams previously (Horgby et al., 2019c). Furthermore, the PLS model results also distinguish influential factors related to the products of weathering (i.e., specific conductivity and sulfate and calcium concentration) or which affect the rate of weathering (i.e., glacier cover, runoff, and total suspended solids). As such, the role of weathering in consuming CO₂ appears substantial.

The importance of geochemical weathering as a CO₂ sink in high-mountain areas is well described (Hilton and West, 2020; Donnini et al., 2016), where rapid weathering of carbonate and silicate rock consumes CO₂. In particular, elevated rates of weathering are expected for subglacial environments, where water flows over recently crushed, fine-grained reactive mineral surfaces (Tranter, 2003; Sharp et al., 1995). This process can continue in proglacial streams, where suspended sediments with high surface areas promote continued CO₂ drawdown (St. Pierre et al., 2019). Indeed, we see the lowest CO₂ saturation at the two most glacially influenced streams (VAU and VAD) within the Valsorey catchment. Thus, glacially enhanced weathering appears significant in this study as well.

Still, with periods of CO₂ undersaturation in all our study catchments, geochemical weathering appears to be relevant regardless of the presence of the glacier. To further constrain weathering as the primary sink of CO₂ in these catchments, we can also assess the potential for carbon fixation via photosynthesis as an alternative cause of undersaturation. With oxygen saturation consistently near or below saturation in all streams, photosynthesis is an unlikely driver of CO₂ undersaturation, as oxygen must inherently be above saturation to balance carbon fixation. Productivity has been shown to be limited in these streams outside of small temporal windows of opportunity (Boix Canadell et al., 2021), further reducing the likelihood. Lastly, the lack of variability in oxygen saturation across streams suggests that photosynthetic rates do not vary significantly across streams and, thus, cannot account for the observed variability in CO₂ saturation.

In contrast, variability within the DIC isotopic data does help explain the contribution of CO₂ to streams derived from the oxidation of organic matter in the terrestrial environment. The effect of organic carbon oxidation on δ¹³C-DIC values is depletion, i.e., more negative values (Pawellek and Veizer, 1994). Thus, it is likely that the depleted δ¹³C-DIC values observed at the Champéry streams are a result of greater rates of organic carbon oxidation, where the pool of organic carbon is evinced by the high vegetation cover and stream DOC concentration. We can more narrowly identify this process as most likely occurring in catchment soils, as the near-equilibrium nature of oxygen and the relatively low concentrations of DOC suggests a minor role for in-stream respiration (Bernhardt et al., 2018). Stream CO₂ is generally supported by external sources of CO₂, such as soil respiration (Hotchkiss et al., 2015; Campeau et al., 2019), and has been shown for mountain streams in particular (Clow et al., 2021; Crawford et al., 2015). Thus, as soils develop and organic carbon accumulates, the potential for terrestrially derived CO₂ inputs to the stream increases and CO₂ saturation increases (Marx et al., 2017). The role of the terrestrial environment in affecting stream CO₂ saturation is reinforced by the PLS model, which selected both vegetation cover and DOC concentration as influential variables. As such, there appears to be a link between increasing CO₂ saturation in these streams and organic matter accumulation and processing in the terrestrial environment.

Given the importance of geochemical weathering to the carbon dynamics of these catchments, consideration of geological variability between sites is necessary. While major lithologies are dissimilar across sites (Table 1), the ratios of major ions are remarkably similar and highlight the importance of carbonate weathering across all sites (Fig. 3). While carbonate-containing lithologies are clear in three of the catchments (Vallon de Nant, Champéry, Val Ferret), the carbonate-like signal in the Valsorey area is likely explained by high levels of calcite in the blue-gray schists (Bucher et al., 2017). This schist is primarily located beneath and in close proximity to the glacier (Burri et al., 1999); thus, glacier-enhanced weathering may disproportionately affect the weathering of this mineral. The prominence of carbonate weathering in these study streams may also indicate that the potential for geochemical weathering to serve as a CO₂ sink is elevated compared with glacier-fed streams globally (Torres et al., 2017). That is, catchments with lower proportions of carbonate-containing lithologies likely have lower potential as geochemical CO₂ sinks (St. Pierre et al., 2019). This elevated weathering as a result of carbonate-rich lithology is exemplified by the tributary stream in the Val Ferret catchment (Fig. 3). Despite no glacier coverage within this subcatchment, the median CO₂ saturation is undersaturated (Table 2). This is likely explained by the abundance of limestone deposits, which weather relatively quickly and geochemically consume CO₂. In expanding these analyses of carbon to other regions and mountain ranges, direct geologic perspectives will be needed to differentiate potential geochemical weathering rates (Hilton and West, 2020) and, hence, the potential for CO₂ consumption.

Figure 6 Conceptual model of processes affecting CO₂ saturation, and thus direction of flux, across glacier, soil, and elevation gradients within glacierized catchments. Geochemical weathering is important across the entire landscape but is enhanced under glaciated conditions and with proximity to the glacier. As vegetation and soil develop at lower elevation, terrestrial inputs add CO₂ through direct inputs from soil respiration and from organic carbon inputs which fuel in-stream respiration. The net balance of these processes determines the CO₂ saturation. In the aerial image of the Valsorey catchment, the transition from glacier to vegetation cover can be seen directly (from Google Earth©, 2023).

4.3 Conceptual model of carbon budgets in glacierized catchments

Altogether, these results provide the basis of a simple conceptual model elucidating contributions to stream CO₂, thereby explaining saturation dynamics across glacier, soil, and elevation gradients in mountain catchments (Fig. 6). Across the entire range of elevation, geochemical weathering acts as a sink of CO₂ (Crawford et al., 2019), where the intensity of this sink is dependent in large part on catchment geology. Where present, glaciers can provide additional weathering potential, whereby higher concentrations of suspended sediment increase mineral surface area greatly (St. Pierre et al., 2019). Moreover, this elevated weathering potential can extend far downstream depending on the suspension and transport of glacial till. Decreasing glacier influence reduces total weathering potential, but CO₂ undersaturation as a result of weathering is not limited to glacierized catchments. With the development of soils within the catchment, inputs of allochthonous organic carbon and CO₂ increases, elevating CO₂ concentrations. This CO₂ likely derives primarily from soil respiration rather than in-stream respiration of organic carbon (Clow et al., 2021; Singer et al., 2012).

Estimated fluxes of dissolved carbon constituents further support this conceptual model and the dominant role of terrestrial processes in determining the relative balance within and between streams. First, the dominance of CO₂ to the absolute total flux emphasizes the significance of gaseous carbon effluxes across within river networks. Our result of a 67 % contribution from CO₂ is similar to a study of a boreal catchment in Sweden, in which CO₂ accounted for 53 % of the net carbon flux (Wallin et al., 2013). Similarly, In a glaciated catchment in Alaska (St. Pierre et al., 2019), the areal rate of CO₂ flux was found to be −0.38 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, an order of magnitude higher than our most highly glaciated system (−0.03 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ at VAU). Following our conceptual model, the difference is explained by the much more heavily glaciated area of the Alaskan catchment (> 40 %) and the limitation of the sampling period to the most intense glacial melt period (June–August). When glacier influence is highest, the potential for weathering is highest as well, driving the consumption of CO₂. Yet, even without glacier influence, the consumption of CO₂ through weathering is still possible within the catchment and should be considered in montane stream carbon budgets.

DIC contributes 29 % to the total carbon flux, whereas DOC contributes the least, generally indicating a greater influence of mineral rather than organic processes (Rehn et al., 2022). The low contribution of DOC differs greatly from the Swedish boreal catchment, where DOC contributed roughly 40 % on average (Wallin et al., 2013). This difference not only highlights the limited soil development within high-mountain systems but also the potential for increased DOC export to the stream with continued soil development. Nonetheless, as the DOC contribution clearly increased with additional vegetation cover across our study systems, the role of the terrestrial landscape in supporting stream organic carbon content is clear.

Our focus on broad relationships across these 12 locations recognizably conceals how local conditions and seasonality may affect site-specific dynamics. Previous analyses have examined CO₂ and DOC individually at these stream locations, and they provide some perspective. On a finer spatial scale within the stream network, local groundwater inputs can disproportionately elevate the CO₂ concentration (Horgby et al., 2019b), which could be a useful tool in more directly quantifying terrestrial inputs to streams. Seasonally, the clearest pattern of pCO₂ indicates elevated contributions of terrestrial CO₂ during the spring snowmelt (Horgby et al., 2019a). Further differentiating temporal patterns over shorter timescales, such as storm events, may be useful in elucidating the contribution of soil respiration to the streams (Marzolf et al., 2022). While there is surely more to be learned at these finer spatial and temporal scales of both organic and inorganic carbon, our focus on broad-scale patterns across catchments allows us to make more generalizable conclusions. Given the strength of the observed relationships within our analyses and their consistency with other studies of high-mountain streams, our conceptual model provides a simple yet important foundation with which to assess carbon dynamics in montane streams globally.

5 Conclusion

The organic and inorganic components of the dissolved carbon pool shift across a glacier and vegetation gradient, driven by the relative balance of geochemical weathering and terrestrial carbon inputs to the stream. Our results also highlight an expanded importance of geochemical weathering in high-mountain ecosystems globally, whereby carbonate and silicate weathering may consume CO₂ across more mountain landscapes than previously considered (Horgby et al., 2019c). The implications for landscape carbon balances are clear, with decreased potential for CO₂ uptake and increased emissions of terrestrially derived CO₂ emerging after glacier retreat and landscape greening. The rate of the transition from carbon sink to source is likely accelerated by climate change (Knight and Harrison, 2014); thus, continued examination of the contributions of these processes to net stream balances is critical to predicting the future role of mountain catchments in the global carbon cycle.