Projecting sediment export from two highly glacierized alpine catchments under climate change: exploring non-parametric regression as an analysis tool

Future changes in suspended sediment export from deglaciating high-alpine catchments affect downstream hydropower reservoirs, flood hazard, ecosystems and water quality. Yet, quantitative projections of future sediment export have so far been hindered by the lack of process-based models that can take into account all relevant processes within the complex systems determining sediment dynamics at the catchment scale. As a promising alternative, machine-learning (ML) approaches have recently been successfully applied to modeling suspended sediment yields (SSYs).

This study is the first, to our knowledge, exploring a machine-learning approach to derive sediment export projections until the year 2100. We employ quantile regression forest (QRF), which proved to be a powerful method to model past SSYs in previous studies, for two nested glaciated high-alpine catchments in the Ötztal, Austria, above gauge Vent (98.1 km²) and gauge Vernagt (11.4 km²). As predictors, we use temperature and precipitation projections (EURO-CORDEX) and discharge projections (AMUNDSEN physically based hydroclimatological and snow model) for the two gauges. We address uncertainties associated with the known limitation of QRF that underestimates can be expected if values in the projection period exceed the range represented in the training data (out-of-observation-range days, OOOR). For this, we assess the frequency and extent of these exceedances and the sensitivity of the resulting mean annual suspended sediment concentration (SSC) estimates. We examine the resulting SSY projections for trends, the estimated timing of peak sediment and changes in the seasonal distribution.

Our results show that the uncertainties associated with the OOOR data points are small before 2070 (max. 3 % change in estimated mean annual SSC). Results after 2070 have to be treated more cautiously as OOOR data points occur more frequently, and glaciers are projected to have (nearly) vanished by then in some projections, which likely substantially alters sediment dynamics in the area. The resulting projections suggest decreasing sediment export at both gauges in the coming decades, regardless of the emission scenario, which implies that peak sediment has already passed or is underway. This is linked to substantial decreases in discharge volumes, especially during the glacier melt phase in late summer, as a result of increasing temperatures and thus shrinking glaciers. Nevertheless, high(er) annual yields can occur in response to heavy summer precipitation, and both developments would need to be considered in managing sediments, as well as e.g., flood hazard. While we chose the predictors to act as proxies for sediment-relevant processes, future studies are encouraged to try and include geomorphological changes more explicitly, e.g., changes in connectivity, landsliding, rockfalls or vegetation colonization, as these could improve the reliability of the projections.

1 Introduction

Fluvial suspended sediment export from glacierized, high-alpine areas can be up to 1 order of magnitude higher (per unit area) than in non-glacierized downstream areas (Hinderer et al., 2013; Beniston et al., 2018). Thus, sediment dynamics in these high-alpine areas and changes therein have important implications for downstream hydropower generation and reservoir sedimentation (Schöber and Hofer, 2018; Guillén-Ludeña et al., 2018; Li et al., 2022), water quality (as well as nutrient and contaminant transport) (Bilotta and Brazier, 2008), aquatic species and riverine ecosystems (Milner et al., 2009, 2017; Gabbud and Lane, 2016; Huss et al., 2017), and also flood hazard (Nones, 2019) and carbon cycling (Tan et al., 2017; Syvitski et al., 2022).

High-alpine areas are particularly sensitive to climate change, experiencing above-average warming (Gobiet et al., 2014) and hence crucial cryospheric changes, such as ongoing and accelerating deglaciation, permafrost melt, and snow cover changes (Huss et al., 2017; Beniston et al., 2018; Abermann et al., 2009). These changes induce changes in discharge volumes, timing and magnitude (Vormoor et al., 2015; Kuhn et al., 2016; van Tiel et al., 2019; Rottler et al., 2020; Hanus et al., 2021). This in turn affects sediment export, and past changes have been observed frequently, e.g., due to enhanced subglacial sediment evacuation and increased sediment availability in expanding erodible landscapes as receding glaciers and melting permafrost e.g., expose glacial till and weaken rockwalls (Micheletti and Lane, 2016; Carrivick and Heckmann, 2017; Lane et al., 2017, 2019; Costa et al., 2018; Delaney and Adhikari, 2020; Li et al., 2020; Vergara et al., 2022).

Nevertheless, future changes in sediment export are understudied (Zhang et al., 2022), and questions such as “Are sediment yields from deglaciating catchments increasing, decreasing or is there no pattern?” or “to what extent is it possible to quantify spatio-temporal patterns of future sediment yields?” (Carrivick and Tweed, 2021) have yet to be answered – although projections of climatological (e.g., Gobiet and Kotlarski, 2020; Gobiet et al., 2014), glaciological (e.g., Stoll et al., 2020; Bolibar et al., 2022; Huss, 2011) and hydrological changes (e.g., Madsen et al., 2014; Hanzer et al., 2018; Hanus et al., 2021; Huss and Hock, 2018; Tecklenburg et al., 2012; Wijngaard et al., 2016), which could serve as a basis for estimating future changes in sediment export, are numerous.

The main reason why answering such questions is challenging is that modeling sediment export at the catchment scale with process-based models remains difficult – if not impossible – because it is determined by a complex system of interconnected processes that is not straightforward to model. For example, the relationship between suspended sediment concentrations and discharge is most often nonlinear in time and space, and univariate models relying solely on discharge are often insufficient (Vercruysse et al., 2017; Zhang et al., 2021). Hence, in addition to variations in discharge, changes in sediment availability, entrainment, transport and deposition would have to be considered; there may be threshold effects and nonlinear responses of geomorphic processes (e.g., triggering of mass movements), correlated influencing factors, hysteresis and seasonality (Huggel et al., 2012; Landers and Sturm, 2013; Vercruysse et al., 2017; Costa et al., 2018; Schmidt et al., 2023; Zhang et al., 2022). Additionally, long-term field observations (i.e., several decades and covering a wide range of conditions) that provide enough training and validation data to develop sediment yield models or to analyze trends are very rare (Zhang et al., 2022; Schmidt et al., 2023).

There are conceptual models on (suspended) sediment export from deglaciating areas (Antoniazza and Lane, 2021; Carrivick and Tweed, 2021; Zhang et al., 2022), which expect an initial increase in sediment export as glaciers begin to retreat and an eventual decrease – after peak sediment – once the glaciers have disappeared and the landscape stabilizes. The timing of peak sediment relative to peak meltwater and the completion of deglaciation is presumed to depend on changes in erosive precipitation; i.e., a negative trend in erosive precipitation implies that peak meltwater and peak sediment may co-occur, while a positive trend or no change in erosive precipitation result in a lag between peak meltwater and peak sediment. However, deducing estimates of future sediment export and implications for individual catchments based on these conceptual considerations is not straightforward or even possible.

Common approaches to model sediment yields at the catchment scale, such as SWAT (e.g., Vigiak et al., 2017), BQART (Syvitski and Milliman, 2007), WBMsed (Cohen et al., 2013), WASA-SED (Mueller et al., 2010) or SAT (Zhang et al., 2021), are mostly empirical or conceptual in their sediment modules, do not consider all relevant erosion processes (i.e., neglecting glacial, gully erosion and landslides in the case of SWAT), and often concentrate on large spatial scales (i.e., sediment fluxes to the oceans for large basins, entire continents or at global scale) and/or large temporal scales (i.e., multiyear averages and long-term fluxes). On the other end of the spectrum, models for individual parts or processes within glacierized catchments exist as, for example, a numerical approach to model subglacial fluvial sediment transport (Delaney et al., 2019) that has also been coupled with models for ice dynamics and bedrock erosion (Delaney et al., 2021) or, e.g., probabilistic or physical models of mass wasting processes, such as landslide or debris flows (Iverson and George, 2014; Hirschberg et al., 2021; Campforts et al., 2022). However, as of yet, there is no all-in-one physical model (fully distributed, incorporating thermal and pluvial drivers of sediment mobilization and transport) to simulate sediment export from cryospheric basins (Zhang et al., 2022) at the catchment scale.

