Sensitivity of headwater streamflow to thawing permafrost and vegetation change in a warming Arctic

Climate change has the potential to impact headwater streams in the Arctic by thawing permafrost and subsequently altering hydrologic regimes and vegetation distribution, physiognomy and productivity. Permafrost thaw and increased subsurface flow have been inferred from the chemistry of large rivers, but there is limited empirical evidence of the impacts to headwater streams. Here we demonstrate how changing vegetation cover and soil thaw may alter headwater catchment hydrology using water budgets, stream discharge trends, and chemistry across a gradient of ground temperature in northwestern Alaska. Colder, tundra-dominated catchments shed precipitation through stream discharge, whereas in warmer catchments with greater forest extent, evapotranspiration (ET) and infiltration are substantial fluxes. Forest soils thaw earlier, remain thawed longer, and display seasonal water content declines, consistent with greater ET and infiltration. Streambed infiltration and water chemistry indicate that even minor warming can lead to increased infiltration and subsurface flow. Additional warming, permafrost loss, and vegetation shifts in the Arctic will deliver water back to the atmosphere and to subsurface aquifers in many regions, with the potential to substantially reduce discharge in headwater streams, if not compensated by increasing precipitation. Decreasing discharge in headwater streams will have important implications for aquatic and riparian ecosystems.


