Regional beryllium-10 production rate for the mid-elevation mountainous regions in central Europe, deduced from a multi-method study of moraines and lake sediments in the Black Forest

Beryllium-10 cosmic-ray exposure (CRE) dating has revolutionized our understanding of glacier fluctuations around the globe. A key prerequisite for the successful application of this dating method is the determination of regional production rates of in situ accumulated ¹⁰Be, usually inferred at independently dated calibration sites. Until now, no calibration site has been available for the mid-elevation mountain ranges of central Europe. We fill this gap by determining in situ ¹⁰Be concentrations in large boulders on moraines and by applying radiocarbon and infrared-stimulated luminescence (IRSL) dating to stratigraphically younger lake sediments in the southern Black Forest, SW Germany. The dating methods yielded concordant results, and, based on age–depth modelling with ¹⁴C ages, the age of a cryptotephra, and IRSL ages, we deduced a regional ¹⁰Be production rate in quartz. Calibrating the Black Forest production rate (BFPR) in the Cosmic-Ray Exposure program (CREp) resulted in a spallogenic sea-level and high-latitude (SLHL) production rate of 3.64±0.11 atoms ¹⁰Be g⁻¹ quartz a⁻¹ when referring to time-dependent Lal–Stone scaling, the European Reanalysis (ERA)-40 atmosphere model, and the atmospheric ¹⁰Be-based geomagnetic database in CREp. The BFPR turned out to be ∼11 % lower than both those at the nearest calibration site in the Alps (4.10±0.10 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL) and the canonical global ¹⁰Be production rate (4.11±0.19 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL) in CREp. A stronger weathering and snow cover bias and a higher impact of forest, soil, moss, and shrub cover at the study site likely explain this discrepancy.

1 Introduction

Beryllium-10 cosmic-ray exposure (CRE) dating is an invaluable tool for age determination in the field of glacial (e.g., Hofmann et al., 2024a), fluvial (e.g., Schoch-Baumann et al., 2022), coastal (e.g., Dawson et al., 2022), and periglacial (e.g., Amschwand et al., 2021) geomorphology, as well as for the dating of mass movements (e.g., Hilger et al., 2018). Calculating CRE ages requires the determination of effective ¹⁰Be production rates at sampling sites. CRE age calculators, such as the Cosmic-Ray Exposure program (CREp; Martin et al., 2017), utilize physical models, such as the Lifton–Sato–Dunai (LSD) scaling scheme (Lifton et al., 2014), to extrapolate ¹⁰Be production rates at calibration sites to sampling sites. The production rate is a key parameter and critically determines the robustness of the ages. At independently dated calibration sites (e.g., Goehring et al., 2012; Claude et al., 2014; Stroeven et al., 2015), geological ¹⁰Be production-rate calibration allows the determination of the rate of the in situ accumulation of cosmogenic ¹⁰Be in rock surfaces with time (cf. Dunai, 2010). Production-rate calibration often relies on radiocarbon ages of organic material that was deposited soon after the exposure of the sampled rock surfaces (e.g., Goehring et al., 2012). However, previous studies (e.g., Small and Fabel, 2016a, b; Lowe et al., 2019) have questioned the robustness of radiocarbon ages that have previously been used for geological calibration. Thus, applying a robust approach, ideally involving multiple independent dating methods, to a potential calibration site is the best practice to overcome dating-method-specific issues. If several lines of evidence converge, the resulting production rate might be more reliable than a production rate based on a single geochronological method.

Particularly thanks to joint efforts, such as the CRONUS (Cosmic-Ray-prOduced NUclide Systematics)-Earth (Phillips et al., 2016) and CRONUS-EU projects (Stuart and Dunai, 2009), the number of ¹⁰Be calibration sites has steadily increased over the last decades (see Martin et al., 2017, for a map of calibration sites). In CREp (Martin et al., 2017), Europe is one of the regions with the most published regional ¹⁰Be production rates, ranging from 3.48±0.10 (Fenton et al., 2011) to 4.46±0.26 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at sea level and high latitude (SLHL; Borchers et al., 2016) when referring to time-dependent Lal–Stone scaling (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008), the European Reanalysis (ERA)-40 atmosphere model (Uppala et al., 2005), and the atmospheric-¹⁰Be-based geomagnetic database of Muscheler et al. (2005). According to Martin et al. (2017), the majority of the European calibration sites are situated at elevations below 500 m above sea level (m a.s.l.). The only exceptions are Maol Chean-dearg in Scotland (Borchers et al., 2016) and the Chironico landslide in southern Switzerland (Claude et al., 2014), located at elevations of 521 and 761 m a.s.l., respectively. However, for rock surfaces at higher elevation (e.g., Le Roy et al., 2017), snow shielding during winter affects effective production rates more strongly than at lower elevation (cf. Ivy-Ochs et al., 2007). Rates of post-depositional weathering or removal of rock (“erosion” or “rock decay” sensu Hall et al., 2012) may also differ from those at production-rate calibration sites. Thus, extrapolating ¹⁰Be production rates from lower-elevation calibration sites to sampling locations at higher elevation may induce an unwanted bias in CRE ages. In addition, soil and vegetation cover (trees, shrubs, and moss) may influence effective production rates. To solve these issues, some studies (e.g., Boxleitner et al., 2019) corrected their ages for these factors, whereas others have presented uncorrected ages and have interpreted them as minimum estimates (e.g., Protin et al., 2019). Other authors performed sensitivity tests to assess the effect of these factors on ages (e.g., Hofmann et al., 2022, 2024b). However, obtaining independent age control is the best strategy to validate proposed correction factors.

The formerly glaciated mid-elevation mountain ranges of central Europe comprise the Jura and the Variscan mountains, i.e., the Vosges, the Black Forest, the Harz, the Bavarian/Bohemian Forest, and the Giant Mountains (Ehlers et al., 2011). During the Late Pleistocene, the ice masses in this region lay in the mostly non-glaciated corridor between the ice sheets over northern Europe and the glacier network of the Alps. As previously discussed (e.g., Hofmann et al., 2020, 2022, 2024b; Hofmann, 2023a), there is an urgent need to date the onset of retreat from their Late Pleistocene maximum positions. Clarifying this issue would help to evaluate the hypothesis that the Alps shielded the ice caps and glaciers from humid air masses from the Mediterranean Sea during the last major advance of piedmont lobe glaciers in the forelands of the Alps (at around 25 ka; e.g., Gaar et al., 2019).

Starting with a pioneering study in the Vosges (Mercier et al., 1999), France, CRE dating became the key method for age determination of ice-marginal positions (e.g., Mercier et al., 2000; Reuther, 2007; Mentlík et al., 2013; see Hofmann et al., 2022, for a compilation), allowing us to overcome three major limitations. Firstly, the lack of organic material prevented radiocarbon dating of moraines. Secondly, luminescence and radiocarbon dating of deposits present in cirque basins only gave minimum deglaciation ages (e.g., Vočadlová et al., 2015). Thirdly, although investigating the sedimentary sequence in the overdeepened Elbe Valley in the Giant Mountains (Czechia) undoubtedly provided valuable insights into the environmental history, applying luminescence and radiocarbon dating to sediments below and above the till of the last glaciation only resulted in a bracketing age for the last glaciation (see Engel et al., 2010). Due to the lack of a regional calibration site, the abovementioned studies have relied on calibration sites in other regions, such as the Chironico landslide in the southern part of the Alps (Claude et al., 2014). It should be noted that each calibration site has its own bias with regard to the amount and timing of snow cover, the type of vegetation, and the rates of weathering. However, the topoclimatic and environmental conditions of the entire Variscan mountains differ from previously established calibration sites (cf. Migoń and Waroszewski, 2022). Introducing a regional calibration site will likely allow us to calculate more accurate ¹⁰Be CRE ages for this region.

In this paper, we present a new calibration site in the southern Black Forest, SW Germany (Fig. 1). We chose the Feldsee cirque (47.9° N, 8.0° E; WGS 1984 coordinate reference system) because (i) we observed multiple large, quartz-bearing boulders on two well-preserved moraines and because (ii) a bog, the Feldsee Bog, is situated in the tongue basin of the former glacier whose sediments are stratigraphically younger than these ice-marginal moraines (Lang, 2005; Hofmann and Konold, 2023). We first measured the concentration of accumulated in situ cosmogenic ¹⁰Be in quartz from moraine boulder surfaces. Obtaining sediment cores from the FSM coring site (“FSM” stands for Feldseemoor, the German name of the bog) on a buried moraine, radiocarbon dating of macrofossils, infrared-stimulated luminescence (IRSL) dating, and establishing an age–depth model with the ¹⁴C ages, the IRSL ages, and the age of a cryptotephra allowed us to derive a minimum age for ice-free conditions at the bog. Based on the ¹⁰Be concentrations in rock samples from the moraine boulders and the modelled basal age of the bog's sediments, we deduced a regional SLHL ¹⁰Be production rate, hereinafter termed the Black Forest production rate (BFPR).