Accordingly, studies that have attempted to project future suspended sediment yields (SSYs) chose rather qualitative approaches, such as comparing sediment yield observations of warmer and colder ablation seasons (Stott and Mount, 2007; Bogen, 2008) and using the responses of SSYs to past predictor changes and applying these to projected changes in the future (Li et al., 2021b) or fitting a multiple regression model to past data (of only 1 year) and increasing the temperature input in the model (Stott and Convey, 2021). However, these approaches may preclude modeling decreases or accounting for interactions between variables.

As a promising alternative, geoscientific machine-learning approaches have emerged and have recently been acknowledged for their potential in applications to Earth system science (Reichstein et al., 2019). Indeed, first studies showed that machine-learning approaches can easily outperform well-known existing models for sediment yield (Gupta et al., 2021; Rahman et al., 2022; Jimeno-Sáez et al., 2022; Schmidt et al., 2023). Such black-box approaches tend to perform well for black-box problems such as high-alpine sediment dynamics, “where the input data and output data are well understood or at least fairly simple, yet the process that relates the input to output is extremely complex” (Lantz, 2019). In a previous study, we have developed and validated a quantile regression forest (QRF) approach to model SSYs in two nested high-alpine catchments and to estimate yields for the past 5 decades (Schmidt et al., 2023). This showed that the QRF model outperformed commonly applied sediment rating curves by about 20 % of explained variance, and other studies found that regression trees and random forest models (which QRF is based on) even outperformed other machine learning approaches in modeling sediment dynamics (Talebi et al., 2017; Al-Mukhtar, 2019).

Thus, the present study is motivated to explore QRF to model future SSYs based on measurement data, emission scenarios and subsequent hydrological model results. We test the approach in two glacierized high-alpine catchments in the Ötztal in Austria, where projections of future climatological and glacio-hydrological conditions from the AMUNDSEN model are available (Hanzer et al., 2018) and where we have successfully trained and applied QRF models to reconstruct past sediment export using records of discharge, precipitation and air temperature (Schmidt et al., 2023).

The goals of the present study are (i) to assess uncertainties of the model due to known limitations of the QRF method in order to identify the limitations of the approach and (ii) to derive estimates of future changes in sediment export with respect to trends in annual yields, shifts in the seasonal distribution and the timing of peak sediment.

2 Methods

In a previous study, we trained and validated quantile regression forest models to retrospectively estimate SSYs at two gauges for the past 5 decades using the available records of turbidity-derived suspended sediment concentrations (4 and 15 years) and long-term records of the predictors, i.e., discharge, precipitation and temperature (Schmidt et al., 2023) (Fig. 2, dashed-line box). In the present study, we use these models and apply them to downscaled and bias-corrected EURO-CORDEX temperature and precipitation projections that were used as input data for the glacio-hydrological model AMUNDSEN, as well as the discharge projections of AMUNDSEN (Hanzer et al., 2018) (Fig. 1). In the following, we outline the quantile regression forest approach including its advantages and limitations with respect to modeling suspended sediment dynamics and the choice of predictors to model sediment dynamics in high-alpine areas. Then, we describe the study area, input data and necessary adjustments, as well as how we analyzed the limitations, sensitivities and the resulting SSY estimates.

Figure 1 Overview of models and resulting projections used in this study. Bias-corrected EURO-CORDEX climate projections and AMUNDSEN model results serve as input data for the QRF models. Q: discharge, P: precipitation, T: temperature, SSC: suspended sediment concentrations, SSY: suspended sediment yields.

Figure 2 Temporal extent of input data, as well as modeled suspended sediment yields (SSYs) in the previous study (dashed-line box, top left) and projections (discharge Q, precipitation P, temperature T) as input data and SSY estimates from this study; 2006–2020 is the training data period at gauge Vent. At gauge Vernagt, the QRF model was trained on the years 2000, 2001, 2019 and 2020, when SSC data were available. Thus, to verify model results in the present study, we use the QRF estimated yields at gauge Vernagt of the years overlapping with the climate and hydrology projections, i.e., the overlap period of 2007–2020 (Sect. 2.4.1).

2.1 Quantile regression forest for suspended sediment concentration modeling

The quantile regression forest (QRF) (Meinshausen, 2006) is a non-parametric regression technique that is based on random forests (RFs) and can be classified as a machine-learning approach. It learns from the training data by growing ensembles of regression trees on random subsets (bootstrap samples) of the training data (Francke et al., 2008a; Schmidt et al., 2023). In each regression tree, the data are recursively partitioned based on splitting rules, where both RF and QRF randomly select the predictors used for splitting. In contrast to random forests, QRF keeps all observations within a node (whereas RF only keeps the mean), which allows us to construct prediction intervals and to assess uncertainty (Schmidt et al., 2023).

The advantages of QRF include that it can handle multiple input variables; makes no assumptions on distributions; and can deal with interactions, nonlinearity and non-additive behavior. In terms of limitations, it does not allow for easy interpretation of the effects of single predictors, and model predictions will always be within the range of observations, i.e., if the predictors in the period of application exceed the range represented in the training data set (hereafter called out-of-observation-range (OOOR) data points), we can expect overestimations or underestimations of the target variable if the respective predictor has a continuing monotonic effect in this range (Schmidt et al., 2023).

With respect to modeling suspended sediment concentrations, studies have shown that QRF is very well suited to modeling sedigraphs and estimating annual SSY (Francke et al., 2008b, a; Zimmermann et al., 2012) and that it performs favorably compared to sediment rating curves and generalized linear models (Francke et al., 2008a; Schmidt et al., 2023). On a related note, RF (which QRF is based on) outperformed support vector machines and artificial neural networks (Al-Mukhtar, 2019) in modeling suspended sediment concentrations.

In a previous study (Schmidt et al., 2023), we trained QRF models at daily resolution on data of the two gauges Vent and Vernagt using the limited available time series of turbidity (4 and 15 years) and long records of the primary predictors, namely discharge (Q), precipitation (P) and air temperature (T) (Fig. 2). These can be seen as drivers or proxies for processes and catchment conditions crucial to sediment dynamics in high-alpine areas: e.g., discharge determines sediment transfer and erosion within the channel; precipitation is key for runoff formation and hillslope erosion, hillslope–channel coupling and the triggering of mass movement events; and air temperature controls the activation of sediment sources (e.g., sub- and proglacial sediments and their transport by glacier meltwaters or hillslope destabilization by permafrost thaw) and whether precipitation occurs as rain or snow. In addition to these primary predictors, we derived ancillary predictors to describe antecedent conditions and the cumulative effects thereof: e.g., longer-term discharge behavior may exhaust sediment sources or lead to sediment storage, long warm periods may deplete snow cover and accelerate glacier melt associated with increased subglacial sediment transport, and high antecedent moisture conditions may amplify surface runoff or promote mass movements in response to precipitation events. To capture these antecedent conditions while keeping the correlation between the derived predictors as low as possible, we computed sums of the primary predictors in non-overlapping windows of increasing sizes (e.g., 24, 24–72, 72 h to 7 d and 7 to 20 d ahead of each time step) (Schmidt et al., 2023; Zimmermann et al., 2012). The lengths of the considered time periods differ between the models as we optimized them for the highest model performance using the Nash–Sutcliffe efficiency (Schmidt et al., 2023). Additionally, we used the day of the year to capture the seasonality and the rate of change in discharge as predictors (Francke et al., 2008b). The models were validated at both gauges. At gauge Vernagt, we trained validation models on the two available 2-year training data sets (2000–2001 and 2019–2020) and compared the model estimates to measurements from the respective period not included in the training data (NSE of 0.73 and 0.6 for daily SSY; for details, see Schmidt et al., 2023). At gauge Vent, a 5-fold cross-validation demonstrated the superior performance compared to sediment rating curves, especially in periods containing threshold effects due to extreme events. For the past 5 decades, OOOR data points (see Sect. 2.4.2) were rare, which strengthened the notion that the available training data covered the majority of typical situations.

