Fluvial Response to Climate Change in the Pacific Northwest: Skeena River Discharge and Sediment Yield

Abstract

Changes in climate affect the hydrological regime of rivers worldwide and differ with geographic location and basin characteristics. Such changes within a basin are captured in the flux of water and sediment at river mouths, which can impact coastal productivity and development. Here, we model discharge and sediment yield of the Skeena River, a significant river in British Columbia, Canada. We use HydroTrend 3.0, two global climate models (GCMs), and two representative concentration pathways (RCPs) to model changes in fluvial fluxes related to climate change until the end of the century. Contributions of sediment to the river from glaciers decreases throughout the century, while basin-wide overland and instream contributions driven by precipitation increase. Bedload, though increased compared to the period (1981–2010), is on a decreasing trajectory by the end of the century. For overall yield, the model simulations suggest conflicting results, with those GCMs that predict higher increases in precipitation and temperature predicting an increase in total (suspended and bedload) sediment yield by up to 10% in some scenarios, and those predicting more moderate increases predicting a decrease in yield by as much as 20%. The model results highlight the complexity of sediment conveyance in rivers within British Columbia and present the first comprehensive investigation into the sediment fluxes of this understudied river system.

1. Introduction

The majority of sediment delivered to coastal oceans arrives via rivers [1]. Riverine sediment influx affects coastal landscape evolution, estuarine and marine ecosystems, and has received increasing attention due its role in moderating global carbon cycles [2,3]. In global estimates of sediment yield, rivers with mountainous catchments deserve particular attention because they deliver a high amount of sediment, disproportional to their basin size in comparison to non-orogenous zones [4]. However, mountainous sub-catchments are typically lacking in climate and hydrological monitoring stations, which complicates estimating sediment yield for these rivers.

The amount of sediment delivered from rivers to coastal areas is shifting from pre-industrial levels due to changes in anthropogenic land use, climate patterns, and sea level rise [1]. There is also growing evidence that rates of warming are amplified with increasing elevations [5], leading to decreased snowpack accumulation and increased glacial melt in mountainous river catchments [6,7]. The 21st century is anticipated to be a critical period for alpine glaciers and their watersheds. Research within the European Alps has predicted that glaciers within the Alps will largely disappear by the end of the 21st century with a strong acceleration in ice loss occurring after 1980 [8]. Along the mountainous coast of British Columbia (BC), rivers are still adapting to the legacy of the last glacial maximum and sediment yield is often increased due to continuous erosion of glacial till deposits within the river valleys [9,10]. Estimates of future changes to sediment yield in this context have not been made. Not only do we therefore need to extend modelling studies of sediment yield to those understudied mountainous catchments, but also study changes thereof in the future. Thus, within this research we examine changes in river discharge and sediment yield in a mountainous catchment in northern BC (Skeena Watershed) through numerical modelling.

The Skeena River is the second largest river in BC and plays a key role in the sediment and nutrient dynamics of the Skeena Estuary. The estuary is of high ecological and socio-economic value to the province (e.g., [11,12]), yet remains understudied. As the Skeena is highly glacially influenced, the timing and nature of glacial melt and subsequent flux over the next century have implications for watershed, estuary, and fisheries management [13,14].

Using the established HydroTrend 3.0 numerical model [15] and climate change scenarios from Global Climate Models (GCMs), we estimate water discharge and sediment yield and changes thereof for the Skeena River throughout the 21st century. Modeling up to the year 2100 covers the period of anticipated largest glacial melt and provides policy makers and managers with future projections of watershed change. We answer the following questions: What is the current amount of water discharge and sediment yield at the Skeena River mouth? How does this compare to the local context? How are water discharge and sediment yield expected to change by 2100?

2. Methods and Study Overview

2.1. The Skeena River, British Columbia, Canada

The Skeena River (Figure 1) empties into the Pacific Ocean after flowing 570 km south and southwest through the Skeena, Hazel, and Coast Mountains [12]. The Skeena River catchment covers an area of ~54,432 km² [12] with 50% of the basin area at elevations above 1000 m. About 1197 km² of the catchment are covered by glaciers, most of these located in the Coast Mountains [12]. The interior river valley is characterized by a wide (5–20 km) floodplain through low-grade metamorphic and sedimentary rocks. The coastal portion of the river valley downstream of Terrace (Figure 1b) is flanked by the steep, high-grade granitoid and gneiss complexes associated with the Coast-Cascade belt [16,17]. River valley fill reflects the recent glacial history of the Skeena watershed and consists of alluvial sediment overlying glacial till, especially in the interior portions of the catchment [18].

The climate governing the watershed varies between the interior and coastal portions of the watershed (Figure 2). The coastal sub-basin (SBC) receives higher precipitation and more moderate winter temperatures than the interior sub-basin (SBI). This is because high fall and winter precipitation occurs frequently in the coastal portion, while the interior portion is located in the rain shadow of the Coast Mountains [12,19]. Mean annual temperature in the Skeena River watershed has increased by 1.1 °C over the 20th century, particularly during the winter months [20]. During the same period, the average precipitation has increased by 10% and the glacial area has decreased by 6% at an estimated volume loss rate of −15.4 km³ yr⁻¹ [20].

The Skeena River is gauged for discharge and water level in the interior portion of the watershed, but hydrometric stations on the main river branch of the river in the coastal area are missing. The hydrometric station furthest downstream is at Usk (for location, see Figure 1), which is located over 100 km upstream of the tidal limit defined by [12] and excludes ~22% of the watershed area. The measured annual mean discharge at Usk is ~900 m³ s⁻¹ (rounded to the nearest 10 m³ s⁻¹ using a 1981–2010 average [21]) and the highest flood event recorded (up to 2020) was in 1948 at 9340 m³ s⁻¹ [21]. Along with a multitude of smaller creeks, roughly fourteen river tributaries join the Skeena River downstream of Usk [12]. These Coast Mountain tributaries are located in an area with substantial rainfall and retreating glaciers and potentially increase discharge and sediment load downstream of Usk significantly. The most comprehensive estimate of mean annual Skeena discharge is 2157 m³ s⁻¹ based on monitoring data at Usk plus the contribution of four downstream gauged tributaries [22]. A past estimate of sediment load for the Skeena River is 11 G kg yr⁻¹ with a yield of 260 t km⁻² yr⁻¹ using a watershed area of 42,000 km² (corresponding to the area above Usk) [23]. An estimate of sediment discharge for the entire watershed is missing.

2.2. HydroTrend 3.0 Model Overview

