Alpine Glacier Shrinkage Drives Shift in Dissolved Organic Carbon Export From Quasi‐Chemostasis to Transport Limitation

The export of dissolved organic carbon (DOC) from catchments is considered as an important energy flux through streams and a major connection between terrestrial and aquatic systems. However, the impact that predicted hydrological changes due to glacier retreat and reduction in snow cover changes will have on DOC export from high‐mountain streams remains unclear. In this study, we measured daily runoff and DOC yield during 1 year in Alpine streams draining catchments with different levels of glacier coverage. DOC yield showed a varied response to runoff across the catchments and varied seasonally as a function of the degree of glaciation and vegetation cover. Using space‐for‐time substitution, our results indicate that the controls on DOC yield from Alpine catchments change from chemostasis to transport limitation as glaciers shrink.


Introduction
Mountain glaciers are shrinking worldwide owing to global warming. As glaciers shrink, runoff generated from ice melt initially increases until a peak is reached, beyond which it diminishes until the ice mass has disappeared (Huss & Hock, 2018). The potential effects of glacier shrinkage on glacier-fed streams have been studied focusing on alterations in hydrological regimes, water availability, geomorphology, and biodiversity (Bakker et al., 2018;Lane et al., 2017;Lane & Nienow, 2019;Milner et al., 2017). Glacier runoff carries sediments and solutes, including inorganic nutrients and organic carbon, and its change is also anticipated to affect the biogeochemistry of glacier-fed streams (Milner et al., 2017).
Dissolved organic carbon (DOC) is a significant component of the land-to-ocean carbon flux (Regnier et al., 2013) and constitutes a large source of energy to stream ecosystem metabolism (Battin et al., 2009). Within streams, DOC can be metabolized leading to respiratory CO 2 release to the atmosphere, transiently stored or exported downstream, and ultimately to the ocean. Below the tree line, terrestrial vegetation and soils constitute a major source of DOC that enters streams via various hydrological flow paths, often during rainfall or snowmelt (e.g., Boyer et al., 1997;Lafreniere & Sharp, 2004). In high-mountain catchments above the tree line, alternative sources of DOC may gain importance because of poorly developed soils. For instance, mountain glaciers can release significant amounts of organic carbon to downstream ecosystems Li et al., 2018), where it can even partially sustain stream food webs .
The transport of DOC across the terrestrial-aquatic interface is thus relevant for stream ecosystem functioning and its export from the catchment is critical for the carbon balance at catchment scale (Tank et al., 2018).
Despite this, our ability to predict the impact of climate change on DOC export is relatively limited at present (Zarnetske et al., 2018). DOC export from catchments can be expressed as the load (mass/time) or yield (mass/area/time). While both quantities are relevant for stream carbon budgets, the comparison of DOC yields between different streams can identify the effects of catchment properties (e.g., drainage area and land use) upon DOC export. Furthermore, the relationship between DOC yield and discharge indicates if the DOC flux is limited by hydrological transport or by sources of DOC within the catchment, or if it is in chemostasis (Creed et al., 2015;Vaughan et al., 2017;Zarnetske et al., 2018). Elucidating these potential drivers is relevant given the expected climate-induced alterations of the hydrological regime of streams and changing land cover.
Despite the understanding that glacier shrinkage may affect the functioning of glacier-fed streams and even downstream biogeochemistry, the drivers of DOC export from glacierized catchments have been rarely studied (e.g., Hood & Scott, 2008;Tockner et al., 2002). Most of our current understanding of DOC export from high-altitude catchments rests on studies of snow-covered catchments and snowmelt, where the latter was reported as a major contributor to annual DOC export (Ågren et al., 2010;Boyer et al., 1997;Hood et al., 2003;Lafreniere & Sharp, 2004). Both ice melt and snowmelt shape the hydrological regime of glacier-fed streams, but the relevance of these periods for annual DOC fluxes is poorly understood. This is also because the ice melt itself can be a DOC source to the stream Li et al., 2018).
The aim of this study was to elucidate how runoff, particularly as mediated by snowmelt and glacier ice melt, affects the DOC yield from non-glacierized and glacierized catchments. We used high-frequency (every 10 min) measurements of discharge, streamwater temperature, turbidity, and FDOM over 1 year to infer DOC concentrations and yields from 12 streams distributed over four catchments in the Swiss Alps. This allowed us to use a space-for-time substitution approach to show that the drivers of DOC yield may shift from quasi-chemostasis to hydrological transport limitation in alpine catchments as glaciers shrink.

