Qualitative transformations of dissolved organic matter along supra-permafrost flow: insights from Arctic subterranean estuaries

Increasing coastal erosion and permafrost thaw along the Arctic shoreline represent major lateral sources of dissolved organic matter (DOM) to the Arctic coastal ocean. One key and still underrecognized flow path is non-point source supra-permafrost groundwater flow discharging through sandy beaches in front of coastal bluffs. Here, we explore the qualitative transformations of DOM and dissolved organic carbon (DOC) along this flow path, with a focus on the beach discharge zone, where fresh supra-permafrost groundwater mixes with recirculated seawater. We sampled meltwater, beach groundwater, and seawater along coastal bluff transects extending up to 0.5–1 km offshore. Along these salinity and redox gradients, DOC and chromophoric DOM (CDOM) declined sharply, partially reflecting dilution. Optical indices (a350, SUVA₂₅₄, HIX) and PARAFAC-based fluorescence analyses revealed a shift from humic-like, high molecular weight (HMW) DOM in meltwaters to protein-like, low molecular weight (LMW) DOM in nearshore waters. The inverse relationship between humic- and protein-like components, combined with reactive Fe-hydroxide coatings on sandy sediments, elevated dissolved iron (Fe_tot) and high dissolved inorganic carbon (DIC) under low oxygen conditions in beach groundwater, suggests that microbial degradation and mineral-organic interactions likely contribute to DOM transformation. Alongside these transformation processes, mixing governs bulk DOC behaviour. These findings suggest that nearshore and intertidal zones, particularly subterranean estuaries (STEs), are zones of DOM transformation along the Arctic coastline. However, whether they act primarily as transient filters, long-term carbon sinks, and/or dynamic biogeochemical reactors remains unresolved, highlighting the need for further research on their role in transferring permafrost-derived carbon to the ocean.

1 Introduction

Permafrost stores around 1300 Pg of organic carbon (OC) within its 13.9 × 10⁶ km² surface area, which represents 60 % of the world's soil organic carbon stored in 15 % of global soil area (Schuur et al., 2015; Obu et al., 2019) . The Arctic permafrost coastline is heavily impacted by global changes, resulting in unprecedented thawing rates and the deepening of the active layer, which, in turn, increases subsurface transport (Lantuit et al., 2012; Jones et al., 2020). Unlithified and ice-bonded permafrost cliffs, such as those along the Beaufort Sea, are highly susceptible to coastal erosion. Over the past twenty years, these cliffs have experienced one of the highest coastal erosion rates recorded in the Arctic, with a recorded rate of 1.1 m yr⁻¹ between 1950 and 2000. This rate has increased by 80 %–160 % in the last two decades (Lantuit et al., 2012; Jones et al., 2020). Accelerating coastline erosion is supplying increasing quantities of terrestrial materials (Kipp et al., 2018), associated nutrients (e.g., nitrogen and phosphorus; Fritz et al., 2017) carbon (including organic carbon and methane; Bristol et al., 2021), and contaminants (such as heavy metals and organic pollutants; Kwasigroch et al., 2018) to the nearshore and coastal ocean. Recent studies indicate that land-derived nutrients from coastal erosion drive ∼ 30 % to 50 % of primary production in the Arctic Ocean, underscoring the critical biogeochemical role of permafrost-derived nutrient and carbon inputs on Arctic shelves (Terhaar et al., 2021; Vonk et al., 2025). Beyond erosion, subsurface pathways also contribute to these fluxes, as thaw depths and coastal hydraulic gradients peak in late summer and remobilize solutes from diffuse, non-point sources into the coastal ocean (Walvoord and Striegl, 2007; Lecher, 2017).

