The role of catchment characteristics, discharge, and active- layer thaw in seasonal stream chemistry across 10 permafrost catchments

High-latitude catchments are rapidly warming, leading to altered precipitation regimes, widespread permafrost degradation, and shifts in stream chemistry across major arctic rivers. In headwater catchments, seasonal deepening of flow paths due to active-layer thaw and declining discharge post-snowmelt drive variability in stream chemistry during the open-water period. In North American permafrost regions, activation of deeper mineral layers as the season progresses typically increases major-ion concentrations while decreasing dissolved organic carbon (DOC) concentrations. Despite decades of research on seasonality in stream chemistry within permafrost regions, the relative influence of active-layer thaw and discharge remains unresolved. Additionally, the role of permafrost extent and topography in driving seasonality of these solutes is poorly constrained. To address these knowledge gaps, we measured discharge and sampled major-ion and dissolved organic carbon (DOC) concentrations across 10 permafrost catchments in Yukon Territory, Canada. We analyzed concentration–discharge relationships using generalized additive models to resolve the distinct influence of discharge and seasonal active-layer thaw on stream chemistry and to identify the role of watershed characteristics on the magnitude and seasonality of solute concentrations. After accounting for seasonal variations in discharge, we found that both major ions and DOC were highly seasonal across all catchments, with DOC declining and major-ion concentrations increasing post-freshet. Seasonal variability in major-ion concentrations was primarily driven by active-layer thaw, whereas DOC seasonality was strongly influenced by flushing of soil organic carbon during freshet. While average major-ion concentrations were geologically mediated, greater permafrost extent led to enhanced seasonality in concentrations. Catchments with strong topographical gradients and thinner organic soils had higher specific discharge and lower DOC concentrations but greater relative seasonality. Our results highlight the important role catchment characteristics play in shaping the seasonal variations in and magnitude of solute concentrations in permafrost-underlain watersheds.

1 Introduction

High-latitude catchments are undergoing rapid warming due to polar amplification (Cohen et al., 2014), which is leading to marked shifts in hydrological and biogeochemical cycling that are primarily attributed to thawing permafrost and altered precipitation regimes (McKenzie et al., 2021; Walvoord and Kurylyk, 2016). Permafrost acts as an impeding layer to water movement, restricting flow and separating supra-, intra-, and sub-permafrost aquifers (McKenzie et al., 2021; Woo, 1986, 1980). Permafrost degradation, which includes deepening of the seasonal thawed zone (termed the active layer) and wholesale loss, influences stream chemistry by lengthening and deepening flow pathways, facilitating greater contact with mineral surfaces (Frey and McClelland, 2009; Walvoord and Kurylyk, 2016). This increased contact is expected to increase weathering-associated solute export, although changes in dissolved organic carbon (DOC) export remain uncertain due to potential increases in adsorption and in-stream processing (Frey and McClelland, 2009; Tank et al., 2023). Studies of large arctic rivers report increased exports of organic carbon (Tank et al., 2016), yet results vary across circumpolar regions (Tank et al., 2023), and long-term studies are limited to extremely large rivers that cross multiple eco-zones and permafrost extent, confounding interpretation.

Seasonal active-layer dynamics and catchment wetness strongly influence runoff and solute export in permafrost regions (Carey, 2003; Carey and Woo, 2001; MacLean et al., 1999; McNamara et al., 1998; Woo and Winter, 1993). During spring, large volumes of meltwater infiltrate organic-rich soils and activate near-surface and preferential flow pathways that are rich in dissolved organic matter (DOM; Carey, 2003; MacLean et al., 1999; Shatilla et al., 2023; Woo and Steer, 1986). At this time, water tables are near the surface and distal areas are connected to the stream, mobilizing solutes from shallow soil layers across large areas when flows are typically at their greatest. As the active layer thaws and water in storage declines following snowmelt, flow pathways descend into the underlying mineral soils, and distal areas become disconnected from the stream network, resulting in greater concentrations of weathering solutes and reductions in DOM concentrations, reflecting longer and deeper flow pathways prior to freeze back (Carey, 2003; MacLean et al., 1999). Flow pathways continue to shift in response to additional summer precipitation, resulting in short-term variability in stream chemistry (Koch et al., 2021; Shatilla et al., 2023; Stewart et al., 2022; Zhi and Li, 2020).

The degree of seasonality in stream chemistry is mediated by catchment characteristics such as permafrost extent and topography, where greater permafrost extent leads to a higher yield of DOM and lower export of weathering-associated ions (Petrone et al., 2006; Webster et al., 2022). Topographic gradients can also play an important role in driving DOC concentrations and export, as higher slopes result in thinner organic layers and faster travel times (Jantze et al., 2015; Lee et al., 2019). In cold-mountain catchments, flushing of thin organic layers during snowmelt has been shown to lead to high DOC concentrations, followed by a rapid decrease in concentrations post-freshet (Boyer et al., 1997; Carey, 2003).

Active-layer thaw and catchment wetness are highly seasonal in permafrost catchments, making it difficult to ascertain the relative influence of each factor on driving seasonal solute dynamics. The potential influences of both discharge (as a proxy for catchment wetness) and active-layer dynamics on seasonal stream chemistry in cold regions are rarely disentangled. Shifts in concentration–discharge (CQ) relationships can be used to detect changes in stream chemistry driven by processes other than discharge, such as changes in the internal structure of the catchment (i.e., ground freeze–thaw) or in the quantity and/or mobility of solute stores. For example, Fork et al. (2020) examined changes in the residuals of water flux versus DOC to infer changes in supply in terrestrial pools of DOC. Biagi et al. (2022) and Ross et al. (2022) utilized generalized additive models (GAMs), which are flexible models that allow for the addition of spline terms to model nonlinear behaviour, to quantify changes in CQ relationships over time in agricultural environments. Similar techniques have yet to be applied to permafrost environments to quantify the influence of seasonal drivers on CQ relationships.