HydroTrend 3.0 (herein referred to as HT or HydroTrend as suitable) uses a water balance model and semi-empirical relationships to estimate suspended sediment load, bedload, and water discharge at a river outlet based on inputs provided by the user that describe the watershed’s physical properties and climate [15]. The model has been successfully applied on studies across the globe, and is capable of reproducing a range of magnitudes in sediment loads and discharges [24,25,26].

Calculations are made over the total watershed area and are not tracked throughout the river network. Thus, results are produced at the drainage outlet. Input values for the model rely primarily on annual and monthly averages of input parameters of basin-wide conditions over any user-defined duration of interest (1–1000 s of years). Annual and monthly standard deviation (for climate inputs), distribution skewness (precipitation), and change per annum (for glacial equilibrium line altitude, or ELA) can be provided to account for temporal variation of the averaged inputs for the calculation of annual, monthly, and daily values.

As described further in Kettner and Syvitski [15], HydroTrend discharge and bedload are calculated at a daily time scale that can be post-processed further into user specified time intervals (e.g., monthly, annual). Discharge is calculated using a water balance that simulates precipitation, storage, reduction (e.g., evaporation), and release through components of rain (𝑄_rain), snow melt (𝑄_nival), glacial melt (𝑄_ICE), groundwater discharge (𝑄_ground), and total evaporation (𝑄_eva) for the entire basin. Hypsometry and lapse rate along with climate inputs control a time-variable freeze line altitude used to estimate the proportion of precipitation stored as snow as well as the volume of glacial ice. Subsurface storage, transpiration, hydraulic conductivity of the dominant lithology cover are also considered as model inputs (see Appendix A section). A power law regression algorithm is applied within the model to estimate more extreme flood [15].

Bedload is calculated using a modified Bagnold’s equation:

$$Q_{B_d}=\left(\frac{{\rho }_s}{{\rho }_s-\rho }\right) \frac{\rho g{Q}_d^{\beta }{S}_l{e}_b}{g tan \lambda }when u\ge \mathrm{when}{u}_{cr}$$

(1)

where daily bedload ($Q_{B_d}$) at the river outlet is calculated using the river bed slope at the river mouth (S_l), daily discharge (𝑄_d), sediment density (${\rho }_s$), fluid density ($\rho$), limiting angle of repose of sediments on the river bed ($\lambda$), and bedload rating ($e_b$) and efficiency ($\beta$) terms under the condition that stream velocity ($u$) exceeds the critical velocity needed to initiate bedload transport ($u_{cr}$) [15].

Unlike the daily discharge and bedload, daily suspended sediment load ($Q_{S_d}$) is calculated by first deriving the long-term suspended sediment load ($Q_S$) for the user-defined simulation duration (1–1000’s of years)) contributed from overland/instream flow ($Q_{S_{BQART}}$) and glaciers ($Q_{S_G}$) [15]. Reservoir retention (trapping and residence time) of suspended sediment is calculated based on reservoir geometry and applied to the long-term calculation prior to deriving daily values. Using values of daily discharge over the simulation duration, a psi model calculates the inter- and intra-annual variability to apply to the long-term suspended sediment load to derive daily $Q_{S_d}$. Monthly values can then be derived within HydroTrend from the daily value calculations in post-processing [15].

Overland and in-stream contributions ($Q_{S_{BQART}}$) to suspended sediment are derived empirically through a basin-wide approach by multiplication of long-term discharge (Q), maximum relief (R), drainage basin area (A), long-term basin-averaged temperature (T), B (a product of the lithology factor, trapping efficiency of reservoirs, and anthropogenic factor), and a coefficient of proportionality $\varpi$ (0.02 kg s⁻¹ km⁻² °C⁻¹) [15]:

$$Q_{S_{BQART}}=\varpi B{Q}^{-0.31}{A}^{0.5}RT$$

(2)

Equation (2) is intended for basins with annual temperatures over or equal to 2 °C (see [15] and [27] for further details on the BQART method).

Glacial contributions ($Q_{S_G}$) are derived through [15]:

$$Q_{S_G}=\left(1-\frac{V_{S_g}}{P_g}\right)\frac{{{\displaystyle \sum }}_{a=1}^n\left(1.93 \ast {10}^{-3}A{\left(9.8\right)}^b-Q_{S_{BQART}}\right) }{n}$$

(3)

where ($V_{S_g}$) is the ratio of precipitation turned to ice (based on the initial temperature, lapse rate, and hypsometry), ($P_g$) the total precipitation fallen on the glacial area, A is the total drainage area, b is the logarithm of the percentage of glaciers in the basin, and n is the total number of simulation years. Through a back calculation, $Q_{S_{BQART}}$ is subtracted from an estimate of annual suspended flux based on watershed area and percentage of glaciers. The subtraction results are then averaged over the number of simulated years and multiplied by a ratio that determines the proportion of precipitation fallen, but not stored on the glacier, and which is available to transport sediment off the glacier. If all precipitation falling over the entire glacial area is stored, the ratio equals 1: the glacier is not melting and there is no discharge of water/sediment from the glacier ($Q_{S_G}$. = 0) [15].

For the purposes of clarity within the text, discharge (Q), suspended load ($Q_S$), bedload ($Q_B$), and suspended sediment concentration ($C_S$) results will be described using different subscript abbreviations depending on whether the presented results are referring to a (d) day (eg: $Q_d$/$Q_{S_d}$/$Q_{B_d}$/etc.), a (m) month (eg: $Q_m$), an (a) annum (eg: $Q_a$), or an ($\overline{}$)average of thirty years (eg: $\overline{Q_m}/\overline{Q_a}/$etc.) of HydroTrend outputs. $S_C$ refers to measured suspended sediment concentration and $C_S$ refers specifically to model-derived suspended sediment concentration. 𝑄_ICE, 𝑄_rain, 𝑄_nival, 𝑄_eva, $Q_{S_{BQART}}$, and $Q_{S_G}$_, are the mean annual discharge derived from glacial melt, rain, snow melt, and total evaporation as well as the mean annual sediment load derived from the BQART and glacial contributions, respectively. Ground water contributions (𝑄_ground) within model outputs are split between base flow contributions that are relatively constant to the river from the ground water table (𝑄_base) and additional subsurface flow contributions (𝑄_sub) that will vary depending on groundwater storage and subsurface storm flow.

2.3. HydroTrend Application for the Skeena: Approach, Inputs and Model Calibration

