A Comprehensive Biogeochemical Assessment of Climate‐Threatened Glacial River Headwaters on the Eastern Slopes of the Canadian Rocky Mountains

Climate change is driving the loss of alpine glaciers globally, yet investigations about the water quality of rivers stemming from them are few. Here we provide an overview assessment of a biogeochemical data set containing 200+ parameters that we collected between 2019 and 2021 from the headwaters of three such rivers (Sunwapta‐Athabasca, North Saskatchewan, and Bow) which originate from the glacierized eastern slopes of the Canadian Rocky Mountains. We used regional hydrometric data sets to accurately model discharge at our 14 sampling sites. We created a Local Meteoric Water Line (LMWL) using riverine water isotope signatures and compared it to collected regional rain, snow, and glacial ice signatures. Principal component analyses of river physicochemical measures revealed distance from glacier explained more data variability than other spatiotemporal factors (i.e., season, year, or river). Discharge, chemical concentrations, and watershed areas were then used to model site‐specific open water season yields for 25 parameters. Chemical yields followed what would generally be expected along river continuums from glacierized to montane altitudinal life zones, with landscape characteristics driving chemical sources and sinks. For instance, particulate chemical yields were generally highest near source glaciers with proglacial lakes acting as settling ponds, whereas most dissolved yields varied by parameter and site. As these headwaters continue to evolve with glacier mass loss, the data set and analyses presented here can be used as a contemporary baseline to mark future change against. Further, following this initial assessment of our data set, we encourage others to mine it for additional biogeochemical studies.

data set collected from glacial river headwaters on the eastern slopes of the Canadian Rocky Mountains • Our findings can act as a contemporary reference for future investigation of glacial headwaters in rapidly evolving alpine regions • The possibilities for data exploration of our 260,000+ measure data set are numerous and we encourage others to mine it for new studies Supporting Information: Supporting Information may be found in the online version of this article.
emerging proglacial landscapes in the forefields of retreating glaciers are known to be hotspots of biogeochemical cycling (Tranter & Wadham, 2013;Wadham et al., 2001), and as a result, meltwaters moving across these landscapes tend to have biogeochemical constituents that differ in both concentration and speciation compared to freshwaters in non-glacierized watersheds (Hood & Scott, 2008;Lafrenière & Sharp, 2004).The unique biogeochemical fingerprint of glacial meltwaters is in part due to the various sources of chemicals to glacial environments.First, atmospheric constituents can be transported to ice via both wet and dry atmospheric deposition and trapped within ice, archiving past environmental conditions (Blais et al., 2001;Boutron, 1995;Pasteris et al., 2014).Second, as glaciers move, through basal ice deformation or the advance and retreat of their terminus, their heavy bodies physically weather underlying bedrock, creating characteristic till that includes fine-grained glacial flour ideally suited for geochemical reactions (S.P. Anderson, 2007;Tranter & Wadham, 2013).As alpine glaciers melt, biogeochemical constituents that were atmospherically deposited and archived within the ice or the result of subglacial grinding could be released downriver.
The downriver impacts of glacially derived biogeochemical constituents can be substantial.For instance, organic carbon and nutrients released from receding glaciers have the potential to drastically alter downstream productivity and riverine food webs (Hood & Scott, 2008;Hood et al., 2015).Contaminants like mercury and polycyclic aromatic hydrocarbons (PAHs) released to the atmosphere via anthropogenic activities and deposited to alpine glaciers could drive in-river toxicity (Beal et al., 2015;Fernández et al., 2003;Schuster et al., 2015).Even glacial flour can contribute not only to suspended sediment loads and turbidity of meltwaters, but to downriver biogeochemical reactions (St. Pierre et al., 2019).For instance, nutrients and other water quality parameters that are integral to freshwater biogeochemical cycling can be associated with sediment, impacting bioavailability and uptake by riverine organisms.Similarly, whether trace elements are associated with sediments helps define whether they are available for microbial use downriver (Hawkings et al., 2020).Yet, as mountain glaciers atrophy and meltwater-fed rivers are increasingly reliant on inputs from spring snowmelt and stochastic summer precipitation events, both the reliability of alpine freshwater resources and their biogeochemical makeup become more uncertain (S.Anderson & Radić, 2020;Immerzeel et al., 2020;Schindler & Donahue, 2006).This makes it critical to investigate contemporary meltwater quality to establish baseline data sets that future researchers can use to monitor biogeochemical changes as mountain glaciers continue to transform under the pressures of climate change (Milner et al., 2017).
The total glacierized region of Western Canada and the US is the seventh smallest of 19 glacierized regions in the world by volume (Radić et al., 2014) and mass (Rounce et al., 2023).Given their relatively small size, these glaciers are therefore particularly sensitive to changes in temperature and atmospheric conditions (Menounos et al., 2019;Rounce et al., 2023).In part because of this, glacierized alpine regions of the Canadian Rocky Mountains have long been a point of scientific interest and inquiry with research generally focused on the rapid diminishment of the icefields and outlet glaciers there (e.g., Marshall et al., 2011;Tennant et al., 2012).For example, between 1919 and 2009, the mean area change of all outlet glaciers stemming from the Columbia Icefield was modeled at −2.4 km 2 (or −34%; Tennant & Menounos, 2013), with corresponding studies highlighting large-scale glacier termini thinning and recession during this century (see also Figure 1 for a visual of Athabasca Glacier in 1918 and2011;Ommanney, 2002;Rippin et al., 2020;Tennant & Menounos, 2013).On the eastern slopes, glacier volume is projected to decline 80%-90% by 2100 with subsequent changes to the volume and timing of meltwater and biogeochemical inputs into downstream rivers (Clarke et al., 2015;Radić et al., 2014).It is no surprise, then, that the glaciers on the eastern slopes have been shown to be past their tipping point (Marshall et al., 2011).
Relatively little has been published on the meltwater-fed rivers that flow from the glacierized regions of the eastern slopes of the Canadian Rocky Mountains, or how the biogeochemistry of these rivers change across the open water season or as they traverse the downstream proglacial environment (e.g., Lafrenière & Sharp, 2005;Staniszewska et al., 2020).As such, we established 14 sampling sites along the Sunwapta, Athabasca, North Saskatchewan, and Bow rivers, which stem from outlet glaciers on the eastern slopes in Jasper and Banff National Parks, Alberta.Between 2019 and 2021, riverine measures were collected from glacier termini up to 100 km downstream across all seasons to investigate spatiotemporal patterns in proglacial freshwater physicochemistry and quality.Given the importance of these headwaters in establishing the initial dissolved and particulate biogeochemical loads of three major Canadian watersheds, and the rapid glacier mass loss that indicates this region is presently on the brink of freshwater resource uncertainty (S.Anderson & Radić, 2020;Schindler & Donahue, 2006), water quality research in these vulnerable headwaters is especially important.(Higgs, 2023).Sampling sites SR1 and SR2 and the Sunwapta Proglacial Lake are labeled in the top right image.For site abbreviations, please see Table S1 in Supporting Information S1.
For each river sampling site, the overall objectives of this study were to: (a) produce biennial hydrographs using measured and modeled discharge data; (b) characterize water isotope signatures of river water, and compare those to water isotope signatures in rain, snow, and glacial ice; (c) quantify an array of 200+ physicochemical parameters, including in situ, basic water chemistry, nutrient, ion, trace element, contaminant, and particulate ion and trace element measures; and (d) combine hydrology and physicochemical measures to estimate chemical loads and yields.Spatial (within river transects and between rivers), seasonal, and interannual patterns are explored for each objective to determine the prominent factors driving data variability at our sampling sites.This interdisciplinary work lies at the intersection of biogeochemistry, limnology, glaciology, and hydrology and is the most comprehensive freshwater assessment ever documented for the headwater reaches of these rivers.As such, the spatiotemporal data set presented here will provide an invaluable baseline reference for those investigating changes to this rapidly changing alpine region in the future, as well as increase our ability to predict the impact of climate change on the water quality of rivers originating in alpine glacierized regions.