In cold catchments, seasonal patterns of ground freeze–thaw, discharge, and solute stores have been identified as drivers of variability in stream chemistry (Boyer et al., 1997; Carey, 2003; MacLean et al., 1999; Petrone et al., 2006). However, the relative importance of each factor and how it is influenced by permafrost extent and other catchment characteristics remain unresolved. To address this, we characterize CQ relationships using discrete major-ion and DOC data across 10 permafrost-influenced catchments in Yukon, Canada. We utilize GAMs with a spline day-of-year (sDOY) term in addition to the traditional CQ approach, along with standard statistical tests, to resolve three hypotheses.

1. Seasonal active-layer thaw activates deeper mineral-rich flow paths, increasing major-ion concentrations and decreasing DOC at all sites, irrespective of seasonal changes in discharge.

2. Permafrost extent is the dominant driver of the magnitude of seasonal changes in concentrations of both major ions and DOC in catchments, where seasonality is higher in catchments with greater permafrost extent (irrespective of seasonality in discharge) due to reduced stream connectivity with sub-permafrost water.

3. Catchments with greater permafrost extent exhibit higher average DOC concentrations and lower major-ion concentrations.

Testing these hypotheses provides new insights into the role of permafrost extent and topography in the seasonality of major-ion and DOC export. As circumpolar regions warm at unprecedented rates, these insights will contribute to a broader understanding of how climate-driven changes in precipitation and permafrost degradation will impact carbon cycling and other geochemical processes.

Figure 1 Map of the Wolf Creek Research Basin (WCRB), highlighting Wolf Creek Outlet (WCO) and its subcatchments: Granger Creek (GC), Buckbrush Creek (BB), and Coal Lake (CL).

2 Data and methods

2.1 Study area

Our study consists of 10 catchments with areas ranging from ∼5 to ∼170 km², all located in Yukon Territory, Canada. Four of the catchments are part of the Wolf Creek Research Basin (WCRB; 60°36^′ N, 134°57^′ W; Fig. 1), with the largest being the WCRB outlet (WCO). Coal Lake (CL) is a ∼ 70 km² subcatchment of the WCRB that contains a ∼1 km² well-mixed lake near its outlet. Granger Creek (GC) and Buckbrush Creek (BB) are alpine headwater subcatchments of the WCRB, each smaller than 10 km². WCO and CL are underlain by sporadic permafrost (Lewkowicz and Ednie, 2004), while high-elevation subcatchments (BB and GC) are underlain by discontinuous permafrost (Lewkowicz and Ednie, 2004).

Low-elevation areas of the WCRB consist primarily of coniferous forests, whereas high elevations have shrub taiga and alpine tundra vegetation (Lewkowicz and Ednie, 2004). Near-surface geology consists of mostly sedimentary rocks such as limestone, sandstone, and siltstone capped by a till mantle and glaciofluvial/glaciolacustrine deposits (Rasouli et al., 2014). The WCRB is part of the traditional territories of the Ta'an Kwach'an Council, Kwanlin Dün, and Carcross/Tagish First Nations. A more detailed description of the WCRB can be found in Rasouli et al. (2014, 2019).

Climate normals (1991–2020) reported near the WCRB at the Whitehorse airport indicate that the average air temperature is 0.2 °C, with an annual precipitation of 279.6 mm (164 mm as rain). Due to the large elevation gradient at the WCRB, average temperatures in the headwater regions are several degrees lower, and precipitation volumes are larger, often exceeding 400 mm yr⁻¹.

Figure 2 Study area map of all Tombstone Waters Observatory (TWO) catchments. Insets show land cover for each catchment.

The other six catchments are part of the Tombstone Waters Observatory (TWO; Fig. 2), which is part of the traditional territory of Tr'ondëk Hwëch'in First Nation. Unlike those of the WCRB, all TWO catchments have areas less than 35 km² and intersect the Dempster highway, a critical north–south road. Named after the kilometre marker where the stream intersects the highway, these catchments are Km 44, Km 71, Km 99, Km 104, Km 175, and Km 185.

Km 44 (64°19^′ N, 138°31^′ W) and Km 71 (64°31^′ N, 138°8^′ W) are the most mountainous catchments, with strong topographic gradients (Table 1). Km 44 primarily consists of bare earth, whereas Km 71 is largely overlain with shrub cover. Km 99 and Km 104 (64°41^′ N, 138°28^′ W) are treeless catchments overlain by a thick layer of peat, except for the high-elevation areas with steep slopes. Km 104 is largely flat, whereas Km 99 has mountainous uplands with more mineral soils. Km 175 (65°12^′ N, 138°26^′ W) and Km 185 (65°16^′ N, 138°18^′ W) are primarily covered by open coniferous forests and are a part of the Ogilvie Mountain Range.

Although there is some discrepancy between various permafrost mapping products regarding permafrost extent, the consensus is that the WCRB sites have less permafrost than the TWO sites (Bonnaventure et al., 2012; Obu et al., 2019; Ran et al., 2022; Yukon Geological Survey, 2023). Surficial geology maps of TWO indicate that Km 71, Km 99, Km 104, and Km 185 are underlain with continuous permafrost (>90 %), whereas Km 44 is underlain with sporadic permafrost and Km 175 is underlain with discontinuous permafrost (Thomas and Rampton, 1982a, b; Yukon Geological Survey, 2023). Climate normals (1991–2020) from Dawson City, the nearest Environment Canada weather station to TWO, report annual precipitation of 319 mm, 208 mm of which is rainfall. The annual air temperature is −3.8 °C. However, local climate can be quite variable amongst the study catchments in TWO.