The results presented in this research are a cumulation of over 48 HydroTrend model runs (per sub-basin, epoch, climate input, and RCP) used to simulate changes in Skeena discharge and load. For the sub-basins, two different model set-ups were used for the calculation of suspended load and bedload of the Skeena River: (1) a one-basin approach (1-basin) was used to estimate bedload and (2) a sub-basin combination approach (2-sub-basins) was used to estimate suspended load. Discharge was calculated through both approaches for means of comparison. The rationale behind the 2-sub-basins approach is that the substantially contrasting climates and lithologies (key elements of Equations (2) and (3)) of the coastal and interior portions of the watershed are better accounted for when splitting the watershed into two basins, as opposed to averaging the contrasting conditions. The basins were divided at Usk, located on the eastern flank of the Coast Mountains, yielding an interior sub-basin (SBI) of 42,360 km² and a coastal sub-basin (SBC) of 12,050 km² (shown in Figure 1B). Results are than added to derive the suspended load at the river mouth. Using the hydrometric station data at Usk [21], validation of the SBI model output is possible. The sub-basin combination approach relies on the assumption that all material present at the sub-basin outlet will be carried to the final drainage sink. This can be plausible for the finer material; yet, it is well known that the coarser grains will change transport modes and eventually come to rest along the river course as slope decreases towards the outlet. With no known constraints on the proportion of bedload that falls out of transport between Usk and the river mouth, a simple addition of SBI and SBC bedload is problematic. Therefore, a 1-basin method, using the entire Skeena River reach and the outlet bed slope, is utilized for the calculation of bedload in this study.

HydroTrend results and simulations are run for the last climate normal (1981–2010) (the ‘reference period’), the current normal (2011–2040), and future climate normals (2041–2070, 2071–2100). Thirty-year epochs were chosen in alignment with the ECCC standard for calculating and releasing measured climate and discharge data. The last ECCC normal was computed and released from 1981–2010 (used during the validation of HydroTrend results) and the current normal including 2011 onward has not been released, as this period is still ongoing.

Each model run (per sub-basin or future epoch) requires a 47-line input file (Hydro.IN) along with basin hypsometry (Hydro.HYPS). For brevity, a list of the input variables, a brief explanation, and sources are given in the Appendix. A more detailed spreadsheet that includes details on the calibration using the long-term discharge data and the inputs chosen for each sub-basin is provided in the Supplementary Materials S1. All HydroTrend input files are also available on a GitHub repository (https://github.com/awild21/SkeenaHydroTrend, accessed on 23 December 2022 ).

Within Hydro.IN, only three input lines (temperature, precipitation (annual and monthly), and ELA) change between the future and reference simulations. For example, ELA change per year is used to calculate an evolving ELA over time at the start of each new epoch based on the reference ELA for each sub-basin (see Appendix and supplementary S1). Aside from climate and ELA inputs, all other variables (~44 input lines) remain constant across time, but vary depending on the sub-basin. For example, changes to groundwater storage over time are unknown (lacking data predictions) and thus, the same starting values for each sub-basin are used for the initial condition, maximum, and minimum ground water storage (described in the Appendix) for running each future scenario. HydroTrend over the simulation duration will depart from the initial storage condition (usual less than 1 month run up time) and will flag if the maximum or minimum groundwater condition were reached (no flags generated). Climate inputs have a dominant influence on the model outputs, are contrasted in the results, and are discussed in the following paragraphs.

For the reference period, Environment and Climate Change Canada (ECCC) [19] provides daily measured climate station data at multiple locations across the Skeena watershed. ECCC precipitation and temperature were averaged spatially across the sub-basin (monthly) or nearest to the river mouth (annually) as required for the model input. All temperatures should be input to the model at the elevation of the river outlet since a lapse rate is applied within the model [15]. Monthly, spatially averaged temperature data were corrected down to mean sea level (SBC) or to the base level of the river at Usk (SBI) using a lapse rate derived for the Skeena (see Supplementary Materials S1 and Appendix for further details). For simulations into the future, when no ECCC data are available, climate data from GCMs were extracted, averaged spatially (river mouth or sub-basin) and averaged into 30-year epochs. A lapse rate correction was performed on (sub)basin-averaged monthly data, following our approach for the ECCC data (see supplementary material S1).

When predicting temperature and precipitation in the current and future epochs (2011–2100), each GCM produces different climate scenarios based on different representative concentration pathways (RCPs) (see [28,29] for further details on RCPs and greenhouse gas trajectories). All of the GCMs presented are coupled earth system models, part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) using the OASIS3 coupler, that differ by using different models and resolutions for simulating the Earth’s atmosphere, oceans, sea ice, continent, and fluxes [30]. The Pacific Climate Impacts Consortium (PCIC) [30]. Provides statistically downscaled daily climate scenarios for Western Canada (including minimum/maximum temperature and precipitation) at a gridded resolution of 300 arcseconds for up to 27 GCMs and multiple RCPs (for detailed methods and comparison of downscaling methods used by PCIC see [30,31]). PCIC orders the GCMs for a given region based on the Giorgi map regions [32], where the first and second GCM listed will offer the most contrast possible with diminishing differences between scenarios further down the list [30]. Twelve GCMs listed for the western North America region by PCIC capture 90% of the variance in temperature and precipitation predictions under RCP 4.5 [33]. The top three GCMs recommended for capturing the climate of the Skeena region were also among the top four GCMs used within a study in southwestern British Columbia to predict changes within the water resources of the Fraser River [34]. These GCMs were: CNRM_CM5-r1(referred to as CNRM), CANESM2-r1 (referred to as CANESM), and ACCESS1–0-r (referred to as ACCESS) (refer to [35,36,37] for details on the GCM model couplings and details for ACCESS, CANESM, and CNRM respectively). We retrieved climate data for these three GCMs and two RCPs, 4.5 and 8.5. A moderate (4.5) scenario was chosen to represent the more likely conditions [38] while the extreme, ‘worst case’ scenario [8.5] is shown for contrast. We compare model results for all three GCMs with model results using measured climate data in the validation section.

Model calibration involved running HydroTrend over SBI using the measured ECCC climate station inputs for the reference period with a range of inputs for a given variable and conducing a two-tailed test of the mean against the long-term 30-year-mean monthly ($\overline{{\mathrm{Q}}_{\mathrm{m}}}$) Usk hydrometric station results. The input variables that yielded the highest p-values were selected, or the most-cited for the region achieving a p-value over 0.9. Within the calibration for the Skeena, next to climate inputs (temperature and precipitation), glacial inputs and groundwater parameters had the largest impact on the discharge results (see input and calibration Supplementary Materials S1). All variables chosen using the long-term discharge yielded a p-value between measured and modelled discharge of more than 0.9. However, some model inputs do not affect the discharge and rather are used for the computation of suspended sediment and bedload. Due to lacking data, sediment was not included in the calibration, but will be compared to the available data from the Canadian Water Office [21] during the validation section and in a calculation of uncertainty/error around the sediment measurements (see Supplementary Materials S2).