2.2 Study area

The two studied gauges Vent Rofenache (hereafter Vent, operated by the Hydrographic Service of Tyrol) and Vernagt (operated by the Bavarian Academy of Sciences and Humanities) are located in the Rofental in the Ötztal Alps, Austria (Fig. 3). The two corresponding nested catchments of 98.1 and 11.4 km² span elevations ranging from 1891 m a.s.l. at gauge Vent and 2635 m a.s.l. at gauge Vernagt to 3772 m a.s.l. The area is characterized by a relatively warm and dry climate (for this alpine setting), with a mean annual temperature of 2.5 ^∘C at gauge Vent (Strasser et al., 2018) and average annual precipitation as low as 660 mm at gauge Vent but a strong precipitation gradient with elevation (Schmidt et al., 2023). Both catchments are heavily glacierized (28 % and 64 % glacier cover in 2015 (Buckel and Otto, 2018)), but accelerating glacier retreat has been observed since the beginning of the 1980s (Escher-Vetter, 2007; Braun et al., 2007; Abermann et al., 2009). Apart from the glaciers, land cover at high elevations is dominated by bare rock or sparsely vegetated terrain, whereas mountain pastures and coniferous forests occupy lower elevations. Geology is dominated by biotite-plagioclase, biotite and muscovite gneisses, variable mica schists, and gneissic schists (Strasser et al., 2018).

Figure 3 Map of the catchment area above gauge Vent, with nested catchment above gauge Vernagt and major glaciers Vernagtferner (VF) and Hintereisferner (HEF). Meteo stations recording precipitation and temperature are located close to the gauges. Map based on 10 m DEM of Tirol (Land Tirol, 2016), glacier inventory of 2015 (Buckel and Otto, 2018) and river network from tiris open government data (Land Tirol, 2021).

The river Rofenache is a tributary stream of the Ötztaler Ache, one of the largest tributaries to the river Inn. The glacial-to-nival hydrological regime shows a pronounced seasonality, with almost 90 % of discharge occurring during snow and glacier melt from April to September (Schmidt et al., 2023). Mean annual suspended sediment concentrations at gauge Vent were the highest in an Austria-wide comparison (Lalk et al., 2014). Annual suspended sediment yields in Vent averaged 1500 t km⁻² a⁻¹ with an even more pronounced seasonality compared to discharge (99 % of the annual SSY transported from April to September) (Schmidt et al., 2023).

2.3 Input data

2.3.1 Climate projections

We used projections of air temperature and precipitation of the European part of the COordinated Regional Downscaling Experiment (EURO-CORDEX) (Jacob et al., 2014), which have been downscaled and bias-corrected for use in their hydrological model by Hanzer et al. (2018). The EURO-CORDEX initiative provides regional climate model results to enable exploration of the impacts of future climate change at a comparatively high horizontal resolution. This is beneficial for modeling future sediment export, for example since regional climate model simulations provide higher precipitation intensities, which are entirely missing in the global climate model simulations (Jacob et al., 2014), and are thus more likely to capture erosion-relevant changes in precipitation. The data used in this study and by Hanzer et al. (2018) were the result of six different regional climate models (RCMs) driven by five different global climate models (GCMs), resulting in a total of 14 different GCM–RCM modeling chains (Table 1). These are forced by three different emission scenarios expressed as representative concentration pathways (RCPs), which correspond to added radiative forcings of 2.6, 4.5 and 8.5 W m⁻² at the end of the 21st century relative to pre-industrial conditions, i.e., RCP2.6 (intervention scenario assuming peak CO₂ concentrations in the middle of the century, followed by slow decline and negative emissions), RCP4.5 (intermediate scenario with peak emissions mid-century followed by strong decline) and RCP8.5 (assuming no implementation of climate mitigation policies, considerably and steadily increasing emissions and greenhouse gas concentrations over time) (Jacob et al., 2014; Hanzer et al., 2018). This results in a total of 31 RCP–GCM–RCM combinations. The horizontal resolution of the original EURO-CORDEX projections is 0.11^∘ (≈ 12.5 km).

