Controls on Sediment Transport From a Glacierized Catchment in the Swiss Alps Established Through Inverse Modeling of Geomorphic Processes

As climate warms, hydrology and geomorphology in glacierized catchments are evolving, changing sediment export from these catchments, thus impacting downstream ecosystems and communities. Currently, much uncertainty persists regarding interactions among geomorphic processes that evacuate sediment from glacierized catchments. Here, we present a catchment‐scale numerical model of subglacial and proglacial sediment transport with debris meltout processes. We apply the model in a Monte Carlo framework to suspended sediment data collected over 2014–2020 from the Fieschertal catchment in the Swiss Alps, assessing possible combinations of geomorphic processes responsible for suspended sediment discharge. The ensemble of model outputs quantifies the interaction of different geomorphic processes, including the trade‐off between bedrock erosion and sediment previously stored in the catchment. Model runs suggest that, at some periods, up to 20% of the sediment leaving the glacier is deposited in the proglacial area relative to the catchment's total sediment discharge. This shows that the model captures the proglacial area change from sediment source to sediment sink in different hydrological and glaciological conditions. Furthermore, the findings highlight the impact of glacier retreat on the catchment's sediment dynamics, which both reduces the proglacial area's slope and introduces additional sediment to the proglacial area through debris meltout. The model outputs and the parameters quantify the interaction among geomorphic processes, yielding key insights into the drivers of sediment transport in alpine regions. The implications of these interactions are discussed in the context of interpreting processes responsible for controlling erosion rates from glacierized regions and modeling sediment transport in glacierized catchments.


Introduction
Climate warming in alpine regions causes glaciers to alter their dynamics (e.g., Dehecq et al., 2019), increase their melt (e.g., Huss & Hock, 2018), and retreat substantially (e.g., Zekollari et al., 2019).These processes will alter the ways that glacierized catchments expel sediment (e.g., Antoniazza & Lane, 2021).Changes in glacierized landscapes also impact downstream ecosystems, by modulating temperature, water, and sediment fluxes in glacierized catchments (Milner et al., 2017;Stibal et al., 2020).Specifically, changing sediment discharge, caused by the evolution of these hydrological, geomorphic, and glaciological processes impacts anthropogenic and ecological systems downstream (Milner et al., 2017).Furthermore, sediment from glaciers can fill in alpine reservoirs (Anselmetti et al., 2007;Ehrbar et al., 2018;Steffen et al., 2022) and abrasion from sediment during high flows can damage hydropower infrastructure (Felix et al., 2016b;Müller-Hagmann et al., 2020).Increased sediment flux caused by increased geomorphic activity in High Mountain Asia, for example, impacts hydropower operations in this region (Li et al., 2022).Glacier retreat also impacts the biodiversity in many catchments, where sediment from glaciers can both bury plants as well as provide a new surface for plants to grow on (e.g., Bosson et al., 2023;Cauvy-Fraunié & Dangles, 2019).Yet, evaluating the processes controlling sediment discharge and the response to climate warming over decadal timescales remains difficult.
From a processes understanding perspective the role of different sediment transport and erosion processes in glacierized catchments remains unknown (Figure 1).Over long timescales, or on glaciers with minimal sediment stored beneath them, abrasion and quarrying provide the means to erode bedrock and produce sediment (e.g., Iverson, 2012;Ugelvig et al., 2018).These erosive processes may be a primary control on sediment discharge in catchments with high gradients and little sediment storage (Herman et al., 2015) and on sediment production (e.g., Brinkerhoff et al., 2017).Additionally, hillslope processes, such as frost cracking, produce sediment that is deposited on glaciers and transported supra-and en-glacially as debris (Andersen et al., 2015;Rowan et al., 2022;Scherler & Egholm, 2020).
On shorter timescales, however, fluvial sediment transport largely influences sediment discharge from glaciers, which is impacted by glacier hydrology and sediment availability (e.g., Beaud et al., 2018;Delaney et al., 2019;Riihimaki et al., 2005).Over yearly to decadal timescales, englacial and supraglacial debris, deposited on or in the glacier, respectively, is advected with the ice.This material can meltout on the lower reaches of glaciers and provide sediment to proglacial areas (Benn et al., 2012;Scherler & Egholm, 2020;van Woerkom et al., 2019).Periodic mass wasting events, such as stream-bank collapse, rockfall, or landslides, may add sediment to proglacial areas, where it can be mobilized by fluvial activity in the proglacial area if the system is supply limited (Battista et al., 2022;Mancini & Lane, 2020;Theler et al., 2010;Warburton, 1990).Proglacial areas may also expose new sediment sources or create depositional environments for sediment leaving the glacier as time passes (Bakker et al., 2018;Delaney et al., 2018;Lane et al., 2017;Mancini & Lane, 2020).Furthermore, sediment transport conditions change with the transition from pressurized flow beneath the glacier (e.g., Beaud et al., 2018) to open channel flow in the proglacial area where sediment discharge is measured (e.g., Church & Ryder, 1972).These processes act, and interact, together to dictate the sediment discharge from a catchment by creating and exporting sediment (e.g., Schmidt et al., 2022).In turn, a multitude of processes may be represented in a sediment discharge record downstream of a glacier, as sediment deposition and mobilization processes occur in the transition zone before sediment arrives at a sink or measurement station (Castelltort & Van Den Driessche, 2003).
Yet, reconciling these geomorphic processes that occur across daily to century timescales with available records of sediment export from glacierized remains difficult.Sediment transport or concentration increases are currently being observed from glacierized regions in the Alps, Greenland, High Mountain Asia, and the Andes (e.g., Bendixen et al., 2017;Costa et al., 2018;Li et al., 2021Li et al., , 2022;;Vergara et al., 2022).In these cases, the measurements are generally collected far from glaciers, near the outlet of the catchments.Closer to glaciers, however, available records of sediment discharge from gauging stations have generally been collected over only a single season or two (e.g., Swift et al., 2005Swift et al., , 2021;;Willis et al., 1996).These short-term data sets have been extremely valuable in understanding some of the key glacial processes responsible for subglacial sediment transport (e.g., Riihimaki et al., 2005;Swift et al., 2005).Yet, a comparison of sediment discharge records across glaciers and seasons shows substantial interannual variability.Furthermore, individual glaciers respond differently to variations in water discharge within the same season (Delaney et al., 2018).In turn, the brevity of the data collection periods and variability in the records pose challenges when establishing the geomorphic processes and effects of climate with available observations (Schmidt et al., 2022(Schmidt et al., , 2023)).
Interpreting the records of sediment transport collected at hydrological stations in glacierized catchments thus requires understanding the role of these different geomorphic processes.Given the multitude of processes mentioned above, modeling this system presents an equifinality problem, where several combinations of processes could result in similar amounts of sediment discharge from the catchment (Figure 1; e.g., Beven & Binley, 2014).Yet, establishing the role of these processes remains important, as each component of the system will evolve differently as glaciers retreat (e.g., Zekollari et al., 2019), hydrology changes (e.g., Brunner et al., 2019), and proglacial processes evolve (e.g., Antoniazza & Lane, 2021).
To effectively describe the interaction of the multitude of processes that serve as a control on sediment discharge in glacierized catchments, we present a lumped element model of subglacial erosion and sediment transport, debris meltout, and proglacial sediment transport in front of the glacier-but not mass wasting processes.We then describe a continuous record of suspended sediment transport collected at the water intake of a hydropower facility at Fieschergletscher in the Swiss Alps from 2014 to 2020.Over this period, the glacier's retreat has caused the proglacial area to grow substantially.The model is applied to sediment transport data from Fieschergletscher to explore the roles and interactions of the sediment transport processes.This is accomplished using a Monte Carlo method that compares the observational suspended sediment discharge data to model outputs with randomly selected parameter values representing sediment storage in the catchment, bedrock erosion, and sediment transportability.Using these data, model outputs, and resultant parameters, we examine the interactions amongst the processes controlling suspended sediment discharge from this glacierized catchment.We discuss the implications of the model's outputs to develop a broader understanding of sediment transport processes in glacierized catchments.Lastly, we evaluate the model's limitations, including its lumped spatial discretization and the omission of mass wasting processes.