Bedrock geology across TWO is highly diverse (Colpron, 2022). Km 44 is primarily underlain by thick-bedded quartz arenite, while Km 71 features dark-grey argillite, slate, and phyllite, along with black graptolitic shale and black chert. The geology of Km 99 and Km 104 includes black shale, chert, dolomitic siltstone, calcareous shale, and buff platy limestone in low-lying areas, with interbedded maroon and green slate, siltstone, and sandstone in the uplands. Km 175 and 185 are primarily underlain by grey and buff-weathering dolostone, limestone, black graptolitic shale, minor chert shale, siltstone, and sandstone (Colpron, 2022).

We delineated catchments using high-resolution digital elevation models (DEMs): a 2 m resolution ArcticDEM product for TWO catchments and a 1 m DEM from a 2018 lidar scan for the WCRB catchments. We calculated topographic parameters, including aspect, slope, and mean elevation, using SAGA GIS. We classified land cover via a supervised approach using a resampled DEM and Sentinel-2 multiband imagery with the RStoolbox package in R (Leutner et al., 2023). A summary of catchment characteristics is provided in Table 1.

Table 1 Table of catchment characteristics for the WCRB and TWO catchments. Permafrost extent classification was determined based on a surficial geology dataset (Yukon Geological Survey, 2023). Topographical gradient was classified using the mean catchment slope.

2.2 Stream chemistry

We collected grab samples for major ions and dissolved organic carbon (DOC) across flow states and open-water seasons at all sites over multiple years. Sampling timelines varied by site: from 2017 to 2022 for WCO and GC; from 2018 to 2022 for Km 99, CL, and BB; from 2021 to 2022 for Km 104; and from 2019 to 2022 for all other sites. A summary of sample numbers for each solute and site is provided in Table S1.

Grab samples were filtered through 0.45 µm syringe filters and stored in 60 mL high-density polyethylene (HDPE) bottles (brown HDPE bottles for DOC samples). Samples were kept cool and dark in the field using ice packs and were immediately refrigerated at 4 °C upon return from the field. Specific conductance (SpC) was measured in the field using a multiparameter sonde (YSI ProPlus or ProDSS).

Major-ion concentrations were analyzed at the University of Waterloo Biogeochemistry Lab using ion chromatography (DIONEX ICS 6000, IonPac AS18, and CS12A analytical columns). DOC concentrations were determined at the University of Alberta Biogeochemical Analytical Service Laboratory using a Shimadzu TOC-5000A total organic carbon analyzer, following the US Environmental Protection Agency test method 415.1 (EPA, 1974).

Out of all the ions measured, only calcium (Ca), magnesium (Mg), and sulfate (SO₄) consistently exceeded detection limits and were included in the analysis. For samples with concentrations below detection limits, we assigned a value of half the detection limit. The number of observations varied across sites, ranging from 24 samples for SpC at CL to 92 samples for major ions at WCO.

2.3 Specific discharge

We instrumented all 10 catchments with Solinst leveloggers to determine the continuous volumetric flow rate via stage–discharge relationships. Whenever grab samples for stream chemistry were collected without concurrent flow measurements, we estimated the discharge using the stage–discharge relationships and continuous logger data. We then used the catchment area to determine specific discharge (mm d⁻¹) for easier comparison among catchments. Most catchments experienced low or zero discharge during winter, and the presence of channel ice often introduced significant uncertainty into winter stage–discharge relationships. Consequently, winter flow data were limited. Flow data were generally available during spring freshet, with the exception of BB, where spring flow data are limited due to logistical access constraints during the snowmelt period. We collected flow data from 2017 to 2022 at the WCRB sites and from 2019 to 2022 at the TWO sites. Km 104 was instrumented in late 2021 and only had data available from August 2021 to September 2022.

2.4 Coefficient of variation ratios

We used ratios of the coefficient of variation of logged concentrations (CV_C) and of logged discharge (CV_Q) to quantify chemodynamic behaviour for solutes at all sites. CV ratios are a valuable tool to broadly compare chemodynamic behaviour across catchments and solutes, particularly in situations where discharge is not the primary driver of stream chemistry variability (Musolff et al., 2015; Thompson et al., 2011).

2.5 Concentration–discharge relationships via generalized additive models

The log–log slope of discharge and solute concentrations can be used to characterize whether a solute displays flushing behaviour (positive slope) or dilution (negative slope; Godsey et al., 2009; Hall, 1970). Positive/negative slopes can also indicate the presence of spatial and vertical heterogeneity of solute sources within a catchment (Stewart et al., 2022; Zhi and Li, 2020), whereas slopes close to zero can reflect homogeneity of solute sources or dominance of in-stream biogeochemical processes (Creed et al., 2015).

Generalized additive models (GAMs) extend traditional linear models by incorporating smooth functions that capture nonlinear relationships (Hastie and Tibshirani, 1986). In this study, we used GAMs to enhance traditional CQ models by adding the day of year (DOY) as a spline term, allowing seasonal variability to be explicitly modelled. The modified CQ equation is expressed as

$$\begin{array}{}\text{(1)}& \log \left(C\right)=b\cdot \log \left(Q\right)+\log \left(a\right)+\text{sDOY},\end{array}$$

where sDOY is a smooth term that accounts for seasonal variation, Q is discharge, C is the concentration of the solute, and b and a are constant coefficients. The effect of the sDOY term on the model can be visualized via partial-effect plots. This equation in exponential form becomes

$$\begin{array}{}\text{(2)}& C={\mathrm{10}}^{\mathrm{sDOY}}a{Q}^b.\end{array}$$