Study Sites
We selected six streams that are not fed by glaciers (that is, draining non-glacierized catchments: La Vièze, Z1, Z2, and Z3; Torrent de la Peule, P1; Tributary of Torrent du Valsorey, T1; Ruisseau du Richard, R1) and six glacier-fed streams under the influence of glacier ice melt dynamics (Avançon de Nant, N1 and N2; Dranse de Ferret, F1 and F2; Torrent du Valsorey, V1 and V2) in the Swiss Alps ( Figure S1 and Text S1 in the supporting information). With this sampling design, we covered non-glacierized and glacierized catchments where the latter had glacier coverage ranging from 2% to 28%. Across these glacierized catchments, glacier coverage has decreased by 35% to 100% since the little ice age in 1850 (Table 1). This selection allowed a space-for-time substitution approach as often done to assess downstream implications of glacier shrinkage (Heckmann et al., 2016). The catchments (1,778 to 2,893 m above sea level, a.s.l.) are dominated by steep slopes and sparsely vegetated (grassland and shrubs) rocky terrain (Table 1 and Table S1). Coniferous forest is present at the Avançon de Nant and the Ruisseau du Richard with less than 16% coverage. Mixed forest is present at the downstream site of La Vièze (Z1, 15% coverage). Spring runoff is dominated by snowmelt in non-glacierized catchments and the combination of both snow and ice melt in glacierized catchments, peaking in spring and summer, respectively. Land cover was assessed based on CORINE© Land Cover Inventory 2012 of the European Union, Copernicus Land Monitoring Service.

Sensor Network and In Situ Measurements
We equipped all 12 streams with sensors measuring streamwater chromophoric DOC (Cyclops-7 CDOM/FDOM, Turner Designs), turbidity (Cyclops-7 Turbidity, Turner Designs), water depth (Odyssey Capacitance Logger, Dataflow Systems Ltd), and temperature. Optical sensors were calibrated in the laboratory before deployment following the manufacturers' specifications. Turbidity sensors were calibrated to nephelometric turbidity units (NTU) using sensor-specific Turbidity 4000 NTU Calibration Standard-Formazin (Sigma-Aldrich®) up to 2000 NTU. Excitation and emission wavelengths of the FDOM fluorometers were centered at 325 (band pass: 120) nm and 470 (band pass: 60) nm, respectively. Such sensors are now being increasingly used to study FDOM, and by conversion DOC concentration, at high temporal resolution in various freshwater ecosystems (e.g., Pellerin et al., 2012, Koenig et al., 2017Tissier et al., 2013;Tunaley et al., 2016;Wilson et al., 2013; Text S3). FDOM measurements were converted from signal voltage to parts per billion (ppb) of quinine sulfate equivalents (fluorescence of 1 ppb quinine sulfate monohydrate (Sigma-Aldrich®) in 0.05 M H 2 SO 4 ) using a two-point calibration curve (0 to 1,000 ppb quinine sulfate equivalents).
Streamwater FDOM, temperature, turbidity and depth were measured at 10-min intervals from 1 October 2016 to the 30 September 2017. We visited all sites at least monthly for sensor maintenance, data downloading, and grab sampling for DOC analysis. Streamwater samples were filtered (precombusted GF/F filters, Whatman) into acid-washed, precombusted glass vials and analyzed using a Sievers M5310c TOC Analyzer (GE Analytical Instruments) (accuracy: ±2%, precision: <1%, detection limit: 22 μg C/L).
Discharge was determined using slug injections of sodium chloride as a conservative tracer (Gordon et al., 2004) and rating curves between streamwater depth and discharge were established for each site (Text S2). Where streamwater depth data were missing because of sensor malfunctioning, we inferred data from gauging stations within the respective catchments (Text S2). We derived instantaneous daily discharge (m 3 /s) from the 10-min frequency time series averaged to daily means. We filled data gaps in daily discharge data using linear interpolation and estimated daily runoff (mm/day).  (168) Note. Bare rocks, glaciers and perpetual snow, and vegetated cover percentage are based on CORINE Land Cover Inventory 2012. Percent vegetation cover in brackets refers to mixed and coniferous forest, moors, and heathlands. Remaining vegetation types include pastures, sparsely vegetated areas, and natural grassland. Glacier coverage in 1850 is from Maisch et al. (2000). Shown are means ± standard deviation. Conservative estimates of dissolved organic carbon yields include in brackets the number of days for which the yield was determined.

