Uncovering the melt: Unmanned Aircraft System (UAS) and in-situ sensor synergies reveal dissolved organic carbon (DOC) pathways in a northern peatland

Spring snowmelt is a critical period for dissolved organic carbon (DOC) export from northern boreal peatlands, yet the spatiotemporal dynamics of this process remain poorly understood. To reveal the spatial patterns, we used a novel combination of high-resolution Unmanned Aircraft System (UAS) snow depth mapping, topographic wetness index, and high-frequency stream monitoring during the peak snowmelt in 2024. The study took place in a small 6 ha fen in northeastern Finland. Our results show that DOC leaching is triggered after widespread snow cover depletion, likely due to thawing of surficial peat layers and increased connectivity in the peatland. High-resolution UAS snow surveys captured the progression of snowmelt from drier, south-facing slopes and forested areas toward wetter fen areas, with the expansion of snow-free areas in high-wetness zones initiating hydrological connectivity and pathways for DOC transport. Event-based hysteresis and flushing analyses enabled by high-frequency stream monitoring revealed transitions from deeper to more surficial flow paths towards the final peak melt. The integration of high-resolution spatial and temporal datasets enabled the detailed identification of DOC transport mechanisms during the snowmelt period. These findings underscore the sensitivity of peatland carbon dynamics to late winter processes and snow conditions, highlighting their potential vulnerability to future shifts in climate.

1 Introduction

Northern peatlands have been recognized as key components in the carbon cycle due to their extensive coverage in northern latitudes (30–90° N) and long-term capacity for carbon accumulation (Alexandrov et al., 2020; Qiu et al., 2020; Xu et al., 2018). In addition to storing carbon, they are also hotspots for dissolved organic carbon (DOC) export on the catchment scale (Croghan et al., 2024; Laudon et al., 2011; Rosset et al., 2019; Tang et al., 2018) and globally, with 58 % of the global peatland DOC exports originating from boreal regions (Rosset et al., 2022). Aquatic carbon export not only contributes substantially to the peatland net ecosystem carbon balance (Nilsson et al., 2008; Tong et al., 2025) but also has a significant role in downstream water quality and biogeochemical cycling (Kritzberg et al., 2020; Oni et al., 2013), and can lead to increased CO₂ emissions caused by mineralization and evasion downstream (Dinsmore et al., 2010; Serikova et al., 2018). The current concern is that climate change and the associated hydrological shift will alter the peatland carbon balance, although the hydrological responses are complex and remain hard to predict (Qiu et al., 2020; Rosset et al., 2022; Zhang et al., 2022). Contrasting responses may occur in one area or even within one peatland complex (Zhang et al., 2022), highlighting the heterogeneous nature and variable hydrological connectivity at peatlands. Still, the spatiotemporal carbon transport processes are poorly understood (Griffiths et al., 2017; Werner et al., 2021), especially during the snowmelt period. Due to strong seasonality in water and carbon cycles, northern peatland ecosystems are sensitive to ongoing climatic change, particularly to snow conditions (Campbell and Laudon, 2019; Croghan et al., 2024; Marttila et al., 2021). Thus, there is an urgent need for documenting and understanding the relevant processes. Particular focus should be given to improved monitoring to facilitate a process understanding during critical periods for ecohydrological connectivity, such as spring snowmelt (Burd et al., 2018; Marttila et al., 2021).

In northern catchments, spring snowmelt dominates annual DOC exports, and studies have found that it contributes 30 %–64 % of annual fluxes (Dyson et al., 2011; Leach et al., 2016; Olefeldt et al., 2013), despite covering only 10 %–15 % of the year (Kortelainen et al., 1997). Climate change is expected to severely modify winter conditions in northern areas, including changes to snow conditions, namely snow accumulation (Pulliainen et al., 2020), duration of snow cover (Luomaranta et al., 2019), and snowmelt timing and rate (Musselman et al., 2017). This will have significant consequences for the catchment hydrological regime, with earlier and less severe spring floods placing more importance on winter streamflow (Teutschbein et al., 2015). As hydrology exerts first-order control on DOC export (Wen et al., 2020), changes in snow conditions and snowmelt timing will directly affect the timing and magnitude of DOC exports (Campbell and Laudon, 2019; Laudon et al., 2013). Changing winter conditions may also impact seasonal soil frost and freeze-thaw cycles, which can modify hydrologic flow paths and soil DOC concentrations (Ågren et al., 2010; Campbell and Laudon, 2019), affecting the amount of DOC that is available to be transported during snowmelt (Ågren et al., 2012; Campbell et al., 2014). Enhanced soil frost during winter has been found to cause higher and delayed DOC peaks during snowmelt (Ågren et al., 2012; Campbell et al., 2014; Tiwari et al., 2018), and more frequent freeze-thaw cycles can alter soil DOC concentrations (Yu et al., 2011), peat hydraulic properties and water flux (Liu et al., 2022). Alternatively, reduced seasonal soil frost could lead to the earlier activation of subsurface flow paths, increased connectivity, and DOC transport during snowmelt (Croghan et al., 2023; Hinzman et al., 2020; Laudon et al., 2011). These complex trajectories highlight the need to understand how snowmelt impacts DOC transport processes.

Peatlands should be considered as mosaics rather than uniform landscape units, as they are characterised by spatial variability in, for example, hydrological connectivity, groundwater influence, vegetation, and microtopography (Isokangas et al., 2017; Rosset et al., 2020; Vitt et al., 2022). Hydrological connectivity, mainly defined by wetness and water levels, represents a key factor regulating DOC transport (Knapp et al., 2022; Prijac et al., 2023), which potentially varies, both spatially and temporally, very quickly during snowmelt in peatlands where the water table is naturally high (Isokangas et al., 2017; Laudon et al., 2011; Peralta-Tapia et al., 2015). To assess the spatial variability, topographical wetness index (TWI) has been proposed as a tool to identify hydrological connectivity and active source zones for DOC transport (Knapp et al., 2022; Richardson et al., 2012; Werner et al., 2021). During snowmelt, activation of sources and carbon transport processes may also be impacted by snow conditions and timing of melt (Campbell and Laudon, 2019; Croghan et al., 2023; Marttila et al., 2021). However, to our knowledge, spatiotemporal variations in snow cover have not yet been combined with TWI to investigate DOC transport during snowmelt.

Snow depth and melt rates can express local or small-scale spatial variation caused by microtopography and vegetation (Litaor et al., 2008; Sannel, 2020; Shirley et al., 2025), and wind processes causing redistribution of snow (Meriö et al., 2023; Shirley et al., 2025). Traditional point-scale snow measurements often fail to capture small-scale spatial variation in snow cover, have limited spatial coverage (tens to hundreds points), and thus can lead to severe over or underestimation, especially during the snowmelt period, as spatial variation in snow depth typically increases towards the end of the snow-covered period (Grünewald et al., 2010; Meriö et al., 2023). Spatial distribution of snow cover can be efficiently monitored with the latest remote sensing technology. However, satellite data is limited by coarse spatial or temporal resolution and issues with cloud cover (Dong, 2018). Thus, novel high-resolution monitoring techniques are needed to assess small-scale snow cover heterogeneity in a resolution relevant to peatland microtopography (< 1 m) (Shirley et al., 2025; Steenvoorden and Limpens, 2023) and to improve spatial representativeness. For this purpose, Unmanned Aircraft Systems (UAS) can serve as a cost-effective solution for data collection at the required spatial and temporal resolution. Many studies have utilized UAS Light Detection and Ranging (LiDAR) and Structure-from-Motion (SfM) techniques for snow depth mapping, achieving better than 10 cm accuracies in open areas using Real-Time Kinematic (RTK) positioning (Bühler et al., 2016; Harder et al., 2016, 2020; Rauhala et al., 2023; Ylönen et al., 2025). Despite the wide application of UAS snow depth mapping, and the crucial role of snow coverage, depth, melting, and the timing of these in ecohydrological processes (Marttila et al., 2021), high-resolution snow cover data is rarely incorporated into biogeochemical studies.

