Deglacial Land-Ocean Linkages at the Alaskan Continental Margin in the Bering Sea

A marine sediment record from the central Bering Sea, spanning the last 20 thousand years (ka), was studied to unravel the depositional history with regard to terrigenous sediment supply and biogenic sedimentation. Methodic approaches comprised the inference of accumulation rates of siliciclastic and biogenic components, grain-size analysis, and (clay) mineralogy, as well as paleoclimatic modelling. Changes in the depositional history provides insight into land-ocean linkages of paleoenvironmental changes. During the finale of the Last Glacial Maximum, the depositional environment was characterized by hemipelagic background sedimentation. A marked change in the terrigenous sediment provenance during the late Heinrich 1 Stadial (15.7–14.5 ka), indicated by increases in kaolinite and a high glaciofluvial influx of clay, gives evidence of the deglaciation of the Brooks Range in the hinterland of Alaska. This meltwater pulse also stimulated the postglacial onset of biological productivity. Glacial melt implies regional climate warming during a time of widespread cooling on the northern hemisphere. Our simulation experiment with a coupled climate model suggests atmospheric teleconnections to the North Atlantic, with impacts on the dynamics of the Aleutian Low system that gave rise to warmer winters and an early onset of spring during that time. The late deglacial period between 14.5 and 11.0 ka was characterized by enhanced fluvial runoff and biological productivity in the course of climate amelioration, sea-level rise, seasonal sea-ice retreat, and permafrost thaw in the hinterland. The latter processes temporarily stalled during the Younger Dryas stadial (12.9-11.7 ka) and commenced again during the Preboreal (earliest Holocene), after 11.7 ka. High river runoff might have fertilized the Bering Sea and contributed to enhanced upper ocean stratification. Since 11.0 ka, advanced transgression has shifted the coast line and fluvial influence of the Yukon River away from the study site. The opening of the Bering Strait strengthened contour currents along the continental slope, leaving behind winnowed sand-rich sediments through the early to mid-Holocene, with non-deposition occurring since about 6.0 ka.


INTRODUCTION
Climate variability during the late Quaternary was characterized by strong insolation-driven contrasts between interglacial and glacial stages, documented by the repeated built-up and decay of northern-hemispheric ice sheets (e.g., Dyke, 2004;Gowan et al., 2021). Ice sheet dynamics controlled global sea level (e.g., Lambeck et al., 2014;Miller et al., 2020) and exerted strong influence on ocean processes (e.g., Khatiwala et al., 2019;Sigman et al., 2021). The time intervals between conditions of glacial maxima and interglacial climate optima are considered as transitions between climate extremes, referred to as terminations. These terminations are of particular interest, as they are often affected by non-linear responses to external climate forcing and show complex environmental interactions with threshold feedbacks in the climate system Lohmann et al., 2020).
A closer look on such interactions is essential for oceans connecting the Arctic with the lower latitudes. As a focal area in the Meridional Overturning Circulation (MOC), the North Atlantic region since long has been identified as an important actor in the global climate system (Broecker et al., 1985;Clement and Peterson, 2008;Ziemen et al., 2019). During the last decades, increasing attention has been given to the paleooceanography of the North Pacific and its marginal seas, where interacting dynamic processes of water-mass formation, ocean stratification, sea-level fluctuations, biological productivity, seaice formation, and fresh-water pulses show global connections as integral part of the global MOC and via atmospheric teleconnections (Cook et al., 2005;Diekmann et al., 2008;Takahashi et al., 2011;Gersonde, 2012;Max et al., 2014;Pelto et al., 2018;Lohmann et al., 2019;Davis et al., 2020;Praetorius et al., 2020).
Here, we present a case study from the Bering Sea, devoted to the characterization of the marine depositional environment and terrigenous sediment supply since the LGM. Our goal was to highlight land-ocean linkages in a marginal sea, which today directly connects the North Pacific with the Arctic Ocean. We have selected sediment cores from site SO202-18 at the Alaskan continental margin (Bering Slope), recovered in the scope of the INOPEX project (Innovative North Pacific Experiment) during cruise 202 of R/V SONNE in 2009 (Gersonde, 2012) (Figure 1). The sedimentary records are suited to decipher late glacial paleooceanographic changes at high temporal resolution (Kuehn et al., 2014;Méheust et al., 2018;Meyer et al., 2019). In this contribution, we focus on glacial melt-water and fluvial runoff, in association with climate warming and sea-level rise with its potential effects on the marine environment. We especially pick up the postulation of an unusual melt-water pulse during late Heinrich Stadial 1 (HS1) in the Bering Sea after 17 ka BP (Khim et al., 2010;Meyer et al., 2019), which implies regional warming during times of northern-hemispheric cooling (e. g., Clark et al., 2012). Our approach is to use information on biological productivity from lithology and bulk sediment geochemistry as well as grain-size and mineralogical data of the non-biogenic fraction to unravel the sedimentary transportation processes and to trace sediment provenance. A modelling experiment is included to constrain climate conditions during HS1.