Relationships Between Streamwater FDOM and DOC Concentrations
We corrected sensor FDOM values for streamwater temperature (FDOM temp ) and turbidity (FDOM corr ) as follows: where FDOM raw is the uncorrected FDOM, Temp m is streamwater temperature, Temp ref is the reference temperature at 10°C, and Turb m is streamwater turbidity (Lee et al., 2015;Watras et al., 2011). We derived the coefficients from laboratory experiments on the sensor-specific effect of a range of temperatures and turbidity on the FDOM signal. Coefficients were comparable with those previously reported (Downing et al., 2012;Lee et al., 2015;Saraceno et al., 2017;Watras et al., 2011). Time series of streamwater DOC concentrations were computed from least-squares linear regressions between FDOM corr sensor data and DOC concentrations measured from grab samples (Text S3).

Calculations of DOC Loads and Yields
We modeled the daily load of DOC (kg C/day) in each stream using LOADEST (Runkel et al., 2004), normalized to catchment area to derive daily DOC yield (kg C·km −2 ·day −1 ). LOADEST uses time series of daily constituent concentration data and the corresponding measurement of stream discharge to build a calibration regression, which is then applied to a complete data series of daily discharge to obtain daily constituent loads (mass/day: Runkel et al., 2004). We created 12 site-specific calibration input files with the available pairs of daily DOC concentration and discharge (Text S4). We derived DOC concentration data from the 10-min frequency time series averaged to daily means. Site-specific regressions were validated using an Adjusted Maximum Likelihood Estimator. We selected the best fit of the predefined regression models available based on Akaike's Information Criterion.
We provide the 95% confidence interval on the predicted DOC yield as provided by LOADEST as a measure of the uncertainty associated with the DOC yields. However, LOADEST does not propagate the error related to discharge rating curves and FDOM/DOC relationships. Therefore, we also calculated DOC yield using the "conservative" approach (i.e., Load = [DOC] × Q) for those days where DOC concentration ([DOC]) and discharge (Q) were available (Text S5). We estimated the uncertainty in DOC yield using a Gaussian error propagation. Calculating these uncertainties makes particular sense for the streams with almost complete time series of both discharge and DOC.