Study Site
Sediment and water discharge data from the Fieschergletscher catchment were collected at a hydropower water intake at Wysswasser, in Canton Valais, Switzerland (46.5°N, 8.12°E) from 2014 to 2020 (Figure 2).Fieschergletscher flows southeast from the top of the Aaremassif, comprised of the Central Aare-Granite.The catchment area is 58.6 km 2 .Water from glacier melt has been collected for hydropower purposes since 1975 with a water intake located at 1,650 m a.s.l.The hydrological and sediment concentration measurements described here were collected at this location.The highest point of the catchment is Finsteraarhorn at 4,274 m a.s.l.As of 2020, the catchment is 60% glacierized with the largest ice mass being the 13 km-long
The glacier has retreated roughly 1.63 km since 1891 (e.g., Linsbauer et al., 2021).In 1975, the glacier terminus was located at the water intake used for hydropower.By 2021, glacier retreat had created a proglacial area spanning 1 km, from the water intake to the glacier terminus (Figure 2).
In this study, we leverage external data sets.Average ice thickness across the glacier was estimated from Grab et al. (2021) to be 225 m.We evaluated the proglacial area's length by measuring the distance from the hydropower intake to the first emergence of a proglacial stream from the glacier from SwissTopo imagery.This length was evaluated from imagery collected in November 2011, fall 2014, and September 2017.We assumed a linear retreat rate among these images (Figure 3).

Methodology
Here we describe the processing of a sediment discharge record from the Fieschergletscher catchment.Then, we present a numerical model of sediment transport and discharge for glacierized catchments and the model's required inputs of catchment and glacier topography, along with water discharge.Lastly, we explain the inversion scheme used to infer the model parameters by comparing the modeled and observed sediment discharge records.

Sediment Discharge Record
From January 2014 until December 2020, nearly continuous records of turbidity and water discharge were collected from the water intake at Fieschertal (Figure 4).The purpose of these measurements was to understand the sediments' impact on hydro-abrasion of Pelton turbines that the hydropower company operates (Abgottspon et al., 2022).A complete description of the installation scheme can be found in Felix et al. (2016bFelix et al. ( , 2021)).
The turbidity data were collected until the 2016 season using an Endress + Hauser Turbimax CUS41 turbidity meter and, from 2016 onwards, using an Endress + Hauser Turbimax CUS5D turbidity meter.A calibration relationship between turbidity and suspended sediment concentration was established through the collection of water samples (Figure S1 in Supporting Information S1).Water samples were placed in an oven to evaporate the water.The dry residual suspended sediment and dissolved load were weighed to yield the sediment concentration of each sample.The data were saved at a 1 min interval with both instruments.The data presented here (Figure S2 in Supporting Information S1) are aggregated over 1 hr intervals, using hourly mean values.
The sediment discharge by mass at a given time, Q m (kg s 1 ; Table 2), is given by where Q w is water discharge (m 3 s 1 ) and SSC is suspended sediment concentration (kg m 3 ).Volumetric sediment discharge is given by, Q = Q m ρ 1 b , where ρ b is the density of bedrock, assumed to be 2,650 kg m 3 (Table 2).Continuous water discharge Q w was made available by the hydropower company in the same location and time interval as the sediment measurements over the study period.
Measured effective suspended sediment export rate from the catchment εm is determined as i=t Q i is the total amount of volumetric sediment discharged from the catchment from the period t to t n , and A is the catchment area.

Forward Model of Sediment Transport, Erosion, and Glacier Change
Here, we describe the forward model for (a) subglacial erosion and sediment transport, (b) englacial and supraglacial debris meltout, and (c) fluvial sediment transport (Figure 4).Each of the three processes exists as an individual element with the glacier components transferring sediment to the proglacial element as a flux and the proglacial element expelling sediment from the domain.The model also considers glacier retreat so that the elements' size changes throughout the model run.The model's lumped nature is designed to be applied to observations, for instance, the data set described in Section 3.1, with the parameter optimization scheme, such as the one presented in Section 3.3.

