Seasonal contrast in rare earth elements concentration in sediment of the Mackenzie Delta

This study reports on the concentration of rare earth elements (REE) along with ancillary geochemical parameters at 11 locations across the Mackenzie River, its delta and coastal waters, both under ice and in open water. Specifically, we analyzed REE, carbon, and redox-sensitive elements (Fe, Mn) in 98 sediment samples and 96 porewater and overlying water samples collected under ice before the spring freshet (April–May) and in open water in early fall (August–September). While sediment REE concentrations remained relatively stable across seasons, results revealed a striking contrast between the two sampling seasons in the porewater, where REE concentrations were nearly two orders of magnitude higher under ice (avg. 216 nmol L⁻¹) than under open water in the fall (avg. 3.20 nmol L⁻¹). Similarly, dissolved organic carbon (DOC) concentrations were approximately one order of magnitude higher under ice than in the fall. Sediment REE concentrations were positively correlated to those of Fe and Mn, particularly under ice, consistent with control by adsorption processes onto their (oxy)hydroxides. In the porewater, winter and fall samples form distinct clusters based on concentration magnitudes. Chromophoric properties of dissolved organic matter (DOM) in the overlying water suggest that under-ice DOM was characterized by low aromaticity, older material compared to the more aromatic, humic-rich DOM measured in open water. We conclude that under-ice conditions, chiefly cold temperature, allow for DOM accumulation in the porewater, which, combined with other possible REE enrichment mechanisms in the porewater, such as REE–carbonate complex formation and exclusion during ice formation, contributes to the elevated winter REE concentrations observed here. To our knowledge, this is the first report of such large seasonal fluctuation in dissolved REE in the fluvial-marine transition zone of the Mackenzie, the largest riverine influence on the Arctic Ocean.

1 Introduction

Arctic coastal environments and deltas represent important biogeochemical hotspots at the interface between terrestrial and marine systems. They serve as conduits and processing sites for trace elements and contaminants. The Mackenzie River system, which drains a large portion of northwestern Canada and discharges into the Beaufort Sea, has received particular attention for its role in transporting contaminants, with recent studies revealing pronounced seasonal and spatial patterns in mercury exposure throughout the Mackenzie watershed, linked to hydrological conditions and biogeochemical processes (Jermilova et al., 2025). These findings highlight the importance of understanding trace element cycling in Arctic river systems, where seasonal variations in temperature, organic matter dynamics, and ice cover can influence contaminant behaviour.

Seasonal ice cover responds strongly to global environmental changes (Gudasz et al., 2010; Landy et al., 2022). Seasonality significantly influences elemental distribution and geochemical processes by affecting water and sediment temperatures, site salinity through ice exclusion or freshwater dilution (Degrandpre et al., 2021; Finlay et al., 2006), photosynthesis activity (Retamal et al., 2008; Semkin et al., 2022), river discharge patterns (Gareis and Lesack, 2017), and carbon input from surrounding terrestrial environments (Holmes et al., 2012; Gareis and Lesack, 2017; Liu et al., 2022). These seasonal dynamics collectively regulate environmental processes and the associated mobility and distribution of reactive elements.

Among the elements of growing scientific interest, rare earth elements (REE) offer insights both into biogeochemical processes due to their coherent geochemical behaviour and sensitivity to environmental conditions (Škerlep et al., 2025; Ye et al., 2019), and as tracers of sediment sources due to their distribution representing geological signatures (Bossé-Demers et al., 2025). Interest in REE thus continues to grow along with their applications in advanced technologies and their increasing utility as proxies for environmental processes (Alonso et al., 2012; Haque et al., 2014; Tostevin et al., 2016; Guo et al., 2024; Grenier et al., 2022).

Generally, REEs tend to behave coherently in the environment due to their similar reactivity patterns and predominant oxidation state, +III (Wall, 2021; Rollinson and Pease, 2021). Organic carbon quantity and quality (i.e., aromaticity, molecular weight and functional group composition) may influence the distribution of REE through organic matter complexation, as dissolved organic compounds can form stable complexes with REE, affecting their mobility and fractionation patterns (Davranche et al., 2004; Pourret et al., 2007a, c; Marsac et al., 2021). Marginson et al. (2024) found positive dissolved organic carbon (DOC)–REE and humic-like dissolved organic matter (DOM)–REE correlations in the sub-Arctic George River, highlighting the influence of organic matter quantity and quality on REE distribution in these systems.

Rare earth elements in aquatic systems partition between dissolved, colloidal, and solid phases through competitive binding to Fe and Mn (oxy)hydroxides and natural organic matter. Iron oxyhydroxides exhibit strong REE affinity, but when REEs are complexed with humic substances, the organic ligand rather than the metal ion interacts with the oxide surface, effectively suppressing direct REE mineral sorption (Davranche et al., 2004). In Fe-rich systems, REE association with Fe–organic matter colloids at circumneutral pH is primarily controlled by the organic component rather than the iron phase itself (Tadayon et al., 2024b), with implications for REE mobility during seasonal redox oscillations. Organic matter binds REEs through heterogeneous sites: at high loading, carboxylic groups, while at low loading, phenolic and chelate groups control binding (Sonke, 2006; Davranche et al., 2015).