During hydrologically dynamic seasons, such as snowmelt, DOC transport processes may vary over very short time scales, making the full dynamics unrecognizable with low-frequency (daily to weekly) sampling methods (Pellerin et al., 2012; Tunaley et al., 2016). This highlights the importance of high-resolution monitoring efforts, particularly in northern peatland headwater streams, where the role of storm events and snowmelt in DOC export is pronounced (Burd et al., 2018). Development in in-situ optical sensors has allowed the extension of the temporal resolution of measurements to match many hydrological and biogeochemical process rates (Bieroza et al., 2023), and there has been growing interest in using high-frequency monitoring to capture seasonal, event-scale, and sub-daily dynamics (Croghan et al., 2024; Gómez-Gener et al., 2021; Shatilla et al., 2023; Tunaley et al., 2018). Analyzing concentration-discharge (C–Q) relationships during storm events can provide insights into the active flow pathways and source areas contributing to DOC transport (Lloyd et al., 2016; Vaughan et al., 2017) as the behavior of the solute concentration relative to discharge can be considered to reflect catchment conditions, hydrological connectivity and event characteristics (Knapp et al., 2022). However, the C–Q relationship is an indirect measurement of these properties, and the typical one-point measurement deployed at the outlet stream can leave spatial processes and active sources unrecognized (Bieroza et al., 2023; Knapp et al., 2022). Thus, more efforts are needed to combine high-resolution spatial and temporal datasets in order to fully understand the link between the two scales.

Although some studies have deployed high-resolution monitoring in peatland headwater streams (Prijac et al., 2023; Rosset et al., 2019; Tunaley et al., 2018), more high-resolution monitoring programmes need to be established closer to the peatland outlets to accurately capture and understand their critical role in landscape carbon exports (Rosset et al., 2019, 2022). In this study, we deployed high-frequency in-situ DOC and hydrological monitoring at small peatland trickle along with UAS mapping during the peak snowmelt period of 2024 at Puukkosuo fen at the Oulanka research station, northeastern Finland. We aimed to (i) document spatiotemporal variations in snow cover melt and hydrological connectivity using frequent UAS surveys; (ii) link the spatial snowmelt processes to DOC export using high-frequency stream monitoring; and (iii) building on the novel datasets, identify the key spatiotemporal hydrological and DOC transport processes occurring during the spring melt period.

2 Methods

2.1 Study area

The research site, Puukkosuo fen, is located in the Oulanka National Park (66.38° N, 29.31° E), northeastern Finland. The mean annual air temperature measured at the Oulanka research station is 0.7 °C, and annual precipitation is 559 mm. The area is characterized by relatively long, snowy winters. The average duration of the snow cover is 194 d, typically ending in May, and the snow depth reaches a maximum of 80 cm on average (1995–2024, Finnish Meteorological Institute (FMI) open data). Puukkosuo is a sloping fen peatland located within the north-boreal aapa mire zone. It receives water inputs as groundwater discharge from the hills on the northwestern and northern edges, as well as surface runoff from the upslope catchment area (Fig. 1b), covering an area of 35 ha and comprising 74 % forest and 24 % peatland. Waters from Puukkosuo eventually drain to the Oulankajoki River. The fen covers an area of approximately 6 ha with peat depths ranging from less than 1 meter near the borders and up to 4.6 m in the center. Due to consistent groundwater discharge, Puukkosuo has a porewater pH close to neutral and provides habitat for several endangered species. Vegetation is mainly open or sparsely tree-covered, dominated by Carex sedges, brown and peat (Sphagnum) mosses (Järvi-Laturi et al., 2025), and some Scots pines.

Figure 1 Map of the study area and measurement locations at Puukkosuo fen (a) and upper catchment boundaries (b). Data © National Land Survey of Finland (NLS) Topographic database 2022.

Puukkosuo fen is part of the EcoClimate environmental monitoring programme maintained by the Oulanka research station, University of Oulu, and is equipped with extensive research infrastructure, including high-resolution hydrological and meteorological monitoring and snow measurements. The groundwater level and temperature at the fen were measured in a dip well with a piezoresistive pressure sensor (Decentlab GmbH) measuring at a 10 min interval. Continuous water level measurements were corrected and converted into water table depth (WTD) using monthly manual water level measurements carried out during the snow-free period as a reference. In addition, air temperature and precipitation were monitored continuously (1 min interval) at the research site, and precipitation samples were collected at the Oulanka research station after rain or snowfall events for isotope analysis. Precipitation sampling is maintained by the Oulanka research station, approximately 800 meters from the Puukkosuo fen.

2.2 In-stream monitoring

Water quality was measured continuously (15 min intervals) at the stream gauging station installed in the peatland outlet, where TriOS OPUS spectral sensor (TriOS GmbH) was installed for water quality measurements. The OPUS sensor measures the full absorption spectrum from 200 to 360 nm, from which spectral analysis is performed to provide readings for the concentration equivalents (eq) of DOC and total suspended solids (TSS). The sensor was cleaned regularly to prevent fouling. Flow measurements (1 min intervals) were obtained with a 1 in. Parshall flume and UGT502 ultrasonic water level sensor (Ifm electronic GmbH), measuring water level at the crest of the flume and converted it into discharge (Q in L s⁻¹) using the equation provided by the manufacturer.

In the same location as in-stream monitoring, weekly grab samples were taken to the laboratory for DOC analyses, which were used to calibrate TriOS OPUS measurements. Samples were first filtered with 1.2 µm GF/C filters (Whatman) and acidified with 2 M HCl before being stored in the dark at 4 °C until analysis. DOC concentration was analysed in the Oulanka research station laboratory (Shimadzu TOC-L Analyzer). As laboratory measurements for TSS were not available, the TSS used in this study will refer to TSS equivalent (TSSeq).

In addition, daily water samples were collected for stable water isotope analysis. Samples were collected in 15 mL plastic tubes and stored in the dark at 4 °C until analysis. Any samples with visible particles or organic matter were filtered with 0.45 µm syringe filters before analysis. The isotopic composition of oxygen (δ¹⁸O) and hydrogen (δ²H) was analyzed at the University of Oulu laboratory with a Picarro 2140-I cavity ring-down spectrometer (CRDS) and were standardized to Vienna Standard Mean Ocean Water. The average measurement precision for all δ¹⁸O samples used in the study was 0.02 ‰. Measurement accuracy calculated based on the mean of measured standards was 99.7 %.

2.3 UAS surveys and snow measurements

The snow data used in the study was collected through UAS flights and sampling in May 2024. Five UAS campaigns were completed during the peak snowmelt period on 1, 14, 15, 16, and 17 May. In addition, UAS data collected on 31 May was used as a snow-free reference in constructing snow depth maps. All UAS surveys were performed using DJI Mavic 3M, at a 60 m altitude with 70 % side and 80 % frontal overlaps, resulting in approximately 1.5 cm ground sampling distance for RGB. Mavic 3M was equipped with a Real-Time Kinematic (RTK) GNSS (Global navigation satellite system) module, producing geotagged images with centimeter-level accuracy. In addition, during UAS campaigns, 5–6 ground control points (GCPs) constructed of painted aluminium plates were surveyed with a Emlid Reach RS3 high-precision RTK GNSS receiver for georeferencing. Snow depth at the site was also monitored as a point measurement with a 10 min interval using a USH-9 ultrasonic snow depth sensor (SOMMER Messtechnik GmbH). In addition, manual snow depth measurements were done during each snow-covered campaign from a snow course transecting the peatland (Fig. 1a).