Site Description
The headwater regions of the Sunwapta, Athabasca, North Saskatchewan, and Bow rivers are contained within the montane cordillera ecozone and span glacierized, alpine, subalpine, and montane altitudinal life zones (Demuth & Horne, 2017).Broadly, they are underlain by calcite and dolomite bedrock, with various plagioclastic feldspars, phyllosilicates, and quartzes intermixed (Gadd, 2009;Pana & Elgr, 2013).In total, we established seven sampling sites along the Sunwapta-Athabasca River continuum, three sites along the North Saskatchewan River, and four sites along the Bow River (Figure 2; Table S1 in Supporting Information S1).
The Athabasca Glacier, a valley glacier of the Columbia Icefield, is the point of origin for the Sunwapta River (SR; Ommanney, 2002).The SR passes through a proglacial lake prior to braiding through a sparsely vegetated glacial outwash plain.The braids converge into one channel and further downstream merge into the Athabasca River (AR), which originates upriver of this juncture at the Columbia Glacier in the northwest margin of the Columbia Icefield (Ommanney, 2002).The large AR flows northeast toward our most northern sampling site 97.8 km downriver, upstream of Jasper Township (Figure 3).
The North Saskatchewan River (NSR) originates at the Saskatchewan Glacier on the southern edge of the Columbia Icefield (Ommanney, 2002).Meltwaters from the Saskatchewan Glacier, the largest outlet glacier of the Columbia Icefield, terminate into a large proglacial lake before flowing downriver through a glacial valley.From there, the NSR meanders through a prominent glacial floodplain area vegetated with low-lying shrubs and grasses before forming a single large channel that flows through montane forest until our last sampling site 46.3 km downriver (Figure 3).
The Bow River (BR) begins at the cirque Bow Glacier in the Wapta Icefield (Ommanney, 2002), flowing from proglacial Iceberg Lake, over a large waterfall, and into a braided forefield.The BR then feeds Bow Lake, a large subalpine lake surrounded by subalpine forest, before it passes through a forested wetland area.Dense montane forest lined the remainder of the headwaters of the BR toward our most southern sampling site 75.4 km downriver, upstream of Banff Township (Figure 3).
For each sampling site, the relative percent of its watershed area covered by major and minor land cover classes was determined with the Alberta Biodiversity Monitoring Institute (ABMI) Wall-to-Wall Land Cover Inventory and the relative percent of its watershed area covered by wetland type was determined with the ABMI Wetland Inventory (Table S2 in Supporting Information S1, Alberta Biodiversity Monitoring Institute, 2010Institute, , 2021)).Of the 10 land cover classes, only three (snow and ice, rock and rubble, and coniferous forest) exceeded 10% cover of the watershed area at some of our sampling sites.Other land cover classes pertaining to diverse vegetation (e.g., grassland, broadleaf forest) or anthropogenic disturbance (i.e., developed, including roads) represented less than 10% of the land cover along our sampling reaches, highlighting the barren and pristine landscape of our study region.Intuitively, the relative percent of cumulative watershed area covered by snow and ice was highest at our sampling sites nearest glaciers (41.5%-54.7%),whereas the relative percent of cumulative watershed area covered by coniferous forest was highest (29.7%-41.3%)at our most distant sampling sites on each river.The proportion of cumulative watershed area covered by rock and rubble was relatively constant (35.7%-55.3%)across river reaches, but interestingly slightly higher at the glacial outwash floodplain sites (Figure 3, Table S2 in Supporting Information S1).

River Discharge
Mean daily river discharge (m 3 s −1 ) was directly quantified at two of our 14 sampling sites (SR2 and BR3) by hydrometric gauging stations maintained by Water Survey of Canada (WSC; Water Survey of Canada, 2021).To estimate daily discharge at the remaining 12 sampling sites, we used day-specific watershed area-discharge linear regression models derived from mean daily discharge measured at nine Albertan mountain river WSC hydrometric gauging stations (Table S3 and Figure S1 in Supporting Information S1).More specifically, we downloaded WSC measured discharge from the nine gauging stations and created simple daily linear regressions for 2019 and 2020 (i.e., 730 linear regressions) in a spreadsheet program before applying the linear regression equations to the watershed areas for each sampling site and for each day of the year(s).The nine gauging stations, including those at our sampling sites SR2 and BR3, were selected because they were all within the SR-AR, NSR, and BR watersheds, had watershed areas (km 2 ) similar to those of our sampling sites, and were operational during our river physicochemistry sampling years of 2019 and 2020.Based on our calibration results (described below), one regression model included data from gauging stations where watershed areas were less than 630 km 2 ; this   S1 and S2 in Supporting Information S1.Distances and elevations are not precisely to scale.
model was used to estimate daily discharge at our study sites where watershed areas were less than 150 km 2 .The second regression model included data from all nine of the gauging stations; this model was used to estimate daily discharge at our study sites where watershed areas exceeded 150 km 2 .When the regression models estimated negative discharges or discharges below 0.1 m 3 d −1 , discharge was set to 0.1 m 3 d −1 .Watershed areas for our sampling sites were obtained directly from WSC (SR2 and BR3) or delineated using the Alberta ArcHydro Phase 2 Data Set on ArcGIS (all other sites; Alberta Environment and Parks, 2018).
To determine the accuracy of our watershed area-discharge models, modeled discharge was compared to discharge obtained at the five WSC gauging stations along the main stem of our study rivers in 2019 and 2020.These calibration sites included SR2 and BR3, as well as three gauging stations within 12.8-41.0km downriver of our sampling reaches (Table S3 in Supporting Information S1).Modeled data was compared to WSC measured data at multiple watershed area thresholds (e.g., 50, 150, 500, or 1,000 km 2 ) during our calibration process.The best fit for all sites was a threshold of 150 km 2 .Thus, the calibration site 07AA007 (Sunwapta River at Athabasca Glacier) with a watershed area less than 150 km 2 was compared to modeled discharge using the first linear model, whereas the other four calibration sites with watershed areas exceeding 150 km 2 were compared to modeled discharge using the second linear model.

Sample Collection and Chemical Concentration Analyses
Sampling sites were visited monthly in 2019 and 2020 during the open water season (OWS), beginning during snowmelt in late May/early June, through peak glacial melt in July/August, then during the receding flow period in September/October.Due to safety concerns of sampling fast flowing rivers in remote areas, all sample collection was done near the shoreline in the main flow of the river.However, turbulence and mixing are high in this high-slope region, and back-eddies, slumping shorelines, and areas where in-river sediment could be disturbed were strictly avoided, ensuring a representative sample.In braided river sections (i.e., SR3 and NSR2), the main channel was selected for sample collection.Additional samples were collected twice in winter (December 2019, January 2021) during base flow, but only at sites where (a) the sampling location was partially ice free, (b) the ice margin shelves were stable enough to support the weight of a person, and (c) a safety harness could be securely connected to an onshore anchor point.
River water sampling and analytical protocols are fully detailed in Table S4 of Supporting Information S1.In general, at each river sampling site and time, we deployed a YSI EXO2 multiparameter sonde to obtain instantaneous measures of optical dissolved oxygen (ODO; % saturation and concentration), conductivity, turbidity, pH, and temperature.Sondes were calibrated at each sampling site for ODO, and prior to each sampling campaign for all other parameters.Low ionic strength standards were not used for the calibration of pH or conductivity, and as a result, possible measurement biases may have occurred in our low ionic glacial meltwaters (Bagshaw et al., 2021).However, we calibrated often, opted for spot sampling over long term deployments, and eliminated the first 5 min of data at each site ensuring in situ equilibrium was reached, all of which help to reduce bias (Bagshaw et al., 2021).Using clean field sampling protocols, we also collected water samples for the analyses of total suspended solids (TSS); total dissolved solids (TDS); nutrients (total and dissolved phosphorus ); trace elements; water isotopes (δ 2 H-H 2 O, δ 18 O-H 2 O); particulate carbon (PC) and nitrogen (PN); dissolved inorganic and organic carbon (DIC, DOC); and the contaminants total and filtered mercury (THg, FHg), methyl-and filtered methyl-mercury (MeHg, FMeHg), and PAHs (2019 only).Henceforth, discussions of "total" chemicals include those that were not filtered prior to analyses, and thus contained non-dissolved (i.e., particulate) and dissolved (i.e., filtered) fractions.We correspondingly collected suspended sediments at each study site for total recoverable particulate cation and trace element concentrations (denoted as "parameter_PTL," for e.g., Ca 2+ _PTL).Depending on protocol, samples were either processed and preserved at the time of sampling or within 24 hr in a clean field laboratory or University of Alberta laboratory.All samples were then analyzed using standard protocols in Canadian Association for Laboratory Accreditation (CALA)-accredited and/or academic laboratories.
Integrated snow cores were collected with a stainless-steel corer once in March 2021 at peak snow accumulation.One fresh snow sample was also opportunistically collected by skimming the surface of the snowpack with a clean Whirl-Pak® bag during a snowfall event that occurred while sampling.Rain samples were collected during a multi-day regional precipitation event in June 2023 by setting out wide-mouthed clean buckets lined with clean plastic bags and placed on elevated surfaces in open areas to prevent ground-level splash back.Precipitation was originally collected for water isotope analyses and to estimate the potential influence of winter wet and dry deposition on chemical contributions to the rivers during snowmelt and rainfall; however, only water isotopes will be discussed in this text (chemical concentrations for precipitation are available in Data Availability).Water chemistry was also determined for a 10.2-m ice core extracted from Mt. Snow Dome, the apex point of the Columbia Icefield, in April 2020 by the Canadian Ice Core Lab (CICL; University of Alberta).In January 2022, the CICL further collected surficial snow samples from longitudinal transects along the Athabasca and Saskatchewan outlet glaciers that were subsequently analyzed for water isotopes (Figure S2 in Supporting Information S1).
Concentration data that fell below analytical detection limits (DL) were modified to half DL values for statistical purposes (Antweiler & Taylor, 2008;Helsel, 1990).Further, all data from four sampling times at site AR1 were eliminated from statistical analyses when inputs from an upstream tributary disproportionately influenced the biogeochemical signal (see starred circles in Figure 4).Unless otherwise stated, physicochemical statistics were performed in base R and physicochemical graphs were built using ggplot2 (R Core Team, 2022; Wickham, 2016).