Regional Setting
The Bering Sea comprises a shallow continental shelf (50-150 m) in its northeastern half and deep basins down to 4,000 m water depths in its southwestern half (Figure 1). The study site SO202-18 is located at 1,105 m water depth at the southern margin of a submarine canyon on the continental slope (Bering Slope) off the wide Alaskan shelf area in the northern Bering Sea (60°07.6′N, 179°26.1′W). To the south, the Bering Sea is connected with the North Pacific through gaps between the islands of the Aleutian Arc (Stabeno et al., 2007). Bering shelf waters escape northward through the Bering Strait to the adjoining Arctic Ocean (Stabeno and Reed, 1994;Roach et al., 1995). Ocean circulation in the deep Bering Sea is dominated by cyclonic contour currents along the continental margins, comprising the Bering Slope Current, the East Kamchatka Current, and the Aleutian North Slope Current (Schumacher and Reed, 1992;Stabeno and Reed, 1994). The Bering Slope Current bifurcates close to the study site to join the outflow to the Arctic Ocean. The vertical water mass structure at the study site is characterized by less saline surface waters (32.8 psu) passing over to increased salinity at 1,000 m (34.3 psu) and on the deep floor at 3,000 m (34.7 psu); the study site lies between 500 and 1,500 m water depths, a level that today reveals suboxic conditions with oxygen less than 1.0 ml/L (Takahashi et al., 2011;Garcia et al., 2019). Sea ice covers the Bering shelf during its maximum seasonal extent and retreats north of the Bering Strait during summer (Danielson et al., 2011). The subpolar climate over the Bering Sea is driven by the variability and seasonal evolution of the Aleutian Low pressure system, controlling sea-ice coverage, cold and warm air-mass distribution, as well as snow fall and melt on surrounding land masses (Overland et al., 1999;Pickart et al., 2009). Continental runoff brings both fresh water and detrital sediments to the Bering Sea (Naidu et al., 1995;Benke and Cushing, 2011). Drainage is basically achieved by the Yukon and Kuskokwim Rivers from Alaska, and the Anadyr River from Chukotka. In addition to the permanent processes of sediment supply, episodic events cause shelf sediment transport to the slope, triggered by storms and turbidity currents. Part of the suspensions are entrained and transported by the Bering Slope Current and transferred to the northern Bering Sea (Naidu et al., 1995;Kim et al., 2015).