During the field campaigns, snowpack samples were collected for stable water isotope analysis. The sample was collected by coring the whole snow profile from the top to the base of the snowpack with a snow tube corer (diameter 3.5 cm). The sample was put into a plastic bag, where it was transferred to 15 mL plastic tubes once completely melted. Also, snow density and snow water equivalent (SWE) were measured at the same sampling location with a snow water equivalent sampler (diameter 10 cm).

2.4 Data Analysis

2.4.1 UAS data processing

Images taken during each survey were imported and processed in Agisoft Metashape 1.7.3, a processing software employing the Structure-from-Motion (SfM) technique with Bundle Adjustment (BA) method to create high-resolution digital elevation models (DEMs) (Westoby et al., 2012). Processing followed the methods and parameterization described in Ikkala et al. (2022). First, all images with motion blur or low quality were identified and removed. Then, images were aligned to create photogrammetric tie points, generating a sparse point cloud (James et al., 2017), and ground control points (GCPs, see Sect. 2.3) were imported as markers and manually verified for all found projections. Although all flights had RTK on board, further georeferencing with one or more GCPs is recommended to correct possible elevation bias (Rauhala et al., 2023). 5–6 GCPs were used, depending on the survey, and for each, two of the GCPs were selected as check points, which are not included in georeferencing but are used to validate the accuracy of the produced model. After georeferencing, poor accuracy tie points were filtered, and camera optimization was performed. Finally, a point cloud was generated with high-quality and moderate depth filtering settings and exported in an LAS format. Additionally, orthomosaics were generated from the point cloud and exported in a TIF format (Fig. S1 in the Supplement). Processing information and final model accuracies for each UAS survey, according to Agisoft Metashape, can be found in Table 1.

Table 1 UAS data processing information and ground control point (GCP) root mean square errors (RMSE) for control and check points according to Agisoft Metashape.

Further processing and data analysis were conducted using R version 4.4.0 in RStudio. First, statistical outlier removal (SOR), together with a cloth simulation filter (CSF), was applied to filter non-ground points from LAS point clouds. CSF is a point cloud filtering method simulating a cloth being dropped on the inverted point cloud surface, and the intersection of the cloth nodes and corresponding points represents the ground (Zhang et al., 2016). As CSF is sensitive to outlier points, it is vital to use an outlier filtering method before implementing ground classification (Zeybek and Şanlıoğlu, 2019; Zhang et al., 2016). Based on multiple trials, the best parameters selected for SOR were k = 20, m = 5, and for CSF were slope processing = TRUE, rigidness = 1, class threshold = 0.1, and cloth resolution = 1. R Packages LidR (Roussel and Auty, 2016) and RCSF (Roussel and Qi, 2018) were used for point cloud processing. Finally, filtered point clouds were rasterized to 10 cm resolution DEMs with package terra (Hijmans, 2020). Snow depth maps were calculated with a method known as DEMs of Difference (DoD) (Rauhala et al., 2023) by subtracting the snow-covered DEM from the bare ground reference. Snow depth mapping with SfM photogrammetry and DoD method is sensitive to errors caused by environmental factors during the survey, as well as underlying vegetation (Harder et al., 2016). Thus, DoDs (i.e., snow depth maps) were aggregated to a 30 cm cell resolution, reducing the noise and outliers but still maintaining sufficient resolution to identify small-scale variation in snow depth. The accuracy of the produced snow depth maps was estimated by comparing the UAS-produced snow depth to manual snow course measurement points. Differences in manual snow course and UAS snow depth measurements were used to calculate root mean square error (RMSE) metrics for each survey. The accuracy of the manual snow depth measurements was not considered when calculating the error metrics.

Using the UAS-generated snow depth maps, we calculated the proportion of snow-covered area (SCA), mean and median snow depth, and the change in snow depth (ΔSD) between surveys. Preliminary analysis of the snow depth maps revealed that the CSF algorithm struggled to identify the ground and non-ground points, especially in the later snow surveys where part of the ground was already snow-free, thus misclassifying vegetation as ground and causing positive artifacts in snow depth maps. This caused some pixels to have a positive change in snow depth (ΔSD) between the surveys. These positive artifacts were interpreted as no change, as they can be assumed to occur in areas where the snow depth had reached zero. All negative artifacts in SD were disregarded in further snow statistic calculations. The estimate of SWE for each raster cell was calculated by multiplying the snow depth by the field-measured snow density. Although manual snow density and SWE were measured as a one-point measurement, the method should provide a sufficient estimate for SWE for the study area since it is fairly homogenous, and variation in snow depth can be expected to override spatial variation in snow density (López-Moreno et al., 2013; Whittington et al., 2012). The produced SWE maps were used to calculate mean and median SWE, and the total volume of stored water in the snow in the study area.

DEM produced from pre-classified 5 points m⁻² LiDAR data surveyed during snow-free period from National Land Survey of Finland (NLS) (2024) was used to delineate the catchment area for Puukkosuo fen and understand the soil moisture and flow path patterns. Although flow paths in peatlands are complex and cannot be completely explained by surface topography, topographical wetness indices (TWI) can serve as a proxy for soil moisture and hydrological connectivity (Knapp et al., 2022; Richardson et al., 2012), particularly under wet conditions (Ikkala et al., 2022; Riihimäki et al., 2021). In this study, the DEM of the whole catchment area was used to calculate the Saga wetness index (SWI) (Böhner et al., 2006), which is a modification of the commonly used TWI. SWI is calculated as:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{SCA}}_M=S_{\max }{\left({\displaystyle \frac{\mathrm{1}}{t}}\right)}^{\mathit{\beta }\exp \left(t^{\mathit{\beta }}\right)}\mathrm{for}S<S_{\max }{\left({\displaystyle \frac{\mathrm{1}}{t}}\right)}^{\mathit{\beta }\exp \left(t^{\mathit{\beta }}\right)}\text{(2)}& {\displaystyle }\mathrm{SWI}=\ln \left({\displaystyle \frac{{\mathrm{SCA}}_M}{{\tan }_{\mathit{\beta }}}}\right)\end{array}$$

where SCA is the specific catchment area, t is a suction parameter, and β is the terrain slope angle. The difference between SWI and traditional TWI is that the SWI algorithm uses iterative modification of the multiple flow direction method (Freeman, 1991) for calculating the specific catchment area for each grid cell (Böhner et al., 2006). The t parameter is defined by the user to control the suction effect, i.e., the capillary attraction in the soil voids, simulating the suction or draw of flow from adjacent pixels, which is expected to increase in flatter terrain (Riihimäki et al., 2021; Winzeler et al., 2022). This approach makes SWI more suitable for flat areas such as peatlands (Böhner et al., 2006; Kopecký et al., 2021). Sink filling, calculation of SWI, and flow accumulation raster based on multiple flow direction algorithm were implemented for a 0.5 m resolution DEM in SAGA GIS 7.8.2 software, following the method described in Ikkala et al. (2022). We used the SWI calculated based on snow-free topography combined with UAS snow depth maps to outline the spatial distribution of areas that are more likely to contribute to the connectivity of DOC transporting flow paths once snow cover disappears. The SWI cells were further classified as low, medium, and high based on quantiles of SWI values to identify changes in snow cover in areas differing in topographically derived wetness conditions. For this purpose, snow depth maps were resampled to a 50 cm resolution to match the SWI raster.

2.4.2 Time-series analysis