Table 1 Overview of EURO-CORDEX scenario simulations used in this study (unaltered from Hanzer et al., 2018, and distributed under CC BY 3.0, https://creativecommons.org/licenses/by/3.0/, last access: 22 December 2023).

Hanzer et al. (2018) used statistical downscaling to represent the local scale. For this, they bias-corrected all RCM outputs using at least 20 years of observations in the period 1971 to 2005 by using quantile mapping, which matches the distributions of the climate model simulations of the current climate to the distributions of observations. This is necessary, especially in Alpine regions, because the 12.5 km spatial resolution of RCMs (despite being comparatively high) cannot sufficiently resolve topographical and climatological heterogeneities (Hirschberg et al., 2021). Hanzer et al. (2018) concluded that the corrected RCM outputs adequately represent the mean and variability of the observed climate. They interpolated the meteorology to a 100 m grid, and we used the projections for the two grid cells that are located closest to the gauges Vent and Vernagt. EURO-CORDEX simulations are provided at daily resolution, and Hanzer et al. (2018), have disaggregated them to 3 h resolution to capture diurnal variability in the energy fluxes. We re-aggregated these data to a daily resolution to match the temporal resolution of our QRF models (see Sect. 2.3.3).

2.3.2 Hydrological projections

We used discharge projections of the physically based hydroclimatological and snow model AMUNDSEN (Hanzer et al., 2018), which is a fully distributed energy and mass balance model including glacier evolution (Δ h method) and that is particularly adapted to high-mountain catchments on a small to regional scale. It comprises a glacier retreat module and has been extensively validated for historic conditions, especially with respect to snow distribution (Hanzer et al., 2016). This is especially beneficial for modeling sediment dynamics since AMUNDSEN can model processes such as changes in glacier melt that govern discharge dynamics and are crucial to sediment fluxes in these high-alpine areas. AMUNDSEN was forced by the downscaled, bias-corrected and temporally disaggregated EURO-CORDEX simulations of precipitation and air temperature described above (as well as relative humidity, global radiation and wind speed) and modeled snow, glaciers and hydrology in the Ötztal Alps until 2100 (Hanzer et al., 2016). The AMUNDSEN model was calibrated and extensively validated for the period 1997–2013 using water-balance-derived mean areal precipitation, snow depth recordings, Landsat- and MODIS-derived snow extent maps, glacier mass balances, and runoff recordings (Hanzer et al., 2016).

The temporal extent of both the meteorological and the hydrological projections is 2006 to 2100, but since data are not available for the entire year of 2006, we use the period of 2007 to 2100. Additionally, three HadGEM-driven models ended in November 2099. The years 2007 to 2020 overlap with observation data at gauge Vent and the results of the previous study at gauge Vernagt (Fig. 2), which we utilize to verify our model results (see Sects. 2.4.1 and 3.2).

2.3.3 Adjustment of input data for the QRF model

As the QRF models were trained at daily resolution, we aggregated the Q and T projections from 3 h resolution to daily means and the P projections to daily sums. However, comparing the AMUNDSEN Q projections to observations in the overlap period (2007–2020; see Fig. 2) showed underestimation of Q during the glacier melt period at gauge Vernagt and substantial overestimation of Q during the snowmelt period at gauge Vent. Hanzer et al. (2018) have acknowledged that but have left the overestimations (percent bias, PBIAS, of up to 23 %) unaltered since “mainly changes rather than absolute values are analyzed; these partial model biases likely do not affect the main conclusions”. However, in our case, it is necessary to correct the discharge data since SSYs are sensitive to discharge amounts, and additionally, unrealistic discharge amounts exceeding the maximum discharge value in the training data represent a challenge (see Sect. 2.4.2). Also, it is necessary to represent discharge seasonality and thus discharge origins as accurately as possible as, usually, more sediment is exported during glacier melt than at similar discharge levels during snowmelt (Schmidt et al., 2023).

For consistency, we applied the same bias correction as Hanzer et al. (2018), i.e., quantile mapping, using the methodology by Gudmundsson et al. (2012) as implemented in the R package qmap (Gudmundsson, 2016). Due to strong season-dependent biases, Hanzer et al. (2018) have performed quantile mapping for each season individually. We followed this approach, yet in order to best represent discharge seasonality, we shifted the limits of the seasons by 1 month (NDJ, FMA, MJJ and ASO instead of DJF, MAM, JJA and SON) as this corresponded better to seasons with similar characteristics of overestimation or underestimation.

2.4 Analyses of model limitations and uncertainties

To analyze model performance and identify the limits of the applicability of the presented QRF modeling approach, we verified the modeled SSYs against measurement data in the overlap period (2007–2020) (Sect. 2.4.1); assessed the frequency of OOOR data points, including by how much the observation range of the predictors is exceeded in the projections; and analyzed whether the modeled SSY are sensitive to changes in these predictors (Sect. 2.4.2).

2.4.1 Verification of model results based on observed data

To determine how well the obtained SSY projections correspond to SSY derived from turbidity measurements, we compared the model results in the overlap period (2007 to 2020) to measurements at gauge Vent. Lacking continuous direct measurements at gauge Vernagt, we used estimated SSYs for the years 2007–2020 from the QRF model trained on measurements for the years 2000, 2001, 2019 and 2020 (Schmidt et al., 2023) for the comparison (see also Fig. 2). To simplify the descriptions in the results, we hereafter refer to these estimates as observations as well. As the hydroclimatic projections (and, thus, the SSY projections resulting thereof) do not mimic the characteristics of single years (let alone months or days) but merely reproduce their distribution, we compared the distributions of observed and simulated annual SSYs. To test for significant differences between these distributions, we used the two-sample Kolmogorov–Smirnov test, which is able to handle the non-normal distribution of some groups, as implemented in the R package stats version 3.5.1 (R Core Team, 2018). Additionally, we assessed whether the seasonality of sediment export is accurately represented in the model results by comparing mean monthly SSY.

2.4.2 Assessment of limits of applicability

As mentioned in Sect. 2.1, a known limitation of QRF is that model bias can result if the predictors in the projection period exceed the range of observed values used as training data. This limitation is a direct consequence of the numerical characteristics of RF and QRF, which are incapable of extrapolation. In order to assess how often and to what extent the model results are affected, we performed a series of analyses (overview in Fig. 4), described in the following.

Figure 4 Overview of analyses performed with respect to OOOR days. Max(x_p,train) denotes the maximum value in the training data, n refers to the number of OOOR days, and e refers to the exceedance extent. Subscripts are as follows: p for predictor, rcp for emission scenario, and ts for time slice. The analysis of OOOR precipitation days (3) is described in Sect. 2.5.

Analysis of out-of-observation-range days

First, we quantified how often OOOR days occurred for each projection proj and predictor p (i.e., Q, P or T) as the mean annual number of OOOR days per year:

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{proj}}}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{years}}}}\cdot {\sum }_{i=\mathrm{1}}^{n_{\mathrm{years}}}{n}_{\mathrm{p},\mathrm{proj},i},\end{array}$$

where n_years is the number of years, and $n_{\mathrm{p},\mathrm{proj},i}$ is the number of OOOR days in a given predictor and projection in a given year i. Additionally, we determined the exceedance extent (percentage of the maximum value in the training data max(x_p,train)), i.e., by how much the maxima in the observations (i.e., the training period) max(x_p,train) were exceeded on the OOOR days j (i.e., in the projection period):

$$\begin{array}{}\text{(2)}& e_{\mathrm{p},j}={\displaystyle \frac{x_{\mathrm{p},j}-\max \left(x_{\mathrm{p},\mathrm{train}}\right)}{\max \left(x_{\mathrm{p},\mathrm{train}}\right)}}\cdot \mathrm{100}\forall x_{\mathrm{p},j}>\max \left(x_{\mathrm{p},\mathrm{train}}\right),\end{array}$$

where x_p,j is the value of the predictor p on OOOR day j.

Although the same limitations of QRF theoretically also apply to predictors falling below the training minima, we did not consider them since Q and P already contained very low or zero values and cannot fall below zero. Likewise, for T, minimum temperatures are already well below zero, and a further decrease (if it occurs) is not further physically relevant to sediment transport. Similarly, we only considered summer precipitation at both gauges (i.e., May–September) to exclude snowfall events that are not directly relevant to sediment dynamics.

Sensitivity analysis

To assess the potential effects of the abovementioned OOOR days on the model results, we performed a sensitivity analysis for the three primary predictors Q, P and T (Fig. 4). For this, we determined the 90th percentiles of the number of OOOR days $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{proj}}}$ and the exceedance extents e_p,j of all projections, i.e., n_p,90 (in d a⁻¹) and e_p,90 (in % of max(x_p,train)), for each predictor p. These values were considered to represent a severe case for possible model deficits due to their lacking extrapolation capability. We created a respective test data set for the sensitivity analysis from the training data in the 14-year overlap period (2007–2020): we selected the corresponding number of days n_ov with the highest values of the respective predictor in the overlap period, as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{n}_{\mathrm{ov}}& =\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{90}}}\cdot t_{\mathrm{ov}}\left(\mathrm{rounded}\mathrm{to}\mathrm{integers}\right),\mathrm{with}{t}_{\mathrm{ov}}=\mathrm{14}\mathrm{a}\\ & \mathrm{as}\mathrm{the}\mathrm{length}\mathrm{of}\mathrm{the}\mathrm{overlap}\mathrm{period},\end{array}\end{array}$$

and altered them by adding or subtracting the respective e_p,90. For example, n_p,90 of Q in Vent is 0.55 d a⁻¹; therefore, we changed n_ov= 0.55 ⋅ 14 ≈ 8 d by the e_p,90 of 9.6 m³ s⁻¹ (Table 2).

Table 2 Amount of reduction or increase in the sensitivity models (average exceedance extent $\stackrel{\mathrm{‾}}{e_{\mathrm{p}}}$ in units of the corresponding predictor) and number of days with a reduction or increase on average per year ($\stackrel{\mathrm{‾}}{n_{\mathrm{p}}}$) and in total in the 14-year period (n_ov).