Introduction
Headwater streams are low-order streams that are closely connected to terrestrial landscapes and provide a strong control on water quality [1,2], biogeochemical cycling, and biodiversity in broader river systems [3]. Given their small size, limited contributing areas, and close connection to terrestrial landscapes, headwater streams have higher hydrologic variability than larger fluvial systems and are vulnerable to hydrologic extremes and desiccation. Climate change and altered ecosystem processes threaten to alter stream flow in headwater catchments globally due to acceleration of the hydrologic cycle, including higher rates of rainfall, runoff, and evapotranspiration (ET) [4]. In the Arctic, permafrost thaw and changing vegetation physiognomy, cover, and productivity may additionally impact headwater streams. Here we use a space-for-time substitution to reveal the potentially severe impacts of warming and altered hydrology on headwater streams in the Arctic, including increased losses to subsurface flow and lower discharge.
High rates of Arctic warming [5] are leading to broad-scale hydrologic changes [6][7][8] including increased ET related to changes in vegetation composition, density, and distribution [9][10][11][12] and increased infiltration related to permafrost thaw [8,[13][14][15]. In addition to increasing ET, replacement of tundra with shrubs and/or trees alters surface  [56] and treeline indicated. The study catchments are located in the Noatak and Kobuk River basins and consist of varying ground cover [37]. (c) Indicates catchments used for the water budget calculations, streambed infiltration and chemistry, and (d) includes catchments from the broader region where only streambed infiltration and chemistry were determined. energy fluxes and can result in higher ground temperatures [6,16,17] especially in winter [18,19], which could accelerate thaw and increase infiltration potential. Increasing groundwater flow has been inferred from hydrologic and biogeochemical trends in large Arctic rivers [20][21][22], but few studies have provided empirical evidence from headwater catchments [7].
Arctic headwater streams are particularly susceptible to drying due to increasing ET and infiltration because of their size and connection to terrestrial hydrologic processes such as runoff and ET. The majority of Arctic stream discharge occurs during the spring freshet when soils are still frozen, leaving a summer baseflow period dependent on soil moisture and rain to sustain flows [23][24][25]. Numerical models indicate that permafrost thaw may reduce runoff to headwater streams, instead routing runoff through subsurface flow paths to higher order streams and rivers [26,27]. Whether water exits catchments through streams or groundwater has major implications for headwater stream discharge and ecosystems, as well as broad-scale implications for fluxes and cycling of carbon and other solutes and nutrients of importance to greenhouse gas emissions and global biogeochemical cycles [28,29].
To assess the impacts of Arctic warming on headwater stream hydrology, we studied sites at the transition from Boreal to Arctic Alaska, where permafrost is warming and thawing, and vegetation is changing. Warming is expected to increase ground temperatures in this region and reduce near-surface permafrost extent by half by 2090 [30]. Dramatic changes in vegetation have been observed including boreal tree line advance, shrub cover expansion, and increased vegetation productivity [11,[31][32][33]. We hypothesize that warmer headwater catchments will display decreased stream discharge due to the increased ET and thaw-induced infiltration. To test this hypothesis, we monitored hydrologic fluxes in catchments with varying ground temperatures, permafrost content, and vegetation cover [34] (figure 1). We compared these hydrologic fluxes to mean annual ground temperature (MAGT), an integrative metric that is impacted by many relevant factors including permafrost, air temperatures, snow, ecosystem properties, and soils [30]. The range of MAGT of these catchments allowed us to use a space-for-time approach to consider likely changes across these diverse catchments. We developed water budgets for the post-snowmelt period and used stream water chemistry and reach-scale discharge trends [35] to extend our findings to the broader region (table 1).

Site description
Our study was conducted in the Brooks Range of Alaska and included headwater catchments in the Agashashok (n = 5) and Cutler River basins (n = 4) in Noatak National Preserve and in the Akillik River basin (n = 2) in Kobuk Valley National Park ( figure 1). This region encompasses the northern limit of the boreal forest [31] and contains a mixture of Arctic tundra and forested (Picea glauca) slopes, with deciduous shrubs (Salix spp., Alnus spp.) along riparian corridors and alpine regions with bare rock and talus. Valleys contain fine-grained glacio-fluvial sediments overlying coarser cobbles. Mean annual air temperatures at the nearest long-term record in Kotzebue, Alaska, averaged −5.6 • C from 1949 to 2005 [36]. MAGT for our study catchments was determined based on permafrost and ground temperature simulated for the year 2000 by Panda et al [30], which incorporates simulations of near-surface permafrost, ecotype, and soil landscape data. Near surface permafrost within each study catchment was quantified based on mapping by Jorgenson et al [37] and determined by summing ecotype classes with permafrost within 2 m of the ground surface.

Data collection and modeling
Water budgets were determined in 2015 through 2017 for the five catchments in the Agashashok River watershed, calculated based on cumulative fluxes beginning after snowmelt, as determined based on the recession of diel signatures in the stream hydrographs, until stream pressure transducers were removed in August or September prior to freeze up. We modeled the catchment water budget as: where dS/dt is the change in soil moisture over the duration, P is cumulative rainfall, Q is cumulative discharge, and I is the model residual, assumed in most cases to equal infiltration. To account for greater hydrologic connectivity and contributing areas in the early summer when the catchments were mostly frozen [24,38,39], we added rain from the 10 d prior to the end of the snowmelt pulse to the cumulative rain term. Rearranging equation (1) We monitored continuous stream discharge in the five headwater catchments in the Agashashok River watershed between 2015 and 2017 [40,41]. In 2016 through 2018, discharge was also measured in the Akillik and Cutler River basins. Stream stage was measured at the mouth of all streams at 15 min intervals between early June and mid-September using vented pressure transducers (In-Situ Inc., Level Logger 500 or 700 s), or absolute pressure transducers (In-Situ Inc., Rugged Troll 100 s) corrected for barometric pressure. Discharge was measured during site visits in the early and late summer using a wading rod and either a mechanical meter, flowtracker 1, or flowtracker 2 acoustic velocity profiler (Sontek, San Diego, CA, USA). Rating curves were developed between stage and discharge measurements and used to calculate continuous discharge. Field measurements of stage and discharge typically bound all but the highest storm peaks, which were further constrained by channel surveys (see supplementary text T1.2 and supplementary figure S1 for more details).
Soil moisture and temperature were monitored on hourly intervals, usually at 0.03, 0.10, 0.25, and 0.50 m depths in tundra and forested soils using Onset, Inc. (Bourne, MA, USA), H21 dataloggers, S-TMB-M002 temperature sensors, and EC5 soil moisture sensors [42]. Mean freeze and thaw dates and thawed duration were calculated for each depth in forested and tundra sites. Differences between forest and tundra freeze and thaw dates were assessed using analysis of covariance to examine the influence of landscape type (forest or tundra) and sensor depth.
Meteorological data were collected at multiple locations within and near the Agashashok River watershed (table S4), using Campbell Scientific (Logan, UT) stations that were placed in the uppermost study catchment (SFT1) and near the confluence of the Agashashok and Noatak Rivers in late summer of 2015 [43]. These data were used to quantify precipitation and parameters necessary for calculating ET. Precipitation records were adjusted for undercatch and compared to determine variability related to location and elevation, resulting in an overall precipitation uncertainty estimate of 20% (see supplementary text T1.1), approximately half of which is due to spatial variability.
ET for each catchment was quantified using an area-weighted calculation considering four general land cover classes and a combination of simulations based on meteorological data and scaled Table 1. Catchment characteristics and monitoring. For the water budget, chemistry, and seepage runs columns, years indicate the summers in which data were collected using these methods. Mean active layer thickness averaged over the catchment area for the year 2000 determined in [30].

Catchment basin
c Near stream permafrost presented in [57], determined based on permafrost mapping [37].
sap flux data [44]. Because of the small distance between the catchments (tens of kilometers), we assumed that meteorological conditions and subsequently ET were spatially constant for each of the four land covers, such that differences in calculated ET are based solely on the proportion of land cover within each catchment. The four land cover classes are based on mapping by Jorgenson et al [37] and included 'barren ground' , 'tundra' , which included dwarf shrub, shrub/scrub, and sedge/herbaceous categories, 'forest' , which included deciduous, evergreen, and mixed forest categories, and 'open water' . Barren ground occurred predominantly on talus slopes and mountain tops and was assumed to have no ET. Tundra and forest floor ET were calculated using Penman Monteith equation [45]. The open water category consisted predominantly of streams and woody wetlands with a few sparse lakes. Open water ET was calculated with the Penman open water equation [46]. Forest ET was determined by summing sap flux measurements scaled to the forest area and forest floor ET. Sap flux estimates were based on thermal dissipation sap flux sensors [47] (TDP30, Dynamax, Houston, Texas, USA) and data from a 3 m tall micrometeorological station installed in a forest on an eastfacing slope along the main stem of the Agashashok River in late May of 2010. Mean hourly sap flux was analyzed in relation to year, hourly air temperature, atmospheric vapor pressure deficit (VPD), PPFD, soil temperature and volumetric soil water content using Random Forest in the randomForest package [48] of R 3.6.1 (R Core Team 2020). The Random Forest analysis explained 91.4% of the variation in mean hourly sap flux and ranked the importance of predictor variables as follows: VPD, PPFD, air temperature, soil temperature, soil water content and year. Sap flux data were upscaled to the two catchments with substantial forest cover (NFT3 and MST4) using WorldView-2 images from late March of 2017 (see supplementary text T1.4). Stream-catchment exchange was assessed by relating water chemistry (specific conductance (SpC) and deuterium (δ 2 H)) and stream seepage to catchment MAGT. Water chemistry was measured at the mouth of each catchment during site visits, which typically occurred in early June and in August or early September of 2016 through 2019 [49]. Seasonal changes in chemistry were calculated as the difference in early and late-season concentrations. For SpC, this difference was then normalized to the mean concentration of the early and late season samples to aid comparison across streams with widely varying ion loads. Seepage runs [35] were conducted to identify exchange between the stream and subsurface by calculating the difference in discharge at the mouth and at a location 0.3-1 km upstream in each headwater stream (see supplementary text T1.5 for more information).
Stream seepage and water chemistry were individually related to MAGT using linear regression, and subsequently a break point analysis was used to relate δ 2 H to MAGT. For each variable, we tested the assumptions of independent values and normally distributed errors with residual plots. Additional details on data collection, analysis methods, and uncertainty estimates are available in the supplement.

Major differences in catchment water budgets between cold and warm catchments
We observed significant differences in catchment water budgets across an approximately 2 • C range in MAGT during 3 years with varying hydrologic conditions ( figure 2 and table 2). Despite the proximity of the study catchments, the magnitude and temporal patterns of stream discharge varied greatly across catchments (figure 2(b) and supplementary figure  S1). In the coldest catchment (SFT1), discharge was high and consistent among years, averaging 199 mm with low interannual variability (coefficient of variation (CV) = 0.27), whereas in the warmest catchment (NFT3) discharge averaged only 39 mm, with greater interannual variability (CV = 1.21) and some days with no discharge (i.e. <0.01 m 3 s −1 , supplementary table S1). In the post-snowmelt period, discharge was the dominant water loss pathway over the three summers in the coldest catchment, accounting for 72%-98% of rain, whereas in the warmest three catchments there were substantial losses to ET (35%-58% of rain) and infiltration (12%-43% of rain). ET was almost always the larger flux in the warm catchments with consistent infiltration, with a mean ET-infiltration ratio of 1.9 ± 0.8 (supplementary figure S2). The coldest catchment displayed a negative residual in 74% of model runs (averaged over 3 years, supplementary table S2), indicating an unaccounted water source, which may be related to melting of perennial snowfields [50] and subsequent drainage through coarse-grained alpine aquifers [25]. Such a flux is most likely in the highest, coldest catchments (SFT1 and SFT2), which had snow patches that persisted into August. Accounting for snowmelt inputs in the water budget model would result in higher infiltration rates than we present.
Discharge, ET, and infiltration displayed significant relationships to MAGT (figures 2(c)-(e) and supplementary table S3). There was a negative relationship between discharge and MAGT in 2015 and 2016 ( figure 2(c)). ET was positively related to MAGT in all years ( figure 2(d)), which is a result of a greater proportion of catchment cover contributing to ET in the warmer catchments: catchment MAGT was positively related to vegetation cover (tundra plus forest cover, r 2 = 0.97, p < 0.05, supplementary figure S3). Per unit ground surface area, forest ET (tree sap flux plus forest floor ET) was 88% of tundra ET, which is similar to previous findings in Alaska [16] and indicates that here, vegetation type has a relatively small direct effect on the magnitude of ET. Infiltration was directly related to MAGT in 2016 (figure 2(e)), consistent with increased infiltration into warmer, thawed soils.

Greater thaw and infiltration potential in forest soils
Terrestrial ecosystems exert indirect controls on soil thermal regimes and subsequently permafrost [51]. Previous work has identified a 25% increase in absorbed net radiation and local atmospheric heating in forests relative to tundra vegetation in Alaska [6], which can result in greater thaw and theoretically greater groundwater flow potential [52]. During winter, forests maintain deeper and more insulative snowpack than the windswept tundra, preventing the development of deeply frozen soils [19] and contributing to rapid springtime thaw. Our results show that forests are an important location of potential infiltration compared to tundra, because forest soils remained thawed for much of the year and exhibited pronounced summer drying (figure 3). Vegetation (i.e. tundra or forest) significantly determined thaw date (F 1,18 = 11.131, p < 0.001) and freeze date  (2) 36 (2) 45 (5) 81 (9) 88 (10) 86 (10) 22 (24) 44 (22) 92 (27) (4) 27 (2) 94 (5) 70 (8) 76 (8) 75 (9) 68 (23) 54 (20) 44 (27) 1.2  (5) 28 (2) 81 (3) 84 (10) 91 (10) 86 (10) 35 (23) 40 (21) 47 (27) 2 (F 1,18 = 5.408, p < 0.05) after controlling for the effect of depth (see section 2). The greatest difference was in thaw date at 50 cm (figure 3(a)), where forest soils remained thawed on average 63 d longer than tundra soils ( figure 3(a)). Active layer soils in tundra remained frozen well into July, limiting summer infiltration. Seasonal patterns of soil moisture provide additional evidence of enhanced infiltration in forested soils ( figure 3(b)). Both forest and tundra soils were nearly saturated at the beginning of the thaw season. Forest soils at 50 cm displayed a significant decrease in soil moisture over the summer (t-test p < 0.05), resulting in a significant decline in water storage in the top 1 m over the summer (p < 0.001, r 2 = 0.67, supplementary figure S4). The forest soil moisture decrease is only a minor component of catchment water budgets (0%-7.5% of precipitation, table 1), but indicates the potential for water to move through these soils, potentially contributing to ET and/or infiltration. Field and laboratorybased hydraulic conductivity tests on unfrozen active layer soils did not differ significantly between forested and tundra soils (supplementary figure S5), suggesting that disparate trends in soil moisture between the two landscapes are related to differences in soil temperature or freeze/thaw state rather than in intrinsic soil properties.

Increased surface water-groundwater interactions in thawed catchments
Stream-groundwater exchange and stream chemistry from the broader region (figure 4) provide additional evidence that warming and permafrost thaw allows greater infiltration into warmer catchments, and that drastic changes in stream-catchment interactions and chemistry may occur within the 2 • C MAGT range assessed in this study. Stream-groundwater exchange in August was significantly related to MAGT ( figure 4(a), natural log transform of positive-shifted data, r 2 = 0.46, p < 0.05), indicating that warm catchment streams are losing water to the subsurface through their stream beds. Although these measurements are representative of the hydrologic conditions at only one point in time, assuming that this flux persists over 1 km of stream length throughout the months of July and August suggests that streambed infiltration could account for 4% and 12% of the total water budget-calculated infiltration for the two warmest catchments (NFT2 and NFT3, respectively). Streams were steady or slightly gaining water in June (supplementary figure S6), indicating that early season streambed infiltration is limited by frozen soils and/or a high water table.
Specific conductance (SpC) and stable isotopes of water are useful tracers of stream-catchment interactions and the mixing of precipitation with water stored in catchment soils [53], especially in the Arctic due to distinct chemistries of shallow and deep soils and shifting water sources with seasonal thaw [38]. We found a negative linear relationship between MAGT and the seasonal change in stream SpC ( figure 4(b), r 2 = 0.41, p < 0.05) and a weak negative relationship between δ 2 H and MAGT for the colder streams (figure 4(c), r 2 = 0.52, p = 0.11, streams with MAGT −4.4 to −3.4), consistent with a seasonal deepening of subsurface flow depths and a shift in water sources in colder catchments [25,54,55]. The negative relationship between SpC and MAGT likely reflects a gradual change or 'press' disturbance associated with warming and thawing impacting the depths of subsurface flow paths. The break point in the relationship between δ 2 H and MAGT indicates a threshold change that may be related to the flushing of stored snowmelt from cold catchment soils (see further evidence in supplementary figure  S7). In warmer catchments, a direct linear relationship between δ 2 H and MAGT (figure 4(c), r 2 = 0.88, p < 0.05, streams with MAGT −3.4 to −2.5) may indicate isotopic enrichment of the stream water due to greater riparian zone evaporation. Together, water chemistry and seepage trends indicate the potential for greater infiltration and loss of stream water in warmer headwater catchments, suggesting that water resources in Arctic headwater streams may be severely impacted by only a small change in ground temperatures.

Drying Arctic streams and implications
We observed substantial differences in water budgets across a set of catchments in northwestern Alaska that vary broadly in vegetation cover and physiognomy, permafrost extent, and ground ice content, but narrowly in terms of soil thermal regime (<2.0 • C temperature range). Our findings highlight the importance of soil thermal regime on catchment processes, as reflected by the significant correlations between water budget elements, stream chemistry, streamgroundwater exchange, and MAGT. We hypothesize that the observations from our study region are likely generalizable for large regions of the Arctic, given projected northern hemisphere MAGT warming of 2 • C-4 • C by mid-century [56]. Approximately 3 Mkm 2 of permafrost in the northern hemisphere (21% of the total permafrost area) falls within a similar temperature range as our catchments (MAGT from −5 • C to −2.5 • C). Another 5.4 M km 2 of the Arctic (36% of total permafrost area) is underlain by warmer permafrost (MAGT from −2.5 • C to 0 • C) [58]. Ground temperature is the landscape property most strongly indicative of permafrost susceptibility to thaw [59], followed by other physical properties, such as soil texture, ground ice, and topography, that further impact thaw vulnerability and also the hydrological consequences of thaw [13,[59][60][61]. For instance, the capacity for increasing infiltration rates following thaw will be mediated by soil hydraulic properties, as determined by soil texture and phase change. Our study region spans a wide gradient in permafrost conditions (thermal state, areal extent, soil properties) and terrain (vegetation cover, topography; table 1), reflecting the broader mosaic of permafrost and hydrological settings across the Arctic. Our study region shows similar thaw vulnerability as broad regions of Russia and the Tibetan Plateau and much of sub-arctic Canada [59,62], suggesting that the changes documented here may occur on similar timescales across the Arctic. Whereas thaw vulnerability metrics provide an assessment of thaw potential, increasing woody vegetation may be a better (albeit conservative) indicator of where thaw and subsequent hydrological changes are already occurring. Treeline advance tends to lag warming and thaw because of slow reproduction and recruitment of trees at their range limit, yet evidence of treeline advance is common in the Arctic (∼50% of treelines investigated in [63]). In the Agashashok River watershed, elevational treeline has risen at a rate of 0.60 m a −1 since 1952 [11], which is very near the circumarctic mean of 0.48 m a −1 [63]. Additionally, shrub expansion, which is widespread even in colder Arctic zones [9,10,12] has been linked to changing hydrology [64], suggesting that hydrologic changes are occurring even farther north of tree line. Together, similarities in vegetation change and thaw potential in our region and throughout the low Arctic suggest that the hydrological trends observed here might also occur over a vast area in the next several decades.
Multiple lines of evidence presented here suggest that Arctic headwater streams may be substantially smaller or ephemeral in the near future, which may fundamentally alter stream-catchment connectivity and ecosystem function [65]. Linear trends between water budget components and MAGT (figures 2(c)-(e) and table S1) suggest that warming could decrease summer discharge by an average of 69 mm per • C of warming. This is far from a robust prediction, given the complexity of factors and interactions that impact the water budget elements-precipitation trends remain uncertain, increased ET will lag until plants and trees grow, and ET and infiltration may compete for catchment water. SpC trends indicate a gradual change in stream chemistry with warming, presumably as a result of integrating catchmentscale deepening of flow paths that are occurring due to vegetation change and permafrost thaw. Concurrently, stream-groundwater exchange trends indicate an exponential increase in stream bed infiltration with warming, possibly as a result of the development of sub-stream taliks and the connection of groundwater flow paths [26,27,66]. Stream hydrographs indicate that warmer Arctic streams not only have lower discharge, but also higher interannual variability. The lack of variability in colder catchment stream discharge may be related to buffering by melt and drainage of alpine snowfields. With warming, rain may become more prevalent than snow, limiting alpine snowmelt during the summer and leading to levels of interannual variability observed in the warmer streams. These factors are likely to vary widely across Arctic headwater catchments, resulting in highly variable rates of discharge decreases.
There is considerable uncertainty about the magnitude and trajectory of precipitation, ET and infiltration rates under projected warming and thawing. We observed small differences between tundra and forest ET, suggesting that continued tree line advance may not result in a large change in ET. Previous work has found that changing physiognomic classes and increasing woody cover will increase arctic ET from 1% to 13% by 2050 [10]. A larger change in ET may be associated with the encroachment of vegetation into barren ground, given its prevalence (9%-56% of catchment cover) in our study area. Whereas changes in ET may be small to moderate, infiltration may change drastically given that warming and vegetation change will lead to continued thaw and enhanced hydrologic connections among thawed regions. Forest soils were thawed for much longer than tundra soils and displayed more seasonal drainage relative to tundra soils despite similarities in soil physical properties. The impact of trees and tall shrubs on catchment-scale infiltration may be outsized relative to their spatial coverage because woody vegetation tends to exist in the wetter regions such as in the stream riparian zones. The potential for greater infiltration has been recently suggested based on a theoretical modeling study in forested soils [52] and observed in tundra sites where woody vegetation is increasing [65]. Our results highlight enhanced subsurface infiltration and flow in forests and parafluvial zones that are often bordered by forests and tall shrubs.

Conclusion
Permafrost thaw and the opening of new subsurface flow paths has been previously inferred from hydrological and geochemical trends in Arctic rivers [20][21][22] and shown conceptually using models [26,27], but few studies have provided empirical evidence of the catchment-scale processes. We show that infiltration into streambeds and thawed, forested soils in headwater catchments is a substantial source of water to the subsurface and potentially to deeper flow paths that bypass small streams and feed large rivers. These changes can occur even in catchments with extensive permafrost, suggesting that these processes may impact large regions of the warming Arctic. The implications for stream chemistry and ecosystems may be severe. Lower water levels in headwater streams reduce resident and anadromous fish habitat, movement, and connectivity among habitats as observed across ecosystem types [67][68][69]. Increased infiltration and decreased inflows will reduce headwater stream discharge if not compensated by greater precipitation, thereby altering stream temperatures [57] and solute fluxes. Altered stream-catchment interactions related to permafrost presence and thaw may also impact carbon sources and cycling in aquatic ecosystems, food webs, and fish [70]. Our findings highlight the potential for a critical loss of discharge from headwater streams that could result from permafrost thaw and vegetation change in a rapidly warming Arctic.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/https://doi.org/10.5066/P9EIX8ET.