Recent streamflow trends across permafrost basins of North America

Bennett, Katrina E.; Schwenk, Jon; Bachand, Claire; Gasarch, Eve; Stachelek, Jemma; Bolton, William R.; Rowland, Joel C.

doi:10.3389/frwa.2023.1099660

ORIGINAL RESEARCH article

Front. Water, 29 June 2023

Sec. Water and Critical Zone

Volume 5 - 2023 | https://doi.org/10.3389/frwa.2023.1099660

This article is part of the Research Topic Women in Critical Zone Science View all 7 articles

Recent streamflow trends across permafrost basins of North America

$\r\nKatrina E. Bennett$ Katrina E. Bennett¹^*

Claire Bachand¹

William R. Bolton²

¹Los Alamos National Laboratory, Earth and Environmental Sciences Division, Los Alamos, NM, United States
²International Arctic Research Center, University of Alaska Fairbanks, Fairbanks, AK, United States

Introduction: Climate change impacts, including changing temperatures, precipitation, and vegetation, are widely anticipated to cause major shifts to the permafrost with resulting impacts to hydro-ecosystems across the high latitudes of the globe. However, it is challenging to examine streamflow shifts in these regions owing to a paucity of data, discontinuity of records, and other issues related to data consistency and accuracy.

Methods: Recent trends for long-term periods (1990–2021, 1976–2021) in observed minimum, mean, and maximum seasonal and annual streamflow were analyzed for a range of watersheds across North America affected by varying degrees of permafrost coverage.

Results: Streamflow trend analysis revealed that areas affected by permafrost are changing variably over the periods in terms of maximum, mean, and minimum seasonal and annual streamflow. These changes indicate a significant shift occurring in the most recent 46 years towards increasing mean streamflow for the dominant (> 50%) permafrost systems. Meanwhile, minimum streamflow increases for all permafrost-dominant systems and many of the other permafrost-affected systems across the seasons and annual periods considered, with the greatest number of significant changes in streamflow over other metrics. Maximum streamflow is shifting variably with significant increases in the permafrost-dominant systems in winter and fall over longer time periods of analysis. Our analysis suggests that streamflow trends are driven by climate (precipitation, followed by temperature), while variables such as permafrost coverage only appear important in the most recent 32-year period.

Discussion: The increases in streamflow trends observed in this study are reflective of deepening active layers and thawing permafrost, indicating that the entire hydrograph is undergoing change within permafrost-dominant streamflow systems as the Arctic moves towards a warmer future under climate change. Despite the many challenges to understanding changing streamflow in cold regions, there are new products and datasets in development that are increasingly allowing researchers to better understand the patterns of change in Arctic and subarctic systems affected by permafrost, offering a range of new tools, which, along with continued observational records, may help in improved understanding of changing Arctic streamflow patterns.

1. Introduction

The globe and in particular, the Arctic and subarctic regions of the world are experiencing enhanced warming and increased movement of atmospheric and terrestrial water through the hydrologic system (Huntington, 2006; England et al., 2021; Chylek et al., 2022; Wendisch et al., 2022), along with accelerated thaw of the permafrost (ground that remains frozen for more than 2 years, Rowland et al., 2010; Lawrence et al., 2012; Chadburn et al., 2017; Nitze et al., 2018; Biskaborn et al., 2019). As these changes occur, complex processes result in deeper soil active layers (the seasonally thawed layer at the top of the permafrost) and thinner permafrost (Hayes et al., 2014), more extensive talik (i.e., holes in the permafrost) formation (Parazoo et al., 2018; Liu et al., 2021), larger areas of discontinuous permafrost (Biskaborn et al., 2019), and less ground ice (Liljedahl et al., 2016).

In the Arctic and subarctic, hydrological intensification also results in increased water flow, both off the land surface and in rivers, deeper snowpack, more summer rain, and more extreme events (Peterson et al., 2002; McClelland et al., 2006; Rawlins et al., 2010; Prowse et al., 2015). Both hydrologic cycle rates and permafrost thaw have positive feedbacks to the vegetation, water quality (i.e., heat, sediment, and carbon/mercury mobilization), and energy exchanges (i.e., water vapor feedback, permafrost extent, snowpack, evapotranspiration, and soil moisture changes) across the Arctic and subarctic landscapes and oceans (Peterson et al., 2002; Rawlins et al., 2010; Bring et al., 2016; Stuefer et al., 2017; Durocher et al., 2019; Zheng et al., 2019; Koenigk et al., 2020; Feng et al., 2021). These changes are anticipated to enhance positive feedbacks to global climate impacts, leading to even more rapid changes and further increasing greenhouse gas concentrations in the atmosphere (Schuur et al., 2015).