Although the Trios OPUS sensor was cleaned automatically and manually, at some time periods, sensor fouling was observed to affect the water quality measurements, resulting in a gradual increase in DOCeq and TSSeq values, followed by a sharp decline after sensor cleaning. We used a systematic method to correct the data values during the fouling periods. First, the periods of fouling were identified from the data using the maintenance visit records and then corrected with linear scaling to match the measurement preceding the sensor cleaning to the post-cleaning measurement. Long-term drift was corrected by using a smooth LOESS fit for the sensor time series and weekly grab samples. The smooth sensor trend attributed as drift was removed from the sensor data and replaced with the grab sample LOESS trend, yielding a drift-corrected DOC time series. LOESS with a very smooth span was used for both sensor and grab samples (0.6 and 0.5, respectively) to ensure that only the long-term drift is corrected while preserving fine-scale variability. A rolling average with a 1-hour window size was then used to smooth out the sensor noise. Local calibration for scaling DOCeq measurements to DOC mg L⁻¹ was then created by plotting DOCeq against laboratory-measured weekly DOC grab samples (DOC = 1.60333 + 0.68178 × DOCeq, n = 83, Fig. S2). The modest R² value (0.42) is likely due to the interference of particulate material and the low number of grab samples, as we noticed the calibration significantly improved with more data points. However, the main focus of this study relies primarily on relative changes in DOC concentration and DOC-discharge relationships rather than absolute concentration values, and therefore uncertainties related to sensor calibration are unlikely to affect the main interpretations. All in-situ monitoring data were harmonized to hourly averages for further analysis.

Any missing values in Q or DOC, covering 5.9 % and 6.3 % of the study period for Q and DOC, respectively, were interpolated to produce a continuous time series for further analysis. Baseflow was first separated from the hydrograph by using the Lyne and Hollick method with R package grwat (Samsonov, 2022). Events were defined from the streamflow data as periods when discharge exceeded baseflow by 20 % (Blaen et al., 2017; Lloyd et al., 2016) for at least a 10 h duration. The 10 h minimum duration appeared to provide a suitable threshold for detecting events in the small peatland outlet stream, where the stream hydrograph can be flashy, especially during the spring freshet. If the detected event was double-peaked, it was split into two separate events, ensuring that each event included a rising and falling limb (Lloyd et al., 2016). Event delineation resulted in 14 events identified during the study period, of which two were disregarded due to the absence of distinguishable rising and falling limbs, resulting in unclear loop geometry. For the remaining 12 events, maximum discharge, mean, minimum, and maximum DOC concentration, total DOC load, and hourly DOC load (DOC load divided by event duration) were calculated.

The hysteresis index (HI) was calculated for each event following the method outlined by Lloyd et al. (2016), where Q and DOC were first normalized for each individual event, and HI was calculated at every 5 % of normalized discharge. HI provides a value ranging between −1 and 1, where negative values correspond to anticlockwise hysteretic patterns and positive values to clockwise hysteretic patterns. HI close to 0 indicates a lack of hysteretic relation, while HI approaching −1 or 1 indicates a strong hysteretic relation (Lloyd et al., 2016). The terms clockwise and anticlockwise refer to the direction of the hysteresis loop when concentration is plotted against discharge. Clockwise hysteresis is characterized by the peak concentration occurring during the rising limb of the event, whereas anticlockwise hysteresis occurs when the concentration peak lags the discharge. In the context of DOC transport, clockwise hysteresis is associated with quickly activating flow paths or proximal sources, which may lead to the exhaustion of carbon supply towards the end of the event, while anticlockwise patterns may suggest distal sources of carbon, the activation of DOC transporting flow paths later in the event, and/or different transit times for water and DOC (Lloyd et al., 2016; Vaughan et al., 2017; Williams, 1989). Thus, the hysteretic pattern can provide insights into the distance of DOC sources or transport time during individual events, which may reflect differences in flow routing, connectivity, or transport pathways (Knapp et al., 2022), but do not uniquely identify specific source areas or their concentrations.

The flushing index (FI) was calculated to explore the behavior of DOC on the rising limb of each event (Vaughan et al., 2017). FI is defined as the difference between the concentration at the point of peak discharge and the concentration at the beginning of the event, making it equal to the slope of the line. Similar to HI, FI also ranges from −1 to 1, where negative values suggest a decreasing DOC concentration on the rising limb (i.e., a diluting effect) and positive values indicate an increase in concentration (i.e., the flushing of DOC), which may be linked to the source supply or its proximity (Vaughan et al., 2017). Positive FI may occur when plentiful proximal sources are transported early during the event and/or from the transport limitation of solute, while diluting can occur when DOC sources are limited, quickly depleted, and/or insufficient to compensate for high runoff (Gómez-Gener et al., 2021; Shatilla et al., 2023; Vaughan et al., 2017).

We used Spearman correlation to explore relationships between event characteristics and mean SWI of the snow-free areas, which was used as a proxy for the wetness state and potential connectivity of the activated areas. Only events after the first UAS survey (1 May) were considered in the mean SWI and event correlation analysis. For the events that did not have UAS snow-cover data collected on the same day as the onset of the events, mean SWI was estimated with linear interpolation. These events took place between the first survey on 1 May and the second survey on 14 May, which had only a small difference in snow-covered area and mean SWI, and thus linear interpolation was expected to produce a reasonable estimate.

The contribution of meltwater in streamflow during the peak melt period (1–20 May) was evaluated by two-component hydrograph separation, which is a common method for quantifying the water sources based on a mass balance approach (Klaus and McDonnell, 2013). Data gaps (max. 2 d) in stream δ¹⁸O samples were interpolated to daily time steps. For the pre-event component, we used an average of weekly δ¹⁸O values sampled from the spring (Fig. 1a) during May to represent the groundwater signal. The spring remains unfrozen during the winter and is assumed to represent local groundwater. For the meltwater component, we used daily δ¹⁸O values from snowpack, sampled 5 times during the peak melt period. The gaps in sampling were filled either by using linear interpolation or using the value of the closest previous sample for the days after the last snowpack sample. Using time variant snowpack isotopic signal during the peak melt period has been shown to provide reasonable estimates of the meltwater contribution (Noor et al., 2023). Also, the following assumptions were made: (i) streamflow is generated from groundwater and meltwater only, (ii) the groundwater signal remains constant in space and time, (iii) no meltwater is stored temporarily in the soil, and (iv) the isotopic signal of groundwater and meltwater is significantly different (Klaus and McDonnell, 2013; Noor et al., 2023).

3 Results

3.1 Spatiotemporal variations in snow depth

Snow cover and melting exhibited spatial variation associated with tree cover, topography, and wetness. On the first survey on 1 May, the highest snow accumulation was observed around the borders of the fen, whereas the central area, areas with higher tree cover, and the south-facing slopes had the lowest snow depth (Fig. 2). During the peak melt period, when daily surveys were conducted (14–17 May), the hillslopes and areas with higher tree cover were the first to become snow-free on 15 May, followed by the central areas of the fen on 16 May. By 17 May, the central area had nearly completely melted, and only residual patches of snow remained, primarily along the boundaries of the fen. Spatial variation in snow depth was greatest during the first survey on 1 May and decreased towards the end of the study period (Fig. S3). In addition to natural snow distribution patterns, some anthropogenic influences were also evident in the UAS snow depth maps. Snowmobile tracks crossing the fen are visible on snow depth maps, creating consolidated and slowly melting tracks in the snow. Based on observations on-site, the snowmobile tracks developed into thick ice covers while melting proceeded. Overall, while UAS-derived snow depth maps demonstrated reasonable accuracy, their performance was highly dependent on survey conditions, as expected. The snow depth maps had a consistent negative bias (Mean Error), i.e., they tended to underestimate the snow depth compared to manual snow measurements (Table S1). When considering all surveys, UAS-produced snow depth maps resulted in a mean RMSE of 9.79 cm.

Figure 2 UAS snow depth maps, showing spatial variation in snow depth (cm) for each survey on 1–17 May.