Figure 1 Topographical map of the southern Black Forest showing the assumed maximum ice extent during the Late Pleistocene (Hemmerle et al., 2016), ice divides (Hemmerle et al., 2016), and outlet glacier names according to the nomenclature of Hofmann et al. (2020). See NASA Jet Propulsion Laboratory (2013) for information on the digital elevation model (DEM) in the background. The inset map shows the location of the Black Forest in Germany.

2 Regional setting

The study site, the Feldsee cirque, is located in the southern part of the Black Forest in SW Germany (Fig. 1). The Feldsee cirque is situated about 2 km ESE of Feldberg (1493 m a.s.l.), the highest summit of the Black Forest. Due to the high abundance of glacial landforms (cf. Liehl, 1982; Metz and Saurer, 2012; Hofmann and Konold, 2023), it is a key site for Pleistocene glaciations of the Black Forest. The cirque has attracted glacio-geomorphological and geological research for almost 2 centuries (Walchner, 1846; Ramsay, 1862; Lang et al., 1984; Schreiner, 1990; Hofmann and Konold, 2023).

2.1 Study site

The Feldsee cirque has an impressive headwall up to 300 m high (Fig. 2; LGRB, 2023). The Feldsee, a moraine-dammed lake up to 33 m deep, covers the tongue basin on the cirque floor (Wimmenauer et al., 1990). The Feldsee Bog (LUBW, 2006), which is about 0.03 km² large, is located downstream from the ice-marginal moraine at the eastern shore of the lake (Fig. 2). The Feldsee cirque leads into the Seebach valley, a well-developed trough valley (Metz, 1985). The study site drains to the Seebach, originating from above the headwall of the Feldsee cirque. This stream pertains to the Rhine drainage network (LUBW, 2022a).

Figure 2 Oblique aerial photograph of the study site. Moraines are marked with dotted lines. The upper limit of the headwall of the Feldsee cirque is marked with a solid line. The Feldsee, a moraine-dammed lake, covers the tongue basin on the cirque floor. Note the prominent moraine at the eastern shore of the lake (position FS-01). The moraine at position FS-02 separates the Feldsee Bog further east into two distinct basins. The bog's sediments have partly buried this moraine. To the east, the bog is bounded by the semi-circular moraine at position FS-03. Also note the initial cirque (marked with a dashed line) on the northeastern flank of Seebuck. Photo: Matthias Geyer.

Quartz-bearing rock of the Variscan basement (age: 380–290 Ma; Geyer et al., 2011), i.e., flaser gneiss, migmatite, porphyry, and paragneiss, dominates the study area (LGRB, 2013). In addition, quartz-rich porphyry outcrops on the cirque's western headwall (LGRB, 2023). Since denudation from about 50 Ma onwards (Eberle et al., 2023) has led to the complete removal of the Permian, Triassic, and Jurassic sedimentary rock on the Variscan basement (Wimmenauer et al., 1990), glacial sediments at the study site (mainly till) originate from quartz-rich rock of Variscan age (Schreiner, 1990). During the Quaternary, the southern Black Forest repeatedly hosted ice caps and glaciers (Metz and Saurer, 2012; Hemmerle et al., 2016; Eberle et al., 2023). For the last (Late Pleistocene) glaciation, previous reconstructions have pointed towards the temporary existence of four interconnected ice caps (Fig. 1; Liehl, 1982; Metz and Saurer, 2012; Hemmerle et al., 2016), covering an area of about 1000 km² (Metz and Saurer, 2012). Figure 1 reveals that the ice cap on Feldberg and the surrounding region was the largest of these ice caps. During the last deglaciation, the ice caps disintegrated into valley glaciers and finally into isolated cirque glaciers (Metz and Saurer, 2012; Hemmerle et al., 2016; Hofmann et al., 2020, 2022, 2024b; Hofmann, 2023a). However, successive phases of ice-marginal stability punctuated the overall trend of glacier recession. According to ¹⁰Be CRE ages for the region NW of Feldberg, formerly covered by the Brugga outlet glacier and the western branch of the Seebach outlet glacier (Fig. 1), repeated phases of glacier recession from moraines occurred no later than 17–16, 14, and by 13 ka at the latest (Hofmann et al., 2022, 2024b; Hofmann, 2023a).

Although the Feldsee cirque glacier at the study site has long since disappeared, snow cover still plays an important role today. Between 1961 and 1990 CE, mean annual precipitation and average annual temperature at the weather station on Feldberg (at 1486 m a.s.l.), situated about 1.5 km west of the study site, amounted to 1909 mm and 3.3 °C, respectively (DWD, 2023). Snowfall accounted for about two-thirds of annual precipitation during this period (Matzarakis, 2012). As will be discussed below, the largely forest-covered Feldsee cirque is thus a challenging site where seasonal snow cover might have considerably slowed down the accumulation of ¹⁰Be in moraine boulder surfaces produced in situ.

2.2 Previous work

From a glacio-geomorphological perspective, the Feldsee cirque is undoubtedly one of the most interesting sites in the entire Black Forest. Both Walchner (1846) and Ramsay (1862) already recognized that glaciers shaped the Feldsee cirque and first described the prominent ice-marginal moraines east of the Feldsee (Fig. 3). These landforms document the Feldsee position, the penultimate period of repeated ice-marginal stability during the last deglaciation of the southern Black Forest (Steinmann, 1902; Schrepfer, 1925; Erb, 1948; Liehl, 1982; Metz and Saurer, 2012; Hofmann et al., 2020; Hofmann, 2023b). Hofmann and Konold (2023) recently revisited glacial landforms at the study site and performed geomorphological mapping at a 1:5000 scale with the aid of a digital elevation model (DEM) derived from light detection and ranging (lidar) data (x–y resolution: 1 m; vertical accuracy: ±0.2 m; LGL, 2015) and field evidence. As this is the most recent study on glacial landforms in the study area, we summarize the findings in the following paragraphs (see Fig. 3a).

As only the ice-marginal moraines of the former Feldsee cirque glacier are relevant for the geological calibration of the regional ¹⁰Be production rate, we only discuss ice-marginal positions of this glacier and not those (positions SH-01, SH-02, SH-03, SH-04, and SH-05) of a transfluent glacier that advanced from a southwesterly direction into the Seebach valley (Fig. 3a; Meinig, 1980).

The ice-marginal moraines of the Feldsee cirque glacier reflect three positions, named, from the oldest to the youngest, FS-03, FS-02, and FS-01 (Fig. 3a; Hofmann and Konold, 2023). A semi-circular and boulder-rich moraine with a sharp crest is located at position FS-03. This moraine partly overlies a moraine at position SH-05. Based on the SW–NE orientation of the latter, this moraine probably formed at the margin of a transfluent ice stream that advanced from the south, i.e., from the Seehalde area, from a southwesterly direction into the Seebach valley. As the moraine at position FS-03 overlies those at position SH-05, the semi-circular moraine at position FS-03 must have formed during a re-advance of the Feldsee cirque glacier. The moraine at position FS-02 subdivides the bog into two parts. Peat and lake sediments partly cover the moraine, and, therefore, it is only visible south of the bog, in the centre of the bog (Lang, 2005), and at its northeastern end. The semi-circular and sharp-crested moraine at position FS-01 dams up the Feldsee. According to Schreiner (1990), this landform is one of the morphologically most distinct and best-preserved moraines in the entire southern Black Forest.

Figure 3 (a) Topographical map of the study site showing the FSM coring site, moraines, and ice-marginal positions (Hofmann and Konold, 2023). “FSM” stands for Feldseemoor, the German name of the bog. For the DEM in the background, see LGL (2015). A diagrammatic eye indicates the photo location of Fig. 8c. (b) Transect from the area west of the Feldsee to the Seebach (4-fold vertical exaggeration) showing ice-marginal moraines and the coring site for a previous study (Lang, 2005). Note the shallower eastern basin of the Feldsee Bog. See panel (a) for the location.

In the 1970s and 1980s, Lang et al. (1984) obtained sediment cores at 13 sites along a W–E transect in the southern portion of the Feldsee Bog and five cores from the Feldsee (summarized by Lang, 2005). This research revealed the presence of the buried ice-marginal moraine at position FS-02 (Fig. 3b). Lang (2005) proposed that this landform subdivides the Feldsee Bog into two distinct basins. From the bottom to the top, Lang et al. (1984) identified the following sediments in the cores: till, gyttja, clayey gyttja/clay, gyttja, and peat. At one coring site (coring site “5” in Lang, 2005), these authors observed a distinct greyish layer at a depth of about 8.1 m below the ground surface (Fig. 4b).

Figure 4 (a) Evolution of the timberline in the Black Forest between 15 ka and 11 ka (redrawn from Lang, 2006). (b) Lateglacial part of the sediment sequence at a coring site in the Feldsee Bog (Fig. 3b), inferred ages of the sediments, and vegetation dynamics, as documented in the pollen record. LPZ: local pollen zones. Redrawn from Lang (2005).