Data Analysis
We approximated hydrologic periods (i.e., winter baseflow, spring snowmelt, summer ice melt, and fall baseflow) typical for alpine streams by partitioning discharge into baseflow and high flows using a recursive digital filter approach and the EcoHydRology package in R (Fuka et al., 2018). The value of the filter parameter (or the hydrograph recession constant) was set at 0.925 (Nathan & McMahon, 1990). This step was followed by manual adjustment of the periods because the onset of the winter baseflow and snowmelt, for instance, vary with altitude.
We quantified the relationships between runoff and DOC yield, rather than DOC concentration, as these are well suited to evaluate the carbon balance at catchment level (Zarnetske et al., 2018). We used a power function to describe the relationship between runoff and DOC yield for each stream throughout the study period of 1 year, where the slope coefficient (b) refers to chemostasis, source or transport limitation of DOC. Chemostasis (b = 1) occurs when the mass of DOC delivered to the stream can keep up with rising discharge; source limitation (b < 1) occurs when DOC in the stream is being diluted due to an inability of the catchment to produce more DOC as discharge rises; transport limitation (b > 1) occurs when new DOC sources are mobilized as discharge increases (e.g., Creed et al., 2015;Zarnetske et al., 2018). The intercept, a, is a DOC concentration normalization coefficient (Zarnetske et al., 2018). We transformed (natural log) the power model and calculated least squares linear regression models on the data to retain the intercept (a) and slope coefficient (b). The goodness of fit and the parameter estimation significance were assessed based on the coefficient of determination (R 2 ) and the P value. As in other studies (e.g., Vaughan et al., 2017;Zarnetske et al., 2018), our models did not account for potential temporal autocorrelation of the data.
We evaluated the effect of runoff on DOC yield across all streams at a daily basis using a linear mixed-effect model with random-effect terms. The linear mixed-effect model was fitted using the lme function (Pinheiro et al., 2018) and included runoff and streamwater temperature as fixed effects and site as random effect (with random intercept and slope). Here we used the temporal correlation structure defined by the residuals (autoregressive process of order 1, or AR(1)). We evaluated the goodness of fit of the model by calculating marginal and conditional R 2 using the r.squaredGLMM function (Barton, 2018). All statistical analyses were conducted in R (version 3.5.1, R Development Core Team, 2018) and on DOC yield estimates from LOADEST.

Runoff Dynamics
Runoff in high-alpine streams is influenced by altitude and latitude, glacier mass balance, and snowpack (Milner et al., 2017) and can vary from seasonal (e.g., snow melt in spring) to daily (e.g., ice melt in summer) scales. Because runoff is a master variable for stream biogeochemistry, we first describe its temporal patterns in our study catchments.
Peak runoff started with snowmelt in spring (glacierized: 748 ± 118 mm, median: 733 mm; non-glacierized: 539 ± 249 mm; median: 512 mm) and contributed on average 58 ± 2% and 58 ± 10% to the annual runoff from glacierized and non-glacierized catchments, respectively. The continuation of peak runoff into summer was particularly pronounced in glacierized catchments where summer runoff (338 ± 83 mm; median: 320 mm) contributed on average 26 ± 5% to the annual runoff compared to 22 ± 8% in non-glacierized (237 ± 199 mm; median: 209 mm) catchments. Summer runoff increased significantly (R 2 = 0.74, P < 0.05) with glacier coverage across the six glacier-fed streams. These differences are likely due to the combined effects of the larger snowpack in glacierized catchments owing to their higher elevations and glacier runoff itself. In fact, it is well established that annual runoff in glacierized catchments can increase due to elevated rates of glacier volume loss (Hood & Scott, 2008;Milner et al., 2017).

Temporal and Spatial Patterns of DOC Yields
First, we explored the partitioning of the annual DOC yield over the hydrological periods. The average contribution of the DOC yield during snowmelt to the total annual yield was 62 ± 3% and 62 ± 13% for streams in glacierized and non-glacierized catchments, respectively (Table S3). This shows that more than half of the total annual DOC export occurred during spring in these high-alpine catchments. Spring contributions to DOC yield tended to be higher (61.7 ± 9.1%) than the equivalent contribution to runoff (57.8 ± 7%), which would imply that reservoirs of organic carbon are being mobilized during snowmelt across all sites. This is in 10.1029/2019GL083424