Water Sources
We used water isotope data from our river, precipitation, and the Mt.Snow Dome ice core samples to define water source contributions to our sampling sites.Analytical precision was ≤0.35 for δ 18 O and ≤0.94 for δ 2 H. Linear regression statistics were performed on seasonal and bulk river water isotope signatures collected at our sampling sites throughout the study period to assess their overall difference relative to the Canadian Meteoric Water Line (CMWL; Gibson et al., 2020), as well as calculated deuterium-excess signatures (d-excess, where d = δ 2 H − 8δ 18 O; Dansgaard, 1964;Pfahl & Sodemann, 2014).Water isotope signatures and d-excess from fresh snow (n = 1) and rain (n = 7) samples were used to represent new precipitation inputs, and integrated snowpack (n = 7) and surficial glacier snow (n = 8) samples were used to represent snowpack contributions, to the mountain region.The distinction between new precipitation inputs and snowpack contributions is important because post-depositional processes are known to enrich water isotope signatures and so do not represent the water isotopes of freshly fallen snow.However, for ease of interpretation, precipitation in the results and discussion refers to all combined rain and snow samples.Water isotope signatures and d-excess obtained from the ice  S1 in Supporting Information S1.Please note different y-axis scales.
10.1029/2023JG007745 9 of 24 core (n = 360) were considered undisturbed firn or glacier ice.Freshly fallen and immediately sampled rain and snow that had not yet undergone post-depositional transformations were compared to riverine water, and elevational trends in precipitation, glacier ice, and river samples were determined to assess the impact of the mountain environment on the water isotope and d-excess signatures in our study region.
To further support the discussion of our water isotope data, water sources to our overall study region were examined by calculating 10-day air mass back trajectories with the National Oceanic and Atmospheric Administration (NOAA) HYSPLIT model using National Centers for Environmental Protection and Atmospheric Research (NCEP/NCAR) Reanalysis I data for 2019-2021 (Kalnay et al., 1996;Stein et al., 2015).Trajectories were launched twice daily (3:00 a.m. and 3:00 p.m. local time) at half the height of the planetary boundary layer at Mt. Snow Dome.Residence time density analysis (e.g., Ashbaugh et al., 1985;Miller et al., 2002) determined precipitation source regions using highly frequented trajectory pathways.Trajectory endpoints over each grid cell were summed to determine total residence time densities, then normalized by multiplying each density value by the Euclidean distance between the center of the grid cell and receptor point (Ashbaugh et al., 1985).See Criscitiello et al. (2016) for further details.

Measures and Patterns in Physicochemical Parameters
Principal Component Analyses (PCA) were used to assess measures of water quality parameters spatiotemporally (across seasons and years, between rivers, and within river transects).Parameters were included in PCAs if sampling of the parameter spanned both 2019 and 2020 and the number of analytical below-detects were less than 25% of the total number of samples taken (Antweiler & Taylor, 2008;Helsel, 1990).The R packages FactoMineR and factoextra were employed for detailed data interpretation and visualization (Kassambara & Mundt, 2020;Lê et al., 2008), with biplots of PC1-PC2 and PC1-PC3 pairs used to capture most of the variability explained by the PCA.
We were primarily interested in summarizing large-scale spatiotemporal patterns across all sampling sites and chemical parameters, so collinear variables were not removed from analyses, and PCAs were specifically chosen to address the high multicollinearity and associated variable inflation factor inherent to the measured physicochemical parameters.Importantly, we were able to exploit the collinearity for the interpolation of missing data.For the PCA run on the main physicochemical parameters that met the criteria outlined above, 1.5% of datapoints were missing (NA), which would have disproportionally reduced the number of sampling records included in the PCA by 17.2% (from n = 151 sampling records to n = 125).Thus, to circumvent the loss of 26 sampling records for use in the main physicochemical parameter PCA (i.e., Figure 6), missing untransformed concentration data was interpolated using the best correlation matrix relationships possible between parameters (Figure S3 in Supporting Information S1; correlation coefficients (r) for data interpolation ranged from 0.50 to 0.96).For example, we used Si to estimate five missing K concentrations based on their correlation relationship (r = 0.79).In total, 10 of the 23 parameters included in the main physicochemical PCAs contained some interpolated data, filling anywhere from 1 to 13 NAs (PN and DOC, respectively).
For PCAs run on cations and trace element chemical concentrations (i.e., Figures S9 and S10 in Supporting Information S1, respectively), only sites that had both complete PTL and dissolved cation and trace element concentrations were included (n = 42).At a given sampling site, if there was not enough particulate matter for analyses, we sometimes combined the particulate matter from several seasonal sampling times at that site prior to digestion and analyses (denoted as such in data set; see Data Availability Statement).Since measured particulate concentrations (e.g., mg g −1 ) were converted to liquid concentrations (e.g., mg L −1 ) using sample-specific TSS data (in mg L −1 ), even if a site had two identical particulate concentrations, the PTL result included in the PCA was unique.