Lang et al. (1984) speculated that this layer is the Laacher See Tephra, having a reported age of 13 006±9 cal. a BP (Reinig et al., 2021). Unfortunately, these authors did not properly document their observation, which would have more convincingly supported the cryptotephra interpretation. Also note that Lang et al. (1984) did not apply radiocarbon dating to their cores. These authors inferred from the succession of lithostratigraphic units and from the position of the apparent tephra and pollen data that a proglacial lake formed between the moraines at positions FS-03 and FS-01 after deglaciation. Sedimentation at the bottom of the lake probably began during the short-lived cooling episode at around 14 ka, often referred to as the Older Dryas (cf. Heiri et al., 2014). Later research revealed the presence of a greyish layer in the sediments of the Feldsee. Lang (2005) interpreted this layer as the Laacher See Tephra and concluded that both the bog and the lake must have been ice-free no later than 13.1 ka.

3 Methods

To the best knowledge of the authors, the cores obtained by Lang et al. (1984) unfortunately do not exist anymore. We thus obtained sediment cores at the FSM coring site during fieldwork in 2021 CE. In addition, we sampled moraine boulder surfaces and determined the concentration of accumulated in situ-produced ¹⁰Be in quartz.

3.1 Coring

Sediment cores at the FSM coring site were obtained with a COBRA vibracorer, allowing us to extract sediment cores with a diameter of 5 cm and a length of 1 m. To prevent daylight exposure of the sediments in the cores, the drilling system was equipped with opaque plastic tubes to store individual sediment cores. The cores were hydraulically extracted, sealed, and opened in the lab with a circular saw in red-light conditions. Sediment samples were obtained from the cores and dried, and loss-on-ignition (LOI) analyses (cf. Heiri et al., 2001) were undertaken. Weighing the samples prior to drying and before the LOI analyses allowed us to determine the sediments' water content.

During opening of the cores in the lab, we noted that all sediment cores were shorter than the penetrated depth; thus, core shortening must have occurred during vibracoring. Generally, core shortening is one of the main limitations of this technique (cf. Glew et al., 2001). As discussed by Glew et al. (2001), sediments with a higher water content are generally more prone to compaction. We assumed that only the clayey and silty lake sediments in the cores (water content: 18 %–85 %), and not the stratigraphically older diamicts (water content: 15 %–17 %), were affected by shortening. Following Glew et al. (2001), we assumed that the sediments in the cores were progressively thinned down-core, i.e., equally affected by compaction. Individual conversion factors were computed for every 1 m long sediment core which then allowed us to adjust the thickness of the lithostratigraphic units to the penetrated depth.

3.2 Radiocarbon dating

To numerically date the sediments at the coring site, radiocarbon dating of macrofossils was undertaken. Macrofossils were hand-picked after core surface cleaning and KOH pre-treatment of bulk samples (see Figs. S11 to S17 for photos of the macrofossils). A total of nine samples were sent to the radiocarbon laboratory in Poznań, Poland (Table 1), to obtain accelerator mass spectrometry (AMS) radiocarbon ages. The ¹⁴C ages were calibrated with the OxCal software (version 4.4; Bronk Ramsey, 2009), available at https://c14.arch.ox.ac.uk/oxcal/OxCal.html (last access: 18 September 2023), and the IntCal20 calibration curve (Reimer et al., 2020).

Table 1 Macrofossils sampled from sediment cores from the Feldsee Bog, conventional ages, and calibrated age ranges. Following the recommendation of Millard (2014), 95 % ranges of calibration are given.

3.3 IRSL dating

To cross-check the radiocarbon ages, five sediment samples from two of the cores (depth: 4–6 m below the ground surface) were sampled for luminescence dating under subdued red light, with two further samples taken to account for potential dose rate inhomogeneity due to the complex stratigraphy (FSM-D1 and FSM-D2; Table 2). The material for equivalent dose (D_e) determination was subsequently treated with 10 % HCl (no reaction) and 30 % H₂O₂ (partly heavy reaction) to remove carbonates and organic matter, respectively. For the lowermost two samples, feldspar from the sand fraction (90–200 µm, CG) and polymineral fine grains (4–11 µm, FG; based on grain-size distributions) were extracted. For the remaining three samples, only fine grains were extracted. Measurements were done on a Freiberg Instruments lexsyg research (Richter et al., 2013), equipped with an ET935QB photomultiplier (Hamamatsu Photonics) and using the combination of a Schott BG39 (3 mm) and a 414/46 BrightLine HC interference filter (3.5 mm) for detection.

Table 2 Samples for IRSL dating. Depth refers to the depth below surface after decompaction. Stim.: stimulation used for equivalent dose (D_e) determination. Given are the measured water content (W) and the sediment moisture (Mois.) used for dose rate calculations. The activity of different isotopes (K, Th, U-238, and Ra-226) is given, revealing disequilibrium in the Uranium decay chain. n: number of replicate D_e measurements. Model: applied approach to extract mean D_e (CAM: Central Age Model; MAM: Minimum Age Model).

It should be noted that optically stimulated luminescence measurements on quartz revealed no suitable signal, similar to reports on other directly bedrock-derived samples (e.g., Preusser et al., 2006) and experience from the nearby Upper Rhine Graben (Preusser et al., 2016, 2021). Therefore, feldspar was selected as a dosimeter. For all samples, a standard IRSL protocol (modified from Preusser, 2003) was applied. This protocol comprised a preheat to 250 °C for 60 s and IRSL stimulation at 50 °C for 100 s (IR-50). For fine grains, a post-IR (pIR) IRSL protocol was additionally tested to potentially overcome the need for fading correction. This protocol involved a preheat to 250 °C for 60 s, IRSL stimulation at 50 °C for 100 s, and a second stimulation at 225 °C for 100 s (pIR). For fine grains, the five replicate measurements were considered sufficient due to the excellent reproducibility, whereas for the sand fraction 20 replicate measurements were performed. For most samples, the average dose was calculated by the Central Age Model (CAM). In two cases, where the overdispersion exceeded 20 %, the Minimum Age Model (MAM) was allied with a sigma_b value of 0.20 (Galbraith and Roberts, 2012). For IRSL and IR-50, a fading correction was undertaken using a g value of 3.7±0.2 g per decade.

The concentration of dose-rate-relevant elements was determined using a high-resolution gamma spectrometer, as outlined in Preusser et al. (2023). An alpha efficiency of 0.07±0.02 and an internal K content of 12.5±1.0 % (Huntley and Baril, 1997) were assumed. Note that the setting is particularly challenging with regard to dose rate determination. While average sediment moisture during burial time was estimated based on the water content measured directly after opening the cores in the laboratory, it should be considered that a certain loss of water may have occurred during and after the coring operations. Furthermore, the deposits will have compacted after initial deposition due to the increase in load with time. Hence, the amount of water used in the calculations is higher than that measured but still only represents an approximation. All age calculations were performed with the ADELEv2017 software (Degering and Degering, 2020). As the different sediment layers are partly quite thin, three-layer models were used for all samples, as implemented in ADELEv2017. Namely, for each layer an individual dose rate was calculated from which the effective dose rate acting on the sample area was determined.

In addition, evidence for significant radioactive disequilibrium in the uranium decay chain was observed for most of the samples. This is common in organic-rich deposits (e.g., Preusser and Degering, 2007; Preusser et al., 2023) and reflects the absorption of uranium by humic substances from water (cf. Ivanovich and Harmon, 1992). This was accounted for in the dose rate calculations for each individual layer, assuming a constant uptake of uranium since the time of deposition. Both the applied layer model and the correction of radioactive disequilibrium represent only approximations to the real situation, as it is unfeasible to determine the detailed 3D geometry of the deposits and the exact timing of uranium uptake.

3.4 Age–depth modelling

To obtain a minimum age for the moraine at position FS-02, age–depth modelling with the online OxCal program and IntCal20 (Reimer et al., 2020) was undertaken by interrogating the ¹⁴C ages of macrofossils, the age of the Laacher See Tephra (13 006±9 cal. a BP; Reinig et al., 2021), and the IRSL ages (converted in cal. a BP). Two IRSL ages were available for the FSM4 and FSM5 samples, i.e., ages for both coarse grains and fine grains. Error-weighted IRSL ages and uncertainties of the error-weighted IRSL ages were computed prior to age–depth modelling. The P_Sequence function was selected in OxCal, as the ¹⁴C ages of the macrofossils and the IRSL ages of sediment samples were expected to increase with depth. In addition, this function allows changes in the sedimentation rate (Bronk Ramsey, 2008). Following the guidelines of Bronk Ramsey and Lee (2013), the model was tuned to find the most suitable value for the k factor, i.e., the number of increments per depth unit (see the Supplement for the input code for OxCal).

3.5 Regional ¹⁰Be production rate calibration and CRE dating of moraines