The mean SD declined by 24.81 cm between the first survey on 1 May and the second survey on 14 May, but at that time, over 95 % of the area remained snow-covered (Table 2). Between 1 and 14 May, snow depth change (ΔSD) (Fig. S4) was higher in the most open areas and south-facing slopes compared to north-facing slopes and the northeastern areas with more tree cover. In addition, the snowmobile tracks crossing the fen showed a slower decline in snow depth compared to the surrounding areas. Considering only the daily surveys between 14 and 17 May, the greatest melting occurred between 14 and 15 May, with a 9.06 cm change in mean SD and 322.29 m³ change in SWE volume. However, the largest change in the snow-covered area occurred between 15 and 16 May, when SCA decreased by 40.71 %, and the largest ΔSD occurred in the central areas of the fen and the area adjacent to the natural spring. Between the last two surveys on 16 and 17 May, greater change occurred in the snow-covered area, which decreased by 17.9 %, with the most substantial changes in SD observed in areas with the greatest snow depth at the beginning of the study period, particularly near the borders and northern and northeastern parts of the fen. At this stage, changes in snow depth and volume were minor compared to previous days.

Table 2 Snow-covered area (SCA), mean and median snow depth (SD), volume of snow water equivalent (SWE), mean and median SWE for each UAS survey, and manual SD, SWE and snow density measurements. Snow density measurement location on 16 and 17 May was changed due to the absence of snow cover at the initial sampling point.

We used the SAGA wetness index (SWI) for the snow-free areas (i.e., areas with SD ≤ 0 cm) (Fig. 3a) to evaluate the interaction between the study site topography and wetness and snowmelt processes. Increasing values were noted throughout the surveys, with mean SWI of the snow-free areas increasing from 7.27 to 9.08 between 1 and 17 May, while the mean SWI for the whole study area was 9.14. This reflects melting progressing from drier to more wet areas. This pattern was also seen by differentiating snow depth by each SWI class (Fig. 3b), with areas classified as having high SWI exhibiting greater snow depth compared to the low or medium SWI classes until 15 May. However, by 16 and 17 May, snow depth in high SWI zones had decreased to levels comparable to or lower than those in the other classes. Based on the pairwise Wilcoxon test, differences in snow depth between SWI classes were statistically significant (p < 0.01), except for snow depth in low and medium classes on 16 May. Between 15 and 16 May, the largest decline in the snow-covered area occurred within the high SWI class (Fig. 3c), where SCA decreased by 49.97 %, while the change in the low and medium classes was 35.44 % and 47.44 %, respectively. Between 16 and 17 May, mean SWI of snow-free areas exhibited a slight increase, from 9.04 to 9.08, and the largest change in SCA occurred in the medium SWI class, where SCA decreased by 16.26 %. High SWI zones had the most persistent snow cover, and by the last survey on 17 May, 17.82 % of high SWI class cells remained snow-covered, while the SCA in low and medium classes was 12.96 % and 13.84 %, respectively.

Figure 3 Map of SAGA Wetness index (SWI) classes for melted areas (snow depth ≤ 0 cm), where areas with snow cover are shown in transparent (a), UAS snow depth (cm) (b), and snow-covered area (%) (c) in SWI classes (low, medium, high) for each survey. In (b), outliers or any values below 0 were considered errors and are not shown.

3.2 High-frequency hydrological and dissolved organic carbon (DOC) time series

The earliest snowmelt events in the study period occurred in mid-April 2024 (Fig. 4). During this time, two snowmelt peaks were recorded in stream flow, both triggered by a rise in air temperature (max temperature 10 °C) and rain events (8.7 mm in total). The next melting event was recorded in early May and was similarly initiated by a rise in air temperature (max temperature 7.8 °C), followed by a rain event (5.2 mm). Before the second melting event in May, low air temperatures (mean temperature 0.1 °C between 6 and 12 May) restrained the melting. The final large snowmelt peak occurred mid-May when rapid melting during the 5 d period caused snow depth to drop from 40 to 0 cm (Fig. 4b). Melting during this period was mainly driven by high air temperatures (mean 7.9 °C between 12 and 18 May).

Figure 4 Water quality and hydrological parameters monitored during the study period April–May. Discharge (Q) and DOC concentration monitored at the stream gauging station (a), air temperature and automatic point snow depth measurement in Puukkosuo (b), precipitation monitored in Puukkosuo and δ¹⁸O isotope samples of the stream, precipitation, and snow (c), DOC load and TSSeq concentration (d), and water table depth (WTD) and water temperature monitored in the GW well in Puukkosuo (e). Dashed lines represent the dates of UAS surveys.

The mean DOC concentration measured at the Puukkosuo outlet stream during the study period was 5.5 ± 0.7 mg L⁻¹. DOC, and TSS concentrations showed a rapid increase during the first recorded snowmelt pulse, which was also the first identified event (Figs. 4a, d, and 6b), during which the highest DOC concentration in the study period (11.3 mg L⁻¹) was recorded. In the following snowmelt events, increases in DOC concentration were smaller and coupled with discharge, causing DOC to display a staircase-like pattern, with each melting event resulting in a slightly elevated base DOC concentration. At the onset of the final snowmelt peak (events 7–12 in Fig. 6) in mid-May, DOC concentrations first increased simultaneously with discharge, reaching a maximum of 6.6 mg L⁻¹ on 13 May or event 8, but on the following peak discharge days, a diluting pattern was observed where DOC decreased with increasing discharge, and the lowest concentration of 4.8 mg L⁻¹ was measured on 15 May (event 10). However, on 16 May (event 11), when the highest discharge (5.2 L s⁻¹) in the study period was recorded, DOC concentration did not display a clear diluting pattern, reaching 6.8 mg L⁻¹ until the discharge increased again on the next day (event 12). TSS concentration remained more stable before the last snowmelt peak, when the concentrations showed a simultaneous increase with discharge.

DOC load also increased towards the end of the snowmelt period, with the highest load (max 118.2 g h⁻¹) observed at the last day of the final snowmelt peak (event 12) (Fig. 4d). Variation in discharge was considerably higher compared to only subtle changes in DOC concentrations, and thus, variation in DOC load was mainly driven by discharge dynamics. This was represented in the cumulative DOC load and discharge during the study period, which exhibited a similar staircase pattern, with each melting event causing an increase in both DOC load and discharge (Fig. 6a). The highest increase in both was observed during the final melt peak, when DOC load increased in parallel with discharge, and the rate of increase slowed after the snowmelt was completed.

Peatland groundwater level response in the GW well was very synchronized with streamflow and precipitation until the snowmelt was exhausted (Fig. 4e). WTD ranged from −4.5 to 5.4 cm, but the water table remained above ground level for most of the study period. Based on the field observations, the GW well was frozen during the study period; thus, high WTD is likely due to meltwater accumulating on top of the ice layer. WTD started to decrease during the final snowmelt peak and reached below the ground surface on May 17, although GW temperatures did not begin to rise until a week later.

Stream δ¹⁸O values in the study period varied between −17.22 ‰ and −13.89 ‰ (Fig. 4c). The mean δ¹⁸O value for precipitation during the study period was −12.48 ‰, with more enriched values occurring during warm periods, i.e., when precipitation occurred as rain, not snow. Snowpack δ¹⁸O values were generally more depleted when compared to stream and precipitation values, ranging from −19.39 ‰ to −16.70 ‰. Each snowmelt pulse was measured as a more depleted δ¹⁸O in stream water, verifying dynamics of snowmelt pulses in the stream. Hydrograph separation revealed that the average daily snowmelt water contribution (SM %) during the peak melt period (1–20 May) was 32 %. The highest contributions were observed on 15 and 16 May, when over 89 % of the streamflow was estimated to originate from meltwater (Fig. 5). This corresponded with the most depleted stream δ¹⁸O values recorded during the last snowmelt peak between 14 May and 17 May, when the mean δ¹⁸O was −16.48 ‰. After May 16, when over 75 % of the fen was free of snow cover and snow depth declined only slightly compared to previous days, SM % began to decrease, indicating increasing GW influence. After the snowmelt was complete, stream δ¹⁸O became less depleted and stabilized quickly, resembling the values before the initial snowmelt pulse in early April (Fig. 4c).