Dissolved REE concentrations in Arctic lakes and rivers vary considerably across sites with otherwise similar geochemical characteristics (MacMillan et al., 2017; Pokrovsky et al., 2016), yet the processes driving this variability remain poorly understood. Insights may be gained from studying sediment porewater – the geochemical compartment that mechanistically connects the solid phase in the sediment to the overlying water column – which to our knowledge have never been studied for REE in Arctic river systems. Porewater REE cycling contributes significantly to REE mobility in marine and coastal sediments (Abbott et al., 2015; Deng et al., 2022), yet systematic sampling of this compartment has been limited, in part by the analytical challenges associated with low concentrations and small sample volumes (Abbott et al., 2015), as well as by the logistical constraints inherent to remote Arctic fieldwork.Here, we set out to better understand the magnitude and controlling factors of seasonal variations in REE concentrations. Previous work at other Arctic sites has hinted at various causes for such variation at other sites, including reduced bacterial degradation of organic matter at low temperatures preserving DOC–REE complexes (Arnosti et al., 1998), and enhanced REE–carbonate complex stability at low temperatures (Marginson et al., 2024). We investigate specifically the linkages between REE concentrations, DOC concentrations and DOM quality. To do so, we analyzed REE, DOC $/$ DOM, and Fe and Mn in sediment, overlying water and porewater samples across 11 sites during two sampling seasons in the Mackenzie River Delta. We use this dataset to improve our understanding of seasonality effects on REE in a cold region and to better understand how these impact REE distribution.

2 Materials and Methods

2.1 Study Region

The Mackenzie River (Northwest Territories, Canada) flows from south to north into the Beaufort Sea through the second-largest delta in the Arctic. The Mackenzie drainage basin encompasses 1.8×10⁶ km² within Canada, including the Rocky Mountains and Mackenzie Mountains. About half of this area is situated in continuous or discontinuous permafrost zones (Abdul Aziz and Burn, 2006; Hill et al., 2001; Holmes et al., 2012). The river is covered by ice from late October to late May (Hill et al., 2001).

The delta exhibits pronounced seasonal and interannual discharge variability that significantly influences the Beaufort Sea (Hill et al., 2001; Mulligan and Perrie, 2019; Nghiem et al., 2014; Emmerton et al., 2008a). During winter, the river discharge diminishes under extensive ice coverage (Goñi et al., 2000; Hill et al., 2001). The ice break-up period in late May and early June triggers the spring freshet, during which snow and ice meltwater cause increased water discharge (e.g., Hill et al., 2001). By early June, freshwater extends more than 100 km from the river mouth onto the shelf (Juhls et al., 2022). Approximately 70 % of annual freshwater discharge occurs between May and September (Leitch et al., 2007; Macdonald et al., 2012), with recent years showing an overall increase in discharge volume (Kopec et al., 2024; Rood et al., 2017). With the contrasting seasonal patterns, the Mackenzie River sets an ideal study site for evaluating seasonality effects on REE distribution in coastal settings. Our two sampling campaigns corresponded to contrasting hydrological conditions (Lizotte et al., 2023). The spring campaign (Leg 1) took place during baseflow (∼4000–5000 m³ s⁻¹), prior to the spring freshet, while the fall campaign (Leg 4) was conducted during the recession limb of the annual hydrograph (∼10 000 m³ s⁻¹).

The Mackenzie River is the principal sediment source to the Arctic Ocean. It delivers an estimated 1.28×10¹¹ kg of sediment annually (Holmes et al., 2002; Stein et al., 2004) and contributes 90 %–95 % of total sediments to the Beaufort Shelf (Vonk et al., 2015). Its sediment sources include the Peel (21×10⁹ kg), Arctic Red (7×10⁹ kg), and Liard (35–45×10⁹ kg) rivers (Vonk et al., 2015). The sediments are mostly composed of fine-grained silt and clay (Hill et al., 2001; Holmes et al., 2012). These sediments are transported through relatively shallow river channels reaching a maximum depth of approximately 10 m (Hill et al., 2001; Mulligan and Perrie, 2019).

2.2 Sample Retrieval and Preservation

The Mackenzie River and its downstream coastal areas were sampled within the context of the EU-H2020 project Nunataryuk. Research teams collected various physicochemical parameters at these sites, with their data and conclusions available in separate publications (Juhls et al., 2022; Bertin et al., 2023; Lizotte et al., 2023). Sampling was performed across two rivers to coast transects along the western and eastern outflow area of the Mackenzie River during two distinct periods: before the spring freshet in April and May 2019, with fully ice-covered sites (labelled W in Table 1), and in early fall in August and September 2019, representative of open-water conditions at the end of the season of high biological productivity (labelled F in Table 1). This corresponds to Leg 1 (17 April to 3 May 2019) and Leg 4 (26 August to 9 September 2019) of the Nunataryuk field campaigns (Lizotte et al., 2023). These sampling campaigns took place before and after large freshwater and sediment discharge events occurring during the summer. Hereafter, we name those sampling periods “winter” and “fall” as they represent those respective conditions. The western sampling transects extended from Inuvik (68.35° N, 133.68° W) towards the western river mouth and the coastal waters of the Mackenzie Bay (69.14° N, 136.85° W). The eastern transects ran from the east river mouth to the coastal water of Kugmallit Bay (69.55° N, 133.41° W), as presented in Fig. 1.

Table 1 Sampling station channel, water column depth, geographical coordinates, bottom water temperature, and salinity.

* (W) Sampled in late winter. (F) Sampled in early fall.

Figure 1 Sampling sites across the Mackenzie Delta. Squares (blue) represent winter sampling sites and circles (brown) represent fall sampling sites. Towns and hamlets are shown by triangles (purple). Blue arrows indicate general flow direction toward the Beaufort Sea along the main channels of the delta. (Source: ESRI | Powered by Esri)