Dissolved organic matter (DOM) represents a fundamental link between terrestrial and aquatic carbon cycles, playing a significant role in the biogeochemistry of aquatic ecosystems (Hedges and Keil, 1995). At the ocean-basin scale, newly mobilized organic matter from Arctic watersheds may enhance emissions of carbon dioxide and other greenhouse gases (Vonk et al., 2015; Kaiser et al., 2017a; Tanski et al., 2019; Lapham et al., 2020; Vonk et al., 2025). It may also modify the productivity of primary producers (Thingstad et al., 2008; McMeans et al., 2015; Nguyen et al., 2022), with potential consequences for the Arctic marine food web and the communities that depend on these resources. This new input affects ocean chemistry by altering carbon and nutrient cycling (Guo et al., 2007; Stedmon et al., 2011; Vonk et al., 2014) and changes optical conditions, reducing light penetration and limiting primary production (Matsuoka et al., 2012; Fichot et al., 2013). A fraction of this terrestrial DOM can be rapidly mineralized through microbial and photochemical processes, affecting nutrient budgets, air-sea CO₂ exchanges, biological productivity, as well as acidification in coastal waters (Kaiser et al., 2017a, b). For example, Kaiser et al. (2017a) showed that ∼ 50 % of the annual dissolved organic carbon (DOC) discharged by Siberian rivers was mineralized in estuaries and on Eurasian continental shelves, indicating significant removal before reaching the open ocean. Therefore, only a small fraction of riverine-derived compounds potentially persists in the ocean over centuries and millennia (Fichot and Benner, 2014; Kaiser et al., 2017a). While the export of terrestrial dissolved organic carbon is known to significantly influence the Arctic marine ecosystem, with an estimated 25–38 Tg C yr⁻¹ discharged into the Arctic Ocean (Holmes et al., 2012), the role and significance of inputs from coastal erosion and permafrost thaw in shaping the biogeochemistry and ecology of nearshore Arctic waters remain poorly understood. This knowledge gap stems not only from the stochastic and episodic nature of erosion- and thaw-related fluxes, which complicates their quantification, but also from the complex interactions between mineral and organic phases that control the stabilization, degradation, transformation, and transport of the DOM pool along flow paths. The mechanisms of organic matter transformations occurring at mineral-organic interfaces are complex (see Li et al., 2023; and references therein), particularly in intertidal sandy sediment, where biological, geochemical, and redox conditions interplay to influence DOM concentrations and molecular compositions. In addition to dilution and bacterial degradation, a key process in these dynamic transitional environments is also the formation of iron (Fe) curtains at the freshwater-seawater redox interface, which may serve as major zones for the storage and transformation of terrigenous organic carbon (Chen et al., 2014; Linkhorst et al., 2017; Sirois et al., 2018; Zhou et al., 2023, 2024; Amoako et al., 2025). At the shoreline, beach subterranean estuaries (STEs), where terrestrial groundwater mixes with recirculating seawater (Moore, 1999), act as biogeochemical hotspots (Anschutz et al., 2009) with the Fe cycle playing a central role in both the remineralization and sequestration of organic carbon (Charette and Sholkovitz, 2002; Sirois et al., 2018; Chaillou et al., 2024). Redox oscillations driven by fluctuations in water levels (e.g., tide amplitude, waves, sea-level rise, and water table variations; Santos et al., 2012) shape the zones where reactive minerals such as Fe-hydroxides precipitate, forming transient or persistent geochemical barriers that influence the mobility of compounds, including DOM and DOC, and thereby modulate their export to adjacent coastal waters (Riedel et al., 2013; Linkhorst et al., 2017; Sirois et al., 2018; Waska and Banko-Kubis, 2024; Amoako et al., 2025). While most studies on beach STEs and their associated biogeochemical processes have focused on lower-latitude systems, our understanding of their occurrence and functioning in Arctic environments remains scarce. The remoteness of Arctic beaches, the logistical challenges of sampling, and the high associated costs limit our understanding of these systems. However, recent research in Alaskan coastal lagoons by Demir et al. (2024) has confirmed that supra-permafrost groundwater flows, originating from the seasonally thawed active layer, can extend beyond coastal bluffs and beaches, even in areas underlain by permafrost. The presence of unfrozen, permeable coastal sediments thus creates pathways for fresh groundwater flow that mixes with recirculating seawater (as in STEs; Moore, 1999) and discharges into coastal zones, as observed in the Laptev Sea (Charkin et al., 2017) and in the Beaufort Sea along the Alaskan and Mackenzie coasts (Bullock et al., 2024; Kipp et al., 2025).

Absorbance- and fluorescence-derived indices are widely used to characterize the origin, reactivity, and transformation of aquatic DOM (Stedmon et al., 2003; Fichot and Benner, 2014; Meilleur et al., 2023). For instance, Fouché et al. (2020) identified permafrost-derived DOM as being predominantly low molecular weight (LMW), proteinaceous, and low in aromaticity, while DOM from the seasonally thawed active layer is characterized by high molecular weight (HMW), aromatic, and condensed organic compounds. In low-latitude STE, selective retention and transport of organic molecules along the flow path led to a fractionation of DOM, thereby shaping the export of terrestrial-derived organic compounds to adjacent coastal waters (Couturier et al., 2016; Linkhorst et al., 2017; Hébert et al., 2022). The transformation and evolution of seasonally mobilized compounds at the land–ocean interface, particularly beyond beaches and coastal bluffs and along supra-permafrost flow paths, remains poorly constrained in Arctic coastal environments. This study aims to characterize the transformation pathways of DOM released from the thawing of coastal Arctic bluffs and to assess the role of sandy beaches as zones of terrigenous DOM transport and transformation, with a focus on compositional changes rather than flux quantification. A fine-scale sampling strategy (<1 km) was deployed in the Kugmallit Bay region (NWT, Canada) to measure the optical characteristics of DOM, including chromophoric dissolved organic matter (CDOM) and fluorescent dissolved organic matter (FDOM), as well as the behavior of DOC along the supra-permafrost flow, from coastal bluffs to the adjacent nearshore waters.

2 Materials and methods

2.1 Site Description

The study area is situated in the Inuvialuit Settlement Region of the Northwest Territories, near the Mackenzie Delta, the fourth-largest river system discharging into the Arctic Ocean (Macdonald et al., 1998). The Mackenzie River is a major source of freshwater, sediments, dissolved organic matter (DOM), and nutrients, exerting a strong influence on the salinity and biogeochemistry of the Beaufort Sea (Holmes et al., 2012). The tidal amplitude in Kugmallit Bay region is small (∼ 0.3 m). Two field campaigns were conducted: the first from 24 July to 6 August 2019, and the second from 22 June to 31 August 2021, coinciding with peak thaw conditions. Approximately 60 samples were collected from four coastal sites characterized by continuous permafrost and active thaw slumps, bordered by sandy and clay-rich beaches: Tuktoyaktuk Island, Peninsula Point, Crumbling Point, and Reindeer Island (Fig. 1).

Figure 1 Map of the four sampling sites (red dots) located in the Northwest Territories, Canada.