Figure 5 The estimated daily meltwater contribution (SM %) during the peak melt period (1–20 May) based on isotope hydrograph separation, mean snow depth (cm), and snow-free area (%) obtained from the UAS surveys.

3.3 Event analysis

Based on the event detection criteria, 12 events were identified during the study period, of which three occurred during the initial melting in early April, three in early May, and the last six during the final melt peak in mid-May (Table 3, Fig. 6a). The average event duration was 23 h. The longest event, number 6, lasted for 60 h, while event number 11, occurring on 16 May during the peak melt period, was the shortest with a 10 h duration. The highest mean DOC concentration was measured in event 12, while event 1 had the highest maximum DOC concentration. However, differences in DOC concentrations between events were small, and greater disparities were observed in DOC loads. When considering event duration, event 12, occurring on 17 May, contributed the largest hourly DOC load, while event 1 had the lowest hourly DOC load. The highest total DOC exports occurred in the longest event, event 6, followed by event 12, while event 4 contributed the lowest DOC exports.

Figure 6 Cumulative DOC load and discharge during the study period (a), events 1–12 identified based on stream discharge (b), quadrant plot of hysteresis index (HI) and flushing index (FI) with numbers representing the event ID (c), and a scatterplot of TSSeq and DOC concentrations with point colors representing the event ID (d).

Table 3 Timing, duration, maximum discharge, hysteresis index (HI), flushing index (FI), minimum, maximum, and mean DOC concentration, hourly DOC load, and total DOC load for events identified during the study period.

Normalized hysteresis loops (Fig. S5) exhibited complex behavior, with the majority being anticlockwise, two clockwise, and three figure-of-eight or otherwise complex-shaped loops. Mean HI was −0.25, ranging from −0.84 in event 7 to 0.6 in event 12. This reflects a dominant pattern of DOC concentrations peaking on the falling limb of the hydrograph. Event 12 was the only one to show a clear clockwise hysteretic pattern. Events 8 and 9 had weak positive HI very close to 0, indicating a lack of hysteretic pattern for these events. Mean FI was 0.15, ranging from −1 to 0.79. FI was mainly positive but underwent a shift to negative in events 9, 10, 11, and 12, indicating a transition from flushing of DOC during the rising limb to a diluting pattern during the events at the peak melt period with increasing daily maximum discharges (Fig. 6b). A quadrant plot of HI versus FI (Fig. 6c) revealed that most of the events represented an anticlockwise hysteresis together with flushing behavior, clustering in the bottom-left corner. However, events 7–12 in mid-May exhibited distinctive flushing and hysteresis patterns, shifting to diluting behavior during the final melt period and from anticlockwise to clockwise between subsequent events.

The Spearman correlation was used to identify potential links between DOC and TSS exported during events (Fig. 6d), and the mean SWI of the snow-free areas and event characteristics (Fig. 7). When examining each event individually, a negative correlation was observed in event 12 (cor = −0.71, p < 0.01), while no significant correlation between DOC and TSS was found in events 2, 7, 9, and 11. Events 1 and 4 had the strongest correlations (cor = 0.9, p < 0.01), indicating a strong link between the DOC and TSS. However, events' mean or maximum DOC concentrations did not seem to have a relationship with any other event characteristics. The mean SWI of the snow-free areas was used as a proxy for the potential hydrological connectivity and correlated with event characteristics. We found that mean SWI was significantly (p < 0.05) positively correlated with SM %, hourly DOC load, and maximum discharge, and negatively with FI (Fig. S6). Also, HI showed a positive, although nonsignificant (p = 0.09) correlation with mean SWI.

Figure 7 Correlation matrix for mean SAGA wetness index (SWI), maximum discharge (Max Q), hysteresis index (HI), flushing index (FI), mean and maximum DOC concentration, hourly DOC load, and daily meltwater contribution (SM %). Significant correlations are flagged, ^* for p < 0.05 and ${}^{***}$ for p < 0.001.

4 Discussion

4.1 Snow cover loss initiates hydrological connectivity and DOC leaching

Our first objective was to document spatiotemporal variations in snow cover melt using high-resolution UAS surveys. In Puukkosuo, spatial variation in snowmelt was initially controlled by environmental/physical factors. The first areas to become snow-free were south-facing slopes and areas with higher tree cover, particularly in the surrounding areas of the fen, which also had lower snow accumulation during the first survey. This pattern can be explained by higher energy received from solar radiation, the formation of tree wells and, to a minor extent, snow interception by trees (Grünewald et al., 2010; Harder et al., 2020). This was also reflected in lower SAGA wetness index (SWI) values for snow-free areas during the first three surveys, and lower snow depths in the low SWI class (Fig. 3a), indicative of drier conditions. Although increasing canopy cover has typically been associated with higher snow depths and delayed melt (Ketcheson et al., 2012; Meriö et al., 2023), areas with higher tree cover showed faster depletion of snow cover in our study. This is likely due to the scale of our study, covering only the fen and its bordering transitional woodland at the very end of the snowmelt period.

During the last two surveys, large snow-free areas first developed in the central parts of the fen, before continuing towards the borders, where the most snow had accumulated. Higher snow accumulation in these areas can be explained by wind distribution processes relocating snow from the open peatland area towards the sheltered edges (Harder et al., 2020; Meriö et al., 2023). Our observations are supported by the findings of Grünewald et al. (2010), who concluded that spatial variation in maximum snow accumulation has a greater impact on snow depletion than variability in snow melt rates. The development of snow-free patches and the subsequent transport of melt energy to adjacent areas are the main controlling factors of snow melt spatial variability (Grünewald et al., 2010). This pattern could be identified in our study area as well, as after the large snow-free areas formed on 16 May, differences in snow depth and snow depth change between the SWI classes decreased, reflecting more uniform melting. However, hydrology had some influence on snow accumulation and/or melting in our study area, as areas with higher wetness identified by SWI classes seemed to have more persistent snow cover and higher snow depths until the very end of the snowmelt. This can be explained by the formation of ice cover in the microtopographical depressions where the water table is above the ground surface, which was also observed in the field. With the UAS snow depth mapping technique, this surficial ice cover is captured as snow (Rauhala et al., 2023). More persistent snow cover of high SWI zones can also be explained by topographical depressions accumulating more snow due to wind-driven snow trapping (Shirley et al., 2025), which could potentially further increase the wetness. In addition, the area immediate adjacent to the spring (Fig. 2) became snow-free earlier in comparison to the surrounding areas, likely due to warmer groundwater enhancing the melt (Van Huizen et al., 2022). However, the impact is very small-scale, and the physical processes described above are likely the main contributors to spatial melting patterns.

The results demonstrated that the decline in snow depth or SWE volume and snow-covered area (SCA) followed different patterns. Snow depth or SWE volume decreased more gradually, with a rate of decline slowing towards the end of the peak melt period, while SCA showed a sharp decline at the end of the study period. The interplay of these variables has important implications, as snow mass defines how much snowmelt runoff can be produced, but the absence or presence of snow cover shapes soil thermal regime, soil moisture and flow paths (Dong, 2018; Webb et al., 2018). This was illustrated by visually comparing the daily meltwater contribution to streamflow with the evolution of snow depth and snow-free area (Fig. 5). First, large decreases in snow depth result in an increase in meltwater contribution, whereas meltwater contribution declines sharply once most of the surface becomes snow-free, indicating meltwater infiltration and mixing with soil water as the soil began to thaw (Eskelinen et al., 2016; Muhic et al., 2023). To explicitly monitor such patterns, UAS-based snow monitoring approach has the potential to address the coverage and resolution issues faced with point snow measurements and satellite data, and improve monitoring of snow cover at local and catchment scales during the period of rapid snowmelt.

