Sub-regions
The study region focuses on 18 sub-regions within the Pacific mountain ranges of North American overlapping with the range of Pacific Salmon with >1.5% glacier cover (Figs. 1 and 2). The term “sub-region” here refers to either a single major salmon watershed or aggregates of small coastal watersheds, which range in area from ~13,000 to ~68,000 km2. For sub-regions within Alaska, USA, we accessed boundary data from the Watershed Boundary Database at the USGS (https://www.usgs.gov/). For sub-regions within British Columbia, Canada, we accessed boundary data from the Freshwater Atlas of British Columbia (https://catalogue.data.gov.bc.ca/). Pacific salmon range data were from the National Center for Ecological Analysis and Synthesis (Fig. 1). The study region covers ~623,000 km2 across British Columbia, Canada and Alaska, USA and ~20% of the total North American range of Pacific salmon.
Glacier outlines
Outlines for the 45,963 glaciers within the study region were obtained from the Randolph Glacier Inventory v6.0 (https://www.glims.org/RGI/; RGI v6.0), which provides a globally complete data set of glacier outlines outside of Greenland and Antarctic ice sheets17. These glaciers cover a total area of ~81,000 km2, which corresponds to 80% of the total glacier area in the Pacific mountain ranges within North America. The glacier outlines refer roughly to the years 2009 ± 2 for Alaska, and 2004 ± 5 for Western Canada17,53. Glacierization for each of 18 sub-regions ranges from 1.5 to 52%.
Present-day streams
Synthetic stream networks were constructed from Digital Elevation Models (DEMs) for each of the 18 sub-regions using Geographic Information Systems (GIS; ArcGIS 10.6 and QGIS 2.18) hydrology tools to represent present-day streams throughout the study region. Specifically, we used Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEMs v2.0 with a spatial resolution of ~30 m54. Open access synthetic stream network datasets such as the National Hydrography Dataset (NHD) from the USGS and the Freshwater Atlas from the British Columbia government are available but were not used due to inconsistencies in spatial resolution across the study region. From our synthetic stream networks, we eliminated all stream segments that overlapped with the RGI glacier outlines because the ASTER global DEMs used to create the synthetic stream networks represent glacier surface elevation rather than estimated deglaciated terrain. All present-day streams within our study region are void of any major dams that inhibit salmon movement based on existing databases of dams55. To summarize present-day stream kms, and all subsequent analyses, we used rstudio: 1.4.1103-4, R: ‘Mirrors’.
Identifying and verifying stream gradient thresholds for migrating salmon and for determining accessible glaciers
We used stream gradient-based thresholds the determine constraints in salmon migration and the number of glaciers that would be accessible and create future streams for migrating adult salmon. Based on the large body of literature suggesting stream gradients (e.g., ranging from <10–20%) suitable for migrating Pacific salmon25,26,27,29, and our own validation results (see below), we applied a conservative stream gradient threshold of 10%, and more inclusive alternative of 15%, to the synthetic present-day stream networks representing streams suitable for migrating adult salmon. The salmon migration constrained synthetic present-day stream network was used to identify which glaciers would be accessible to adult salmon. The 15% stream gradient threshold represents accessibility for salmon that are capable of swimming up steeper gradients (e.g., Chinook salmon). We calculated stream gradient within the present-day stream network by breaking the network into ~500 m segments, then determining the slope of each stream segment by extracting elevation values from the ASTER global DEMs for both ends of each segment, then dividing the elevation difference by the segment length. The ~500 m stream segment length was the minimum distance possible given the spatial resolution of our study region and memory limitations with the ArcGIS software (see Uncertainties). The upper limits of salmon migration within present-day streams were identified by selecting contiguous stream segments in the direction from the river mouth to headwaters (i.e., glacier tongue) that were below the different stream gradient thresholds (10 and 15%).
The migration stream gradient thresholds selected for this analysis were determined based on our validation analysis and those presented in the literature. For the validation analysis, we used Pacific salmon location data from the Anadromous Waters Catalogue (AWC; www.adfg.alaska.gov/sf/SARR/AWC) and segment slope values from our synthetic present-day stream network, for a sample watershed, the Susitna River watershed. The AWC is maintained by the State of Alaska and contains observation data for all Pacific salmon species with associated spatial reference (latitude and longitude). However, given the difficulties inherent in surveying for fish throughout Alaska, the AWC is incomplete, and Pacific salmon are likely to present in many unsurveyed streams. Thus, our analysis of fish presence within streams represents a minimum of the extent of fish occurrence and therefore is conservative. We used the Susitna River watershed for our validation analysis given that it supports good coverage of all five species of Pacific salmon, with each species having a minimum of ~200 individual observations (Supplementary Fig. 1). In addition, the Susitna River watershed is large in size (Supplementary Table 1), containing a range of watershed types (i.e., rain, snow, and glacier dominated) representing diversity in salmon habitats56.
We ran the GIS network analysis tool to assess the maximum stream gradient each Pacific salmon could migrate beyond. Using our present-day stream network (segmented into ~500 m stream lengths containing slope values), we obtained the maximum stream gradient value crossed by each salmon species (obtained from the AWC) from river outlet to observation location (Supplementary Fig. 1). For each salmon species, apart from coho, ~75% of the salmon observations cross a maximum stream gradient threshold of at least 10%, with Chinook, chum, and pink salmon only infrequently being observed in stream segments beyond ≥10% (Supplementary Fig. 1). Given the migration thresholds identified in the literature25,26,27,29 and this verification analysis, we confirmed that the migration stream gradient thresholds of 10% represent a conservative estimate, and 15% represent a more inclusive estimate. Our analysis covering 623,000 km2 should be viewed as the best-available predictions that require field validation at specific locations.
After determining the salmon migration thresholds, we selected the glaciers that were butted against the salmon migration constrained present-day stream network using the 10% and 15% thresholds. In other words, the selected glaciers that feed into streams that have the potential to currently be colonized by migrating salmon. Thus, as the glaciers retreat, they present new habitat that salmon can colonize. We define these glaciers as “accessible”.
Future salmon-accessible streams created from estimated deglaciated bedrock
From each of the accessible glaciers, we created future stream networks from estimated deglaciated bedrock. We estimated the deglaciated bedrock terrain by subtracting gridded ice thickness data from ice surface DEMs for every accessible glacier within our study region. Ice thickness distribution was calculated at a grid resolution of 25–200 m (depending on glacier area) using a simple dynamic model that considers glacier mass turnover and ice flow mechanics, and by inverting the glaciers’ surface topography57 (described below). The data set uses glacier outlines from RGI v6.0 and is in close agreement with the recently released global consensus glacier ice thickness product (see Methods)58. Additionally, we compared the ice thicknesses and overall ice volumes used in this analysis59, with a more recent ice thickness data set58, and our assessment indicates that the differences are small (<3%) with respect to the other uncertainties that we account for. Inferred ice thickness was validated against a set of ice-penetrating radar observations for 300 glaciers from most glacierized regions of the world57. Ice surface DEMs were obtained from the Shuttle Radar Topography Mission DEMs for glaciers below 60°N (resolution of ~90 m), and from ASTER global DEMs for glaciers above 60°N (resolution of ~30 m). Both DEMs had an elevation uncertainty of ~±10–20 m for mountain areas57.
We built a future stream network derived from estimated deglaciated bedrock terrain beneath accessible glaciers using the ArcGIS hydrology toolbox. We extracted slope values from each ~500 m segment using the same methods as used for the present-day stream network. To each segment, we assigned a stream order to each stream segment using the program RivEX to determine each stream size. Stream order is a metric used to measure the relative size of streams, where the smallest tributaries are referred to as first-order streams, which flow into larger streams that combine to form streams of higher order. All future stream segments derived from deglaciated bedrock terrain beneath accessible glaciers below the adult migration thresholds (<10 and 15%) of all sizes (stream orders First–Fourth in this case) were termed “future salmon-accessible streams”. The absolute future salmon-accessible stream kms were summed for each of the 18 sub-regions (Supplementary Table 3) to determine the extent of new streams suitable for migrating adult salmon. Last, we calculated the relative increase in stream kms for each of the 18 sub-regions by dividing the total future salmon-accessible stream kms, projected to be created in 2100, by the total present-day stream kms below either 10% (Fig. 2) or 15% (Supplementary Fig. 2) stream gradient threshold.
Glacier retreat modeling and exposure of future salmon-accessible streams
To project the retreat of glaciers, we applied the Global Glacier Evolution Model (GloGEM) to each of the accessible glaciers that requires the use of DEMs and glacier outlines (RGI v6.0)23. GloGEM computes glacier mass balance and associated geometry changes for each individual glacier in the study region23. To calculate glacier surface mass balance, as a difference between accumulation (snowfall and refreezing) and ablation (glacier surface melting), GloGEM was forced with a monthly time series of near-surface air temperature and precipitation. Glacier geometry changes (e.g., thinning and/or shrinking) were assessed from the surface mass balance coupled with the empirically derived functions of glacier thinning along the glacier centerline60. For annual mass losses at marine- or lake-terminating glacier fronts, glacier retreat was approximated by accounting for glacier front height and width61. The total mass changes for each glacier were used to adjust surface elevation and extent on a yearly basis. Previous work gives further details of the model, its calibration, and downscaling procedures23,39.
To project the glacier retreat for the benchmark years 2050 and 2100, we forced the glacier model with temperature and precipitation time series from an ensemble of five GCMs. The five GCM models were selected as they showed better performance in simulating climatology over North America relative to other Coupled Model Intercomparison Project 5 GCMs:62 CanESM2, CSIRO-Mk3-6-0, GFDL-CM3, MIROC-ESM, and MPI-ESM-LR. The GCMs are subjected to a range of specified climate forcings that correspond to plausible scenarios for the rate of change in the concentration of atmospheric CO2 and other greenhouse gases. For the IPCC AR5, these scenarios are referred to as Representative Concentration Pathways (RCPs) and the standard emission scenarios are RCP2.6, RCP4.5, RCP6.0, RCP8.563. For several of the GCMs we used, the RCP6.0 was omitted, and the RCP2.6 is likely to not be reached64, therefore we selected the RCP4.5 and RCP8.5 for the glacier modeling. The glacier model is forced by each GCM in the ensemble, whereas the projections of glacier retreat are presented as the ensemble mean for each RCP.
We used the glacier retreat projections to assess when each future salmon-accessible stream segment, below either the 10% or 15% stream gradient thresholds, would become exposed from glacier ice based on 10-year averages of modeled glacier extent centered around the years 2050 and 2100 for both the RCP4.5 and RCP8.5 and each of the five GCMs (Fig. 2 and Supplementary Fig 2), and the ensemble-mean. We then summed the stream kms that would be available by the projected 2050 and 2100 years for each of the 18 sub-regions (Supplementary Table 3), and the total stream kms by stream order after complete deglaciation (Supplementary Table 6).
Defining salmon habitat requirements using stream gradient and order
Based on a strong body of literature on habitat suitability (Supplementary Table 2), we determined how many of the future salmon-accessible stream kms could be used specifically for salmon spawning and rearing habitat by selecting stream segments with gradients using two scenarios: either a 0–2% or 0–4%. We also define salmon spawning and rearing habitat as having a stream order greater than first order. Although we acknowledge that salmon can use first-order streams, we focus on streams greater than first order, which includes second-, third-, and fourth-order streams. Thus, our analysis indicates that suitable spawning and rearing habitat from the future salmon-accessible stream networks had segments with stream orders that ranged from second to fourth, with stream gradients ranging from either 0–2% or 0–4%.
From the future salmon-accessible streams, we determined the amount of salmon habitat by summing the stream kms based on two scenarios (i.e., 0–2% spawning/rearing with a < 10% stream gradient threshold; and, 0–4% spawning/rearing with a <15% threshold), then determined when the stream habitat would become available for the years 2050 and 2100 for each of the 18 sub-regions (Fig. 3). The two scenarios help capture the fact that different salmon species have different tendencies in terms of stream gradients associated with spawning and rearing (Fig. 3)32.
Uncertainty estimates
Given that our study integrated models, data inputs, and analytical approaches for determining future salmon habitat created by glacier retreat across a large study region, it is important to consider potential uncertainties. All values presented throughout the manuscript with uncertainty estimates are propagated from the uncertainty derived from the (1) GCM projections, and from the mean uncertainty extracted from the results of our sensitivity analyses of (2) glacier ice thickness estimates, and (3) stream segment length, as presented below. We used the root sum of squares from the three individual uncertainty estimates to determine the global uncertainty. Uncertainty estimates are presented throughout the paper as ± one standard deviation. Given the time and computation power required to run some of these sensitivity analyses, we use sample sub-regions, North Southeast, AK, and Taku River, BC, that represent diversity in watershed size, terrain complexity, and extent of future stream network then applied the error to other sub-regions (Supplementary Table 1).
Other sources of uncertainty that were considered and discussed below, but not included in our uncertainty estimates, are the (4) glacier retreat model (GloGEM) used in our analysis, and (5) Unknown potential changes in habitat and landscape variables important to salmon spawning and rearing.
GCM projections
There are uncertainties in the IPCC GCM that project future greenhouse gas emissions, which appears to be a dominant source of error when considering rates of glacier retreat19,24. However, climate models similar to those used in this analysis have shown to be quite accurate at predicting forecasted temperature changes65. To illustrate the uncertainty of the climate models, for each of the 18 sub-regions, we present the individual five GCMs projections as one standard deviation around the ensemble mean.
Glacier ice thickness estimates
There is considerable uncertainty in modeled ice thickness distribution. An intercomparison of ice thickness models showed a ±25.9% uncertainty in inferred ice thickness58. Thus, for our analysis, there might be additional uncertainty due to the estimated bedrock elevation in the deglaciated topography. However, ice thickness models show a good performance regarding the patterns of thickness distribution (thin/thick parts of the glacier), even when the average estimated ice thickness of an individual glacier may be too high or low. Moreover, errors in estimated elevation of the bedrock are likely similar at either end of a ~500 m segment59. Hence, the errors in the calculation of slopes are likely reduced due to the spatial correlation of errors in the ice thickness estimates.
To determine the uncertainty of the future stream kms derived from deglaciated bedrock topography, obtained by subtracting ice thickness from ice surface DEMs, we ran a sensitivity test on the ice thickness data. In a conservative approach, we systematically increased/decreased all ice thicknesses by 25.9% according to the uncertainty stated in Farinotti et al. (2019) and re-computed bedrock topography. We then re-ran the ~500 m segment analysis to determine stream gradient, and applied the salmon migration constraints of 10%, and calculated the total number of future stream kms. We applied this analysis to the sample watersheds, North Southeast, AK, and Taku River, BC. Our sensitivity analysis resulted in estimates of total future stream kms for these sub-regions that were 13% greater (given error in increased ice thickness) or 12% less (decreased ice thickness) than our predictions. Thus, the error associated with the ice thickness data used to determine total future stream kms is comparable to the error from uncertain climate forcing for a given CO2-emission pathway presented in the main text.
Stream segment length
The stream gradient thresholds applied were derived from ~500 m segment lengths and a DEM with a 30 m resolution. We selected a ~500 m segment length as it was not feasible due to limited computation resources to shorten our segment length. However, there is some imprecision in using a ~500 m segment length, such as it is impossible to know the exact stream characteristics present (e.g., longer series of riffles vs. a single large waterfalls), and Pacific salmon upstream migration is generally restricted by certain stream features such as waterfalls. In addition, we could not validate our stream gradient thresholds for adult salmon migration against known barriers (e.g., waterfalls) because there are no known data sets of salmon migration barriers between southern British Columbia and Alaska. We also acknowledge that some migration barriers can vary seasonally depending on stream flows. Therefore, some segments with a 10% stream gradient may contain a migration barrier whereas others may not. For comparison, other studies have used similar reach lengths of 200 m66 and 500 m67.
To assess the uncertainty arising from the assumptions on the segment lengths used throughout the study, we ran two sensitivity analyses, (1) on the present-day stream network, which determines the number of accessible glaciers, and (2) on the future stream kilometres. First, we broke the present-day salmon migration stream network into different segment lengths (~250 m, ~400 m, ~600 m, and ~750 m segment length with a 10% stream gradient threshold), as the number of accessible glaciers is determined by the present-day stream network. We applied this sensitivity analysis on the two sample sub-regions, North Southeast, AK, and Taku River, BC. For the North Southeast, AK sub-region the number (and size) of salmon-accessible glaciers changed within the range of −16 to 9% (~90–~350 ha), depending on the chosen segment length (Supplementary Table 7). For the Taku River, BC sub-region there was the same number of salmon-accessible glaciers in all segment length scenarios (Supplementary Table 7). Second, we estimated the uncertainty of the future stream networks segment length derived from each accessible glacier’s future deglaciated terrain (see below). For this analysis, we used a ~250 m and ~750 m segment lengths, and 10% stream gradient threshold, for each of the 18 sub-regions to understand variations in total future salmon-accessible stream kms given different segment lengths (Supplementary Table 8). On average for all 18 sub-regions, the total future salmon-accessible stream kms when using the ~250 m segment length was 14% less, and when using the ~750 m segment length, was 11% more. The ~250 m segment length is more precise and therefore there are more opportunities for the slope calculation to be above the 10% stream gradient threshold, whereas the opposite is true of the ~750 m segment length. In addition, the stream network using the ~750 m segment lengths extends further into the upper reaches of the stream network, leading to more salmon-accessible stream kms. Given that for some glaciers the spatial resolution of the ice thickness data was only 200 m, that there were computational limitations, and to be consistent with the present-day segment lengths for comparison purposes, we chose to use the ~500 m segment length, a good trade-off between capturing the accuracy and precision of the stream network. The associated error in determining segment length is less than or equal to the error in using the different GCMs as presented in the main findings (Supplementary Table 3).
Glacier retreat model
We predicted rates of glacier retreat using GloGEM, as described above. A recent intercomparison of global glacier models showed that GloGEM predicts somewhat faster rates of glacier retreat compared with some other models19. However, the most recent GloGEM runs, based on a new glacier-specific data set of mass balances 2000–202020, are consistent with the results used in the present study for the two relevant regions Alaska and Western Canada. This indicates that the model correctly captures glaciers changes in the past and is, thus, likely to yield reasonable projections for the next decades.
Unknown potential changes in habitat and landscape variables important to salmon spawning and rearing
Many habitats and landscape variables (e.g., channel width and confinement, streamflow, stream temperature, riparian forest development) can be important to Pacific salmon throughout their life cycle29,68,69,70,71, but it was not possible to consider these explicitly in this study. Field measurements of variables such as sediment supply, grain size, bankfull discharge, and channel slope are not readily available over large geographic areas72,73, and are impossible to obtain for sub-glacier environments. Many studies have determined ways to extract important environmental variables for salmon using DEMs66,74. However, this type of analysis was not possible given the uncertainties and errors in estimating the deglaciated bedrock topography (see above). Therefore, only stream order and gradient were used in determining salmon streams and suitable habitats for spawning and rearing75. These metrics are useful and accurate in identifying suitable habitats for Pacific salmon, as shown by other studies75,76,77. Further, we were not able to incorporate natural hazards associated with glacier retreat78, such as landslides or river piracy42, into predictive models, but acknowledge that such stochastic events can alter the accessibility and suitability of rivers following glacier retreat.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