Sediment was sampled with a gravity corer (UWITEC, Austria) with a 9 cm diameter core liner. Winter, under ice sampling was conducted through a hole in the ice, with the gravity corer maintained using a tripod. Sites were accessed via helicopter or snowmobile. Fall sampling was conducted on a boat. At each site, two sediment cores were taken, with one used for sediment subsampling at 1 cm intervals and the other for porewater extraction (see below). Subsampling was conducted at the Aurora Research Institute (ARI) in Inuvik for western sites and at the Tuktoyaktuk Learning Centre for eastern sites. The subsamples were placed in Falcon cups (Corning, USA) and frozen. Back in the laboratory at Université Laval, frozen sediment samples were freeze-dried and homogenized using an agate mortar and pestle. Aliquots of these homogenized samples were used for subsequent solid phase analyses.

Porewater and overlying water were sampled from pre-drilled core liners with holes covered by tape at 1 cm intervals. Upon retrieving the core, acid-washed 5 cm Rhizon samplers (Rhizosphere Research, Netherlands) with 0.15 µm PES membranes were inserted into the core. Porewater was retrieved using acid-washed syringes (VWR, Canada) by creating a vacuum. This allowed collection of 7–10 mL of porewater (Seeberg-Elverfeldt et al., 2005), starting at 1 cm above the sediment-water interface (the overlying water) and continuing down to the depth where bottom clays prevented Rhizons insertion (5–20 cm). Porewater replicates could not be recovered due to the limited volume of porewater available for analysis in the samples. The collected porewater was distributed into different vials for preservation. For DOC concentration measurements, 2 mL of porewater was delivered to amber glass vials fitted with Teflon-lined caps. The vials were previously cleaned in hydrochloric acid (HCl) (10 %) baths followed by sodium hydroxide (NaOH) (0.1 M) baths and ultrapure water baths for 24 h each, then combusted at 450 °C overnight. Teflon-lined caps were washed separately in an ultrapure water bath. For major and trace elements, samples were transferred to 15 mL clear polyethylene terephthalate (PET) Vacutainer tubes for the winter samples. Due to the slow vacuum retrieval rate of the Vacutainer, we switched to acid-washed syringes combined with 15 mL high-density polyethylene (HDPE) centrifuge tubes for the fall samples. Both containers were acid washed in 4 % nitric acid (HNO₃) baths and rinsed with ultrapure water baths prior to use. Approximately 4 mL of porewater sample was amended with 320 µL of 50 % double-distilled Omnitrace HNO₃ (Fisher Scientific, Canada) to obtain a final concentration of ≈4 % HNO₃ suitable for preservation. Field blanks could not be collected because the ultrapure water available at the Tuktoyaktuk field laboratory had degraded during storage and was found to be unsuitable. This issue was identified post-campaign during data analysis, after both sampling seasons had been completed, precluding any corrective action. We refer to field blanks from porewater campaigns conducted the same year in the same laboratory using identical handling and preparation procedures for indirect evidence against systematic contamination (see Table A1).

Prior to major and trace metal analysis in the sediment, 50 mg of homogenized sediment was placed in a Teflon microwave reaction vessel (EasyPrep, CEM Corporation, Canada) amended with 6.5 mL of Aristar Plus 12 M HCl (VWR Canada) and 3.5 mL of doubled-distilled 16 M Omnitrace HNO₃ (Fisher Scientific, Canada). A procedural blank and a certified reference material sample (CRM MESS-4, NRC Canada) were included in the microwave carousel. The samples were digested in a MARS 5 microwave (CEM Corporation, Canada) with a 20 min ramp to 1600 W (setting of 800 maximum PSI, 240 °C), followed by a 20 min hold time and a 30 min cooling period. After digestion, the samples remained in the vessels overnight to cool completely. The vessels were opened and placed on a temperature-controlled Digiprep block (SCP Sciences, Canada) for evaporation at 120 °C for approximately 3.5 h, until nearly dry. The samples were then transferred under a laminar flow hood (ESCO, USA) for recovery in 2.87 mL of doubled-distilled 16 M Omnitrace HNO₃ (Fisher Scientific, Canada). This solution was then transferred to an acid-washed 50-mL HDPE centrifuge tube (VWR Canada). The Teflon vessels were rinsed three times with ultrapure water, and the solution was added to the centrifuge tube. The final solution was brought to a volume of 50 mL to achieve ≈4 % HNO₃.

2.3 Instrumental Analysis