Glacier Erosion and Sediment Transport Element
This element aims to model the storage and transport of subglacial sediment to the proglacial area and quantify the subglacial erosion that occurs.Following Delaney et al. (2019), the element simulates the evolution of the thickness of a subglacial till layer H g , which we define as fluvial transportable sediment that originated from glacier erosion of bedrock-but in a single element representing the whole glacier that retreats and changes in area (Figure 4).The density of the till layer is assumed to be the same as bedrock density (ρ b ) so that no density conversion occurs as the bedrock is transformed into sediment.
Erosive processes such as abrasion and quarrying add material to the subglacial till layer.Conversely, fluvial sediment transport in supply-and transport-limited regimes mobilizes or deposits sediment, thus removing or adding material to the till layer, respectively (Brinkerhoff et al., 2017;Delaney et al., 2019).To quantify these processes, we modify the Exner Equation (Exner, 1920a(Exner, , 1920b;;Paola & Voller, 2005), a mass conservation relationship, so that it evolves the till layer height given the changing size of the glacier, along with the erosive and sediment transport conditions where H g is the till layer height, l g (w g ) is the length (width) of the glacier, ṁ is a sediment production rate, representing the addition of sediment to the till layer through erosion, and the retreat rate is given by r r .The last term, Hg , is a factor that accounts for the conservation of subglacial and proglacial sediment as the glacier changes in size as it advances or retreats (see Section 3.2.2).
The sediment flux Q g represents fluvial sediment transport below the glacier and is calculated using the mobilization equation modified from Delaney et al. (2019) (4a) (4b) where Qg is the sediment transport capacity from the glacier, that is, the maximum amount of sediment that could be transported given the water's flow velocity and channel width, as discussed below.σ is a sigmoidal function of that enables a smooth transition from supply-limited to transport-limited sediment transport conditions over the till height range: transition from transport-to supply-limited regime during periods with high sediment transport capacity, relative to the till layer height.
Condition 4a in Equation 4 represents the case where bedrock erosion exceeds sediment mobilization and, thus sediment transport exists in a transport-limited regime.Condition 4b allows sediment mobilization to smoothly transition between transport-and supply-limited regimes.Here, when H g is small, sediment mobilization is limited to the sediment production term ṁ.
We calculate sediment transport capacity from the glacier Qg using the total sediment transport relationship by Engelund and Hansen (1967) where f r is the Darcy-Weisbach friction factor (Table 2), ρ b (ρ w ) is the density of the bedrock or sediment (water; Table 2), D 50 is the median sediment grain size, ŵg denotes the width of the channel floor that integrates the sediment transport rate across the width of the subglacial channel, and τ represents the shear stress between the water and the channel bed (Table 1).We determine the shear stress τ through the Darcy-Weisbach formulation where v g = Q w S is the water velocity and S is cross sectional area of the channel.
To evaluate the size of the subglacial channel S and the width of the channel floor ŵg , we leverage the subglacial hydraulics model presented in Delaney et al. (2019) to which we refer for further details.This model is based upon the assumption that the cross-sectional area S of the subglacial channel can be calculated using characteristic water discharge and hydraulic gradient of the channel over a sufficiently long period.In this way, the model reconciles the time-scale differences between variations in water discharge and the size of the subglacial channel, controlled by creep closure and melt opening of the subglacial channel (Röthlisberger, 1972).The model then uses the channel shape to evaluate a channel floor width ŵg with respect to the channel's area S.
The source term ṁ is defined as, where ε is the sediment production rate of the glacier, and H max is a till height beyond which no further bedrock erosion may occur (Table 1).

Debris Meltout and Glacier Retreat Element
This element of the model represents the meltout of sediment from englacial sources, as well as supraglacial debris, which originates through processes such as headwall erosion (e.g., Andersen et al., 2015;Krautblatter et al., 2013).Both Equations 3 and 12 contain terms that conserve sediment and account for the meltout of sediment as the glacier retreats.To capture these processes in Equation 3, we define the term Hg so that the first condition results in no change in the till layer height as the glacier retreats at rate r r .H p is the sediment height in the proglacial area (see Section 3.2.3),l p (w p ) is the proglacial area length (width), and l g (w g ) is the glacier's length (width).
The second term captures the absorption of the proglacial sediment of height layer H p into the subglacial till layer as the glacier advances.This term also represents the conservation of mass in the subglacial environment through the redistribution or thinning of subglacial sediment as the subglacial area grows.The glacier and proglacial area widths are represented by w p and w g , respectively.
Similarly, the term Hp in Equation 12captures processes related to glacier advance and retreat in the proglacial area so that the first condition conserves the mass of sediment in the proglacial area and represents the emergence of subglacial and effective thickness of englacial sediment (H g and H e , respectively) and deposition in the proglacial area as the glacier retreats.
We define the emergent sediment flux Q e , or the flux of emergence of englacial and subglacial sediment, through glacier retreat and debris advection as where H e is the effective thickness of debris, that is the integrated quantity of debris in the ice column, including basal ice, englacial ice, and surficial debris on the ice surface.v t represents the near-terminus velocity of ice so that debris can leave the glacier when it is not retreating, effectively creating a moraine (e.g., Rowan et al., 2022).Note, that the near-terminus velocity of this glacier during the study period is close to zero and thus we set v t = 0 (Millan et al., 2022).

Fluvial Sediment Transport Element
To describe fluvial transport of sediment, we implement an approach similar to the subglacial sediment transport element (Section 3.2.1).To account for fluvial sediment transport in the proglacial area, we use Equation 3, adjusted for the fluvial environment H p is the sediment layer height in the proglacial area, Q p is the sediment flux in the proglacial area (Table 1), and Q e is the flux of englacial sediment (see Section 3.2.2).Note that the source term has been removed from Equation 3, as we assume that the rivers in proglacial areas create or detach negligible amounts of bedrock through processes such as hydro-abrasion.The equation's last term represents the transfer of sediment to the proglacial area by mass conservation of sediment and meltout of englacial sediment and subglacial sediment Hp ) from the glacier as it retreats at a rate of r r (see Section 3.2.2).As above, the density of this sediment layer H p is assumed to be that of bedrock (ρ b ).
To describe the evolution of the divergence of the flux in Equation 12, we implement a similar scheme as in Equation 4 that represents sediment transport Q p in both supply-and transport-limited regimes.The sediment flux Q p is calculated using a modification of the mobilization equation from Delaney et al. ( 2019) Qp is the sediment transport capacity in the proglacial area, σ(H p ) is determined with Equation 5 and controls the transition from supply-to transport-limited regimes.
We use the sediment transport relationship in Equation 6 to establish Qp .To evaluate shear-stress driving this relationship is implemented from the fluvial hydraulics model following Tucker and Slingerland (1997).This formulation assumes that the channels are wide compared to their depth, and there is uniform flow.The model evaluates shear-stress τ p , to evaluate sediment transport capacity, at the river bed as Water Resources Research f p is a friction factor and Ψ p is the slope of the proglacial area.We describe channel width ŵp empirically in a power law, where k c is a scaling factor (Leopold & Maddock, 1953).Equipped with τ p , we replace channel width ŵg in the subglacial area with channel width ŵp in the proglacial area in the sediment transport relationship (Equation 6) to solve for sediment transport capacity in the proglacial environment Qp ( ).

Discretization, Initial Conditions, and Model Inputs
Equations 3 and 12 are lumped in space so that the two evolving sediment layer heights H g and H p are each represented by a single element.The length of each element may change through the model run.Sediment fluxes from the subglacial and englacial elements of the model are passed to the proglacial element (Figure 4).
To evolve sediment layer heights H g and H p in time as they respond to erosion and sediment transport conditions, along with changing glacier extent, in Equations 3 and 12, the model uses the VCABM differential equation solver (Variable step, variable-Coefficient Adams-Bashforth-Moulton method; Hairer et al., 1992;Rackauckas & Nie, 2017;Radhakrishnan & Hindmarsh, 1993).We impose a maximum time step of 6 hr, ensuring that diurnal variations in melt input are captured.In practice, the solver typically uses a time step of roughly 2.5 hr, which varies depending on sediment transport conditions and solver tolerance.Longer time steps, for instance, occur over periods when glacier melt ceases along with sediment transport (i.e., over winter months).The numerical parameters used in the solver are presented in Table S1 in Supporting Information S1.
Initial conditions in till height must be selected for the subglacial and proglacial areas (H g0 and H p0 ).Previous work has shown the sensitivity of model outputs to the initial subglacial till height condition H g0 (Delaney & Adhikari, 2020;Delaney & Anderson, 2022).However, we note that the lumped nature of the model means that only one value of each parameter must be selected and, thus, a range of initial conditions can be examined (Section 3.3).
The model requires inputs of glacier surface and bed topography, proglacial topography, and water discharge at the catchment's outlet.We assume all water discharge at the catchment outflow comes from the glacier and moves through the proglacial area.Thus, no additional water is added to the proglacial area from pluvial sources or lakes.
In other applications, however, the model could include water discharge inputs separately into the different elements of the model, representing hydrological processes such as rainfall at lower elevations or drainage from proglacial lakes that occur solely in the proglacial environment.The glacial surface and proglacial topography evolve across a model run in response to observed glacier retreat.We assume that the glacial area is a rectangle (Figure 4), with glacier length l g that evolves in response to glacier retreat: Here l g at time t is a result of the glacier length at the beginning of the model run t 0 minus the cumulative amount of glacier retreat at time t given the retreat rate r r that varies in time with the observed retreat rate (Figure 3).
The proglacial area size must also be evaluated in response to glacier retreat.Its length l p is given as where l is the catchment length, which is constant throughout the model run (Figure 4).
Glacier thickness h, as well as bed and surface slope, must also be given due to the need to evaluate the hydraulic gradient in the subglacial sediment model (Section 3.2.1;Delaney et al., 2019).Ice thickness at a time t is given by Water Resources Research where t r is the thinning rate at time t, assumed to be 2 m a 1 , based upon observed glacier-wide mass balance estimates from several Swiss glaciers (Bauder et al., 2022).The slope of the glacier bed is given by max z s (t) is the maximum elevation of the glacier, min z s (t) is the minimum elevation of the glacier, h(t) is the glacier thickness at time t.
Note that as glacier length l shrinks, the glacier slope generally steepens (Huss et al., 2010;Zekollari & Huybrechts, 2015).The gradient of the glacier surface is given by The gradient of the proglacial area Ψ p is evaluated as where z bp is the base level of the proglacial area, assumed to be fixed and representing the water intake in this case (Figure 2).z b (t) is the minimum elevation of the glacier bed at time t, allowing the gradient to evolve as the glacier retreats.z b (t) and z bp are established from topographic analysis of the proglacial area at different times, and they are used to establish the proglacial gradient in Figure 3.We add the sediment height in the proglacial area H p so that sediment deposition (mobilization) can steepen (reduce) the proglacial area's gradient.The sediment height H p is added to the observed surface elevation (z b (t) and z bp ).Experience with the model in this application shows that the effect of adding H p is very small, as its value is on the order of millimeters (see Sections 4.2 and 4.3).Yet, its implementation means that in other applications an equilibrium characteristic of rivers between sediment transport and supply can be reached, if deposition occurs, therefore steepening the channel and thus increasing sediment transport capacity (Equation 14; e.g., Wickert & Schildgen, 2019).
The glacier width w g is 2,200 m, and the width of the proglacial area w p is 200 m, based upon examining their respective widths from aerial imagery.

Parameter Optimization
To understand the role of the model parameters and their sensitivity in the system, we test model outputs with different randomly selected combinations of parameter values against suspended sediment discharge data from Fieschergletscher using a Monte Carlo framework (Section 3.1).Running the model in this way allows us to examine the interaction amongst parameters, as well as to identify the parameter space with the most favorable outcomes.The independent sampling, where a set of parameter values does not depend on the outcome of previous parameter sets, also allows us to examine the equifinality of the system, whereby different combinations of parameters, and thus contributions of geomorphic or hydrological processes, may yield equally acceptable results (Beven, 1996;Beven & Binley, 2014).
To implement this scheme, we create a vector of randomly selected parameter values θ to run in the forward model Median sediment grain size is D 50 , representing a control on the transportability of sediment in both the subglacial and proglacial environments as applied to Equation 6.The glacier's sediment production rate is ε (Equation 8) and is constrained by the minimum erosion rate presented in Hallet et al. (1996) and the maximum erosion rate observed from the nearby Aletschgletscher (Table 3; Delaney et al., 2018;Meile et al., 2014).H g0 and H f 0 are the initial till and sediment heights below the glacier representing subglacial and proglacial sediment existing before the model initialization, respectively (Section 3.2.4).They are estimated to range from 5 × 10 4 (smaller than the Δσ parameter used in Equation 5) to 5 cm.We estimate 5 cm to be the maximum thickness of sediment, averaged over the glacier bed.This value is chosen as it results in adequate availability of sediment over the model run, while also maintaining the characteristic reduction in sediment production with sediment layer thickness in Equation 8. H e is the total thickness of englacial and supraglacial debris (Equation 11) representing the quantity of sediment melting out of the glacier.Its upper bound is set by reasonable values presented in McCarthy et al. (2022).Values and sampling distributions are given in Table 3. Log-uniform sampling distributions are used when a potential value of a parameter spans several orders of magnitude.
Model outputs are tested against data using the absolute misfit (ξ) of model outputs (Q p ) and observations (Q), given by We sum model outputs and data over 5 days time periods to reduce the impact of lags between the transport of water and sediment through the catchment, as the model does not contain a mechanism to consider hysteresis in sediment transport following variations in water velocity through the catchment (Williams, 1989).We tune the friction factors, source quantile, and response times (Table 2; f r , s p , and Δt; Section 3.2.1;Delaney et al., 2019, and f p in Equation 14) such that values of subglacial and proglacial water velocities in the observed range of water discharges align with measured values (subglacial: ∼1.2 m s 1 , proglacial: ∼0.6 m s 1 ; e.g., Magnusson et al., 2012;Werder et al., 2010).The proglacial channel shape factor k c is selected so that the resulting channel widths are roughly 80 m at maximum water discharge, similar to observed channel widths at Fieschertal from aerial photos.We also make no density conversion between the bedrock and the till and sediment layers, such that they both have the density of bedrock (ρ b ).The near-terminus ice velocity near v t is close to zero (c.f.Millan et al., 2022), so the flux of sediment leaving the glacier is solely through meltout in Equation 11.
The forward model is run 3 × 10 6 times with randomly selected parameter combinations of θ (Table 3).We present the model results below, from parameter combinations resulting in the top 0.02% (600 retained runs) of model performance given by misfit ξ (Equation 23).

Measured Sediment Discharge Trends
From spring 2014 to winter 2020, we measured 4.98 × 10 8 kg of suspended sediment leaving the catchment.
During the same period 1.01 × 10 9 m 3 of water left the catchment.The mean sediment concentration in the water was 0.49 kg m 3 , and the total effective sediment export over the study period was 3.24 mm, equal to a mean annual effective sediment export rate of 0.42 mm a 1 assuming a bedrock density of 2,650 kg m 3 .
The largest quantity of suspended sediment discharge occurred in 2017, when 37,100 m 3 of sediment left the catchment.The greatest mean suspended sediment concentrations over the season occurred in that year (Table 4).Total water discharge was highest in 2018 (Table 4 and Figure 5a).Across the 7 years, there is a positive correlation between the annual amount of water discharged and the corresponding quantity of sediment discharge (Pearson r: r = 0.78, p-value: p = 0.04, number of samples: n = 7).No significant correlation persists between the  yearly average suspended sediment concentration and sediment discharge (r = 0.6, p = 0.15, n = 7).This lack of correlation suggests that overall sediment discharge is mainly controlled by sediment transport capacity from the water discharge (i.e., transport-limited regime), and not by sediment access (i.e., not supply-limited regime), represented in part by the sediment concentration.The smallest water discharge and sediment discharge volume occur in the first years of the study period, increase toward 2018, and decrease thereafter (Figure 5a).
The correlation between sediment discharge and water discharge over 5-day periods ranges from 0.40 to 0.83 across the individual years.Across the 7 years, Pearson's r correlation between water discharge and suspended sediment concentration is lower compared to the relationship between water discharge and sediment discharge (Figure 5a).The exception is the year 2016 when sediment discharge correlates more strongly with suspended sediment concentration (r = 0.45) compared to water discharge (r = 0.40).This is possibly due to reduced access of meltwater to sediment over this year.The highest correlation between water discharge and sediment discharge occurs in 2018 (r = 0.83) and occurs with the highest annual water discharge (Figure 3, Table 4).The strong correlation in 2018 may have resulted from an abundance of sediment caused by increased glacier retreat and sediment discharge the year before, that is, in 2017.It could be that the high melt in 2018 could have transported this sediment in increased quantities.
Concave-up (down) relationships between cumulative amounts of water and sediment discharge result from sediment being expelled from the catchment at a greater (lower) rate than water (Figure 5b).The convexconcave nature of these relationships varies substantially across the years (Figure 5).In 2015 the relationship was concave-up or linear 34% of the time, suggesting that decreasing sediment availability occurred through the melt season.In 2016, the relationship between cumulative water and sediment discharge was more linear or concave-up with sediment transport increasing at a faster rate relative to water discharge 53% of the time, more than any other year.This trend suggests that sediment availability was greater through the melt season in 2016.

Model Behavior
Model outputs from the 600 retained model runs capture annual and seasonal variations in sediment discharge over the study period (Figure 6).The model can capture many of the peaks and cessations in sediment discharge over the 5-day periods (Figure 6f).Over the study period, glacier retreat caused between 2,400 and 9,500 m 3 of emergent sediment to be transferred to the proglacial area (Q e ; Figure 6c), depending on the thickness of the subglacial till (H g ) and englacial sediment layers (H e ).This quantity represents a maximum of 5% of the observed suspended sediment discharged from the catchment throughout the study period (Table 4).The emergent sediment flux is greatest during the period of greater glacier retreat, that is, prior to fall 2017 (Figure 6c).The remaining sediment discharged into the proglacial area resulted from the fluvial transport of subglacial sediment (Q g ).
The retained model runs capture most interannual variations in sediment discharge, with 542 model runs having yearly sediment discharge sums with Spearman Rank correlation of 0.96, when compared to sediment discharge data.Spearman Rank correlation of 0.96, in this case, means that 1 year of the 7 years was out of sequence when the model's annual sediment yields were ranked against the observed ones.The remaining 58 model runs have a Spearman Rank correlation of 0.83.
Throughout the study period, many model runs exhibit periods of sediment deposition in the proglacial area that result from reduced sediment transport conditions in the proglacial environment (Figure 7).The deposition occurs as the water transitions from pressurized flow below the glacier to open channel flow in the proglacial area (Section 3), since the water may lose sediment transport capacity at the glacier margin.Over some periods through the study, 20% of the sediment leaving the glacier is deposited in the proglacial area and not discharged from the proglacial area.Slightly more deposition, and less mobilization, occurred in the proglacial from 2018 onwards, compared to the beginning (Figure 7b).This change is mostly due to the reduced slope, and thus reduced sediment transport capacity, in the proglacial area as the glacier retreated (Figure 3b).
The proportion of sediment evacuated from the catchment that was mobilized in the proglacial area also decreases throughout the model run, from up to 70% in some peak melt periods before 2018 to up to about 30% later in the study (gray area in Figure 7a).This occurs because of a steeper proglacial slope and increased glacier retreat prior to 2018 since both increase sediment transport capacity and thus introduce more sediment to the proglacial area (Figure 6c).Furthermore, by the later years of the model run, sediment mobilization in the proglacial area is reduced by the smaller sediment layer thickness (Figures 7a and 7c).
Decreasing sediment layer height in the ensemble of model outputs shows reduced sediment availability in both the proglacial and subglacial areas over the model run (Figure 7).Till heights (H g ) decrease every year, while in the proglacial area, the sediment height (H p ) decreases throughout each season to close to 1 mm for the last years of the study period.The observed concave-down relationships between sediment and water discharge also suggest seasonal sediment exhaustion throughout the melt seasons (Figure 5b).

Model Parameters
Posterior distributions of all parameters θ are substantially different from their prior sampling distributions and generally lay within their parameter bounds (Figure 8; Table 3).(Log-) Gaussian posterior distributions emerge from (Log-) uniform posterior of initial till height (H g0 ) and sediment grain size (D 50 ), demonstrating their strong influence on sediment discharge from the catchment (Figure 8).The sediment production rate ε ( ) develops a Gaussian posterior distribution as well.However, the distribution is bound on one side by the minimum sediment production rate of 0.06 mm a 1 .The mean sediment production rate of the posterior distribution (0.4 mm a 1 ) aligns very closely with the measured sediment export rate for the catchment (ɛ m = 0.42 mm a 1 ).However, because the glacier covers roughly 60% of the catchment, this sediment production parameter causes sediment to not be replenished in the subglacial till layer, resulting in subglacial sediment exhaustion (Figure 6c).
The initial sediment height condition in the proglacial area (H p0 ) develops a wider distribution than the initial till height below the glacier (H g0 ; Figure 8).This could be due to the smaller proportion of the catchment which is occupied by the proglacial area.While a log-normal posterior distribution does not occur from the log-uniform sampling distribution in the effective englacial sediment height (H e ; Table 3), a log-uniform distribution develops with values less than 6 mm of sediment (Figure 8).This suggests that greater values of H e result in an excess flux of sediment into the proglacial area to meet the criterion of an adequate model run.Furthermore, the potential for small values of H e shows that viable model runs are possible with nearly no debris input into the proglacial area, despite the visual evidence of debris cover at the glacier (Figure 2).
Negative relationships also emerge amongst the three sediment layer height conditions (H g0 , H p0 , H e ).These occur as retained model runs require a range of sediment to be available in the catchment for an adequate parameter combination to occur.Should any one of these sediment sources increase, then the remaining sources would need to decrease to accommodate this condition.
The strongest relationship between variables occurs between the sediment production rate ε ( ) and the initial till height condition (H g0 ), with a negative correlation of nearly 1 (Figure 8).This results in a trade-off between bedrock erosion and sediment availability when the model quantifies sediment transport.All other parameters that correlate positively (negatively) with the sediment production rate ε, correlate negatively (positively) with the initial till height condition (H g0 ).For instance, sediment grain size (D 50 ) correlates negatively with sediment production rate ( ε; Pearson's r: r = 0.21).In this case, smaller sediment grain sizes D 50 , or more transportable sediment, are required for increasing sediment production rates ε as the sediment is being replenished or continually added to the till-layer throughout the model run.Conversely, sediment grain size (D 50 ) correlates positively with initial till height condition (H g0 ; Pearson's r: r = 0.15), indicating that sediment transport capacity must be reduced with greater sediment availability for the model to yield a favorable outcome.This correlation suggests that for greater initial till height conditions (H g0 ), larger sediment grain sizes (D 50 ) are needed to maintain the availability of sediment throughout the model run.
As with the trade-offs between sediment grain size (D 50 ) and initial till height condition (H g0 ), a positive relationship persists between effective englacial debris thickness (H e ) and sediment grain size (D 50 ; Figure 8).This suggests that as greater debris thickness (H e ) increases sediment availability, larger sediment grain sizes (D 50 ) are required to reduce transport capacity and modulate the amount of sediment mobilized in the catchment.No such relationship exists between the initial proglacial sediment height (H p0 ) and grain size (D 50 ; Figure 8).lack of correlation between the two parameters (H p0 and D 50 ) likely occurs due to the large amount of deposition and erosion that occurs in the proglacial area over the model run (Figure 7b).Large variability in the initial proglacial sediment layer height parameter (H p0 ) minimizes the impact of the initial condition on the transportability of sediment leaving the catchment.The variability in proglacial sediment layer height (H p ) over the model run lies contrary to the emergence sediment (Q e ), which only varies with the retreat rate of the glacier (Figure 6c).

Erosion Rates and Sediment Transfer in Glacierized Catchments
Model results show that a multitude of processes interact together to both produce and transfer sediment through the catchment.Of particular interest is the presence of both erosional and depositional processes occurs in the proglacial area, as sediment leaves the glacier and is observed at the water intake (Figures 1, 6, and 7).Erosional and depositional processes in the proglacial area occur as water leaves the pressurized subglacial environment and enters the open channel environment (Church & Ryder, 1972;Mancini et al., 2023;Perolo et al., 2018).In many model runs, substantial amounts of erosion and deposition can occur in the proglacial area over a single season.As a result, the model shows that the changes to sediment dynamics measured downstream of glaciers (Costa et al., 2018;Lane et al., 2019) could be influenced by depositional processes in proglacial areas and glacier-fed river systems (e.g., Blöthe & Korup, 2013;Mancini et al., 2023;Tofelde et al., 2021).Yet, fluxes from the glacier drive the broader sediment transport in the catchment, especially as the model suggests that the proglacial area is in a supply-limited regime by the end of the model run (Figure 7).
The deposition of sediment in the proglacial area may actually increase in the future due to reduced water flow and sediment transport capacity from about 2055 onwards (Figure S3 in Supporting Information S1; Compagno et al., 2021), along with potentially increased sediment flux from the glacier (e.g., Delaney & Adhikari, 2020).The potential abundance of sediment in the proglacial area could make large amounts of sediment available for transport during extreme events, such as meteoric floods, intense rain events, or outburst floods (e.g., Cook et al., 2018;Korup & Tweed, 2007;Zhang et al., 2022).Such extreme events could become more frequent in the future (e.g., Brunner et al., 2019;CH2018, 2018 and may well result in excess sediment transport capacity that mobilizes sediment in the proglacial area (Figure 7).Especially high sediment transport could occur when the water discharge intensity and variability exceed the water discharge that has historically resulted from snow and glacier melt (Lane & Nienow, 2019;Schmidt et al., 2022).These events could evacuate large amounts of sediment at high concentrations (e.g., Cook et al., 2018), negatively impacting the hydropower operations in the catchment (Felix et al., 2016a;Patro et al., 2022).
The modeled flux of emergent sediment (Q e ) that has melted out of the glacier is relatively small, at about 5% of the catchment's yield.Yet, it impacted the catchment's geomorphic processes (Figures 6 and 8), by increasing sediment availability in the proglacial area before 2018 (Figure 7).The effective englacial sediment thickness parameter (H e ) represents the impact of erosional and climatic conditions that produced this sediment from rockfall and hillslope processes over past decades and centuries (e.g., Andersen et al., 2015;Pohl et al., 2015;Rowan et al., 2022;Scherler, 2014).Thus, the impact of englacial debris represents the legacy of past erosional conditions on present fluvial sediment transport, as the debris has been advected through the ice over time.

Omission of Mass Wasting Processes
Despite being visually evident in Fieschertal (Figure 2), we omit mass wasting processes that introduce sediment to the proglacial area from the valley sides or stream channel banks (e.g., Legg et al., 2014).We have chosen not to include these processes in the absence of data quantifying specific events, their highly intermittent nature, and difficulties in accurately modeling them (e.g., Hirschberg et al., 2022).Furthermore, many mass wasting events are above the glacier margin, such as those visible in Figure 2, and are disconnected from the proglacial area, depositing sediment directly onto the glacier.Thus, these processes might be partially captured in the effective englacial sediment height parameter (H e ).Mass wasting events might have a greater impact on the catchment's sediment discharge from 2018 onwards when potential sediment exhaustion results in less sediment being mobilized from the proglacial area (Figure 7).Reduced sediment availability over this period could allow for additional sediment transport from mass wasting processes when there is an excess of sediment transport capacity in the catchment.
In the future, fluxes from mass wasting processes could increase as the glacier retreats, increasing sediment connectivity from the valley sides with glacier debuttressing of steep slopes (e.g., Carrivick & Heckmann, 2017).In 2020, near the end of the study period, a gullied area became visually evident on the east side of the proglacial area (Figure 2b).Sediment deposited from this area covers around 0.04 km 2 of the proglacial area.Assuming sedimentary deposition rates from gullied areas between 21 and 300 mm a 1 (Mancini & Lane, 2020) and a bulk density of sediment of 1,500 kg m 3 , sediment fluxes into the proglacial area could be between 480 and 6,800 m 3 a 1 .This resultant sediment input would have occurred over the last year of the study period in 2020.The quantity is comparable to the emergent sediment flux introduced to the proglacial area by the emergent sediment flux over the study period (Q e ; Figures 6c and 6f).Yet, it remains uncertain if this depositional zone is connected to the proglacial river, where the sediment can be mobilized and evacuated from the catchment (e.g., Mancini & Lane, 2020).

Model Limitations and Modeling Alpine Sediment Dynamics
The model's character allows for the temporal evolution of geomorphic processes to be considered in describing the suspended sediment evacuation from the catchment.However, spatially distributed sediment transport processes are not represented here because of the model's lumped nature, despite their potential importance (Delaney et al., 2023).Spatially distributed subglacial sediment transport processes have been attributed to higher reaches of glaciers experiencing substantial glacier melt, so that meltwater may access increasing quantities of previously stored sediment to evacuate (Delaney & Adhikari, 2020;Li et al., 2021;Vergara et al., 2022).Despite the potential increase in sediment discharge from glaciers (Delaney & Adhikari, 2020;Vergara et al., 2022), the periods of deposition and erosion in the proglacial area suggest that despite the large proportion of glacier coverage, both glacial and proglacial processes should be considered in evaluating the controls of sediment supply to glacially-fed river systems (e.g., Costa et al., 2018;Schmidt et al., 2022;Steffen et al., 2022).The model's simple spatial discretization likely reduces its representation of sediment transport processes over extended periods.However, a more complex framework might be limited in its ability to the interacting geomorphic processes, examine the sensitivity of the multitude of different parameters, or be applied to observations as done here (Beven, 1996).The increased computational time and the multitude of parameters mean that a spatially distributed model constrains its capacity to be run iteratively.Amongst the parameters most difficult to select in a spatially distributed model are the initial sediment layer height conditions, the selection of which can impact outputs throughout a model run (Delaney & Adhikari, 2020).This effect was circumvented here by selecting individual initial till and sediment layer height conditions in the inversion (H g0 and H p0 ; Table 3).
Note that the assumed seasonally steady glacier retreat causes sediment deposition in the proglacial area even during winter months, which increases sediment availability early in the following melt season (Figures 6  and 7).We have chosen to keep the retreat rate constant throughout the year in the absence of further data regarding its seasonal variations in length, however.Additionally, the lumped nature of the model means that the width of the glacier is 2,200 m, when in fact the glacier terminus is roughly 200 m wide, similar to the proglacial area (Figure 2).This means that the effective englacial sediment layer height (H e ) is underestimated in the inversion of the model to arrive at an optimum flux (Figure 8).When the lumped glacier width is scaled by the actual width of the glacier terminus, then the ensemble mean of the effective englacial sediment layer heights approaches 2 cm, which lies within the range of established debris thicknesses on many debris-covered glaciers (McCarthy et al., 2022), neglecting the englacial and basal debris components of (H e ).
The dominant role of the grain-size parameter D 50 points to the importance of hydrological processes in mobilizing sediment from the catchment (Sections 3.2.1 and 3.2.3).Yet, the performance of water discharge alone in predicting sediment discharge from the catchment varies substantially from year to year and glacier to glacier (Figure 5; e.g., Delaney et al., 2018;Swift et al., 2005;Willis et al., 1996).Furthermore, the changing availability of sediment, particularly in the proglacial area, indicates that supply-and transport-limited sediment transport conditions should be considered in evaluating sediment transport through the catchment.For example, less sediment mobilization persists in the proglacial area at the end of the study period compared to the beginning (Figure 6).
The model results suggest that bedrock erosion rate and sediment storage below the glacier (the legacy of erosion prior to the model initialization) can be conflated over the annual timescales examined here ( ε and H g0 ; Figure 8).Yet, the strong trade-off between sediment storage below the glacier and bedrock erosion in achieving favorable model outputs demonstrates the need for their individual consideration.In another application, the use of glacier velocity data sets or a more sophisticated subglacial erosion model to temporally evolve the bedrock erosion rates may help to separate these processes in driving variations in sediment discharge (e.g., Herman et al., 2015;Ugelvig et al., 2018).However, more complex erosion parameterizations are not applied here, as they require the addition of a large number of parameters, which are also poorly constrained, such as glacier sliding velocity and erosional exponents and constants required in "erosion rule" formulations (Herman et al., 2015;Koppes et al., 2015).Furthermore, the model's lumped nature means that a single or mean bedrock erosion rate can be used, unlike in a distributed model, where spatial heterogeneities in sediment production and transport capacity impact outputs (e.g., Delaney & Anderson, 2022).

Conclusions
We present a numerical model of sediment transport in glacierized catchments which considers the interacting processes related to glacier retreat, subglacial sediment transport, proglacial sediment transport, bedrock erosion, existing sediment in the catchment, and meltout of englacial material.Fluvial subglacial and proglacial sediment transport is considered in supply-and transport-limited regimes.We run this model iteratively with different parameter combinations and compare the outputs to suspended sediment discharge data collected over 7 years from the highly glacierized Fieschergletscher catchment in the Swiss Alps.Retained model runs capture seasonal variations in sediment discharge and many peak events observed in the data.The model outputs show that the proglacial area can act as both a sediment source and sink.By examining the interaction of parameters in the retained model runs and the behavior of the model's representation of geomorphic processes, we quantify the interaction among different processes controlling the sediment discharge from the catchment.
In other applications, inverting the model with additional data from different locations of the catchment may help to better evaluate geomorphic processes and their fluxes.For instance, applying the model to additional measurements of sediment discharge directly from the glacier snout, or to height changes in the proglacial area evaluated from DEMs, might yield more nuanced insights into a catchment's sediment dynamic.Implementing similar modeling frameworks in other catchments will help to understand the difficult-to-observe processes driving variability in sediment transport from source to sink.

Figure 2 .
Figure 2. Maps and photo of Fieschergletscher modified from Felix et al. (2022) (aerial images: Swiss Federal Office of Topography).(a) Overview of catchment boundary (orange) and waterway (dotted blue) of the hydropower plant Fieschertal (image collected in 2021).(b) Glacier terminus with historical glacier outlines (Linsbauer et al., 2021) and water intake, that is, measurement location in the present study.The shaded area represents the gullied section of the proglacial area.(c) Picture of the proglacial area (VAW, 18 July 2014), camera location and direction indicated with red diamond and arrow in (b), and (d) keymap.

Figure 1 .
Figure 1.Sketch of geomorphic processes acting in a catchment.The model presented in this study lumps individual glacial and proglacial processes over their respective areas to describe sediment transport from the catchment.Note that rockfall and mass wasting processes are omitted from the model although shown in the sketch.

Figure 5 .
Figure5.Cumulative amounts of measured water discharge and sediment discharge from the measurement station at the Fieschertal water intake are calculated from daily totals (a).The table denotes Pearson's correlation between sums of water discharge Q w and sediment discharge Q over 5-day periods in the first column and Pearson's correlation between sums of sediment discharge Q and suspended sediment concentration SSC over 5-day periods in the second column.The third column denotes the fraction of the time over which sediment transport is increasing with respect to water discharge (i.e., concave-up curvature), also over 5-day periods.Normalized cumulative sediment and water discharge over the season to 1 (b).The black line represents a 1-to-1 relationship between sediment discharge Q and Q w .

Figure 6 .
Figure 6.Observed water flux and modeled sediment fluxes in the catchment resulting from the accepted model runs (Section 3.3).(a) Water discharge from measurement station (i.e., Figure S2 in Supporting Information S1).(b) Subglacial sediment discharge (Q g ) in red with arrows denoting values exceeding axis limits.(c) Modeled sediment flux out of the proglacial area in blue (Q p ) and measurements in orange (Q).Emergent sediment flux (Q e ) into the proglacial area in green, has a different scale than the other sediment fluxes.Colored areas represent the maximum and minimum quantities of the ensemble of retained model outputs at each timestep.Panels (d)-(f) show magnified values for 2017.

Figure 7 .
Figure 7. Deposition of sediment dynamics in the model.(a) Ratio of sediment entering the proglacial area from the glacier to sediment leaving the proglacial area.The gray area below 1 denotes greater sediment transport from the proglacial area (Q p ) compared to the glacier (Q g ), resulting from the mobilization of sediment in the proglacial area.White area denotes reduced sediment transport in the proglacial area, resulting in sediment deposition.Periods with less than 50 m 3 of sediment discharge per day have been removed from the plot to avoid extreme ratios when minimal sediment transport occurs.(b) Sediment layer height in the proglacial area.(c) Sediment layer height in the subglacial area with median values in black lines.Shading represents the percentage of runs that fall into bins of till height at each time step.

Figure 8 .
Figure 8. Correlation plot with parameters tested (Table3).The lower left of the figure is scatter plots of parameters that resulted in retained model runs against each other.On the diagonal lays histograms of parameters in retained model runs.The red line denotes the mean retained parameter value, and gray area is the prior distribution.The upper right panels of the plot are Pearson's correlation for the parameters against each other.Parameter correlations that are significant at the 2σ confidence level are shades of red or blue, with color intensity increasing with absolute correlation.The axis limits span the sampling range presented in Table3.The median sediment size (D 50 ) parameter represents the transportability of sediment.ε controls the production of sediment.Initial till condition (H g0 ), initial sediment condition (H p0 ), and effective englacial debris thickness (H e ) all represent the storage of sediment in parts of the catchment prior to the model run.

Table 1
Model Variables

Table 2
Physical Model Parameters and Constants Note.Parameters with blank values are used in the parameter optimization and presented in Section 3.3.aFordescription, refer to the subglacial hydraulic model presented inDelaney et al. (2019).

Table 3
Model Parameters Used for Inversion

Table 4
Annual Values of Sediment Transport and Hydrological Characteristics at Fieschertal Note.Complete time series are given in Figure S2 in Supporting Information S1. a Maximum hourly values of measured variables.DELANEY ET AL.