Holocene sea ice and paleoenvironment conditions in the Beaufort Sea (Canadian Arctic) reconstructed with lipid biomarkers

The Beaufort Sea region in the Canadian Arctic has undergone substantial sea ice loss in recent decades, primarily driven by anthropogenic climate warming. To place these changes within the context of natural climate variability, Holocene sea ice evolution and environmental conditions (sea surface temperature, salinity, terrestrial input) were reconstructed using lipid biomarkers (HBIs including IP₂₅, OH-GDGT, brGDGT, C_16:0 fatty acid, phytosterols) from two marine sediment cores collected from the Beaufort Shelf and slope, spanning the past 9.1 ka and 13.3 cal. kyr BP, respectively. The Early Holocene (11.7–8.2 ka) is characterized by relatively higher sea surface temperature, lower salinity and no spring/summer sea ice until 8.5 ka on the Beaufort Sea slope. Around 8.5 ka, a peak in organic matter content is linked to both increased terrestrial input and primary production and may indicate increased riverine input from the Mackenzie River and terrestrial matter input from coastal erosion. Following this period, terrestrial inputs decreased throughout the Mid-Holocene in both cores. A gradual increase in IP₂₅ and HBI-II concentrations aligns with relatively higher salinity, lower sea surface temperature and rising sea levels, and indicate the establishment of seasonal (spring) sea ice on the outer shelf around 7 ka and on the shelf around 5 ka. These patterns suggest an expansion of the sea ice cover beginning in the Mid-Holocene, influenced by decreasing summer insolation. During the Late Holocene (4.2–1 ka), permanent sea ice conditions are inferred on the slope with a peak during the Little Ice Age. After 1 ka, seasonal sea ice conditions on the slope are observed again, alongside an increase in salinity and terrestrial input, and variable primary productivity. Similar patterns of Holocene sea ice variability have been observed across other Arctic marginal seas, highlighting a consistent response to external climate forcing. Continued warming may drive the Beaufort Sea toward predominantly ice-free conditions, resembling those inferred for the Early Holocene.

1 Introduction

Sea ice is a critical component of the Arctic climate system, influencing ocean–atmosphere interactions, modulating surface albedo (Kashiwase et al., 2017), regulating heat fluxes (Lake, 1967), and influencing ecosystem structure through its control on light penetration and nutrient cycling (Lannuzel et al., 2020). Its high sensitivity to temperature makes it both a driver and indicator of Arctic climate change. Since the late 1970s, satellite observations have revealed a significant decline in Arctic sea ice extent, sparking interest in the mechanisms that govern sea ice variability over multiple timescales (Stroeve et al., 2012). The Canadian Beaufort Sea is a marginal sea of the western Arctic Ocean which exhibits strong seasonal and interannual variability in sea ice cover. Characterized by landfast ice on the shelf and mobile pack ice offshore, this region has experienced significant sea-ice loss in recent decades due to rising atmospheric and oceanic temperatures (Carmack et al., 2015; Comiso et al., 2008).

Understanding natural variability of sea-ice prior to the industrial era is critical for contextualizing recent trends. Throughout the Holocene, Arctic sea ice has responded to changes in orbital forcing, ocean circulation, and ice sheet dynamics (Park et al., 2018; Stein et al., 2017). During the Late glacial, abrupt climatic events such as Bølling-Allerød (∼ 14–12.8 ka) and Younger Dryas (∼ 12.8–11.7 ka), contributed to the instability of the Arctic cryosphere. In the Canadian Arctic, enhanced meltwater discharge and re-routing following the retreat of the Laurentide Ice Sheet (LIS) contributed to oceanographic shifts and transient cooling events (Broecker et al., 1989). Lipid biomarker records and climate simulations from the Arctic suggest reduced sea ice during the Early Holocene thermal maximum (11–6 ka in the Canadian Arctic) (Kaufman et al., 2004), followed by expansion through the Middle to Late Holocene, consistent with declining summer insolation (Stranne et al., 2014; Wu et al., 2020a). Numerous studies on Arctic sea ice variability have focused on offshore locations highlighting heterogeneity in sea-ice cover history and the importance of local currents (Belt et al., 2010; Detlef et al., 2023; Fahl and Stein, 2012; Hörner et al., 2016, 2018; Stein et al., 2017; Stein and Fahl, 2012; Vare et al., 2009; Wu et al., 2020a). However, these studies often neglected the spatial extent of sea ice cover toward the coast and the migration of the marginal ice zone.

Sea-ice cover can be reconstructed from microfossil and lipid biomarker evidence preserved in marine sediments. Remains of sea-ice organisms such as dinocysts (de Vernal et al., 2013) and diatoms, the latter producing a specific biomarker known as IP₂₅ (Belt et al., 2007), provide valuable records of past sea-ice conditions. This highly branched isoprenoid (HBI) and its isomer HBI diene (HBI-II) are used to trace the presence of spring sea-ice in modern and geological settings (e.g., Belt et al., 2010; Fahl and Stein, 2012; Hörner et al., 2016, 2018; Vare et al., 2009). However, because the absence of these two HBIs may reflect either a permanent sea-ice condition (due to the absence of light) or completely sea-ice free waters, the PIP₂₅ ratio was developed (Müller et al., 2011). This ratio includes a phytoplankton biomarker (typically dinosterol, brassicasterol or HBI-III) that represents open-water productivity. PIP₂₅ values have been used to distinguish between seasonal sea-ice (> 0.5) and permanent sea-ice cover (> 0.75).

Sea-ice cover is influenced by both the atmospheric and the oceanic forcing, including salinity and temperature, parameters that are often difficult to reconstruct in polar regions. Lipid biomarker and their (isotopic) ratios are a useful toolkit for this purpose. In particular, compound-specific hydrogen isotopes of phytoplankton biomarkers have shown promise for reconstructing past salinity (e.g., Lattaud et al., 2019; Sachs et al., 2018; Weiss et al., 2019). However, in the Arctic Ocean, the generally low abundances of lipid biomarker restrict this approach to the dominant lipid biomarkers, such as palmitic acid (C_16:0 fatty acid; Sachs et al., 2018). In contrast to salinity, several established proxies exist for reconstructing seawater temperature, including microfossil assemblages (e.g., dinocyst; Richerol et al., 2008), inorganic ratios (e.g., Mg $/$ Ca of foraminifera; Barrientos et al., 2018; Kristjánsdóttir et al., 2007) and lipid biomarkers (Ruan et al., 2017; Varma et al., 2024). Among lipid biomarker proxies for cold water (< 15 °C) environments, hydroxylated glycerol dialkyl glycerol tetraether (OH-GDGT) are particularly useful, with RI-OH^′ and TEX-OH identified as promising seawater temperature indices (Lü et al., 2015; Varma et al., 2024). Nevertheless, even the latest calibrations (Varma et al., 2024) reveal substantial uncertainty at low temperatures, partly due to the limited representation of polar core-top sediments, which are biased toward the European and Russian Arctic. These limitations highlight the need to further develop and test Arctic-specific proxies for both salinity and seawater temperature.