Tuktoyaktuk Island, the primary sampling site (N=34; Fig. 2A), features a coastal bluff roughly 9 m high, 1.5 km long, and 100 m wide (Ouellette, 2021). Located southeast of Kugmallit Bay and opposite the Hamlet of Tuktoyaktuk, the island is undergoing rapid coastal retreat, losing approximately 1.8 m of shoreline annually due to storm events and permafrost thaw. This erosion rate has increased by 22 % over the past 15 years (Berry et al., 2021; Whalen et al., 2022; Tanguy et al., 2023), and the entire site is projected to vanish within the next 20–30 years (Jones et al., 2020). In front of the bluff lies a 50 m-wide beach composed of fine to medium sands (0.3–0.5 m deep) overlying a frozen clay horizon (Fig. 2A). The exact submarine extension of this permafrost layer is unknown. Peninsula Point (N=7), southwest of Tuktoyaktuk Island within the Pingo Canadian Landmark, features a complex retrogressive thaw slump system underlain by a massive ice body ranging from 5 to 20 m in thickness (Mackay, 1986). Large, muddy lobes of thawed permafrost material and meltwater flow downslope toward the coast, intermittently covering the sandy intertidal zone (Fig. 2B). The shoreline has retreated at an average rate of ∼ 3.4 m yr⁻¹ between 1985 and 2018 (Hayes et al., 2022). Crumbling Point (N=5), located at the northwest edge of Kugmallit Bay, is another retrogressive thaw slump system with geomorphological characteristics similar to Peninsula Point. Reindeer Island (N=11), north of Richards Island, is part of a lagoonal complex formed by thermokarst lakes and surrounded by coastal bluffs. Although specific erosion rates for this site are not yet available, they are likely comparable to regional averages for the Canadian Beaufort Sea (∼ 0.5 m yr⁻¹; Solomon, 2005). However, localized retreat rates can be substantially higher as seen in nearby Pullen Island, where erosion exceeds 12 m yr⁻¹ (Berry et al., 2021). The width of the intertidal zone along these beach systems varies from a few meters to several tens of meters, depending on oceanographic conditions (e.g., wind regime and seasonal beach morphology).

Figure 2 Schematic cross-sections and photographs illustrating the coastal permafrost landscape at Tuktoyaktuk Island (A) and Peninsula Point (B; not to scale) as examples of a retrogressive thaw slump system. Black arrows show theoretical surface (solid lines) and subsurface (dashed lines) flow paths, while blue arrows represent recirculating seawater driven by tidal and wave inputs into the intertidal zone of sandy beaches. Blue dots in panel (B) mimic melt ponds. The locations of the different sample types (“massive ice”, “meltwater”, “beach groundwater”, and “seawater”) are indicated. These conceptual models are adapted from Kipp et al. (2025) based on Hayes et al. (2022) for Peninsula Point and Whalen et al. (2022) and Coote et al. (2025) for Tuktoyaktuk Island.

2.2 Water and Sediment Sampling

At each site, we carried out localized sampling along a transect, from coastal thawed bluffs, through beaches and intertidal zones, to offshore seawater as presented in Kipp et al. (2025). Meltwater and beach groundwater (here defined as porewater in sandy sediment) samples were collected from coastal permafrost slumps and the adjacent sandy shore, respectively, while seawater samples were collected 0.5 to 1.0 km offshore from the coastline at each study site. Surficial sandy sediments from beaches in front of permafrost bluffs were also collected during the 2019 field campaign at Tuktoyaktuk Island (N=3) and Peninsula Point (N=2). Meltwater was sampled directly in puddles formed on the thaw slumps using a submersible pump. Beach groundwater was collected using push-point piezometers inserted ∼ 50 cm deep into the sandy ground, above the frozen clay layer, within the intertidal zone, with continuous pumping via a Solinst^® peristaltic pump. A massive ice sample was collected from a permafrost core collected at Richards Island, which was sectioned and thawed in a hermetically sealed, acid-cleaned bucket before collecting the meltwater via a peristaltic pump. Finally, seawater samples were collected from a small vessel in front of the slump systems positioned ∼ 0.5–1.0 km offshore using a submersible pump placed 0.5–1.0 m below the surface. At each site, water samples were pumped into an online flow cell where salinity (S), temperature, pH (at NBS scale), and oxygen saturation were monitored with a daily calibrated multiparametric probe (600QS, YSI Inc.). Note that these stations were likely under the influence of the Mackenzie River plume, which mainly controls the water chemistry and optical signature in the region (see Lizotte et al., 2023, and references therein).

Water samples were collected for CDOM/FDOM into acid-washed 60 mL glass amber bottles after on-line filtration through a 0.22 µm Millipore Opticap^® XL4 cartridge with a Durapore^® membrane. The samples were stored in the dark at 4 °C. Total dissolved Fe (Fe_tot) samples were collected in 60 mL metal-free Falcon^® tubes after filtration through the same 0.22 µm Millipore Opticap cartridge. The samples were acidified with 3 drops of 70 % nitric acid to prevent the re-oxidation of reduced trace metals and stored at 4 °C. DOC samples were taken using 60 mL acid-cleaned polypropylene syringes and rapidly filtered with pre-combusted (450 °C for 5–6 h) 25 mm Whatman GF/F (0.7 µm) and stored in pre-combusted and acid-washed 12 mL borosilicate EPA tubes with PTFE caps. The same day, the DOC samples were acidified to pH < 2 with high-purity HCl 2N and stored in the dark at 4 °C until analysis. During the 2019 campaign, samples were also collected in 30 mL scintillation vials, hermetically sealed for further water stable isotope analysis.

2.3 Chemical and Optical Analysis

Total organic carbon (TOC) content was measured on freeze-dried, ground and homogenized sediment using a CHNSO analyzer (Costech ECS 4010 CHNSO analyzer^®; precision < 0.3 %). Reactive Fe-hydroxide content was determined on the same sediment samples, extracted using the citrate-dithionite-bicarbonate (CDB) reduction method of Mehra and Jackson (1960), as modified by Lalonde et al. (2012), and previously applied to sandy sediments by Sirois et al. (2018). Briefly, reactive Fe-hydroxides are reductively dissolved with dithionite at circumneutral pH (bicarbonate buffer) using citrate as a complexing agent to keep dissolved Fe in solution. Fe in solution was analyzed using a Microwave Plasma Atomic Emission Spectrophotometer (4200 MP-AES; Agilent Technologies). The detection limit of the method is 2.4 µg L⁻¹ for Fe concentration at a wavelength of 391 nm and analytical precision was <5 %.

