Multi-scale water balance analysis of a thawing boreal peatland complex near the southern permafrost limit in northwestern Canada

Permafrost thaw profoundly changes landscapes in the Arctic-boreal region, affecting ecosystem composition, structure, function and services and their hydrological controls. The water balance provides insights into water movement and distribution within a specific area and thus helps understand how different components of the hydrological cycle interact with each other. However, the water balances of small- (<10¹ km²) and meso-scale basins (10¹–10³ km²) in thawing landscapes remain poorly understood. Here, we conducted an observational study in three small-scale basins (0.1–0.3 km²) of a thawing boreal peatland complex. The three small-scale basins were situated in the headwater portion of Scotty Creek, a meso-scale low-relief basin (drainage area estimated to range between 130–202 km²) near the southern permafrost limit in the Taiga Plains ecozone in northwestern Canada. By measuring water losses (discharge and evapotranspiration [ET]), inputs (rainfall [R] and snow water equivalent [SWE]) and storage change (ΔS), and by calculating runoff (Q), we (1) aimed to quantify the growing season water balances (May–September, 2014–2016) of the three small-scale headwater sub-basins. After (2) comparing monthly sub-basin and corresponding basin water losses through ET and Q, we aimed to (3) assess the long-term (1996–2022) annual basin water balances using publicly available observations of discharge (and thus calculated Q), R and SWE in combination with simulated ET. (1) Growing season water balance residuals (RES) for the sub-basins ranged from −81 to +122 mm. The monthly growing season water balance for the sub-basin for which all water balance components throughout the 3 year study period were recorded exhibited large positive RES for May (+117 to +176 mm), since it included late-winter SWE routinely estimated in late March right before snowmelt. In contrast, lower monthly and negative RES were obtained for June–September (−41 to 0 mm). For two sub-basins, we provide two different drainage area estimates, highlighting the challenges associated with automated terrain analysis using digital elevation models (DEMs) in low-relief landscapes. Drainage areas were similar for one sub-basin, but they exhibited a fivefold difference for the other. This discrepancy was attributed to the high degree of landscape heterogeneity and resulting hydrological connectivity, with implications for Q calculations and RES. (2) Spring freshet contributed 41 % to 100 % (sub-basins) and 50 % to 79 % (basin) of the April–September Q. Spring freshet peaks were comparable, except for the driest year (2014), when the basin Q was more than 10 times lower than in the sub-basins. At both scales, ET was the dominant source of water loss, more than twice Q. (3) Over the long term (1996–2022), the increase in the basin runoff ratio (the ratio of runoff to precipitation) from 1996 to 2012 (0.1 to 0.5) has been attributed to the increased connectivity of wetlands to the drainage network due to permafrost thaw. However, the smaller mean and more variable runoff ratio from 2013 to 2022 may be due to wetland drying and/or changes in precipitation patterns. Overall, we demonstrate how the hydrological responses of rapidly thawing boreal peatland complexes – at both the sub-basin and basin scales – are shaped by complex factors that extend beyond year-to-year changes in precipitation and ET. Long-term hydrological monitoring is crucial to identify and understand potential threshold effects (e.g. changes in land cover and hydrological connectivity) and ecohydrological feedbacks affecting the local (e.g. subsistence activities), regional (e.g. water storage) and global ecosystem services (e.g. carbon storage) provided by thawing boreal peatland complexes.

1 Introduction

A large portion of the Arctic-boreal region is characterised by permafrost (perennially frozen ground). Understanding the interactions between permafrost-thaw-induced landscape changes and hydrological processes is critical for predicting changes in ecosystem composition, structure, function and services in response to climate change (Walvoord and Kurylyk, 2016). Permafrost coverage varies widely across the Arctic-boreal region and increases with latitude and/or altitude (Gruber, 2012). The maximum thickness of the seasonally thawed and hydrologically active layer above the permafrost generally decreases from the southern permafrost limit northwards (Ran et al., 2022). Active layer thickness, which is partly controlled by local climate, ecosystem characteristics and ground properties (e.g. porosity, water content) ranges approximately from more than 1 m (∼60° N) to less than 0.5 m (∼70° N) across Canada (Ran et al., 2022). Higher water content, by simultaneously increasing the latent heat of fusion during thaw and enhancing thermal conductivity, has an opposite effect on active layer thickness. The latent heat of fusion exerts a stronger control on active layer thickness, leading to a thinner active layer (Clayton et al., 2021). For example, in saturated peat deposits with a porosity of about 80 % at 61° N latitude, active layer thickness did not exceed 0.8 m (Connon et al., 2018).

In recent decades, the Arctic-boreal region has experienced a rapid increase in air temperature, up to four times greater than that seen on a global scale (Rantanen et al., 2022). This atmospheric warming has led to accelerated permafrost thaw (Biskaborn et al., 2019; Smith et al., 2022). Additional factors, including natural (e.g. wildfires) and anthropogenic disturbances (e.g. extractive activities; Foster et al., 2022; Klotz et al., 2023), were shown to increase ground heat flux, accelerating permafrost warming and thaw (Gibson et al., 2018; Li et al., 2021). Recent scientific advances have provided insights into the multifaceted and interdependent ecological, hydrological, atmospheric, and biogeochemical consequences of permafrost thaw (e.g. Burd et al., 2018; Carpino et al., 2021; Gordon et al., 2016; Quinton et al., 2019; St. Jacques and Sauchyn, 2009; Torre Jorgenson et al., 2013). In addition, permafrost thaw presents a substantial socio-environmental challenge in the 21st century (Pi et al., 2021; King et al., 2018). For example, accelerated permafrost thaw threatens local communities, infrastructure, and Indigenous livelihoods and cultural practices across the northern circumpolar permafrost region (Gibson et al., 2021; Langer et al., 2023).

From the hydrological and biogeochemical perspectives, permafrost thaw has the potential to cause changes in land cover and hydrological connectivity, and thus in how water and matter moves across and through the changing landscapes of the Arctic-boreal region (Box et al., 2019; Walvoord and Kurylyk, 2016; Wright et al., 2022). For example, thaw-induced changes in land cover and hydrological connectivity can potentially affect the composition and export of both particulate and dissolved organic carbon (Burd et al., 2018; Vonk et al., 2015), mercury methylation (Gordon et al., 2016) and sulfide oxidation and weathering (Kemeny et al., 2023). Additional complexity is added through changes in precipitation regimes, projected to shift from snow- to rainfall-dominated at least in parts of the Arctic-boreal region (He and Pomeroy, 2023; Thackeray et al., 2022). A better hydrological understanding of thawing landscapes in the Arctic-boreal region is crucial to predict the permafrost–carbon feedback strength at the global scale (Ramage et al., 2024; Schuur et al., 2022; Treat et al., 2024).

In the Taiga Plains ecozone of northwestern Canada, permafrost coverage increases from south to north, ranging from isolated (<10 % in areal extent), through sporadic (10 %–<50 %) and discontinuous (50 %–<90 %), to continuous (90 %–100 %) (Ecosystem Classification Group, 2007; Wright et al., 2022). There, a large portion of the low-relief landscape comprises boreal peatland complexes, including permafrost peat plateaus dominated by black spruce (Picea mariana) and permafrost-free, treeless wetlands resulting from surface subsidence due to ground ice melt (i.e. thermokarst; Wright et al., 2022). Such thermokarst wetlands form depressions and receive water from surrounding permafrost peat plateaus. Some thermokarst wetlands are connected to the drainage network and basin outlet through channel fens. Since the 1970s, the faster thaw rate of ground ice-rich permafrost has resulted in the expansion of thermokarst wetlands at the expense of permafrost peat plateaus, especially near the southern permafrost limit in the southern Taiga Plains (Chasmer and Hopkinson, 2017; Wright et al., 2022). There, permafrost thaw was found to be an equal driver of boreal forest loss as wildfire (Helbig et al., 2016a). For example, from 1970 to 2010, permafrost peat plateaus transformed into thermokarst wetlands at rates ranging from 6.9 % to 11.6 % across 10 sites, each covering 10 km² and spanning from 59°58^′ N to 61°18^′ N (Carpino et al., 2018). This prominent thaw-induced land cover change has increased hydrological connectivity across boreal peatland complexes (Connon et al., 2014, 2015; Quinton et al., 2019) and modified the water balances of small- and meso-scale basins, with areas of <10¹ km² and 10¹–10³ km², respectively (Carey et al., 2010; Uhlenbrook et al., 2004).