Geophysical Research Letters
line with the DOC flushing from near-surface soil horizons observed in an alpine catchment in Colorado (Boyer et al., 1997). Snowmelt was also reported to be critical for the DOC export from high-latitude catchments suggesting that carbon production during the cold season is important for regulating arctic DOC transport (Finlay et al., 2006). We evoke a similar relevance of DOC production during winter for high-alpine streams, most likely originating from the sparse vegetation from the antecedent summer degrading underneath the snowpack as reported from high-elevation catchments in the Rocky Mountains (Brooks et al., 1999). Furthermore, direct leaching of DOC from terrestrial plant material contained within the snowpack (Lafreniere & Sharp, 2004) or atmospheric deposition (Mladenov et al., 2012) may also contribute to the streamwater DOC..
DOC yield during summer contributed on average 21 ± 3% and 22 ± 11% to the annual DOC yield from glacierized and non-glacierized catchments, respectively. In the glacier-fed streams, DOC yield increased with the degree of glacier coverage (R 2 = 0.92, P < 0.01) during summer, which is likely attributable to elevated glacier runoff. DOC yield during autumn and winter contributed on average 8 ± 5% and 9 ± 8% to the annual DOC yield from glacierized and non-glacierized catchments. We suggest that this is attributable to the DOC poor groundwater predominantly feeding high-alpine streams during these periods (e.g., Milner et al., 2017).

Drivers of DOC Yields
For most of our study streams, the power function describing the relationship between runoff and DOC yield resulted in a slope coefficient (b) higher than 1. Interestingly, b values in non-glacierized streams (0.99 ± 0.05 to 1.51 ± 0.03) were significantly (t test: P < 0.05) greater than 1, whereas b values in glacier-fed streams (0.92 ± 0.02 to 1.13 ± 0.03) were not significantly (t test: P > 0.05) greater than 1 (Figure 2a and Figure S2). These results indicate that DOC yield in non-glacierized catchments was largely controlled by transport limitation. Peak runoff in these streams occurred during snowmelt and sometimes in summer during storm events (i.e., in Z2 and Z3; Figure 1). Overall values of b in non-glacierized catchments show that the hydrological mobilization of DOC and its supply to the stream was guaranteed and likely enhanced by the activation of terrestrial DOC reservoirs as runoff increases during these events. On the other hand, b values closer to 1 for glacier-fed streams suggest chemostasis where DOC production and mobilization are nearly proportional to changes in runoff (Godsey et al., 2009). Additionally, ice melt from glaciers in summer contributes directly to stream runoff, and concomitantly and proportionally to the DOC yield in glacier-fed streams. In fact, upon melting, ice-locked DOC can contribute to the downstream DOC pool Li et al., 2018;Singer et al., 2012). We argue that ice melt as a common source of water and DOC underlies the quasichemostasis of the DOC yield in glacier-fed streams.

Geophysical Research Letters
The slope coefficient, b, was positively related to vegetation coverage, including mixed or coniferous forest, moors, and heathlands, typically indicative of later successional stages following glacier recession in alpine catchments (Figure 2b). In addition, the intercept, a, was positively related to vegetation coverage including all types (e.g., grass land) of vegetation (R 2 = 0.42; P < 0.05), further supporting the notion of DOC source availability but transport limitation. We found the highest b (1.51 ± 0.03) associated with the spring-fed stream T1 where sources of DOC may be contained within fringing wetlands. This is in agreement with the observation that wetland coverage increases the likelihood of transport limitation in a wide range of catchments (Zarnetske et al., 2018).
Vegetation coverage was highest in the catchments that experienced highest glacier loss by area since 1850 (R 2 = 0.57; P < 0.05) ( Figure S3) and we assume, therefore, that the buildup of soils as a source of DOC to the streams was higher in these catchments. It is reasonable to assume that glacierized catchments are out of equilibrium in terms of the relationships between climate and soil development as they have undergone ice retreat over the last 120 years at least. Because soil development on deglaciated terrain is a complex function of temperature, and hence altitude, moisture limitation and weathering (Miller & Lane, 2019), for instance, we argue that our observed patterns of DOC yield between glacierized and non-glacierized catchments may not be exclusively linked to the degree of glacier coverage. Rather, they could be linked to the extent to which the surface regolith has weathered to allow soil and associated organic matter stocks to develop as glaciers recede.
The streams (Z1, Z2, and Z3) draining the Champéry catchment had highest annual DOC yields and, along with T1, the highest slope coefficient b as well. The Champéry catchment has not been glacierized since the little ice age at least and could be considered as an end-member toward which the presently glacierized catchments may evolve in terms of DOC yield.
This assumes that the increase in streamwater DOC concentration with the encroachment of vegetation has a stronger positive effect on DOC yields than the decrease in runoff associated with glacier shrinkage. We acknowledge that the local climate, changing with altitude and exposure, may interfere with this relationship as indicated by the DOC yield in P1, also glacier-free since the little ice age at least. In fact, the difference between the P1 and the streams in the Champéry catchment may be attributable to the differences in altitude between these systems and its impact on climate and hence on the development of vegetation and soils.
In overall agreement with the power function analyses, the mixed effect model, including a temporal autocorrelation structure, revealed higher slopes for DOC yield and runoff relationships for streams draining nonglacierized than glacierized catchments (Tables S4 and S5). We recognize that the results from the power function analyses and the linear mixed effect model are not directly comparable. Nevertheless, we interpret the similar results from both analyses as evidence of different drivers of DOC yield in glacierized and non-glacierized catchments.