Stable isotopes of water (δ¹⁸O, δ²H) were analyzed by EA-IRMS during the following year after the collection. Uncertainties are ±0.05 ‰ and ±1 ‰ for δ¹⁸O and δ²H, respectively. Reference materials were used throughout the isotopic water analyses, and isotopic analyses are reported compared to the international Vienna Standard Mean Ocean Water (VSMOW). DOC samples were analyzed a few weeks after data collection by Total Organic Carbon analyzer (TOC-VCPN Shimadzu) based on the method of Wurl and Tsai (2009). Trueness was verified against certified reference solutions (EMSL secondary standard and DSR Batch 21 Lot 04-21, University of Miami), with analytical uncertainty less than 4 % and a detection limit of 5.8 µmol L⁻¹. To ensure instrument stability, fresh acidified deionized water (blank) was regularly analyzed. The concentration of total dissolved iron (Fe_tot) was measured according to the ferrozine method proposed by Stookey (1970) and adapted by Viollier et al. (2000). The detection limit of the method was 0.4 µmol L⁻¹ and the reproducibility was better than 0.3 %.

Absorbance and fluorescence spectroscopy were used for the measurement of the chromophoric and fluorescent fractions of DOM (CDOM and FDOM) a few weeks after sampling. The CDOM absorbance was measured using a Lambda 850 UV-VIS Perkin Elmer spectrophotometer with 1 cm path-length quartz cuvettes. Measurements were taken from 220 to 800 nm at 1 nm intervals with a scanning speed of 100 nm min⁻¹ and a 5 nm slit width. Blanks and references were measured using fresh Milli-Q water. The FDOM was measured concomitantly using a Varian Cary Eclipse spectrofluorometer. Fluorescence spectra were measured within the emission wavelengths of 220 to 600 nm and the excitation wavelengths of 220 to 450 nm at 5 nm intervals as described by Couturier et al. (2016). Similarly, fresh Milli-Q water was used as a blank to rinse the cuvette in between samples. Fresh deionized water was used as a blank, and absorbance measurements were used to correct the inner-filter effect. Fluorescence spectra were further corrected for instrument-specific excitation and emission bias using manufacturer-provided correction files, and the dataset was corrected for Rayleigh and Raman scattering, according to the method described by Pucher et al. (2019).

2.4 Optical-derived indices and PARAFAC model

Absorbance and fluorescence indices were extracted using the staRdom toolbox on the R Studio Software (Pucher et al., 2019). Different indices were explored, including a350, BIX, HIX, FI, SUVA₂₅₄ and SR. However, only three were reported here to characterize the DOM pool, as they showed distinct variations among the different categories of samples (i.e. beach groundwater, meltwater, massive ice, and seawater). The spectral absorption coefficient at 350 nm (a350) was used as a proxy for CDOM concentration in the samples (Helms et al., 2008). It was calculated as 2.303 times the absorbance at the wavelength λ=350 nm divided by the path length of the cuvette (m). The specific UV absorbance (SUVA₂₅₄ in L mg C⁻¹ m⁻¹) was calculated as the absorbance at the wavelength λ=254 nm normalized to DOC concentration and is commonly used as an indicator of DOM aromaticity (Weishaar et al., 2003; Helms et al., 2008): greater SUVA₂₅₄ values correspond to a greater degree of aromaticity (Helms et al., 2008). It has also been shown to be positively correlated with the molecular weight of the DOM compounds. In the FDOM pool, the humification index (HIX) corresponds to the peak area under emission of 435–480 nm divided by the peak area under emission of 300–345 nm, at an excitation of 254 nm. HIX is an indicator of humic substances and the extent of humification of DOM compounds (Ohno, 2002; Hansen et al., 2016): higher HIX indicates a greater humification of the DOM source and HMW compounds.

In combination with absorbance and fluorescence indices, a PARAFAC model based on over 250 Arctic and subarctic coastal water samples was developed to further investigate the composition and sources of FDOM across samples (Bro, 1997; Stedmon et al., 2003; Murphy et al., 2013). Five components were validated using split-half analysis, following the method adapted from Pucher et al. (2019) in R, with the model explaining more than 94 % of the total variance in the dataset. The components were also matched with the literature for identification and external validation, using OpenFluor (Murphy et al., 2014). The five fluorescent components (C1–C5) identified are presented in Fig. 3, and their characteristics based on the literature are summarized in Table 1. Briefly, components C1 and C4 are associated with terrestrial humic-like DOM, with C1 representing the dominant terrestrial humic signal, while C4 contributes a weaker humic fraction. Component C5 corresponds to a minor microbial humic-like signal. In contrast, components C2 and C3 correspond to protein-like DOM, with C2 representing the dominant tyrosine-like signal and C3 a weaker tryptophan-like component.

Figure 3 Excitation-emission matrix (EEM) fluorescence spectra of the 5-component PARAFAC model. Fluorescence is expressed in Raman Unit (RU).

Table 1 EEM-PARAFAC components and their characteristics, based on literature assignments.

Pearson correlation coefficients were calculated between PARAFAC components, expressed as absolute fluorescence intensities (RU). Spearman rank correlations were used when relating PARAFAC components to ratio-based optical indices (HIX). Correlations were considered significant at p<0.001. Differences in geochemical parameters and PARAFAC component intensities across sample types were assessed using Kruskal-Wallis tests, followed by pairwise Wilcoxon tests with Bonferroni correction for post-hoc comparisons.

3 Results and discussion

3.1 Physico-chemical characteristics along the permafrost-nearshore continuum