This study presents a multi-proxy reconstruction of Holocene sea ice and oceanographic variability from two sediment cores (PCB09, PCB11) collected from the Beaufort outer shelf and shelf slope. Lipid biomarkers, including highly branched isoprenoids (HBIs), glycerol dialkyl glycerol tetraethers (GDGTs), the hydrogen isotopic compositions of algal-derived fatty acids and terrestrial sterols, are used to reconstruct sea-ice cover, sea surface temperatures (SSTs), salinity and terrestrial organic matter input. Additionally, a set of surface sediments is used to assess the applicability and calibrate salinity and SST proxies in sediments of the Beaufort Sea.

The primary objectives are to (1) reconstruct past variations in sea-ice cover on the Beaufort Shelf throughout the Holocene, (2) explore the potential roles of insolation changes, meltwater input, and oceanic conditions in shaping regional sea-ice variability, and (3) place the Beaufort Sea record within a broader Arctic context to provide insights into past and present climate variability.

2 Materials and Methods

2.1 Study area

The study focuses on the Canadian Beaufort Sea, one of the marginal seas of the Arctic Ocean (Fig. 1), bounded by the glacially excavated Amundsen Gulf to the east, Mackenzie Trough to the west, and the Mackenzie River delta to the south (Carmack et al., 2004). The shelf is a large estuarine setting at the interface between the Arctic Ocean and the Mackenzie River (Omstedt et al., 1994) (Fig. 1). The Mackenzie River is a significant source of freshwater to the Beaufort Sea, with an annual water discharge of 316 km³ yr⁻¹ (Holmes et al., 2012) and is considered the largest Arctic river in terms of sediment flux (124–128 Mt yr⁻¹) (Stein et al., 2004). At the same time, permafrost coastal erosion adds another 8–9 Mt yr⁻¹ of sediment into the Beaufort Sea, including carbon and nutrients (Wegner et al., 2015). Surface water circulation in the Beaufort Sea is primarily characterized by the clockwise Beaufort Gyre, which drives offshore currents towards the west and traps the majority of the Arctic Ocean's freshwater in the Canada Basin (Serreze et al., 2006). There is also a eastward flowing shelf-break current at depths beneath 50 m, which transports Pacific Water (coming from the Bering Strait) along the slope (Pickart, 2004). Sea ice cover on the Beaufort Shelf north of the Mackenzie River Delta varies from year to year, but generally begins to form in mid-October, persisting until ice break up in April–May (Fig. S1 in the Supplement). During ice break up, an open water flaw lead occurs along the outer edge of the landfast ice establishing a spring marginal ice zone over the shelf.

Figure 1 Combined topographic and bathymetric map of the Beaufort Shelf region (Canadian Arctic) displaying the cores from this study as black squares (PCB09 and PCB11) and key records discussed in the text (ARA04C/37 from Wu et al., 2020a; JPC15/27 from Keigwin et al., 2018). The 2021 summer sea-ice extent is outlined by the orange dotted line, the winter sea ice extent follows the coastline (Meier et al., 2018). Map generated using Ocean Data View (Schlitzer, 2025).

2.2 Material

Two sediment cores were analyzed in the study. They were recovered as part of the Permafrost Carbon on the Beaufort Shelf (PeCaBeau) project during the 4th Leg of the 2021 CCGS Amundsen expedition (Bröder et al., 2022, Fig. 1). At station PCB09 (71.1° N, 135.1° W) at a water depth of 675 m on the Beaufort shelf slope, a piston core (PC, length of 420 cm) and multicore (MC, 30 cm) were retrieved (Fig. 1). At station PCB11 (70.6° N, 136.0° W) on the outer Beaufort shelf (74 m water depth) a giant gravity core (GGC, length of 290 cm) and MC (32 cm) were recovered (Fig. 1). PCB09 is found within the modern Atlantic bottom water mass, while PCB11 lies within the Pacific summer water (Fig. S2), the water masses were defined as in Matsuoka et al. (2012). The core tops (0–1 cm) from 22 multicores collected during PeCaBeau, were used to ground truth the hydrogen isotope ratio of C_16:0 fatty acid proxy for reconstructing salinity and test the applicability of the SST reconstructions (Fig. S1).

2.3 Methods

2.3.1 Core processing

During the PeCaBeau expedition, all cores were scanned shipboard on a Geotek multi-sensor core logger (MSCL). Bulk density and magnetic susceptibility were measured at a 1 cm downcore resolution on the piston and gravity cores, (Bröder et al., 2022). The cores remained unopened and were shipped to AWI Potsdam following the expedition. They were subsequently split in the fall/winter and the working halves sampled using 2 × 2 cm u-channels, before being cut into 1–2 cm thick slices which were frozen and freeze dried.

An ITRAX XRF-core scanner was used to measure relative elemental abundances at Stockholm University, Sweden. Measurements were performed on u-channel samples at a downcore resolution of 2 mm. Analyses were made using a Mo tube at a voltage of 55 kV, a current of 50 mA and an exposure time of 20 s. Here we present the ratios of Ca $/$ Ti, reflecting detrital carbon inputs regionally elevated by meltwater delivery from either the Mackenzie or Amundsen Gulf during deglaciation (Klotsko et al., 2019; Swärd et al., 2022; Wu et al., 2020a). Zr $/$ Rb was used as a proxy for grain size variations as Zr content is elevated in coarse minerals, while Rb is associated with clay minerals (Wu et al., 2020b) and Br $/$ Cl as a proxy for marine organic matter as Br usually correlates well with OC content (Wang et al., 2019).

2.3.2 Age model

The chronology of the piston and gravity cores (Fig. 2) were determined by ¹⁴C dating of foraminifera (n = 13, PCB09) and bivalve shells (n = 7, PCB11) (Table S1). The MSCL data was used to stratigraphically correlate PCB09 with JPC15/27 (Keigwin et al., 2018) allowing us to integrate existing radiocarbon ages (n = 8) from this record with our new data (n = 5) (Fig. S3).