Chemical Yields
Measured chemical concentrations and modeled discharge were combined to model OWS (May 1 to October 31) loads (total mass OWS −1 ) at each of our sampling sites in 2019 and 2020.Chemical parameters were included in load modeling only if the number of analytical non-detects were less than 25% of the total number of samples taken for that parameter.Site AR1 was not included in chemical load modeling due to the removal of dissolved concentration data at four time points (see starred circles in Figure 4), leading to fewer sampling records (n = 5) and the overfitting of load models.
Data from two parameters of differing concentrations (DIC in mg L −1 and NO2 − + NO3 − in μg L −1 ) and two river sites of differing morphology and size (SR1, single-channel proglacial stream and SR3, braided and mid-river) were initially run through the U.S. Geological Survey's (USGS) loadflex package in R to test which modeling approach would be best applied to all sites (n = 13) and parameters of interest (n = 25) (Appling et al., 2015).In loadflex, we generated continuous load through four built-in and customizable models: (a) rectangular interpolation, (b) linear regression, (c) LOAD ESTimator (LOADEST) regression, and (d) composite, which applied a mathematical correction to a regression method of choice (in this study (c) (Appling et al., 2015)).The four loadflex models produced comparable mean load results across parameters and sites (relative standard deviation ranged from 0.06% to 5.1%; Table S5 in Supporting Information S1), indicating little overall difference between models.However, rectangular interpolation consistently produced the largest mean standard errors when compared to the other three loadflex models and was thus eliminated from consideration.Manually calculated load datapoints (i.e., point sample concentrations multiplied by daily flow; n = 11 for the parameter-site combinations) were then overlayed on modeled loads plotted against time for the remaining three loadflex models.The linear regression model consistently resulted in the most deviation between the manual and modeled load and was therefore excluded from further load calculations.Finally, to select between the LOADEST and composite models, the Durbin Watson statistic was performed on the regression residuals to determine the degree of autocorrelation (Appling et al., 2015).Autocorrelation was a requirement of the composite model and was determined to be weak to non-existent in the regression residuals (Durbin Watson ranged from 1.5 to 2.6).Because the LOADEST regression model was determined to be the overall best model to use for load calculations across all sites and parameters of interest in this study, loads were modeled in the USGS rloadest package (Lorenz et al., 2017;Runkel et al., 2004).
We followed conservative model elimination criteria to reduce the likelihood of overfit load models.As we did not collect time-series data in this study, we did not consider models where trigonometrically transformed time was not a component (i.e., LOADEST models 3, 5, 7-9).Out of LOADEST models 1, 2, 4, and 6, we chose the one with the lowest Akaike Information Criterion (AIC), which indicated the best regression fit.We then determined if the R 2 and percent bias (Bp) for the model was greater than 50 and lower than 25, respectively.If the model with the lowest AIC did not meet the statistical criteria or it resulted in a model error, it was eliminated from use (Table S6 in Supporting Information S1).In these cases, the next best models were assessed in the same way.Once models passed these criteria, we graphed manually calculated and LOADEST modeled loads against time to determine if any large deviations in load occurred in the model.All the selected models resulted in good or moderate fits (Figure S4 in Supporting Information S1).In total, we surveyed the statistics and time versus calculated/modeled load graphs of 345 LOADEST models to determine the best loads possible in this study (13 sampling sites × 25 parameters + 20 additional models from best model elimination; Tables S6 and  S7 in Supporting Information S1).Open water season chemical yields (mass OWS −1 km −2 ) for each site were then calculated by dividing the LOADEST-derived chemical load (mass OWS −1 ) by the watershed area (km 2 ).While annual OWS chemical loads and yields were calculated for each sampling site, only yields will be reported and discussed to allow for standardized results between sites of differing watershed areas.To determine whether specific reaches along the rivers were sources or sinks of chemicals, we also calculated reach-specific yields between two neighboring sampling sites using the equation:

River Discharge
There were few overall discrepancies between WSC-measured discharge and our modeled discharge at the five calibration sites in 2019 and 2020 ( 2 adj ranged from 0.774 to 0.998; Figure S5 in Supporting Information S1).Measured and modeled discharge had the weakest statistical relationship at the calibration site with a watershed area less than 150 km 2 , station 07AA007 (Sunwapta River at Athabasca Glacier, R 2 adj = 0.774, Figure S5a in Supporting Information S1), compared to the four calibration sites with watershed areas exceeding 150 km 2 (R 2 adj ≥ 0.912, Figures S5b-S5e in Supporting Information S1).During the OWS, however, modeled discharge visually deviated from measured WSC station discharge at two of the large watershed calibration sites: (a) station 05DA009 (North Saskatchewan River at Whirlpool Point) located approximately 20 km downstream of our sampling site NSR3; and (b) station 05BB001 (Bow River in Banff) located approximately 40 km downstream of our sampling site BR4.Station 05DA009 recorded average to above average discharge from March to September 2019 (106%-127% basin-wide) and 2020 (115%-155% basin-wide) relative to long-term discharge records (Alberta Environment and Parks, 2019a, 2020a).During 2019 and 2020, our model underestimated discharge at this station (Figure S5c in Supporting Information S1), possibly due to higher-than-expected snow and glacier melt hydrologic inputs from the Howse and Mistaya, two large glacial rivers, that entered the NSR upstream of station 05DA009, although the Mistaya River was hydrometrically gauged by WSC upstream of NSR3 and included in our model (station 05DA007, Table S3 in Supporting Information S1).The southernmost WSC hydrometric gauging station in our model, station 05BB001 recorded average discharge from March to September 2019 (81%-103% basin-wide) and 2020 (87%-118% basin-wide) relative to long-term discharge records (Alberta Environment and Parks, 2019a, 2020a).During 2019 and 2020 our model overestimated discharge at station 05BB001 (Figure S5e in Supporting Information S1), potentially because of the abnormally large variability of flow noted throughout the downstream reaches of the BR basin during these 2 years (Alberta Environment and Parks, 2019a, 2020a).Regardless, because stations 05DA009 and 05BB001 were high volume downstream stations, the deviations between modeled and actual discharge described here were small (R 2 adj = 0.966 and 0.912, respectively) and were not considered an issue when applying our discharge models to the 12 ungauged sampling sites to calculate chemical loads and yields.
Biennial (January 2019-January 2021) WSC-measured and modeled hydrographs for the SR, AR, NSR and BR sampling sites are shown in Figure 4.The hydrographic profiles of our 14 sampling sites generally displayed steep rising limbs in spring, notable spring and mid-to late-summer discharge peaks, gradual receding limbs in early fall, and low winter baseflow.The roughly two-peaked, right skewed distribution of the annual hydrographs in our study system was the result of a rapid snowmelt freshet in the spring associated with the onset of warmer temperatures, and the combined influence of summer precipitation and expulsion of glacier meltwaters from supra-, en-, and sub-glacial meltwater channels that were hydrologically active in the summer and early fall (Demuth & Pietroniro, 2003;Marshall & White, 2010).The hydrology of our sampling sites is thus reflective of a characteristic mid-latitude glacierized alpine hydrograph (Comeau et al., 2009;Demuth & Pietroniro, 2003;Marshall & White, 2010).Even so, some spatiotemporal trends did emerge in our biennial hydrographs (Figure 4).Total discharge increased downstream in conjunction with the greater capacity of larger watershed areas to funnel winter melt and summer precipitation to the downstream sites.Meltwater peaks appeared slightly later in 2020, yet it was an overall higher flow year.For example, summing the total discharge of the most downstream sampling sites (i.e., AR3, NSR3, BR4), total OWS Q 2020 was 4.22 km 3 compared to total OWS Q 2019 at 3.79 km 3 .The difference in the timing and magnitude of flow in 2019 and 2020 may have been partially due to the prevailing regional meteorological conditions during those years.Winter 2020 experienced a mild February (Table S8 in Supporting Information S1; Environment and Climate Change Canada, 2019, 2020) and above normal precipitation (approximately 115%-200% when compared to the mean winter precipitation over the previous 30 years; Alberta Environment and Parks, 2020b).More winter precipitation in 2020 compared to 2019 could have delayed the onset of flow as deeper, better insulated snowpacks and glaciers required more energy to melt (Alberta Environment and Parks, 2019b).Once landscape melt was triggered, the larger snowpacks combined with higher temperatures recorded in Jasper and Banff National Parks in July and August 2020 compared to July and August 2019 (Table S8 in Supporting Information S1) likely contributed to the larger flow volumes seen in the 2020 hydrographs.
At sampling sites close to their source glacier, namely SR1, SR2, NSR1, and BR1 (Figure 2, Table S1 in Supporting Information S1), the difference in discharge between the height of annual snowmelt and glacier meltwater discharge peaks was small.Moving downstream, glacier melt peaks decreased relative to snowmelt peaks, a reflection of the glacierized area making up progressively less of the total watershed area at those sites (Table S2 in Supporting Information S1; Marshall & White, 2010).Even though glacier meltwaters make up only a small proportion of the overall annual riverine flow at downstream sites in glacierized alpine basins, that proportion increases in the peak glacial melt months of July to September.For example, on an annual basis, the mean proportion of glacier meltwaters to BR at the WSC hydrometric gauging station 05BA001 (Bow River at Lake Louise) between 1975 and 1998 was 4.2% whereas the mean July to September proportion of glacier meltwaters was over double that at 8.7% (Comeau et al., 2009).The same study also noted strong positive correlations between the total glacierized area of a watershed and the glacial melt contributions to river flow across the eastern slopes of the Canadian Rocky Mountains (Comeau et al., 2009).This correlation supports the relative decrease in summer discharge peaks (driven by glacial melt) compared to spring discharge peaks (driven by snow melt) with increasing distance downstream observed along each study river (Figure 4; also see Marshall, 2014), and suggested that the proportional glacier meltwater outputs from SR-AR and NSR watersheds, where nearly 77% of Alberta's glacierized area is located (Marshall & White, 2010), are likely far greater than glacier meltwater outputs from the less glacierized BR watershed (Comeau et al., 2009;Demuth & Pietroniro, 2003;Marshall & White, 2010;Marshall et al., 2011).
The hydrology of this region is rapidly evolving.Conservative climate models (RCP 4.5) predict that glacial mass loss in western Canada may reach 85% by 2100 (Radić et al., 2014).As glacier melt contributions to headwater river discharge continue to decrease, the relative importance of snow and summer precipitation to sustaining flow will increase (Comeau et al., 2009;Marshall & White, 2010;Marshall et al., 2011) and the right skewed distribution of hydrographs produced in deglaciated alpine regions could become more exaggerated.Yet, despite the research regarding future hydrological changes to the eastern slopes described above, far less has been documented regarding the consequent impact of waning glacier flow to riverine water physicochemistry.