To assess the uncertainty of the input variables concerning sediment loads, tests were conducted to examine the influence of these variables on computed suspended and bed load. One of the most important input variables with uncertainty controlling suspended sediment load derived from the BQART equation is the lithology factor [15,27]. According to the map by [27] (shown in the input Supplementary Materials S1), the Skeena watershed lies within a lithological factor of 1 (coast) and 2 (interior), due to the differing geology (in agreement with [17] descriptions). We adapted the values that are dominant on the [27] map for SBC and SBI, respectively. However, since there are limited sediment data for validation and it is unclear how representative they are, in the results, we will indicate what influence changing the factor by 0.5 (the resolution of the [27] lithological map) would have on our estimates of suspended load on the river mouth (shown in greater detail in the validation and uncertainty Supplementary Materials S2).

One of the most important input variables controlling bedload estimates is the delta plain gradient or bed slope at the river mouth. One estimate exists of 0.5 m km⁻¹ [12] for the bed slope closest to the SBI outlet (near Terrace). This estimate, however, is located a distance upstream from the tidal limit [12] and from the actual river mouth in the tidally influenced estuary (see Supplementary Materials S1 for a map of the tidal limit). Measurements in Google Earth near the Skeena River tidal limit yielded slopes of ~0.333 m km⁻¹. We have elected to use this value for our model simulations, but will indicate in our results what influence reducing or increasing the slope by 0.1 m km⁻¹ has on our final bedload estimates (presented in detail in Supplementary Materials S2).

3. Validation

3.1. SBI Discharge and Sediment: Reference Period Model Predictions vs. Usk Measurements

Following the calibration of the model, 30-year measurements of discharge at Usk are compared to simulations of discharge for SBI using i) ECCC climate station data and ii) GCM input data. Very limited sediment data are available for the Skeena River. Measurements of $S_C$ at Usk station were conducted sporadically over only 13 days in various seasons in 1988–1992 by the Water Office Canada and are publicly available [21]. On these days, measurements were taken at 12 min intervals at five depths within the channel. To derive one estimate of ‘daily’ $S_C$, all measurements were averaged across depths for a given day. Measurements at Usk were collected during a range of discharges between 156 and 3290 m³ s⁻¹, representing freshet and non-freshet conditions providing a range of validation data suitable to test the model. In order to compare these values to daily modelled $C_S{}_d$, all model values for the same month over the reference period under similar discharge (+/−50 m³ s⁻¹) were extracted, averaged, and the standard deviation determined. For example, if one measurement was collected on May 15, 1990 at a discharge of 1000 m³ s⁻¹, all $C_S{}_d$ modelled for days in May 1981–2010 with a discharge of 1000 m³ s⁻¹ +/−50 m³ s⁻¹ were extracted. This approach was used because 30-year average climate data was input to the model and a psi model used to produce variability between months and years (based on the 30-year mean), which prohibits an exact timing of 𝐶𝑠𝑑 to a given measurement day.

All model runs produce higher $\overline{Q_a}$ than that measured at Usk (~904 m s⁻³) during the same period (1981–2010) (Figure 3A). Estimates of $\overline{Q_a}$ increase for runs ACCESS, CANESM, to CNRM. This order of increase holds for model runs of the coastal sub-basin (SBC), for which no measurement data are available (Figure 3A). Comparison of measured and modelled $\overline{Q_m}$ (Figure 3C) shows overprediction of the annual river freshet by all model runs, except run ECCC, which predicts mean peak flows on average in May, rather than June. GCM Model runs overpredict $\overline{Q_m}$ from April–October and underpredict $\overline{Q_m}$ in the winter months (Figure 3C). A seven-year subset of measured and modelled $Q_m$ (Figure 3D) shows that the underprediction may result because the GCM model runs are failing to predict a secondary peak in $Q_m$ in the fall. Run ECCC on the other hand, which follows measured $Q_m$ during this time period more closely, is at times capturing the secondary fall peak, and is also capturing the accurate timing of the freshet in some years (Figure 3D). A comparison of the frequency of daily discharge (𝑄_d) events from measured and modelled data (Figure 3B) shows that the general overprediction of $\overline{Q_m}$and $Q_m$, results from an overprediction of intermediate daily discharges (2000–5000 m³ s⁻¹, as seen by the ‘hump’ in the otherwise exponential curves) and prediction of daily discharges > 8000 m³ s⁻¹ that do not occur in the measurements. All models underpredict the frequency of daily discharges under 2000 m³ s⁻¹.

A two-tailed t-test comparing measured and modelled $\overline{Q_m}$showed that p-values were 0.57, 0.5, and 0.38 for ACCESS, CANESM and CNRM, respectively. Out of the model runs, ECCC is closest to the measurements, with a p-value of 0.98 (Figure 4A). ECCC overestimated $\overline{Q_a}$ by only 1%. GCMs overpredicted the Usk discharge by 29% (260 m³ s⁻¹), 37% (330 m³ s⁻¹), 51% (460 m³ s⁻¹) for ACCESS, CANESM, and CNRM respectively (see Supplementary Materials S2 for values). Since the CNRM scenario overpredicted the discharge by over 50%, we do not present future scenarios using CNRM in the results.

The lack of sediment concentration data is a limitation of this study and thus model results of daily $C_S{}_d$ are presented with error bars within one standard deviation to see if the measured values lie within this range (Figure 4B). Since data are limited, having the highest p-value capturing the mean does not mean that this is the most accurate representation of the long-term river data. However, the model should be able to capture the measured data within one standard deviation. All measured values fell within one standard deviation of the ECCC and ACCESS inputs using a lithology factor of 2. Modelled $C_S{}_d$ from run ECCC within one standard deviation reached into ranges well above the observed measurements, while the GCM runs tended to underpredict $S_C$ (Figure 4). The GCM models performed more poorly under the lithology factor of 1.5 whereas the ECC data achieved a higher p-value (see table in Figure 4B). In addition, the ECCC and GCM model results were better able to capture the measured values within one standard deviation range under factor 2 (see Supplementary Materials S2). Despite a lower p-value for the ECCC scenario, results for SBI using a lithology factor of 2 will be presented for the current and future 2011–2100 epochs since it matches more closely with the lithology factor map for SBI [27]. The difference between this, and a lithology factor of 1.5, is displayed in brackets as error/uncertainty within the historical results.