In this study, we refer to samples collected at the nearshore as “seawater”; nevertheless, we acknowledge that these samples more accurately represent a brackish environment, with salinity S values ranging from 0.5 to 20, mostly related to the Mackenzie River plume occurrence. The massive ice and meltwater samples exhibited the lowest salinities with S<1, whereas the salinity of beach groundwater samples ranged between 0.2 and 5.4. The S range measured in beach groundwater samples mostly reflected the tidal and wave pumping effect and the recirculation of the infiltrated seawater within the permeable sediments (Santos et al., 2012). The higher salinities (S>5.4) were only measured in seawater samples. The temperature varied between 8.1 and 16.7 °C (with a mean value of 13.0 ± 2.6 °C), with the higher temperatures measured in meltwater and some beach groundwater samples. Oxygen saturations ranged from 6 % to 141 % of saturation. The nearshore surface seawater and the meltwater samples were all oversaturated due to mixing and photosynthetic processes. However, the low-salinity beach groundwater samples exhibited low oxygen saturation (6 %–48 %) despite the infiltration and recirculation of well-oxygenated seawater. Redox oscillations and transient oxygen-depleted conditions are observed in micro- and meso-tidal sandy STEs (Sirois et al., 2018; Waska et al., 2021; Hébert et al., 2022), where oxygen supplied through tidal exchange is rapidly consumed by heterotrophic microbial processes (Chaillou et al., 2018, 2024; Moore et al., 2024). The total dissolved iron (Fe_tot) concentrations ranged from 0 µmol L⁻¹ (or below detection limits) in well-oxygenated seawater samples to more than 1400 µmol L⁻¹ in a deoxygenated meltwater sample. In beach groundwater samples, most Fe_tot concentrations were on the order of a few µmol L⁻¹, but one sample reached values above 600 µmol L⁻¹, corresponding to lower oxygen saturation levels. Regarding pH (NBS scale), the highest median values were observed in both meltwater and seawater (7.9), while beach groundwater exhibited a slightly lower median pH of 7.5, ranging from 7.0 to 8.0, with one sample dropping to 6.5 (also associated to the highest Fe_tot concentrations of ∼ 600 µmol L⁻¹). Although pH data for massive ice samples were not available, their values likely fall within the lowest range of measured pH. This assumption is based on the origin of the massive ice in the region, which is thought to result from intrasedimental processes (e.g., ice intrusion and segregation) rather than from relict glacial ice (Mackay and Dallimore, 1992). Surficial beach sediments (<50 cm below the beach floor) are relatively poor in particulate organic carbon (TOC < 0.45 %) and mainly composed of clastic sediment coated with reactive iron-oxide (Table 2).

Table 2 Particulate fraction characteristics of sandy beach sediment collected at Tuktoyaktuk Island and Peninsula Point.

3.2 Origin of the subsurface water flow in the intertidal zone

Little is known about the hydro(geo)logical connectivity between supra-permafrost groundwater and Arctic coastal seawater (Lecher, 2017). However, recent studies and conceptual models, mostly developed for lower-latitude regions, suggest that the transport of water and permafrost-derived solutes to the coastal zone is governed by three main hydrological pathways (Walvoord and Kurylyk, 2016): (i) riverine transport, which includes both active layer runoff and contributions from deeper groundwater (Dittmar and Kattner, 2003; Demir et al., 2024); (ii) coastal erosion, which directly delivers organic-rich permafrost soils and particulate matter into the marine environment (Kipp et al., 2018); and (iii) submarine groundwater discharge (SGD) through subterranean estuaries (STEs), where supra-permafrost groundwater mixes with recirculating seawater before entering nearshore waters (Demir et al., 2024; Kipp et al., 2025). The δ¹⁸O and δ²H results support the theoretical model in which beach groundwater represents a mixture of supra-permafrost groundwater flowing seaward and recirculating seawater advected landward (Fig. 4). The δ¹⁸O and δ²H values measured in water samples collected in 2019, mostly at Tuktoyaktuk Island, Peninsula Point sites, and Crumbling Point (Fig. 1), ranged from −28 ‰ to −10 ‰ and from −215 ‰ to −82 ‰, respectively, the massive ice sample (N=1) presents the most depleted signature (Fig. 4). These depleted values are largely explained by low air temperatures and are typical of permafrost hydrology reported in the western Arctic (Fritz et al., 2011; Utting et al., 2012). The samples are well aligned along the local meteoric water line (LMWL; δ²H = 7.39 × δ¹⁸O – 6.70; Fritz et al., 2022), regardless of their salinity values, except for the three (3) meltwater samples that are slightly below it, likely due to evaporation processes at the surface. The alignment of the different samples along the LMWL regression line suggests a common meteoric origin, probably from the permafrost watershed. The supra-permafrost groundwater flow that transits through the beach sediment does not seem to be affected by surficial processes (e.g. evaporation process), as observed in a few meltwater samples, likely limiting photochemical degradation of the concomitantly transported CDOM/FDOM. However, the high range of isotopic values in beach groundwater samples, from −25 ‰ to −12 ‰ and −192 to −112 ‰ for δ¹⁸O and δ²H, respectively, suggested a mixing between depleted thawed waters (as surficial meltwater and supra-permafrost groundwater) and enriched seawater, as expected in a STE system. These results are consistent with the recent study by Kipp et al. (2025), which reported elevated activities of radon isotopes (²²²Rn) in the same beach groundwater samples. As a noble gas, radon escapes rapidly upon exposure to the atmosphere, indicating a groundwater origin and minimal atmospheric interaction. These authors proposed that thaw water derived from massive ground ice is transported into the supra-permafrost groundwater system, where it is intermittently diluted by radon-poor recirculating seawater (e.g. offshore water from the Mackenzie River) percolating through the beach sediments, depending on tide conditions. Using stable isotope signatures, a first-order estimate of the relative contribution of fresh groundwater to beach groundwater samples can be derived using a simplified two-endmember mixing approach. However, this estimation remains uncertain due to the variability in the isotopic composition of the nearshore seawater endmember, which reflects a wide range of salinities ($\mathrm{5}<S<\mathrm{20}$), as well as the potential mixing of multiple freshwater sources (local and regional supra-permafrost groundwater, and surface meltwater runoff, see Fig. 2) within the supra-permafrost system. To constrain the endmembers, we used the mean isotopic composition of the most saline sample to represent seawater (−74.5 ‰ δ²H; −9.1 ‰ δ¹⁸O; S∼20), and the massive ice sample to represent the fresh supra-permafrost water endmember (−224 ‰ δ²H; −32 ‰ δ¹⁸O). Based on these assumptions, the contribution of fresh supra-permafrost water ranges from approximately 20 % to 76 %. These values should be interpreted as indicative ranges rather than precise quantifications, given the simplified assumptions and the dynamic mixing conditions within the system. Although these estimates do not constrain the total SGD flux through the beach system, they indicate that freshwater inputs can represent a substantial fraction of beach groundwater, consistent with previous SGD studies in the Beaufort Sea region (e.g., Connolly et al., 2020; Demir et al., 2024; Bullock et al., 2024). Our results indicate that supra-permafrost groundwater mixes with recirculated seawater in a nearshore zone, highlighting the dynamic nature of this system along the shoreline, where stable isotopes, salinity, oxygen saturation, and pH likely vary in response to seawater infiltration and dilution.