Major elements in the sediments (Fe and Mn) were analyzed by ICP-OES (Thermo Scientific iCAP 7400, USA) using iridium (Ir) as an internal standard. Analysis of the digested CRM MESS-4 yielded accuracy greater than 99 % for Fe (n=11) and 96 % for Mn (n=10). REE in the sediment were measured by triple-quadrupole ICP-MS (Agilent 8900, Agilent Canada) using iridium (Ir) as an internal standard. The method for REE in sediment was validated using MESS-4 reference material, yielding accuracy greater than 91 % for La, Ce and Eu (n=21), greater than 75 % for Nd (n=12), and 97 % for Sm (n=12). Lu could not be assessed, likely due to low recovery rates. The good recovery of Yb (72 %), another heavy REE, suggest that the low recovery of Lu is not due to an instrumental bias. Dissolved REE as well as dissolved Fe and Mn were also analyzed by ICP-QQQ-MS, using Te and Rh as internal standard. The method for dissolved REE was validated using the CRMs SLRS-6 (NRC Canada; a river water reference material with certified dissolved REE concentrations; Yeghicheyan et al., 2019) and TM-DWS.3 (Environment and Climate Change, Canada), with accuracy exceeding 92 % (n=11 for REE and n=3 for Fe and Mn) for all analytes but Eu, with an accuracy of 85 % (n=11). The dominant isotope was monitored for each REE. The analyses were performed in O₂-mode with a collision cell to minimize isobaric interferences for REE and Fe, while Mn was performed in single-quad mode. Matrix separation prior to ICP-MS analysis was not required given the predominantly freshwater nature of the samples. (Yeghicheyan et al., 2019). A supplementary table reporting measured concentrations, certified values, and percent recoveries for MESS-4 and SLRS-6 is provided in the Appendix (Tables A2 and A3). Of the 1344 total measurements, 335 (25 %) were below the limit of detection (LOD) and substituted with 0.5× LOD to enable statistical analysis. Antweiler and Taylor (2008) demonstrated that half-LOD substitution is an unbiased and adequate technique for datasets with less than 70 % censored values; our 25 % censoring rate falls well within this acceptable range. We acknowledge that this approach is inferior to non-parametric methods such as Kaplan–Meier and note this as a limitation.

Solid-phase carbon content in the sediment was determined using a CHN analyzer (Flash 2000, Thermo Scientific, Canada). Homogenized sediment samples (3 mg) were weighed and placed in tin capsules for analysis. The method was validated using reference materials (cystine and sulfanilamide), yielding accuracy greater than >99.3 % (n=10). Porewater DOC analysis was performed using a total organic carbon and nitrogen (TOC/TN) analyzer (Vario Cube, Elementar, Germany). In-house carbon standards were used to validate the method, achieving accuracy of >90 % (n=17).

Fluorescent DOM (FDOM) from Matsuoka et al. (2021b) were measured using a spectrofluorometer (Aqualog, Horiba, Japan) with corrections for inner-filter effects and Raman-Rayleigh scattering. Fluorescent components were identified via a parallel factor analysis (PARAFAC) and compared with the OpenFluor database. Chromophoric DOM (CDOM) absorption spectra (200–722 nm) were measured in triplicate within 12 h of collection using an UltraPath liquid waveguide system (Juhls et al., 2021; Matsuoka et al., 2021a; Matsuoka et al., 2021b) to derive specific UV absorbance at 254 nm (SUVA₂₅₄). From the absorption and fluorescence measurements described above, SUVA₂₅₄ was calculated as the absorption coefficient at 254 nm normalized to DOC concentration (L mg C⁻¹ m⁻¹). We used the measured humification (HIX) and biological (BIX) indices, as well as PARAFAC components. From the latter, we extracted the following excitation (ex) $/$ emission (em) wavelength pairs (Hansen et al., 2016): C (340_ex:440_em), A (260_ex:450_em), T (275_ex:304_em), and M peaks (300_ex:390_em).

2.4 Data Analysis

Spatial interpolation of dissolved and sedimentary REE, DOC, Mn and Fe concentrations was performed using Inverse Distance Weighted (IDW) interpolation in QGIS (v. 3.40.9). To facilitate visual comparison across compartments with different concentration ranges, data were normalized within each compartment (dissolved concentration, sediment) by dividing all values by the compartment-specific maximum concentration, resulting in relative concentrations ranging from 0 to 1.

Spearman rank correlations were performed to assess relationships between geochemical variables without assuming linear distributions using the ranks of Pearson's correlation calculated in Grapher v.26 (Golden Software, USA). Correlations were considered statistically significant at p<0.05. The remaining data treatment was performed in R (R Core Team, 2025). Principal component analysis (PCA) was performed on seven variables: dissolved ΣREE, DOC, dissolved Fe, dissolved Mn, SUVA₂₅₄, temperature and salinity. Six variables (REE, DOC, Fe, Mn, SUVA₂₅₄ and salinity) were tested for normality using Shapiro-Wilk tests. Temperature was not tested for normality as it directly represents the seasonal variations. Variables exhibiting right skewness (REE: W=0.479, p<0.0001; DOC: W=0.865, p<0.0001; Fe: W=0.532, p<0.0001; Mn: W=0.552, p<0.0001; Salinity: W=0.305, p<0.0001) were log10-transformed (Reimann and Filzmoser, 2000). SUVA₂₅₄ (W=0.781, p<0.0001) and temperature were not transformed. Finally, MANOVA and ANOVA were performed on six variables, excluding temperature as it defines the seasonal grouping, using Pillai's trace as a criterion for seasonal separation. This resulted in biplots with loading vectors scaled ×4 with 95 % confidence ellipses. PCA, MANOVA and ANOVA were calculated using the base stats package from R. 95 % confidence ellipses were calculated using ellipse v0.5.0 package. Data normalization prior to PCA followed z-score standardization, centring each variable to a mean of zero and scaling to a standard deviation of 1 (Reimann and Filzmoser, 2000).

3 Results