The flexibility of the sDOY is user defined to prevent overfitting. GAMs provide the effective degrees of freedom (edf) for the spline term and provide p values for both the sDOY term and the log–log slope (b). The standard error in sDOY varies with the day of the year, often increasing during freshet due to interannual variability in its timing and magnitude. We calculated the mean standard error to estimate the average uncertainty in seasonal effects. To quantify seasonality, we used the range of the sDOY fit. This range is logarithmic, where a value of 1 represents a 10-fold change in concentration after accounting for discharge variability.

We fit a GAM to all 5 reported solutes across all 10 sites using the restricted maximum likelihood (REML) estimation. We used the mgcv package in R for all analyses involving GAMs (Wood, 2000). We reported the b term, R², sDOY range, p values for both the b term and the sDOY, and the mean standard error in sDOY for each model.

The presence of heterogeneity in sources within the catchment and changes in connection to these sources due to ground thaw is indicated by statistical significance in both the log–log slope and the sDOY term. A nonsignificant log–log slope with significant seasonality suggests other processes (i.e., in stream biogeochemical processes) are driving seasonality.

2.6 Influence of topography and permafrost extent on stream chemistry

Accurately measuring permafrost extent remains challenging, so we categorized it as a nominal variable for statistical analyses: continuous, discontinuous, and sporadic permafrost. To evaluate the influence of permafrost extent on solute seasonality (sDOY range), we performed ANOVA and Tukey tests. Given the limited number of sites within each permafrost category, we combined the sDOY ranges for weathering-associated solutes (Mg, Ca, SpC, and SO₄) to assess significant differences.

We did not group DOC with the other ions for ANOVA testing, as it exhibits flushing behaviour. To increase the sample size for statistical testing, we grouped sites with discontinuous and sporadic permafrost into a single category. We then performed t tests comparing the sDOY ranges of DOC between continuous and discontinuous/sporadic permafrost categories. Similarly, we used t tests to examine the relationship between permafrost extent and median concentrations of all solutes, applying the same grouping strategy.

We conducted Pearson's correlation tests to evaluate the relationship between mean catchment slope and both sDOY range and median solute concentrations of all five solutes.

3 Results

3.1 Discharge and concentrations across catchments

To account for sampling gaps at extreme low and high flows, we used the median specific discharge of grab samples for comparisons. The highest median specific discharge was observed at Km 44 and Km 71 (2.2 and 2.03 mm d⁻¹, respectively), while the lowest occurred at Km 99, WCO, and Km 104 (0.57, 0.64, and 0.65 mm d⁻¹, respectively). Notably, Km 44 and Km 71 were characterized by steep topographical gradients, whereas Km 104 had the lowest topographical gradient.

Solute concentrations were generally lower at all WCRB sites compared to the TWO sites (Fig. 3). Km 71 typically had the highest concentrations of major ions, particularly Mg (median concentration of 87.6 mg L⁻¹). Km 175 and Km 185 (underlain by discontinuous and continuous permafrost, respectively) also had relatively high concentrations of major ions in all seasons (Fig. 4). Km 104 and Km 99 had the highest median concentrations of DOC, 15.2 and 8 mg L⁻¹, respectively. In contrast, the lowest DOC concentrations generally occurred at the WCRB sites and at Km 44. Concentrations in spring were typically more variable than in the late season across sites, likely due to high variability in discharge (Fig. 4). Detailed time series information on when samples were collected is represented in Fig. S1.

Figure 3 Solute concentrations plotted against flow for all sites. Winter samples were lumped in with spring. Analogous figures where samples are colour coded by permafrost extent and topographical gradients instead of site names are available in the Supplement (Figs. S2, S3).

Figure 4 Boxplots showing log solute concentrations for spring and summer for all catchments. Winter values were grouped with spring, and grab samples where flow was not available are also included in the boxplots. Samples taken on a DOY between 130 and 280 were classified as summer samples.

3.2 GAMs

3.2.1 Seasonality

GAM results revealed that seasonality, defined by the sDOY term (p<0.05), varied across sites and among solutes (Fig. 5). The sDOY range, indicating the magnitude of relative changes in seasonal concentration changes after accounting for discharge, is detailed in Table S1.

For major ions, sDOY generally increased post-freshet, reflecting rising concentrations as the season progressed. In contrast, DOC exhibited decreasing concentrations over time. Sites within the WCRB, characterized by less permafrost, showed weaker seasonality (lower sDOY range) for major ions compared to the more permafrost-rich TWO sites.

Model performance was generally lower for DOC than for other solutes, particularly for CL (a large lake subcatchment of the WCRB with an outlet lake underlain with sporadic permafrost), Km 175, and Km 104, all of which had an R² less than 0.2. The models exhibited poor performance (R²<0.25) in predicting SO₄ and SpC at CL, as well as in predicting Mg at BB (an alpine headwater subcatchment of the WCRB; Fig. 6). It is important to emphasize that BB had limited winter/spring samples, which may have led to a lower sDOY range.

DOC seasonality was significant at most sites except CL, Km 104, and Km 175. Among sites with significant seasonality, BB and Km 99 exhibited the lowest sDOY ranges (0.23 and 0.33, respectively). Conversely, Km 44 and Km 71, which have steep topography, had the strongest seasonality (sDOY ranges of 0.69 and 1.02, respectively). The standard error for sDOY was generally higher for DOC compared to other solutes.

SO₄ exhibited significant seasonality at all sites. WCO (sporadic permafrost) had the weakest seasonality (sDOY range of 0.18), whereas Km 104 (continuous permafrost) displayed the strongest (sDOY range of 1.78), with high uncertainty (mean standard error of 0.19). Limited data along with some values being below detection limits may explain this high uncertainty. Km 185, also underlain by continuous permafrost, had the second-highest sDOY range (0.81), with a lower standard error (0.05).