Figure 4 Isotopic composition of massive ice, meltwater, beach groundwater and seawater samples collected in 2019, the global meteoric water line (GMWL; Craig, 1961) and the local meteoric water line for Inuvik (LMWL, Fritz et al., 2022) are also reported. Samples are colour-coded by salinity (S). Note that only one massive ice sample (N=1) was collected in 2019.

3.3 Behaviour of the DOC and DOM pool

The DOC concentrations dropped from 2360 µmol L⁻¹ in meltwater samples to 236 µmol L⁻¹ in the saltiest seawater sample (S=20; Fig. 5A), with concentrations declining sharply along the supra-permafrost flow path, reaching values below 400 µmol L⁻¹ in beach groundwater samples. These differences among sample types were statistically significant (Kruskal-Wallis, p<0.001). Within the beach groundwater–seawater subset, DOC concentrations showed a nearly linear decrease with increasing salinity and clustered closely around the theoretical mixing line connecting the freshest and saltiest endmembers, with only limited scatter around the predicted conservative trend. This mixing line, based on the most enriched and depleted stable isotope signatures in the subset comprising beach groundwater and seawater samples, corresponds to salinity values of 20.08 (N=1 seawater sample) and 0.18–0.30 (N=8 beach groundwater samples), respectively. These observations indicate that conservative mixing and dilution become the dominant controls on bulk DOC concentrations once supra-permafrost groundwater enters the beach STE and mixes with recirculating seawater. This interpretation is consistent with the relatively small volumetric contribution of fresh supra-permafrost groundwater compared with recirculating seawater along Arctic coastlines (Sect. 3.2). Dilution of the high-concentration, low-volume freshwater endmember therefore contributes substantially to the observed decline in bulk DOC concentrations. Absorption coefficients at 350 nm (a350) ranged from 2.0 to a maximum of 134.0 m⁻¹, the latter observed in a meltwater sample (Fig. 5B). As with DOC concentrations, the a350 values declined sharply along the flow path, with a median value of 24.0 m⁻¹ in meltwater samples and median values of 8.4 and 7.8 m⁻¹ in beach groundwater and seawater, respectively. However, unlike DOC, a350 values in beach groundwaters and seawater show substantial scatter around the theoretical mixing line rather than following a simple conservative dilution trend. At comparable salinities, samples exhibited a broad range of a350 values, with deviations occurring on both sides of the predicted mixing relationship. The variability of the endmembers, particularly in beach groundwaters, limits the robustness of a single mixing model and likely contributes to the observed dispersion, pointing to processes in addition to dilution. This interpretation is further supported by the decoupling between DOC and a350 (Fig. 6), as discussed below. The HIX indices showed a wide range of values across all sample types, with median values decreasing along the continuum from 3.5 in meltwater samples to 1.4 in beach groundwater and 0.8 in seawater (Fig. 5C). Linear mixing models do not apply to ratio-based indices, and non-linear mixing curves could not be reconstructed, as the underlying fluorescence components used to calculate HIX were not retained during data processing. The decrease in HIX along the salinity gradient is consistent with a shift toward less humified DOM. However, as for a350, the variability of HIX values, particularly in seawater samples, limits the strength of this interpretation. In contrast to other parameters, SUVA₂₅₄ increased slightly along the flow path, from 2.4 L mg C⁻¹ m⁻¹ in the meltwater to 2.9 in beach groundwater and 3.4 in seawater (Fig. 5D) and generally followed a dilution pattern along the salinity gradient. However, the variability observed in seawater samples limits the strength of this interpretation.

Figure 5 Distribution of (A) DOC, (B) a350, (C) HIX, and (D) SUVA₂₅₄ across the salinity gradient (S) and among the different sample types (i.e., beach groundwater, massive ice, meltwater, and seawater). Note that only one massive ice sample is included. The zoomed-in panels display only beach groundwater (yellow dots) and seawater (brown triangle) samples, along with a theoretical mixing line (dashed black line) connecting the “fresh” and “salty” endmembers for DOC and a350. These endmembers were defined based on the most enriched and depleted stable isotope signatures in the dataset, corresponding to salinity values of 20.08 (N=1) and 0.18–0.30 (N=8), respectively. As the underlying fluorescence and absorbance components used to calculate HIX and SUVA₂₅₄ indices were not retained, non-linear mixing curves could not be reconstructed. In the boxplots, black lines indicate the medians, whiskers represent the full range of data, and individual points denote outliers.

Figure 6 Relationship between absorption coefficients at 350 nm (a350, m⁻¹) and DOC concentrations (µmol L⁻¹) across all samples along the permafrost-to-nearshore aquatic continuum.