Thirty-year mean annual suspended sediment load $(\overline{Q_S{}_a}$) (using a lithology factor of 2) are very close (205 kg s⁻¹ from ECCC to 221 kg s⁻¹ from CNRM in Supplementary Materials S2) among all model runs, likely because an overprediction of the frequency of intermediate discharges, paired with an underestimation of concentrations, leads to overall similar load estimates. However, there are no long-term suspended load data to offer comparison. Out of the runs using GCM climate data, run ACCESS is seen to perform best, followed after by run CANESM and run CNRM. Between the ECCC and GCM inputs, reducing the SBI lithology factor from 2 to 1.5 produced a reduction in suspended load of 30 (ECCC)-34 (CNRM) kg s⁻¹ (Supplementary Materials S2). Between all model runs, the lithology factor of 1.5 produces a 17–18% decrease in sediment load from the lithology factor 2 results. SBC has only a 7% higher difference between the lithology factor 1 and 1.5 scenarios (19 kg s⁻¹ for 245–264 kg s⁻¹ using ECCC) (Supplementary Materials S2).

3.2. Skeena Watershed GCM Climate Predictions

Basin-averaged, monthly temperature and precipitation normals (30-year averages) from GCMs for the reference (1981–2010), current (2011–2040) and future (2041–2070, 2071–2100) periods are shown in Figure 5 for the entire Skeena watershed. Also shown are the measured climate normals as derived from ECCC stations. During the reference period, temperature as derived from the ACCESS model appears consistently lower except for the winter months, while CANESM temperatures agree well with measured data (Figure 5, left). This is reflected in p-values of 0.64 and 0.87 for ACCESS and CANESM compared to the ECCC reference climate data, respectively. Precipitation from CANESM is, at times, up to 15 mm higher than ECCC data, while ACCESS precipitation fluctuates around the ECCC data, resulting in a higher p-value (0.94) than CANESM (0.64).

Both ACCESS and CANESM models predict increasing temperatures under RCP 4.5 and 8.5 scenarios, with maximum summer temperature increasing by 2–4° (RCP 4.5) or 5–6° (RCP 8.5) and minimum winter temperature increasing by 1° (RCP 4.5) or 3–4° (RCP 8.5). Both GCMs predict generally drier summers (1–5% drier for CANESM and 16–29% for ACCESS), while fall/winter precipitation increases (31–44% in CANESM and 10–17% in ACCESS). The CANESM model predicts the development of a more prominent peak in mid-summer precipitation as well as a more variable fall/winter precipitation compared to the ACCESS model (Figure 5). Differences between model outputs using RCP 4.5 or RCP 8.5 are generally higher temperatures, and more pronounced precipitation peaks under RCP 8.5, though shifts in timing also occur.

4. Results

4.1. Skeena River Mouth Reference Discharge and Sediment Load

Modelled discharge and sediment yield estimates at the Skeena River mouth for the reference (1981–2010) period are summarized in

Table 1. Modelled reference (1981–2010) discharge and sediment yield (30-year averages) at the Skeena River mouth using measured ECCC climate normals to drive HydroTrend simulations. Results are presented for sub-basin simulations and whole basin approach. Bracketed values are used to show error/uncertainty around estimates with the greatest uncertainty. Results in the Table 1 are displayed to the nearest 10 m³ s⁻¹ and 0.1 G kg yr⁻¹. Values used to calculate the best estimate for the Skeena mouth are italicized (SBC and SBI combined for suspended load 1-basin used for bedload).

	Basin Area	$\mathbf{Mean}\mathbf{Annual}\mathbf{Discharge}(\overline{{\mathit{Q}}_{\mathit{a}}})$	$\mathbf{Mean}\mathbf{Annual}\mathbf{Glacial}\mathbf{Suspended}\mathbf{Load}(\overline{{\mathit{Q}}_{{\mathit{S}}_{{\mathit{G}}_{\mathit{a}}}}})$	$\mathbf{BQART}\mathbf{Mean}\mathbf{Annual}\mathbf{Suspended}\mathbf{Load}(\overline{{\mathit{Q}}_{{\mathit{S}}_{\mathit{B}\mathit{Q}\mathit{A}\mathit{R}{\mathit{T}}_{\mathit{a}}}}})$ (Lithology 2 ( −0.5) for SBI & 1 (+ 0.5) for SBC)	$\mathbf{Total}\mathbf{Mean}\mathbf{Annual}\mathbf{Suspended}\mathbf{Load}(\overline{{\mathit{Q}}_{{\mathit{S}}_{\mathit{a}}}})$	$\mathbf{Mean}\mathbf{Annual}\mathbf{Bedload}(\overline{{\mathit{Q}}_{{\mathit{B}}_{\mathit{a}}}})$ (Outlet Slope 0.333 m/km +/−0.1 m/km)	Total, Mean Annual Sediment Load (S)	Total Sediment Yield (Ys)
Units	km2	m3 s−1	G kg yr−1	G kg yr−1	G kg yr−1	G kg yr−1	G kg yr−1	t km−2 yr−1
Interior Sub-basin (SBI) Model Run	42,360	910	2.8	3.9 (to 2.9)	6.7 (to 5.7)	3.6	10.3	240
Coastal Sub-basin (SBC) Model Run	12,050	580	6.8	1.2 (to 1.8)	8.0 (to 8.6)	1.5 (1.1 to 2)	9.5	790
Skeena whole basin Model Run	54,410	1570	22.4	4.1 (2.7 to 5.4)	26.5 (25 to 28)	4.2 (5.5 to 2.9)	30.7	560
Best Estimate for the Skeena Mouth	54,410	1490–1570	9.6	5.1 (+0.6/−1)	14.7 (+0.6/−1)	4.2 (+/−1.3)	18.9 (+1.9/−2.3)	350 (+30/−40)