The seasonality patterns for Ca, Mg, and SpC were similar. WCO and CL, the largest catchments, were the only sites without significant seasonality of Ca (p values of 0.051 and 0.16, respectively). Mg seasonality was significant only at WCO among the WCRB sites and at most TWO sites except Km 175. SpC seasonality was significant at all sites except CL (p=0.37).

Km 104 and Km 185 showed the strongest seasonality of Ca (sDOY ranges: 0.93 and 0.74), Mg (0.97 and 0.75), and SpC (0.74 and 0.64). In contrast, weaker seasonality occurred at the WCRB sites. CL, GC, and Km 44 had the lowest significant sDOY ranges for Ca (0.17), while WCO exhibited the weakest Mg seasonality (0.19). SpC seasonality was weakest at GC and WCO (0.14 and 0.12, respectively). Complete GAM results are reported in the Supplement.

Figure 5 Partial-effect plot for sDOY term from GAMs for all solutes and sites where the sDOY term was significant (p value <0.05) and mean SE was less than 0.15. The plots represent how the value of the sDOY term in Eqs. (1) and (2) changes depending on the day of year when a sample was taken. Larger sDOY shift over time represents greater relative change in solute concentrations driven by processes other than seasonal discharge. The sDOY value is logged; thus, a shift of 1 (i.e., Mg for Km 104) would result in a 10× increase or decrease in solute concentrations if discharge is constant.

3.2.2 Log–log slopes, sDOY ranges, and CV ratios

For major ions, BB and GC (alpine headwater subcatchments of the WCRB) exhibited the steepest negative slopes, particularly for Mg, with values approaching −0.6 (Fig. 6). Km 104 showed low CV ratios and log–log slopes near zero, indicating chemostatic behaviour for major ions. This is likely due to the presence of a thick organic layer at Km 104, which largely extends to the permafrost table. CV ratios for major ions were generally below 0.7 across all sites, with most sites showing values below 0.3. Km 104 and CL had relatively high CV ratios for SO₄. However, these sites also showed large standard errors in their log–log slopes, reflecting greater uncertainty in their estimates.

For DOC, log–log slopes were either positive or close to zero at most sites. The strongest positive slopes were observed at WCO (0.20) and GC (0.31), indicating high heterogeneity in the soil profile. Conversely, Km 71 (strong topographical gradient), CL (lake-influenced subcatchment), and Km 104 (flat with continuous permafrost) showed negative or near-zero slopes (insignificant b term). Km 104 had a relatively low CV ratio for DOC, whereas Km 71 and CL had higher CV ratios. This indicates that the near-zero log–log slopes for DOC reflect chemodynamic behaviour driven by processes other than discharge at Km 71 and CL. In contrast, DOC is chemostatic at Km 104, likely due to the largely homogeneous soil profile in the catchment.

Figure 6 Log–log slope plotted against 10^sDOY range from the GAMs (a–c). Points are colour coded by site (a), permafrost extent (b), and topographical gradients (c). The 10^sDOY range can be interpreted as the magnitude of change in solute concentrations if discharge remains unchanged. For example, a 10^sDOY range value of 10 (or sDOY range of 1) indicates a 10× increase in concentrations of weathering-associated solutes (decrease for DOC concentrations) if discharge is constant. The sDOY range for SO₄ for Km 104 was large and was not included in these plots for visualization purposes. The CV ratio plotted against log–log slopes from GAMs (d). Gam R² was plotted against sDOY range (e). Error bars for sDOY range represent the mean standard error in sDOY range (mean sDOY SE $\cdot \sqrt{\mathrm{2}}$). Several samples were below the detection limit for SO₄ at Km 104, potentially leading to high uncertainty in the model.

3.3 Influence of topography and permafrost extent on seasonality

ANOVA revealed significant differences in sDOY ranges for weathering-associated solutes based on permafrost extent (Fig. 7). Tukey tests showed that sites underlain by continuous permafrost exhibited significantly greater sDOY ranges compared to those with discontinuous or sporadic permafrost (adjusted p<0.01). The difference between discontinuous and sporadic permafrost was significant (adjusted p=0.061). In contrast to major ions, no significant differences in sDOY ranges for DOC were found between the continuous and discontinuous/sporadic permafrost sites (t test, p>0.1).

Median solute concentrations did not differ significantly between continuous and discontinuous permafrost sites for any solutes. However, DOC concentrations showed a marginally lower p value (p=0.156). Notably, three of the four sites with the highest median DOC concentrations were classified as continuous permafrost, with Km 71 being the only continuous permafrost site exhibiting low DOC concentrations.

Pearson's correlation analysis revealed that sDOY ranges for DOC were positively correlated with mean catchment slope (r=0.59, p=0.07), while all other solutes showed negative but insignificant correlations between sDOY ranges and catchment slope (p>0.1). Median DOC concentrations were negatively correlated with mean catchment slope ($r=-\mathrm{0.62}$, p=0.057). For SO₄, a marginally significant positive correlation (r=0.55, p=0.097) was observed with catchment slope, but no other solutes exhibited significant relationships between median concentrations and slope.

Sites such as Km 71, characterized by both continuous permafrost and strong topographical gradients, can limit statistical significance. High permafrost extent has been shown to lead to higher DOC concentrations, but strong topographical gradients are associated with lower median DOC concentrations.

Figure 7 Boxplots of sDOY ranges for weathering-associated solutes grouped by permafrost extent.

4 Discussion