Water temperatures in winter were ∼0 °C at all stations, while they ranged from 7.3 to 10.6 °C ($\stackrel{\mathrm{‾}}{x}=\mathrm{9.0}\pm \mathrm{1.6}$ °C) during fall (Lizotte et al., 2023). Salinity was similar across sites, with values ranging from 0.002 ‰ to 0.199 ‰. One anomalous value of 14.26 ‰ was measured at site 810(F), located in Kugmallit Bay (Table 1). This elevated salinity persists despite higher fall river discharge (∼10 000 m³ s⁻¹) as a remnant of summer marine intrusion: in late summer, reduced river flow allows Beaufort Sea water to intrude into the delta, with surface salinities of up to 25 ‰ recorded in Kugmallit Bay (Juhls et al., 2022). The same site was also slightly brackish in winter (810(W)), consistent with incomplete flushing before freeze-up. Figure 2 illustrates the marked differences in concentrations of dissolved ∑REE between both sampling seasons in the water column (hereafter referred to as overlying water, corresponding to the river water column above the sediment-water interface) and porewater. Similar plots for carbon, Fe and Mn data are presented in supporting information (Figs. A1, A2, and A3 respectively). Although our spatial coverage is admittedly too sparse to allow quantitative interpolation across the Delta, the figure usefully highlights the seasonal contrast in ∑REE concentrations, and the absence of such contrast in the sediment phase, both for ∑REE and carbon. ∑REE concentrations in sediments exhibited minimal seasonal variation, with winter values ranging from 588 to 1088 µmol kg⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{816}$ µmol kg⁻¹) and fall values ranging from 591 to 1250 µmol kg⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{952}$ µmol kg⁻¹). Similarly, winter carbon concentrations in the sediment ranged between 2.89 % and 4.87 % ($\stackrel{\mathrm{‾}}{x}=\mathrm{3.66}$ %) and fall carbon concentrations ranged between 2.40 % and 4.91 % ($\stackrel{\mathrm{‾}}{x}=\mathrm{3.41}$ %).

Figure 2 Spatial interpolation of normalized average ΣREE concentrations in the Mackenzie River Delta for the overlying water (water column, WC, top panels), porewater (PW, middle panels) and sediment (Sed, bottom panels) for the winter (left) and fall (right) sampling seasons. Concentrations were normalized for dissolved and solid compartments with the compartment maximum, resulting in relative values between 0 (blue) and 1 (red). Grey circles show sampling sites. This figure is intended to illustrate the dominance of seasonal over spatial variability in REE concentrations; absolute concentrations are reported in the text and in Figs. 3–5. (Source: ESRI | Powered by Esri)

In contrast to the sediment, porewater and overlying water samples exhibited strong seasonal trends for both ∑REE and DOC (Figs. 2, A1). To further illustrate the seasonal contrast, we plotted the distribution of concentrations across defined intervals (Fig. 3), which shows that winter samples are within the highest concentration ranges for both ∑REE and DOC. These analytes exhibited lower concentrations in fall. During winter, porewater ∑REE concentrations ranged from 2.98 to 1624 nmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{217}$ nmol L⁻¹), and overlying water concentrations ranged from 19.8 to 926 nmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{196}$ nmol L⁻¹). In fall, porewater ∑REE concentrations ranged from 0.26 to 67.96 nmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{3.12}$ nmol L⁻¹) and overlying water concentrations ranged from 0.36 to 0.62 nmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{0.50}$ nmol L⁻¹). This represents a difference of almost two orders of magnitude in average ∑REE concentrations from winter to fall. Similarly, average DOC concentrations are lower by about an order of magnitude in fall compared to winter. In winter, porewater DOC concentrations ranged from 13.71 to 58.32 mmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{22.66}$ mmol L⁻¹) and overlying water concentrations ranged from 10.40 to 25.03 mmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{19.41}$ mmol L⁻¹). In comparison, fall concentrations ranged from 0.15 to 20.66 mmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{2.86}$ mmol L⁻¹) for porewater, and from 0.29 to 19.36 mmol L⁻¹ ($\stackrel{\mathrm{‾}}{x}=\mathrm{6.08}$ mmol L⁻¹) for overlying water.

Figure 3 Distribution of samples across logarithmic concentration intervals (n=96) for ∑REE (mmol L⁻¹) and DOC (mmol L⁻¹) for the winter (blue) and fall (brown) sampling seasons.

Spearman correlation analysis revealed significant relationships in sediment samples for both seasons (Fig. 4). In winter sediments, Fe and Mn, Fe and Ce as well as Mn and Ce were positively correlated (Fig. 4a), while Ce and Nd were strongly positively correlated (Fig. 4d). In the dissolved phase (overlying water and porewater), significant correlations were only found among redox-sensitive metals. Fe and Mn were moderately correlated both in winter and fall (Fig. 5a) or Ce (Fig. 5b). Notably, no significant correlations were detected between REE or metals and DOC in either season (Fig. 5c–d). Winter and fall samples form instead distinct clusters based on concentration magnitudes (Fig. 5c), with no overlap between seasons. This seasonal separation reflects the approximately two-order-of-magnitude difference in REE concentrations and one-order-of-magnitude difference in DOC concentrations between seasons (Figs. 2, A1). Results from the MANOVA on winter (n=53) and fall (n=36) samples support this (Pillai's trace = 0.883, F-Stat = 103.14, degree of freedom = 6 and p<0.0001). The combination of large effect size (Pillai's trace = 0.883) and high statistical significance indicates that winter and fall samples have distinctly different multivariate profiles across the measured variables. Finally, PCA analysis (Fig. 6; Tables A4 and A5) of the porewater Mn, Fe, DOC, ∑REE concentrations, SUVA₂₅₄ and temperature data shows significant multivariate separation, capturing 79 % variance in two dimensions. It should be noted that PCA components are descriptive and reduce dimensionality; they do not establish causality. Other unmeasured variables may contribute to the variance structure, and the attribution of mechanisms to PCA axes should be interpreted with this in mind. Finally, SUVA₂₅₄ and HIX indices are higher in the fall than the winter, while values for BIX are slightly lower in fall than in winter. The C : T and A : T peak ratios are higher in fall compared to winter, as is the C : M peak ratio, but to a lesser extent than the two other peaks (Table 2).