Figure 2 Core description for (A) PCB09 and (B) PCB11 presenting lithostratigraphy, bulk density, magnetic susceptibility, X-Ray Fluorescence (XRF) results including Ca $/$ Ti, Zr $/$ Rb and Br $/$ Cl ratios. Total organic carbon content (%) and abundance of N. pachyderma (N. p., number per gram of sediment) as well as the age models of PCB09 and PCB11 generated using the Bacon Rpackage (Blaauw and Christen, 2011). The red bands illustrate the 95 % confidence level around the modelled median age (blue line). Symbols indicate the calibrated age of the dated material before age-modelling. The errors on these are lower than the symbol size. The numbers are the median modelled ages and 95 % error at the location of each sample. Blue radiocarbon ages originated from nearby core HLY13-15JPC (Keigwin et al., 2018).

Bivalve shells were either picked from the split cores when sampling, or later from the freeze-dried sediments. Foraminifera were picked from the > 45 µm fraction of the wet-sieved samples following organic extractions. Foraminifera samples consisted of either planktonic (Neogloboquadrina pachyderma), benthic, or a combination of both in horizons when specimens were extremely rare. Care was given to pick well preserved foraminifera to avoid age bias (Wollenburg et al., 2023). Foraminifera and mollusk samples were prepared for Accelerator Mass Spectrometry (AMS) analyses at the Laboratory for Ion Beam Physics at ETHZ using procedures described in Missiaen et al. (2020) which include sieving and acid cleaning to remove impurities from the shells.

Radiocarbon-based age models were generated using the BACON package in R (Blaauw and Christen, 2011) and the Marine20 calibration curve (Heaton et al., 2020). A reservoir age of 330 ± 41 years was applied to the Holocene-age mollusc samples in PCB11 as determined by West et al. (2022) for Pacific waters entering the Arctic Ocean in the Chukchi Sea. For PCB09 we applied the approach used for JPC15/27 (Keigwin et al., 2018) and for ARA04C/37 (Wu et al., 2020a) but updated for Marine20 as described by Lin et al. (2025). A reservoir correction of −150 ± 100 years was applied to Holocene planktic foraminifera, and a larger reservoir correction (50 ± 100 years) for the bottom 4 samples (Table S1) that fall within the Younger Dryas. In our age model we also incorporate samples of benthic foraminifera that were calibrated using a reservoir correction of 206 ± 67 years, determined by West et al. (2022) for Atlantic waters near the Chukchi Sea. Samples containing mixed planktic and benthic foraminifera were calibrated using an average of these values (28 ± 85 years).

2.3.3 Bulk organic matter

For the determination of total organic carbon (TOC) content and stable carbon isotope composition (δ¹³C) at the University of Basel, about 12 mg of freeze-dried sediment was weighed into each silver capsule and 1–2 drops of distilled water were added. The samples were exposed to fuming hydrochloric acid (HCl, 37 %) in a desiccator for 24 h to remove inorganic carbon. Samples were dried (48 h, 50 °C) and analyzed using an elemental analyser coupled to an isotope mass spectrometer (Sercon, Integra 2). The standards used to calculate TOC was Ethylenediaminetetraacetic acid (EDTA, Sigma Aldrich) and for δ¹³C were USGS40 (−26.389 ± 0.042 ‰, IAEA), USGS64 (−40.81 ± 0.04, IAEA), and USGS65 (−20.29 ± 0.04, IAEA). The analytical precision, defined as the standard deviation of the measurement of the USGS standards for the δ¹³C sequence was ±0.03 ‰.

2.3.4 Lipid biomarkers

Lipid biomarkers were analysed from 42 samples (every 10 cm for the first 143 cm, then every 20 cm) for PCB09 and 21 for PCB11 (every 20 cm). Core top samples from the MC's were also analysed for lipid biomarkers. For each sample, 5 g of homogenized freeze-dried sediment was extracted using an Energy Dispersive Guided Extraction (EDGE) following (Lattaud et al., 2021a). Briefly, after extraction with dichloromethane (DCM) : methanol (MeOH) (9 : 1, $v/v$), the total lipid extracts (TLE) were saponified at 70 °C for 2 h (with KOH in MeOH at 0.5 M). The neutral phase was collected by liquid-liquid extraction with 10 mL of hexane, three times. The leftover TLE was acidified to pH 2 and the acid phase was recovered by liquid-liquid extraction adding 10 mL hexane : DCM (4 : 1, $v/v$), three times. The acid compounds were methylated by adding MeOH : HCl (95 : 5, $v/v$) and heated at 70 °C overnight. The methylated fatty acids were recovered by liquid-liquid extraction (three times) with 10 mL hexane : DCM (4 : 1, $v/v$). Internal standards were added to the neutral fraction prior to silica chromatography: 7-hexylnonadecane (7-HND, provided by Simon Belt), 9-octylheptadec-8-ene (9-OHD, provided by Simon Belt), C₂₂ 5,16-diol (Interbioscreen), C_36:0 alkane (Sigma Aldrich) and C₄₆ GDGT-like compound (C₄₆ glycerol trialkyl glycerol tetraether, GTGT; Huguet et al., 2006). The neutral phase was separated into three fractions (F1, F2, and F3) through silica column (combusted and deactivated 1 %) using hexane : DCM (9 : 1, $v/v$), DCM, and DCM : MeOH (1 : 1, $v/v$), respectively.

The F1 containing HBIs was analyzed on a GC-MS (Agilent 7890-5977A) operating in Selective Ion Monitoring (SIM) mode at the Institute of Polar Sciences (ISP), Bologna, Italy, following (Belt et al., 2014). The column used was a J&W DB5-MS (length 30 m, id 250 µm, 0.25 µm thickness). Integrations were done in SIM mode for IP₂₅ ($m/z$ = 350), HBI II ($m/z$ = 348) and HBI III and HBI IV ($m/z$ = 346). Concentrations of IP₂₅ were corrected for $m/z$ 348 influence (4 %) and instrumental response factor. 9-OHD ($m/z$ 350) was used to quantify HBIs. A reference sediment containing known amount of IP₂₅ was run in parallel to correct IP₂₅ concentration.

F3 was split in two with one aliquot that was filtered using a polytetrafluoroethylene filter (PTFE, 0.45 µm pore size) and analyzed for GDGTs with high performance liquid chromatography (LC)/atmospheric pressure chemical ionization–MS on an Agilent 1260 Infinity series LC-MS according to Hopmans et al. (2016) with the adaptation of Lattaud et al. (2021b). GDGTs were quantified using the C₄₆ GTGT internal standard assuming the same response factor.

The other F3 aliquot was silylated with bis(trimethylsilyl)trifluoroacetamid (BSTFA) (70 °C 30 min) and analysed at the ISP for sterol concentration (brassicasterol, stigmasterol, β-sitosterol, campesterol) on a GC-MS. The C₂₂ 5,16-diol is used to quantify sterols. Specific $m/z$ ratios have been extracted from chromatograms in order to identify each biomarker according to their respective mass spectra, sterols were quantified on the total ion current.