Contextualizing Annual DOC Yields From High-Alpine Catchments
LOADEST and our conservative approach yielded reasonably close estimates of annual DOC yields, at least for the sites with rather complete time series (Table 1). Uncertainty associated with the second approach was relatively high for some streams, which we attribute to the error related to discharge measurements. This is not unexpected though as it is inherently difficult to get proper measurements of high discharge in steep mountain channels with high roughness. Nevertheless, our estimates are closely bracketed by those from similar systems or lower than those from more productive systems ( Figure S6).
The annual DOC yield averaged 189 ± 46 kg C·km −2 ·year −1 (median: 186 kg C·km −2 ·year −1 ) in glacierized and 401 ± 308 kg C·km −2 ·year −1 (median: 253 kg C·km −2 ·year −1 ) in non-glacierized catchments ( Table 1). As expected from our analysis of DOC yield at the daily scale, the response of annual DOC yields to runoff also differed between glacierized and non-glacierized catchments (Figure 2c). In fact, per unit annual runoff, the annual DOC yield from non-glacierized catchments was on average 3.32 times larger than the DOC yield from glacierized catchments. These differing relationships are comparable to those reported from catchments with and without wetlands (Mulholland, 1997), for instance. They emphasize the various sources (e.g., soils, glacier ice) of organic carbon and their relative extent in glacierized and non-glacierized catchments. Annual DOC yields estimated in this study are lower than those reported from catchments with wetlands and forests but comparable to those from temperate grasslands (Hope et al., 1994). For instance, the DOC yield for Z2 (937 kg C·km −2 ·year −1 ) in the La Vièze catchment with its high vegetation cover is close to the yield reported from a high-altitude stream in the Rocky Mountains (960 kg C·km −2 ·year −1 ; Pacific et al., 2010). On the other hand, the annual DOC yields from our glacierized catchments are not comparable to the high DOC yields from Alaskan catchments, likely attributable to the very high runoff from Alaskan glaciers (Hood & Scott, 2008).

Conclusions
Mountain glaciers are shrinking rapidly worldwide, but impacts on downstream biogeochemistry, such as on the DOC export fluxes, remain poorly understood (Milner et al., 2017). Using a space-for-time substitution approach, we found that a quasi-chemostasis of DOC yield in glacier-fed streams may shift to a limitation by hydrologic transport as glaciers shrink. Underlying processes are complex and likely entail shifts at the level of catchment hydrology, soil development, and weathering. For instance, we propose that as glaciers shrink, runoff production from ice melt directly feeding into the stream channel will diminish, whereas hydrological flow paths through catchments may diversify potentially delivering DOC from various sources to the streams. At the same time, soil organic carbon and its spatial distribution within catchments will be impacted by changes in plant communities and soil properties due to warmer temperatures (Mayor et al., 2017). Shifts in the seasonal patterns of runoff and DOC yield will impact the timing and magnitude of the lateral DOC fluxes from terrestrial to aquatic ecosystems, and further to downstream ecosystems.