Water Sources
The water isotope signatures of riverine samples ranged from −23.1 to −19.1‰ δ 18 O and −176.0 to −143.4‰ δ 2 H (n = 151), with the lowest signatures recorded in spring and the highest signatures recorded in autumn.A Local Meteoric Water Line (LMWL) based on the riverine samples was calculated to be δ 2 H = 8.0*δ 18 O + 8.9 (r Pearson = 0.98) and was nearly identical to the CMWL at δ2H = 8.0*δ18O + 8.5 (Figure 5a; Gibson et al., 2020), suggesting no major in-river fractionations such as evaporation or condensation were occurring during our sampling period (Beria et al., 2018;Mahindawansha et al., 2022).Minor deviations from our riverine LMWL and the CMWL were observed for river samples binned by season with the highest slope in winter (8.8, r Pearson = 0.97), followed by autumn (8.5, r Pearson = 0.98), spring (7.6, r Spearman = 0.96), and summer (7.4,r Pearson = 0.96).There were also no discernible patterns in d-excess signatures across seasons (Figure 5b).Riverine d-excess ranged from 6.2 to 11.1‰ except for NSR3 in spring of 2019 which produced a d-excess value of 4.2‰, potentially due to instantaneous evaporation during a major spring rain event (Pfahl & Sodemann, 2014).The similar slopes and d-excess signatures across seasons were not surprising given the considerable overlap and small isotopic range of seasonal riverine water isotope measures (Figures 5c and 5d).
We further collected precipitation (rain and snow) and glacial ice core samples to identify the isotope signatures of water that fed our rivers.Isotope signatures spanned a range from −25.5 to −7.3‰ δ 18 O and −201.8 to −73.2‰ δ 2 H for precipitation (n = 22) and from −30.7 to −13.8‰ δ 18 O and −231.0 to −118.9‰ δ 2 H for the 10 m glacial ice core collected on Mt.Snow Dome (n = 360).Snow isotopes tended to sit near or below the CMWL with a slope of 8.7 (r Pearson = 0.98), typical of cold environments (Mahindawansha et al., 2022), but adding rain samples to the analysis changed the precipitation slope to 7.3 (r Spearman = 0.98; Figure 5e).This slope was more similar to that produced using glacial ice core samples, which sat near or above the CMWL with a slope of 7.5 (r Spearman = 0.98; Figure 5g).The difference between the slope produced using snow samples compared to the slope produced using ice core samples may have been due to post-depositional transformation processes (Penn et al., 2023).The two most enriched water isotope signatures, which diverged from the other precipitation samples (Figure 5e), were in rain collected near the terminus of the Athabasca Glacier at the onset of a severe thunderstorm with strong katabatic winds.We also observed negative d-excess signatures in precipitation.The negative d-excess value of the snow sample could have been the product of sublimation or freeze-thaw events and the negative d-excess signatures of the rain samples could have been due to sub-cloud evaporation (Bershaw et al., 2020), as rain samples were collected and immediately processed during rain events (Figure 5f; Beria et al., 2018;Dansgaard, 1964;Hu et al., 2022;Penn et al., 2023).In contrast, the d-excess signatures of the ice core skewed more positively than the river or precipitation data and contained signatures up to 20.5‰ (Figure 5h).Detailed descriptions of atmospheric water transport, including air mass back trajectory residence time density analysis, for our study region can be found in the SI.We also explored water transformations pre-and post-deposition and how water isotopic signatures changed with elevation at our study sites and at sampled sites from atop the Columbia Icefield in the SI.Ultimately, given the complexity of water transport and transformation processes to and within mountain watersheds (Sinclair & Marshall, 2008), our limited number of precipitation samples relative to the large range of air mass sources, and the large overlap in isotope signatures of river Water source is color-coded and the gray dashed line is the Canadian Meteoric Water Line (CMWL), while the black long dashed line is the line of best fit for all the data on each graph.Seasons are abbreviated as SPR (spring; May 14-June 25), SUM (summer; July 14-September 3), AUT (autumn; October 9-14), and WINT (winter; December 20-January 29).Season specific lines of best fit are colored by season as teal (SPR), pink (SUM), brown (AUT), and blue (WINT).Rivers are abbreviated as SR, Sunwapta River; AR, Athabasca River; NSR, North Saskatchewan River; BR, Bow River.waters, precipitation, and glacier ice at our study sites, presenting a water isotope mixing model to apportion river isotopic signatures to specific water sources is outside the scope of this study.Regardless, our combined riverine, precipitation, and ice core water isotope data set greatly expands on the existing literature on, and the characterization of, water isotopes in our study region (Arendt et al., 2015;Lafrenière & Sinclair, 2011;Niu et al., 2017).

Measures and Patterns in Physicochemical Parameters
The physicochemical measures explored in the main PCA are summarized as yearly means in Table S9 of Supporting Information S1.Dissolved cations were included in the main physicochemical PCA, but trace elements were not because only few, correlated trace elements met the criteria for inclusion.The cation-and trace element-specific PCAs included only sampling points when both PTL and dissolved concentration data were available.Further, while the BR is represented in the main physicochemical PCA, it was wholly excluded from the cation and trace element PCAs because at no sampling times was enough suspended particulate material obtained for PTL digestions and analysis.For each PCA, we visualized PC scores by season and year, river, and distance from glacier (km) to assess if these spatiotemporal factors drove the variation in measures of the physicochemistry presented on the loading plots.Please see Supporting Information S1 for a further discussion on parameters included in the PCA analyses.