We used the resulting altered time series of the primary predictors to compute the corresponding ancillary predictors (that describe antecedent conditions; see Sect. 2.1), ran the QRF model with them and compared mean annual SSC after the alterations to the original data set. Thus, we performed six individual runs for the sensitivity assessment at each gauge, two (one where the predictor was increased and one where it was reduced) for each of the three primary predictors Q, P and T. We chose to compare mean annual SSC instead of annual SSY as discharge is needed to compute SSY so that the alterations in Q would have affected the estimated SSY twice.

Classification

Third, we assessed whether the sensitivity analysis was informative for the different RCPs, time slices and predictors, i.e., if the sensitivity analyses contained sufficiently extreme conditions to represent the projections. For this, we determined the mean exceedance extent per predictor p, emission scenario rcp and time slice ts $\stackrel{\mathrm{‾}}{e_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}$ and the mean number of OOOR days per year $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}$. We compared these to the e_p,90 and n_p,90 and marked the respective predictor–time-slice–RCP combination with N or S in Table 5 if ≥ $\mathrm{1}/\mathrm{3}$ of projections $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>n_{\mathrm{p},\mathrm{90}}$ or $\stackrel{\mathrm{‾}}{e_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>e_{\mathrm{p},\mathrm{90}}$ and with bold N or S if this applied to ≥ $\mathrm{2}/\mathrm{3}$ of the projections (see Fig. 5 and Table 5).

Figure 5 Example of the classification based on the OOOR analysis. The boxplots show the distribution of exceedance extents $e_{\mathrm{p},\mathrm{rcp},j}$ per RCP and predictor on all days j within the respective time slice ts. Gray numbers denote the average $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}$ of all projections within the respective RCP, time slice and predictor. Colored numbers indicate the number of projections with the mean exceedance extent $\stackrel{\mathrm{‾}}{e_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>e_{\mathrm{p},\mathrm{90}}$, i.e., projections more extreme than the sensitivity analysis. Black numbers indicate the number or projections with the mean number of OOOR days per year $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>n_{\mathrm{p},\mathrm{90}}$, i.e., projections with more OOOR days than the sensitivity analysis.

2.5 Analysis of model results

We analyzed the model results, i.e., estimated annual yields in the projection period for trends, as well as shifts in seasonality. To assess trends, we used two methods implemented in the R package FUME (Santander Meteorology Group, 2012): the Mann–Kendall test, which is a non-parametric tool to detect linear trends (specifically, we used a version that was modified to detect trends in serially correlated time series; Madsen et al., 2014; Yue et al., 2012) and Sen's slope estimator (Sen, 1968) to assess trend magnitude. Further, we compared the estimated yields in three time slices (near future: 2011–2040, intermediate future: 2041–2070, and far future: 2071–2100; see also Fig. 2), comparable e.g., to Jacob et al. (2014) and Hanzer et al. (2018). To assess changes or shifts in the seasonality of sediment export, we compared the mean monthly yields of the observations and the projections.

To assess whether the detected trends are sensitive to the potential underestimation of yields on OOOR precipitation days, we multiplied the daily yields estimated by our QRF model on these days by a factor of 5, i.e., assuming a very severe underestimation in the original estimates (Fig. 4). We chose this factor as it is close to the most severe exceedance extents in precipitation at both gauges, which are 456 % at gauge Vernagt and 442 % at gauge Vent (see also Sect. 3.3). We then compared the trends in annual SSY of the altered time series to the trends in the original QRF estimates. All analyses were conducted with the statistical software R (R Core Team, 2018).

3 Results

3.1 Verification of bias-corrected discharge for the present climate (2007–2020)

The bias-corrected discharge data yield more adequate representations of measured monthly discharge amounts and their seasonal distribution (Fig. 6), as well as mean annual discharge volumes (Table 3). At gauge Vernagt, for example, maximum mean monthly Q in the observations and the bias-corrected data is in August, whereas the original AMUNDSEN simulations suggested a maximum in July. Nevertheless, some underestimation of August discharge remains at gauge Vernagt. At gauge Vent, the original AMUNDSEN simulations substantially overestimated discharge amounts in April to July, i.e., the snowmelt period, which was successfully corrected by the bias correction.

Table 3 Comparison of mean annual discharge volumes based on the original AMUNDSEN output, observations and bias-corrected AMUNDSEN estimates in the overlap period (2007–2020).

Figure 6 Comparison of mean monthly discharge [mm] at gauge Vernagt (a) and Vent (b) derived from measurements, unaltered multi-model means of the original AMUNDSEN output (solid lines) and multi-model means of the bias-corrected AMUNDSEN output (dashed lines) in the overlap period (2007–2020).

3.2 Verification of modeled SSY for the present climate (2007–2020)

We find good agreement between observed and modeled annual SSY at both gauges (Fig. 7), and the Kolmogorov–Smirnov test does not yield significant differences between the observations and model results in mean annual sediment yields. Nevertheless, years with extreme annual yields (both lower and higher) occur in the model results, especially under RCP4.5 and RCP8.5 (e.g., for Vent, max. 3250 t a⁻¹ in RCP4.5 vs. 2120 t a⁻¹ in the observations), likely due to the higher sample size in the projections (42 or 196 years in the projections compared to 14 years in the observations; see also description of Fig. 7).

Figure 7 Comparison of annual specific SSY in the overlap period (2007–2020) derived from measurements (obs, n=14 years) and QRF modeling results per RCP (n=42 for RCP2.6 and n=196 for RCP4.5 and RCP8.5, respectively) at gauges Vernagt (a) and Vent (b).

Similarly, the seasonality of sediment export is well represented in the overlap period of the projections (Fig. 8), and the Kolmogorov–Smirnov test does not yield significant differences in comparison to the seasonal distribution of the measurements. Monthly SSYs tend to be slightly lower in the projections in August at gauge Vernagt and in July and August at gauge Vent. Similar patterns had already become apparent in the comparison of mean monthly discharges at gauge Vernagt (Fig. 6).

Figure 8 Comparison observations (see also Fig. 2) to QRF model forced by climate projections (multi-model means per emission scenario) during the overlap period (2007–2020).

3.3 Assessment of limits of applicability

3.3.1 Out-of-observation-range days

Generally, we find more frequent OOOR days and higher exceedance extents in later time slices and in the higher-emission scenarios (Fig. A1 in the Appendix). At both gauges, OOOR days in Q are relatively rare, but in the higher-emission scenarios and later time slices, OOOR days become more frequent, and individual exceedance extents of more than 100 % occur (Fig. A1 and Table 4). Exceedances in temperature are more frequent at gauge Vernagt, especially under RCP8.5 and after 2040. At gauge Vent, there are only a few OOOR days in T, except under RCP8.5 after 2070.

Table 4 Mean and maximum exceedance extents e_p of the three primary predictors, namely discharge (Q), precipitation (P) and air temperature (T), across all emission scenarios and time slices in percent of the maximum during the training period max(x_p,train) and original units.

OOOR data points in summer precipitation are rather rare at both gauges. However, precipitation shows very high exceedance extents of up to ca. 450 % (RCP4.5 after 2070 and RCP8.5 before 2040 at gauge Vent, RCP4.5 after 2070 at gauge Vernagt; Fig. A1 and Table 4). This corresponds to daily precipitation sums of approx. 280 and 240 mm d⁻¹ at gauge Vent and Vernagt, respectively, and is equivalent to over a third of the current mean annual precipitation at gauge Vent (687 mm (Hydrographic Yearbook of Austria, 2016)). Yet, even without the most extreme cases, exceedance extents in precipitation can be quite severe, which corresponds to very heavy precipitation events.

3.3.2 Sensitivity analysis

Figure 9 shows the results of the sensitivity analysis, which indicates the extent to which annual SSY estimates may be affected by underestimations on days with OOOR observations for the different predictors. P is the most sensitive predictor at gauge Vernagt. Yet, although precipitation amounts were altered quite substantially on the 6 d with the highest precipitation in the overlap period (by almost 29 mm; see Table 2), the effect on mean annual SSC is small (≤ 3 %). The same applies to Q, the most sensitive predictor at gauge Vent, which was altered by almost 10 m⁻³ s⁻¹ on 8 d. Temperature is the second most sensitive parameter at gauge Vernagt, while at gauge Vent, the temperature alterations had little effect on mean annual SSC. At gauge Vent, P is the second most sensitive parameter but with a maximum effect of < 2 % on mean annual SSC.

Figure 9 Results of the sensitivity analysis for gauges Vernagt (a) and Vent (b) with respect to mean annual SSC of the years with altered days. Subscripts signify an increase (inc) and decrease (dec) in the predictors by the respective average exceedance extent e_p,90 and frequency n_p,90 as identified based on the exceedance analysis. Gray numbers represent the number of altered days in the test data set (2007–2020).

The results of the sensitivity analysis also give an indication of the behavior of the QRF model in response to OOOR data points: as expected, we generally observe a decrease in mean annual SSC if we decrease the predictor values on the selected days (Q, P and T) and vice versa. However, for most predictors, the decrease is more pronounced than the increase (although the same days were altered by the same extent). We presume that this is due to the described incapability of QRF to extrapolate. Thus, we can expect to underestimate the additional effect, e.g., of precipitation exceeding max(P_train).

3.3.3 Classification

Table 5 shows that, until 2070, all predictors and RCPs fall within the conditions covered by the sensitivity analysis (with the exception of Q after 2040 at gauge Vent under RCP8.5). This implies that the results of the sensitivity analysis are informative in these cases and that we expect similar or smaller effects of OOOR days on mean annual SSC or SSY in the projection period. After 2070, exceptions occur at both gauges in two out of three RCPs and for several predictors, which implies that the uncertainty is higher than in the results of the sensitivity analysis.

Table 5 Results of the classification per emission scenario, predictor and time slice. S denotes more than one-third (italic) and two-thirds (bold, does not occur) of the projections with $\stackrel{\mathrm{‾}}{e_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>e_{\mathrm{p},\mathrm{90}}$ (assumed in sensitivity); N denotes more than one-third (italic) and two-thirds (bold, does not occur) of the projections with $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}>n_{\mathrm{p},\mathrm{90}}$.

The OOOR analysis showed that very high exceedance extents occur in precipitation and that precipitation is a sensitive parameter at both gauges (although the effect on mean annual SSC was small). Additionally, we find that heavy summer precipitation becomes more intense (and only slightly more frequent) (Fig. A2 in Appendix): the 99.5th​​​​​​​ percentile of the summer precipitation projections increases over time, which suggests that precipitation events become more intense. At the same time, the number of precipitation events that exceed the 99.5th percentile determined from the precipitation observations in the overlap period (2007–2020) hardly increases on average, which suggest that precipitation events of a certain strength do not become (much) more frequent. We also find an increase in daily SSY associated with heavy precipitation events (Fig. A2 in Appendix).

Thus, we additionally assessed whether the trends in annual yields were sensitive to changes in yields on days with OOOR precipitation (Sects. 2.5 and 3.4.1) as extreme precipitation can be very important for sediment dynamics (e.g., by triggering mass movements).

3.4 Projections of future sediment export: changes in annual yields, timing of peak sediment and changes in seasonality

3.4.1 Changes in annual yields and timing of peak sediment

The resulting projections suggest an overall decrease in mean annual SSY for both gauges and each of the three emission scenarios, which is more pronounced at gauge Vernagt (Fig. 10 and Table 6). Accordingly, we consistently find significant negative trends in the projections (2007–2100) in mean annual SSY (Table 6). The differences between the RCPs are small and are smaller than the spread within individual RCPs (Fig. 10). Accordingly, trends in mean annual SSY are only slightly more negative in the high-emission scenarios. With respect to the 99th percentile of annual SSY estimates, trends are less strong than for mean SSY estimates at gauge Vent, while at gauge Vernagt, the trends in the 99th percentile are even stronger than for mean annual SSY estimates.

Table 6 Trends in mean and 99th​​​​​​​ percentile of annual specific SSY projections (2007–2100) given as Sen's slope [t km⁻² a⁻²] for the original estimates and the altered estimates (5-fold-increased SSY on days with OOOR precipitation). Significance levels: ^* indicates 0.05, ${}^{**}$ indicates 0.01, and ${}^{***}$ indicates 0.001.

Figure 10 (a, c) Mean annual suspended sediment yields per RCP, with the entire range from minima to maxima of the individual projections indicated by the colored envelopes. The period after 2070 is faded due to lower confidence in model results, as pointed out in Sect. 4.3. (b, d) Annual SSY of all respective years and projections within the three time slices.

Negative trends were detected for all individual projections as well: at gauge Vent, 26 (out of 31) are significant (α= 0.05, Sen's slope ranging from −10.8 to −3.8 t km⁻² a⁻²), and at gauge Vernagt, 30 of 31 are significant (α= 0.05, Sen's slope ranging from −15.2 to −6.1 t km⁻² a⁻²).

The trend in the altered time series (with 5-fold-increased daily yields on days with OOOR precipitation; see Sect. 2.5) hardly differs from the trend in the original time series (Table 6). Specifically, at gauge Vernagt, trend characteristics are basically unchanged. The only trend reversal occurs in the 99th percentile at gauge Vent under RCP8.5, where the trend is slightly positive (and significant) instead of negative. We conclude that the overall trend characteristics remain very robust, even if we assume severe underestimation of the model on days with OOOR values in the predictors. Thus, the overall future sediment budget seems to be governed by their mean behavior rather than solitary extreme events.

The synopsis with estimates of annual SSY for the past decades shows that we find increases in annual yields at both gauges up until sometime between 2000 and 2020 and decreases afterwards, which are much more distinct at gauge Vernagt (Figs. 10 and 11). This suggests that peak sediment has already been reached or is underway at both gauges and occurs simultaneously with peak water (Fig. 11).

Figure 11 Estimated timings of peak water and peak sediment. Black lines indicate past mean annual Q from measurements and mean annual SSY estimates of QRF model; colored lines correspond to different RCPs (compare to Fig. 10). Underlying data have been smoothed using a 15-year moving average.

3.4.2 Changes in seasonality

Mean monthly SSY is projected to decrease substantially during the glacier melt period in August in all RCPs and at both gauges (Fig. 12). As a result, the highest mean monthly SSY shifts from August to July or even to June under RCP8.5 after 2070 at both gauges. Additionally, the spring onset of sediment export is projected to occur at a higher rate and slightly earlier in the year under RCP8.5 after 2070. This represents a decrease in the importance of glacier melt for sediment export. After 2070, only relatively minor further changes are projected under RCP2.6 and RCP4.5, whereas RCP8.5 experiences further decreases in mean monthly SSY throughout the year.

Figure 12 Seasonality of mean monthly SSY in three time slices and emission scenarios (as the monthly mean of all respective modeling chains).

At gauge Vent, a slight increase in mean July SSY is projected after 2070 under RCP2.6. This is likely related to an increase in discharge since this increase is not visible in mean monthly concentrations (Fig. A3 in the Appendix). It also has to be considered that only three projections are averaged for RCP2.6 (as compared to 14 in the other RCPs), which makes it less robust with respect to outliers.

A comparison to the seasonal distributions determined from the altered time series (5-fold-increased SSY on days with OOOR precipitation) showed only very slight differences, which indicates that the seasonal distribution is also insensitive to underestimations of SSY on days with heavy precipitation.

4 Discussion

Testing new methods to estimate future suspended sediment export from glacierized high-alpine areas can provide important information, e.g., with regard to flood hazard, sediment or water quality management, since estimating such changes has so far been limited to relatively rough approximations. This study represents the first attempt, to our knowledge, to derive SSY projections using a machine-learning approach and to investigate them in synopsis with reconstructed past SSYs.

4.1 Projected changes in sediment export and location of peak sediment

The presented SSY projections in the Ötzal, Austria, suggest an overall decrease in annual SSY (Fig. 10). This is consistent across emission scenarios (i.e., with respect to the ensemble means), as well as all individual projections (i.e., based on the 31 different RCP–GCM–RCM chains). This is linked to shrinking of glaciers and thus decreases in glacier melt as temperature is projected to increase in all models (between +1.1 (RCP2.6) to +3.8 ^∘C (RCP8.5) in the annual mean) (Hanzer et al., 2018). Accordingly, snow cover is projected to decline, especially in low elevations. With respect to precipitation, there is no clear general trend with respect to annual sums; both decreases and increases are projected by individual models (−14 % to +24 %), and multi-model averages are close to zero for all scenarios, yet with a general shift in precipitation from summer to winter

The modeled decrease in SSY is much more distinct at gauge Vernagt, where a trend analysis in the previous study (Schmidt et al., 2023) showed significant positive trends in the period until 2020. At gauge Vent, significant positive trends were detected if all data points since the 1970s were considered (Schmidt et al., 2023). However, if only the years after the distinct increase around 1980 were considered, the trend was slightly negative (Schmidt et al., 2023). This suggests that peak sediment has already been reached or is underway at both gauges and occurs simultaneously with peak water.

These findings match the expectations of conceptual models that sediment yield from deglaciating basins will initially increase (due to increases in glacial erosion, sediment supply accessibility, transport capacity and occurrences of extreme floods) and subsequently decrease as glacier masses decline, meltwater volumes and freeze–thaw weathering decrease, and vegetation colonizes (Antoniazza and Lane, 2021; Zhang et al., 2022). It is expected that peak sediment may lag behind peak meltwater, with a lag that can be up to decades or centuries (Delaney and Adhikari, 2020). This lag is hypothesized to be scale dependent – i.e., it will be shorter for areas closer to the glacierized regions – and to depend on the changes in erosive precipitation: if erosive precipitation decreases, peak sediment occurs simultaneously with peak water, while increasing or stable erosive rainfall scenarios are associated with a lag (Zhang et al., 2022). Indeed, for the study area of this study, a decrease in summer precipitation sums (i.e., June to August, which is the time of minimum snow cover and thus maximum erodibility) is projected (Hanzer et al., 2018). At the same time, heavy precipitation events are projected to become more intense (and only slightly more frequent; see Fig. A2 in the Appendix). However, the negative trend in discharge appears to prevail as our estimates suggest that peak sediment coincides with peak meltwater (Fig. 11).

Sediment export projections differed only slightly (if at all) between emission scenarios; i.e., the spread between projections within one emission scenario is much larger than differences between ensemble means of the three RCPs (Fig. 10). This is due to irreducible internal climate variability, which has been acknowledged in climate modeling (Deser et al., 2012), as well as in erosion modeling (Hirschberg et al., 2021). It should be noted that comparisons to RCP2.6 need to be treated with care as it comprises fewer GCM–RCM combinations (only 3 as compared to 14 in the higher-emission RCPs). Nevertheless, the absence of major differences between RCPs is in accordance with findings by Gobiet and Kotlarski (2020) that “until the middle of the 21st century […] it is projected that climate change in the Alpine area will only slightly depend on the specific emission scenario”. Accordantly, Hanzer et al. (2018) projected glacier volumes to decline by 60 %–65 % until 2050, “largely independent of the emission scenario”.

4.2 Changes in seasonality and response to heavy precipitation events

Despite the overall decrease in SSYs, our results suggest that high annual SSYs are possible, especially at gauge Vent and towards the end of the century (Fig. 10). Additionally, yields on days with heavy precipitation may increase at both gauges – in absolute terms and in relation to the annual export (Fig. A2 in Appendix). This is reasonable, given that increases in heavy precipitation intensity and/or frequency in the European Alps have been detected in measurement data from the past (e.g., Hiebl and Frei, 2018; Scherrer et al., 2016; Gobiet and Kotlarski, 2020), as well as future projections (Gobiet and Kotlarski, 2020; Jacob et al., 2014; Kotlarski et al., 2023) – despite the overall decrease in summer precipitation mentioned above. As a result, we can expect an increase in sediment-related harmful events triggered by heavy precipitation, such as flash floods and gravitational mass movements (i.e., debris flows, landslides) (Huggel et al., 2012; Savi et al., 2020; Gobiet and Kotlarski, 2020). Similar expectations, i.e., increasing high-magnitude transport events in the context of an overall decrease, have e.g., been expressed with respect to bedload in South Tyrol (Coviello et al., 2022). Such a development would have important implications e.g., for sediment management in hydropower production and flood hazard management as the flood risk can increase if cross-sections are reduced after sedimentation, and potential backwater effects need to be considered (Nones, 2019) as much of the damage is associated with transported solids rather than the water itself (Badoux et al., 2014).

Our findings suggest a shift in sediment export seasonality since the highest mean monthly SSY shifts from August to July (or even June) due to substantial reductions in sediment export in July, August and September at both gauges (Fig. 12). This is linked to the projected distinct reductions in glacier melt (Hanzer et al., 2018) and appears to be reasonable given that glacier melt has so far been the dominant transport medium of suspended sediments at these gauges (Schmidt et al., 2023). These results are not sensitive to potential underestimations of SSY on days with very heavy precipitation. Such shifts in seasonality and the concomitant overall reduction in fluvial sediment transport will likely have severe effects on biodiversity, i.e., flora and fauna of glacier-fed streams (Milner et al., 2009, 2017; Gabbud and Lane, 2016; Huss et al., 2017).

4.3 Limitations

As a potential limitation to the presented quantile regression forest approach, out-of-observation-range data points in the predictors can lead to underestimates in SSY on the affected days. Yet, the analysis of such incidents in synopsis with the results of the sensitivity analysis showed that, before 2070, the effect on annual yield estimates is ≤ 3 % Fig. 9). This is very small given the overall high variability in SSC (Vercruysse et al., 2017; Delaney et al., 2018; Schmidt et al., 2023). On a similar note, even assuming rather generous increases in yields on days with OOOR precipitation altered the trends only marginally, which shows that underestimations on individual days with OOOR precipitation have little effect on long-term annual averages. However, we have less confidence in the model results after 2070 for two reasons. First, more frequent and severe OOOR incidents occur during this time, especially in the high-emission scenarios, and fewer projections are covered by the assumptions of the sensitivity analysis. We can therefore expect a higher uncertainty in the model results. Specifically, the effect of underestimation for single large events will aggravate. Second, more than a few glacier simulations suggest that glaciers could have disappeared almost entirely by 2070 (Hanzer et al., 2018), which implies a major shift in the hydro-geomorphic functioning of these catchments. While our QRF models were able to model threshold effects better than sediment rating curves (Schmidt et al., 2023) (likely because they are not bound to linear or monotonous relationships), this is only true for effects that are represented in the training data. Thus, the results for the period after 2070 need to be treated with caution. We have indicated this in the presentation of our results by using transparency or dashed lines instead of solid lines.