Lipid δ²H values were analyzed by GC coupled to an isotope ratio mass spectrometer (IRMS) on all acid fractions having adequate compound abundance. Samples were analyzed using splitless injection with a split/splitless inlet at 280 °C and a Restek Rtx-5MS GC column (30 m × 0.25 mm × 0.25 µm) with helium carrier gas at 1.4 mL min⁻¹. The GC oven was held at 60 °C for 1.5 min, ramped to 140 °C at 15 °C min⁻¹, then to 325 °C at 4 °C min⁻¹, and held for 15 min. Column effluent was pyrolyzed at 1420 °C, and δ²H values were measured on a Thermo Delta V Plus IRMS. The H${}_{\mathrm{3}}^{+}$ factor was evaluated with each measurement sequence to confirm stability. Values were always lower than 3 ppm mV⁻¹. Reference standards with known isotopic compositions (Mix A7, USGS71, C_30:0 FAME; provided by Arndt Schimmelmann, Indiana University, USA) were analyzed alongside samples to normalize values to the Vienna Standard Mean Ocean Water-Standard Light Antarctic Precipitation (VSMOW-SLAP) scale. Standards were injected at a range of concentrations so that peak size effects could be assessed and corrected for. Quality control samples with known δ²H values were measured as unknowns to check precision and accuracy (C_16:0 FAME in mix F8-40, C_30:0 FAME; Arndt Shimmelmann), which were 4.2 ‰ or better, and 1.0 ‰ or better, respectively (n = 13–16). Final fatty acid δ²H values of C_16:0 were corrected for added hydrogen during methylation following Eq. (1).

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{C}\mathrm{16}:\mathrm{0}}={\displaystyle \frac{\begin{array}{c}\left({\mathrm{nH}}_{\mathrm{FAME}}+{\mathrm{nH}}_{{\mathrm{CH}}_{\mathrm{3}}}\right)\times {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FAME}}\mathrm{measured}\\ -{\mathrm{nH}}_{{\mathrm{CH}}_{\mathrm{3}}}\times {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{CH}}_{\mathrm{3}}}\end{array}}{{\mathrm{nH}}_{\mathrm{FAME}}}}\end{array}$$

Where nH${}_{{\mathrm{CH}}_{\mathrm{3}}}$ = 3, nH_FAME=32.

2.3.5 Biomarker ratios and environmental indicator reconstructions (sea ice, salinity, SSTs)

In order to describe sea ice variability in the Holocene, the PIP₂₅ index is used (Müller et al., 2011). The PIP₂₅ index (Eq. 2) uses additional phytoplankton biomarkers (i.e. brassicasterol, dinosterol, and HBI-III) which indicate open water conditions to compare with the abundance of IP₂₅ (Belt et al., 2007):

$$\begin{array}{}\text{(2)}& {\mathrm{PIP}}_{\mathrm{25}}={\displaystyle \frac{{\mathrm{IP}}_{\mathrm{25}}}{\left[{\mathrm{IP}}_{\mathrm{25}}\right]+\left[\mathrm{Phytoplankton}\mathrm{biomarker}\right]\times c}}\end{array}$$

HBI-III was used in this study (Belt et al., 2015; Kolling et al., 2020; Köseoğlu et al., 2018; Smik et al., 2016) as a reference for pelagic phytoplankton to derive P_IIIIP₂₅ index (afterward called PIP₂₅). Dinosterol was not detected in the samples which is common in the Beaufort Sea (Fu et al., 2025), and brassicasterol has been shown to derive mainly from terrestrial input in the region (Wu et al., 2020a). The c value represents the ratio of the mean concentration of IP₂₅ over the mean concentration of HBI-III of all samples for each core (4.79 for PCB09 and 5.55 for PCB11).

Surface salinity was reconstructed using the calibration between δ²H of C_16:0 fatty acid (palmitic acid) and salinity from the test study of Sachs et al. (2018) (Eq. 3), adding data from Allan et al. (2023) and after testing surface sediments from multicores from the region (Fig. S4):

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{PA}}=\mathrm{4.22}(\pm \mathrm{0.6})\times S-\mathrm{338}(\pm \mathrm{15})\end{array}$$

where S is salinity in practical salinity units (psu). Based on the known calibration errors (4 ‰ for the δ²H measurement), reconstructed salinity should have an associated error of ±7 psu.

To reconstruct SST, hydroxylated GDGTs (OH-GDGTs) were used as the hydroxyl group in these GDGTs is suggested to be an adaptation feature to regulate permeability in cold waters (Liu et al., 2012). In this study, the RI-OH^′ (Eq. 4) and TEX-OH (Eq. 5) indexes were calculated:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\text{RI-OH}}^{\prime }={\displaystyle \frac{\left[\text{OH-GDGT-1}\right]+\mathrm{2}\times \left[\text{OH-GDGT-2}\right]}{\begin{array}{c}\left[\text{OH-GDGT-0}\right]+\left[\text{OH-GDGT-1}\right]\\ +\left[\text{OH-GDGT-2}\right]\end{array}}}\text{(5)}& {\displaystyle }\text{TEX-OH}={\displaystyle \frac{\text{GDGT-2}+\text{GDGT-3}+\text{Cren isomer}}{\begin{array}{c}\mathrm{\text{GDGT-2}}+\text{GDGT-3}+\text{Cren isomer}\\ +\text{OH-GDGT-0}+\text{GDGT-1}\end{array}}}\end{array}$$

For the conversion from RI-OH^′ and TEX-OH to SST, the recent calibration of Varma et al. (2024) is used (Eqs. 6 and 7):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\text{RI-OH}}^{\prime }=\mathrm{0.04}\times \mathrm{SST}+\mathrm{0.003}\text{(7)}& {\displaystyle }\text{TEX-OH}=\mathrm{0.021}\times \mathrm{SST}+\mathrm{0.08}\end{array}$$

Several organic proxies have been used to interpret terrestrial organic matter input such as branched glycerol dialkyl glycerol tetraether (brGDGTs), long chain n-alkanes, and plant sterols. BrGDGTs are membrane lipids synthesized by bacteria and are known to be ubiquitous in terrestrial environments (Schouten et al., 2013). The BIT index (Hopmans et al., 2004) (Eq. 8) is a common indicator of terrestrial input into the marine realm:

$$\begin{array}{}\text{(8)}& \mathrm{BIT}={\displaystyle \frac{\left[\mathrm{BrGDGT}-\mathrm{Ia}+\mathrm{IIa}+\mathrm{IIIa}\right]}{\left[\mathrm{BrGDGT}-\mathrm{Ia}+\mathrm{IIa}+\mathrm{IIIa}\right]+\mathrm{Crenarcheol}}}\end{array}$$