Main Physicochemical PCA
The first three principal components (PC) explained 67.1% of the variation in measures for the main physicochemical parameters among all our river sampling sites and times (PC1 = 37.1, PC2 = 19.2,PC3 = 10.9; Figure 6 with loadings summarized in Table S10a of Supporting Information S1).There was a distinct grouping of dissolved parameters (e.g., dissolved carbon, ions, and nutrients) and particulate and unfiltered parameters (e.g., TSS, PC, and THg) along the PC1 axis (Figure 6a).The bulk of dissolved chemical parameters were negatively correlated with pH, as might be expected due to carbonate or silicate geochemical weathering consuming CO 2 and driving up pH (Brown et al., 1996).In situ temperature was negatively correlated with ODO concentrations along PC2 in line with the known inverse relationship between oxygen solubility and temperature in aquatic systems (Figure 6a; Dodds & Whiles, 2020).PC3 further separated the nitrogen species from most of the dissolved parameters and accentuated the negative correlation between temperature and ODO concentrations (Figure 6b).
The spatiotemporal factor that captured the most variation in this PCA was distance from glacier.While mid-to far-downstream sampling site (approximately 25-100 km) scores tended to overlap, the river sampling site scores closest to their source glacier (approximately 0-25 km) clearly differed on the PC1 axis (Figures 6g and 6h).River sampling sites closest to the source glacier tended to have the highest particulate and total concentrations and the lowest dissolved concentrations, with less obvious differences among downstream sites (Figures 6g  and 6h).This aligns with the understanding that untransformed glacial meltwaters tend to have high sediment loads but are dilute in dissolved chemistry (S.P. Anderson et al., 1997), whereas various sources and sinks to the meltwater as it flows downriver across differing landscape types can alter the relative proportion of the particulate, total, and dissolved constituents in less predictable ways (Deuerling et al., 2018;Kohler et al., 2022).We also found that samples collected during spring, when river chemistry was largely driven by snowmelt runoff, stood out from samples collected in the other seasons on PC2 (Figure 6c).During spring 2019, river samples had higher concentrations of dissolved chemical parameters than during the rest of the year, whereas during spring 2020 when snowmelt runoff volume was greater than in 2019 (Figure 4), river samples had higher concentrations of particulates and particle bound chemicals, possibly because there was heavy rainfall and overland flooding during our late June 2020 sampling trip (Figure S8 in Supporting Information S1).The rivers' contact with, and mobilization of, sediments during this extreme flow event indicate that an erosion threshold was crossed and that is likely the reason we see the highest particulate and total chemical concentrations during that period.By comparison, river (Figures 6e and 6f) poorly explained differences in variation in physicochemical measures on all the PCs.

Particulate and Dissolved Cation PCA
The first three PCs explained 93.5% of the variation in PTL and dissolved cation concentrations (PC1 = 57.3,PC2 = 29.4,PC3 = 6.8; Figure S9 in Supporting Information S1 with loadings summarized in Table S10b of Supporting Information S1).Along PC1, the PTL cations grouped at a nearly 90° angle from the dissolved cations, indicating a lack of correlation between the two groups (Figure S9a in Supporting Information S1).Their similar deviation from zero along the PC1 axis (∼−5 for PTL and ∼+5 for dissolved) and similar vector lengths on PC2 suggested relatively equal contributions of these groups to the variation captured by this ordination (Figure S9a in Supporting Information S1).Yet, the PTL and dissolved cations showed a pattern along PC3 that was not captured by PC1 or PC2 (Figure S9b in Supporting Information S1).The PTL cations grouped close together on PC3, which suggested that PC3 was a poor explainer of the PTL cation concentration variability, while the dissolved cations split between Ca 2+ /Mg 2+ and Na + /K + in the top right and bottom left quadrants of the loading plot, respectively (Figure S9b in Supporting Information S1).
The bulk of PTL and dissolved data were from summer when glacial melt resulted in rivers being turbid enough to obtain enough material for PTL analysis.A sample from the outwash floodplain along the NSR, NSR2, during the high flow spring of 2020 drove much of the variation in PTL cation concentrations along PC1 and PC2 (Figures S8, S9c and S9d in Supporting Information S1).Along PC1 and PC2, the dissolved cation concentrations were more influenced by river (Figures S9e and S9f in Supporting Information S1) and distance from glacier (Figures S9g and S9h in Supporting Information S1) as the lowest concentrations occurred at SR sampling sites closest to the source glacier.Along PC3, all spatiotemporal factors clustered closely around zero and it was consequently not considered important in explaining cation concentration variation (Figure S9d, S9f, and S9h in Supporting Information S1).

Particulate and Dissolved Trace Element PCA
The trace element PCA was conducted the same way as the cation PCA with both PTL and dissolved constituents.In the PCA of PTL and dissolved trace elements, the first three PCs explained 87.2% of the variation in concentrations (PC1 = 65.9,PC2 = 11.2,PC3 = 9.4; Figure S10 in Supporting Information S1 with loadings summarized in Table S10c of Supporting Information S1).The variation captured by PC1 was largely accounted for by the short-vectored PTL samples as the dissolved trace elements did not deviate far from zero on the PC1 axis (Figures S10a and S10b in Supporting Information S1).The opposite trend was observed on PC2 where PTL trace elements remained clustered around zero on the PC2 axis and the dissolved trace elements Cr, Mo, Sr and Al, Ba, Mn were separated by PC2 (Figure S10a in Supporting Information S1).PC3 also captured more variation in dissolved species, separating Cr from the other dissolved trace elements (Figure S10b in Supporting Information S1).Hence, PC1 captured more variation in PTL trace element concentrations while PC2 and PC3 captured similar variations in dissolved trace element concentrations.
Though both river (Figure S10e in Supporting Information S1) and distance from glacier (Figure S10g in Supporting Information S1) appeared to capture some variation in dissolved trace element concentrations, particularly at AR sampling sites far from the source glacier, PC2 contributed only approximately 1/6th the amount of PC1 to the explanation of variation in our data.Thus, it was season and year that explained the most variation in PC1 versus PC2 with spring 2020 dominating the PC1 axis (Figure S10c in Supporting Information S1).Spring 2020 particulate samples along the PC1 axis again drove the bulk of variation in PC1 versus PC3 (Figure S10d in Supporting Information S1).Apart from Cr, dissolved trace elements vectored away from SR sampling sites (Figure S10f in Supporting Information S1) close to the source glacier (Figure S10h in Supporting Information S1) where the lowest dissolved trace elements would exist.This trend was presumably due to the high quantity of snow and glacier meltwaters diluting dissolved chemistry at glacier forefield sampling sites (S.P. Anderson, 2007;S. P. Anderson et al., 1997).