As a more general limitation, there are several other factors with the potential to substantially alter and influence sediment dynamics in the study area, which we cannot consider in our models. This concerns geomorphological changes, such as increased paraglacial erosion: debuttressed slopes may trigger landslides and rockfalls, and indeed, increased debris flow and rockfall activity have been shown in response to warming in other areas, likely associated with intensified alpine permafrost thaw (Savi et al., 2020; Hartmeyer et al., 2020; Huggel et al., 2012). Additionally, sediment availability and accessibility increase as erodible landscapes expand (Li et al., 2021a), and subglacial sediment availability might also increase (more subglacial sediment can be accessed by meltwaters as the equilibrium line altitude retreats upslope) until the glacier size becomes smaller than a critical size (Delaney and Adhikari, 2020; Zhang et al., 2022). Although these processes are likely already partially reflected in the observations used for the model training, their intensity may still be too low to be sufficiently learned by the model. Thus, future intensification of these processes could lead to higher sediment export rates than our estimates suggest and might thereby affect the estimated location of peak sediment. Notwithstanding, there are also several factors that could lead to decreases in sediment export, such as decreases in connectivity (such as the formation of supra-, sub- or proglacial lakes or outwash fans which act as sediment traps) or decreasing glacial erosion as glaciers recede (Zhang et al., 2022). Additionally, freeze–thaw weathering may decrease (Hirschberg et al., 2021), and it is not clear how quickly the deglaciating landscapes will stabilize, e.g., through eluviation of fine materials and fluvial sorting of sediment, which progressively increases the resistance to entrainment and transport (Ballantyne, 2002; Lane et al., 2017), or vegetation colonization (Haselberger et al., 2021; Altmann et al., 2023; Musso et al., 2020; Eichel et al., 2018). Many of these processes are ultimately governed by temperature and/or precipitation, and we have chosen the predictors to act as proxies (e.g., antecedent moisture and temperature conditions could be crucial for mass movements). While this is out of the scope of the presented study, we encourage future studies to work towards including more advanced proxies for geomorphological changes.