In addition, in situ marine production of brGDGT can occur in coastal sediments between 50 and 300 m water depth (Peterse et al., 2009; Sinninghe Damsté, 2016). To assess the potential for brGDGT to be produced in situ we calculated the number of cyclopentane rings present in the tetraether brGDGT (#ring_tetra) (Eq. 9):

$$\begin{array}{}\text{(9)}& \mathit{\#}{\mathrm{ring}}_{\mathrm{tetra}}={\displaystyle \frac{(\mathrm{Ib}+\mathrm{2}\times \mathrm{Ic})}{(\mathrm{Ia}+\mathrm{Ib}+\mathrm{Ic})}}\end{array}$$

2.3.6 Micropaleontology

Extracted sediments were wet sieved using a 45 µm mesh. The > 45 µm fraction was dried in the oven (40 °C) and picked for foraminifera using a stereoscopic microscope. Planktonic foraminifera species are identified (https://www.mikrotax.org/pforams/, last access: 1 August 2023) using the morphological descriptions compiled in Microtax and counted for each sample.

3 Results

3.1 Core chronology, lithostratigraphy and bulk organic matter

The age models of cores PCB09 (Fig. 2a) and PCB11 (Fig. 2b) show that they cover the last 13350 ± 690 and 9120 ± 720 cal. yr BP, respectively. PCB11 had a mean sedimentation rate of 35 ± 10 cm kyr⁻¹ with slightly higher sedimentation rates in the Late Holocene (< 4 ka). Conversely, PCB09 had an average sedimentation rate of 50 ± 27 cm kyr⁻¹, with substantially higher sedimentation rates before the Late Holocene.

The upper 300 cm of PCB09 (0–11 ka) display a gradual downcore increase in bulk density, reflecting porosity loss in largely homogenous silty-clay sediments (Fig. 2a). Below 300 cm, slightly higher variability in the bulk density, elevated magnetic susceptibility and higher Zr $/$ Rb all point towards a transition to slightly coarser-grained sediments. Ca $/$ Ti tended to increase downcore becoming more variable below 105 cm (7.6 ± 0.6 cal yr BP). There was a notable stepwise increase in Ca $/$ Ti at 345 cm (11.7 ± 0.4 ka, Fig. 2a). Discrete peaks in Zr $/$ Rb and bulk density, indicative of sediment coarsening, co-occurred with elevated Ca $/$ Ti ratios at depths of 276 cm (10.8 ± 0.4 cal. yr BP), 345–352 cm (11.8 ± 0.4 cal. yr BP) and 402 cm (12.8 ± 0.4 cal. yr BP).

A similar pattern is seen in PCB11, where below 240 cm (7.4 ± 0.6 cal. yr BP) there was an abrupt increase in bulk density, magnetic susceptibility, Zr $/$ Rb and Ca $/$ Ti. TOC concentrations and the Br $/$ Cl ratio (which mirrors small scale changes in the TOC) also decreased notably through this interval (Figs. 2b, 3a). This lithologic transition post-dates the Late glacial and Early Holocene detrital carbonate inputs in cores recovered from the Beaufort Sea slope (Klotsko et al., 2019). It is likely that this coarser basal facies is related to the inundation of the shelf during transgression.

Figure 3 Bulk characteristics and biomarker concentrations (in ng g${}_{\mathrm{sediment}}^{-\mathrm{1}}$) for core PCB09 (brown circles) and PCB11 (red squares) with (a) total Organic Carbon (TOC), (b) δ¹³C, (c) IP₂₅, (d) HBI-II, (e) HBI III, (f) HBI IV, (g) isoprenoid GDGTs (isoGDGT), (h) hydroxylated GDGTs (OH-GDGT), (i) branched GDGTs (brGDGT) and (j) terrestrial sterols (sum of brassicasterol, stigmasterol, β-sitosterol, campesterol).

Bulk sediment δ¹³C of PCB09 was lowest in the Early Holocene at approximately −26.3 ‰ until 8.7 ± 0.4 ka, before increasing during the Middle and Late Holocene to −24.7 ‰ (Fig. 3b). The trend in δ¹³C in PCB11 is similar to PCB09 showing a steady increase over time from −26.5 ‰ to −25.8 ‰.

3.2 Lipid biomarkers

IP₂₅ and HBI II (C_25:2) concentrations were generally low (< 2 ng g⁻¹) in the Early Holocene (Fig. 3c, d). IP₂₅ in both cores increased throughout the Middle to Late Holocene. During the Late Holocene, IP₂₅ and HBI II concentrations dropped in PCB09 around 1.8 ± 1.5 ka. Concentrations of both biomarkers were higher in PCB11 than in PCB09 after 3 ka, reaching modern values of 40 and 150 ng g⁻¹ (IP₂₅ and HBI II). PIP₂₅ values in both cores increased from the Early to the Mid-Holocene (Fig. 4a).

Figure 4 Reconstructed environmental parameters for PCB09 (brown circles) and PCB11 (red squares) (a) P_IIIIP₂₅ (PIP₂₅ calculated with HBI III as phytoplankton biomarker), (b) δ²H of C_16:0 fatty acid and corresponding reconstructed salinity (Sachs et al., 2018), (c) BIT index (Hopmans et al., 2004), (d) SST derived from TEX-OH (Varma et al., 2024) shaded areas indicate the uncertainties of the reconstruction, modern (1995–2004) annual and summer SSTs were extracted from (Locarnini et al., 2024) (e), global sea level estimates derived from Lambeck et al. (2014) and (f) 21 June insolation at 70° N (Laskar et al., 2004). (g) Illustrative examples of plaeoshorelines at 8 and 4 ka compared to the modern. These were generated by adjusting the sea-level using the modern bathymetry portrayed in IBCAO V. 5 (Jakobsson et al., 2024). Relative sea-level adjustments were taken from ICE 6G_C (Peltier et al., 2015) for the grid cell encompassing the position of PCB11. The sea-level adjustments for this location were 46 m at 8 ka and 12 m at 4 ka (Fig. S7).

HBI III (C_25:3) and HBI IV (C_25:4) were low in both cores with values below 8 ng g⁻¹ (Fig. 3e, f). Concentrations were higher in PCB11 than in PCB09 after 4 ka.

The concentration of isoGDGTs and OH-GDGTs followed a similar pattern throughout the Holocene (Fig. 3g, h). IsoGDGTs and OH-GDGT concentrations in PCB09 were stable during the Early Holocene at around 400 and 25 ng g⁻¹, respectively. At around 8.5 ka, the isoGDGTs and OH-GDGT amounts doubled. Throughout the Mid-Holocene, isoGDGTs and OH-GDGT concentrations were variable but above 500 ng g⁻¹. A drop in PCB09 to almost below detection limits occurred between 1–1.5 ka. IsoGDGTs and OH-GDGTs in PCB11 showed a steady increase from around 100 and 10 ng g⁻¹, respectively, in the Early Holocene to > 500 and > 50 ng g⁻¹ (Fig. 3g, h).

BrGDGTs concentrations in PCB09 were below 100 ng g⁻¹ throughout the cores except for peaks during the Early and Mid-Holocene at 11.2 ± 0.3, 8.2 ± 0.5 and 5.7 ± 0.5 ka, the latter was also seen in PCB11 (albeit concentrations were higher in PCB11) (Fig. 3i). BIT values varied from 0.1 to 0.4 in PCB09 and from 0.1 and 0.5 in PCB11. The #ring_tetra values were always < 0.7 for both cores. Terrestrial sterol concentrations in PCB09 were relatively stable throughout the core except for short-lived peaks at 9.5 ± 0.4, 8.2 ± 0.5 and 5.0 ± 0.9 ka (Fig. 3j). In PCB11, the concentration remained stable throughout the core after an initial increase at 6.9 ± 0.6 ka and a peak in the surface sediment.

Although biomarker concentrations expressed per gram of sediment may be influenced by mineral dilution or post-depositional degradation, TOC contents in our cores varied only slightly (0.9 %–1.3 % in PCB09 and 1.1 %–1.3 % in PCB11; Fig. 3a). Such limited variability indicates that TOC normalization would not substantially alter the observed trends. The distinct downcore patterns among lipid biomarker classes therefore likely reflect differences in source input and preservation rather than a uniform effect of degradation or dilution.

3.3 Salinity, SST and terrestrial input inferred from biomarker ratios

Surface sediments δ²H C_16:0 values from the Beaufort Sea (Figs. S1, S4) range from −275 ‰ to −200 ‰, comparable with the values obtained by the preliminary study of Sachs et al. (2018). δ²H C_16:0 values of all sediments correlate with summer salinity (r² = 0.63, p < 0.001) and the calibration equation is the same as the one obtained by Sachs et al. (2018). It is to be noted that the uncertainty associated with the analysis and calibration reaches 7 psu, which is quite large for salinity changes over glacial-deglacial timescales. This is to the contrary to what Allan et al. (2023) and Wu et al. (2025) observed in a set of surface sediments in Baffin Bay and a downcore record of the Beaufort Sea, where no relationship with salinity is observed. This contrast for the surface sediments could come from the different environment as Baffin Bay is a much more enclosed basin compared to the Beaufort Sea, or that δ²H C_16:0 values encompassed a too small range of salinity (28–34 psu). The downcore Beaufort Sea record yielded no notable salinity response of the δ²H C_16:0 when large freshwater influence were expected, which was attributed to the production of C_16:0 by heterotrophs (Wu et al., 2025).

Sea surface salinity inferred from δ²H C_16:0 values at PCB09 increased from 27 ± 7 psu during the Early Holocene to 30 ± 7 psu during the Mid-Holocene, and remained stable until 3 ka, increasing to 32 ± 7 psu during the Late Holocene (Fig. 4b). Reconstructed salinity at PCB11 was more stable during the Mid-Holocene and Late Holocene (28 ± 7 psu) until an increase to 30 ± 7 psu in the last centuries. (Fig. 4b). Although uncertainties associated with reconstructed salinity are large (±7 psu) the salinity trend between both locations agree with modern observations of lower salinities at PCB11 than further offshore at PCB09 (Fig. S2).

Two different sets of SSTs were reconstructed using the OH-GDGT only (RI-OH^′) or a combination of OH- and isoGDGT (TEX-OH) (Figs. 4d, S4, S5a, b). SSTs were only reconstructed when the BIT index was below 0.3 (Fig. 4c) as both calibrations are sensitive to terrestrial input (Varma et al., 2025). RI-OH^′ in the surface sediments varies from 0.05 to 0.17 while TEX-OH varies from 0.08 to 0.32. Both indexes plot in the global calibration curves from (Varma et al., 2024) and the reconstructed SST varies from 0.9 to 4.0 and −0.1 to 11.6 °C, respectively. TEX-OH reconstructed SSTs in PCB09 varied between 7 ± 2.6 °C in the Early Holocene, remained stable during the Mid-Holocene (∼ 3 ± 2.6 °C), decreased to 0 ± 2.6 °C between 1–1.5 ka and after which they increase to 5 ± 2.6 °C, close to modern summer SST (Locarnini et al., 2024). The reconstructed SST in the multicore is −1 ± 2.6 °C, similar, within the calibration error, to modern annual SST (Locarnini et al., 2024) (Fig. 4d). RI-OH^′ reconstructed SSTs in PCB09 produce high temperatures (∼ 12 ± 2.6 °C) between 7 and 10 ka, maybe due to non-thermal impact on the production of OH-GDGT (Harning and Sepúlveda, 2024). For both cores, reconstructed SST using RI-OH^′ is stable around 3 °C (Fig. S5b). PCB11 TEX-OH reconstructed SSTs are stable during the Early to Mid-Holocene (∼ 5.0 ± 2.6 °C). However, after 1 ka the reconstruction produces unrealistically large variation (5–15 °C, Fig. S5a).

The BIT index showed a steady decrease in PCB11 throughout the Holocene and until 3 ka in PCB09 (Fig. 4c). In PCB09, this decrease was interrupted at 9 and at 1.5 ka with BIT index values reaching 0.3 and 0.4, respectively. The increase at 9 ka was likely due to a relative decrease in crenarchaeol concentration (Fig. S5a) whereas the 1.5 ka increase was likely due to a decrease in brGDGT concentration (Fig. 4i).

3.4 Micropaleontology

Almost all of the planktonic foraminifera (99 %–100 % in abundance relative to other species) observed in PCB09 are N. pachyderma, formerly N. pachyderma sinistra. This is expected since this species has been found to dominate polar water masses (e.g. Eynaud 2011; Moller, Schulz, and Kucera 2013). Planktonic foraminifera are mostly absent in PCB11, consistent with data from plankton tows indicating that planktic foraminifera are rare on the Canadian shelf where surface waters are influenced by Mackenzie River discharge (Vilks, 1989). In PCB09, the foraminiferal shells appear white and fragmented in sections with abundant light-colored and sand-sized ice-rafted debris and other detrital materials (Fig. S6). Foraminifera are more abundant in samples that have relatively more mud aggregates than sand-sized debris (Fig. 2a, b). There is almost zero accumulation rate (per mm yr⁻¹) of N. pachyderma within the shelf slope from 10 ka.

4 Discussion

This study aims to reconstruct Holocene paleoenvironmental conditions in the southeastern Beaufort Sea focusing on spatial variability between the shelf slope (> 500 m water depth) and the outer shelf (< 100 m water depth). By analyzing the abundance and ratios of sea ice biomarkers (IP₂₅, HBI II), phytoplankton and heterotrophic archaeal productivity markers (HBI III, HBI IV, iso- and OH-GDGT), terrestrial inputs (brGDGTs, terrestrial sterols), and reconstructed environmental indicators (salinity, SST) this study aim to highlight spatial environmental difference between a shallow (PCB11) and deep (PCB09) site. In the following sections, we interpret biomarker records in a chronological framework, highlighting the dynamic relationship between freshwater inputs, ocean circulation, and sea ice conditions.

4.1 Late glacial to Early Holocene (12–8.2 ka)

The Late glacial to Early Holocene is only recorded at the shelf slope location (PCB09, Fig. 1). This period is characterized by low concentrations of sea ice biomarkers (Fig. 3a, b) resulting in low PIP₂₅ values (Fig. 4a). The low concentration suggests that this area had intermittent sea ice coverage during the Late glacial to Early Holocene, but the presence of HBI III and HBI IV (Fig. 3e, f) indicate that the region was only under seasonal ice cover until spring allowing late spring/summer open-water diatom primary production (Belt et al., 2015). Heterotrophic production in the shelf slope region during this period is relatively low (as suggested by the presence of ammonium oxidizer Thaumarchaea-derived isoGDGTs; Schouten et al., 2013; Fig. 3g, h) but increased and peaked at 8.2 ± 0.5 ka. During 12–8.2 ka, SSTs were elevated in comparison with the rest of the Holocene (Fig. 4d) which coincided with peak 21 June insolation (Fig. 4f) (Clemens et al., 2010; Laskar et al., 2004). The warmer surface waters might have inhibited the development of sea ice over the Beaufort Shelf.

During the Late glacial to Early Holocene, large freshwater inputs to the Beaufort Shelf, inferred from the low reconstructed salinity (Fig. 4b) likely originated from the decaying Laurentide Ice Sheet. Such water masses derived from drainage regions that had undergone minimal weathering would have released low amounts of nutrients. The influx of low-salinity freshwater may have intensified salinity-driven stratification on the shelf, reducing the upwelling of nutrient-rich saline Pacific waters to the surface which also limited nutrient availability. This stratification and lower nutrient availability likely limited primary productivity and the presence of ammonia-oxidizers on the Beaufort Shelf. It is important to note that sea level on the Beaufort Shelf was > 60 m lower in the Early Holocene than it is today (Fig. 4e). Implying that between 10–12 ka, the Beaufort Sea was a shallow estuarine environment (Fig. 4g; Hill et al., 1993).

The concentration of brGDGTs and terrestrial sterols in the shelf slope location during the Early Holocene peaked at 11.3 ± 0.3 and 8.2 ± 0.5 ka (Fig. 3i, j), which agrees with meltwater outflow from the LIS and freshly deglaciated surfaces as seen in nearby cores (Klotsko et al., 2019; Wu et al., 2020a). Additionally, increased freshwater input may have transported more detrital calcium (Ca), as indicated by elevated Ca $/$ Ti ratios (Fig. 2a), which could have enhanced the preservation of foraminifera by buffering the water column and limiting carbonate dissolution, in sediments along the shelf slope. Murton et al. (2010) used optically stimulated luminescence (OSL) dating to identify two major meltwater pulses through the Mackenzie River system between 13 and 11.7 ka and between 11.7 and 9.3 ka. This timing is supported by sedimentary and isotopic records from the Beaufort Sea indicating a major Lake Agassiz flood route through the Mackenzie system (Keigwin et al., 2018; Klotsko et al., 2019). These meltwater events coincide with events (11.3 ± 0.3, 8.2 ± 0.5 ka) in the biomarker records from this study (Fig. 3), and one event at 10.1 ± 0.4 ka is recorded in the reconstructed salinity (Fig. 4b), suggesting enhanced freshwater forcing impacted ocean circulation and increased sea ice extent. The massive meltwater discharge from the LIS (at least ∼ 9000 km³) into its surrounding oceans has been the major cause for eustatic sea level rise from 10 to 6 ka (Moran and Bryson, 1969).

4.2 Middle to Late Holocene (8.2–0 ka)

After 8.2 ka, environmental conditions at the shelf slope (PCB09) and the outer shelf (PCB11) diverged markedly, reflecting their contrasting depositional and oceanographic settings. At PCB09, SSTs cooled from ∼ 6 to 3 °C (Fig. 4d), while steadily increasing sea-ice biomarker concentrations led to PIP₂₅ > 0.5 by 7–6 ka, indicating the establishment of stable ice-edge or polynya conditions at the slope. The higher salinity (Fig. 4b) and greater distance from the coast at this site suggest enhanced influence of offshore Pacific-derived waters and reduced terrestrial input.

In contrast, at the outer shelf site PCB11, sea-ice biomarkers also increased after 8.5 ka but were accompanied by persistently high concentrations of open-water diatom markers, implying continued seasonal sea-ice and productive flaw-lead conditions, i.e., an open-water or newly formed sea-ice zone between landfast ice and sea ice. The flaw-lead today occurs about 80 km away from shore (Fig. S1) (Carmack et al., 2004). The proximity of PCB11 to the coastline before 6 ka (Fig. 4g) likely favored landfast-ice diatom assemblages and greater sensitivity to freshwater discharge. Only after 4 ± 0.5 ka did PCB11 reach PIP₂₅ values comparable to the slope, indicating a delayed transition to stable seasonal sea ice, approximately 2 kyr later than at PCB09. Thus, while both sites record a Middle-Holocene trend toward increasing sea-ice cover, the slope experienced earlier stabilization and reduced productivity linked to offshore cooling and stratification, whereas the outer shelf remained a dynamic, seasonally open-water environment likely sustained by coastal flaw-lead formation, and strong riverine influence.

In the Late Holocene (4.2 ka to present), the two Beaufort sites again show distinct sea-ice and productivity trajectories reflecting their different hydrographic settings. At the slope site, perennial sea ice cover developed after 3 ka as indicated by higher PIP₂₅ values (> 0.8), increased reconstructed salinity (Fig. 4a, b), as well as a decline in open-water diatom biomarkers (Fig. 3e, f). A sharp decrease in PIP₂₅ and open-water diatom biomarkers at 1.5 ± 1.5 ka suggests a shift to year-round ice cover coincident with low SST (∼ 0 °C, Fig. 4d). The parallel reduction in heterotrophic archaeal production may reflect strengthened water-column stratification or reduced nutrient supply, consistent with restricted shelf-break upwelling (Schulze and Pickart, 2012). These changes are broadly consistent with the timing of the regional cooling associated with the Little Ice Age (Mann et al., 2009), although the resolution of the biomarker record does not allow precise attribution to centennial-scale events. In contrast, on the outer shelf, seasonal sea-ice conditions persisted longer and sea ice cover expanded gradually and became well established after about 2 ± 0.6 ka. Even as sea-ice biomarkers increased, open-water diatom markers remained relatively abundant, implying continued flaw-lead or marginal-ice-zone productivity sustained by intermittent open-water formation and coastal influence.

By the last 0.4 ka, both sites reached modern configurations, but with clear spatial differences: PCB09 records persistently thick, multi-year ice and limited primary production, whereas PCB11 retains higher concentrations of both sea-ice and open-water diatom biomarkers, consistent with modern flaw-lead dynamics observed near the outer shelf (Fig. S1).

4.3 Comparison with other Arctic marginal seas

Previous studies using IP₂₅ to reconstruct sea ice variability in Arctic marginal seas have reported largely open-water conditions with significant freshwater influence during the Late glacial to Early Holocene.

Figure 5 Arctic sea-ice records (PIP₂₅ and IP₂₅ concentration) covering the Holocene: PCB (this study), ARAC04 (Wu et al., 2020a), ARA2B (Stein et al., 2017), PS51-154 and PS51-159 (Hörner et al., 2016), PS5428 (Fahl and Stein, 2012), BP0036 (Hörner et al., 2018), Ryder19-12 (Detlef et al., 2023). Average thickness of sea ice (in meters) for September 2023. Data from the Arctic Ocean Sea Ice Reanalysis (Williams et al., 2021) obtained via the E.U. Copernicus Marine Service Information MyOcean Viewer (https://marine.copernicus.eu/access-data/myocean-viewer, last access: 1 October 2025). The shaded rectangles indicate the limit of the proxy for no sea-ice in red, seasonal sea ice in blue (0.5 < PIP₂₅ < 0.75) and permanent sea ice above PIP₂₅ > 0.75.

The nearby cores JPC15 (Keigwin et al., 2018), ARAC20 (Wu et al., 2020a) (Fig. 1) recorded similar environmental changes (sea ice cover, freshwater input) as in PCB09 but different from those recorded in the shallow PCB11 site, highlighting the differences between shelf slope and outer shelf and the spatial variation of the polynya position. Aside from the close-by cores (Keigwin et al., 2018; Klotsko et al., 2019; Wu et al., 2020a), other Arctic records in the Canadian Archipelago (Vare et al., 2009), East Siberian (Dong et al., 2022), Kara (Hörner et al., 2018) , Chukchi (Stein et al., 2017), Laptev (Fahl and Stein, 2012; Hörner et al., 2016), and Lincoln (Detlef et al., 2023) Seas and along the Lomonosov Ridge (Stein and Fahl, 2012), report minimum sea-ice cover during the Early Holocene (centered around 10 ka) (Fig. 5). Detlef et al. (2023) reconstructed sea ice conditions from a sediment core covering the last 11 ka, showing that while the Lincoln Sea currently experiences perennial sea ice cover (PIP₂₅ = 0), it underwent a shift to seasonal sea ice during the Early Holocene (around 10 ka) due to significantly warmer conditions (PIP₂₅ > 0.5). This period of reduced sea ice cover is associated with increased marine productivity and meltwater input indicated by biomarker and sedimentary facies. The Northern Greenland (Detlef et al., 2023) and the site offshore the Laptev Sea on the Lomonosov Ridge (Fahl and Stein, 2012; Hörner et al., 2016) are the first regions to record permanent sea-ice cover beginning around 9 ka after the Early Holocene minimum. The Beaufort Sea (this study; Wu et al., 2020a) showed permanent sea-ice cover over the slope after 3 ka. Seasonal sea-ice cover in the shallower regions of the Laptev and Beaufort Seas (PS51-159 and PCB11) was developed after 5 and 3 ka, respectively. The Chukchi Sea (ARA2B) had seasonal sea-ice from the Mid-Holocene onwards, with an increase after 4.5 ka (Stein et al., 2017). The variations in sea ice cover and primary production in the Chukchi Sea were attributed to differences in solar insolation and variability in Pacific water inflow, which brought increased heat flux and episodic declines in sea ice cover. In the Canadian Archipelago, a record that did not include the Early Holocene (Belt et al., 2010; Vare et al., 2009), increased sea-ice cover exists from from 7 to 3 ka.

In contrast, studies using dinocyst assemblages from around the Arctic Ocean (see the review of de Vernal et al., 2013) report constant sea ice cover for the Early to Mid-Holocene with a clear decrease around 6 ka, followed by a return to pre-6 ka conditions until an increase toward modern times. This could be due to a warm-bias in the dinocyst estimate or a non-representative training set (de Vernal et al., 2013).

Together, many of the biomarker studies provide a consistent narrative of spring sea ice development during the Holocene across the Arctic Ocean following a period of high insolation in the Early Holocene. The transition from largely open-water and freshwater-influenced conditions during the Late glacial to Early Holocene to increasing sea ice cover from the Mid-Holocene onward is a shared feature across the Arctic shelf seas, although regional variations in sea ice dynamics and productivity are observed due to local freshwater input and oceanographic conditions. Evidence from areas with permanent sea ice, such as the Lincoln Sea, shows that the minimum ice cover during the Late glacial extended even into the high Arctic, offering insights into the extent of sea ice reduction during this time.

5 Conclusions

Analysis of two sediment cores from the outer Beaufort Shelf and shelf slope help elucidate the region's paleoenvironmental variability throughout the Holocene. The shelf slope experienced ice-free to minimal sea ice extent during the Late glacial and Early Holocene. During the Early Holocene, the Beaufort Shelf was ∼ 60 m shallower than today, and experienced large freshwater influxes due to the decaying LIS. The following sea level rise brought the core sites further away from the river mouth and eroding permafrost coasts, lowering the input of terrestrial organic matter. Cooling during the Mid-Holocene drove an increase in sea ice cover for the Beaufort Shelf and other Arctic marginal seas. Sea-ice cover and its impact on local upwelling and regional Pacific inflow impacted local primary production, concentrating the phytoplankton production in open-water flaw-lead or polynya conditions. Open water conditions substantially decreased during the Late Holocene as extended sea ice cover developed at the shelf slope, which caused primary productivity to further decline. This study highlights broadly consistent patterns of sea ice variability across Arctic marginal seas, implying similar controls on sea-ice dynamics, and underscoring the vulnerability of perennial sea-ice as contemporary climate conditions reach or surpass those from the Early Holocene.