In permafrost-influenced catchments, solute concentrations are typically highly seasonal (Carey, 2003; Koch et al., 2021; MacLean et al., 1999; Shatilla and Carey, 2019; Shogren et al., 2021; Townsend-Small et al., 2011). Concentrations of DOC generally decrease and concentrations of weathering-derived ions generally increase over the course of the spring-to-fall transition. Seasonal active-layer dynamics, depleting solute stores, and seasonality in discharge have been used to explain these patterns previously (Carey, 2003; Shatilla and Carey, 2019). However, little work has been done to resolve the relative influence of these drivers of variability. Here we provide insights into the mechanisms of solute export in permafrost-underlain catchments by assessing the influence of seasonality on CQ relationships using GAMs.

Results largely confirm our first hypothesis, which states that after accounting for the seasonality in discharge, concentrations of DOC would decrease and the concentrations of major ions would increase post-freshet due to seasonal ground thaw. However, our results also suggest that depletion of soil organic matter during spring and not ground thaw is the primary driver of seasonality in DOC concentrations.

Our second hypothesis states that greater permafrost extent would lead to greater seasonality of major ions and DOC. Although we found evidence to support this hypothesis for major ions, our results indicate that permafrost is less important than other catchment characteristics (i.e., topographical gradients) when considering seasonality of DOC.

Our third hypothesis states that greater permafrost extent would lead to higher median DOC and lower median concentrations of weathering-associated ions. We found some evidence for higher DOC concentrations in catchments with continuous permafrost. However, we found stronger evidence for high topographical gradients leading to lower median DOC concentrations. We found no evidence of permafrost extent influencing median concentrations of weathering-derived ions.

4.1 Seasonal drivers of DOC and major ions

4.1.1 Spatial/vertical heterogeneity in soil chemistry is evident in most catchments

Log–log CQ slopes were predominantly negative (dilution) for all major ions and positive (flushing) or not significantly different than zero for DOC. This corresponds to the work of others in permafrost-underlain catchments, who reported typically positive log–log slopes for DOC (MacLean et al., 1999; Skierszkan et al., 2024) and negative ones for major ions (MacLean et al., 1999; Shatilla et al., 2023). Negative log–log slopes for major ions indicates the presence of more mineral rich soils at depth. Log–log CQ slopes for major ions were closest to zero at Km 104 (a relatively flat catchment with uniform catchment characteristics), indicating high homogeneity in terms of the spatial and vertical distribution of solute sources within the catchment.

A near-zero log–log slope suggests that mobile organic carbon (OC) may be homogeneous vertically/spatially or that in-stream/near-stream biogeochemical processes drive variability in DOC concentrations (Creed et al., 2015; Zhi and Li, 2020). Although many catchments had significant positive log–log slopes for DOC, CL, Km 44, Km 71, Km 104, and Km 175 all had nonsignificant log–log DOC slopes. Km 104 has near-zero log–log slopes with a very low CV ratio. The low log–log slopes likely reflect the largely homogeneous vertical/spatial soil profile due to thick organic soils, a relatively thin active layer (∼40 cm), and weak topographical gradients. Additionally, the stream is slow moving, and in-stream processes may be more important in this catchment than in the other TWO catchments. CL is the second-largest catchment in this study, with a ∼1 km² lake just above the sampling outlet. The nonsignificant log–log slope, along with the high CV ratio at CL, indicates the dominance of in-stream/surface water processes (i.e., photosynthesis, respiration, photo-oxidation), as the high heterogeneity of catchment characteristics in the WCRB rules out homogeneity in vertical and spatial soil characteristics as a potential driver of low log–log slopes. Both Km 104 and CL had insignificant sDOY terms and near-zero log–log slopes, indicating the lack of seasonality in DOC concentrations.

Despite the near-zero log–log slopes, the high seasonality of Km 44 and Km 71 for DOC (Sect. 3.2) suggests that seasonal processes other than discharge and active-layer thaw are the primary drivers of variability in DOC in these catchments. Km 44 and Km 71 have strong topographic gradients, a thin organic layer, and a high proportion of bare ground, particularly upslope of riparian areas. This result suggests that as the source area expands during high flows, water with low DOC concentrations from upslope rocky areas (which do not have a defined organic layer) moves quickly through the organic layer near riparian areas (limiting contact time) due to the steep gradients and dilutes the stream DOC signal.

4.1.2 Seasonal depletion of soil stores is important for DOC export

Our analysis suggests that flushing of organic soils during freshet is likely the primary driver of seasonal declines in DOC, yet active-layer thaw may be a secondary driver depending on the spatial/vertical distribution of organic matter concentrations in the soil. At most sites, high seasonality of DOC (as indicated by the sDOY range), along with nonsignificant log–log slopes, indicates that changing flow paths due to active-layer thaw are likely not the primary driver of the seasonality of DOC concentrations (Fig. 6). The presence of nonsignificant log–log slopes indicates a lack of vertical heterogeneity in soil OC concentrations or the dominance of in-/near-stream processes. Thus, seasonality in CQ relationships for DOC in catchments where log–log slopes are near zero must be primarily driven by processes other than changing flow paths due to active-layer thaw or discharge. Previous work in snow-dominated mountains has shown that the flushing of organic matter in soils during freshet can reduce the soil reservoir of DOC (Boyer et al., 1997; Hornberger et al., 1994). For example, Boyer et al. (1997) attributed declines in DOC as melt progressed to source depletion in the upper soil horizons in a mountain headwater catchment in Colorado, USA. Certain sites (GC, WCO, BB, Km 99, and Km 175) did have positive log–log slopes for DOC, suggesting that the depletion of DOC stores and thawing of the active layer may have a combined effect on controlling DOC concentrations. The strong decline in the sDOY term during freshet at several sites (Fig. 5) indicates that flushing of DOM in soils during freshet rapidly depletes finite labile-OC stores within catchment soils, as has been observed elsewhere (Boyer et al., 1997; Hornberger et al., 1994). Previous work at GC attributed a strong hysteresis of DOC concentrations during freshet and summer events to both a rapid decline in organic matter and greater emphasis on deeper flow pathways (Carey, 2003; Shatilla and Carey, 2019). Townsend-Small et al. (2011) compared DOC concentrations at the same discharge level in the upper Kuparuk River in Alaska and saw a decrease in DOC concentrations over time, suggesting that seasonality is at least partially driven by active-layer thaw and/or depleting solute stores. Skierszkan et al. (2024) used a mixing model and found seasonality in DOC concentrations to primarily be driven by changing flow paths in the Dawson Range, Yukon, Canada, which contrasts with our findings. This discrepancy is potentially a result of strong heterogeneity in spatial and vertical soil OC concentrations in the Dawson Range and a greater presence of organic carbon in soils.