4.4 Uncertainties

The presented results are associated with uncertainties, which are a combination of uncertainties inherited from the underlying climatological and hydrological projections and uncertainties inherent in the QRF approach. Climate model uncertainty represents a combination of uncertainties in assumptions of future anthropogenic greenhouse gas emission, GCM uncertainty (different GCMs produce different responses to the same radiative forcing) and RCM uncertainty (different RCMs forced by same the GCM produce different regional responses) (Evin et al., 2021; Gobiet et al., 2014). It has been found that EURO-CORDEX simulations may be biased towards “too cold, too wet, too windy”, but these uncertainties are mostly within the observational uncertainties, and simulations “reproduce fairly well the recent past climate despite some biases” (Vautard et al., 2021). To address this, it was recommended to carry out bias correction, which has been performed by means of quantile mapping for the precipitation and temperature projections (Hanzer et al., 2018). The hydrological model results are also associated with uncertainties, such as the tendency to overestimate spring runoff, winter snow accumulation and glacier mass balances. We have addressed this through bias correcting the discharge projections, which resulted in a more adequate representation of discharge seasonality and volumes. Certainly, bias correction methods such as quantile mapping in turn introduce uncertainties, e.g., by assuming that the biases are stationary, i.e., do not change over time (Gudmundsson et al., 2012). Hydrological simulations that do not show the necessity for this correction could eliminate this issue. Uncertainties in the QRF approach have been addressed in a previous study and include the tendency to underestimate rare, high-magnitude daily SSYs (albeit with small effects on the respective annual yields); the underestimations on days with OOOR values (which had small effects until 2070, as discussed in detail above); and the choice of temporal resolution (i.e., daily compared to hourly resolution involves some loss of information, e.g., on precipitation intensities, but the effect was also found to be small) (Schmidt et al., 2023). Since QRF is a data-driven approach, the quality of the estimates hinges on the underlying training data set, as well as the choice of predictors, i.e., a large and varied enough data set in combination with predictors that meaningfully represent the most important processes improves the quality of the estimates (Schmidt et al., 2023). Thus, future studies are recommended to explicitly sample extreme events and/or verify the representativity of the training data set.