Understanding the water balances of small- and meso-scale basins is essential for assessing hydrological responses at broader, regional scales (Evenson et al., 2018; Zhang et al., 2018). In the southern Taiga Plains and in other boreal regions in Canada, several studies have focused specifically on evapotranspiration (ET; Helbig et al., 2016b; Isabelle et al., 2018; Warren et al., 2018) or runoff (Q; Connon et al., 2014; Mack et al., 2021; St. Jacques and Sauchyn, 2009). In some studies, Q or water storage changes (ΔS) were obtained as water balance residuals (RES), or ET was estimated using a hydro-chemical method or an empirical equation (Barr et al., 2012; Bolton et al., 2004; Carey et al., 2010; Hayashi et al., 2004). However, studies that investigate the full water balance of small- to meso-scale basins in thawing boreal peatland complexes, with all water balance components measured, are lacking.

Here, we provide a multi-scale water balance analysis using field observations made in three small-scale basins of a thawing boreal peatland complex in the headwater portion of Scotty Creek, a meso-scale, low-relief basin near the southern permafrost limit in the Taiga Plains. The goal was to constrain the headwater sub-basin water balances in a basin context. Specifically, our three objectives were to

1. estimate daily sub-basin water losses (runoff and evapotranspiration), inputs (rainfall and snow water equivalent) and storage change to quantify sub-basin water balances over three growing seasons (May–September, 2014–2016),

2. examine sub-basin hydrological responses in a basin context by comparing monthly sub-basin and corresponding basin water losses through evapotranspiration and runoff, and

3. assess the long-term (1996–2022) annual basin water balance in relation to changes in land cover and hydrological connectivity.

2 Methods

2.1 Study site