The peak melting documented through UAS surveys proceeded rapidly due to high daily temperatures (mean 7.9 °C) and was potentially further accelerated by rainfall (2 mm) on the evening of 15 May (Fig. 4). During the peak melt period, an increase in mean SWI, reflecting both decreasing snow cover and higher potential connectivity, was linked to higher discharge peaks, snowmelt water contribution, and hourly DOC load. This states that discharge dynamics and increasing connectivity are the dominant drivers of the DOC load, which was also seen as a clear response in cumulative load and discharge to the melting events, and more constant DOC export between the events, resulting in an almost linear increase (Fig. 6a). However, the hourly DOC load increased almost linearly in relation to mean SWI, whereas maximum discharge and meltwater contribution showed a saturating nonlinear relationship. This indicates that once hydrological connectivity in the peatland was established, further DOC export was governed by increasing activation of DOC-rich flow paths rather than by runoff magnitude. This was also supported by the observation that on 16 May, when the first large uniform snow-free areas formed and melting occurred in the wettest areas, DOC concentration remained almost invariant despite the highest maximum discharge and over 89 % meltwater contribution. Exposure of the peat surface, particularly in high wetness areas, can be associated with increased hydrological connectivity, activating DOC-transporting source areas and soil layers (Knapp et al., 2022; Prijac et al., 2023), which might have offset the diluting effect of meltwater and resulted in relatively stable DOC concentrations even during peak discharge.

Subsurface hydrological connectivity during snowmelt is restricted by soil frost or ice layers at the top of the peat, forcing water to be transported as overland flow or as deeper preferential flow paths (Laudon et al., 2007; Peralta-Tapia et al., 2015). This process can lead to a dilution effect in high flows, as the streamflow consists of snowmelt water with low DOC concentrations (Eskelinen et al., 2016; Leach et al., 2016). This pattern was also observed in our study during the peak melt period with high daily maximum discharges (Fig. 4). However, when the soil begins to thaw and meltwater is able to infiltrate, new flow paths and pools of DOC can activate for transport (Croghan et al., 2023; Rose et al., 2023). Thus, depletion of snow cover likely facilitated DOC transport by allowing the thawing of surficial peat layers to begin as the insulating snow cover was removed. This was supported by the simultaneous decrease in WTD once large snow-free areas were formed, indicating soil thawing. Interestingly, the decrease of WTD below the ground surface after the final melt peak did not seem to have a clear impact on DOC concentrations. This could be explained by the consistently high water table, which remained less than 5 cm below the surface until the end of the study period. These surficial peat layers can be efficiently flushed during the snowmelt (Birkel et al., 2017; Rosset et al., 2020) or can potentially have limited DOC due to the freeze-out effect (Ågren et al., 2012). Typically, an increase in water table is presumed to reflect hydrological connectivity (Knapp et al., 2022; Prijac et al., 2023; Rosset et al., 2020), but in Puukkosuo, the water table remains high and fluctuations are small, and we expect the soil frost to be the main factor restricting subsurface hydrological connectivity during spring snowmelt.

The depletion of snow cover and subsequent thawing of peat layers could also boost biological activity and DOC production due to warming soil temperatures (Campbell and Laudon, 2019; Wen et al., 2020). However, production rates of DOC during cold periods are typically restricted by low temperatures, and it is unlikely that the production rate during spring snowmelt could compensate for DOC exports (Campbell et al., 2014; Wen et al., 2020), often resulting in source limitation (Gómez-Gener et al., 2021). Thus, prevailing conditions during the antecedent winter and the previous year's DOC exports largely define the amount of DOC accumulated and available for transport during snowmelt (Ågren et al., 2010, 2012; Tiwari et al., 2018). High DOC export during the preceding summer and autumn can decrease the DOC pool available for transport during the next spring unless production is sufficient to compensate (Ågren et al., 2010). Cold autumns and winters, associated with severe soil freezing, have been found to cause higher DOC peaks during snowmelt (Tiwari et al., 2018), as severe soil freezing and recurring freeze-thaw cycles during winter can increase peat DOC concentrations and lead to delayed leaching of more labile DOC due to physical disruption of soil aggregates (Campbell et al., 2014; Liu et al., 2022; Yu et al., 2011). However, it is the timing and degree of lateral hydrological connectivity that eventually defines when and whether DOC is transported (Knapp et al., 2022; Prijac et al., 2023; Wen et al., 2020).

On a catchment scale, initial activation of new DOC flow paths during snowmelt has been linked to the completion of melting and soil thawing in open peatland areas (Croghan et al., 2023), as these areas are typically first to melt (Meriö et al., 2023). Still, it is possible that part of the exported DOC in our study area originated from the upper forested catchment, especially at the end of the snowmelt period (Ågren et al., 2008, 2012; Laudon et al., 2011; Tang et al., 2018). Although we did not have data on the progress of snowmelt in the upslope forested area, the research station has a forest measurement site with point snow depth monitoring located approximately 600 m downhill from the Puukkosuo study site. Average snow depth in this forest site was 28 cm on 14 May and 12 cm on 17 May, which is on average 6.7 cm higher than the mean snow depth obtained from UAS surveys. Corresponding snow depths could be assumed in the forested upslope catchment of the study area. Delayed melting in the upper catchment can enhance hillslope connectivity and increase the amount of DOC transported from the forested upslope areas during the final phases of snowmelt (Croghan et al., 2023; Laudon et al., 2011). At the same time, increasing overland flow at the peatland can reduce the contribution of DOC transported from peat layers (Ågren et al., 2008; Laudon et al., 2011). Still, peatland can be expected to have a significant contribution to DOC exports, given that 24 % of the catchment area is covered by peatland and the stream gauging station is located immediately at the fen outlet, forcing all incoming water to pass through the organic layers (Laudon et al., 2011; Rosset et al., 2019). Although the gauging station captures runoff draining from the fen, some flow bypassing is likely due to multiple outlet pathways. We assume that this bypassing flow does not differ substantially in volume or water quality from the measured flow, however, this remains a source of uncertainty.

4.2 Event analysis reveals DOC transport dynamics and flow path activation during snowmelt

The measured DOC concentrations (mean 5.5 mg L⁻¹) in Puukkosuo were overall slightly lower compared to other studies carried out in boreal peatland catchments during snowmelt (Dyson et al., 2011; Gómez-Gener et al., 2021; Leach et al., 2016) but similar to concentrations measured in the subarctic Pallas catchment (Croghan et al., 2023, 2024). DOC concentration in the peatland outlet stream was highly influenced by the onset of snowmelt and the related increases in snowmelt water inputs. Initially, DOC concentration and discharge were low, but the first snowmelt pulse facilitated a peak in DOC and TSS, which was considerably higher than the events occurring during the rest of the study period (Fig. 4). A similarly high initial pulse in turbidity was observed by Croghan et al. (2023), who suggested it reflected DOC and particulate organic carbon (POC) mobilization from sedimented layers. After the first event, increases in DOC with discharge were more subtle, but sources remained sufficient to increase DOC concurrently with discharge until the final melting peak, suggesting the increasing activation of flow paths and transport of DOC as melting proceeds (Marttila et al., 2021). While most events exhibited positive correlations for TSS and DOC, the absence of correlation in event 11 and the negative correlation in event 12 suggest a decoupling of particulate material and DOC sources or contrasting transport mechanisms during late melt stages.

Analysing events' concentration-discharge responses enabled identifying the relevant DOC transport processes. Most of the identified events during the study period exhibited anticlockwise hysteresis together with positive flushing indexes (FI) (Fig. 8c), which has been attributed to rapid flushing of DOC and reconnection of carbon-rich subsurface flow paths later in the event in peatland ecosystems, causing concentrations to peak in the falling limb of the event (Prijac et al., 2023; Tunaley et al., 2016). The dominant C–Q patterns observed reflect predominantly transport-limited DOC export, yet surprisingly small dilution effects, even with high contributions of meltwater, suggest the system being nearly chemostatic (Shatilla et al., 2023; Wen et al., 2020). Previous studies in peatland catchments have found both anticlockwise (Burd et al., 2018; Prijac et al., 2023; Tunaley et al., 2016) and clockwise (Ågren et al., 2008; Dyson et al., 2011) hysteretic patterns, and often a combination of both in spring (Croghan et al., 2023, 2024; Rose et al., 2023), highlighting the complex transport dynamics and sources during the snowmelt season (Shatilla et al., 2023).