4.1.3 Seasonal active-layer freeze–thaw is important for major-ion export

Our study largely agrees with others on permafrost environments that attributed the seasonality of major ions to seasonal active-layer thaw. Similar to other studies, we found that CQ relationships are variable among seasons for major ions in permafrost-underlain catchments (MacLean et al., 1999; Shatilla et al., 2023). Non-zero log–log slopes can indicate the presence of heterogeneity in soil chemistry, whereas significant sDOY terms indicate seasonal changes in ion concentration after accounting for variability in discharge. The negative log–log slopes and strong seasonality at multiple sites for multiple major ions suggest that the activation of different flow paths due to seasonal active-layer thaw along with variability in discharge drives concentrations of weathering-associated solutes in these regions. Seasonal increases in ion concentrations suggest that the seasonality in CQ patterns is not primarily driven by depletion of ion stores (Fig. 5). Shatilla et al. (2023) attributed ground thaw and the increasing connectivity of deep groundwater to seasonal increases in ion concentrations at GC. Lehn et al. (2017) attributed a seasonal increase in major-ion concentrations in Alaskan watersheds to the seasonally thawing active layer and cryoconcentration in soils during later fall and winter of the previous years. Our results are generally in agreement with these findings, although we did not have the data to assess the influence of cryoconcentration on soils.

Km 104 had significant but relatively low log–log slopes, with the highest seasonality of major ions. This suggests that other drivers (aside from discharge and active-layer thaw) may play an important role in seasonal ion concentrations at this site. Km 104 and Km 99 were found to have highly seasonal transit time distributions, where the fraction of young water is much lower in the summer than in the spring (unpublished data). Greater transit times in the summer can lead to greater contact time between water and mineral soils (Benettin et al., 2015), potentially driving seasonality in this system. However, calculating transit times requires estimates of input tracer concentrations, which was difficult to do in many of our study catchments, particularly during snowmelt.

4.2 Catchment characteristics drive the magnitude of seasonality

4.2.1 Stronger topographical gradients lead to greater DOC seasonality

We found median DOC concentrations to be generally higher at sites with a greater permafrost extent, which supports findings from other permafrost catchments (Frey and McClelland, 2009; MacLean et al., 1999; Petrone et al., 2006). However, we did not find evidence of greater permafrost extent leading to a greater relative change in concentrations of DOC. For example, DOC was highly seasonal at WCO (sporadic permafrost) and had relatively low seasonality at Km 99 (continuous permafrost). Our results contrast with those of MacLean et al. (1999) and Petrone et al. (2006), who observed stronger seasonality of DOC/DOM in catchments with greater permafrost extent. MacLean et al. (1999) reported stronger CQ relationships and greater model performance for DOC in a high-permafrost catchment compared to the low-permafrost catchment in central Alaska during the summer. However, the authors found the CQ relationships to be similar between the catchments during snowmelt. The authors attributed this to a lower variability in discharge during the summer for the low-permafrost catchment. Petrone et al. (2006) observed larger increases in DOC export in catchments with greater permafrost extent. However, Petrone et al. (2006) did not assess the relative influence of seasonal discharge.

At our study sites, catchments with steeper topographic gradients were more important for DOC seasonality than permafrost extent was. Correlation analysis revealed that catchments with stronger topographical gradients, which are often associated with thin organic soils, had higher seasonality in DOC. For example, the two most mountainous sites (Km 44 and Km 71) had the highest seasonality of DOC, indicating the greatest relative seasonal changes in concentrations after accounting for discharge. Both Km 44 and Km 71 had high specific discharge, strong topographical gradients, and the lowest average concentrations of DOC, supporting the idea that organic matter stores were limited in these catchments. Nonsignificant log–log slopes imply that nearly all the seasonality at these sites was due to the depletion of DOC stores and not due to active-layer thaw, as was conceptually outlined in Fig. 8. Although other seasonal processes (water age, stream temperature, etc.) may be important, flushing of soil OC has been observed to be an important mechanism in multiple mountain catchments (Boyer et al., 1997; Shatilla and Carey, 2019). The limited DOC supply and rapid flushing due to high flows presumably led to a more rapid depletion of DOC stores in these catchments.

Conversely, Km 99 and Km 104 had high average DOC concentrations and high permafrost extent yet exhibited lower seasonality than Km 71 and Km 44, which are underlain by sporadic and continuous permafrost but have stronger topographical gradients. Additionally, the CV ratio at Km 104 was also very low, indicating that DOC was chemostatic. Km 99 and Km 104 had much weaker topographical gradients and specific discharge relative to other sites, potentially leading to limited flushing of soils (low specific discharge), which may have played an important role in the weak seasonality of DOC. One caveat is that Km 104 had a very low R² value (<0.2) for DOC and did not have significant seasonality or log–log slopes, potentially due to the lack of multiyear sampling. Km 99 is similar to Km 104 but had a significant positive log–log slope for DOC. Unlike Km 104, Km 99 is not entirely overlain with peat in the active layer, as mineral soils are present in the headwaters. The significant seasonality at Km 99 may be driven by thawing of the active layer, leading to deactivation of shallow DOC-rich flow paths, or due to flushing of the organic layer during freshet. We cannot disentangle the influence of depletion of solute stores and thawing of the active layer on seasonality of DOC in catchments that do have significant log–log slopes. It is likely that in these catchments, both flushing of organic soils and active-layer thaw drive seasonality in DOC concentrations.