5 Conclusions

We found decreasing trends in annual SSY at both gauges regardless of the emission scenario, which suggests that peak sediment was already reached between 2000 and 2020. This is linked to substantial decreases in discharge volumes, especially during the glacier melt phase in late summer, as glaciers continue to melt as a result of increasing temperatures. These findings persist even if yields on days with projected heavy precipitation are inflated by a factor of 5. Despite the projected overall decrease, high(er) annual yields are possible, likely in response to heavy summer precipitation, which could become more intense according to the projections. This discrepancy has important implications for sediment management but also e.g., for flood management.

To our knowledge, this study represents the first attempt to combine machine learning for suspended sediment modeling with climate and hydrological projections in order to derive projections of sediment export in high-alpine areas. It demonstrated that a quantile regression forest can be a valuable tool for this application. We addressed known issues of the QRF, i.e., underestimations on days where predictors in the projection period exceed the range represented in the training data. The influence of such underestimations on the results were shown to be negligible until 2070. We conclude that the presented results are much more uncertain after 2070, partly because of more frequent and severe out-of-observation-range data points, but mainly since a major shift in the functioning of the hydro-geomorphic system can be expected as deglaciation is quasi completed.

However, while the chosen predictors represent proxies for crucial processes controlling sediment transport in these high-alpine environments, several potentially crucial geomorphological factors that could increase or decrease sediment export (and thereby change the projected trends and location of peak sediment) could not be taken into account. These include increases in rockfalls and landsliding, changes in connectivity, or vegetation colonization. Future studies are encouraged to include these factors more explicitly.

Appendix A

Figure A1 Results of the classification based on the sensitivity and OOOR analyses. The boxplots show the distribution of exceedance extents for each RCP and predictor on all days j within the respective time slice ts. Gray numbers denote the average $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{rcp},\mathrm{ts}}}$ of all projections within the respective RCP, time slice and predictor. Colored numbers indicate the number of projections with mean exceedance extent $\stackrel{\mathrm{‾}}{e_{\mathrm{p},\mathrm{ts}}}>e_{\mathrm{p},\mathrm{90}}$ (as used for the sensitivity analysis). Black numbers indicate the number or projections with the mean number of OOOR days per year $\stackrel{\mathrm{‾}}{n_{\mathrm{p},\mathrm{ts}}}>n_{\mathrm{p},\mathrm{90}}$ (as used for the sensitivity analysis).

Figure A2 Analysis of summer precipitation projections (a–d) and SSY projections (e–h) at gauges Vernagt and Vent. (a, c) Frequency of days with heavy summer precipitation (> 99.5th percentile of the training data, i.e., 33.7 mm d⁻¹, Vernagt, and 28.88 mm d⁻¹, Vent). (b, d) Intensity of heavy summer precipitation events over time, expressed as the 99.5th percentile. (e, g) Sediment export on days with precipitation > 99.5th percentile of the training data. (f, h) Sediment export on days with precipitation > 99.5th percentile of the training data relative to the respective annual yields.

Figure A3 Seasonality of mean monthly SSC in three time slices and emission scenarios.