. The difference between using a slightly (within 0.5) lower (SBI: 2–0.5) or higher (SBC: 1 + 0.5) lithology factor based on the nearest values representing the area on the [27] map are shown in brackets for suspended load. Increasing or decreasing the river outlet gradient by 0.1 m km⁻¹ is shown as a potential uncertainty in brackets surrounding the bedload estimate. Overall, the model results suggest that the larger, interior basin contributes more to mean annual discharge ($\overline{Q_a}$), but less to mean annual suspended load ($\overline{Q_{S_a}}$) than the smaller coastal basin. This is because of the modelled high contribution to suspended load from glacial sources ($\overline{Q_{S_{G_a}}}$) in the coastal basin. Mean annual bedload ($\overline{Q_{B_a}})$ derived from the 1-basin approach contributes to about ~22% of the total, mean annual sediment load (S) at the Skeena River mouth as per model result. Total Sediment Yield (Y_s) for the watershed is estimated to be ~350 t km⁻² yr⁻¹. Values used to compute the best Skeena estimate within Table 1 are italicized and consist of the combination of SBI and SBC results/uncertainty for suspended load and the one whole basin results/uncertainty for bedload.

4.2. Skeena River Mouth Future Discharge and Sediment Load

4.2.1. Mean Annual Discharge and Sediment Load

Mean annual discharge and sediment load for all four epochs are shown in Figure 6. The largest contributions to discharge in the Skeena in the reference period stem from snowpack melt ($\overline{Q_{niva{l}_a}}$) (28–39% across all runs). All model simulations predict a decrease in $\overline{Q_{niva{l}_a}}$ by the end of the century following an initial mid-century increase, except the ACCESS RCP 4.5 run, which predicts an increase in $\overline{Q_{niva{l}_a}}$, if small (3%). Of those models predicting a decrease, RCP 8.5 models predict a higher decrease (by 15% (CANESM) and 18% (ACCESS)) than RCP 4.5 (7% (CANESM)). The lowest contributor to discharge is discharge from glacial melt ($\overline{Q_{Ic{e}_a}}$) (4–15% across all runs). This contributor sees a rapid increase in the period 2011–2040 (current) compared to the reference period and a subsequent decrease in all simulations. The initial peak is more pronounced using the RCP 8.5 scenario than using RCP 4.5. By the end of the century $\overline{Q_{Ic{e}_a}}$ remains higher than in the past (57% (ACCESS) and 54% (CANESM)) for RCP 4.5 and approaches reference values (only 3% higher (ACCESS and CANESM)) for RCP 8.5.