Arctic and subarctic river outflows to the Arctic Ocean have been noted to be increasing over the past 50 years (Peterson et al., 2002; McClelland et al., 2006; Bring et al., 2016; Durocher et al., 2019). A particular focus has been on the annual discharge of the largest Eurasian rivers (e.g., Lena, Ob, and Yenisei) as these rivers input a large percentage (55%) of fresh water annually to the Arctic Ocean compared with 28% from North American (NA) systems (Shiklomanov et al., 2021). While Eurasian rivers are unquestionably increasing in annual mean streamflow (Shiklomanov et al., 2021), if and how NA Arctic river discharges are changing remains in question. For example, in the Mackenzie River basin, no obvious trends were found in terms of changing annual flow magnitudes from 1973 to 1999, although variability was noted to be increasing (Zhang et al., 2001; Woo and Thorne, 2003). Updated work showed increases over a longer time frame (early-to-mid 1900's-2013; Rood et al., 2017). On the other hand, Déry and Wood (2005) found trends of decreasing annual streamflow discharge from 1964 to 2003 at coastal stations in Arctic Canada. An updated paper in 2009 extended the time series by 4 years (1964–2007), which revealed a trend reversal toward increasing annual streamflows and streamflow variability (Déry et al., 2009). NA river basin annual discharges to the Arctic Ocean were noted to be decreasing from 1964 to 2000 by McClelland et al. (2006) while Zhang et al. (2001) found increasing trends in some basins from 1967 to 1996. More recent work by Durocher et al. (2019) showed that annual river discharges around the Hudson, James, and Ungava bays in northern Canada increased significantly while other NA systems increased as well (although results were not significant at p-values < 0.05). Déry et al. (2016) noted increasing annual flows in northern Canada (18%) over the latter part of the 1964–2013 (1989–2013) time period. Annual discharge in 2020 from the two large NA rivers contributing to the Arctic (Mackenzie and Yukon) was ~25% greater than the 1981–2010 average, largely owing to a very wet summer in 2020 (Druckenmiller et al., 2021). Overall, the consensus is toward increasing NA river flows on the order of 1 km³ per year from the larger NA river systems (Druckenmiller et al., 2021; Shiklomanov et al., 2021).

Previous research has focused on the causal relationships for changing streamflow patterns. Some authors have linked changing streamflow to specific climate variability indexes, such as the Pacific Decadal Oscillation (Rood et al., 2017) or the Arctic Oscillation (Déry and Wood, 2004), although no significant relationships were found by Durocher et al. (2019) for either climate variability index Peak streamflow within high latitude systems is largely driven by snowmelt processes, a seasonal phenomenon that is highly influenced by climate warming (Wang et al., 2016). Peak streamflow carries with it important constituents, including sediment, carbon, oxygen, and nutrients from river systems into the ocean, which may strongly affect nearshore environments and deltas (Holmes et al., 2008, 2012; Rawlins et al., 2010). Low flows, which occur during summer and fall in Arctic and subarctic systems, are largely driven by summer rainfall and baseflow conditions. Deeper active layers and subsidence from permafrost thaw may change low flows by increasing the storage of water in the active layer and the movement of water through the soil profile and into rivers (Walvoord and Striegl, 2007; St. Jacques and Sauchyn, 2009; Walvoord and Kurylyk, 2016). Low flow increases in permafrost regions also have important ramifications for exports of nutrients from Arctic systems (McNamara et al., 2008; Frey and McClelland, 2009; Neilson et al., 2018). Finally, streamflow may be changing in the Arctic due to shifting climate and water balance relationships over long timeframes. Rivers in northwestern Canada were observed to have declining magnitudes of annual snowmelt peaks (Cunderlik and Ouarda, 2009). Further, precipitation minus evaporation (P-E) relationships in these basins are likely changing (Serreze et al., 2000; Zhang et al., 2000; Huang et al., 2022).

Challenges in defining streamflow trends and linkages to increased precipitation may be due in part to sparse measurements and cold season biases in precipitation observations (Adam and Lettenmaier, 2003; Yang et al., 2005; Cherry et al., 2007). Indeed, sparse spatial and temporally short observations, low-quality data, and intermittency in the records limit the ability of researchers to conclusively determine the direction and magnitude of changes and to properly attribute these changes (Ahmed et al., 2020). While many studies on changing Arctic streamflow use established methods to infill data (Pigott, 2001; Teegavarapu and Nayak, 2017; Dembélé et al., 2019) and new techniques utilizing machine learning and modeling methods to gap fill records show promise (Zhang and Post, 2018; Arriagada et al., 2021), no current review of infilling in Arctic systems data has been conducted. Thus, determining trends in streamflow records accurately using high quality streamflow gage records along with process-based and statical tools to model streamflow is considered a grand challenge for cold regions hydrological sciences (Blöschl et al., 2019).

Our objective in this work is to consider NA streamflow trends for recent and longer-term periods within gaged basins affected by varying amounts of permafrost coverage. We aim to (1) quantify differences and trends in annual and seasonal streamflow discharges for mean, minimum, and maximum streamflow and (2) assess whether these vary with permafrost coverage across their basins. We look at 74 gaged sites in North America, ranging from 1 to 90% permafrost coverage. We consider the seasonal and annual response in these differences and trends for varying amounts of permafrost coverage and use multiple tools along with newly available datasets and machine learning methods to examine the drivers of change within the basins. Our goal is to highlight changes in streamflow, identify the drivers of shifts within these systems, and to discuss the challenges associated with identifying trends and their drivers in Arctic and subarctic hydrologic systems.

2. Materials and methods

2.1. Study site and streamflow gage selection

We considered streamflow data from across the North American (NA) Arctic, focusing on river basins that contained at least some permanently frozen ground (i.e., permafrost is ground that is frozen for more than 2 years). Thus, our study basins span Alaska and Canada, and range from those draining to the Chukchi Sea in the west to Hudson Bay, Quebec, Canada in the east with latitudes spanning 55° N to 70° N and basin sizes ranging between 27 to 50,200 km² (Figure 1).

FIGURE 1

Figure 1. Streamflow gages (observed, n = 74) by permafrost coverage percentages with infilled gages (n = 55, open circles), and GloFAS v3.1 watersheds (gray circles, n = 958) shown over (Obu et al., 2019) permafrost zones for North America.

To compute basin characteristics needed for gage selection and analysis (e.g., permafrost coverage, topographic slope, aspect, and others) we used the River and Basin Profiler (rabpro, Schwenk et al., 2022). Rabpro automatically delineates watersheds by snapping gage points to appropriate MERIT-Hydro (Yamazaki et al., 2019) flowlines and performing digital-elevation model-based pixel-tracing to identify all pixels draining to the snapped location. After basins were delineated, we used rabpro to compute basin attributes that were necessary to further downselect gages and characterize their basins.

We selected streamflow gages based on several criteria aimed at obtaining a high-quality network of gages that would be reliable for trend analysis. Most importantly, we selected gages affected by frozen ground with mostly complete data records. Gage selection was performed through the Veins of the Earth data platform (https://www.essoar.org/doi/10.1002/essoar.10509913.2). Our initial set of gages (n = 967) included those in the United States Geological Survey (USGS) and Environment and Climate Change Canada (ECCC)'s Hydat databases that contained a minimum of 15 years of continuous data, excluding gages that could not be automatically snapped to flowlines via rabpro. Further downselection was performed using the Reference Hydrometric Basin Network (RHBN 2020 Update, https://www.canada.ca/en/environment-climate-change/services/water-overview/quantity/monitoring/survey/data-products-services/reference-hydrometric-basin-network.html) and GAGES-II (Falcone, 2011) to exclude basins with significant anthropogenic impacts. Daily streamflow time series were downloaded using the tidyhydat (Hydat; Albers, 2017) and dataRetrieval (USGS) R packages. Additionally, real time station data was downloaded manually from ECCC for Canadian stations and merged into a matrix (1975–2022). From this set, we chose stations that had more than 75% of their records complete within the most recent 32-year period, 1989–2021. We then further removed stations that had significant upstream lakes, reservoirs, or glaciers. We also excluded gages located below 55°N whose basins were comprised of <10% permafrost coverage (Obu et al., 2019), and gages located on the same river system and in the same permafrost group (see below) to arrive at the final subset of 74 streamflow gages and basins. This process intentionally selected gages where winter snow was the predominant driver of peak streamflow.

2.2. Data synthesis

We considered the most recent 32-year period (1990–2021, referred to herein as recent-32) and divided the period into two equal length records, 1990–2005, referred to herein as early, and 2006–2021, referred to herein as late, for comparative trend analysis. We also separately considered a longer period (1976–2021, referred to herein as recent-46) with records of 1976–1998 (early), and 1999–2021 (late). Because we downselected stations using the recent-32 period, a few gages were missing larger percentages of data in the recent-46 period. Seasonal streamflow values were calculated from daily data for Winter (November-April), Spring (May-June), Summer (July-August), and Fall (September-October). Annual results are calculated based on the hydrological water year from November-October. We considered minimum, maximum, and mean streamflow for each of the time periods and seasons described above. All data sets associated with this publication can be found in Bennett and Schwenk (2023).

We separated out the stations by percentage of their catchments underlain by permafrost using three groups: 1–25% (non-dominant), 25–50% (moderate), and >50% (dominant) permafrost coverage for our analysis. We did not use the normal classifications of isolated (10% or less), sporadic (10–50%), discontinuous (50–90%), and continuous (90–100% of area) because the watersheds that we examined (particularly the sites located in southern regions) were comprised of a mixture of these zones. Thus, we used rapbro to characterize average present-day permafrost coverage (%) based on the permafrost extent map from Obu et al. (2019) for our study watersheds.

Multiple approaches were used to synthesize our results. We looked at average monthly hydrographs for each permafrost group (km³ yr⁻¹) and discuss the first part of the period (referred to as records, early) compared to the second part of the periods (late). Further, we examined 5-year rolling streamflow means were normalized and centered around three (to remove any confusion around negative value being misinterpreted, unitless). Trends were calculated for periods only using Sen's slope and Mann-Kendall significance tests with the zyp.R software package (Bronaugh and Werner, 2009). Streamflow autocorrelation, or the relationship in the streamflow trends to its lagged values, which can be an issue in trends analysis, was addressed in the trends by following methods of Yue et al. (2002). Streamflow trends are reported in km yr⁻¹ per 30 year (where we calculated the trends per year and then multiply that value by 30, which placed all trends values on the same time interval) and are square root transformed to adjust for outliers for some figures. Trends were considered detectable based on a threshold of the Mann-Kendall p-value of 0.1 for each unit of interest (i.e., season or annual minimum, mean, or maximum streamflow or climate variable), with p-value ≤ 0.05 noted in some cases. These trend analysis and data synthesis methods were also applied to the climate data analysis described below.

To consider impacts of seasonal climate factors on annual streamflow trends, we extracted ERA5-Land data, a global reanalysis data product developed by the European Center for Medium-Range Weather Forecasts (ECMWF) at hourly time steps for surface variables at a resolution ~9 km and is based on the model H-TESSEL (Muñoz-Sabater et al., 2021). The ERA5-Land data used here was extracted using Google Earth Engine (Gorelick et al., 2017) and rabpro. We retrieved monthly rainfall, snow water equivalent, and temperature time series for each basin and calculated seasonal and annual trends.

Finally, average regional streamflow trends across the three permafrost groups for the recent-46 period were estimated at 55 of the sites with the most complete data. GF3.1 was used to infill gage station data and to reduce biases caused by missing data. Average regional annual and seasonal trends were considered for each of the permafrost groups at p-values ≤ 0.1, and p-values ≤ 0.05.

2.3. Machine learning

We used the caret R package (Kuhn, 2008) to run tens of machine learning algorithms to determine the best fit model between annual streamflow trends and a suite of climate and physiographic features, including seasonal mean rainfall, snow water equivalent, and air temperature. We used R² and RMSE to select the best-performing machine learning approach against a test data set and refined the suite of features for further modeling by removing highly correlated features. A single, best performing algorithm (random forest, Breiman, 2001; adjusted R² > 70%) and 10 features were selected, tuned, and used to develop models of maximum, mean, and minimum annual streamflow trends.

2.3.1. Machine learning input data

In addition to trend analysis, we also investigated the probable features that account for observed trends and their changes within permafrost-dominant systems. This analysis was focused only across permafrost-dominant systems. For more robust statistical analyses, we used reanalysis streamflow data from the Global Flood Awareness System (GloFAS) v3.1 (GF3.1, Harrigan et al., 2020), available globally from 1979-present. We pseudo-randomly selected basins from GF3.1 (ensuring a minimum of 1/2° spacing between each sample), resulting in 958 new basins after using the same physiographic down-selection criteria that we applied to the gage data (see Section 2.1). We used these basins as training data for our machine learning model. To ensure the GF3.1 accurately represented trends observed across the permafrost-dominated gaged basins, we compared biases and trends in GF3.1 daily streamflow records with those of the permafrost-dominated gages (see Section 3).

2.3.2. Random forest modeling

Random forests are based on decision trees, which consist of a series of yes/no questions that lead to a predicted value (Breiman, 2001). Random forests can be applied to either classification or regression problems. In our case, we use the random forest for regression to predict a continuous variable (trend). In a random forest, individual decision trees are applied as an ensemble, each of which creates a vote for some value and the average vote across all trees is applied for the final prediction (Liaw and Wiener, 2002). Random forests generally are not impacted by reasonable levels of correlation between input predictor variables and can be developed to avoid over-fitting due to built-in feature randomization and tunable hyperparameters. However, because of the complexity of these models, it can be difficult to understand and control the output of a random forest.

Prior to running the random forest model, we scaled and centered the data and used a grid-search function to optimize the number of variables used in the splits (mtry = 1–10) and the number of trees in the forest (ntree = 50). We applied a 10-fold, repeat (three times) cross-validation, which randomly divided the datasets into 10 non-overlapping parts and then trained on nine of those and tested on one, repeated 10 times to test each subset, and completed this process three times. This approach was used to resample the dataset during hyper parameter selection to avoid overfitting. Ultimately, our final models accurately predicted trends in minimum, mean, and maximum streamflow over both periods (R² > 75%) and was used to determine feature importance (Section 2.3.3).

2.3.3. Machine learning feature importance

Feature importance using caret's varImp function was assessed for these models to determine the most important features for predicting mean, minimum, and maximum annual streamflow trends across the two different periods (recent-32 and recent-46). varImp measures feature importance using the MeanDecreaseAccuracy approach. With this method, a variable's importance is measured by randomizing that variable and then measuring the decrease in out-of-bag accuracy for each individual decision tree in the random forest. Decrease in accuracy is measured as percent change in mean squared error. These decreases in accuracy are then averaged across all trees in the forest to determine the mean decrease in accuracy, which becomes the feature importance of that variable. A larger mean decrease in accuracy means that the variable is more important.

3. Results

Table 1 includes an overview of the different sites utilized in this study, based on rabpro and gage metadata analysis and GF3.1 reanalysis data, and Supplementary Table 1 includes all the observational gages used in this work, ordered by permafrost coverage (%). Figure 1 illustrates the site locations across the NA Arctic/subarctic. Gages in the non-dominant permafrost group (n = 28) were generally at the lowest latitudes (average 59.06° N) with the lowest average drainage areas (26,859 km²) and an average elevation of 781.4 m. Gages in the moderate permafrost group (n = 23) were at moderate latitudes (average 61.51° N), at an average elevation of 884.2 m, and had average drainage areas similar in size to the non-dominant systems (26,493 km²). Gages in the dominant-permafrost group (n = 23) were highest in latitude (average 65.52° N), at an average elevation of 555.6 m, and had the largest drainage basin areas (44,291.3 km²). Average annual temperatures range from −2.6, −4.4, and −7.6°C, while precipitation (rain) values range from 16.8 (9.8), 15.8 (9.4), and 21.1 (12.2) km³ yr⁻¹ in the non-dominant, moderate, and dominant permafrost groups, respectively.

TABLE 1

Table 1. Different streamflow basins sources, number of basins, figures, and analyses used in this study.

3.1. Observed monthly and seasonal streamflow over periods

Average maximum, minimum, and mean seasonal and annual streamflow hydrographs for the recent-32 and recent-46 periods are given in Figure 2 (with values below expressed as the difference between the late and the early record for the recent-32 period. Shaded values in Figure 2 are the average maximum and minimum ranges for monthly streamflow, with the first part of the record (early) showing in gray and the second part of the record (late) shown in pink.

FIGURE 2

Figure 2. Average maximum and minimum ranges for monthly streamflow hydrographs (shaded areas, km³ yr⁻¹), from the non-dominant permafrost (1–25%, n = 28), moderate (25–50%, n = 23), and permafrost-dominant (>50%, n = 23) systems for the early and late records (see legend) in 1990–2021 (left) and 1976–2021 (right). Average mean monthly streamflow is shown with a solid gray (early) and dashed pink (late) lines. Specific months are demarcated with vertical lines, solid black (June), dashed black line (May, permafrost-dominant systems).

The hydrographs show the variation between the late and early records of streamflow averages within each permafrost group. For example, the non-permafrost-dominant systems had, on average, higher annual maximum streamflow during the late record (1.95 km³ yr⁻¹), with little to no change in the peak timing (see top left panel, Figure 2). The moderate permafrost systems, on average, showed peak flow in June in the early records. The late record shows a reduction in the June peak flows, with an overall decrease in maximum spring flows (−2.27 km³ yr⁻¹) and minimum flows (−0.22 km³ yr⁻¹) compared to the early record. Permafrost-dominant systems peak in June in both early and late records but exhibited an earlier and higher May streamflow (+4.77 km³ yr⁻¹) with a lower June maximum (−3.55 km³ yr⁻¹) and minimum (−0.99 km³ yr⁻¹), on average. Minimum streamflow was lower in spring for the < 50 % permafrost groups (recent-32) and moderate permafrost systems (recent-32 and recent-46) in the late compared to the early record, while it was higher for the permafrost-dominant group (0.56 km³ yr⁻¹), and generally exhibited increases in all other seasons and annually.

The standardized and scaled 5-year rolling means for seasonal and annual mean streamflow over the recent-46 years (minus the first 5 years, 1980–2021) are provided in Figure 3, with the recent-32 period indicated by orange dashed lines. For most winter and fall seasons and annual, lower streamflow transitions toward much higher streamflow in the end of the period. Winter and fall mean streamflow patterns indicated that the permafrost-dominated systems were lower earlier in the period compared with other permafrost groups. After 1990, the >50% permafrost systems followed similar trajectories and transitioned to much higher streamflow toward the end of the period. Spring and summer showed coherence across the permafrost groups with substantial interannual variability but no visible patterns. Annual signals showed the permafrost-dominant systems responded differently compared to the other permafrost groups from ~1998 to 2015, and then increased toward the end of the period while the < 50% permafrost groups decreased. Minimum streamflow (Supplementary Figure 1) showed increases across the periods in these systems overall, particularly for winter and fall seasons, and annually. Dominant permafrost systems showed increased flows earlier in time compared to the < 50% permafrost systems. Supplementary Figure 2 shows maximum streamflow rolling means, where strong annual decreases were observed toward the end of the period for all permafrost groups.

FIGURE 3

Figure 3. Five-year rolling means for mean streamflow (scaled and centered around 3, unitless) across permafrost groups for seasons and annually. Time periods for the rolling means are the 1980 to 2021 (42 years). Boundaries for the shorter period (recent-32) are shown with orange dashed lines. Minimum and maximum streamflow are shown in Supplementary Figures 1, 2.

3.2. Streamflow trends

Trends in minimum, mean, and maximum streamflow for individual systems are shown in Figure 4 (with average values for all periods provided in Table 2) for recent-32 and recent-46 periods, using a p-value threshold of 0.1 to identify detectable trends. Overall, a lot of variability exists across the seasons, periods, and permafrost groups. Minimum streamflow trends are generally increasing, with some decreases for a few seasons and periods, such as during spring (non-dominant permafrost group) for the recent-32 (Table 2). The permafrost-dominant group exhibited the strongest positive annual minimum streamflow trend for the recent-46 period (0.25 km³ yr⁻¹ 30 yrs for annual trends, Table 2), with almost half (49%) of the systems exhibiting significant positive seasonal and annual trends (Figure 4). The permafrost-dominant group also exhibited the strongest positive annual mean streamflow trends for the recent-46 period (1.82 km³ yr⁻¹ 30 yrs, Table 2), with ~1/3 (34%) of the systems exhibiting significant seasonal and annual trends (Figure 4). Trends in maximum annual streamflow were mixed, with negative annual trends for most seasons in both periods (Table 2). The strongest decline in maximum streamflow occurs in the permafrost-dominant group during the recent-32 period (−9.51 km³ yr⁻¹ 30 yrs, Table 2). Maximum streamflow had the least number of significant records, with some of the moderate permafrost systems having no significant trends (Figure 4, Table 2).

FIGURE 4

Figure 4. Maximum, minimum, and mean streamflow trends for winter, spring, summer, fall, and annual, recent-32 and recent-46 periods, for all 74 systems, ordered by permafrost % (separate groups shown by a gray line). Negative and positive values are shown for trends and large/small circles show trend significance at a p-value threshold of 0.10.

TABLE 2

Table 2. Permafrost groups, record/period type, seasonal and annual trends for minimum, mean, and maximum streamflow.

The regional streamflow averages across 55 systems were examined for permafrost groups for the recent-46 period in Figure 5 and trend results for minimum, mean, and maximum streamflow are given in Table 3, with significant results at p-value = 0.1 (0.05) denoted in italics (bold italics). Overall, only increasing streamflow trends are significant. Minimum and mean streamflow have more significant trends than maximum streamflow (n = 9, 6, 2, respectively). Minimum and mean streamflow are significantly increasing in winter and fall in almost all permafrost groups (with the exception of the non-dominant mean streamflow in fall). Minimum streamflow trends are significantly increasing in permafrost-dominant systems in all seasons. The largest significant streamflow trend is for maximum streamflow in permafrost-dominant systems in the fall (69.7 km³ yr⁻¹ 43 yrs), followed by mean streamflow in permafrost-dominant systems in the fall (65 km³ yr⁻¹ 43 yrs). Maximum streamflow is increasing significantly only in permafrost-dominant systems in winter and fall. Annually, minimum streamflow trends are significantly increasing for all permafrost groups, while mean streamflow is only increasing in permafrost-dominant systems, and maximum streamflow are not changing significantly.

FIGURE 5

Figure 5. Regional mean annual streamflow averaged across 55 gages. Non-dominant (red) and moderate (green) are not significant, while only dominant permafrost systems are significantly increasing. The non-dominant and dominant permafrost systems have a similar trend (~32 km³ yr⁻¹ for 46 years).

TABLE 3

Table 3. Trends for permafrost groups, seasons, and annual for minimum, mean, and maximum regionalized streamflow (n = 55).

3.3. Drivers of change in streamflow trends

Figure 6 compares observed and GF3.1 mean annual streamflow values for the 23 permafrost-dominant systems. GF3.1 showed strong coherence with the observations of streamflow at these sites (adjusted R² ≥ 0.9) annually and for all seasons (Supplementary Table 2). This finding justified the use of randomly selected permafrost systems from the GF3.1 reanalysis (n = 958) as training data for machine learning models, which we use to consider the drivers of changes in streamflow. Features used in these machine learning models are given in Supplementary Table 3. Out of the 10's of machine learning models we tested, random forest models performed best in the testing period (i.e., maximized R², Supplementary Table 4).

FIGURE 6

Figure 6. Comparison of the observations against the GloFas v 3.1 data sets for annual mean streamflow within the 23 permafrost-dominant systems (km³ yr⁻¹). Orange line are two times the standard deviation.

Ultimately, we trained six random forest models using trends in seasonal ERA5 climate variables and static features as input data (nine inputs total). We trained one model each for predicting trends in annual minimum, maximum, and mean streamflow for the two different periods (recent-32 and recent-46). After training these models, we examined feature importance (Figure 7).

FIGURE 7

Figure 7. Feature importance for random forest predictor variables (Qmin = minimum streamflow, light blue, Qmean = mean streamflow, green, and Qmax = maximum streamflow, dark blue), the recent-32 period (left), and the recent-46 period (right). Climate variables are trends over the period of interest and are derived from ERA5. Features are defined in Supplementary Table 3.

Generally, trends in seasonal climate features such as summer rain and winter and spring snow water equivalence (SWE), and temperature in summer for one model (maximum streamflow, recent-46), rank above static features such as permafrost coverage in feature importance (Figure 7). Additionally, precipitation trends were generally ranked as more important than temperature trends.

We found that features importances changed depending on the random run we considered, an important caveat of this work, therefore we present only the most important variables and general findings herein. For annual minimum streamflow trends, which represent the low flow season, i.e., fall/winter, the top two predictor variables are summer rain and spring SWE for both periods. Permafrost is ranked as the third most important variables in the recent-32 period, while temperature variables are ranked as third, fourth, and fifth most important in the recent-46 period, with other static variables being generally ranked low in importance. For annual mean streamflow, we see again that climate trends are ranked higher than other features, with summer rain, winter SWE (recent-32), and spring SWE (recent-46) ranked as the top variables. Static variables such as permafrost is ranked as fourth, while other static variables are lower in importance. In maximum streamflow trends, spring SWE (recent-32) and rain in summer (recent-46) are the most important variables, followed by winter SWE (recent-32) and temperature in summer (recent-46). This is the only time that temperature variables are ranked in the top two in terms of importance. Permafrost coverage is ranked as third in the recent-32 period, but second to lowest in the recent-46 period.

Overall, permafrost seems to be most important for predicting streamflow trends for the recent-32 period compared to the recent-46 period.

4. Discussion

Seasonal and annual streamflow patterns within permafrost-affected systems shifted in unique ways which can be examined using various methods, as considered in this work. Minimum streamflow increased across all permafrost systems regardless of the percentage of permafrost coverage, particularly for winter seasons, and annually; trends that held whether we considered individual systems or regional patterns. Annual trends in maximum streamflow were decreasing in the recent-46 (for all permafrost-affected systems) and over the recent-32 period (for the >25% permafrost systems). However, for regional trends, the only significant changes (increases) in maximum streamflow was for permafrost-dominant systems during the winter and fall, a new finding that may be reflective of increasing minimum streamflow patterns, owing to increases in storage within deepening active layers, as presented by other authors (St. Jacques and Sauchyn, 2009; Walvoord et al., 2012; Bring et al., 2016; Walvoord and Kurylyk, 2016; Cooper et al., 2023). The changes in annual maximum streamflow, which is indicative of changing peak streamflow, likely reflect a lower snowpack during the spring snowmelt period, and a shift toward an earlier or “flatter” spring streamflow peak, with the non-dominant permafrost systems experiencing no change in timing and a higher peak (Bokhorst et al., 2016; Rixen et al., 2022). As temperatures continue to warm in high latitude regions and peak streamflow moves earlier in the spring season, the streamflow hydrograph may become flatter as observed in our systems (increasing minimum flows and decreasing maximum flows in the spring), leading to a less dramatic peak flow and a flatter hydrograph overall (Stadnyk et al., 2021).

Overall, mean streamflow increased in our permafrost systems and most prominently in permafrost-dominant systems (regional and individual systems), with more significant results compared to maximum streamflow trends across the metrics we considered. This finding is well-documented for permafrost systems across the Arctic (McClelland et al., 2006; Durocher et al., 2019; Ahmed et al., 2020; Shrestha et al., 2021). The strong shift toward increasing flows may be indicative of accelerated streamflow runoff conditions in permafrost dominated systems. Permafrost-dominant systems are very likely changing with regards to their active layers and surface conditions (Makarieva et al., 2019; Wang et al., 2021; Wright et al., 2022), a feature that appears to be already leading to increased mean streamflow now, with stronger increases anticipated in the future (Stadnyk et al., 2021). Interestingly, our overall trend results for permafrost-dominated systems align very closely with recent estimates of increased streamflow for the Mackenzie and Yukon River systems (Druckenmiller et al., 2021; Shiklomanov et al., 2021).

Aside from differences in permafrost extent, the unique changes in permafrost-dominant systems compared to non-dominant permafrost systems could be in part due to differences in elevation and precipitation between permafrost systems. More of the non-dominant permafrost systems are headwater (higher elevation) vs. lowland systems and snowpack change under warming is elevation-dependent (Solander et al., 2018; Rottler et al., 2019), leading to unique changes in streamflow between permafrost systems. Further, permafrost-dominant systems considered here also had higher precipitation compared to the other systems.

The drivers of streamflow changes within our permafrost-dominated systems can be examined using reanalysis tools and models such as ERA5-Land and ERA5-GloFAS datasets in conjunction with machine learning techniques, although these results will be affected by biases within the climatic drivers and modeled streamflow at ungaged basins. As well, our results come with a caveat that the methodology is susceptible to data availability, randomness in the tools we applied, and methodological approach overall. Nevertheless, we saw notable differences in the grouping of climate trends and static features driving streamflow trends in our models. Overall, climate appeared to be a much stronger feature in the streamflow trends examined in these systems compared to the basin characteristics such as permafrost percent, with precipitation being ranked as more important than temperature features. Further, in the different records, basin characteristics did not usually outrank temperature variables, except for a few cases. Only for trends in streamflow in the most recent time period did we see permafrost percent start to become an important variable, which may be indicative of enhanced thawing of permafrost in more recent years, as this process begins to play an increasingly focal role in changing streamflow trends across all components of the streamflow hydrograph.

An important issue in streamflow trend analysis for Arctic and subarctic systems is the high variability that confounds the signal-to-noise within these systems. Part of this is due to the lack of streamflow gages, data intermittency, and lack of long term (i.e., 100+ years) data sets. Intermittency arises from the wintertime freezing of rivers and difficulties associated with maintaining and recording remote streamflow gage data. Longer records were associated with significant increases in statistically significant trends in our data sets; but at the same time, the more recent period showed shifts that may be indicative of enhanced warming and thawing within permafrost basins of the high latitudes. Further, as temperatures rise, precipitation and rain/snow partitioning are shifting, along with changes in the permafrost conditions, all of which can cause increases and decreases in streamflow that obscure changes in the systems. Thus, considering observations alone may be insufficient to clearly evaluate streamflow trends in a statistically robust manner. However, these limitations will likely diminish in the future given the work being done on new datasets, infilling techniques, and improved remote sensing and modeling tools (Muskett and Romanovsky, 2009; Biancamaria et al., 2016; Solander et al., 2017; Harrigan et al., 2020; Feng et al., 2021).

5. Conclusions

Our conclusions for NA streamflow system affected by permafrost are:

• Minimum streamflow significantly increased for permafrost-dominant systems over both periods analyzed, with the strongest annual increases occurring for the longer-term trends over the past 46 years.

• Mean streamflow changed with significant increases in mean streamflow trends in permafrost-dominant systems over the past 46 years.

• Maximum streamflow is changing variably, with the strongest declines in the past 32 years for the permafrost-dominant systems. The only significant longer-term changes occur in the permafrost-dominant systems in fall and winter (increases).

• Statistically significant trends were observed most clearly when regional streamflow trends longer than ~46 years.

• Increases in mean and maximum flows in permafrost-dominant systems are reflective of increasing low flows, deepening active layers, and thawing permafrost, indicating that the entire hydrograph is undergoing change within these streamflow systems.

• Climate drivers of streamflow trends are largely precipitation (rain, SWE), followed by temperature, and then by static drivers such as permafrost coverage.

• Summer rain trends were the most important predictor variable in five out of six of the variables/periods (otherwise ranked fourth), followed by SWE in spring (in the top 2 predictors for 3/6) and SWE in winter (in the top 2 predictors for 2/6).

• Percent of the basin covered by permafrost was not ranked highly (i.e., top 2) overall as a driver of streamflow trends but was the third and fourth most important in streamflow trends for the most recent 32-year period.

• Despite the many challenges to understanding changing streamflow in cold regions, there are new products and datasets in development that are increasingly allowing researchers to better understand the patterns of change in Arctic and subarctic systems affected by varying degrees of permafrost, offering a range of new tools, which, along with continued observational records, may help in improved understanding of changing Arctic streamflow patterns.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers can be found at: https://data.ess-dive.lbl.gov/view/doi:10.15485/1985921.

Author contributions

KB conceived of the work, performed data collection, analysis, original writing, and editing of the paper. JSc and JSt curated all datasets used in analysis except streamflow time series and contributed to writing and editing of the original text. CB provided original writing and edits to the text. EG provided data analysis, figure development, and editing of original text and through the review phases. WB provided feedback on concepts, data consultations, and review of original text and revisions. JR provided feedback on work concepts, review of original text and revisions, and funding as project PI. All authors contributed to the article and approved the submitted version.

Funding

This research was funded as part of the Interdisciplinary Research for Arctic Coastal Environments (InteRFACE) project through the Department of Energy, Office of Science, Biological and Environmental Research Earth and Environment Systems Sciences Division's Regional and Global Model Analysis (RGMA) program. Awarded under contract grant: #89233218CNA000001 to Triad National Security, LLC (Triad). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

Acknowledgments

We thank Alessandro Cattaneo for providing important review and discussions of results during the process of conceptualization, analysis, and writing of this manuscript. We also thank a reviewer for providing valuable and insightful comments on an early iteration of this work.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2023.1099660/full#supplementary-material

Supplementary Figure 1. Five-year rolling means for minimum streamflow across permafrost groups for seasons and annually. Time periods for the rolling means are the 1980 to 2021 (42 years). The shorter record transitions are shown with orange dashed lines. Values are scaled and centered around 3.

Supplementary Figure 2. Five-year rolling means for maximum streamflow across permafrost groups for seasons and annually. Time periods for the rolling means are the 1980 to 2021 (42 years). The shorter record transitions are shown with orange dashed lines. Values are scaled and centered around 3.

References

Adam, J., and Lettenmaier, D. P. (2003). Adjustment of global gridded precipitation for systematic bias. J. Geophys. Res. 108, 4257–4271. doi: 10.1029/2002JD002499