MATERIAL AND METHODS
Sediment cores were collected during the SONNE Expedition 202 in 2009 (Gersonde, 2012). We used the combined core sections of SO202-18-3 and SO202-18-6, which were retrieved at 1,105 m water depth from the continental slope off the wide Alaskan shelf area in the northern Bering Sea (60°07.6′N, 179°26.1′W) ( Figures  1, 2). Inter-core correlation and splicing between piston core SO202-18-3 and kasten core SO202-18-6 was based on the correlation of elemental scanning data (Kuehn et al., 2014). The age model for the composite sediment cores SO202-18-3 and SO202-18-6 was adopted from former publications of the sediment record, inferred from radiocarbon dating and tuning approaches (Kuehn et al., 2014;Meyer et al., 2019).
Non-destructive X-ray fluorescence (XRF) elemental scanning was performed using an AVAATECH XRF core scanner with a Rh X-Ray tube at 1 mA, at 10kV for major elements. Count times were set to 10 s. Peak area intensities for K and Ti were measured in counts per second (cps) and used to calculate K/Ti ratios.
Discrete samples for further sedimentological analyses were taken at 10-cm intervals, corresponding to time intervals of 50-350 years, depending on the sedimentation rates ( Figure 2). Bulk samples were used for the determination of sediment water content (w, in wet mass%), densities, organic carbon, carbonate, and biogenic silica, as well as the calculation of mass accumulation rates. About 10 cm 3 of sediment material were taken for these analyses, freeze-dried and milled. Density of the dried sediment samples (ρ d ) were measured with Micromeritics AccuPyc 1,330 gas-pycnometers. Wet bulk density (WBD), porosity (φ, vol%), grain density (GD), and dry bulk density (DBD) were calculated according to equations 1-4 that we compiled after Niessen et al. (2013) and Kuhn et al. (2017): Here, as well for all other analyses, a correction for the pore water salt within the dried samples was applied (Kuhn, 2013), where ρ s is the density of salt (2.1 g/cm3), r is the mass ratio of salt to water ( 0.036, at 3.5% salinity), and ρ w is the density of the pore water ( 1.024 g/cm3). Measurements of the total carbon (TC) and organic carbon (TOC, after CaCO 3 removal with HCl) contents were analyzed using an Elementar Vario EL III carbonnitrogen-sulphur analyzer and ELTRA CS-2000, with analytical errors of 1 and 3%, respectively. The CaCO 3 content was calculated using equation [CaCO 3 (TC-TOC) • 8.3]. Molar ratios of TOC and total nitrogen (C/N ratio) as well were calculated from these analyses. Biogenic opal (BSi•10H 2 O) was FIGURE 1 | Map of the Bering Sea with adjoining oceans and land masses and their mountain ranges as well as the location of the study site (red dot) at the Bering Slope. The map also shows bathymetry (Schlitzer, 2002), ocean currents as black arrows (Ren et al., 2014;Wang et al., 2016), seasonal sea-ice limits as colored dashed lines (Danielson et al., 2011), and main rivers as blue lines (Naidu et al., 1995). Sediment cores SO202-18-3 and SO202-18-6 were collected during expedition 202 of research vessel SONNE in 2009 (Gersonde, 2012). Abbreviations: ACC-Alaska Coastal Current; ANSC-Aleutian North Slope Current; AS-Alaska Stream; BSC-Bering Slope Current; EKC-East Kamchatka Current; AR-Anadyr River; YR-Yukon River; KR-Kuskokwim River.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 analyzed by the sequential leaching method (Muüller and Schneider, 1993). Total sediment mass-accumulation rate (MAR, in g/cm 2 /ka) was calculated from the product of sedimentation rate and the DBD. The salt-corrected component concentration in weight % was applied for specific MARs. The analysis of the non-biogenic sediment fraction, also referred to as the lithogenic fraction, was used to characterize terrigenous sediment supply, in terms of provenance and modes of transportation. Terrigenous sediment is defined as the detrital siliciclastic material carried from land to sea. Analyses included measurements of grain-size distributions, which were conducted on a 10 g split from each sediment sample. Organic materials, as well as biogenic carbonate and biogenic silica, unrelated to our primary focus on clastic transport and sedimentation, were removed using H 2 O 2 (3%), HCl (10%) and NaOH (1 M). The samples were then rinsed with deionized water and, to facilitate dispersion, treated with 100 ml of a 0.05 M (NaPO 3 ) 6 solution on a rotating hybrid instrument for 24 h. The coarse materials (>2 mm diameter) were removed by wet sieving, dried and weighed. Grain-size measurements were achieved by laser diffraction particle size analysis, using a Malvern Mastersizer 3,000, following the method described in Asikainen et al. (2006); data were thus provided on 93 subclasses between 0.375 μm and 2 mm. Because the data showed multimodal distribution patterns, an end-member statistical analysis using the MATLAB modelling algorithm of Dietze et al. (2012) was used to assess the likely relative importance of different sedimentary processes, such as fluvial, eolian, or hemipelagic transport (c.f. Wang et al., 2016). In this approach, individual grain-size populations (end-member loadings, vol%) as well as their contributions to bulk grain-size composition (scores, %) are derived from eigenspace analysis, weight transformation, varimax rotations, and various scaling procedures (Dietze et al., 2012Wang et al., 2015. Mineralogical analyses were undertaken by X-ray diffraction (XRD) on both bulk sediments and the clay fraction of nonbiogenic sediments. Bulk mineralogy was measured on freezedried and milled sediment samples, analysed with a X-ray diffraction device (XRD), using a Philips PW1820 goniometer at 40 kV, 40 mA, with Co. k-alpha radiation. XRD data were analysed using MacDiff 4.0.7 software (Petschick et al., 1996). Here we use the 3.34/(3.18 + 3.24) Å ratio as quartz/feldspar ratio (Qz/Fsp), the 3.34/10 Å ratio the quartz to illite ratio (Qz/Illite), and the half-height width at 10 Å as indicator of illite crystallinity (Illite HHW).
Semi-quantitative XRD mineralogy of the clay fraction was determined on texturally oriented aggregates of the carbonatefree clay fraction (< 2 µm), following standard methods at the Alfred Wegener Institute (Ehrmann et al., 1992;Petschick et al., 1996). Sediment subsamples for this approach (approximately 10-15 g dry) were taken along the entire length of the core. The clay and silt fractions were obtained by removing sand FIGURE 2 | Composite section of sediment cores SO202-18-3/-6, showing downcore and temporal variations in lithology, linear sedimentation rates (LSR), weight percentages (wt%) of the sediment components of carbonate, biogenic opal, total organic carbon (TOC), the molar ratio of total organic carbon versus total nitrogen (C/N ratio). Mineralogical X-ray diffraction data comprise quartz/feldspar ratios (Qz/Fsp), quartz/illite ratio (Qz/Illite) of bulk sediment, illite crystallinity (Illite HHW, and mineralogical composition (%) of the clay fraction: Smectite, Illite, Chlorite, Kaolinite).
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 (250-2000 μm, 125-250 μm, and 63-125 μm diameters) and gravel (>2 mm diameter) by wet sieving. The coarse materials were then dried and weighed. Clay (< 2 μm) was separated from silt (2-63 μm) by application of Stokes' law of settling; the clay was subsequently treated with H 2 O 2 (3%) to remove organic material and HCl (5%) to remove carbonate. Mineral proportions were calculated semi-quantitatively from weighted XRD peak areas in the glycolated state (Biscaye, 1965): Sum of expandable clays (referred to as smectite) at 17 Å, Illite 10 Å (*4), kaolinite and chlorite at 7 Å (*2). Distinction of the proportions of kaolinite and chlorite was achieved by the deconvolution of the 3.54/3.57 Å double peak. Relative analytical precisions are 6-9% and 8-14% for the major and minor clay components, respectively (Ehrmann et al., 1992). We employ the Earth system model COSMOS, which was mainly developed at the Max Planck Institute for Meteorology in Hamburg (Jungclaus et al., 2006). The model incorporates the ocean-sea ice model MPIOM (Marsland et al., 2003), the ECHAM5 atmosphere model at T31 spherical resolution (∼3.75°× 3.75°) with 19 vertical levels (Roeckner et al., 2003) and the land surface model JSBACH including vegetation dynamics (Brovkin et al., 2009). The ocean model is resolved at 40 unevenly spaced vertical layers and takes advantage of a curve-linear grid at an average resolution of 3°× 1.8°in the horizontal dimension, which increases towards the grid poles at Greenland and Antarctica (∼ 30 km). Net precipitated water over land which is not stored as snow, intercepted water or soil water is redirected towards the ocean via a high-resolution river routing scheme (Hagemann and Dümenil, 1997). Our version (COSMOS-landveg r2413) has no flux correction and has been successfully applied to test a variety of paleoclimate hypotheses, including the Pliocene (Stepanek and Lohmann, 2012), glacial (Gong et al., 2013;Zhang et al., 2013Zhang et al., , 2017Maier et al., 2018), interglacial (Wei and Lohmann, 2012;Lohmann et al., 2013;Pfeiffer and Lohmann, 2016) as well as future (Gierz et al., 2015) climates. Here, we present results obtained from model set-ups encompassing the pre-industrial (PI) (Wei and Lohmann, 2012), LGM and glacial hosing experiments (Gong et al., 2013;Zhang et al., 2013).
A lithologic description of sediment cores SO202-18-3/-6 is provided in Figure 2. The sediment section is mainly composed of high proportions of lithogenic mud (>70%) with lower amounts of biosiliceous opal (diatoms, radiolarians). Biogenic carbonate (foraminifera/coccoliths) ranges below 5%. The Holocene section is enriched in sand, diatomaceous sand dominated the weakly laminated sediment until 10.2 ka. Maxima of opal appear in the late glacial interval between 16.5 and 10.2 ka, which were deposited at high sedimentation rates (up to 250 cm/ka) in accordance with maxima in organic carbon. A prominent feature is the occurrence of several, meter-thick laminated sequences in the intervals of the B/A interstadial and the PB stage. The laminae are composed of mm-to submm thick couplets of pure light-colored diatomaceous ooze, alternating with dark diatom-bearing lithogenic clayey silt, frequently including remains and fibers of marine macro algae. From 10.2 ka, two thicker, well laminated diatomaceous ooze sequences followed downward until 14.6 ka, interrupted and followed by diatomaceous mud between 11.7-12.9 ka and 14.6-16.5 ka. Lithogenic mud was visible down to the bottom of the core. Detailed core descriptions are provided by Bernhard Diekmann and co-workers in Gersonde (2012) as well as by Kuehn et al. (2014).

Grain-Size Distribution
The SO202-18 sediment record exhibits marked variations in the grain-size distribution of the non-biogenic fraction ( Figure 3). The mean grain size ranges between 4 and 53 μm ( Figure 3A), indicating the dominance of silt fractions throughout the core. The downcore wet sieving and setting weight percentages (>63 μm sand, 2-63 μm silt, and < 2 mm clay) ( Figure 3A) show dominant and gentle silt from 36 to 74 wt%. In addition to a dominance of silt, sand constitutes up to 20 wt% and increases to 30-60 wt% in the Holocene part ( Figure 3A). Clay content generally varies moderately between 4-20 wt% through the core and increases up to 31 wt% at 15 ka. The K/Ti ratio of the XRF elemental data strongly correlates with clay concentrations, reflecting the dominance of K-rich clay minerals in the very fine-grained sediments ( Figure 3C). Statistical unmixing of grain-size distributions yielded an optimal model with three end members (EMs), explaining 97% of the data variance, of which EMs one to three explain 82, 12, and 3%, respectively ( Figure 3B). The grain-size end members EM1 to EM3 reveal peak modes at 3 µm (very fine silt), 83 µm (very fine sand) and 43 µm (very coarse silt), respectively. EM1 and EM3 are respectively dominant at the lower and upper core sections and their scores vary from 0 to 99%, while EM2 accounts for more than 50% throughout the rest sections. The down-core distribution of end-member scores ( Figure 3A) confirm the dominance of sand in the Holocene part, of silt throughout the core, and a pronounced clayey interval between 16.5 and 11.5 ka BP. Secondary subdued peaks of EM1 and EM2 coincide with the major peaks of each other, may represent artifacts of statistical noise or two-phase sedimentary dynamics.

Bulk and Clay Mineralogy
Stratigraphic variations in core SO292-18-3/-6 clay minerals and the related intensity ratios of Quartz/Feldspar (Qz/Fsp) ratio, Quartz/Illite ratio and Illite HHW are shown in Figure 2. The mineralogy of the lithogenic sediment fraction, as inferred from XRD analysis, is quite uniform throughout the studied sediment Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 section. Smectite, Illite, chlorite and kaolinite were the major clay minerals. Smectite overall varied within a rather low range (3-17%), only intermittently does smectite account for more than 20% (e.g., at 8.6 ka BP). Illite content ranged from 40 to 53%, showed a generally featureless distribution throughout the core. Chlorite, a detrital product of high-latitude weathering conditions, contributes only occasionally during Preboreal (>40%). Kaolinite content ranged from 4 to 19%, being Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 6 generally somewhat constant (approximately 7%) but with relatively high values (>10%) during 16.5 and 14.5 ka BP, the upper part of HS1. In bulk mineralogy, an increase in quartz/ feldspar ratios (>5%) characterize the sand-rich sediments of the Holocene. Poor illite crystallinity marks the upper HS1 interval with the >0.5 values.

Modeling Results
The results of our hosing experiment are shown in Figure 6. To explore the HS1-to-LGM changes of surface air temperatures and surface winds associated with the Aleutian Low-pressure system, we have re-analyzed results obtained from model set-ups encompassing the pre-industrial (PI) (Wei and Lohmann, 2012), LGM and glacial hosing experiments (Gong et al., 2013;Zhang et al., 2013). Here, details of the simulated PI climate have been well documented and compared to climate data in Wei and Lohmann (2012), while the COSMOS-modelled LGM and HS1 conditions have been shown comparable with the simulations by most PMIP3 models (c.f. Gong et al., 2013;Zhang et al., 2013).

DISCUSSION
The lithostratigraphy of the sediment record from site SO202-18 confirms descriptions, datings, and compositional data from other sediment records along the Bering Slope (Sancetta et al., 1984;Cook et al., 2005;Pelto et al., 2018). In the following section, we will consider modes of sedimentary processes and terrigenous sediment provenance at site SO202-18 on the Alaskan continental margin of the Bering Sea. Changes of both will be discussed in the context of marine processes as well as regional and climatic changes.

Modes of Terrigenous Sediment Supply
Modes of sedimentary processes in the marine environments of the Bering Sea and nearby Arctic Ocean can be inferred in their relative intensity from grain-size distributions and their subpopulations that reflect the intensity, transport capability, and consistency of glacial, sea-ice, aquatic, and eolian sediment sources and transport (Darby et al., 2009;Wang et al., 2016;Deschamps et al., 2018;Murdmaa et al., 2019).
Of the three grain-size end members, identified in the SO202-18 record (Figure 3), EM3 dominates and has broad bimodal distribution in the coarse silt fraction between 20 and 50 µm with major peak at 44 um and minor peak at 20 µm. In recent surface sediments of the North Pacific and the Bering Sea, this end member is widely distributed in the Bering Sea (Wang et al., 2016). Silt particles >10 µm represent the so-called "sortable silt", which is transported within nepheloid layers of focused ocean currents (McCave et al., 1995;Diekmann et al., 2008;Tegzes et al., 2015;Wang et al., 2016. In the temporal context, we see EM3 as an indication of hemipelagic background sedimentation, which prevailed during the pre-Holocene except HS1 period at the study site. In other sediment cores from the Bering Slope, taken during IODP expedition 323, silt-rich sediments, consistent with our grainsize EM3, have been found to characterize glacial periods (Lund et al., 2021). Their interpretation of silt-rich sediments as drift deposits underneath prevailing sea ice during glacial stages is consistent with our explanation of dominant background sedimentation (Lund et al., 2021). Also, on the Shirshov Ridge in the western Bering Sea, silt-rich sediments, similar to our EM3, are interpreted as hemipelagic sediments under the influence of bottom currents (Murdmaa et al., 2019). According to the strong variation through the deglacial process (e.g., shelf flooding, and intensification of the Bering Slope Current), EM3 in our study does not follow a clear glacial-interglacial variation pattern. As hemipelagic background sedimentation, explaining 3% of the data variance, EM3 showed low values or even absence during the peaks of other end members.
The grain-size population of EM2 peaks in the fine sand fraction at 80 µm and is the coarsest grain size observed in our record. It also characterizes the surface sediments at the study site (Wang et al., 2016). The common interpretation is that sand-rich sediments, which generally characterize interglacial periods at the Bering Slope, originate in coarse materials that become mobilized during interglacial times of flooded shelves (Lund et al., 2021). In the pelagic Bering Sea, sand-rich sediments point to sea-ice rafting of anchor ice from the shelves (Murdmaa et al., 2019). We favour our interpretation that the enrichment in sand possibly points to strong winnowing effects, causing residual sand enrichment by the vigor contour currents of the Bering Slope Current that branches off towards the Bering Strait at the study site (Figure 1). Moreover, the rough local topography of the continental slope might promote strong turbulence, enhancing winnowing effects. Strong current activity might also explain the absence of sediments younger than 6 ka at the site. In nearby core ARC6-B10 from the Bering Slope of deeper water depth (2,493 m), grain-size coarsening characterize the time since the mid-Holocene, when our sediments disappear (Sun et al., 2021). This is consistent with the invigoration of currents, winnowing in our shallow site and coarsening depositions in the deeper site. If other sedimentary processes are responsible, such as those by seaice transportation, we should see an uniform grain-size signal at both sites.
EM1 has a broad mode from clay to fine silt (1-4 µm) and represents the finest observed grain-size population, appearing in the upper part of HS1 (16.2-14.7 ka). The down-core variation of EM1 is also clearly displayed by the curves of mean grain size and K/Ti ratios (Figures 3, 4). As for the other end members, according to information from literature, several ways of interpretation of EM1 are possible. Such clay-rich sediments were also recognized in other sediment records from the Bering Slope for the same time (Pelto et al., 2018). In modern sediments, this grain-size population cannot be found in the Bering Sea (Wang et al., 2016). Today, fine silt and clay delineates the trajectories of the westerly winds and the spread of fine dust from Central Asia to the pelagic North Pacific (Muhs, 2013;Wang et al., 2016). In modern depositional environments at continental margins, clays usually flocculate at the fresh-water/ ocean-water interface, in the distal offshore part of deltas and river mouths (Coleman and Wright, 1975;Diekmann et al., 2008). Otherwise, peaks in clay at the Bering Slope can be seen as an indication of coastal erosion and inundation as well as particle transport by dust storms in early deglacial times (Pelto Frontiers (Lambeck et al., 2014); (e, f, g) loading of grain-size end member 1, the K/Ti ratio, and clay proportions as proxies of fluvial runoff; (h, i) illite crystallinity (HHW) and kaolinite proportion as indicators of provenance change during late HS1; (j, k, l) mass-accumulation rates (MAR, g/cm 2 / ka) of carbonate (Carbo), terrigenous sediment supply (litho), and biogenic opal. The blue shade marks the timing of the glaciofluvial melt-water pulse during late HS1. Light brownish shades indicate periods of high biological productivity with dark brown shades in conjunction with laminite preservation and anoxic events (Kuehn et al., 2014).
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 8 et al., 2018). In polar and subpolar regions, another factor is fine sediment trapped in sea ice and released by melting, an important process in the Arctic Ocean (Nürnberg and Tiedemann, 2004;Darby et al., 2009;Wang et al., 2021). In the pelagic Bering Sea, clay-rich sediments are related to hemipelagic sedimentation under sea ice and weak bottom currents (Murdmaa et al., 2019). From the different interpretations, we conclude that the finest observed grain-size population, represented by EM1, might be best related to the influence of fluvial suspension load in a near-coastal setting. However, this interpretation can only be validated, when considering a marked change in sediment provenance and the dynamics of deglaciation in Alaska, as outlined in the next sections.
In summary, we see clay-rich sediments of EM1 at late HS1 as products of nearby fluvial runoff of aquatic suspensions. Sea-ice cover and sea level rise suspension may support undisturbed settling of clays. Silt-rich sediments of EM3 represent hemipelagic background sedimentation under the influence of contour currents throughout the studied section. Sand-rich sediments of EM2 originate from remobilization of shelf sediments and winnowing of fines under turbulent bottom currents during the Holocene. The recognition of sediment provenance underpins these interpretations.

Terrigenous Sediment Sources
With the exception of the clay-rich interval during late HS1, mineralogical provenance data at site SO202-18 show a consistent downcore pattern (Figure 2), pointing to quite persistent sources of terrigenous sediments through time. They basically originate from the surrounding land masses of Siberia, Alaska, and the Aleutian Islands and are dispersed by the counterclockwise ocean currents of the Bering Sea (Naidu and Mowatt, 1983;VanLaningham et al., 2009;Wang et al., 2016). Furthermore, the assemblage varies within the limits of the mineralogical composition of sediment loads in modern rivers around the Bering Sea (Naidu and Mowatt, 1983). Such a pattern of low variability and common sediment sources through time has also been inferred from the pelagic Bering Sea for the late Quaternary (Murdmaa et al., 2019) and for the last 2.4 million years from IODP Site U1343, situated in the central part of the Bering Slope (Kim et al., 2015). However, the latter record has no sufficient temporal resolution to be comparable with our record. A former study over the last 60 ka from a nearby sediment core, nonetheless, exhibits a marked kaolinite spike of >15% around 16 ka, though only represented by one sample (Kim et al., 2011). In our record, however, a clear change in mineralogy between 15.7 and 14.5 ka BP, is well displayed by increased kaolinite concentrations and the presence of poorly crystallized illite at high sampling resolution (Figures 2, 4). The clay mineralogy of another nearby sediment core has been studied over the time interval from 13.7 ka to present (Sun et al., 2021). It shows kaolinite concentrations around 10% in the pre-Holocene part, declining to lower values between 3 and 7% in the Holocene part, but never reaching such high values around 20%, as in the HS1 interval of SO202-18.
Today, such an assemblage with increased kaolinite proportions cannot be found in any of the modern surface sediment samples of the Bering Sea (Wang et al., 2016). Clues for a regional kaolinite source arise from patches of clays with high kaolinite/chlorite ratios in modern sediments adjacent to the Alaskan coast, extending farther on the shelf off the modern Yukon Delta (Naidu and Mowatt, 1983). In neighboring seas, kaolinite concentrations in marine sediments are generally low in the North Pacific (Wang et al., 2016) and in Siberian shelf surface sediments of the Arctic Ocean (Viscosi-Shirley et al., 2003). In any case, higher kaolinite concentrations appear in marine sediments north of the Bering Strait in the Chuchki and Beaufort Seas (Naidu et al., 1971;Wang et al., 2021) and are also evident as detrital components in Tundra soils along the northern Alaskan coast of the Arctic Ocean (Borden et al., 2010). A common source for Yukonderived sediments and marine sediments along the Arctic coast of Alaska are kaolinite-bearing Mesozoic sedimentary rocks and paleosols as well as altered Paleozoic tuff beds, that prevail in northern Alaska and form part of the mountain chain of the Brooks Range (Naidu et al., 1971;Naidu and Mowatt, 1983;Dalrymple and Maass, 1987;Borden et al., 2010;Darby et al., 2011). Kaolinite in the course of the Yukon River may also originate from diagenetic clays in Cretaceous-Cenozoic coal deposits of the Yukon-Tanana River catchments of central Alaska (Triplehorn, 1976;Merritt, 1986;Levitan et al., 2013).
The kaolinite maximum in our record indicates that the sediment source in northern and central Alaska was temporarily relevant for the depositional environment at site SO202-18 after the LGM, when the Bering land bridge blocked the connection to the Arctic Ocean and the Alaskan coastline was close (Bond, 2019). The Bering Slope either became directly connected to fluvial runoff of the Yukon River or older sediments from the catchment became remobilized. The following interpretation of the regional paleoenvironmental history from our sedimentary proxy record gives support to test this hypothesis.

Land-Ocean Linkages of Marine Depositional Processes Through Time
The interpretation of the marine depositional environment at the Bering Slope has to consider climate-driven changes in sea-level (transgression), glacial dynamics in the continental catchment as well as oceanographic processes.

Time Interval 20-16.5 ka BP: Background Sedimentation at Low Sea Level
This period around the end of the LGM, the ED and early HS1 was dominated by hemipelagic background sedimentation of siliciclastic muds with subordinate biogenic compounds (Figures 2, 4, 5). The Laurentide and Cordilleran Ice Sheets reached their maximum extent (Dyke, 2004), including the Brooks Range, Alaska Range and southern coast of Alaska (Briner and Kaufman, 2008). In contrast, glaciation of eastern Siberia was restricted to isolated mountain patches (Barr and Clark, 2012). Global eustatic sea-level was up to 135 m lower than in the Holocene (Lambeck et al., 2014). Maximum sea ice in spring covered the whole Bering Sea and subpolar North Pacific (Méheust et al., 2018). Though the emerged Alaskan shelf caused site SO202-18 to be located closer to the coast, fluvial discharge under such ice-age conditions was diminished, documented by Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 low clay input. Frozen ground in the catchment trapped sediments in the permafrost (Meyer et al., 2019). The marine sediments with high concentrations of sortable silt (EM3) point to an active slope current at site SO202-18.

Melt-Water Discharge
The beginning of this time interval in the younger part of HS1 is marked by an abrupt shift in sediment provenance (indicated by kaolinite, poorly crystallized illite). During this time, freshwater supply with increased suspension load brought clays to the Bering Slope. Mass-accumulation rates of both biogenic remains and terrigenous sediments started to rise (Figure 4). Terrestrial processes influenced the marine depositional environment. The runoff event started after the initiation of postglacial climate amelioration and associated northern-hemispheric glacial melt and sea-level rise since about 19 ka BP (Carlson and Winsor, 2012;Gowan et al., 2021). During HS1, sea-level rise reached a level of FIGURE 5 | Paleoenvironmental scenarios in the Bering Sea area during deglaciation, as inferred from this study. Map drawn with Ocean Data view (Schlitzer, 2002). Glacial extent and chronology in Alaska according to Dyke (2004) and Briner and Kaufman (2008). The timing of deglaciation on the Siberian side is still unclear and shows the maximum extent during the LGM ( Barr and Clark, 2012). The configuration of river systems follows reconstructions by VanLaningham et al. (2009) andBond (2019). Modern surface currents shown modified after Kim et al. (2015).
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 10 about 80 m lower than today, still leaving a well exposed shelf area close to the study site ( Figure 5).
In the Bering Sea catchment, this time was characterized by general ice-sheet shrinkage in Alaska and saw the deglaciation of the Brooks Range (Dyke, 2004;Briner and Kaufman, 2008;Pendleton et al., 2015). The glacial meltwaters produced the marked signal in provenance change through the tapping of the Yukon River catchment to the receding glaciers, which before were mostly frozen and untouched by rivers. A big proglacial lake in the Old Crow Basin (roughly 60 × 80 km in diameter), situated at the southeastern foot of the Brooks Range, was drained to the upper tributaries of the Yukon River and enhanced fluvial discharge to the Bering Sea (Schweger et al., 1989;Zazula et al., 2004). Meltwater from the Brooks Range and incipient erosion and remobilization of exposed bed rocks, sediments, coal deposits and soils in the valleys of central Alaska increased clay fluxes from distal sources to the Bering Sea. This process also explains the enrichment of terrigenous organic matter (Khim et al., 2010;Meyer et al., 2019) and petrogenic carbon (Meyer et al., 2019) in marine sediments of the Bering Slope, as also indicated by elevated C/N ratios ( Figure 2). The fluvial supply of nutrients apparently fertilized the sea water and stimulated biological productivity, as suggested by increased accumulation rates of biogenic silica and carbonate together with clay (Figures 4, 5). Possibly the processes were enhanced by pronounced seasonal runoff strongly during spring and weakened during summer, permitting biological productivity without light-limiting sea-ice coverage. Clay settling may have appeared during winter time underneath the protecting cover of sea ice. However, direct clues for such seasonal processes are first preserved in the overlying sediments of the B/A interstadial.
We are aware that the interpretation of fresh-water discharge should be documented by stable oxygen isotope proxies of salinity change, as for instance documented in the North Pacific (Maier et al., 2018). Such direct clues so far were not recognized in the Bering Sea (e.g., Cook et al., 2005;Caissie et al., 2010;Riethdorf et al., 2013) and not investigated in our record. However, changes in regional oceanic conditions are evident at the central Bering Slope, where biomarker records indicate a change towards marginal sea-ice conditions during late HS1 after 16 ka BP (Detlef et al., 2020). At the southern Bering Slope a strong pulse in the appearence of the diatom species Odontella aurita marks the time interval between 16 and 15 ka BP, indicating the onset of deglacial open-water diatom blooms (Caissie et al., 2010). Moreover, in the Bering Sea, O. aurita often is associated with salinity fronts (Sancetta et al., 1984) and usually stands for a high-productive, neritic environment in a marginal ice-zone setting.

Time Interval 14.5-12.9 ka BP: First Biological Bloom Event
This time is consistent with late glacial climate amelioration during the B/A interstadial Lohmann et al., 2020, Max et al., 2014, Kuehn et al., 2014. After the global and regional meltwater pulses, sea level reached −60 m below present (Lambeck et al., 2014) (Figures 4, 5). Fluvial runoff from the nearby coast remained active (Figure 5), as indicated by still high values of the clay proxies. Fluvial discharge was driven by the modern hydrological configuration of a catchment under warmer climate conditions, but did no longer include glacial melt waters from central Alaska. The latter is indicated by a return in sediment provenance towards the modern Bering Sea mineral assemblage with low amounts of kaolinite (Wang et al., 2016). Support for this interpretation also comes from the observation of cessation of petrogenic carbon supply (Meyer et al., 2019). Terrigenous young carbon, now was released by permafrost thaw and the formation of tundra wetlands in the catchment (Meyer et al., 2019). The interstadial moreover was characterized by enhanced biological export production, shown by an increase in the biogenic sediment components: diatoms, calcareous plankton, and marine organic carbon at high accumulation rates (Kuehn et al., 2014;Pelto et al., 2018;Gorbarenko et al., 2019) (Figures 4, 5). Substantial parts of the time interval were characterized by anoxic conditions that led to the preservation of laminites (Cook et al., 2005;Caissie et al., 2010;Kuehn et al., 2014). This time interval covers the YD stadial (12.9-11.7 ka BP), which went along with a reduction in fluvial runoff (EM1 and clay contents) and biological productivity (Figures 4, 5), during a period of northern-hemispheric cooling (Broecker et al., 1985;Clark et al., 2012;Lohmann et al., 2020). The transition of the YD to the PB was accompanied by a dramatic change in the depositional environment. Between 11.5 and 10 ka, postglacial climate warming was accompanied by the second biological bloom period including repeated periods of anoxic conditions (Figures 4, 5). However, fluvial runoff no longer reached the study site, because the global Meltwater Pulse 1B shifted the sea level to 40 m below present (Lambeck et al., 2014) and caused a dramatic lateral shift of the coast line ( Figure 5). This shift possibly also initiated the opening of the Bering Strait (Figures 4, 5) and established the modern vigorous current system at the study site. This switch is well documented by residual sandy sediments, left behind by winnowing processes. The assumption of a Bering Strait opening fits latest estimates from geological evidence around 11 ka (Keigwin et al., 2006;Jakobsson et al., 2017), but contradicts simulation experiments related to glacial melt and sea-level rise that suggest a timing prior to the YD stadial (Pico et al., 2020).

Paleoenvironmental Implications and Outlook
The outcomes of our study cast a light and leave open questions on two specific aspects of the deglacial environment in the Bering Sea: 1) incipient climate amelioration already during late HS1, and 2) a prominent role of fresh-water fluxes to the Bering Sea with possible effects on oceanographic processes and biological productivity.
The recognition of ice-sheet retreat in Alaska (especially in the Brooks Range) calls for early postglacial climate warming that allowed for glacial melt. Warming in Alaska during HS1 thus appeared during a time, when temperatures dropped in the North Atlantic region (Clement and Peterson, 2008;Praetorius and Mix, 2014;Zhang et al., 2017;Gong et al., 2019;Lohmann et al., 2019;Ziemen et al., 2019). Warming in Alaska and NW Canada is confirmed by paleolimnological proxy records, showing continuous summer warming through the ED, HS1, to B/A (Kurek et al., 2009;Fritz et al., 2012). In the Bering Sea and North Pacific, changes in sea-Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 surface temperatures in detail revealed spatiotemporal heterogeneities, but with a trend towards warming during HS1 in the northwestern Pacific and the Bering Sea (Meyer et al., 2016;Davis et al., 2020;Praetorius et al., 2020). The evidence of late HS1 warming points to both insolation and greenhouse gas forcing of regional climate ( Figure 4). Model experiments were applied to check the rationale for regional climate warming ( Figure 6). Glacial hosing experiments into the northern North Atlantic mimic HS1 conditions (Gong et al., 2013). With the exception of Alaska and the northern Bering Sea, the HS1 temperature anomalies with respect to the LGM give evidence of widespread cooling through all seasons in the North Pacific region and adjoining continents, with the exception of Alaska and the northern Bering Sea, where winter and spring temperatures increased by about 0.8°C. Such a regionalized warming can be regarded as a response to the cooling and increased sea-level pressure over the North Atlantic, which strengthened the Aleutian Low over the subarctic Pacific area. Once a stronger AL, surface wind anomalies in anti-clockwise direction will prevail over the subarctic Pacific region in glacial winter. More specifically, this gives anomalous winds from the Siberian region to the Okhotsk Sea and the western Bering Sea, thus a stronger transport of landsourced, relatively drier and colder air mass southwards. In the contrast, over the Bay of Alaska, the stronger AL exhibits itself in form of northward wind anomalies, thus enhancing a delivery of warmer and more-moisture air mass from lower latitudes ocean. Overall, change in the AL strength acts to shift the climate over the western and northern subarctic Pacific in the opposite directions (c.f. Gong et al., 2019). In consequence, less snow fall during winter and an earlier start of glacial melt in spring is detected in Alaska. The balance between ice accumulation and ablation thus was shifted towards the loss of ice.
Regional warming and incipient deglaciation enhanced continental runoff at low sea level to the submarine canyon at the Bering Slope and affected the depositional environment at our study site. Nutrient supply from land-derived aquatic suspensions might have fertilized the Bering Sea during late HS1 as well as during the B/A and PB interstadials as important prerequisite for marine biological productivity. With the exhaustion of glacial melt water supply after HS1 (indicated by our sedimentological data), fluvial terrigenous sediment supply was fed by the inundation of the Bering Shelf and permafrost thaw in the catchment, causing loosening and erosion of formerly frozen ground (Meyer et al., 2019).
Land-ocean linkages during times of enhanced fresh-water supply are well known from the North Atlantic during MIS 2-3, but also have been recognized for the wider North Pacific region (Maier et al., 2018;Praetorius et al., 2020). A prominent case is related to the deglaciation of the Cordilleran ice sheet in southern Alaska, where the supply of glacial debris by icebergs and meltwater likely spurred biological productivity in the Gulf of Alaska during HS1 (Maier et al., Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 712415 2018), and particularly with the onset of B/A warming (Davies et al., 2011;Müller et al., 2018). As at our study site, another decisive factor there was the remobilization of sediments by transgression of the shelf (Davies et al., 2011). The clarification of the prominent role of fresh-water flux on biological productivity at the Bering Slope, suggested by our study, has to invoke other important processes that cannot directly be constrained from our sedimentological proxy data. One important factor is the dynamics of sea-ice coverage with its effects on light availability, which so far has not been fully reconstructed over the deglacial time interval in our record (Méheust et al., 2018). Episodes of suboxic to anoxic bottomwater conditions at site SO202-18 between the B/A and PB contributed to a better preservation of biogenic and organic matter (Kuehn et al., 2014), restricting the interpretation of the intensity of primary biological productivity. The availability of oxygen in the North Pacific and its marginal sea depends on ocean stratification and the configuration of North Pacific Intermediate Water and North Pacific Deep Water and vertical convection controlled by sea-ice dynamics . Ocean stratification on the other hand is also related to fresh-water supply (Lam et al., 2013;Praetorius et al., 2020). Further proxy studies and quantitative modelling approaches are needed to disentangle all these interrelated processes.