Chemical Yields
Chemical concentrations in themselves only represent a sample-specific mass of chemical in a given volume of water and thus inherently reduces the discussion to specific sampling sites and times.Because snow and glacier melt-fed systems are characterized by high summer flow that drive downstream mobilization of chemicals (S.P. Anderson et al., 1997), to determine the watershed-scale load of our sampled chemical parameters, we now discuss chemical yields.
Year-specific 2019 and 2020 OWS chemical yields for each sampling site (except AR1, please see Section 2.6 for details) and parameter of study are summarized in Table S11 of Supporting Information S1.Basic chemical parameters TSS, TDS, PC, and DIC (Table S11a in Supporting Information S1) and ions Ca 2+ , Mg 2+ , and SO 4 2− (Table S10c in Supporting Information S1) had the highest chemical yields (reported in Mg OWS −1 km −2 ), while trace elements Cr and Mo (Table S11d in Supporting Information S1) and contaminants THg and FHg (Table S11e in Supporting Information S1) had the lowest (reported in g OWS −1 km −2 ).All other chemical parameters, including nutrients (Table S11b in Supporting Information S1), had moderate chemical yields (reported in kg OWS −1 km −2 ).Water yields were reported in m 3 OWS −1 km −2 (Table S11e in Supporting Information S1).With only the exceptions of Na + at BR2 (97%) and Cr at SR3 and SR4 (93% and 95%, respectively), all chemical yields were greater in 2020 than in 2019.Because year had little overall influence on the variation of the concentrations assessed in the main physicochemical PCA as indicated by the overlap of data in temporal space (Figures 6c-6d), we can likely account for the majority of the difference in yearly chemical yields from the difference in yearly discharge.Total OWS discharge of the most downstream sampling sites in 2020 exceeded that of 2019 by approximately 0.44 km 3 (see Section 3.1), resulting in larger loads and yields of chemical constituents in our study rivers.To identify spatial differences in chemical yields between our study rivers and sampling sites, mean 2019 and 2020 (±SD) OWS chemical yields for each sampling site and parameter of study were graphed (Figure 7).We also graphed mean 2019 and 2020 (±SD) OWS chemical yields between two neighboring sampling sites to assess whether discrete river reaches were sources or sinks of chemicals in our study rivers (Figure S11 in Supporting Information S1).Please note that because AR1 was removed from analysis, SR4 and AR2 are neighboring sites in Figure S11 of Supporting Information S1 and its discussion.
Along the SR-AR and BR, chemical yields of all parameters were generally highest at the sampling sites closest to the source glacier (i.e., SR1 and BR1) then decreased or stabilized at downstream sampling sites (Figure 7).Glaciers of the Canadian Rockies are warm-based, which allows for basal sliding in addition to subglacial abrasion, and potentially high sediment formation (S. P. Anderson, 2007).It is thus unsurprising that the strongest source of particulate (e.g., TSS, PN, PC) and total (e.g., TP, THg) chemicals was at the most upstream sampling sites along the SR-AR and BR, dwarfing sources of the same chemicals at downstream sites (Figure S11 in Supporting Information S1).This pattern was not consistent along the NSR which tended to have the highest chemical yields at the mid-river glacial outwash plain site, NSR2 (Figure 3).Here, river waters increased their surface area with extensive braiding and had the potential to entrain more underlying and adjacent sediments, primarily creating a  S11 of Supporting Information S1.Rivers are abbreviated as SR, Sunwapta River; AR, Athabasca River; NSR, North Saskatchewan River; BR, Bow River.For parameter abbreviations, please see Section 2.3.spike in particulate and total chemical yields (Goff & Ashmore, 1994).Along the NSR, the area between NSR1 and NSR2 was consistently the greatest source of particulate and total chemicals and contributed to the elevated chemical yields along the river, highlighting the importance of glacial outwash plains to riverine sediment loads and distribution (Figure S11 in Supporting Information S1).
The greatest decline of particulate and total chemical yields typically occurred along the upper reaches of the SR-AR and this loss was not recovered at downstream sampling sites.Between SR1 and SR2 sits Sunwapta Lake, a proglacial lake fed directly by sediment-laden glacial meltwaters flowing from the Athabasca Glacier (Figure 3).Sunwapta Lake has been shown to have extremely high sedimentation rates typical of proglacial lakes (Gilbert & Shaw, 1981;Wankiewicz, 1979), and accordingly acted as an initial settling pond for particulate and total chemicals suspended in meltwaters (Carrivick & Tweed, 2013).Importantly, the AR, NSR, and BR similarly have substantial proglacial lakes.However, the effects of sedimentation on particulate and total chemical yields are not explored there as our first sampling sites along those rivers occurred downstream of their proglacial lakes due to sampling access limitations (Figure 3).Given the dramatic decline of particulate and total chemical yields due to sedimentation in Sunwapta Lake, reason follows that the particulate and total chemical yields of the AR, NSR, and BR were likely highest between the glacier toe and their respective proglacial lakes and that this study consequently only captured some of the more chemically subdued glacial meltwaters that are possible along those rivers.For instance, Geilhausen et al., 2013 found that a proglacial lake in Austria reduced the TSS of a glacial river by over 85% during normal flow conditions and Bogen et al., 2015 found a similar 80% loss of TSS in proglacial lakes in Norway.River reach-specific chemical yields confirmed that not only did particulate and total yields decrease from SR1 and SR2, but that the area between SR1 and SR2 acted as a large sink of those chemicals (Figure S11 in Supporting Information S1).
The loss of a large portion of particulate and total chemistry in proglacial lakes may be most pronounced within the NSR.The Saskatchewan Glacier was the largest of the source glaciers in this study with an area of 38.3 km 2 and a length of 12.2 km (Tennant & Menounos, 2013).Despite NSR1 being sampled 5.6 kms downstream of the glacier tongue (Table S1 in Supporting Information S1), the relative percent snow and ice land cover was the highest of all sites at 54.7% (Table S2 in Supporting Information S1).This compares to SR1 which was sampled only 0.2 kms from the Athabasca Glacier tongue (Table S1 in Supporting Information S1) yet had a slightly lower relative percent snow and ice land cover of 50.7% (Table S2 in Supporting Information S1).Theoretically, with more glacier area in direct contact with underlying bedrock, more subglacial grinding and abrasion could occur, and more sediment-laden meltwaters could be expelled from the system.Bedrock hardness is a key factor in how much abrasion, and consequent sediment export, can occur in a subglacial system.The bedrock underlying our study glaciers is relatively consistent, with Middle Cambrian limestone and dolostones, including those of the Pika Formation (Arendt et al., 2016), underlying the Athabasca, Saskatchewan, and Bow glaciers (Pana & Elgr, 2013), with Bow Glacier also partially underlain by the Lower Cambrian Gog Group consisting of limestone and dolostone beds with some quartzose sandstone (Pana & Elgr, 2013).Yet, the chemical yields along the NSR consistently fall low-to mid-range compared to the other rivers of study, suggesting that the deposition of material in the North Saskatchewan proglacial lake was likely the highest of all the proglacial lakes in this study and that the chemical signature of the NSR downstream of the proglacial lake was markedly different from what would be measured coming from the glacier.Further, the proglacial lakes of the NSR and AR are more than twice the surface area of the proglacial lakes of the SR and BR, so size differences could play a role in sedimentation rates as well.
Bow River sampling sites BR2, BR3, and BR4 generally had the lowest particulate and total chemical yields of all rivers (Figure 6), following the low concentrations of particulate matter at these sites (Figures S9 and S10 in Supporting Information S1).The decrease of those same chemical yields from BR1 to BR2 may have been due to the presence of Bow Lake, a large subalpine lake (2.8 km 2 surface area; Blais et al., 2001) between the two sampling sites.The reduction in particulate and total chemical yields is not as apparent in Bow Lake as it was for Sunwapta Lake, likely due to its position downstream of BR's proglacial Iceberg Lake which would have been the site of the most substantial glacial sediment settling (Figure 3).The BR1-BR2 river reach-specific yields for particulate and total chemicals were neutral to slightly positive, suggesting that the area between BR1 and BR2, including Bow Lake and a small wetland, had very little impact on the chemical load into lower reaches of the BR despite a decrease in chemical yield values from BR1 to BR2 described above (Figure S11 in Supporting Information S1).Yet, our chemical yields suggest that Bow Lake impacted the dissolved chemistry.Most dissolved (e.g., TDS, DOC, TDN, all trace elements, FHg) chemical yields decreased from BR1 to BR2 and never recovered to their pre-lake levels at downstream sites (Figure 7).This trend was especially evident as the BR tended to have the highest dissolved chemical yields compared to the other rivers of study, in part due to the elevated dissolved loads estimated for BR1.As our most southern study river and the only one not to originate from the prominent Columbia Icefield, the BR was arguably our most atypical glacial river.Thick subalpine forests began to line the waters midway along Bow Lake, a relatively more upstream locale for the onset of substantial vegetation than the other rivers in this study, and wetland zones were present at downstream sites (Figure 3, Table S2 in Supporting Information S1).Nitrogen species yields decreased from BR1 to BR2 suggesting lake and/or wetland processes may have been using TDN and NO 2 − + NO 3 − faster than they could be replenished (Figure 7b and Figure S11b in Supporting Information S1), while DOC (Figure 7a and Figure S11a in Supporting Information S1) and Na + (Figure 7c and Figure S11c in Supporting Information S1) yields peaked at BR2, a site with wetland regions along the river's edge, suggesting site-specific inputs.
When chemical yields were considered by parameter grouping, few anomalies stood out.Water yields were greatest at the most upstream site and leveled out at downstream sites along each river (Figure 7e).Following this, yields of basic chemical parameters mostly decreased along the SR-AR and BR, but DOC increased from the most upstream sites to the most downstream sites along all three rivers (Figure 7a).The increase in DOC yields with distance from glacier can potentially be explained by the transition from a glacierized altitudinal life zone with fast-flowing cold, turbid meltwaters to an alpine altitudinal life zone with freshwaters more hospitable to in situ biological productivity and/or more additions of allochthonous terrestrial organic matter (Robison et al., 2023).Dissolved Si yields also increased with distance downstream and was the only nutrient to do so (Figure 7b).An increase in water temperature and terrestrial-sourced organic acids can have a stimulatory effect on the erosion of silicates from surrounding soils and minerals, which may account for the higher Si yields at downstream sites (S.P. Anderson, 2005).Further, amorphous silica bound to suspended sediment has been found in glacial regions (Pryer et al., 2020), and the dissolution of bound silica along our rivers could have also contributed to higher Si yields at downstream sites.And while related ions such as Ca 2+ and Mg 2+ paralleled each other's downstream yield patterns (Figure 7c), trace element yields were found to be more element-and site-specific potentially due to the influence of geology or in situ ion exchange processes, for example (Figure 7d).Ultimately, sources and sinks of dissolved chemicals varied with parameter and site (Figure S11 in Supporting Information S1).Apart from the area between SR1 and SR2 acting as a minor sink for Ba, Mo, and FHg, river reaches were sources of dissolved chemicals.For instance, every river reach was a source of DOC, though as we would expect, downstream watershed DOC contributions were greater than upstream ones (Figure S11a in Supporting Information S1); Ca 2+ contributions across all rivers and river reaches were steady (Figure S11c in Supporting Information S1); and the strongest sources of Sr were midstream for the SR-AR, mid-and down-stream for the NSR, and upstream for the BR, and thus river-specific (Figure S11d in Supporting Information S1).In summary, the patterns in yields explored here were intuitively impacted by the diverse landscape features of our sampling reach (Figure 3;Dixon & Thorn, 2005).Proglacial lakes acted as settling ponds for turbid meltwaters coming from glacier drainage systems (Gilbert & Shaw, 1981;Wankiewicz, 1979), while braided glacial outwash floodplains contributed more particulate and total chemicals to the rivers (Goff & Ashmore, 1994).As you move downriver, the warmer climate, development of soils, and introduction of more montane vegetation in adjacent off-channel areas increased riverine yields of key dissolved parameters.