DOC concentration and a350 are commonly used as proxies to assess both the quantity and optical properties of DOM. The relationship between these two parameters provides insights into the biogeochemical sources and transformation processes that shape the composition of the DOM pool along environmental gradients (Stedmon et al., 2003; Fichot and Benner, 2011; Spencer et al., 2013; Massicotte et al., 2017). A strong correlation between DOC and CDOM typically reflects a stable DOM pool, a linear relation frequently observed in freshwater systems (Frenette et al., 2012; Massicotte et al., 2017 and references therein). However, deviations from this relationship, referred to as DOC-CDOM decoupling, suggest shifts in DOM composition and point to selective processing of specific DOM fractions. Microbial degradation, photo-oxidation and organo-mineral interaction mechanisms selectively degrade specific organic compounds, and they lead to shifts in CDOM composition or a loss of CDOM along the land-ocean continuum (Nelson et al., 1998; Del Vecchio and Blough, 2004; Weyhenmeyer et al., 2012; Nelson and Siegel, 2013). In subarctic siliciclastic STE, Couturier et al. (2016), for example, reported a strong decoupling between CDOM-DOC, where HMW CDOM compounds were largely retained within the STE sediments, limiting their export to nearshore adjacent seawater. In such environments, the Fe cycling at the oxic-anoxic interface acts as a geochemical barrier, promoting the sequestration of terrestrial organic compounds (Sirois et al., 2018). As in soils and cohesive marine sediments, molecular fractionation of DOM within STEs tends to promote the in situ degradation and/or preferential export of non-Fe-stabilized, protein-like derived DOM to coastal waters (Linkhorst et al., 2017; Sirois et al., 2018; Waska et al., 2024; Zhou et al., 2023).

The strong decoupling between DOC and a350, particularly in seawater, massive ice, and beach groundwater samples (Fig. 6), indicates that these compounds follow different trajectories along the supra-permafrost flow path, leading to a decrease in humification and a shift toward lower molecular weight material as revealed by HIX indices. Given the light-limited conditions in beach groundwater environments, microbial degradation and mineral–organic interactions likely control the fate and composition of the DOM pool. Despite this general decline in humification and molecular weight of the DOM pool, the impact on DOM aromaticity (as indicated by SUVA₂₅₄) remains relatively weak. SUVA₂₅₄ values tend to increase slightly along the flow path, likely reflecting the dilution effect with the offshore seawater (Fig. 5D). Together, these patterns indicate that DOM composition evolves along the supra-permafrost flow path, resulting in a more decomposed and refractory DOM pool in nearshore waters. This interpretation should be considered in light of the high variability observed in seawater samples, which span a broad range of DOM compositions and reflect the complexity of the system.

3.4 Drivers of DOC – DOM decoupling along the supra-permafrost flow path

In addition to the transformations inferred from optical properties, other sources and processes may contribute to DOM composition along the flow path. The in situ production of a new DOM pool along the supra-permafrost flow path cannot be ruled out despite the low TOC content in the surficial sandy sediment (TOC < 0.50 % dry sediment; Table 2). Similar patterns have been observed in transgressive subarctic STEs (Hébert et al., 2022) and in areas with buried peat lenses (Waska et al., 2021), where ancient sedimentary organic matter preserved in paleosols acts as a hotspot for terrestrial DOM production along the groundwater flow path. Such local production can significantly influence the chemical and optical properties of the beach groundwater DOM pool. Moreover, the presence of permafrost extending below the STE suggests that additional particulate organic matter reservoirs may be mobilized under thawing conditions along the transport. These deeper in situ inputs may explain the very depleted stable water isotopic signature of a few beach groundwater samples (Fig. 4) and the high radon activities (Kipp et al., 2025). Further efforts are needed to improve our understanding of Arctic beach hydrogeology by using geochemical and molecular proxies.

Despite their low TOC content, sandy sediments can act as hotspots of redox activity (Anschutz et al., 2009). In beach groundwater, the absence of light combined with fluctuations in pH, salinity, and redox conditions creates a dynamic environment conducive to a complex network of redox reactions, including the remineralization of organic matter and the cycling of iron between its oxidized (Fe(III)) and reduced (Fe(II)) states. Field observations support the occurrence of these reactions. The presence of reactive Fe-hydroxide coatings on sandy sediment surfaces (Table 2), along with intermittently high total dissolved iron concentrations (median Fe_tot of 1.6 µmol L⁻¹, reaching 671.0 µmol L⁻¹; Fig. 7A) and high Fe:DOC ratio observed in beach groundwater (median of 0.00844, ranging from 0.00353 to 0.92404 in more acidic samples; Fig. 7B), together suggest the potential for the formation of Fe-OM associations through oxidative precipitation of reactive Fe-hydroxide. However, at the Fe:DOC molar ratios observed here, Fe-DOM interactions are expected to drive selective molecular fractionation of the DOM pool rather than substantial bulk OC removal, preferentially scavenging aromatic, high molecular weight humic compounds from solution while leaving bulk DOC concentrations largely unaffected (Linkhorst et al., 2017; Riedel et al., 2013). This is consistent with the sharper decline of a350 relative to DOC along the salinity gradient, while bulk DOC follows more closely conservative mixing behaviour. The concurrent decrease in HIX further supports a progressive loss of humic, high molecular weight compounds along the flow path, independent of simple dilution. The chemical nature of the DOM pool plays a critical role in these interactions: higher molecular weight, aromatic, oxygen-rich terrigenous DOM tends to bind more strongly to Fe(III) than aliphatic, marine-derived DOM (Linkhorst et al., 2017), and slightly acidic conditions further promote Fe(III)-DOM coagulation (Nierop et al., 2002; Amoako et al., 2025). Further evidence consistent with active microbial remineralization within the STE is provided by DIC concentrations measured in beach groundwater in 2019 at the same sites, which exceeded 3000 µmol L⁻¹ (Lizotte et al., 2023), together with a low O₂ saturation. The relative influence of the different mechanisms controlling the molecular fractionation of the DOM pool in supra-permafrost groundwater remains to be determined, and further studies are required to explore the role of mineral–organic interactions and Fe-OM associations in regulating the export of thaw–derived material.