The presence of a lake near the outlet of CL can lead to dampening of the terrestrial DOC signal and may be dominated by in-lake biogeochemical processes. This may explain the lack of significance in the sDOY term in CL. Low seasonality at BB relative to GC is likely a function of limited sampling during freshet at BB. Seasonality was also low in Km 175, likely due to high median Fe concentrations (1.8 mg L⁻¹) in the stream (unpublished data), which may lead to the dominance of in-stream processes in Km 175 as Fe interacts strongly with DOM (Gu et al., 1994; McKnight and Duren, 2004).

Figure 8 Seasonality in DOC concentrations is largely influenced by flushing of the organic layer during freshet. This conceptual figure depicts changes in chemistry purely driven by seasonal processes other than discharge. The change in DOC concentrations is less dramatic in low-gradient catchments due to a thicker organic layer. This also results in greater overall DOC concentrations.

4.2.2 Permafrost extent is important for major-ion seasonality

Unlike DOC, permafrost extent and winter groundwater contributions do influence the seasonality of major ions, although the overall concentration and export are more strongly influenced by the regional variability in surficial geology, which is considerable in this study. ANOVA and Tukey tests revealed higher seasonality in major ions for sites underlain with continuous permafrost than for sites underlain with sporadic and discontinuous permafrost. The lower seasonality in catchments with lower permafrost extent can be explained by significant surface–groundwater connection throughout the ice-free season, as outlined conceptually in Fig. 9. Our results are similar to those of Webster et al. (2022), who observed seasonal increases in NO₃ in Alaska, which the authors attributed to a thawing active layer in a high-permafrost-extent catchment but not in low-permafrost catchments. The authors argued that the dominance of groundwater was responsible for the lack of a seasonal trend in the low-permafrost catchments. This supports our hypothesis for major ions, where the presence of deep-groundwater contributions in low-permafrost catchments can dampen seasonality. However, these results should be interpreted with caution, as permafrost extent in mountain headwater catchments is uncertain.

Figure 9 Conceptual diagram depicting the influence of permafrost extent on the relative seasonal change in ion concentrations and the mean concentration of DOC after accounting for seasonal discharge.

5 Conclusion

In this study, we assessed the role of permafrost extent and topography in seasonality of DOC and major-ion concentrations in permafrost-underlain catchments. We obtained new insights into chemical transport through hydrochemical data collected across multiple catchments and seasons, a task rarely undertaken across such a diverse range of remote, permafrost-underlain environments. We utilized GAMs to assess changes in connectivity to solute sources after accounting for seasonally changing discharge. We largely confirmed our first hypothesis, where ion concentrations increased and DOC concentrations decreased (irrespective of seasonality in discharge) after the start of freshet. However, we found evidence of differing drivers of seasonality between DOC and major ions: seasonality in DOC was primarily caused by flushing of OC from soil stores during freshet, while seasonality in major ions was driven by seasonal active-layer thaw. For major ions, our second hypothesis, suggesting that permafrost extent drove seasonality in chemical concentrations, was found to be largely true. However, although permafrost extent may be important for the average concentrations of DOC, topography was a more important driver in terms of seasonality, as mountainous catchments with thin organic soils and high specific discharge led to rapid flushing of limited soil OC stores. We found some evidence for higher DOC concentrations in catchments with continuous permafrost, which supported our third hypothesis. However, we found stronger evidence for high topographical gradients leading to lower median DOC concentrations. We found no evidence of permafrost extent influencing median concentrations of weathering-derived ions, which is presumably due to the diverse geology in the WCRB and TWO.

Our results suggest that climate-change-induced increases in active-layer thickness and greater connectivity of sub-permafrost, intra-permafrost, and suprapermafrost water will lead to elevated concentrations of weathering-derived ions, reduced concentration of DOC, and a decrease in seasonality of major ions in high-latitude catchments. The role of permafrost extent in the seasonality of DOC concentrations remains unclear. This work highlights the need for long-term stream chemistry sampling across a biogeophysical range of high-latitude watersheds. Greater sampling resolution could provide a means of quantifying the influence of secondary drivers (i.e., stream temperature) on stream chemistry in these environments. Currently, much of the literature assesses the role of catchment characteristics in solute exports, but limited work has been done on the drivers of seasonality in cold regions, largely due to the lack of sampling in the shoulder seasons. This study is one of the few to assess the drivers of seasonality of stream chemistry in permafrost-underlain catchments.

6 Future directions

GAMs offer a valuable approach for assessing trends in solute concentrations driven by permafrost thaw, while simultaneously accounting for variability in discharge. Long-term-trend analysis of solute concentrations in permafrost regions is often complicated by concurrent trends in specific discharge, which can lead to significant data omissions (Keller et al., 2010). By leveraging partial-effect plots within a CQ framework, GAMs provide a powerful technique to disentangle the effects of discharge variability from long-term changes in stream chemistry.

Water transit times may be a critical factor influencing stream chemistry, as they govern the contact time between water and soil substrate, while also reflecting flow path length (Benettin et al., 2015; Hrachowitz et al., 2016). This is particularly relevant in permafrost environments, where water transit times can exhibit pronounced seasonality (Piovano et al., 2019). As of yet, the role of seasonal ground freeze–thaw processes in water ages is poorly characterized in permafrost environments. Future work should examine the influence of water transit times on seasonality of stream chemistry in permafrost catchments.