Conclusions
The objectives of our research in the glacial headwaters of the Sunwapta-Athabasca, North Saskatchewan and Bow rivers on the eastern slopes of the Canadian Rocky Mountains were to: (a) produce biennial hydrographs; (b) characterize water isotope signatures; (c) provide a general analysis of a 200+ parameter physicochemical data set; and (d) estimate chemical yields for each sampling site and assess the spatiotemporal factors that drove data variability.We found that we could accurately model the hydrology at our sampling sites using existing regional discharge data sets, and that the hydrographs we produced showed the roughly bimodal shape typical in a snow-and glacier-melt hydrological regime during the open water season.Back trajectories from 2019 to 2021 showed that air masses traveled to our study region from western North America, circumpolar Arctic, and the North Pacific Ocean, explaining some of the diversity in precipitation and glacial ice water isotope signatures that we measured.The LMWL produced from our riverine water isotope signatures was similar to the CMWL (Gibson et al., 2020) and reflected a mix of inputs.Principal component analyses of physicochemical measures revealed that of all the spatiotemporal factors examined, distance from glacier had the overall greatest impact on the variability of physicochemistry in our system.Yet, that may rapidly change as the glaciers disappear from the landscape and extreme precipitation events, such as those experienced in June 2020, become more prevalent.Chemical yields followed what would generally be expected for rivers flowing through glacierized, alpine, subalpine, and montane altitudinal life zones.For example, particulate and total yields were highest at sampling sites nearest their source glacier and along glacial outwash plains.And apart from DOC and Si, whose yields increased with distance downstream, dissolved yields were parameter and site specific.
As mountain glaciers disappear due to climate change on a global scale, proglacial landscapes are expanding past their current areal boundaries and evolving biogeochemically.Further, the hydrology of the rivers that once stemmed from glaciers will become increasingly reliant on snowmelt and summer precipitation for their flow.Hence, chemical yields will likely become more reflective of non-glacierized systems, making it especially important to provide detailed scientific evaluations in these regions.We provided an initial high-level assessment of such a data set here, at the intersection of biogeochemistry, limnology, glaciology, and hydrology.Thus, the analyses explored here can also be considered a contemporary baseline reference for those investigating river headwaters stemming from the eastern slopes of the Canadian Rocky Mountains in the future.However, the data set is expansive with over 260,000 datapoints, and as such, we encourage other to mine and analyze it in ways we did not here.For example, water apportionment via isotope mixing models could be investigated using the water isotope signatures of our river water, precipitation, and glacier ice.Further, data on dissolved phosphorus concentrations, MeHg, and PAHs were not discussed in detail here, yet they could provide new biogeochemical information potentially important for gauging productivity and toxicity in glacier rivers.We believe the possibilities for data exploration are numerous and support further investigations given the lack of biogeochemical studies in these river systems and the collective effort in producing this data set.

Figure 1 .
Figure 1.Side-by-side comparisons of historic and contemporary photographs of Athabasca Glacier, Jasper National Park.All 1918 scenes were captured by Arthur Oliver Wheeler during an Interprovincial Boundary Survey and archived by the Mountain Legacy Project.All 2011 photographs were captured by the Mountain Legacy Project team(Higgs, 2023).Sampling sites SR1 and SR2 and the Sunwapta Proglacial Lake are labeled in the top right image.For site abbreviations, please see TableS1in Supporting Information S1.

Figure 2 .
Figure 2. Map of source glaciers, our sampling sites on the Sunwapta (SR), Athabasca (AR), North Saskatchewan (NSR), and Bow (BR) rivers, and Water Survey of Canada (WSC) hydrometric gauging stations that were used as hydrology calibration sites in this study (see TablesS1 and S3in Supporting Information S1).Sampled watersheds are color-coded and match the watersheds depicted on the inset map of Alberta, Canada.All sampling was conducted in the headwaters of these major river systems in Jasper and Banff National Parks (Park boundaries are not shown).Base map imagery from Esri, Government of Canada (Alberta Project NA1) 1983 Transverse Mercator.

Figure 3 .
Figure 3. Profiles of the (a) Sunwapta, (b) Athabasca, (c) North Saskatchewan, and (d) Bow rivers with sampling sites, distance from glacier, watershed area, elevation, and notable river features included.Relative percent watershed area covered by major land cover classes in our study system (snow and ice [snowflake], coniferous forest [tree], and rock and rubble [mountain]) are listed for each site.For further site descriptions, including river abbreviations, exact sampling locations, and relative percent land covered by minor land cover classes, please see TablesS1 and S2in Supporting Information S1.Distances and elevations are not precisely to scale.

Figure 4 .
Figure 4. Modeled (solid line) and measured Water Survey of Canada (WSC; dotted line) discharge (Q; m 3 s −1 ) and physicochemical sampling dates (colored circles) at the 14 sampling sites along the (a) Sunwapta (SR), (b) Athabasca (AR), (c) North Saskatchewan (NSR), and (d) Bow (BR) rivers for 2019 through early 2021.Stars in the orange circles for AR1 symbolize sampling dates where dissolved concentration data was eliminated from all data analyses.For site abbreviations, please see TableS1in Supporting Information S1.Please note different y-axis scales.

Figure 5 .
Figure 5. (a and c) The δ 18 O and δ 2 H of water for river, (e) precipitation, and (g) the Mt.Snow Dome ice core samples (left column), and (b and d) d-excess and δ 18 Oof water for river, (f) precipitation, and (h) the Mt.Snow Dome ice core samples (right column).Water source is color-coded and the gray dashed line is the Canadian Meteoric Water Line (CMWL), while the black long dashed line is the line of best fit for all the data on each graph.Seasons are abbreviated as SPR (spring; May 14-June 25), SUM (summer; July 14-September 3), AUT (autumn; October 9-14), and WINT (winter; December 20-January 29).Season specific lines of best fit are colored by season as teal (SPR), pink (SUM), brown (AUT), and blue (WINT).Rivers are abbreviated as SR, Sunwapta River; AR, Athabasca River; NSR, North Saskatchewan River; BR, Bow River.

Figure 6 .
Figure 6.Principal Component Analyses (PCA) of the main physicochemical parameter measures quantified among all our river sampling sites and times.Shown are biplots of PC1 versus PC2 (left column) and PC1 versus PC3 (right column), with parameters and scaled vectors (a and b), and sampling sites color-coded spatiotemporally by season and year (c and d), river (e and f), and distance from glacier (km; g and h) for visual comparison.The first three PCs account for 67.1% of the variation in measures of the main physicochemical parameters.Seasons are abbreviated as SPR (spring; May 14-June 25), SUM (summer; July 14-September 3), AUT (autumn; October 9-14), and WINT (winter; December 20-January 29).Rivers are abbreviated as SR, Sunwapta River; AR, Athabasca River; NSR, North Saskatchewan River; BR, Bow River.For parameter abbreviations, please see Section 2.3.