Our study site is located within the headwater portion of the Scotty Creek basin (61°18^′ N, 121°18^′ W), which has been estimated to cover 130 km² (this study) to 202 km² (Water Survey of Canada, https://wateroffice.ec.gc.ca, last access: 31 May 2024). The basin lies approximately 50 km south of Fort Simpson, NT, in the sporadic permafrost zone of the southern Taiga Plains (Fig. 1a and b). The continental, subarctic climate of the Fort Simpson region is characterised by long, cold winters and short, dry summers. The climate normals (1981–2010), mean annual air temperature (T_air) and mean annual total precipitation (P), are −2.8 °C and 388 mm, respectively, of which 40 % falls as snow (data from the Fort Simpson Climate station, WMO ID: 71365, were gap-filled with data from the Fort Simpson A station, WMO ID: 71946, Environment and Climate Change Canada, https://climate.weather.gc.ca, last access: 31 May 2024). No significant difference in the snow water equivalent (SWE) between Fort Simpson and observations made in the headwater portion of Scotty Creek were found, suggesting that the Fort Simpson station is a good proxy for SWE at Scotty Creek (Connon et al., 2021). The snow-covered season usually begins in mid- to late October and lasts until mid- to late April or early May. The snow-covered season duration has shortened by 35 d between 1998 and 2014 (Chasmer and Hopkinson, 2017). It is estimated that the permafrost loss rate across the basin has increased from 0.19 % yr⁻¹ (1970–2000) to 0.58 % yr⁻¹ (2000–2015) since the 1970s (Chasmer and Hopkinson, 2017).

Figure 1 (a) Scotty Creek basin location within the southern Taiga Plains ecozone. (b) Study site location within the headwater portion of the Scotty Creek basin. (c) Landscape (i.e. boreal peatland complex) and wetland (i.e. collapsed permafrost-free wetland) eddy covariance towers: 2014 flux footprint climatology (90 % contribution) (Helbig et al., 2016b). Contours of the three small-scale basins, i.e. the West, East, and South sub-basins, derived from automated terrain analysis using a digital elevation model (DEM; West_DEM, East_DEM and South_DEM) and of the East sub-basin derived from field observations (East_FIELD, Connon et al., 2015). The land cover map is from Chasmer et al. (2014). The two outlets, South1 and South2, were located approximately 10 m apart, appearing as a single point.

Underlain by various glacial tills, silts, and clays deposited during the last glacial retreat (Aylesworth and Kettles, 2000), the relatively flat (mean slope: 0.3 %; Quinton et al., 2003) study site is dominated by low-lying peatland ecosystems, which are interspersed with well-drained mineral uplands. The forested mineral uplands are covered by trembling aspen (Populus tremuloides) and white spruce (Picea glauca). The low-lying peatland ecosystems include spatially extensive forested permafrost peat plateaus (“forests”), and permafrost-free thermokarst wetlands (“wetlands”) and lakes (Fig. 1c). Separated from the forests by narrow (a few metres), actively thawing forest–wetland transitions, the topographically lowered (0.5–1 m) wetlands and lakes receive some lateral inflow from the surrounding forests. The wetlands occur mainly as saturated treeless collapse features. Channel fens (a few 10s m in width) connect some of the wetlands to the drainage network and thus route water to the Scotty Creek basin outlet (Quinton et al., 2019; Fig. 1b).

The forest overstorey is dominated by black spruce (Picea mariana) interspersed with tamarack (Larix laricina). The forest understorey and ground cover is dominated by birch shrubs (Betula spp.), bog Labrador tea (Rhododendron groenlandicum), bog rosemary (Andromeda polifolia), reindeer lichen (Cladina spp.), feather moss (Pleurozium schreberi) and Sphagnum spp., respectively (Garon-Labrecque et al., 2016). Abiotic conditions (e.g. soil water content and temperature) change abruptly within a few metres across the transition from “drier and cooler” forests to “wetter and warmer” wetlands (Baltzer et al., 2014; Helbig et al., 2016b, 2017a, 2017b). Wetland vegetation in the collapse mostly includes Sphagnum spp. and ericaceous shrubs such as leatherleaf (Chamaedaphne calyculata) and pod-grass (Scheuchzeria palustris) in the wettest sections. The channel fens are dominated by herbaceous species, including scattered tamarack and glandular birch (Betula glandulosa), abundant seaside arrowgrass (Triglochin maritima) and bog buckbean (Menyanthes trifoliata), and some dense patches of Cyperaceae species. Channel ground cover is dominated by woolly feathermoss (Tomenthypnum nitens) and ribbed bog moss (Aulacomnium palustre).

Peat thickness across the headwater portion of Scotty Creek is generally >3 m, and the mean (± one standard deviation [std]) organic carbon (C) stock was estimated as 167 ± 11 kg C m⁻² (n=3; Pelletier et al., 2017). Forest permafrost thickness is <10 m (McClymont et al., 2013; Quinton et al., 2009), with a maximum active layer thickness in late August/early September of <1 m (Devoie et al., 2021). The mid to late growing season (June to late August/early September) wetland water table position (WTP) usually ranges between 0.1 and 0.2 m below the ground surface, respectively (Helbig et al., 2016b). Table A1 shows a list of all the variables and expressions used in this study, alongside the corresponding abbreviations and acronyms.

2.2 Sub-basin water balance: eddy covariance and supporting measurements

Boreal peatland complex (ET_LAND; 2014–2016) and wetland evapotranspiration (ET_WET; 2014–2016) were obtained from “nested” turbulent energy flux measurements using the eddy covariance technique (Baldocchi, 2014). Identical eddy covariance instrumentation was mounted at the top of a 15 m “landscape flux tower” (AmeriFlux-ID: CA-SCC) and at 1.9 m on a nearby (100 m) 2 m “wetland flux tower” (AmeriFlux-ID: CA-SCB; Fig. 1c). The instrumentation on each tower included a three-dimensional sonic anemometer (CSAT3A; Campbell Scientific Inc., Logan, UT) and an open-path carbon dioxide (CO₂)/water vapour (H₂O_(g)) infrared gas analyser (EC150; Campbell Scientific Inc.) to measure high-frequency fluctuations (10 Hz) in vertical wind velocity and sonic temperature, and CO₂ and H₂O_(g) molar densities, respectively. Due to instrument failure in the landscape flux tower, CO₂ and H₂O_(g) molar densities were measured with an enclosed CO₂/H₂O_(g) infrared gas analyser (LI7200; LI-COR Biosciences Inc., Lincoln, NE) between March and August 2015. Further details on the instrumentation set-up; the calibration and maintenance procedures; the data acquisition, processing and quality control; and the flux footprints calculation for the landscape and wetland flux towers are provided in Helbig et al. (2016b).

Supporting measurements on or near the landscape and wetland flux towers included incoming and outgoing short- and long-wave radiation (CNR4; Kipp & Zonen B.V., Delft, the Netherlands), rainfall (TR-525USW; Texas Instruments Inc., Dallas, TX), T_air and relative humidity (HC2-S3; Rotronic AG, Basserdorf, Switzerland), soil temperature and moisture along vertical profiles, and relative wetland WTP (OTT PLS; OTT Hydromet GmbH, Kempten, Germany; Levelogger Gold F15/M5, Solinst Canada Ltd., Georgetown, ON; HOBO U20 Water Level Data Logger, Onset Computer Corporation, Bourne, MA). Wetland volumetric water content at 5 cm depth was measured with water content reflectometers (CS616; Campbell Scientific Inc.) at a wetland location in each of the three sub-basins. The different low-frequency ancillary data streams were stored as 30 min block averages in an external storage device connected to additional data loggers (CR1000, CR3000; Campbell Scientific Inc.). Forest and wetland SWE were obtained from snow depth (metal ruler) and density measurements (Eastern Snow Conference [30 cm² cross-sectional area] snow tube or snow sampler) along several representative forest and wetland transects during late March (i.e. late winter) snow surveys in 2014–2016 (Connon et al., 2015, 2021).

2.3 Sub-basin boundary delineation

The headwater portion of the Scotty Creek basin was studied based on three small-scale basins (“sub-basins”): West (two outlets, West1 and West2), East (one outlet) and South (two outlets, South1 and South2), together draining approximately 48 % of the landscape flux tower footprint area (Fig. 1c). The wetland flux tower footprint area was located within the West sub-basin. Delineating low-relief basin boundaries and thus drainage areas using automated terrain analysis remains challenging, and estimates tend to vary depending on the level of topographic detail in the digital elevation model (DEM) and the algorithm used (Al-Muqdadi and Merkel, 2011; Datta et al., 2022; Keys and Baade, 2019; Moges et al., 2023). In boreal peatland complexes, differences between “potential” and “effective” drainage areas may arise due to the presence of isolated wetlands disconnected from the drainage network and the basin outlet (Connon et al., 2015). We delineated the boundaries of potential drainage areas for the sub-basin outlets using a LiDAR-derived 1 m DEM using terrain analysis techniques implemented in the ArcGIS Hydrology toolset from the Spatial Analyst toolbox (version 10.2; Environmental Systems Research Institute, 2014; Chasmer et al., 2014). Considering the low-relief landscape, we verified the resulting sub-basin boundaries' plausibility (West_DEM, East_DEM and South_DEM) through visual interpretation of 2010 WorldView-2 imagery (Chasmer et al., 2014). Questionable boundary sections were surveyed using a differential global positioning system (Leica SR530; Leica Geosystems, St. Gallen, Switzerland) in post-processing kinematic mode (centimetre accuracy). Based on a decision-tree land cover classification (Chasmer et al., 2014), West_DEM, East_DEM and South_DEM were dominated by forests (including forest–wetland transitions) and wetlands (combined >95 %). The resulting drainage areas and wetland-to-forest ratios are 0.105 km² (West_DEM), 0.328 km² (East_DEM) and 0.099 km² (South_DEM), and 1.06 (West_DEM), 0.84 (East_DEM) and 1.24 (South_DEM), respectively (Fig. 1c).

Focusing on the hydrological connections between individual wetlands and the sub-basin outlets, the boundaries of the effective drainage areas for the West and East sub-basins were delineated previously (West_FIELD and East_FIELD; Connon et al., 2015). These delineations were based on visual inspection of the same DEM and the 2010 WorldView-2 imagery used in the potential drainage area delineation described in the previous paragraph followed by extensive field observations. Permafrost ridges acting as barriers to water flow, and permafrost-free hydrological connections to channels around and between wetlands and the sub-basin outlets, were identified using a frost probe. All wetlands in the West sub-basin were hydrologically well-connected to the drainage network, resulting in similar drainage area estimates for West_FIELD (0.090 km²) and West_DEM (0.105 km²). In the East sub-basin, several isolated wetlands were not connected to the drainage network, resulting in a fivefold smaller drainage area estimate for East_FIELD (0.068 km²) compared to East_DEM (0.328 km²). We used both drainage area estimates for the East sub-basin, East_DEM and East_FIELD, to calculate sub-basin Q. The South sub-basin contained one individual wetland directly connected to the two outlets (Fig. 1c); thus, we expect the difference between effective and potential drainage area to be negligible (South_FIELD≈South_DEM).

2.4 Sub-basin water balance: discharge measurements

We estimated daily discharge (L d⁻¹) as open water flow at five narrow (1–8 m in width) stream channel locations (= sub-basin outlets) in the vicinity of the landscape and wetland flux towers using rectangular cutthroat flumes (Fig. S1 in the Supplement). The flumes were constructed following open-source design plans (Siddiqui et al., 1996; Skogerboe et al., 1972) and were installed 0.8 m above the channel bottom on wooden damming structures to divert the flow of water through the flumes. At each flume, WTP was measured every 5 min and averaged and recorded every 30 min from April to late August/early September in 2014–2016 using vented pressure transducers (DCX-38 VG; Keller AG, Winterthur, Switzerland). Twelve rating curves to convert WTP to half-hour discharge estimates were obtained from manual discharge and WTP measurements made during and shortly after snowmelt in late April to early May (spring freshet) and late May (baseflow) in 2014–2016, respectively. For the West sub-basin, we used one rating curve per year for each of the two outlets (West1 and West2), thus using six rating curves in total. For the South sub-basin, we used one rating curve per year at the South1 outlet (thus using three rating curves in total) and a single rating curve at the South2 outlet, created in 2015 and used for all 3 years. The East sub-basin consisted of one outlet, which was monitored in 2014 and 2015, with one rating curve per year (no data was available for 2016).

Gaps in the half-hour discharge time series were filled in two steps. First, the half-hour WTP recorded at nearby upstream wetland locations within the respective sub-basin (Haynes et al., 2018) were used to construct monthly and growing season (May–September) proxy rating curves, based on non-gap-filled half-hour discharge for each flume in 2014–2016. At the West sub-basin, 79 % of the discharge data were gap-filled using the wetland WTP method. The mean coefficient of determination, R², (± std) of the monthly linear relationships between the wetland and outlet WTP, calculated for the months with available data between May and September over the 3 years, was 0.70 ± 0.33 (n=8). In contrast, discharge gap-filled using the wetland WTP method accounted for only 3 % and 14 % of the discharge data for the East and South sub-basins, respectively. These monthly rating curves were then used to gap-fill the half-hour discharge time series. Growing season rating curves were used in case of insufficiently strong monthly proxy rating curves. Second, any remaining gaps (0 %, 9 % and 19 % of the data for the West, East and South sub-basins, respectively) due to missing upstream relative wetland WTP were gap-filled using linear regression analysis based on a mean 2014–2016 growing season proxy rating curve. Gap-filled half-hour discharge was summed to obtain daily discharge for the three sub-basins, which was converted to daily sub-basin runoff (Q_WEST, Q_EAST-FIELD, Q_EAST and Q_SOUTH; mm d⁻¹) using the corresponding effective (East_FIELD only) and potential drainage areas (West_DEM, East_DEM and South_DEM).

2.5 Basin water balance: data sets

We obtained several data sets for Scotty Creek spanning 27 hydrological years (October–September 1996–2022). Instantaneous discharge for the Scotty Creek basin outlet (Fig. 1b) along the Liard Highway (61°24^′ N, 121°26^′ W) is publicly available (Scotty Creek at Highway No. 7, 10ED009; Water Survey of Canada, wateroffice.ec.gc.ca). Daily P (mm d⁻¹; R and SWE) are publicly available for the nearest weather station in Fort Simpson (Fort Simpson Climate station, WMO ID: 71365; Fort Simpson A station, WMO ID: 71946, see above). We obtained daily ET (mm d⁻¹) for Scotty Creek (21 hydrological years: October–September 2001–2022) from the Breathing Earth System Simulator (BESS; Jiang and Ryu, 2016), a global biophysical model using seven atmosphere and land products from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument at a spatial resolution of 0.05° (Fig. S2 in the Supplement). We used the mean daily ET for the 2002–2022 period as the daily ET for the 1996–2001 period, i.e. the pre-MODIS era.

We delineated a drainage area for the Scotty Creek basin outlet from the publicly available 90 m DEM of the Shuttle Radar Topography Mission (SRTM, hole-filled SRTM for the globe version 4; Jarvis et al., 2008) using automated terrain analysis implemented in the ArcGIS Hydrology toolset from the Spatial Analyst toolbox (Environmental Systems Research Institute [esri], 2014). The terrain analysis derived potential drainage area was 130 km², which is smaller than previously published drainage area estimates for the Scotty Creek basin outlet: 134 km² (Burd et al., 2018), 139 km² (Chasmer and Hopkinson, 2017), 150 km² (Quinton et al., 2004), 152 km² (Connon et al., 2014) and 202 km² (Water Survey of Canada). For reproducibility and methodological consistency with the sub-basin drainage areas, the BESS model estimates of ET were averaged across Scotty Creek using the terrain derived drainage area (130 km², this study). All data sets were temporally aggregated to monthly and annual (hydrological year: October–September) runoff (Q_BASIN), precipitation (P_BASIN), SWE and rainfall (SWE_BASIN, R_BASIN) and evapotranspiration (ET_BASIN). We used the lower (130 km² (this study) and upper basin drainage area estimates (202 km², Water Survey of Canada) to calculate Q_BASIN (Q_{BASIN_130} and Q_{BASIN_202}, respectively).

2.6 Multi-scale water balance analysis

We calculated monthly (mm month⁻¹; West sub-basin), growing season (mm per growing season; West, East and South sub-basins denoted as “WEST”, “EAST” and “SOUTH” subscripts) and annual (hydrological year: October–September, mm yr⁻¹; Scotty Creek basin denoted as “BASIN” subscript) water balances as:

$$\begin{array}{}\text{(1)}& R+\text{SWE}=\text{ET}+Q+\mathrm{\Delta }S\end{array}$$

where ΔS is water storage change, rainfall (R) plus snow water equivalent (SWE) is total precipitation (P), and ET is evapotranspiration. Groundwater discharge from permafrost thaw was expected to be negligible (Connon et al., 2014; Quinton et al., 2019).

For simplicity, we loosely defined the growing season as the May–September period when actual measurements for all water balance components (Eq. 1) for the complete months were available. For example, the wetland WTP measurements start in May, because before then, the wells were frozen. Water table position was used to calculate ΔS_SUB-BASIN, as we assumed that Q occurs from forests to the topographically lower wetlands (Wright et al., 2022). Therefore, we calculated ΔS_SUB-BASIN for the West, East, and South sub-basins from the wetland ΔS, scaled by the sub-basin specific wetland area coverage (A_WET). ΔS_WET was calculated based on saturated and unsaturated peat layers using WTP variation, volumetric water content at 5 cm depth and peat porosity values at 3 cm (=0.92) and 15 cm (=0.86) from Isabelle et al. (2018).

Precipitation (P_SUB-BASIN), including R and SWE in late March just before the start of snowmelt (i.e. SWE_MAX), was obtained from rain gauge measurements ($R_{\text{WEST}}=R_{\text{EAST}}=R_{\text{SOUTH}}$) and calculated as a weighted mean for each sub-basin (SWE_{MAX_SUB-BASIN}) according to sub-basin-specific cover areas (i.e. wetland [A_WET] and forest areal coverage [A_FOR]) and associated measured SWE from late-winter snow surveys (i.e. forest [SWE_{MAX_FOR}] and wetland SWE [SWE_{MAX_WET}]), respectively:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\text{SWE}}_{\text{MAX\_SUB-BASIN}}=\\ & {\displaystyle \frac{A_{\text{FOR}}\times {\text{SWE}}_{\text{MAX\_FOR}}+A_{\text{WET}}\times {\text{SWE}}_{\text{MAX\_WET}}}{A_{\text{SUB-BASIN}}}}\end{array}\end{array}$$

where A_SUB-BASIN denotes the sub-basin area. We added SWE_SUB-BASIN to R in May, as we assumed that the main contribution of snow to the P_SUB-BASIN, and thus to the growing season and annual water balances, occurred primarily through complete snowpack melting.

Mean energy balance closure fractions at the landscape and wetland flux towers were 0.70 (0.67, 0.72 and 0.72 from 2014 to 2016) and 0.67 (0.65, 0.69 and 0.68 from 2014 to 2016), respectively. To account for sensible (H; W m⁻²) and latent heat underestimation (LE; W m⁻²), we applied the closure fraction correction by preserving the Bowen ratio (H LE⁻¹) to obtain the corrected LE (i.e. ET) (Barr et al., 2012; Isabelle et al., 2020). The closure fraction correction was calculated using 30 min average fluxes for the months of July–September, when the most complete energy flux data were available. Mean growing season forest and wetland flux footprint area contributions to ET_LAND (corresponding to ET_WEST) measured at the landscape flux tower were approximately 50 % each (Helbig et al., 2017a, 2016b; Warren et al., 2018; Fig. 1c). In contrast, the mean growing season footprint for ET_WET consisted solely of wetland surrounding the tower (Helbig et al., 2016b; Warren et al., 2018). For the South and East sub-basins, we calculated forest ET (ET_FOR, Eq. 3) using ET_LAND and ET_WET as:

$$\begin{array}{}\text{(3)}& {\text{ET}}_{\text{FOR}}={\displaystyle \frac{\left({\text{ET}}_{\text{LAND}}-\frac{A_{\text{WET}}}{A_{\text{SUB-BASIN}}}\times {\text{ET}}_{\text{WET}}\right)}{\left(\frac{A_{\text{FOR}}}{A_{\text{SUB-BASIN}}}\right)}}\end{array}$$

Evapotranspiration for the South and East sub-basins was calculated as weighted means as for SWE_SUB-BASIN (Eq. 2). Sub-basin runoff (Q_SUB-BASIN) was obtained from daily discharge measurements and the corresponding sub-basin areas.

Annual basin water balances (mm yr⁻¹, Eq. 1) were calculated using temporally aggregated precipitation (P_BASIN) and rain (R_BASIN) measurements from Fort Simpson (Fort Simpson Climate station, WMO ID: 71365; Fort Simpson A station, WMO ID: 71946), with snow water equivalent (SWE_BASIN) simply calculated as P_BASIN minus R_BASIN, and ET estimates from the BESS model (ET_{BESS_BASIN}). The ΔS_BASIN was calculated as the difference between the water inputs (P_BASIN) and outputs (Q_BASIN and ET_{BESS_BASIN}) of Eq. (1). A positive value indicated an increase in the water stored in the basin, and vice versa.

We compared growing season monthly Q_SUB-BASIN and ET_SUB-BASIN, both calculated as the means of the corresponding West, East and South sub-basin estimates, with Q_BASIN and ET_{BESS_BASIN}, respectively, using ordinary least squares (OLS) regression analysis. The Q_BASIN and Q_SUB-BASIN used for this comparison were obtained from the drainage area derived from the automated terrain analysis of a DEM in this study. Similarly, we compared monthly ET_LAND with headwater ET estimates from the BESS model (ET_{BESS_HEAD}) using OLS regression analysis. The OLS regressions uncertainty was estimated using bootstrapping with 1000 iterations. The ET headwater estimates from the BESS model were calculated as the mean of four pixels, i.e. the pixel containing the landscape flux tower and three adjacent pixels representing the Scotty Creek headwater portion (Fig. S2). We examined the annual (hydrological year: October–September) hydrological balance components, i.e. Q_BASIN, P_BASIN, R_BASIN and SWE_BASIN time series (1996–2022), and calculated the annual ratio of runoff to precipitation (the runoff ratio).

3 Results

3.1 Meteorological conditions

The annual mean T_air of the Fort Simpson region over the 3 year study period fell in the range of (2014) the 27 year mean (1996–2022, Table 1), or was higher (2015, 2016) than it. The first year of the 3 year study period (2014) was much drier, with less snow and rainfall compared to both the other 2 years and the 27 year study period. The annual total P in 2016 and 2015 was lower and higher than the 27 year mean, respectively, but within one std. The start and end of the snow cover period were consistent throughout the 3 year study period.

Table 1 Annual mean air temperature (T_air), total precipitation (P), snow water equivalent (SWE) and rain (R) at Fort Simpson (data from the Fort Simpson climate station, WMO ID: 71365, were gap-filled with data from the Fort Simpson A station, WMO ID: 71946, Environment and Climate Change Canada, https://climate.weather.gc.ca, last access: 31 May 2024), dates of snowmelt end and start of a spatially continuous snow cover and snow-free season length at Scotty Creek.

3.2 Sub-basin growing season water balances

The hydrographs of the West, East and South sub-basins were dominated by the spring freshet, caused by the rapid melting of the snowpack starting in late April (Fig. 2a–c). Each year, the peak in Q_SUB-BASIN occurred within 2–4 d after the start of snowmelt. For each sub-basin, the spring freshet (April–May) Q was the lowest in 2014 (15, 44, 27 and 130 mm for the West, South, East_DEM and East_FIELD sub-basins, respectively) and the highest in 2016 (83 and 104 mm for the West and South sub-basins, respectively), with intermediate values in 2015 (54 and 77 mm for the West and South sub-basins, respectively). A maximum in daily Q of 12 mm d⁻¹ was observed in the South sub-basin (highest wetland-to-forest ratio) in 2016, coinciding with a heavy rainfall event (>30 mm d⁻¹) 10 d before the start of snowmelt (Fig. 2c). The spring freshet accounted for 99 % and 100 %, 73 % and 87 %, and 83 % and 89 % of Q over the April–September period in 2014, 2015 and 2016, for the West and South sub-basins, respectively. In contrast, the spring freshet for the East sub-basin accounted for 41 % and 47 % of Q over the April–September period in 2014 and 2015, respectively. Once the spring freshet ceased, only the East sub-basin sustained continuous Q throughout the remainder of the growing season (baseflow) in 2014 (drier than normal conditions; Fig. 2a). All three sub-basins sustained a continuous Q after the spring freshet in 2015 (wetter than normal conditions) but not in 2016 (drier than normal conditions; data only for West and South sub-basins in 2016). All post-spring freshet variations in Q were in response to individual rainfall events, reaching amounts of up to 30 mm d⁻¹.

Figure 2 Basin and sub-basin hydrographs in (a) 2014, (b) 2015 and (c) 2016 – daily rainfall ($R_{\text{EAST}}=R_{\text{SOUTH}}=R_{\text{WEST}}$, mm d⁻¹), boreal peatland complex evapotranspiration (ET_LAND) approximately corresponding to ET from the West sub-basin (ET_LAND≈ET_WEST, mm d⁻¹), runoff (Q, mm d⁻¹) from the Scotty Creek basin, and Q (mm d⁻¹) from the East, South, and West sub-basins approximately draining the landscape flux tower footprint area (Fig. 1c). East_DEM and East_FIELD drainage areas are used to compute the lower and upper Q range contours (DOY = day-of-year).

Over the 3 year study period, mean daily ET_LAND was 2.9 ± 1.1 mm d⁻¹ (ranging from 0.6 to 5.5 mm d⁻¹) and ET_WET was 3.3 ± 1.5 mm d⁻¹ (ranging from 0.4 to 8.1 mm d⁻¹). The daily ET of the boreal peatland complex (ET_LAND≈ET_WEST) increased continuously from 0.3 mm d⁻¹ in early April to 2.5 mm d⁻¹ in late May, coinciding with the rapid melting of the snowpack. From late May until late September, daily ET ranged between 2.0 and 4.0 mm d⁻¹ for most of the time (Fig. 2a–c). At 366 mm, total ET in April–September was the lowest in 2014 (mean T_air was 11.1 °C). In contrast, total ET and mean T_air from April–September were similar in 2015 and 2016 (447 and 458 mm, and 11.5 and 11.6 °C, respectively). Comparatively, total Q_WEST was 15, 75 and 101 mm for the April–September period in 2014, 2015 and 2016, respectively. Thus, total ET_WEST was approximately 24, 6 and 5 times greater than total Q_WEST in 2014, 2015 and 2016, respectively.

Differences in growing season (May–September) water input (P_SUB-BASIN) and combined losses (ET_SUB-BASIN and Q_SUB-BASIN) ranged between −211 mm (net loss: 2016, South sub-basin) and +21 mm (net gain: 2015, West sub-basin), resulting in ΔS_SUB-BASIN of similar magnitudes (ranging from −250 mm [2016, South] to +3 mm [2015, East]) across sub-basins and years (Fig. 3a–c; Table S1 in the Supplement). However, the difference between water input and combined losses for the East_FIELD sub-basin was −354 and −311 mm, in 2014 and 2015, respectively (Fig. 3b).

Considering the variations in ΔS_SUB-BASIN, growing season water balance residuals, RES_WEST and RES_SOUTH (Eq. 1), were positive for the West (+114, +122 and +34 mm in 2014, 2015 and 2016, respectively) and South sub-basin (+38 mm in 2016) (Fig. 3, Table S1). In contrast, growing season water balance residuals for the East sub-basin, RES_EAST and RES_EAST-FIELD, were negative in 2014 (−81 and −287 mm, respectively) and 2015 (−30 and −285 mm, respectively). In the West sub-basin, all the water balance components were recorded over the 3 year study period, enabling us to calculate the monthly water balance during the growing season.

Figure 3 Growing season (May–September, 2014–2016) water balances (mm per growing season) for the (a) West, (b) East and (c) South sub-basins: rainfall ($R_{\text{EAST}}=R_{\text{SOUTH}}=R_{\text{WEST}}$), snow water equivalent (SWE_EAST, SWE_SOUTH and SWE_WEST), evapotranspiration (ET_EAST, ET_SOUTH and ET_LAND≈ET_WEST), runoff derived from the terrain analysis drainage area (Q_EAST, Q_SOUTH and Q_WEST) and water storage change (ΔS_EAST, ΔS_SOUTH and ΔS_WEST). The black dot symbol indicates the water balance residual (RES_EAST, RES_SOUTH and RES_WEST) resulting from Eq. (1). (b) For the East sub-basin, ET_EAST-FIELD, Q_EAST-FIELD and ΔS_EAST-FIELD are estimated from the effective drainage area derived from field observations (EAST_FIELD, Connon et al., 2015). SWE_EAST-FIELD is comparable to SWE_EAST. The white dot indicates RES_EAST-FIELD. Wetland water table position and discharge data to calculate ΔS_SOUTH and Q_EAST are not available for 2014 and 2015 (not measured), and 2016 (instrument failure), respectively.

3.3 Sub-basin monthly growing season water balance – West sub-basin

The negative ΔS_WEST in May indicates a large reduction in the water stored in the West sub-basin, even though the total water input (R_WEST plus SWE_WEST) exceeded losses by 20 % (2016) to 50 % (2014 and 2015) (ET_WEST plus Q_WEST, Fig. 4a–c, Table S2 in the Supplement). This discrepancy is reflected in the large positive monthly water balance residuals (RES_WEST) in May each year (149, 176 and 117 mm in 2014, 2015 and 2016, respectively), reaching almost twice the magnitude of ΔS_WEST in 2014 and 2015 (Fig. 4a and b). In contrast, monthly RES_WEST from June–September in all 3 years were an order of magnitude lower than those in May (from −41 to 0 mm with a mean of −14 mm, Table S3 in the Supplement). In the 3 year study period, ET_WEST was similar during the early- to mid-growing season (June–August: mean monthly ${\text{ET}}_{\text{WEST}}\pm \text{one std}=\mathrm{95}$ ± 9 mm). Mean monthly ET_WEST during the late growing season (September) was 45 ± 8 mm. For the June–September period, the 2014 total R_WEST (188 mm) was lower than the total ET_WEST (291 mm), and ΔS_WEST was −69 mm. Similarly, in 2016, ET_WEST (361 mm) largely exceeded R_WEST (185 mm), and ΔS_WEST was −110 mm. In contrast, during the June–September period in 2015, R_WEST (291 mm) was closer to ET_WEST (336 mm), and ΔS_WEST was −10 mm.

Figure 4 Growing season monthly (May–September, 2014–2016) water balances (mm month⁻¹) for the West sub-basin: rainfall (R_WEST), snow water equivalent (SWE_WEST), evapotranspiration (ET_LAND) approximately corresponding to evapotranspiration from the West sub-basin (ET_LAND≈ET_WEST), runoff (Q_WEST) and water storage change (ΔS_WEST). The black dot symbol indicates the monthly water balance residual (RES_WEST) resulting from Eq. (1).

3.4 Comparison between sub-basin and basin evapotranspiration and runoff

Comparable spring freshet peaks were observed between the basin and sub-basins, except for the driest year (2014), when Q in the basin hydrograph (<0.6 mm) was substantially lower than in the sub-basin hydrographs (from 1.6 to 9.4 mm, Fig. 2). At the basin scale, the spring freshet contributions (April–May) to Q varied between 50 % and 79 % over the April–September period in 2014 to 2016, i.e. in the range observed for the three sub-basins (from 41 % to 100 %). Monthly Q between the sub-basins (using the drainage area obtained with terrain analysis techniques) and the basin were comparable (Fig. 5a). The greatest absolute difference was twofold in May (from 1.6 to 2.3; Fig. 5a). Total ET from the BESS model over the April–September period ranged from 237 to 252 mm for both the basin and its headwater portion, while the values measured from the landscape flux tower ranged from 366 mm (2014) to 458 mm (2016). Consequently, the comparison of monthly ET shows underestimation of the modelled ET (BESS) at both the basin and headwater scales compared to the ET obtained from flux tower measurements (Fig. 5b and c). Higher growing season water losses (ΔS_SUB-BASIN) in 2014 and 2016 observed for the sub-basins (Fig. 3) are consistent with the annual (hydrological year: October–September period) basin response, i.e. ΔS_BASIN (Fig. 6a).

Figure 5 Monthly comparisons of growing season (May–September 2014–2016, mm month⁻¹) water losses (evapotranspiration [ET] and runoff [Q]) between the Scotty Creek basin (x axis) and the sub-basins (y axis). (a) Q_BASIN and average (vertical error bar corresponding to minimum and maximum) Q estimates for the East, South and West sub-basins (Q_SUB-BASIN). Estimates of Q are obtained for the drainage area derived from automated terrain analysis using a digital elevation model. The symbol shapes (i.e. dot or triangle) indicate the number of months available to calculate mean sub-basin Q. No discharge data to calculate Q_SUB-BASIN are available in September 2016. (b) Basin and (c) headwater ET estimates obtained with the BESS model (ET_{BESS_BASIN} and ET_{BESS_HEAD}, respectively) compared with corresponding (y axis) landscape flux tower estimates of ET (ET_LAND). For (a)–(c), the continuous black line is the ordinary least square (OLS) regression. The OLS regression uncertainty (grey coloured band) is estimated using bootstrapping with 1000 iterations. The stippled black line is the 1:1-line.

Figure 6 (a) Annual (hydrological year: October–September, 1996–2022) water balances (mm yr⁻¹) for the Scotty Creek basin obtained from daily precipitation (P_BASIN) and rainfall measurements (R_BASIN) resulting in snow water equivalent (${\text{SWE}}_{\text{BASIN}}=P_{\text{BASIN}}-R_{\text{BASIN}}$), daily runoff (Q_BASIN) and evapotranspiration (ET) estimates from the BESS model (ET_{BESS_BASIN}). ET for the 1996–2001 period (dashed green line) corresponds to the 2002–2022 average period. Basin scale water storage change (ΔS_BASIN) is the difference between incoming and outgoing water fluxes. (b) Annual ratio of runoff to precipitation (i.e. the runoff ratio). The hashed area corresponds to the mean runoff ratio over the temporal period considered (1996–2002; 2003–2012; 2013–2022). For panels (a) and (b), the range of values for Q_BASIN, ΔS_BASIN and the runoff ratio corresponds to the lowest and highest basin drainage area estimates, i.e. 130 and 202 km².

3.5 Basin annual water balance

Over the 27 year (1996–2022) study period, annual water inputs were dominated by R, ranging from 111 to 324 mm (mean ± std, 243 ± 63 mm), while SWE_BASIN ranged from 81 to 181 mm (mean ± std, 112 ± 24 mm, Fig. 6a, Table S3). For water losses, annual ET estimated with the BESS model ranged from 223 to 311 mm (mean ± std, 261 ± 22 mm) over the 2002–2022 period (Fig. 6a). In comparison, annual Q_{BASIN_130} and Q_{BASIN_202} ranged from 26 to 317 mm (mean ± std=164 ± 81 mm) and from 17 to 204 mm (mean ± std=105 ± 52 mm) for the 2002–2022 period, respectively. Thus, annual ET was between 2.2 and 3.5 times higher than annual Q, given the range of drainage area estimates.

ET_{BESS_BASIN} and SWE_BASIN were relatively stable over time (261 ± 22 mm and 112 ± 24 mm, respectively, Fig. 6a). ΔS_BASIN, R_BASIN and Q_BASIN experienced higher between-year variability from 1996 to 2022 ($\mathrm{\Delta }{S}_{\text{BASIN\_130}}=-\mathrm{60}$ ± 75 mm, $\mathrm{\Delta }{S}_{\text{BASIN\_202}}=-\mathrm{5}$ ± 63 mm; R=243 ± 63 mm; Q_{BASIN_130}=155 ± 76 mm; Q_{BASIN_202}=100 ± 49 mm) than ET_{BESS_BASIN} and SWE_BASIN.

ΔS_{BASIN_202} and ΔS_{BASIN_130} ranged from −172 to 105 mm and from −95 to 121 mm, respectively. ΔS_BASIN decreased from ∼120 to 0 mm over 1996 to 2001, while Q increased from ∼30 to ∼140 mm. ΔS_BASIN was negative (∼100 mm) over the 2004–2014 period. Then, ΔS_BASIN was either positive or negative from 2015 to 2022 for both drainage area estimates (Fig. 6a).

The annual ratio of runoff to precipitation (i.e. the runoff ratio, Fig. 6b) ranged from 0.1 to 0.5 (runoff ratio_{_202}) and from 0.1 to 0.8 (runoff ratio_{_130}). The runoff ratio strongly increased from 1996 to 2002 (from ∼0.1 to 0.4–0.6 the runoff ratio_{_202} and runoff ratio_{_130} average were 0.2 and 0.3, respectively) followed by a period of higher and more stable values until 2012 (runoff ratio_{_202} and runoff ratio_{_130} average for 2003–2012 are 0.4 and 0.6, respectively). For the 2013–2022 period, the runoff ratio was more variable, but on average, lower (runoff ratio_{_202}=0.2 and runoff ratio_{_130}=0.4) than for the 2003–2012 period.

4 Discussion

4.1 Growing season water balance components in three small-scale basins of a boreal peatland complex: Objective 1

From mid-May until the end of September, the growing season water balances were dominated by water input and loss through rainfall and ET, respectively. Growing season daily ET ranged among values commonly observed elsewhere across the boreal biome, with higher wetland than forest ET (Arain et al., 2003; Isabelle et al., 2018; Nakai et al., 2013; Volik et al., 2021; Wu et al., 2010). For example, higher wetland (2.9 ± 1 mm d⁻¹) than forest ET (1.7 ± 0.6 mm d⁻¹) at Scotty Creek was reported for June–mid July 2013 (Warren et al., 2018), with transpiration from black spruce and tamarack accounting for only approximately 6 % to 12 % of forest ET (Perron et al., 2023).

The spring freshet contribution to growing season water losses was the lowest for the East sub-basin. Despite the uncertainty in the East sub-basin drainage area, the range of the wetland-to-forest ratio for the East sub-basin (0.34 to 0.84) was lower than for the two other sub-basins (South: 1.24 and West: 1.06). The larger forested portion in the East sub-basin compared to the other two sub-basins could lead to more post-spring freshet runoff, as the gradually deepening frost table can promote subsurface runoff (Sjöberg et al., 2021). In contrast, during the mid-growing season, wetlands can act as “gatekeepers”, reducing hydrological connectivity (Connon et al., 2015; Phillips et al., 2011). Land cover control over runoff dynamics has also been observed in other permafrost-affected basins, including a mountainous permafrost landscape where differences in vegetation types were shown to affect the rainfall–runoff relationship (Genxu et al., 2012).

Regarding monthly water balances, the high residuals observed in May for all 3 years (Fig. 4a–c) might be explained by the inclusion of snowmelt input through SWE that month. Due to limited data availability, SWE_MAX, estimated in late March just before the onset of snowmelt, served as a proxy for snowmelt input in the May water balance, highlighting the difficulty of appropriately accounting for the spring freshet in growing season water balance observations. To shed light on this challenge, we estimated the amount of snowmelt at the end of April using a simple temperature index model (Fig. S3 in the Supplement). The estimated snowmelt amounts (median [25th–75th percentiles] from 10 000 Monte Carlo simulations) at the end of April were 105 [78–136] mm in 2014, 187 [138–238] mm in 2015, and 125 [92–159] mm in 2016. These ranges correspond closely to the SWE_MAX measured each year (102, 167 and 128 mm in 2014, 2015 and 2016, respectively), suggesting that only a small portion of SWE_MAX contributed to the May water balance. This would reduce the high residuals in the May water balance in the West sub-basin, estimated as +149 (2014), +176 (2015) and +117 mm (2016).

Despite the observational challenges, particular attention should be paid to this snowmelt period, which is profoundly influenced by climate warming. Firstly, the spring freshet in recent years has been shown to occur earlier in the Arctic-boreal region compared to previous decades (Chasmer and Hopkinson, 2017; Mack et al., 2021; Pohl et al., 2007; Woo et al., 2008). At Scotty Creek, an earlier snowmelt of 16 d was observed during the 2000–2009 period compared to the 1970–1979 period (Chasmer and Hopkinson, 2017). The Scotty Creek basin hydrograph analysis consistently reveals an earlier increase in discharge (∼15 d) during the 2009–2022 period compared to the 1995–2008 period (Fig. S4 in the Supplement). Secondly, earlier snowmelt leads to a longer snowmelt period, as projected for the Liard River watershed, resulting in a more gradual snowmelt (Woo et al., 2008). However, an increase in wetland extent due to forested peat plateau collapse can contribute to shorter snowmelt period since snow melts faster in wetlands than in forest stands (Connon et al., 2021; Quinton et al., 2019). Shorter snowmelt periods can result in higher spring freshet peaks, as observed at Scotty Creek and the adjacent Jean Marie River meso-scale basin (Connon et al., 2021).

Except for May, the remainder of the growing season showed reasonably well closed monthly water balances with low residuals (Fig. 4), suggesting that obtaining water storage from measured wetland WTP and water content appears to be appropriate in low-relief landscapes such as the thawing boreal peatland complex in this study. To better understand the hydrological response of small- and meso-scale basins, we compared the hydrographs and monthly average runoff and ET estimates of the three headwater sub-basins, with corresponding estimates obtained at the basin scale, as described in the following section.

4.2 Small-scale basin evapotranspiration and runoff from a boreal peatland complex in a meso-scale basin context: Objective 2

The annual basin water balance (Fig. 6) had higher water losses in 2014 and 2016 than in 2015 (Fig. 3), similar to the growing season sub-basin water balances. Using independent data sets (i.e. sub-basin measurements and publicly available and modelled data for the basin), we observed that ET is the dominant annual water loss at both the sub-basin and basin scales, averaging more than twice the runoff. The hydrographs at both scales were comparable, i.e. they were dominated by the spring freshet peak, typical for regions with a subarctic nival regime (Gandois et al., 2021; Woo et al., 2008). However, an exception occurred during the driest year (2014), when the peak in basin runoff peak was more than 10 times lower than that for the sub-basins (Fig. 2). This difference might be partially explained by the higher proportional coverage of wetlands in the headwater sub-basins (∼40 %) compared to the entire basin (∼20 %) and the high coverage of mineral uplands in the basin (∼40 %; Chasmer et al., 2014). More water is expected to be stored in saturated wetlands than in mineral uplands (McCarter et al., 2020; Price, 1987), which may help sustain a higher runoff ratio during years with low late-winter SWE, as observed in 2014. The higher degree of saturation in wetlands compared to mineral uplands can favour surface runoff over water infiltration during the spring freshet, as observed in small-scale basins in Sweden (Jutebring Sterte et al., 2018, 2021). Dry conditions in 2013 (annual total P=259 mm, Fig. 6) may have further exacerbated the drying of mineral uplands compared to wetlands, thereby enhancing infiltration at the basin scale during the 2014 snowmelt.

For the concurrent monitoring period at both scales (2014–2016), sub-basin runoff agreed well with basin runoff (Fig. 5a). May showed the greatest difference, with values differing by a factor of two, highlighting difficulties in adequately measuring discharge during the spring freshet. Runoff discrepancies between the sub-basin and basin scales may also be partly attributed to time lag effects, e.g. the spring freshet peak was delayed (∼2–4 d) between the headwater sub-basin and the basin outlets (Fig. 2). Thus, a certain portion of headwater sub-basin runoff in late April might have been accounted for in May at the basin scale. Additionally, the observed runoff difference in May may partly reflect differences in snow depth and melt dynamics between the basin scale (130 and 202 km²) and the finer sub-basin scale (<1 km²), as snowpacks are often heterogeneous in forests and tend to melt more rapidly in wetlands (Connon et al., 2021; Nousu et al., 2024).

Our results show that modelled ET obtained with the BESS model at the basin scale underestimated (annually ∼100 mm) observed ET (Fig. 5b). Given that wetland ET is higher than forest ET (Helbig et al., 2016b), the underestimation of ET might be related to land cover heterogeneity at the basin scale. The northern, i.e. downstream, portion of the basin is dominated by mineral uplands with better drainage and is mainly covered by deciduous or mixed forest stands (Chasmer et al., 2014). ET estimations from a chemical method at the Scotty Creek basin scale consistently ranged from 280 to 300 mm yr⁻¹ for the 1999–2002 period (Hayashi et al., 2004). However, the modelled ET was lower than the observed ET at the sub-basin scale, probably underestimating the contribution of wetlands (Fig. 5c). Although this difference may stem from the tendency of the BESS model to underestimate the spatial variability of ET in wetland-rich landscapes such as boreal peatland complexes near the southern permafrost limit, we cannot disentangle the extent to which it reflects a general underestimation of ET versus a specific underestimation in wetlands.

4.3 Annual basin water balance in relation to changes in land cover and hydrological connectivity: Objective 3

Understanding the long-term runoff dynamics of thawing boreal peatland complexes remains challenging due to strong ecohydrological feedbacks (Shirley et al., 2022; Song et al., 2024; Walvoord and Kurylyk, 2016). Variability in precipitation regimes may have influenced runoff ratio dynamics, as suggested by the peak in the runoff ratio in 2020 (0.5–0.8), the rainiest year in the 1996–2022 period (Fig. 6b). Under wet conditions, ephemerally connected wetlands can increase the effective drainage area (Connon et al., 2015), whereas during dry periods, some wetlands become hydrologically disconnected, thereby reducing the runoff ratio. The lower runoff ratio in summer compared to the spring freshet is consistent with intensified wetland drying during the summer months (Figs. S5 and S6 in the Supplement).

The weak correlation between current-year effective precipitation (precipitation minus ET) and runoff (R²=0.2, Fig. S7 in the Supplement) suggests that other processes such as rapid changes in land cover and hydrological connectivity may have played a more dominant role in controlling runoff. Additionally, the cross-correlation analysis showed that current-year effective precipitation provides the best linear correlation with runoff, while antecedent wetness offers no explanatory power (Fig. S8 in the Supplement).

From a landscape perspective, in the headwater portion of the Scotty Creek basin, Haynes et al. (2022) estimated a 1.4 % forest loss between 2010 and 2018. Rapid permafrost thaw and the connection of wetlands to the drainage network are expected to increase hydrological connectivity, leading to an increase in permanent and transient runoff (Connon et al., 2014, 2015; Haynes et al., 2018). However, despite these changes, the average runoff ratio over the 2013–2022 period was lower than during the 2003–2012 period (Fig. 6b). Meanwhile, the drying of hydrologically connected wetlands has been reported at Scotty Creek between 2010 and 2018 (Haynes et al., 2018), which facilitated the development of individual hummock landforms, indicative of drier near-surface peat layers (Haynes et al., 2022). Wetland drying limits the saturation of permeable near-surface peat layers, which can coincide with a decrease in drainage efficiency. For example, the high hydraulic conductivity of the near-surface peat layer promotes more effective drainage compared to deeper peat layers (Ingram, 1978; Morris et al., 2011; Quinton et al., 2008). Vegetation succession occurring within approximately a decade following wetland initiation can lead to vertical peat accumulation above the water table. This elevation of the peat surface can contribute to reduced saturation in near-surface peat layers by hydrologically decoupling them from the saturated underlying peat (Errington et al., 2024). Therefore, the decrease in runoff ratio observed after 2012 can be attributed to reduced drainage efficiency resulting from wetland drying.

Our analysis indicates that the competing influences of wetland expansion (which increases hydrological connectivity) and wetland drying (which reduces hydrological connectivity) are key drivers of long-term runoff variability in boreal peatland complexes near the southern permafrost limit. Accordingly, sustained long-term hydrological monitoring is essential to disentangle the impacts of precipitation and land cover changes on the runoff ratio within such rapidly changing landscapes. A sub-basin-scale modelling study suggests that replacing 50 % of forested peat plateaus with wetlands leads to a reduced runoff ratio following wetland drainage. This change was attributed to increased surface storage capacity, reduced runoff efficiency, and higher landscape evapotranspiration, assuming no increase in precipitation (Stone et al., 2019).

In the boreal biome, wetlands exhibit higher mid-day ET than adjacent forests during the growing season (Helbig et al., 2020b). Projections for the 21st century indicate that wetland ET will exceed forest ET by more than 20 % across approximately one-third of the boreal biome under the Representative Concentration Pathway (RCP) 4.5 scenario and up to two-thirds under the RCP 8.5 scenario (Helbig et al., 2020a). While wetland expansion at the expense of forested peat plateaus increases ET, it may significantly influence both the water balance and the regional climate (Helbig et al., 2016b). These findings highlight the continued need for long-term measured and modelled ET comparisons, as ET is likely to play a crucial role in shaping the future water balance of boreal peatland complexes near the southern permafrost limit.

4.4 Effective versus potential drainage area: implications for water balance studies

Defining basin and sub-basin boundaries and drainage areas in low-relief landscapes such as vast swaths of the Taiga Plains using automated terrain analysis is challenging and estimates tend to vary, at least partly, depending on the DEM used (Al-Muqdadi and Merkel, 2011; Datta et al., 2022; Keys and Baade, 2019; Moges et al., 2023). Although difficult to apply across large regions, field observations are crucial in low-relief landscapes for accurately defining the effective drainage area (Connon et al., 2015). Our comparison of effective and potential drainage areas – based on field observations and automated terrain analysis of a DEM- showed that both estimates are consistent for the sub-basin almost entirely composed of connected wetlands (factor 1.2, West sub-basin, Fig. 1c). However, the two drainage areas exhibit important differences for the sub-basin with a high proportion of isolated wetlands (East sub-basin). There, the potential drainage area is five times higher than the effective drainage area. Field observations may lead to a more precise delineation of the effective drainage area contributing to the drainage network (Connon et al., 2015). However, regarding the growing season water balance for the East sub-basin (Fig. 3b), the water balance residual is 3.5 to 9.5 higher using the effective drainage area. In this case, the automated terrain analysis derived drainage area is more adequate to close the water balance. Subsurface water flows can occur at greater depths in permafrost-free basins (Sjöberg et al., 2021). Unobserved subsurface flows, such as through taliks, defined as perennially thawed ground below the active layer (Devoie et al., 2019), potentially lead to an underestimation of the effective drainage areas from field observations.

At the basin scale, automated terrain analysis produces different drainage areas (Burd et al., 2018; Chasmer and Hopkinson, 2017; Connon et al., 2014; Quinton et al., 2004; Water Survey of Canada), with the two most distinct estimates being used in this study (i.e. 130 and 202 km²). The increase in wetlands hydrologically connected to the effective drainage area due to permafrost thaw is expected to be captured by the substantial increase in the runoff ratio from 1996 to 2012 (Fig. 5). Delineating drainage areas at the sub-basin and basin scales remains a challenge, with proportionally larger errors in smaller areas such as the East, West and South sub-basins. Thus, minor differences in landscape heterogeneity (e.g. the hydrological connectivity of wetlands to the drainage network) may lead to large variations in the drainage area. Given the rapid permafrost thaw and associated land cover changes occurring near the southern permafrost limit (Quinton et al., 2019), improving the constraints on the hydrological connectivity of small low-relief (sub-)basins is essential for accurately quantifying and modelling water and carbon losses (Gao et al., 2018; Wei et al., 2024).

4.5 Constraining the water balance in thawing boreal peatland complexes: broader implications and perspectives

Non-linear hydrological responses, such as changes in the runoff ratio, ET and water table position, to variations in precipitation and hydrological connectivity driven by permafrost thaw are linked to shifts in soil physical properties, microbial communities and vegetation composition and structure. These interconnected changes collectively influence ecosystem services at multiple scales, including the local (e.g. subsistence activities), regional (e.g. water storage) and global levels (e.g. carbon storage as reflected in the net ecosystem carbon balance [NECB]; Camill et al., 2001; Chapin et al., 2006; Ernakovich et al., 2022; Jones et al., 2022; Li et al., 2023; Shirley et al., 2022). Assessing whether thawing boreal peatland complexes act as a net source or sink of carbon (NECB), once both vertical and lateral fluxes are considered, is therefore an important avenue of research (Song et al., 2024). For example, a recent review showed that the dissolved organic carbon concentration can be elevated in sporadic and discontinuous permafrost areas and tend to increase with permafrost thaw (Heffernan et al., 2024). Thus, understanding the mechanisms driving runoff, such as the spring freshet, is essential for quantifying lateral carbon exports to NECB (Chapin et al., 2006; Gandois et al., 2021; Laudon et al., 2004).

Long-term hydrological monitoring is also essential for understanding how gradual changes (e.g. vegetation shift or increasing T_air) are interlinked with more frequent and intense pulse disturbance events (e.g. weather extremes, abrupt permafrost thaw or wildfires) (Li et al., 2023). Wildfires have been shown to accelerate permafrost thaw (Gibson et al., 2018), posing an increasing threat to ecosystem services. The year 2023 set a record for surface burned across Canada (MacCarthy et al., 2024; Wang et al., 2024). As water table position and moisture can constitute an indicator of fire risk, understanding the water balance dynamics of peatland-dominated basins may help in managing fire risk (Kartiwa et al., 2023; Mortelmans et al., 2024). In October 2022, the Scotty Creek basin was impacted by a late-season wildfire. While the wetland flux tower and several cutthroat flumes remained intact (Fig. S1), the landscape flux tower was destroyed and rebuilt in March 2023. Our work, which contributes to understanding the hydrological response of a rapidly thawing boreal peatland complex, can serve as a baseline for understanding the combined effects of permafrost thaw accelerated by wildfire.

5 Conclusions

This study contributes to a better understanding of the hydrological response of small-scale basins (here: “sub-basins”) within the headwater portion of a meso-scale basin (here: “basin”) in the Taiga Plains in northwestern Canada. We provide insights into how the hydrological responses of rapidly thawing boreal peatland complexes – at both the sub-basin and basin scales – are shaped by complex factors (e.g. changes in land cover and hydrological connectivity) that extend beyond year-to-year changes in precipitation and ET. Specifically, we find that:

• determining runoff in low-relief landscapes such as thawing boreal peatland complexes is challenging because

  • –

    sub-basin and basin boundaries and the resulting drainage areas must be approached with caution, since permafrost ridges act as barriers isolating wetlands from the effective drainage network, and

  • –

    of difficulties in integrating spring freshet runoff into the growing season water balance.

• the small-scale headwater portion is representative of the corresponding meso-scale basin. At both scales, our analysis shows that

  • ET is the dominant water loss, on average more than double that of runoff,

  • the temporal dynamics of the growing season (sub-basin) and annual water balance components (basin) are similar,

  • spring freshet peaks are similar, except for the driest year, when the basin runoff is more than 10 times lower than the sub-basin runoff, and

  • spring freshet contributions to runoff in the April–September period are similar.

• long-term changes in the basin-scale runoff ratio cannot be explained by precipitation and ET alone. The increase in the runoff ratio from 1996 to 2012 likely reflects enhanced hydrological connectivity and wetland drainage. In contrast, the shift to a lower mean runoff ratio from 2013 to 2022 may be attributed to wetland drying following loss of connection to the drainage network. We propose that wetland drying, observed at the headwater sub-basin scale, accounts for the declining runoff ratio at the basin scale. At the same time, new isolated wetlands are forming, and additional wetlands may be becoming connected to the drainage system. Thus, the observed changes in runoff ratio likely reflect the competing influences of wetland drying and the emergence of new hydrological connections.

Appendix A

Table A1 List of all the variables and expressions used in this study (left column), alongside the corresponding abbreviations (right column).