For establishing the BFPR, we collected surface rock samples (Table 3) from (i) six gneiss boulders on the moraine at position FS-03 and (ii) four gneiss boulders on the ice-marginal moraine at position FS-02. We also sampled the surface of the FS-01a gneiss boulder on the moraine at position FS-01a for age calculations. However, we did not include the sample in the calibration dataset, as moraine formation might have post-dated the onset of deposition of lake sediments at the FSM coring site (Fig. 3b). The modelled basal age of the Feldsee Bog then offered the opportunity to derive a regional SLHL spallogenic ¹⁰Be production rate.

3.5.1 Fieldwork, sample preparation, and measurements

Since the study of Tomkins et al. (2021) demonstrated that landform stability mainly influences the scatter in age distributions from moraines (and thus in ¹⁰Be concentrations), only well-embedded boulders were selected to avoid underestimated ¹⁰Be concentrations due to boulder rotating and to post-depositional and post-stabilization exhumation. As the moraine at position FS-03 consisted of clast-supported diamicts and some parts of the moraine were solely composed of boulders, identifying stable and large boulders proved to be straightforward. The same was true for the moraine at position FS-02, although this landform consisted of matrix-supported diamicts. Identifying large and stable boulders on the moraine at position FS-01 turned out to be difficult, as this landform consisted of matrix-rich diamicts and exhibited only a few boulders; thus we only sampled one large boulder. Hofmann and Konold (2023) mapped a kettle on the proximal side of the moraine at position FS-03 and on the moraine at position FS-02 in the centre of the Feldsee Bog (Fig. 3), pointing to paraglacial reworking and delayed moraine stabilization (cf. Porter et al., 2019). To minimize the risk for paraglacial reworking issues, we avoided sampling boulders in the vicinity of these landforms.

Table 3 Location of the sampled moraine boulders and their height, the thickness of the rock samples, and topographic shielding at the sampling sites.

Previous sampling guidelines suggested that flat-topped boulders should be selected for CRE dating of moraines (e.g., Ivy-Ochs and Kober, 2008). However, none of the moraine boulders in the study area were flat-topped; thus rock samples (Table 4) were obtained with an angle grinder, a chisel, and a hammer from dipping surfaces with a constant angle. See the detailed sample documentation (Tables S1 to S11 in the Supplement) for the strike and dip of the sampling surfaces, measured with a geological compass. Masarik and Wieler (2003) argued that the cosmic ray flux on non-flat rock surfaces is lower than on flat rock surfaces, since the number of scattered neutrons is larger on non-flat rock surfaces. Due to these potential edge effects, rock surfaces near boulder edges were avoided during fieldwork. As scaling of the ¹⁰Be production rate to sampling sites required coordinates of the sampling sites, x–y coordinates of the boulders were recorded with a global navigation satellite system (Leica CS20 controller and Leica Viva GS14 antenna).

The elevation of the sampling sites was later retrieved from the abovementioned high-resolution digital elevation model of the study site of the Baden-Württemberg State Agency for Spatial Information and Rural Development (LGL, 2015). Obtaining pairs of azimuth and elevation angles of the horizon around the sampling sites for topographic shielding factor calculations was not possible at all sampling sites due to the presence of dense mixed forests. Therefore, the ArcGIS toolbox of Li (2018) was chosen for shielding factor calculations, considering (i) self-shielding of dipping surfaces and (ii) shielding by topographical obstructions around the sampling sites. Following the guidelines of Hofmann (2022), the high-resolution DEM of the study site (LGL, 2015) was resampled to an x–y resolution of 30 m for shielding factor calculations. The mass and thickness of each rock fragment in the samples were determined to compute the mass-weighted average of the sample's thickness. Since the sampled boulders were rich in quartz, the density of quartz (2.65 g cm⁻³) was assumed for the thickness correction. For photos of the boulders, dimensions of the boulders, and height above ground of the sampling surfaces, see the detailed sample documentation in the Supplement.

Table 4 Characteristics of the samples: ⁹Be added, the results of AMS measurements, blank-corrected ¹⁰Be concentrations, and scaling factors for the samples in the calibration dataset (geomagnetic database: Muscheler et al., 2005). ¹⁰Be concentrations in the samples were corrected with the ¹⁰Be concentration in a batch-specific blank.

^a Analytical blank. ^b The ¹⁰Be concentration was classified as an outlier with Peirce's criterion (Peirce, 1852). ^c Due to low ⁹Be currents during AMS measurements, the reported ¹⁰Be concentration should not be regarded as reliable. ^d This concentration was not used for production rate calibration.

Quartz separation and ¹⁰Be extraction from purified quartz were accomplished in the laboratory facilities of the University of Freiburg (Germany) and the Laboratoire National des Nucléides Cosmogéniques (LN₂C) in Aix-en-Provence (France) according to the protocol described in Hofmann et al. (2024b). After crushing and wet sieving (target grain size: 0.25–1 mm), the samples were passed through a magnetic separator (S.G. Frantz Co.). The samples were then treated with mixtures of 37 % HCl and 35 % H₂SiF₆ to further isolate quartz and remove feldspars. Because the samples still contained feldspar after this treatment, the samples were subsequently etched with diluted 5.5 % HF, dried, spiked with magnetite powder (325 mesh), and passed through a magnetic separator. Meteoric ¹⁰Be was removed with 48 % HF in three steps, with 10 % of the quartz dissolved in each step. A dose of about 150 mg of a ⁹Be carrier solution (3025±9 µg ⁹Be g⁻¹) was added to each sample before total dissolution with 48 % HF. Chromatography with anionic and cationic exchange resins (DOWEX 1X8 and 50WX8) and precipitation stages were then performed to further separate and purify beryllium. To thermally decompose the final Be(OH)₂ precipitate to BeO, the samples were subsequently heated to 700 °C. The samples were finally mixed with Nb and pressed into copper cathodes for AMS measurements.

The ¹⁰Be concentrations in the samples (Table 4) were deduced from AMS measurements at ASTER (Accélérateur pour les Sciences de la terre, environnement, Risques; Arnold et al., 2013), the French AMS national facility, at CEREGE (Centre Européen de Recherche et d'Enseignement des Geosciences de l'Environnement) in Aix-en-Provence (France). The measured ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios were normalized with respect to the in-house standard STD-11 using an assigned ${}^{\mathrm{10}}\mathrm{Be}/{}^{\mathrm{9}}\mathrm{Be}$ ratio of (1.191±0.013) × 10⁻¹¹ (Braucher et al., 2015) and the Be half-life of (1.387±0.012) × 10⁶ years (Chmeleff et al., 2010; Korschinek et al., 2010). The ${}^{\mathrm{10}}\mathrm{Be}/{}^{\mathrm{9}}\mathrm{Be}$ ratio uncertainties comprise (i) the measurement uncertainty (counting statistics), (ii) the error of average standard measures, and (iii) the systematic error of ASTER (0.5 %; Arnold et al., 2010). The ¹⁰Be concentrations in the samples were corrected for the ¹⁰Be concentration in a batch-specific chemical blank (Table 4).

3.5.2 Production rate calibration

The calibration of the spallogenic ¹⁰Be SLHL Black Forest production rate followed the workflow of Martin et al. (2017; their Fig. 3). Following Eq. (4) in Martin et al. (2017), the ¹⁰Be concentration in each sample was first corrected for the sample thickness and for topographic shielding to compute a theoretical ¹⁰Be concentration in a sample with null thickness. The ¹⁰Be concentrations were then scaled to the mean latitude, longitude, and elevation representative for all sampled boulders (Martin et al., 2017). Following the guidelines of Ross (2003), the ¹⁰Be concentrations were subsequently evaluated with Peirce's criterion (Peirce, 1852), and a weighted ¹⁰Be concentration was computed after the exclusion of outliers. The mean squared weighted deviation (MSWD) for the ¹⁰Be concentrations turned out to be higher than 1. As recommended by Martin et al. (2017), the standard error of the weighted mean ¹⁰Be concentration (calculated with 1σ uncertainties of the ¹⁰Be concentrations) was multiplied with $\sqrt{\text{MSWD}}$ to obtain the uncertainty of the average ¹⁰Be concentration.

The regional ¹⁰Be production rate was calculated with the aid of CREp, available at https://crep.otelo.univ-lorraine.fr (last access: 6 November 2023). For calibration, the modelled basal age (15 640±420 cal. a BP) of the lake sediments above the moraine at ice-marginal position FS-02 was converted into ka before 2010 CE. The spallogenic SLHL ¹⁰Be BFPR in quartz was computed for the scaling schemes in CREp, i.e., time-dependent Lal–Stone (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008) and LSD (Lifton et al., 2014) scaling, and all geomagnetic databases in CREp, i.e., the atmospheric ¹⁰Be-based virtual dipole moment (VDP; Muscheler et al., 2005, and references therein), the LSD framework (Lifton et al., 2014), and the Lifton 2016 VDM (Lifton, 2016, and references therein). Following the approach in a previous calibration study (Fenton et al., 2011), a “baseline” production rate was first calculated, i.e., a production rate that accounts for the site-specific bias induced by snow cover, vegetation cover, soil cover, and post-depositional weathering and by other factors, such as changes in atmospheric circulation.