Apart from ACCESS run (RCP 4.5) which shows little overall change, model simulations predict an increase in $\overline{Q_{Rai{n}_a}}$ by the end of the century when compared to reference values. CANESM runs predict a more drastic increase (by 49% for RCP 4.5 and 112% for RCP 8.5) than ACCESS runs (63% for RCP 8.5). While both simulations using RCP 4.5 decrease in the last period (2071–2100), simulations using RCP 8.5 predict a continued increase in $\overline{Q_{Rai{n}_a}}$. Discharge from subsurface flow ($\overline{Q_{su{b}_a}}$) increases in all model simulations by the end of the century. CANESM runs predict a more drastic increase (27% (RCP 4.5) or 40% (RCP 8.5)) than ACCESS runs (4% (RCP 4.5) or 16% (RCP 8.5). The relative contribution of both $\overline{Q_{Rai{n}_a}}$ and $\overline{Q_{su{b}_a}}$ to overall discharge range between 15–30%, with $\overline{Q_{Rai{n}_a}}$ experiencing a larger increase in relative importance than $\overline{Q_{su{b}_a}}$. As initial model input values for groundwater are held constant, discharge from groundwater flow ($\overline{Q_{bas{e}_a}}$) shows no variability throughout the century.

Overall, mean annual discharge $\overline{(Q_a}$) is predicted to increase from reference values by the end of the century. CANESM predicts a 14% (RCP 4.5) to 22% (RCP 8.5) increase, while ACCESS predicts a 4% (RCP 4.5) to 7% (RCP 8.5) increase. Under RCP 4.5, $\overline{Q_a}$ appears to be on a decreasing trend towards the end of the century, while values seem to stabilize under RCP 8.5.

Model results for sediment load show that the contribution of suspended sediment from glaciers ($\overline{Q_{S_{G_a}}}$) decreases throughout the century, from approximately half of the total load to between 16% (CANESM) and 18% (ACCESS) of the total load under RCP 4.5 or 5% (CANESM) and 7% (ACCESS) under RCP 8.5. Suspended sediment derived from erosion in the catchment ($\overline{Q_{S_{BQAR{T}_a}}}$) increases in all model simulations throughout the century. CANESM runs predict an increase in $\overline{Q_{S_{BQAR{T}_a}}}$ by 33% (RCP 4.5) and 148% (RCP 8.5), while ACCESS runs predict an increase by 72% (RCP 4.5) and 138% (RCP 8.5). Combined, model simulations suggest that suspended load ($\overline{Q_{S_a}}$) decreases by about ~30% by the end of the century, except for CANESM run RCP 8.5, in which the increase in $\overline{Q_{S_{BQAR{T}_a}}}$ in the period 2071–2100 is able to compensate for the loss in $Q_{S_{G_a}}$ and no change is observed.

All model simulations predict a rapid increase in bedload ($\overline{Q_{B_a}}$) in the period 2011–2040 (current) compared to the reference period and a subsequent decrease, resulting in an overall small increase in $\overline{Q_{B_a}}$ by the end of the century. The predicted increase is more drastic in RCP 8.5 than RCP 4.5 and in CANESM simulations compared to ACCESS. Overall, total load (combined suspended load and bedload, $\overline{S_a}$) is predicted to increase in the period 2011–2040 (current) compared to the reference period driven by the rapid increase in bedload and only moderate decrease in suspended load. Following this peak, $\overline{S_a}$ decreases by the end of the century. Whether this leads to an overall decrease compared to the reference period differs depending on the GCM used. CANESM runs predict a contrasting ~15% decrease in $\overline{S_a}$ for RCP 4.5 and ~10% increase under RCP 8.5, driven by the increase in bedload and increasing $\overline{Q_{S_{BQAR{T}_a}}}$ under this scenario. ACCESS runs predict a decrease in $\overline{S_a}$ by up to 20% driven by the decrease in $\overline{Q_{S_a}}$.

4.2.2. Mean Monthly Discharge and Sediment Load

Analyzing the seasonal trends can provide additional details on shifts in the timing of future river discharge and load. Figure 7 depicts modelled mean monthly discharge ($\overline{Q_m}$) and suspended sediment load ($\overline{Q_{S_m}})$for the four time periods. All model runs predict a shift towards a more bimodal pattern in $\overline{Q_m}$and $\overline{Q_{S_m}}$by the end of the century due to an increase in the secondary fall peak and decrease in mid-summer values. The fall peak is more pronounced in CANESM simulations than in ACCESS simulations for both variables. Also observed in all simulations is a shift in timing of the primary peak from June to May by the end of the century. The magnitude of peak flow fluctuates under RCP 4.5 with no clear trends, but decreases for both GCM runs under RCP 8.5 following an initial increase in 2011–2040. Peak $\overline{Q_{S_m}}$does not display a similar trend.

5. Discussion

5.1. Skeena River Sediment Yield in Context

Remote sensing along the British Columbian coastline has shown that the highest sediment concentrations occur in sediment plumes at the mouths of the Fraser, Skeena, and Stikine rivers [39]. Our best model estimates give a mean annual discharge of ~1500 m³ s⁻¹ and a mean annual sediment load of 18.9 (+1.9/−2.3) G kg yr⁻¹ (Table 1). This discharge estimate is below previous estimates of 2157 m³ s⁻¹ [22], while sediment load is larger than the only previous estimate available (11 G kg yr⁻¹ [23]), likely due to the previous neglection of the coastal sub-basin.

For rivers in mostly non-glaciated watersheds in BC, [9] predicted an upper and lower envelope or ‘main trend’ of specific suspended sediment yield per basin area, invoking the importance of the legacy of glacial deposits and erosion in river valleys over modern hillslope processes. Accordingly, specific suspended sediment yield increases with drainage area up until 30,000 km² and then decreases for larger basins. For context, the well-studied Fraser River to the south (Figure 1) drains the largest basin area within British Columbia (250,000 km²) with a mean annual discharge of 3410 m³ s⁻¹, a mean annual suspended load of 17 G kg yr⁻¹ [40] (and a resulting specific suspended sediment yield of 68 t km⁻² yr⁻¹). Our modelled Skeena River (basin area of 54,410 km²) specific suspended sediment yield is 250–280 (270 best estimate) t km⁻² yr⁻¹ and lies well within the envelope of specific suspended sediment yields of this size (95–400 t km⁻² yr⁻¹, [9]). This may be surprising given the large amount of estimated modern glacially derived sediment in the Skeena (~65% of the total suspended load, most of which from the coastal sub-basin) that is usually thought to cause deviations from this main trend [41]. Calculating specific sediment yield for non-glacial sources only, yields 96 t km⁻² yr⁻¹ and is at the lower limit of the range specified by [9]. The Fraser River similarly falls at the lower end of the range specified for rivers of its size (70–200 t km⁻² yr⁻¹, [9]). Recall that an attempt was made to consider the Pleistocene glacial legacy of till deposits within the Skeena by using the higher lithology factory (more available/erodible sediment) within the interior sub-basin that may aid in ensuring that the Skeena estimates fall within the Church and Slaymaker [9] envelope. Given the topographic complexity of the Skeena River (the river travels through a secondary mountain range shortly before reaching the Pacific), these results are encouraging in that our predictions for the Skeena fall within the bounds described by [9] both when including and excluding the modern input from glaciers.

Total specific sediment yield (including bedload) is 350 t km⁻² yr⁻¹ according to our model results. Bedload accounts for ~22% (+/−6%) of this total yield, a larger proportion than the 10% commonly attributed to bedload in many rivers worldwide [42,43]. In the Pacific Northwest, the similarly sized Susita River (Alaska), for example, has a total sediment yield of 560 t km⁻² yr⁻¹ with a bedload contribution of 29 t km⁻² yr⁻¹ (5%) [44]. Susita River may be less steep in its final reaching due to a wide flood plain (visible in Google Earth) in comparison to the Skeena that is closely constrained (almost fjord-like) by mountains on either side, but this is uncertain. Higher bedload proportions do exist, particularly in mountainous catchments (>30%, [45,46]), and the Skeena River runs through the steep topography of the Coast Mountains in its final reaches. Still, our model results may be an overestimation of bedload yield and very greatly depend on where one defines the river mouth and delta slope. Thus, bedload estimates should be treated with caution in the light of uncertainties around the bed slope and tidal influence in the final river reaches.

5.2. Future Changes in Discharge and Sediment Load

Projected changes in temperature and precipitation by the end of the century are expected to influence the hydrological regime in the Skeena River watershed. Generally, changes in temperature are thought to influence the timing of runoff (e.g., timing of freshet peak), while changes in precipitation are thought to affect maximum snowpack accumulation and runoff volume [47]. The Skeena River today exhibits a hybrid hydrological regime, with a freshet discharge peak in late spring due to snow-melt and a secondary peak in late summer/fall due to increased fall precipitation in addition to glacial melt [21]. Hybrid hydrographs (opposed to a snow melt dominant hydrograph with a single prominent discharge peak) show two distinct discharge peaks, one that is snow melt (freshet) derived (typical of the interior of British Columbia) and one that is rain (fall or winter storm) driven (typical of coastal British Columbia catchments). Generally, our results predict an increase in mean annual flow by the end of the century, caused by an increase in glacial ice melt, and an increase in discharge from rain and subsurface flow. This concurs with an earlier onset of the freshet peak. These observations are consistent among all GCMs and RCPs chosen. Several modelling studies that have examined changes in the hydrological regime of the Fraser River (southern BC) have shown similar predictions where mean annual discharge increases throughout the century, accompanied by an earlier onset and decrease of the freshet peak [48,49].

Shifts in timing of the freshet peak lead to/reflect a reduction in the length of snow accumulation season. However, results vary, depending on rising precipitation and temperature fluctuations, as to how this is expressed within the discharge hydrograph and mean annual flow. Initially, Skeena discharge predictions from snow increases mid-century, associated with an increase in winter precipitation, particularly in CANESM simulations. This appears to moderate the effect of rising temperatures on snowpack accumulation. Such observations have also been made for other high elevation mountainous regions where the projected increase in precipitation is large [6,50,51]. Near the end of the century, snowpack melt follows a decreasing trend, suggesting that increased precipitation may no longer be able to offset rising temperatures [52]. The magnitude of the spring melt freshet peak is predicted to decrease only under one, more extreme RCP 8.5 scenario. Results for model runs using RCP 4.5 show no consistent trend among the epochs. Such lack of consistency is indicative of the variability in spring runoff that is strongly influenced by fluctuations in annual spring temperature and precipitation. The highest floods can be expected when high temperature anomalies in spring melt season coincide with heavy rainfall causing rain on snow events. For example, the latest flood warnings in the Skeena region of BC were issued in 2022, when a cool and wet winter/early spring season with heavy snowpack accumulation was followed by heavy rainfall and a rapid increase in temperature [53].

Skeena predictions show a shift towards increased bimodality, with an increase in secondary fall peak magnitude and a decrease in summer discharge. This is related to an increase in fall precipitation as well as on average drier and hotter summers, particularly in the CANESM simulations, for which these bimodal patterns are more distinct. The matches with recent events in the summers of 2021 and 2022 in BC have shown that this average increase in summer temperature is frequently accompanied by high temperature anomalies, or heat waves [54], that can induce additional variability in the hydrological conditions and add stress to fisheries (e.g., [55]). Within the Fraser River, winter and spring runoff are predicted to increase with ongoing climate changes while summer runoff will decrease, leading to a shift from snow-dominant to hybrid regimes in the Coast Mountains and shifts from hybrid to rain-dominated regimes in the central plateau [48]. Within the Skeena, a possible shift towards a rain-dominant regime such as observed by [48] may occur beyond the century, but is not yet prevalent throughout this century based on our model results.

Examining the contributions to discharge, over the four time periods studied, illustrates that the largest increase in discharge for most models occurs from the reference period (1981–2010) to the current epoch (2011–2040), due to the sudden change of discharge contributions from ice or glacial melt. Such changes are well underway and have already been documented in many studies [6,50]. From 2011–2040 onwards, both snowpack and glacial melt follow a decreasing trend in the model results.

Future changes in sediment load are more complex, especially in the Pacific Northwest, where glacial legacy may at times overshadow recent fluvial activity [41]. The large amount of sediment derived from glacial sources, however, makes the Skeena susceptible to climate change. All our model simulations predict a rapid decrease in glacial sediment load, a trend that will likely continue beyond this century as glaciers diminish in size. Predicted changes to the non-glacial sediment load, are not consistent among our simulations, which depend largely on differing precipitation predictions between the RCPs. Increased precipitation and warmer temperatures throughout the century lead to an increase in basin-derived suspended sediment, which has been attributed to higher runoff and higher sediment availability in other river basins in the Pacific Northwest (e.g., [56,57]). Analysis of the daily model outputs suggests that this increase is accompanied by a tendency towards more episodic rain events and more extreme fluctuations in daily discharge. These predictions are consistent with observations within the past few decades or predictions of more intense and episodic storms and atmospheric rivers followed by periods of drought in BC (e.g., [58,59,60]). The increase in basin-wide-derived suspended sediment predicted by our model may in part be driven by the more episodic nature of discharge events. This increase cannot compensate for the loss of glacial sediment, except for the simulation using CANESM GCM forcing and RCP 8.5, which predicts the largest amount of increase in rainfall, compared to the other runs. A slight increase in the overall suspended load from 2041–2070 to 2071–2100 in the simulation using ACCESS GCM forcing and RCP 8.5 may point to a delayed compensation of the loss of glacial sediment beyond this century. Bedload initially increases, but follows a decreasing trend towards the end of the century concurrent with our predictions of discharge.

These results are carefully scrutinized, a variety of scenarios are given, and a validation of results provided wherever possible; however, a modelling approach does come with assumptions and some inconsistencies which are explained here. Though the validation showed that models using GCM forcing overpredict discharge, particularly from snow for 1981–2010, the trends observed are consistent among model runs. This overestimation could lead, however, to shifts in timing of the real-world scenario compared to our observations. Glacial melt, for example, could peak and subsequently decrease earlier or later in the century and the snowpack contributions could more rapidly or more slowly decline. In addition, there are a number of factors that have not been considered in our model. Firstly, groundwater processes are directly affected by changes in temperature and precipitation, which may cause changes in the contribution of groundwater to discharge, particularly in the summer [61]. Few data are available regarding groundwater processes (except for immediately around select communities) and future changes therein for the Skeena watershed. The model predicts that base flow and subsurface flow account for ~0.3–0.35% of the current discharge and warrants further study within the watershed. Changes in discharge are also caused by changes in vegetation or soil characteristics (human- or climate-induced) and can at times moderate effects of changing temperatures and precipitation [62]. Changes in vegetation or land use will also affect riverine sediment load and need to be considered for a holistic understanding of future change [63,64,65]. For example, increased mining, forest fires, landsliding, deforestation, or agriculture could all increase $S_C$ concentrations in a river [64], but are not accounted for by the model aside from potentially changing lithology or anthropogenic constants. Beyond the century, based on the observed trajectory of a decrease in glacial load and bedload, hillslope and valley-fill erosion driven by increasing precipitation and temperature will likely become the dominant drivers of Skeena River sediment load.

6. Conclusions

• This study examined fluvial discharge and sediment load and changes thereof due to climate change of a mountainous river in the Pacific Northwest (Skeena River, BC). For 1981–2010, using measured ECCC climate station inputs, mean annual discharge of the Skeena River is predicted at 1500 m³ s⁻¹ with a mean annual suspended load of 14.7 + 0.6/−1 G kg yr⁻¹ and a sediment yield of 350 t km⁻² yr⁻¹. The largest contributor to suspended sediment is derived from glacial sediment.
• Future changes in the climate of the Skeena Basin using GCM inputs predict an increase in summer and winter temperatures, and a development towards drier summers and wetter fall/winter periods by 2100. As a consequence, mean annual discharge is projected to increase rapidly during the period 2011–2040, driven by increased ice melt. Discharge then plateaus or decreases toward the end of the century, as contributions by snow and ice decline and those from rain and subsurface flow increase. An increase in seasonality of discharge occurs, with a distinct late summer secondary (in addition to the primary freshet in spring) flow peak developing due to high fall precipitation. A shift towards earlier peak flow is observed in all model simulations. Sediment load contributed by glaciers is projected to decrease by 64–69% (RCP 4.5) to 87% (RCP 8.5) by the end of the century. Rising precipitation increases the basin-wide mobilization of sediment from overland and instream contributions, particularly towards the end of the century. Bedload increases initially but continues towards the end of the century with a downward trend matching the discharge trend.
• Our findings highlight the sensitivity of mountainous river basins to climate change and stipulate a transition towards rain-dominant fluvial regimes with a loss of glacial ice and snow pack at high altitudes. Our results do not show that such a transition is completed within the Skeena by the end of the century (some contributions from glacial ice and snow pack melt remain by the end of the century). Future work should examine whether this shift is accompanied by changes in transported grain size and stored alluvium, which in turn can alter the stability of riverine and coastal habitats.