Figure 4 Correlation between solid-phase concentrations of Mn and Fe (a), Ce and Fe (b), Ce and Mn (c), and Ce and Nd (d) for the winter (blue) and fall (brown) sampling seasons. Spearman's rank (r) are indicated only for significant correlation (p<0.05).

Figure 5 Correlation between porewater Mn and Fe (a), Ce, Mn and Fe (b), ΣREE and DOC (c), and Fe (circles) and Mn (triangles) over DOC (d) for the winter (blue) and fall (brown) sampling seasons. Spearman's rank (r) are indicated only for significant correlation (p<0.05). Highlighted zones on panel (c) indicate the two distinct concentration ranges observed between the two sampling seasons.

Figure 6 PCA biplot of log-transformed seasonal dissolved Mn, Fe, REE and DOC concentration, as well as SUVA₂₅₄, temperature and salinity data. Points represent individual samples colored by season (winter in blue, n=53 and fall in brown, n=36) with 95 % confidence ellipses. Arrows indicate variable scaled loadings. Black circle represents the equilibrium descriptor contribution (EDC). Data shown are for porewater and overlying water samples.

Table 2 Average chromophoric properties of DOM in the water column in the Mackenzie Delta at sampling sites during winter and fall. HIX, BIX, C : T, A : T and C : M are dimensionless ratios; SUVA₂₅₄ is expressed in L mg C⁻¹ m⁻¹.

4 Discussion

4.1 Sediment REE stability and Fe-Mn control

Our results show that ∑REE concentrations in the sediments are slightly higher in the fall (952 µmol kg⁻¹) than in the winter (816 µmol kg⁻¹). The strong positive correlations between sediment ∑REE concentrations and redox-sensitive elements (Fe and Mn) (Fig. 4) are consistent with Fe and Mn redox-driven recycling controlling REE distribution in the solid phase. Both field (Ye et al., 2019; Toyoda et al., 1990; Takahashi et al., 2015) and laboratory evidence (Bau, 1999) show that REE can be sequestered by Fe and Mn (oxy)hydroxides. The fact that correlation between Fe, Mn and ∑REE in our study are stronger under-ice suggests that near-zero temperatures are particularly favourable for the control of Fe and Mn on REE. Reduced microbial respiration along with the thermal suppression of microbial activity likely create stable conditions near the sediment–water interface (Bouvet et al., 2025; Arnosti et al., 1998) for REE to accumulate onto Fe and Mn (Dang et al., 2022). The covariance between Ce and Nd (Fig. 4d) indicates that REE co-vary together, allowing the use of ΣREE concentrations as the key variable in this study. Ce and Nd were selected as representative REEs on Fig. 4 because they are among the most abundant in natural waters. Their co-presentation allows detection of any Ce anomaly relative to the smooth lanthanide trend defined by Nd; the absence of such an anomaly in our dataset supports treating the REE series as a coherent group.

4.2 Dissolved Organic Matter Quantity versus Quality in REE Complexation

Porewater and overlying water (i.e., water column) dissolved ∑REE concentration data shows a seasonal trend (Figs. 2 and 3), with two distinct populations of porewater ∑REE and DOC measurements (F-statistics of 409.13 and 220.09 respectively, p values <0.0001; Table A7) (Fig. 5c). PCA shows that DOC and ∑REE co-vary along the PC1 axis (53 % variance), inversely correlated with temperature and SUVA₂₅₄, due to the 180° angle between these vectors. Meanwhile, Fe and Mn co-vary along PC2, independently of the temperature, DOC concentrations and DOM quality or of ∑REE concentrations, given the 90° separation between the vectors. Salinity is represented similarly, with an angle that is not representative of a strong co-variation. This seasonal variability, where REE and DOC co-vary and are inversely related to temperature and DOM aromaticity (SUVA₂₅₄), is consistent with OM complexation of REE during colder periods being controlled by DOM quantity. The next largest variance is explained by the second principal component (PC2) encompassing redox-sensitive processes, for which dissolved Fe and Mn concentrations are reasonable proxies in the sediment porewater (Schulz and Zabel, 2006).

Collectively, those findings align with recent work in the subarctic George River, Canada (Marginson et al., 2024) and the temperate Sleepers River, USA, in which dissolved water-column REE and DOC exhibited positive correlations (Norton and Shanley, 2025). Seasonal variability in REE dynamics in rivers has also been evidenced for tropical systems, where dissolved REE concentrations showed strong seasonal variation driven by watershed runoff during wet seasons versus scavenging by organic-rich particles during dry seasons (Dang et al., 2023). While our Arctic system operates under different temperatures and biological productivity regimes, both studies demonstrate that seasonal changes in organic matter quality and quantity influence dissolved REE dynamics in river systems with pronounced seasonal variability. Several studies point to organic matter quality playing a crucial role in REE complexation and distribution in aquatic systems (Catrouillet et al., 2020; Pourret et al., 2007c; Marsac et al., 2021; Tadayon et al., 2024a). Overall, it is generally accepted that complexed REE account for up to 95 % of the dissolved phase in freshwater (Tang and Johannesson, 2003; Johannesson et al., 1995; Pourret et al., 2007b; Revel et al., 2025).