Figure 7 Distribution of total dissolved Fe (Fe_tot, µmol L⁻¹) (A), and Fe-DOC ratio (B) across the salinity gradient (S), and by sample type. The y-axis in each panel is shown on a logarithmic scale to accommodate the broad range of Fe concentrations. Note that there is no Fe_tot data from the massive ice samples and Fe:DOC ratios were calculated only for Fe_tot concentration higher than the detection limit.

While the present dataset does not allow a quantitative separation between mixing-driven dilution and in-situ transformation processes, the convergence of multiple independent lines of evidence discussed above, suboxic conditions, elevated DIC, Fe-OM associations, and the deviation of a350 and HIX from the theoretical mixing line, collectively points to significant in-situ biogeochemical transformation of the DOM pool within the beach STE, operating alongside the mixing dynamics that govern bulk DOC behaviour.

3.5 PARAFAC components in the DOM pool

Among the five fluorescent components identified by the PARAFAC model, the protein-like (tyrosine-like) component C2 dominated the FDOM pool across all sample types and salinities (Table 1; Fig. 8B). While no trend was observed along the salinity gradient, clear shifts in DOM composition were evident along the continuum. Component C2 was strongly negatively correlated with both humic-like components C1 and C4 ($r=-\mathrm{0.87}$ and $r=-\mathrm{0.77}$, respectively; p<0.001), as well as with HIX ($\mathit{\rho }=-\mathrm{0.90}$; p<0.001). In contrast, humic-like components were highly positively correlated with each other (C1 vs C4: r=0.90, p<0.001) and with HIX (ρ=0.96 and ρ=0.93 for C1 and C4, respectively; p<0.001).

Figure 8 Distribution of PARAFAC components (A) C1, (B) C2, (C) C3, (D) C4, and (E) C5 along the salinity gradient (S) and for the different types of collected samples (i.e., beach groundwater, massive ice, meltwater and nearshore seawater samples). Note that only one massive ice sample is reported here. For the boxplots, the black lines represent the median values, the whiskers represent the extent of the data, and the dot points represent the outlier values. Note that there is only one massive ice sample.

Consistent with optical indices (HIX, a350) and DOC patterns previously described, humic-like components (C1 and C4) decreased along the continuum, with the lowest median values observed in seawater (Fig. 8A, D), while the absolute intensity of protein-like DOM (C2) increased significantly across sample types (Kruskal-Wallis, p=0.017), with values significantly higher in seawater than in meltwater (Wilcoxon, p=0.013). Massive ice and meltwater samples were predominantly composed of humic-like, HMW, and terrestrially derived DOM (C1 and C4), consistent with active layer-derived FDOM described by Fouché et al. (2020). The tryptophan-like component C3 and microbial humic-like component C5 showed comparatively low and stable intensities across sample types and salinity gradients, consistent with their minor contributions to the overall DOM fluorescence signal.

As DOM is transported along the supra-permafrost flow path, the terrestrial humic-like signal (C1 and C4) decreases, coinciding with an increase in protein-like, LMW DOM that is consistent with microbial processing. This shift is coherent with the trends observed in optical indices (decreasing HIX and variable a350), and supports the decoupling between DOC and CDOM, indicating that bulk DOC concentrations and DOM composition are influenced by distinct processes. Together, these observations suggest that meltwater and massive ice act as precursors of supra-permafrost solutes exported toward the coastal ocean through Arctic sandy beaches. Along this flow path, the decline in humic-like DOM pool is accompanied by an increase in biologically derived, tyrosine-like FDOM (component C2), typically associated with microbial production (Table 1 and references therein). These patterns are consistent with transformations reported by Fouché et al. (2020), where microbial processing converts DOM into LMW protein-like material. In this context, the presence of non-Fe-stabilized DOM and redox oscillations likely promotes microbial reprocessing and transformation. These processes are consistent with the production of LMW and protein-like compounds, and they contribute to the observed compositional changes along the flow path .

4 Conclusion

This study reveals a rapid decline in DOC and CDOM concentrations along the diffuse and non-point source supra-permafrost groundwater flow from coastal permafrost bluffs through sandy beach sediments and to adjacent nearshore seawater. This decline is likely largely influenced by the dilution of a high-concentration, low-volume supra-permafrost groundwater endmember with volumetrically dominant recirculating seawater. At the same time, our results provide evidence consistent with organo-mineral interactions and microbial reprocessing contributing to a concurrent compositional transformation of the DOM pool within the STE. Rather than acting as a quantitative DOC sink, the Arctic beach STE can be described as a DOM quality filter, where these processes are associated with a relative decrease in humic-like, high molecular weight compounds and a shift toward more protein-like, microbially-derived DOM before export to coastal waters. Our results suggest that nearshore and intertidal zones, particularly subterranean estuaries where supra-permafrost groundwater mixes with recirculated seawater, contribute to shaping DOM composition before its export to coastal waters. Similar transformation processes have been reported in temperate and tropical STEs, suggesting that Arctic beach systems share common biogeochemical features despite the influence of permafrost. However, whether Arctic beaches primarily act as transient filters, long-term carbon sinks, or dynamic biogeochemical reactors remains unresolved, underscoring the need for further research into their contribution to the transfer of carbon from permafrost to the ocean. Given ongoing permafrost thaw and increasing coastal erosion, improving our understanding of these nearshore processes is important for constraining Arctic carbon cycling and predicting their potential role in climate feedback.