As many geomorphologists prefer to work with version 3 of the online exposure age calculator, formerly known as the CRONUS-Earth online exposure age calculator (available at http://hess.ess.washington.edu/math/v3/v3_cal_in.html, last access: 8 November 2023), the “baseline” production rate was additionally determined with this calculator. See the Supplement for the input sheet.

Fenton et al. (2011) argued that snow cover and post-depositional weathering of sampling surfaces probably affected their regional spallogenic SLHL ¹⁰Be production rates at two rock avalanches in Norway and proposed post-depositional-weathering- and snow-cover-corrected regional spallogenic SLHL production rates. The sampling sites in the southern Black Forest were situated in a heavily forested and sheltered area. A layer of organic soil, mosses, and shrubs, most notably blueberry bushes (Vaccinium myrtillus), a few centimetres thick, covered the sampled boulders. As soil cover on boulders leads to enhanced chemical weathering rates (cf. Ahnert, 2009), the boulders probably underwent significant post-depositional weathering. In addition, snow cover during winter has probably slowed down the accumulation of ¹⁰Be in quartz in the sampled boulders. During field surveys in 2022 CE, we observed a substantial snow cover in winter on the boulders, being up to a few decimetres thick.

Therefore, proposing a post-depositional-weathering- and snow-cover-corrected regional ¹⁰Be production rate was deemed mandatory. On the FS-02b boulder, we observed a protruding quartz vein with a height of 1 cm. The height of the quartz vein and the abovementioned modelled basal age of the lake sediments at the FSM coring implied a post-depositional weathering rate of 0.06 cm ka⁻¹ for the weathering-corrected BFPR. It should be noted that this rate of weathering and removal of rock is lower than a previous rate for a similar setting in the Bavarian Forest (0.24 cm ka⁻¹; recalculated from Reuther, 2007). According to data from the weather station of the DWD in the municipality of Schluchsee (at 990 m a.s.l.), situated about 10 km to the SE of the study area, seasonal cover lasted for 4 months in the 1961–1990 CE period (average snow depth: 0.3 m; DWD, 2023). Assuming similar conditions during the whole duration of exposure, a snow density of 0.3 g cm⁻³, and an attenuation length for fast neutrons in snow of 109 g cm⁻² (Zweck et al., 2013), the commonly used Eq. (3.76) in Gosse and Phillips (2001) yielded a snow shielding factor of ∼0.974. This factor was adopted for production rate calibration.

Modelling has revealed up to a 7 % reduction in the intensity in cosmic ray flux in old-growth temperate forests (Plug et al., 2007). As the sampled moraine boulders were situated in forested areas, a vegetation-corrected production rate was calculated. According to palynological data from several bogs and fens in the Black Forest, the timberline rose to an elevation of ≥1100 m a.s.l. during the last interstadial of the Late Pleistocene (Fig. 4a; Lang, 2006), commonly referred to as the Allerød interstadial (13.9–12.8 ka; Heiri et al., 2014). During this period, open pine and birch forests probably covered the study site (Fig. 4b). The presence of Pinus in the pollen record, however, does not necessarily indicate that the study site was forested. As discussed by Andrieu et al. (1997), Pinus generally tends to be overrepresented in pollen records due to high pollen production and winds easily transporting the pollen over long distances. Unfortunately, Lang (2005) did not observe Pinus or Betula macro-remains that would have more convincingly supported the idea of open pine and birch forests. Lang (2006) suggested that the timberline shifted back to about 750 m a.s.l. during the subsequent cold phase (Fig. 4a), often referred to as the Younger Dryas stadial (12.8–11.7 ka; Heiri et al., 2014). However, with regard to the northern Alps where the timberline lowered by about 200 m during this cold phase (see Stojakowits et al., 2014), the reconstructed shift of the timberline in the Black Forest (ca. 350 m) appears quite strong. It is therefore possible that the study site remained a sparsely forested area (Philipp Stojakowits, personal communication, 25 October 2023). Since we could not conclusively clarify the Lateglacial vegetation history of the study site, we assumed for the vegetation correction that forests covered the study area throughout the duration of exposure. We are, however, aware that this correction might overestimate the impact of the vegetation on the cosmic ray flux at the sampling sites. For the vegetation correction, we assumed a shielding factor of 0.98 for boreal forests (Plug et al., 2007), probably the best analogues for the forests at the study site.

3.6 Assessment of the impact of the new production rate

To assess the impact of the newly calibrated BFPR on CRE ages, CRE ages, internal (analytical) uncertainties, and external uncertainties (i.e., analytical uncertainties plus the error of the ¹⁰Be production rate added in quadrature) for the sampled moraine boulder surfaces were calculated with the Chironico landslide spallogenic SLHL production rate and the BFPR. The following parameters were chosen to scale the SLHL production rates to the sampling sites: time-dependent “Lm” scaling (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008) and the ERA atmosphere model (Uppala et al., 2005), as recommended by Martin et al. (2017). The ages were corrected for changes in the Earth's magnetic field with data from the atmospheric ¹⁰Be-based geomagnetic database (Muscheler et al., 2005). Since the sampled moraine boulders were derived from quartz-rich lithologies, the density of quartz (2.65 g cm⁻³) was assumed for the thickness correction.

3.7 Statistical assessment of age datasets

Multiple ages were available for most of the moraines. The assessment of the ages followed the guidelines of Balco (2023). If three or more ages were available from the same landform, reduced chi-squared (${\mathit{\chi }}_{\mathrm{R}}^{\mathrm{2}}$) was computed. ${\mathit{\chi }}_{\mathrm{R}}^{\mathrm{2}}$ was then compared with a critical value from a standard χ²-table (degree of freedom: n−1). The confidence interval was set to 95 %. If ${\mathit{\chi }}_{\mathrm{R}}^{\mathrm{2}}$ turned out to be lower than the critical value, the hypothesis that the data formed a single population was at 95 % confidence (cf. Balco, 2011). In this case, the landform age was determined by computing the error-weighted mean. If ${\mathit{\chi }}_{\mathrm{R}}^{\mathrm{2}}$ turned out to be higher than the critical value, implying that measurement uncertainties did not fully account for scatter in ages from the same landform (cf. Balco, 2011), the age that was furthest from the average age in relation to its measurement uncertainty was considered an outlier and removed from the dataset. This procedure was repeated until the pruned dataset yielded an acceptable ${\mathit{\chi }}_{\mathrm{R}}^{\mathrm{2}}$ value. No more than half of the ages from the same landform were excluded. The procedure was also stopped if three ages were remaining. If the original or pruned datasets yielded a p value >0.05, the internal landform uncertainty was determined by calculating the standard error. If the original or pruned datasets gave a p value <0.05, the standard deviation of the ages was selected as the internal landform age uncertainty. The external landform age uncertainties were calculated by adding the internal landform age uncertainties and the error of the selected production rate in quadrature.

4 Results

Vibracoring at the FSM coring site allowed us to obtain sediment cores with a total length of 8 m. Borehole FSM recovered the sedimentary succession of the Feldsee Bog and the uppermost 0.32 m of the partly buried moraine at position FS-02 (Fig. 5). The percentage of sediment recovery increased from 67 % (depth: 4–5 m) to 81 % (depth: 5–6 m). Decompaction was thus undertaken with correction factors of 1.49 (4–5 m) and 1.39 (5.00–5.68 m). Beryllium-10 concentrations in a total of 10 moraine boulder surfaces were successfully determined, allowing us to calibrate the production rate.

Figure 5 Sediment sequence at the FSM coring site, calibrated ages (95 % ranges of calibration), IRSL ages, the age–depth model, and LOI. Photos of the sediment cores from the Feldsee Bog were acquired with the methodology of Gegg and Gegg (2023). 95 % ranges of the modelled ages are marked with dashed lines. The solid line represents the mean modelled age. Agreement indices (A) for individual ages are given in parentheses.

4.1 Lithostratigraphy of the sediment sequence

In the lowermost 0.32 m, this succession consisted of greenish-grey, medium-densely packed silty diamicts that were mostly clast-supported and occurred in a fining-upwards package. Individual clasts were angular to subangular and of diverse crystalline lithologies of local origin (gneiss, migmatite, and granite porphyry; LGRB, 2013). With 1.2 % to 1.9 %, LOI values in samples from the sandy diamicts at a depth of 5.80–5.68 m were very low. At a depth of 5.68 m, the diamicts transitioned into crudely laminated light-grey to olive-grey fines, silt, and clay that contained dispersed coarse sand to fine gravel clasts at their base. LOI values varied between 1.7 % and 10.5 %. Dark-brown to grey fines, silt, and clay overlay these deposits. LOI values ranged from 12.3 % to 35.7 %. Macroscopic plant fragments and a distinct 2–3 mm thick light-grey lamina occurred in this unit. Tephra observed in non-visible ash beds (cryptotephra; cf. Krüger and van den Bogaard, 2021) at a decompacted depth of 4.91 m were glassy, generally <100 µm in size, and rich in fluid inclusions and vesicles (Fig. 5). They were isotropic and often conchoidally fractured. Oftentimes the tephra shards featured elongated fluid inclusions that were directionally aligned. Light-beige grey fines (sand, silt, and clay) with low LOI values (4.9 %–5.6 %) overlay the dark-brown to grey fines. The basal contact of this unit was sharp. Starting at a depth of 4.75 m, organic-rich (LOI: 30.5 %–54.8 %) gyttja-type deposits that were crudely bedded on a centimetre scale occurred (Fig. 5).

The succession of FSM (Fig. 5) reflects the glacial–postglacial transition of the study area. The FSM borehole was situated in elongation of a geomorphologically distinct ice-marginal moraine (Fig. 3a) and recovered its continuation into the Feldsee Bog. The corresponding diamicts (6.00–5.68 m) were medium-densely packed and contained washed-out sections with decreased fines content. This indicates an origin from melt-out and subaqueous mass movements (e.g., Schlüchter, 1997). The decreasing gravel content in this cycle further suggests a progressive shift from an immediately ice-proximal to a slightly more distal position due to glacier retreat. Well-sorted basin fines abruptly overlay the basal diamicts, reflecting a cessation of direct glacial input, although dispersed outsized clasts near the bottom have likely been deposited as ice-rafted debris. Frequently occurring plant fragments and fibres, and higher LOI values, indicate that a milder climate prevailed during the emplacement of the sediments. This study confirmed the assertion of Lang (2005) that a cryptotephra occurs in the sedimentary sequence of the Feldsee Bog. Future work will have to determine the major element composition of the tephra shards to unequivocally link it to the eruption of the Laacher See volcano. A mostly inorganic package at 5.00–4.75 m followed the first organic-rich package from 5.28 to 5.00 m, possibly due to cooler climatic conditions. The emplacement of the gyttja-type sediments above 4.75 m depth occurred during persistently temperate conditions.

4.2 Radiocarbon ages

A total of nine macrofossils in sediment cores from the FSM coring site were radiocarbon dated. Table 1 summarizes conventional ¹⁴C ages and calibrated radiocarbon ages (see also Fig. 5). Following the recommendation of Millard (2014), we hereinafter refer to 95 % ranges of calibration.

Calibrating ¹⁴C ages of unidentified plant fragments sampled at decompacted depths of 5.61, 5.56, and 5.54 m resulted in overlapping ages of 15 600–15 210, 15 480–15 040, and 15 680–15 240 cal. a BP, respectively. Sphagnum moss at a decompacted depth of 5.53 m yielded a modern ¹⁴C age. Leaf fragments at a depth of 5.43 gave an age of 13 730–13 170 cal. a BP. Sphagnum moss at a decompacted depth of 5.26 m yielded a ¹⁴C age of 720±30 a BP, whereas unidentified plant fibre at a decompacted depth of 5.24 m gave a modern ¹⁴C age. In contrast to the FSM-563, FSM-560, and FSM-558 samples, the appearances of sphagnum moss (FSM-553 and FSM-538 samples) and the plant fibre (FSM-536) were comparably fresh; thus the macrofossils are probably modern samples that were squeezed into older sediments during vibracoring. As these macrofossils were apparently modern samples, their ¹⁴C ages were not included in the age–depth model. Macrofossils sampled at a depth of 4.28 m gave almost identical ages of 9750–9530 and 9750–9490 cal. a BP, respectively (Table 1).

4.3 Luminescence ages

The IRSL ages of seven samples from the sediment cores from the FSM coring site ranged from 10.2±0.7 to 15.8±1.7 ka (Fig. 5). Table 4 summarizes the luminescence data for the FSM coring site.

While in the lower part of the sequence (samples FSM3, FSM4, and FSM5), the pIR ages are some 25 % higher than the fading corrected IRSL and IR-50 ages, there is a good agreement between the different approaches for the two samples taken just above the cryptotephra (samples FSM1 and FSM2). The mean age for these latter samples is 11.4±0.4 ka, which is some 10 %–15 % lower than expected for sediments above the cryptotephra, provisionally classified as Laacher See Tephra (13 006±9 cal. a BP; Reinig et al., 2021). FSM3 is also slightly to moderately underestimating the reported tephra age with estimates of 10.2±0.7 ka (IRSL) and 11.7±1.0 ka (IR-50). The ages for sample FSM4 (FG: IRSL = 12.7±1.0 ka and IR-50 = 13.2±1.0 ka; CG: 12.7±0.9 ka) overlap with the radiocarbon age of ca. 13.4 ka. The same applies to the basal sample FSM5, with ages of 14.0±1.2 ka (FG: IRSL), 14.2±1.2 ka (FG: IR-50), and 15.8±1.7 (CG: IRSL). Overall, there is a tendency of slightly lower IRSL/IR-50 ages compared to independent age control through radiocarbon dating (Fig. 5). This is possibly explained by an assumed sediment moisture that is slightly too low due to compaction of the deposits with time, uncertainties related to the applied layer model, and/or correction of radioactive disequilibrium (timing of uranium uptake). On the other hand, some the pIR ages appear slightly overestimated, which might be related to small residual signal levels at the time of deposition. However, the consistent results determined for sample FSM4 indicate that fading correction is appropriate. Note that Fig. 5 only displays the IRSL ages for simplification.

4.4 Age–depth model

Figure 5 shows the age–depth model for the depth between 5.68 and 4.00 m. Including all ¹⁴C and IRSL ages, the reported age of the cryptotephra gave an agreement index of 0.1 %, being well below the critical threshold of 60 % advocated by Bronk Ramsey (2008). The IRSL ages for the FSM2, FSM3, and FSM4 samples were excluded during the second and subsequent runs, as individual agreement indices turned out to be lower than 60 %. The final age–depth model gave an agreement index of 76.3 %. The modelled age (mean age and 1σ uncertainty) and the 95 % modelled age range for a decompacted depth of 5.68 m turned out to be 15 640±420 cal. a BP (15.70±0.42 ka before 2010 CE) and 16 470–15 230 cal. a BP, respectively. Hence, the deposition of lake sediments commenced at 16 470–15 230 cal. a BP, thus providing the minimum age of ice-free conditions which will just post-date glacier recession from positions FS-02 and FS-03.

4.5 Beryllium-10 concentrations in rock samples from moraine boulders

The blank-corrected beryllium-10 concentrations in moraine boulder surfaces ranged from 92 300±4200 to 161 000±13 000 atoms g⁻¹ quartz. Table 4 summarizes the results of AMS measurements and blank-corrected ¹⁰Be concentrations. Figure 6 displays normalized ¹⁰Be concentrations with respect to the normalized error-weighted mean ¹⁰Be concentration. During AMS measurements on the FS-03b sample, ⁹Be currents were very low, decreasing from 0.1 to 1 µA during AMS measurements. Hence, the ¹⁰Be concentration in the sample from the boulder should be considered with greatest caution.

Figure 6 Normalized ¹⁰Be concentrations in rock samples from the study site. The ¹⁰Be concentrations have been corrected for thickness and topographic shielding, scaled to SLHL using time-dependent “Lm” scaling, and normalized to the error-weighted mean ¹⁰Be concentration utilized for production rate calibration. The purple curves are Gaussian approximations of ¹⁰Be concentrations in individual samples. The vertical purple band and the vertical light-purple band represent 1σ and 2σ uncertainties of the error-weighted mean ¹⁰Be concentration, respectively. The vertical purple line corresponds to the error-weighted mean ¹⁰Be concentration.

4.6 Regional ¹⁰Be production rate

Calibrating the spallogenic BFPR with the ¹⁰Be concentrations in the FS-02a, FS-02b, FS-02c, FS-03a, FS-03c, FS-03d, FS-03e, and FS-03f moraine boulder surfaces resulted in SLHL production rates between 3.61±0.11 and 3.65±0.11 ¹⁰Be g⁻¹ quartz a⁻¹ for the different scaling schemes and geomagnetic databases in CREp (Table 5).

The ¹⁰Be concentration in the FS-03b moraine boulder surface was not used for production rate calibration due to low ⁹Be currents during AMS measurements. Correcting the ¹⁰Be concentrations for topographic shielding and the sample thickness led, in most cases, to minor shifts in the ¹⁰Be concentrations in samples. During the evaluation of the set of ¹⁰Be concentrations with the criterion of Peirce, the ¹⁰Be concentration in the sample from the FS-02d moraine boulder turned out to be an outlier. Hence, weighted mean concentrations were computed with the ¹⁰Be concentrations in the remaining samples. In all cases, MSWD exceeded the critical value of 1 (see Martin et al., 2017, for further discussion). Following Martin et al. (2017), the standard error of the weighted mean and $\sqrt{\text{MSWD}}$ were added in quadrature to obtain the uncertainties of the weighted mean ¹⁰Be concentrations. See Table 4 for the ¹⁰Be concentration in moraine boulder surfaces that have been corrected for topographic shielding and sample thickness and have been adjusted to SLHL.

The “baseline” production ranged from 3.61±0.11 to 3.65±0.11 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL for the different scaling schemes and geomagnetic databases in CREp (Table 5). These production rates resembled the production rate for “Lm” (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008) and “St” (Stone, 2000) scaling (Table 6), calculated with version 3 of the online exposure age calculator, formerly known as the CRONUS-Earth online exposure age calculator (Balco et al., 2008).

Removing the snow cover bias at the study site led to 2.7 % higher production rates between 3.70±0.12 and 3.75±0.12 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL. The snow-shielding- and post-depositional-weathering-corrected production rates ranged from 3.73±0.12 to 3.78±0.12 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL and were thus 6 % higher than the “baseline” production rates. With 3.81±0.12 to 3.86±0.12 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL, the post-depositional-weathering-, snow-cover-, and forest-cover-corrected production rates were 8 % higher than the “baseline” production rates (Table 5).

Table 5 The spallogenic ¹⁰Be SLHL BFPR (atoms ¹⁰Be g⁻¹ quartz a⁻¹) for different scaling schemes and geomagnetic databases available in CREp (Martin et al., 2017).

Table 6 The spallogenic ¹⁰Be SLHL BFPR for two scaling schemes in version 3 of the online exposure age calculator, formerly known as the CRONUS-Earth online exposure age calculator (Balco et al., 2008).

4.7 Beryllium-10 CRE ages

Calculating the ¹⁰Be CRE ages of moraine boulder surfaces with the BFPR resulted in ages between 10.2±0.5 and 17.6±1.5 ka. CRE ages calculated with the Chironico landslide spallogenic SLHL production rate lay between 9.1±0.5 and 15.7±1.3 ka. See Table 6 and Fig. 7 for CRE ages of moraine boulders, landform ages, and reduced χ² values.

Table 7 CRE ages of sampling surfaces on moraine boulders computed with the Black Forest (this study) and Chironico landslide (Claude et al., 2014) production rates.

^a Classified as an outlier. ^b Due to low ⁹Be counts during AMS measurements, the ¹⁰Be concentration in the sample is not reliable; therefore, no CRE age was computed.

Table 6 reveals that the ages calculated with the BFPR turned out to be, on average, 12.3 % higher than those computed with the Chironico landslide production rate. We hereinafter only comment on the ages and 1σ uncertainties derived with the BFPR. External uncertainties are given in parentheses.

Figure 7 Map of moraines at the study site and the immediate surroundings. CRE ages and associated external uncertainties of moraine boulder surfaces are given in ka before 2010 CE and ka, respectively. Bold ages are landform ages and landform age uncertainties (i.e., external uncertainties). Ages in grey are classified as outliers. See NASA Jet Propulsion Laboratory (2013) for a description of the DEM in the background.

The FS-03a, FS-03c, FS-03d, FS-03e, and FS-03f boulders on the ice-marginal moraine at position FS-03 yielded ages of 17.59±1.39 (1.48) ka, 15.55±0.48 (0.66) ka, 14.75±0.50 (0.66) ka, 16.34±0.59 (0.76) ka, and 16.34±0.68 (0.83) ka, respectively. As mentioned above, the ¹⁰Be concentration in the sample from the FS-03b boulder should not be regarded as reliable due to low ⁹Be currents during Be AMS measurements; thus, no age was calculated. The age of the FS-03d boulder (14.75±0.50 (0.66) ka) was classified as an outlier. The remaining ages gave a landform age of 16.06±0.42 (0.64) ka (Fig. 7 and Table 7).

The FS-02a, FS-02b, FS-02c, and FS-02d boulders on the moraine at position FS-02 gave ages of 15.52±0.59 (0.74) ka, 15.58±0.57 (0.72) ka, 15.62±0.48 (0.67) ka, and 10.24±0.45 (0.54) ka, respectively. After the exclusion of the outlying age of the FS-02d boulder (10.2±0.5 (0.5) ka), the remaining ages yielded a landform age of 15.58±0.03 (0.47) ka (Table 6 and Fig. 7).

The only sufficiently large and stable boulder on the ice-marginal moraine at position FS-01 was exposure dated to 14.94±0.50 (0.66) ka (Table 6 and Fig. 7).

Overall, the presented ages were, apart from a few outliers, internally consistent, and the landform ages complied with the stratigraphy. With 15.4±0.8 ka, 15.4±0.8 ka, and 15.4±0.7 ka, the ages of the FS-02a, FS-02b, and FS-02c boulders on the moraine at position FS-02 were remarkably consistent (reduced χ²: 0.0; Table 6).

5 Discussion

The determination of concentration of in situ-accumulated cosmogenic ¹⁰Be in boulders on moraines in the Feldsee cirque and applying radiocarbon dating, tephrochronology, and IRSL dating to stratigraphically younger lake sediments allowed us to establish the internally robust BFPR production rate, the first ¹⁰Be production rate for the mid-elevation mountain ranges of central Europe. In Sect. 5.1, we discuss the internal robustness of this production rate. We then evaluate previously suggested corrections of CRE ages for post-depositional weathering and snow shielding (Sect. 5.2). In Sect. 5.3, we finally discuss the newly calibrated BFPR in the European context. A detailed discussion of the implications of this study for the Late Pleistocene glacial history of mid-elevation mountain ranges would go beyond the scope of this paper and will be discussed in future work.

5.1 Internal robustness of the BFPR

The modelled basal age of the lake sediments above the moraine at the ice-marginal position FS-02 (mean age and 1σ uncertainty: 15 640±420 cal. a BP) underpins the SLHL BFPR. Radiocarbon dating of macrofossils near the base of the lake sediments, i.e., at decompacted depths of 5.61, 5.56, and 5.54 m, resulted in consistent ages of 15 600–15 210, 15 480–15 040, and 15 680–15 240 cal. a BP, respectively. These ages agree with the concordant IRSL ages of the FSM5-CG and FSM5-FG samples, luminescence dated to 15.8±1.7 and 14.0±1.2 ka, respectively. Due to consistent ages obtained with two independent dating methods, we regard the modelled basal age of the lake sediments (15 640±420 cal. a BP) as reliable.

The calibration of the BFPR was based on two implicit assumptions. Firstly, we deemed the modelled basal age at the FSM coring site a minimum age for the ice-marginal moraines at positions FS-03 and FS-02. The time lag between glacier recession from the moraines at these positions and the deposition of the macrofossils in the former proglacial lake remains unknown. Unfortunately, no bog exists downstream from the position FS-03 that would have allowed us to determine a maximum radiocarbon or luminescence age of the moraines at position FS-03. We were thus unable to derive the lowest-possible ¹⁰Be production rate at the study site.

Secondly, we implicitly assumed that all sampled moraine boulder surfaces had already been exposed to cosmic radiation when the macrofossils were deposited at the coring site. Previous work (Heyman et al., 2016) has revealed that larger boulders on moraines are more likely to have been continuously exposed to cosmic radiation than smaller boulders. Plotting the boulder height versus the normalized ¹⁰Be concentrations, i.e., the scaled ¹⁰Be concentrations, revealed the absence of a relationship between these variables (Fig. 8a). The low and statistically insignificant coefficient of determination (R²=0.03; p=0.84) contradicts the hypothesis that post-depositional and post-stabilization exhumation is a serious issue at the study site.

Figure 8 (a) Boulder heights versus normalized ¹⁰Be concentrations. (b) Protruding quartz veins on a random boulder on the ice-marginal moraine at position FS-03. (c) Boulder-rich ice-distal side of the moraine at position FS-03 where the large FS-03b and FS-03c boulders were sampled. All photos: Felix Martin Hofmann.

Previous work has also demonstrated that landform stability critically determines CRE ages of moraine boulder surfaces and that the position of the boulders on the landform, i.e., the moraine crest, the ice-proximal slope, or the ice-distal slope, does not matter (cf. Tomkins et al., 2021). Figure 8c shows a photo of the ice-distal side of the moraine at position FS-03 where the FS-03a, FS-03b, and FS-03c boulders have been sampled. Note that this part of the moraine almost entirely consists of boulders. We therefore deem significant post-depositional reworking of the moraine unlikely. It should be noted that, despite their high age, the moraines at positions FS-03 and FS-02 still have distinct crests. We interpret the comparably high landform stability as a further indication that the sampled moraine boulders had already been exposed to cosmic radiation when the deposition of lake sediments at the FSM drilling site began.

5.2 Evaluation of previously suggested corrections for snow shielding and post-depositional weathering

Due to the lack of a regional calibration site for the mid-elevation mountains of central Europe, CRE ages presented in previous studies were determined with SLHL ¹⁰Be production rates at localities outside this region and subsequently corrected for snow shielding and weathering (e.g., Mercier et al., 1999; Reuther, 2007; Hofmann et al., 2022, 2024b). Note that this approach is technically wrong, as SLHL production rates at calibration sites have their own bias with regard to weathering, snow shielding, and other factors. Previous studies thus corrected, to some extent, twice for snow shielding and weathering.

Irrespective of this discussion, the “baseline” BFPR offers the unique opportunity to critically assess whether these corrections were sufficient. For example, Hofmann et al. (2022, 2024b) calculated CRE ages for the southern Black Forest with the SLHL ¹⁰Be Chironico landslide production rate (Claude et al., 2014), as this calibration site was the closest calibration site to the southern Black Forest. These authors corrected CRE ages for the supposedly stronger snow shielding bias in the Black Forest with an additional snow shielding factor (∼0.974), deduced from meteorological data obtained between 1961 and 1990 CE at the weather station near Lake Titisee, located at about 850 m a.s.l. and about 10 km to the NE of the study area. The correction for the supposedly stronger post-depositional weathering bias was based on a weathering rate for homogenous paragneiss from the Bavarian Forest (0.24 cm ka⁻¹; recalculated from Reuther, 2007). Correcting the ages resulted in an age shift of about 5 % (Hofmann et al., 2022, 2024b).

Figure 9 shows the uncorrected and snow-shielding- and weathering-corrected ¹⁰Be ages for the study site (Chironico landslide SLHL production rate), as well as those computed with the BFPR. Figure 9 reveals that the uncorrected CRE ages, apart from the ages for the FS-03a boulder surface, do not overlap with those calculated with the BFPR within external uncertainties. Three possible explanations for the age differences appear possible: (i) previous studies have corrected the ages sufficiently for snow cover and post-depositional weathering but have ignored other factors that do affect CRE ages, such as soil cover on the boulders at the study site. (ii) Previous studies might have simply corrected the ages insufficiently for snow cover and post-depositional weathering, and other factors do not play an important role. (iii) The correction for snow cover and post-depositional weathering was insufficient, and previous studies have ignored other important factors. As the CRE ages computed with the BFPR align with independent age control, we propose that the BFPR should be used in future dating studies on the mid-elevation mountains of central Europe and that previous age estimates should be reconciled with the BFPR.

Figure 9 CRE ages of moraine boulder surfaces and external CRE age uncertainties, calculated with the Chironico landslide SLHL ¹⁰Be production rate (Claude et al., 2014) and the SLHL ¹⁰Be BFPR. Note that all ages were computed in CREp (Martin et al., 2017) with the following parameters: time-dependent Lal–Stone scaling (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008), the ERA atmosphere model (Uppala et al., 2005), and the atmospheric ¹⁰Be-based geomagnetic database of Muscheler et al. (2005).

5.3 Comparison with previously published European production rates

We hereinafter discuss the SLHL ¹⁰Be BFPR in the context of previously published European production rates (Table 7). Note that all SLHL production rates in CREp mentioned in the following paragraphs refer to time-dependent Lal–Stone scaling (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008), the ERA atmosphere model (Uppala et al., 2005), and the atmospheric ¹⁰Be-based geomagnetic database (Muscheler et al., 2005).

The internally robust “baseline” SLHL ¹⁰Be BFPR (3.64±0.11 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL), apart from the SLHL ¹⁰Be production rates at the Russenes and Grøtlandsura rock avalanches (Fenton et al., 2011), does not overlap at the 1σ level with other European SLHL ¹⁰Be production rates in CREp (Fig. 8a). Most notably, the SLHL ¹⁰Be BFPR is 11.2 % and 11.4 % lower than the ¹⁰Be Chironico landslide production rate (4.10±0.10 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL) and the canonical global production rate (4.11±0.19 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL) in CREp, respectively. However, the BFPR overlaps with the global production rate and several regional production rates at the 2σ level (Fig. 8b). The agreement at the 2σ level with the global SLHL ¹⁰Be production rate in CREp is due to the relatively large 1σ error of this production rate (0.19 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL).

Table 8 The BFPR and previously published European ¹⁰Be production rates in CREp (Martin et al., 2017). Note that all SLHL production rates refer to time-dependent Lal–Stone scaling (Nishiizumi et al., 1989; Lal, 1991; Stone, 2000; Balco et al., 2008), the ERA atmosphere model (Uppala et al., 2005), and the atmospheric ¹⁰Be-based geomagnetic database (Muscheler et al., 2005). We did not include the Rannoch Moor production rate, as this production rate is not commonly accepted (cf. Lowe et al., 2019, and references therein). See Fig. 8b for the location of the sites.

Figure 10 (a) European ¹⁰Be production rates available in CREp (Martin et al., 2017). We did not include the Rannoch Moor ¹⁰Be production rate in this compilation, as this production rate is not commonly accepted (cf. Lowe et al., 2019, and references therein). (b) Map of the calibration sites mentioned in panel (a).

5.4 Explanations for the anomalously low BFPR

Overlapping SLHL production rates (at the 1σ level) for the Chironico landslide and the Black Forest would have additionally required that (i) changes in the geomagnetic field simultaneously affected both sites; (ii) post-depositional weathering rates were similar; (iii) the snow cover on the sampling surfaces during winter was comparable; (iv) shielding by trees reduced the cosmic ray flux at both sites at a similar magnitude; and (v) other factors, such as soil cover on sampling surfaces, played a similar role during the exposure of the boulders.

Due to the comparably short distance between the Chironico landslide and the Feldsee Bog (∼170 km), changes in the geomagnetic field probably do not account for the offset. The post-depositional weathering bias at the study site was probably stronger due to the presence of soils, mosses, and shrubs on most of the sampled boulders that slowed down the production rate (cf. Dunai, 2010) and resulted in higher weathering rates (cf. Ahnert, 2009). It should be noted that the sampled boulder surfaces on the Chironico landslide “did not show any obvious signs of damage or spalling” (Claude et al., 2014, p. 277). As we observed protruding quartz veins on three moraine boulders at the study site, we argue that the higher rate of post-depositional weathering at the study site likely partly explains the difference in production rates. The snow cover bias at the study site was probably stronger than at the Chironico landslide, as the Feldsee Bog is situated in a sheltered position, at a much higher elevation, and at a higher latitude. To the best knowledge of the authors, no data on the duration and thickness of the snow cover during winter are available from a weather station close to the Chironico landslide; thus we were unable to confirm this assumption with quantitative data. The sampled boulders at both the Feldsee Bog and the Chironico landslide were situated in forested areas; thus, the forest canopy probably slowed down the effective production rates at both sites at a similar order of magnitude. As the landslide is located at both a lower latitude and elevation, the forest cover bias at the latter site should have been stronger. In contrast to the boulders at the study site, the boulders on the Chironico landslide sampled by Claude et al. (2014) were devoid of soil, mosses, and shrubs. We thus conclude that stronger snow-cover- and post-depositional-weathering bias in the study area, and the presence of soil, mosses, and shrubs on the moraine boulder surfaces at the study site, explain the relatively low BFPR.

6 Conclusions

Applying a multi-method approach to ice-marginal moraines and stratigraphically younger lake sediments in the southern Black Forest yielded consistent geochronological data and allowed us to derive the SLHL ¹⁰Be BFPR, the first ¹⁰Be production rate for the mid-elevation mountain ranges of central Europe. The BFPR turned out to be lower than most European production rates and about 11 % lower than the canonical global ¹⁰Be production rate in CREp. We attribute this offset to a strong influence of snow shielding; post-depositional weathering; and soil, moss, and shrub cover at the study site. Assessing previously suggested corrections for snow shielding and post-depositional weathering revealed that these corrections were probably too weak. Due to a lower 1σ error (0.11–0.12 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL) than the 1σ uncertainty of the canonical global production rate in CREp (0.19 atoms ¹⁰Be g⁻¹ quartz a⁻¹ at SLHL), the application of BFPR will likely allow us to calculate more accurate CRE ages for settings that resemble the study site. Among other factors, these localities must be similar in terms of weathering rates, snow cover biases, and tree shielding.

This study shows that ages generated with three independent geochronological methods, namely radiocarbon dating, IRSL dating, and tephrochronology, can successfully be interrogated in a Bayesian approach for establishing an age–depth model for lake sediments. If such a model relies on multiple independent geochronological methods, ages with higher uncertainties, such as IRSL ages, can also be included. The integration of these ages then leads to improvements to the age–depth model and increases the precision of the production rate. However, integrating a comparably imprecise geochronological method only works if this method is applied in tandem with a precise dating method, such as radiocarbon dating. If this is not the case, the low precision of the unmodelled ages results in a large uncertainty in the age model, which is reflected in a large error in the production rate.