CONCLUSION
Our reconstruction of the late Quaternary depositional and paleoenvironmental history at the Bering Slope points to landocean linkages, strongly related to fresh-water influx and postglacial transgression with impacts on marine biological productivity and oceanography since 20 ka: • The period 20-16.5 ka around the end of the LGM to early HS1 was dominated by hemipelagic background sedimentation under the influence of contour currents, dominated by siliciclastic muds with subordinate biogenic compounds. • During late HS1, between 16.5 and 14.5 ka BP, hemipelagic sedimentation was supplemented by strong runoff from distal sediment sources, mainly by the Yukon River, indicated by a change in terrigenous sediment provenance and a high proportion of clays from aquatic suspensions. Fluvial sediment loads were fed by meltwaters from the Brooks Range, which lost its ice sheet during that time. Glaciofluvial runoff stimulated the onset of postglacial biological productivity. • The time interval between 14.5-12.9 ka BP covers the B/A interstadial and was also dominated by fluvial runoff of clay suspensions, now from near-coastal sources, in the course of late glacial climate amelioration and sediment reworking by sea-level rise and thawing permafrost. The marine depositional environment was dominated by high biological productivity and the temporary preservation of laminites under oxygenpoor conditions.
• A cooling episode between 12.9 and 11.5 ka BP during the YD stadial, slowed down biological productivity, which was rejuvenated during the PB, the earliest stage of the Holocene. The opening of the Bering Strait around 11 ka strengthened contour currents and turbulence at the Bering Slope, leaving behind residual sand-rich sediments, heralding the modern depositional environment at the study site. • Wider implications of the study are related to the recognition of ice-sheet retreat in Alaska (especially in the Brooks Range). This is consistent with geomorphological findings and land-based proxy records of regional warming, which thus appeared during a time, when temperatures dropped in the North Atlantic region during HS1. Our modelling results indicate possible teleconnections to the North Atlantic with effects on the dynamics of the Aleutian Low, that gave way to regional warming during winter and spring over Alaska and parts of the Bering Sea.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
RW: Main author of the paper, which is part of his doctoral thesis. Responsible for interpretations, text, and figures. Raised the grain-size data and mineralogical data. GK: Provided data on sediment composition and accumulation rates and contributed to the text and the paleoenvironmental interpretations. XG, GL: Did the modelling, provided text and paleoclimatic interpretations. 41706213). GL, XG, and BD received funding through PalMod. We thank the scientific and technical staff on shipboard. We are also grateful to laboratory and computer staff at Alfred Wegener Institute.