The strong seasonal contrast in DOC concentration (Fig. 2) can be investigated in the light of the chromophoric properties of DOM, namely low SUVA₂₅₄ and HIX values in the winter relative to the fall (Table 2). Those values point to winter OM samples having a lower aromaticity and a lower proportion of humic acid than the fall sample (Hansen et al., 2016). Similarly, BIX values point to a higher proportion of autochthonous dissolved organic matter in the winter and the C : T, A : T and C : M peak-ratio point to a lower humic character in winter. Collectively, those indicators of autochthonous DOM in the winter are consistent with the enzymatic breakdown of solid-phase OM to aqueous DOM maintaining a significant activity at low temperatures (Davidson et al., 2006), while the subsequent mineralization of aqueous DOM remains inhibited. This differential temperature sensitivity may explain in DOM accumulation in the porewater in the winter (German et al., 2012; Schädel et al., 2016). This lower winter OM degradation is consistent with the markedly lower bacterial abundance observed in winter samples compared to fall samples for the Mackenzie Delta (Lizotte et al., 2023). In fall, increased DOM degradation in the porewater suppresses the autochthonous signal, leaving catchment-derived terrestrial DOM. The shift in DOM character between seasons likely extends to molecular weight distributions, with autochthonous winter DOM, characterized by enzymatic breakdown products, typically exhibiting lower molecular weight compared to the terrestrial, humic-rich DOM observed in fall (Hansen et al. 2016). However, the literature present conflicting evidence regarding the role of molecular weight in REE complexation, with some studies showing that REE generally preferentially complex with high molecular weight OM (Tang and Johannesson, 2003; Catrouillet et al., 2020; Tang and Johannesson, 2010), and other preferentially complex with low molecular weight OM (Zilber et al., 2024). Without size fractionation data, we cannot evaluate which mechanism dominates in our system. Future studies should consider pyrolysis-GC-MS to directly characterize DOM functional groups and molecular composition beyond what optical indices can resolve, and size-fractionation approaches (e.g., dialysis or ultrafiltration) to identify which DOM size fraction dominates REE complexation. Explicit characterization of REE fractionation patterns (LREE/HREE ratios) would also strengthen mechanistic interpretation. Thus, we find two distinct porewater DOM pools depending on the season. Fall DOM has higher-affinity binding sites (phenolic groups, log K=4.93) compared to winter DOM (carboxylic groups, log K=3.29) (Marsac et al., 2011). Using carbon-normalized binding site density for two contrasting DOM types as proxies (Mueller et al., 2023), a higher site density (0.023 µmol mg C⁻¹) representative of our fall DOM and a lower site density (0.0105 µmol mg C⁻¹) representative of our winter DOM, and applying these to our measured seasonal DOC concentrations (∼267 mg C L⁻¹ in winter vs. ∼41 mg C L⁻¹ in fall), we estimate the absolute concentration of binding sites: ∼0.94 µmol L⁻¹ in fall vs. ∼2.80 µmol L⁻¹ in winter. Despite fall DOM having higher binding site density per unit carbon, the higher DOC concentration in winter results in ∼3-fold more binding sites in absolute terms, demonstrating that DOC quantity dominates over binding site quality in controlling seasonal REE mobility. This interpretation is consistent with Sonke (2006), who showed that humic complexation dominates dissolved REE speciation across pH 4–10 in average world river conditions, outcompeting carbonate even under alkaline conditions. Having established that DOC ligand abundance dominates REE complexation in both seasons, we note that inorganic ligand concentrations for the same campaigns (Lizotte et al., 2023) show limited capacity to alter this pattern. Phosphate concentrations do vary seasonally (winter: ∼0.04–0.06 µ mol L⁻¹; fall: up to 1.26 µmol L⁻¹) yet remain well below DOC in terms of REE complexation capacity. Silicate (∼55–71 µmol L⁻¹) and nitrate (∼4–11 µmol L⁻¹) show limited seasonal variation at concentrations unlikely to compete with DOC for REE complexation. The seasonal variations in organic matter quantity and quality observed in the Mackenzie Delta parallel those documented in the Lena Delta, another major Arctic river system (Juhls et al., 2020). In both systems, winter DOM is characterized by older, more degraded material with reduced overall fluxes. These seasonal organic matter dynamics are particularly important because coastal environments and deltas serve as crucial carbon reservoirs and processing hotspots (Bianchi and Allison, 2009), collectively accounting for approximately 80 % of total marine carbon burial (Hedges and Keil, 1995). This significant contribution results from high sedimentation rates and proximity to continental organic matter sources (Bianchi et al., 2018; Goñi et al., 2000). Recent modelling work on the Mackenzie River plume demonstrated that seasonal variations in terrestrial CDOM export affect coastal light attenuation, phytoplankton phenology, and sea-surface temperature, switching the coastal zone from a CO₂ sink to a source (Bertin et al., 2025). We note that logistical constraints resulted in partially non-overlapping sampling locations between seasons, particularly in the western transect. However, eastern sites were sampled at identical locations in both seasons and confirm the winter enrichment pattern, while the complete separation of seasonal samples in multivariate space (Fig. 6) demonstrates that seasonal processes dominate over site-specific spatial variability. Our findings that seasonal dynamics shape REE–organic matter interactions therefore have broader implications for understanding coupled carbon-trace element cycling in these biogeochemical hotspots.

4.3 Inorganic Controls on Dissolved REE Mobility