Figure 8 Conceptual diagram summarizing the observed and hypothesized snowmelt-driven DOC transport processes during the snow-covered and final melting period in Puukkosuo fen. Blue arrows represent the main flow paths contributing to DOC transport.

We observed a change in hysteretic behavior during the final melting peak contrasting the dominant anticlockwise pattern, as the hysteresis index (HI) shifted first to slightly positive but was close to zero, then back to negative, and reached the highest value during the last detected event. A similar shift from anticlockwise to slightly clockwise hysteresis during snowmelt was observed by Croghan et al. (2023), who attributed it to a change from slower acting subsurface flow paths to the activation of dominant DOC transporting pathways or the depletion of DOC sources in the catchment. We also found changing flushing behavior during the last snowmelt peak. A shift to negative FI was observed in events 9, 10, 11, and 12 when daily discharge increased, which can be attributed to a diluting effect of snowmelt water when DOC sources are inadequate to keep up with increasing runoff (Gómez-Gener et al., 2021). However, particularly in event 11, almost negligible changes DOC concentration indicate that a fraction of runoff continued to interact with DOC-rich flow paths, rather than depletion of DOC sources or flow paths bypassing the peat entirely. The shift to negative HI in event 11 suggests progressive activation of DOC-transporting pathways later in the event (Prijac et al., 2023). The activation of these hydrological pathways likely primed the rapid transport of DOC for the final event, event 12, (Croghan et al., 2023), when clockwise hysteresis was observed and DOC peaked on the rising limb of the event. However, slightly stronger dilution in the final event might have resulted from depletion of the newly activated sources and/or high water table facilitating overland flow (Croghan et al., 2023; Leach et al., 2016).

By integrating high-resolution stream monitoring and spatial snow cover mapping, we were able to link the observed DOC dynamics to the snowmelt-driven DOC transport processes (Fig. 8). The event analysis results suggest a shift in dominant sources or flow paths as the snowmelt progressed, as early melting events exhibited hysteretic and flushing patterns typical for lateral subsurface flow but shifted to more surficial flow path behavior during peak melting (Croghan et al., 2023; Prijac et al., 2023). This is consistent with previous studies showing that most of the exported DOC during spring originates from deeper flow paths, but shifts to surficial carbon pools as melting proceeds (Leach et al., 2016; Rose et al., 2023). With the snow surveys, this shift can be attributed to the depleting snow cover. However, during the final melt peak, the relative stability of DOC concentrations despite the high contribution of meltwater reflects changes in transport efficiency and hydrological connectivity rather than major changes in sources. This is consistent with the observed relationships between mean SWI of snow-free areas and event metrics during peak melt, where higher SWI was associated with lower FI and a weakly positive relationship with HI, reflecting increasing, although subtle, dilution alongside the activation of quick, well-connected flow paths (Blaen et al., 2017; Vaughan et al., 2017). The observed hysteretic behavior during snowmelt is complex and varies between subsequent events, which may be linked to depletion of accumulated DOC sources (Shatilla et al., 2023) or, more likely in our case, to spatially and temporally heterogeneous thaw conditions (Rose et al., 2023). As indicated by the snow depth maps, snow cover, and by inference soil frost, is not spatially uniform, which increases the flow path complexity and connectivity as melt progresses.

The role of catchment conditions and spatial variation in hydrological connectivity has often been overlooked in concentration-discharge analysis, and hysteretic relations have been expected to represent the activation of different source areas (Knapp et al., 2022). The integration of novel monitoring technologies holds potential to bridge this gap by explicitly linking spatial and temporal patterns. Although the topographically derived wetness index is not able to capture complex peat-specific flow paths, it can provide a reasonable estimate to describe the spatial distribution of areas that are more likely to contribute to hydrological connectivity, especially during periods of high water table when stronger topographic control of lateral flow can be expected (Richardson et al., 2012). In peatlands, preferential flow paths may represent an important pathway, especially in the presence of soil frost (Laudon et al., 2007; Peralta-Tapia et al., 2015). Based on the hydraulic conductivity measurements done at the study site, a layer with higher conductivity exists at approximately 30–50 cm depth, which could act as a preferential flow path and support DOC supply during the snowmelt. Although C–Q analysis can provide insights into DOC transporting flow pathways (Lloyd et al., 2016), an integrated analysis of spatial and vertical connectivity gradients could provide a more detailed understanding of subsurface flow routing. In addition, further information on variation in pore water DOC concentrations is needed to provide a better assessment of spatial variability, mobilization, and transport processes (Knapp et al., 2022; Vitt et al., 2022).

5 Conclusions

This study used a novel combination of high-resolution UAS snow depth mapping, topographical wetness index, and high-frequency stream monitoring to reveal spatiotemporal snowmelt dynamics, hydrological connectivity, and DOC transport in the northern Puukkosuo fen. By coupling spatial and temporal high-resolution monitoring and multiple analytical techniques, we found that snow cover depletion, particularly in high wetness areas, triggers hydrological connectivity and DOC transport. Event analysis also suggested the activation of new sources over time, likely due to the thawing of surficial peat layers. Overall, shifting hysteresis and flushing behavior were observed during the study period, reflecting changes in dominant flow paths from deeper sources early in the melt to more surficial contributions as the snowmelt period progressed and snow cover started to deplete. However, the stability of DOC concentrations regardless of high meltwater contributions to streamflow suggests nearly chemostatic behaviour and continuous supply of DOC.

The combination of high-resolution temporal stream monitoring and spatial snow cover mapping allowed detailed insights into the connection between snowmelt progress and stream Q–DOC response. As shown by our results, these processes exhibit rapid variation even on a daily timescale, emphasizing the need for high-resolution monitoring to adequately capture the processes in play and estimate DOC export during snowmelt. UAS snow monitoring proved to be a cost-efficient method to acquire spatially and temporally resolved information on snow cover changes during the snowmelt period. UAS-based information on snow depth, volume and snow-covered area could also provide more explicit information on these variables valuable for ecohydrological and biogeochemical modelling and its development. To our knowledge, this study is the first to use UAS surveys to document and link spatial variation in snowmelt-induced hydrological activation to event Q–DOC characteristics. A similar UAS-based monitoring method could provide insights into spatial variation in hydrological connectivity and provide a more robust link between the catchment conditions and events C–Q dynamics during other seasons as well. Introducing spatial mapping is needed to address the spatial hydrological processes associated with DOC transport within peatlands, and we encourage further studies to combine spatially resolved information on catchment hydrological conditions with event C–Q analysis to provide more robust interpretations on activation and contribution of different source areas. To further improve the understanding of the dynamic nature of peatland ecosystems, more information on how different parts of peatland contribute to DOC transport and production is needed.

Our identified processes can be expected to be relevant also in other fen ecosystems within northern boreal regions with comparable seasonal and snowmelt patterns. Our results indicate that hydrological activation is sensitive to late winter processes, as the fen is already activated in the early phase of melting. Thus, earlier or more rapid snowmelt and warmer periods in winter can directly impact hydrological activation and DOC leaching. Also, DOC transport processes may be impacted by any changes in snow, soil frost, and thawing conditions, thus Puukkosuo and similar fens in boreal and arctic regions are likely to be sensitive to climate change impacts. These findings emphasize the importance of snowmelt timing and spatial patterns in controlling peatland DOC export and suggests that water-carbon dynamics may undergo substantial changes under future climate conditions.