Temperature is a potential factor controlling REE release from the sediment to the porewater. Marginson et al. (2024) documented a significant correlation between lower temperatures and elevated REE concentrations in the George and Koroc Rivers, in their tributaries, and in thermokarst lakes in northern Québec. Their proposed mechanism invokes a higher partial pressure of CO₂ in colder waters (Young et al., 2025) which subsequently increases the formation and stability of REE–carbonate complexes in solution (Marginson et al., 2024). Our results are consistent with this previous report.

DOM ligands can mobilize REE from solid to dissolved phase in natural waters (Wen et al., 2024). However, the ability of DOM to mobilize REE can be hindered when dissolved Fe and Mn occupy DOM binding sites (Tang and Johannesson, 2003; Pourret et al., 2007b; Marsac et al., 2011; Neweshy et al., 2022). This competition varies seasonally. During fall, the season of high biological productivity and of high microbial respiration produced soluble Fe²⁺ via the reductive dissolution of Fe (oxy)hydroxides, which accumulates to levels that have been shown to displace REE from DOM binding sites (Neweshy et al., 2022). Indeed, Fe³⁺, Fe²⁺ and Mn²⁺ all have binding constants (pK_MHA for humic acid of respectively 0.8, 2.1 and 3.4) that are similar or higher than those of REE (average pK_MHA for humic acid of 1.53) (Tang and Johannesson, 2003), with Fe acting as one of the key competitors to REE in complexation by DOC (Takahashi et al., 1997). The increased competition for DOM ligands may thus partly explain the lower dissolved ∑REE concentrations in the fall, leaving in solution the poorly soluble free REE³⁺ ions or carbonate complexes (Tang and Johannesson, 2003; Johannesson et al., 1995; Pourret et al., 2007b).

5 Conclusions

Our study reveals the sensitivity of REE cycling to seasonal environmental changes in Arctic coastal systems. Sedimentary ∑REE concentrations remain remarkably stable across seasons, with positive correlations to Fe and Mn (oxy)hydroxides that strengthen in winter. In contrast, porewater ∑REE concentrations exhibit pronounced seasonal variability, with significantly higher values in winter than in fall. The nearly two-order-of-magnitude enrichment of dissolved ΣREEs in winter porewaters reflects the combined effects of seasonal changes in dissolved organic matter quantity and quality, with secondary contributions from carbonate speciation and Fe and Mn redox cycling. DOC ligand abundance emerges as the primary control: under ice-covered, hydrologically isolated conditions, winter porewaters accumulate DOC in sufficient quantity to generate approximately three times more absolute binding sites than fall DOM, despite its lower per-site affinity. This hierarchical interpretation is consistent with the PCA structure, in which DOC and associated optical indices (SUVA₂₅₄, HIX, BIX) load most strongly on PC1, which accounts for the dominant seasonal separation in the dataset. These seasonal dynamics have implications for understanding REE transport in Arctic coastal systems, which are experiencing intensifying riverine influence as warming accelerates (Emmerton et al., 2008b; Opsahl et al., 1999; Kipp et al., 2020).

Appendix A

Figure A1 Spatial interpolation of normalized average carbon concentrations in the Mackenzie River Delta for the overlying water (water column, WC, top panels), porewater (PW, middle panels) and sediment (Sed, bottom panels) for the winter (left) and fall (right) sampling seasons. Concentrations were normalized for dissolved and solid compartments with the compartment maximum, resulting in relative values between 0 (blue) and 1 (red). Grey circles show sampling sites. (Source: ESRI | Powered by Esri)

Figure A2 Spatial interpolation of normalized average Fe concentrations in the Mackenzie River Delta for the overlying water (water column, WC, top panels), porewater (PW, middle panels) and sediment (Sed, bottom panels) for the winter (left) and fall (right) sampling seasons. Concentrations were normalized for dissolved and solid compartments with the compartment maximum, resulting in relative values between 0 (blue) and 1 (red). Grey circles show sampling sites. (Source: ESRI | Powered by Esri)

Figure A3 Spatial interpolation of normalized average Mn concentrations in the Mackenzie River Delta for the overlying water (water column, WC, top panels), porewater (PW, middle panels) and sediment (Sed, bottom panels) for the winter (left) and fall (right) sampling seasons. Concentrations were normalized for dissolved and solid compartments with the compartment maximum, resulting in relative values between 0 (blue) and 1 (red). Grey circles show sampling sites. (Source: ESRI | Powered by Esri)

Table A1 Field blank concentrations (µmol L⁻¹) from porewater campaigns conducted the same year in the same laboratory using identical handling and preparation procedures.

Table A2 Measured concentrations and certified reference values for the solid-phase certified reference material MESS-4 (NRC Canada), with accuracy and precision for each analyte.

* MESS-4 certificate reports those values as information only without uncertainties.

Table A3 Measured concentrations and certified reference values for the aqueous-phase certified reference material SLRS-6 (NRC Canada), with accuracy and precision for each analyte.

Table A4 Principal component analysis summary of log-transformed seasonal dissolved Mn, Fe, REE, DOC and salinity data, SUVA₂₅₄ and temperature data for winter (n=53) and fall (n=36).

Table A5 Principal component variable loadings of seasonal dissolved Mn, Fe, REE, DOC and salinity data, SUVA₂₅₄ and temperature data for winter (n=53) and fall (n=36).

Table A6 MANOVA test for differences between winter (n=53) and fall (n=36).

Table A7 Univariate ANOVA of seasonal dissolved Mn, Fe, REE, DOC, SUVA₂₅₄ and salinity data for winter (n=53) and fall (n=36).