A chronological and palaeoenvironmental re‐evaluation of two loess‐palaeosol records in the northern Harz foreland, Germany, based on innovative modelling tools

The continuing development of analytical methods for investigating sedimentary records calls for iterative re‐examination of existing data sets obtained on loess‐palaeosol sequences (LPS) as archives of palaeoenvironmental change. Here, we re‐investigate two LPS (Hecklingen, Zilly) in the northern Harz foreland, Germany, being of interest due to their proximity to the Scandinavian Ice Sheet (SIS) and the position between oceanic climatic influence further west and continental influence towards the east. First, we established new quartz OSL and polymineral IRSL chronologies. Both methods show concordant ages in the upper part of the Hecklingen profile (~20–40 ka), but in the lower part IRSL underestimates OSL ages by up to ~15 ka for the period 40–60 ka. Interpretations hence refer to the OSL data set. Second, we applied Bayesian age‐depth modelling to data sets from Hecklingen to resolve inversions in the original ages, also reducing averaged 1σ uncertainty by ~19% (OSL) and ~12% (IRSL). Modelled chronologies point out phases of increased (MIS 2, early MIS 3) and reduced (middle and late MIS 3) sedimentation, but interpretation of numerical rates is problematic because of intense erosion and slope wash particularly during MIS 3. Finally, previously obtained grain‐size data were re‐investigated by end member modelling analyses. Three fundamental grain‐size distributions (loadings) explain the measured data sets and offer information on intensity and – combined with modelled OSL ages – timing of geomorphic processes. We interpret the loadings to represent (i) primary loess accumulation, (ii) postdepositional pedogenesis and/or input of aeolian fine fractions, and (iii) input of coarse aeolian material and/or slope wash. The applied modelling tools facilitate detailed understanding of site‐formation through time, allowing us to correlate a strong peak in mean grain size at ~26–24 ka to the maximum extent of the SIS and increased influence of easterly winds.

The northern Harz foreland, Germany, is close to the former southern border of the Weichselian Scandinavian Ice Sheet (SIS) and comprises large areas covered by loess. Research on the distribution and characteristics of loess-palaeosol sequences (LPS) in the region is thus of special interest for reconstructing Pleistocene palaeoclimatic conditions in ice-marginal environments (Reinecke 2006;Lehmkuhl et al. 2016). Atmospheric dust dynamics, palaeowind velocity and change in main palaeowind directions during individual phases of a glacial cycle in the vicinity of the SIS margin can be retrieved from LPS proxy data in that area. As such, the Harz foreland also represents an essential link for comparative Late Pleistocene palaeoenvironmental reconstructions from the western and oceanic parts of Europe to the more continental areas in eastern Europe.
Furthermore, due to rain-shadow effects of the Harz Mountains the eastern part of the northern foreland only receives modern precipitation ≤500 mm a -1 (D€ oring 2004; Fig. 1). This small-scale gradient from an oceanic climate in the west to a continental one in the east offers the possibility to study the environmental effects of such variations for the past.
However, most studies north of the Harz Mountains focused on the transitional zone and its silt-and sandsized Weichselian aeolian sediments (Poser 1951;Brosche & Walther 1978;Gehrt 1994;Gehrt & Hagedorn 1996;Hilgers et al. 2001). One of the few recent investigations on LPS in the northern Harz foreland by Reinecke (2006) addressed two LPS outcrops near the villages Hecklingen and Zilly (Fig. 1). While the Zilly profile preserves only Lateglacial loess, the LPS at Hecklingen encompasses a more continuous record back to MIS 3 at least (Reinecke 2006). This is in stark contrast to observations from Saxony, Poland and western Ukraine, indicating a strong erosional phase causing large depositional gaps in most places (e.g. between~65 size, geochemical (XRF, CNS), colour and rock magnetic analyses). For chronological correlation, they used re-calculated infrared stimulated luminescence (IRSL) ages by Reinecke (2006). However, this study raised several questions on the chronology of different layers: for instance, data of the assumed MIS 3 soil complex indicate more than one phase of soil development and intercalated erosional events , implying that the period displayed in the preserved material likely exceeds 10 ka in Hecklingen (Blume et al. 2016). However, IRSL ages show an inversion between the layer below and the assumed MIS 3 complex and therefore only a short period available for pedogenesis. Additionally, Krauß et al. (2016) obtained data indicating wind direction changes towards stronger easterly winds during the local Last Glacial Maximum (lLGM). They further assumed age underestimation of Reinecke's (2006) ages by 4-6 ka (Frechen & Schirmer 2011;Schmidt et al. 2011;Zens et al. 2017) and a timing of the LGM from 26.5 to 19.0 ka (Clark et al. 2009). These discrepancies and assumptions call for further reliable chronometric control on loess deposition. Moreover, previous IRSL ages were not corrected for anomalous fading, an effect causing age underestimation (e.g. Wintle 1973;Huntley & Lamothe 2001).
When sedimentation rates as proxies of dust activity and loess depositional dynamics are to be inferred from chronostratigraphical data, age inversionsalthough to  Karger et al. (2017). B. Glacial extent during the Last Glacial Maximum (modified from Ehlers et al. 2011). Study sites Zilly and Hecklingen (Gleina for comparison) in the North European loess belt. The loess distribution is according to Lehmkuhl et al. (2018). In the Netherlands the data are based on the geological map (scale 1:600 000; Zagwijn & Van Staalduinen 1975) and for Poland (Dobrza nski et al. 1974) national soil maps were digitized. [Colour figure can be viewed at www.boreas.dk] BOREAS be expected for geomorphological and statistical reasons and a common phenomenonpose major problems in interpretation. The IRSL ages by Krauß et al. (2016) contain such an age inversion in the lower part of the Hecklingen profile. Here, Bayesian age-depth modelling represents a viable tool because (i) age inversions do not exist in modelled data and (ii) uncertainty can be reduced considerably. The method was first applied to radiocarbon dates mostly in archaeological contexts (e.g. Buck et al. 1991Buck et al. , 1996Bronk Ramsey 1995), but comparatively few methodological studies on Bayesian modelling of stratigraphically ordered luminescence samples (Rhodes et al. 2003;Millard 2004;Comb es & Philippe 2017;Zeeden et al. 2018;Lanos & Dufresne 2019) and only a small number of applied investigations exist (Noppradit et al. 2018;Peri c et al. 2019;Fenn et al. 2020). Such age-depth modelling, however, can potentially resolve the problem of a too short period of pedogenesis of the MIS 3 soil complex at Hecklingen, resulting from the age inversion in the previous IRSL data set.
Another constraint inherent to the proxy data by Krauß et al. (2016) is the selection of specific grain-size fractions taken to be indicative of source area, transport pathway, intensity and characteristics of aeolian geomorphic as well as (post)depositional processes. In general, grain size is a sensitive proxy of the listed factors (e.g . Doeglas 1968;Vandenberghe 2013;Ujv ari et al. 2016;Vandenberghe et al. 2018) and the ISO 14688 standard separation into sand, silt and clay based on the logical mathematical classification of Atterberg (1905) is often used to describe sediment characteristics including transport energies and modes (Assallay et al. 1998;Blott & Pye 2012). However, natural processes do not always follow such imposed categorizations, and the reconstruction of source area and driving transport regimes should thus avoid predefined grain-size classifications. Moreover, grain-size populations characteristic for the source material and transport process can be multi-modal or overlapping, further complicating an interpretation based on fixed grain-size classes (Dietze & Dietze 2019). Deterministic (e.g. Gan & Scholz 2017) and robust (probabilistic; e.g. Dietze et al. 2012;Paterson & Heslop 2015;Dietze & Dietze 2016) end member modelling analyses (EMMA) allow extraction of informative and unbiased subpopulations from bulk grain-size data, indicative of involved materials, transport modes and postdepositional alteration. It seems therefore worthwhile applying EMMA to the grain-size data by Krauß et al. (2016) for the Hecklingen and Zilly LPS to explore if unconstrained and unbiased extraction of grain-size subpopulations allows for new findings beyond existing interpretations on palaeowind characteristics.
The aim of this paper is to tackle the mentioned challenges and inconsistencies by producing new quartz OSL and polymineral IRSL ages for the LPS at Hecklingen and Zilly. The new chronometric datawill allow re-evaluation of part of the palaeoclimate proxy data by Krauß et al. (2016) and resolution of the discrepancies arising with their IRSL data set. We apply Bayesian agedepth modelling to circumvent age inversions and increase the precision of the luminescence chronology of the Hecklingen profile. In addition, we perform EMMA with the grain-size data previously obtained for both sites and place median and modelled grain-size data as a proxy for palaeowind velocity onto the modelled OSL chronology to estimate sedimentation rates for the Hecklingen record. Finally, we re-evaluate the previous chronometric and proxy data for the study sites and discuss the implications within the wider context of dust flux and palaeoenvironments in the North European loess belt.

Study sites and previous research
The profiles Hecklingen (latitude 51°50.451 0 N, longitude 11°31.580 0 E, altitude 106 m a.s.l.) and Zilly (51°56.280 0 N, 10°50.643 0 E, 182 m a.s.l.) are located in the northern and northeastern Harz foreland in Germany, respectively ( Fig. 1). Rain-shadow effects of the Harz Mountains influence the current climate conditions in this region, whereas in the west a stronger suboceanic influence is registered ( Fig. 1; Haase et al. 1970;D€ oring 2004;Reinecke 2006;Fabig 2007;Krauß et al. 2013Krauß et al. , 2016. Pleistocene deposits within the northern Harz foreland are mainly loose sediments and vary in thickness, composition and sedimentary origin. The loess-covered areas are part of the northern European loess belt . During the Elsterian and Saalian glacial cycles, the SIS covered the northern Harz foreland, and during the last glacial this region belonged to the periglacial area with the ice margin further north (Fig. 1). Therefore, preserved loess sediments only date back to the last glacial period and form a 40-50 km wide belt of Weichselian loess (Fig. 1). Generally, the loess cover rarely exceeds 2.5 m of thickness.
According to Krauß et al. (2016), the Hecklingen section (Fig. 2, left) comprises from top to bottom a grey tundra gley (supposedly MIS 4) at the base followed by a loess layer and a truncated soil (MIS 3). This is overlain by a presumably MIS 3 soil complex of 2 m thickness and a stony loess layer containing large amounts of gravel, followed by a~4-m-thick loess complex with four intercalated tundra gley soils (Gelic Cryosols, MIS 2). The uppermost 1.2 m represents the recent soil complex.
In the Zilly section (Fig. 2, right), a~4-m-thick loess complex with three intercalated tundra gleys (MIS 2) is preserved. The uppermost 1.2 m of the profile represents the recent soil complex. Based on proxy data, both profiles were divided into different units (units I to VI in Hecklingen and units I to III in Zilly), as displayed in Fig. 2 and explained in detail in Krauß et al. (2016).
Photographs of the cleaned profiles are compiled in Figs S1-S5.
In 2003/2004, A. Hilgers (University of Cologne) performed IRSL measurements for polymineral siltsized samples from the Hecklingen and Zilly sections. Since Reinecke (2006) published these ages without additional information concerning the dating procedure, this data set was re-evaluated by Krauß et al. (2016) using different signal integration intervals and statistical processing. Main procedural differences between the previous IRSL data set ) and the here presented new polymineral IRSL ages refer to the preheat conditions (270°C for 10 s vs. 270°C for 120 s), the IRSL reading temperature (50 vs. 125°C), an inserted 'hot-bleach' at the end of each SAR cycle (only new data set) and the methods used for dose rate assessment (c-ray spectrometry vs. thick-source a-counting and ICP-OES). A summary of the IRSL ages by Krauß et al. (2016) along with relevant dose rate data is given in Table 1.

Luminescence dating
During fieldwork, neither suitable organic or carbonbearing material for radiocarbon dating nor tephra was found that could provide numeric or relative age control for the Hecklingen and Zilly LPS. Therefore, in 2014 eight new luminescence samples from the Hecklingen profile and two new luminescence samples from the Zilly profile were taken by hammering steel tubes into the cleaned profilewalls. Doserate sampleswere takenfrom the matrix in the vicinity of the luminescence sampling holes relevant for the contribution of c-radiation. Unfortunately, it was not possible to obtain samples from identical sampling locations as used for the previous IRSL data set (Reinecke 2006;Krauß et al. 2016) due to erosion of the profile sections in between the two sampling campaigns.
The preparation of luminescence samples was conducted following common procedures for the luminescence fine grain (~4-11 µm) technique (Berger et al. 1980; refer to Data S1 in Supporting Information for details, also on instrumentation for luminescence measurements).
For OSL equivalent dose (D e ) determination on quartz, a single aliquot regeneration (SAR) protocol ; Table 2) was applied, using the signal integrated from the first 0.6 s of the decay curve minus a background averaged from the last 7.5 s. The dose points (signal integral vs. dose) were fitted with a single saturating exponential function. The D e of polymineral samples was determined using a modified IRSL SAR protocol, building on previous experience by Preusser (2003), Rother et al. (2010) and Faust et al. (2015). While employing a preheat at 270°C for 120 s prior to IRSL readout for 300 s at 90% LED power and 125°C, it additionally allows for a resting period of 20 min between preheat and IR stimulation (Table 2). This pause is intended to minimize anomalous fading (e.g. Wallinga et al. 2001;Balescu et al. 2003), which frequently causes D e underestimation of feldspar-bearing samples (Wintle 1973). Except for sample BT1543, a 'hot-bleach' of 300 s IR stimulation at 270°C was inserted at the end of each SAR cycle to reduce the signal carry-over to the next cycle. The dose points, obtained from the signal integral of the first 2 s of the IRSL decay curve corrected for a background averaged from the last 50 s, were fitted with a single saturating exponential function. For details on fading and residual dose experiments, see Data S1, Table S1 and Fig. S6.
Dose rate samples were analysed by thick-source acounting (Aitken 1985;Z€ oller & Pernicka 1989) and ICP-OES to quantify U, Th and K concentrations. Prior to radioelement analyses, sediment samples were homogenized to ensure that minute, cm-scale variations in radioactivity were averaged. Moisture estimation was mainly based on measured water contents of the samples but includes an additional absolute uncertainty of 5% to take into account moisture fluctuations during burial time. Samples from the Zilly site showed unusually low water content of 1%, certainly not reflecting a representative value for the burial period. We therefore adopted an assumed moisture content of 12.5AE5% for the two samples. The sediment cover was estimated to 50% of the present value for cosmic dose rate calculation. Following Mauz et al. (2006) and Schmidt et al. (2018), an a-value of 0.030AE0.003 (OSL) or 0.08AE0.02 (IRSL) was assumed for a-dose rate calculation. Dose rates and luminescence ages were calculated with DRAC (v1.   Krauß et al. (2016;re-calculated from Reinecke 2006). Dose rate calculations are based on radionuclide concentrations from laboratory c-ray spectrometry and equivalent dose (D e ) determination on single aliquot regeneration (SAR) measurements following  and Wallinga et al. (2000). Main differences between this previous approach and our new data set include the preheat conditions (270°C for 10 s vs. 270°C for 120 s), the IRSL reading temperature (50 vs. 125°C), an inserted 'hot-bleach' at the end of each SAR cycle (only new data set) and the methods used for dose rate assessment (c-ray spectrometry vs. thick-source a-counting and ICP-OES). Accepted in relation to measured number of aliquots.

Bayesian age-depth modelling
Age-depth modelling was done using the ADMin method (Zeeden et al. 2018), which is especially tailored for luminescence data. The age-depth model inferred from the new OSL ages was then taken to place data on an age scale through linear interpolation between modelled ages. To this aim, the function tune() in the R (R Core Team 2020) package 'astrochron' (v0.9; Meyers 2014) was used. For comparative purposes, the same procedure was performed with the grain-size data in conjunction with the 'old' IRSL data set by Krauß et al. (2016) to reveal the changes in interpretation associated with the new OSL chronology.

Unmixing grain-size data by end member modelling analyses
Details on raw grain-size data generation are given in Krauß et al. (2016). The R package 'EMMAgeo' (Dietze et al. 2012;Dietze & Dietze 2016 served to extract grain-size populations representing most closely source material, pedologic and geomorphic processes. Initially, the number of useful end members was established, thereafter end member analysis was performed without specifying weight transform limits or enforcing specific sums for scaling. However, these parameters have little influence over the result and different settings would not result in different interpretations.

Luminescence dating
The OSL signal of all samples decays to instrumental background level within <10 s (Fig. 3A). Hence, it can be assumed that the signal integral used to construct the dose-response curve is dominated by the thermally stable fast component. Dose-response curves can be well fitted with a single saturating exponential function (Fig. 3A 3B). A preheat temperature of 220°C (middle of the preheat plateau) was therefore adopted for the DRTs, which were successful with a dose reproducibility better than 2% for sample BT1534 (given dose 52.4 Gy), while a slight tendency for dose underestimation becomes apparent for sample BT1537 (7% underestimation for 196.3 Gygiven dose). The recovered dose of sample BT1541 is 11% lower than the given dose of 130.9 Gy (Table 3). With a recovery ratio of 0.98AE0.01 a given dose of 179.8 Gy could be accurately reproduced for polymineral sample BT1537 with the IRSL SAR protocol. A residual dose of 6.7AE0.4 Gy was subtracted from all IRSL D e values and considered in the DRT. IRSL fading tests indicate generally negligible to low fading rates for all measured samples. While for samples BT1535 and BT1543 all of the five aliquots yielded gvalues not significantly different from zero, results from sample BT1539 gave g-values between 0.51AE1.50 and 2.50AE1.44% per decade. Averaged g-values of samples BT1535, BT1539 and BT1543 amount to 0.62AE0.65, 1.37AE0.72 and 0.67AE0.39% per decade, respectively (n = 5 in all cases, uncertainty given as 1r standard deviation; Table S1 and Fig. S6). However, due to the additional resting period of 20 min between irradiation and measurement an additional 'fading' time delay of only~1.5 decades could be achieved. Fading correction was applied to all IRSL ages using the average of all measured g-values (0.89AE0.14% per decade) and the function calc_FadingCorr() in the R package 'Luminescence' (Kreutzer et al. 2012a;Kreutzer 2020), following the approach by Huntley & Lamothe (2001). The IRSL ages corrected for fading and residual dose are statistically indistinguishable from the uncorrected IRSL ages.
Results from radioelement analyses, moisture determination and cosmic dose rate calculation are compiled in Table 4, while D e values, total dose rates and ages are summarized in Table 5. A correlative overviewof the new Table 2. OSL and IRSL protocols used to determine the equivalent dose.
Step 1 is skipped for the first cycle (measuring the natural OSL or IRSL signal).

OSL IRSL
Step Procedure Signal Step Procedure Signal The stratigraphically lowermost sample in Hecklingen produced OSL and IRSL ages overlapping again at 1r (52.8AE3.3 and 49.1AE5.6 ka, respectively). Both OSL and IRSL ages do not significantly increase with depth below~7 m, and, unlike the OSL chronology, the IRSL data set does not contain significant age inversions. Considering recent methodological investigations (see further discussion below), OSL and IRSL ages in the lower part of the profile might suffer from underestimation.

Bayesian age-depth modelling
Calculated probability density functions (PDFs) of random uncertainties in luminescence ages for both the OSL and IRSL data sets from the Hecklingen profile are plotted in Fig. 6. They visualize the relative degree of random uncertainty that is able to explain the observed age inversions in the data set (Zeeden et al. 2018). Due to the high number of age inversions in the OSL data set, a high random uncertainty must be assigned, and the PDF curve starts to rise only above~50% random uncertainty. By contrast, the IRSL data set can be explained reasonably well with random uncertainties starting from 10-20%. A comparison of original and modelled ages along with their 1r uncertainties is shown in Fig. 4B and Table 5.
The derivation of sedimentation rates and the reinterpretation of palaeoenvironmental proxy data will, however, be based on modelled OSL ages alone. The reason is that signal fading has not been reported for OSL from non-volcanic quartz but cannot be entirely excluded as a reason for the underestimated IRSL ages from the lower part of the Hecklingen profile.

Sedimentation rates inferred from the modelled OSL ages
While it is not possible to derive trends in sedimentation rate from the original ('un-modelled') OSL or previous IRSL  ages, the modelled OSL ages suggest three intervals of rather stable mean sedimentation rates for the Hecklingen LPS (Fig. 4B). While generally the delineation of changes in sedimentation rate from age-depth plots has an inherently subjective element to it, the changes in dust accumulation rates for the Hecklingen site can be clearly constrained to the lower and upper boundary of the MIS 3 soil complex. Until~54 ka, sedimentation rates are in the order of  0.44 m ka -1 . From~54-28 ka, sedimentation rates are lower and at 0.11 m ka -1 on average. Thereafter, they increase again to a level of~0.45 m ka -1 . Note that the model does not make any initial assumptions on the sedimentation process and that it cannot differentiate between intervals including erosion and those characterized by low sedimentation rates. These reservations apply in particular to units V.I and V.II, where phases of soil erosion and subsequent deposition of soil sediment were inferred from lithology and geochemical data (Ba/ Sr and Rb/Sr ratios; Krauß et al. 2016). In comparison to LPS formed at plateau positions, the topographic situation at Hecklingen may favour distortions of the calculated sedimentation rates due to erosion and/or accumulation of material from further upslope. This should be kept in mindwhen assessing the validityof such sedimentation rates and contrasting them with values obtained at LPS that formed in different genetic settings.

End member modelling analysis (EMMA) of grain-size data
Applying the EMMAgeo algorithm results in three multi-modal grain-size populations ('loadings') as the best description of the grain-size data set of the Hecklingen profile, characterized by maxima at~25,~35 and 400 µm. The coarsest loading has additional peaks at 55 and~200 µm of roughly the same relative intensity. The EMMAgeo algorithm suggested two to seven loadings for the Zilly grain-size data set, but for reasons of consistency we decided to use three loadings here as well. Maxima of the loadings are at~15,~30 and~50 µm. Generally, the loadings obtained for the two sites show two to three additional sub-modes towards larger grainsize values (Fig. 7A, B). The end member scores as the relative contribution of each loading are shown against sampling depth for the two profiles in Fig. 7C, D. Grain size class-wise explained variance amounts to 74 and 77% for Hecklingen and Zilly, and the sample-wise explained variation averages to 89 and 83%, respectively (Fig. S7). Figures S8 and S9 show the modelled loadings and scores at both sites along with the 'classically' obtained grain-size contributions (sand, silt, clay). A more detailed description of the modelled grain-size data sets will be provided in the next section when they are referenced to age by placing them onto the Bayesian agedepth model.

Placing proxy data onto the Bayesian age model
To demonstrate the effect of Bayesian age-depth modelling on the interpretation of sedimentary proxy data for palaeoenvironmental reconstruction, we first compare the old 'linear IRSL' with the new 'Bayesian OSL' chronology for the proxy 'mean grain size' (see Krauß et al. 2016Krauß et al. , 2017 for the data set). The 'Bayesian OSL' model in combination with the mean grain size shows a shift compared to the 'linear IRSL' model ( Fig. 8) and two main changes in the timing of enhanced mean grain size. (i) The increased mean grain size documented for the MIS 3 soil complex is extended time-wise. While according to the 'linear IRSL' age model grain size reduces to values around 30-40 µm at~38-37 ka and thus long before termination of MIS 3, the new age model implies a roughly contemporaneous reduction of grain size with the end of MIS 3 at~29 ka. (ii) The threefold and strong peak in mean grain size during MIS 2 is 'squeezed' time-wise according to the new age model, approximately by a factor of 3. It now represents a comparatively short period between~26 and~24 ka. The results from EMMA as placed onto the Bayesian OSL age model are shown in Fig. 9 for the three loadings derived for the Hecklingen profile. The most characteristic feature is the drop to zero of loading 1 during MIS 2, which is considerably shortened in duration with the new age model (from~4 to~1-2 ka). Furthermore, the peak in relative contribution of loading 1 is shifted from~36-34 ka in late MIS 3 to~28-27 ka in early MIS 2 (Fig. 9A). The period of complete absence of loading 2 in the data set is extended time-wise with the new age model from a duration of~15 to~20 ka and now reaches into early MIS 2, while it ended at~35 ka (MIS 3) with the old age model. Additionally, the four peaks in relative abundance of loading 3 previously occurring between ~20 and~27 ka are condensed in time with the new model to a time window at~26-24 ka (Fig. 9B). Finally, loading 3 (Fig. 9C) represents grain sizes >55 µm and most closely resembles the mean grain-size data set shown in Fig. 8 and described earlier. The sediment-age relationship of the grain-size data plotted with thin lines in Figs 8 and 9 is based on extrapolation beyond the youngest and oldest available luminescence ages and must therefore be interpreted with care.

Luminescence methodology
Contrasting the new OSL against the IRSL data set generated on identical samples, it appears that in the upper part of the profile (<6 m) IRSL ages are systematically underestimated by~1-5 ka as compared to OSL ages, although four out of the five samples in this section (BT1538-BT1543) overlap within 1r uncertainty. If this observation is related to fading that was not adequately corrected for, remains open at this point. Slight inaccuracies in the chosen dose rate parameters (e.g. a-value) could likewise cause an offset between OSL and IRSL ages.
A somewhat different picture emerges for the lower part of the profile (>6 m), encompassing samples BT1535-BT1537. Here, two out of the three samples yield OSL and IRSL ages significantly different from each other (at 1r), while IRSL ages underestimate OSL ages by~4-15 ka. Both OSL and IRSL ages do not significantly increase with depth here, which could be explained either by high sedimentation rates at the base and below the soil complex of units V.I and V.II or by onset of dose saturation for both OSL and IRSL signals. Although OSL D e values are well below the 2D 0 threshold, the commonly accepted criterion for reliable SAR dose estimation, age underestimation of OSL fine grain samples cannot be ruled out. Such behaviour has frequently been reported from loess records and other environments across Eurasia (e.g. . Consequently, these OSL ages could be regarded as minimum ages. Additionally, for reliable estimation of D 0 , dose-response curves should be constructed up to full saturation (Timar-Gabor & Wintle 2013). For IRSL data, signal saturation is not considered a major problem with D 0 in the range 300-350 Gy.

Comparison with the previous IRSL chronology
The new OSL and IRSL chronologies of the upper part (<5 m) of the Hecklingen profile conform within uncertainties with the re-calculated IRSL ages by Krauß et al.  (2016) (Tables 1, 5), with the new ages generally appearing to be slightly older. However, since the sampling positions of the two sets of samples (Fig. 5) are not identical, a direct numerical comparison including significance tests is not possible. In the lower part of the profile, the upper IRSL ages (BT1536 and C-L1309; unit VI.I) assigned to MIS 4 appear to be underestimated, while the lowermost IRSL age by Krauß et al. (2016) roughly corresponds to stratigraphical expectations. The OSL age of sample BT1535 (52.8AE3.3 ka; unit VI.II) might be underestimated, as indicated by the DRTresults of an adjacent sample (Table 3). A systematic underestimation of the Krauß et al. (2016) IRSL data cannot be confirmed, because some of the ages seem plausible with regard to the new OSL chronology (e.g. C-L1310) or show even an older age (C-L1308).
Comparing the dose rate data from Reinecke (2006) and Krauß et al. (2016) to those of the new OSL and IRSL results, it becomes clear that determined radioelement concentrations and thus the total dose rate for the more recently taken samples are systematically higher, especially for the Hecklingen profile. In particular, the K and U values measured by a-counting and ICP-OES significantly exceed those determined by c-ray spectrometry (K by~16% and U by~50%). Therefore, dose rates calculated for the new IRSL samples are on average~25% greater than those of the previous samples. Despite intercalated zones of pedogenesis potentially enriched in radioactivity, loess is an overall homogeneous sediment, and within analytical uncertainties we could not establish a correlation between the presence of a palaeosol (tundra gley or the MIS 3 soil complex) at a sampling spot and the measured radioelement concentration. Since some of the samples were taken in approximately comparable stratigraphical positions, a systematic difference of the dose rate values of both data sets is hence unexpected and points to inaccuracies in some of the analytical procedures for dose rate determination. Consequently, some of the mismatch between old and new luminescence chronologies could be attributed to issues in dose rate assessment.
OSL ages of samples BT1536 (59.9AE3.8 ka) and BT1537 (57.6AE3.6 ka) support the previous assumption ) that the truncated transition zone of a former MIS 3 soil (unit V.II) and the loess layer below (unit VI.I) have a similar age and therefore the former soil must have developed in the loess layer. Furthermore, the age of sample BT1538 (40.6AE2.7 ka) suggests indeed that the MIS 3 soil material (unit V.I) was re-located and pedogenically overprinted after re-deposition as previously proposed.
The new OSL ages at Zilly of~15 ka are indistinguishable (within 1r uncertainty) from the previous IRSL ages, irrespective of the sampling position. This seems to rule out substantial anomalous fading in the latter ages, being consistent with the observation that the newly generated IRSL ages show significant age underestimation only beyond 40 ka.
Identical OSL ages of samples BT1534 and BT1533 support the results of geochemical and grain-size proxies discussed in Krauß et al. (2016), indicating high accumulation rates of predominantly silt-sized (~80%), unaltered material during MIS 2.

Bayesian age-depth modelling
Generally, Bayesian age-depth models use the initial age information together with the stratigraphical principle of superposition commonly leading to a better understanding of the true age and reduced uncertainty (Bronk Ramsey 2009;Zeeden et al. 2018). Outliers can be problematic for such modelling approaches and affect ages and credible intervals for the entire modelled sequence. Some methods allow the setting of outlier   (Parnell 2018), which we did not include here. The 1r credible interval assumes that all ages are true (and no outliers) and that the uncertainty derived from luminescence dating is correct. Most models assume uncertainty to be fully random, which allows a maximum in uncertainty reduction. The ADMin model applied here (Zeeden et al. 2018) assumes the reported uncertainty to be correct, but deliberately does not make an initial assumption about which proportions of the reported uncertainty are systematic or random. Rather it derives this information from the data set itself in a probabilistic manner. This is expected to lead to less increase in accuracy, which, however, is more realistic for luminescence ages. Yet, an unsolved issue is the treatment of data points that might be subject to a systematic uncertainty larger than those for the remaining data points, or which are characterized by an asymmetric systematic uncertainty. For instance, this becomes relevant for luminescence ages probably underestimated as indicated by diagnostic laboratory tests, or for luminescence ages derived from poorly bleached sediment such as soil that has undergone gelifluction. For those cases, a statistically valid procedure is still pending, but beyond the scope of this contribution.
For the Hecklingen profile, the modelled ages (both OSL and IRSL) show a consistent stratigraphy with increasing age corresponding to increasing depth. For most samples, original and modelled ages show overlapping 1r uncertainty, and several samples hardly change in age. The modelled uppermost and lowermost samples are not obviously younger or older than expected, which is a good model outcome, as outermost samples can indicate modelling problems (Bronk Ramsey 2000;Steier & Rom 2000). In our data set, the age-depth modelling reduces 1r uncertainty by 19 and 12% on average for OSL and IRSL ages, respectively, ranging from 1 to 37% and from À3 to 18% for individual samples. This corresponds to results by Zeeden et al. (2018), who found uncertainty reduction of 21-30% for a loess data   set with more samples. The youngest IRSL sample showed an increase in uncertainty following age-depth modelling, which might be an issue of boundary effects in the modelling process. Finally, it is possible that the age inversions in the OSL and IRSL data sets are not of pure statistical origin but related to geomorphic processes. The location of the two LPS close to the former ice margin favours periglacial processes such as gelifluction, leading to spatio-temporal variations in bleaching conditions and therefore potentially to age inversions. Likewise, temporal changes in dose rate induced by varying moisture/ice content in the sediment or secular disequilibria in the U decay series could also produce reversed ages (Degering & Degering 2020). The change in grain-size composition with a higher proportion of loading 1 at <4.70 m might indicate a shift in source area with material of different luminescence properties, hence potentially triggering age inversions as well.
Age-depth modelling and grain-size unmixing for an improved interpretation of palaeoenvironmental proxy data The end member modelling results (Fig. 7) demonstrate how complex LPS are in their function of documenting palaeoenvironmental developments over time. The loadings alone (Fig. 7A, B) would not allow a detailed evaluation of the proxy data since they represent the entire profile without distinguishing between different phases. This is especially problematic when LPS include different phases of soil development. Contextual knowledge about the investigated systems is therefore vital to interpret the results correctly (Dietze & Dietze 2019). Thus, the interpretation of the sequences attempted here additionally takes into account the results of Krauß et al. (2016).
In Hecklingen, loading 1 (Fig. 7A) very likely includes postdepositional soil formation and aeolian deposition of finer fractions without distinguishing between the two with depth (Fig. 7C). Vandenberghe et al. (2018) state that clay minerals might reach similar grain sizes in their pristine shape as dust sediments and hence recommend including all available site information (cf. Vandenberghe 2013; Schulte & Lehmkuhl 2018). Loading 3 (Fig. 7A) seems to include portions of aeolian but also slope wash deposits of coarser material. However, loading 2 (Fig. 7A) shows a single peak at 35 µm, lying within the range of European loess. Vandenberghe et al. (2018) and Dietze & Dietze (2019) suggest that it is possible to unmix grain-size composition and identify the primary windblown signal after alteration and potentially filter out the characteristic imprint of reworking processes. Thus, we interpret loading 2 as the primary loess end member representing dominant westerly wind directions (e.g. Renssen et al. 2007), since the two other loadings might display different geomorphic and pedogenic processes. The bottom of the sequence is mainly composed of loadings 1 and 2 (Fig. 7C). Following the penultimate argumentation, results display a phase of pronounced deposition of aeolian material at the bottom, which decreases continuously towards the top in favour of finer material, indicating successive soil development in the deposited loess material. The mean grain size combined with the Bayesian OSL age-depth model (Fig. 8) shows a substantial improvement in reasonable timing of the loess accumulation phase and first pedogenesis during MIS 3. Figure 9 gives a more detailed picture of the continuously lowering contribution of loading 2 towards pedogenesis, matching the timing of MIS 3 interstadials and soil development well.
Between 7 and 6.5 m depth, a sudden switch in the contribution of different loadings indicates an erosional event that might have truncated the assumed soil (Fig. 7C). Up to~5 m depth loadings 1 and 3 contribute alternating but similar scores. This suggests shallow slope wash and re-deposition of soil sediment but also pedogenic overprinting by formation of finer grained secondary clay minerals. During weathering processes, iron oxides are released and material is coloured brownish to reddish (brunification), the primary minerals are transformed, and secondary clay minerals form as well (loamification, e.g. Blume et al. 2016). Geochemical and colour data in Krauß et al. (2016) support the assumption of pedogenically overprinted soil sediment. This phase is bordered by a 100% contribution of loading 1 at~4.7 m, indicating either erosion followed by sedimentation or a rapid burial of the soil complex. Since no intact A-horizon was left and almost all palaeosols are more or less truncated, we propose the EMMA results to reflect erosion followed by loess accumulation. This agrees with the Bayesian OSL age model putting the end of the MIS 3 soil complex at~29 ka (Figs 8,9). The varying proportions between loadings 1 and 2 at~4.5-3.5 m might represent tundra gleys (E0, E1, E2, Fig. 5).
At~3.5 m loadings 2 and 3 solely make up the grainsize composition (Fig. 7C). Considering that loading 1 is completely missing, loading 3 might in this case represent rather aeolian deposition of coarser material than reworking of material. Further, the uniqueness of the layer (~3.4-2.8 m) also argues for a change of source contribution and therefore wind direction. Geochemical and colour data of Krauß et al. (2016) support this assumption and indicate a shift of wind direction towards stronger and more frequent easterly winds during a postulated time frame confined more or less to the LGM. The strong influence of the SIS and its last maximum extent should receive more attention as a major cause for these shifts than the overall influence of the Greenland Ice Sheet. The timing and the number of ice advances of the SIS are still a matter of debate and under investigation (L€ uthgens et al. 2020).
The Bayesian OSL age model puts this accumulation phase in a narrower time window of 26-24 ka as compared to the linear IRSL age model (Figs 8,9). Considering the findings of Hughes et al. (2016) andB€ ose (2011), it is thus possible that this accumulation phase in MIS 2 is connected to the most southerly extent of the SIS and its stronger influence on atmospheric circulations in the region, causing easterly winds with strong dust storms. However, the most pronounced east wind activity in the Eifel (western Germany) appears to fall into the period 23-20 ka with a short interruption near 21 ka (R€ omer et al. 2016). As the time frame preserved is different at Zilly (Fig. 5), loadings show a different composition than at Hecklingen (Fig. 7A, C). However, loading 2 could be interpreted as the primary loess end member displaying regular westerly wind induced dust storms. Loadings 1 and 3 also seem to include aeolian components but might as well show a mixture of process signals (Fig. 7B). At the bottom of the sequence, EMMA suggests a different input source through high relative contributions by loading 3 (Fig. 7D), indicating a shift towards a coarser aeolian component. However, at Hecklingen the scores of loadings 2 and 3 behave similarly during the most southerly extent of the SIS, while they behave contrarily at the bottom of the Zilly profile (Fig. 7D). This suggests that at Zilly a short phase during the deposition of the last cover loess and post-LGM ice sheet retreat initiated a switch of source material assuming that loading 2 represents dominantly westerly winds. At~3-2.5 m the sediment composition is almost solely explained by loading 2; loading 1 takes the lead from~1.2 m towards the top. Above 2.2 m, loading 2 represents a phase of windblown sediments, in which the Holocene soil successively developed. As loading 3 is missing for depths <1 m, the modern soil was not ploughed deeply, unlike in Hecklingen.

Comparison with neighbouring loess areas
The following comparison of our results from the northern Harz foreland with those from neighbouring loess areas is restricted to well-studied LPS southeast and east of our study area. Further to the west or southwest well-studied LPS appear only quite far away along the Lower Rhine Embayment with different actual climate and much greater distance from the LGM ice advance (see Fig. 1). Thick loess occurrences in Saxony and Lower Silesia (western Poland), however, are also situated at the periphery of the Central European Lowland towards the Upland (Ore Mountains, Sudetes) and affected by the drier and more continental climate in the rain-shadow of the Harz Mountains and the northern ridges of the Bohemian Massif.
The 'Gleinaer Bodenkomplex' (Gleina Pedocomplex) was supposed to represent cumulative interstadial pedogenesis from the upper part of the Middle Weichselian (Lieberoth 1963). More recent studies (Meszner et al. 2013(Meszner et al. , 2014 on Saxonian LPS and a composite profile demonstrated, however, that based on OSL dating this complex marks an important gap in loess sedimentation between 60 and 30 ka in central Europe. Therefore, this sequence is designated as the 'Gleina Complex' (without the word 'Boden', i.e. 'soil'). As far as loess-like material of the Gleina Complex (III) is preserved at all, it consists of brownish and gleyic remnants of interstadial soils and is characterized by repeated re-deposition. The Gleina Complex is estimated to be some 30 ka old, following a phase with strong erosion. The thick loess above the Gleina Complex (~30-18 ka) looks foliated in its lower part (IIb) and like homogenous loess in its upper part (IIa). So far, a strong similarity between loess of the Harz foreland and of Saxony is apparent.
The composite profile of Saxonian LPS reaches further back in time than the Hecklingen profile, and a comparison of both below the Gleina Complex is questionable with respect to the present database. The drier position of Hecklingen compared to Gleina may explain the preservation of the Hecklingen soil complex despite some evidence of re-deposition (coarser embedded material, loading 3 in Figs 7C ice advance, this erosional unconformity is twofold and~34 ka long.

Conclusions
We established new OSL and IRSL chronologies for two LPS (Hecklingen, Zilly) in the northern Harz foreland, which are of particular interest for palaeoenvironmental research due to their proximity to the maximum extent of the SIS. Application of a Bayesian age-depth model was successful in resolving age inversions in the Hecklingen OSL and IRSL records and produced a consistent agedepth relationship. Thereby, we significantly increased precision of individual ages (by~19 and~12% on average for OSL and IRSL). Although fading-corrected, IRSL underestimates OSL in the lower part of the Hecklingen profile, and further inferences were based on the OSL data set. While the modelled age-depth relationship technically facilitates the calculation of sedimentation rates, its validity has to be carefully assessed owing to the non-continuous nature of the LPS at Hecklingen, featuring erosive phases and re-location through slope wash.
The effect of the new Bayesian OSL age model for the Hecklingen site as compared to previous IRSL ages was demonstrated with the grain-size data set compiled in a precursor study (Krauß etal.2016):inputofcoarseaeolian material extends over a much longer period in MIS 3 until its termination at~29 ka and the sharp peakof mean grain size in MIS 2 can be constrained to a much shorter interval (~26-24 ka). This distinct accumulation of relatively coarse aeolian sediments in MIS 2 can thus be more confidently correlated to the maximum extent of the SIS and strong easterly winds prevailing during that period.
Finally, EMMA of the grain-size data by Krauß et al. (2016) allowed us to identify three fundamental end members for the Hecklingen and Zilly sections that can be interpreted to represent characteristic geomorphic processes such as the deposition of primary loess, slope wash of relatively coarse material or the breakdown of grain sizes due to pedogenesis. Combining the EMMA results with the new Bayesian OSL age model created a high-resolution record of environmental dynamics tracing changes in dominating geomorphic processes independent of imposed grain-size classifications.
We conclude that the new chronometric data along with the Bayesian age-depth and grain-size end member modelling offered a series of new insights into palaeoenvironmental dynamics in the study area. As such, these modelling approaches might represent an incentive to reevaluate existing LPS or other sedimentary data sets to fully unfold the potential of these records.
Acknowledgements. -The work undertaken in this study was financed by the German Science Foundation within the project 'Harsh periglacial climate at the MIS3/MIS2 transition reflected in Central European Loess-Paleosol sequences' (project-ID ZO51/39-1) and is affiliated to the CRC 806 'Our Way to Europe', subproject D1, also generously funded by the German Science Foundation (DFG, project-ID 57444011 -SFB 806). We thank Manfred Fischer for valuable help with routine work in the Bayreuth luminescence laboratory, as well as Prof. Frank Preusser, Prof. Regina DeWitt, two anonymous reviewers and the Editor for helpful and critical comments that helped substantially to improve the quality of this paper. We declare that we have no conflict of interest.
Author contributions. -LK, CZ and FL carried out fieldwork and took the samples. CS was responsible for luminescence dating and related analyses. CZ, LK and CS did the Bayesian age-depth modelling, EMMA and the interpretation of results. FL, LK and CZ prepared part of the figures. LZ and FL helped to discuss the data and to edit the manuscript. CS wrote the manuscript and prepared part of the figures with significant input from all co-authors.
Data availability statement. -The data that support the findings of this study are available from the corresponding author upon reasonable request.         Table 2 in the main text) for the samples BT1535, BT1539 and BT1543.    Table S1. Summary of fading rates measured for samples BT1535, BT1539 and BT1543 with the IRSL protocol outlined in Table 2. The shown g-values are normalized to a delay time of 2 days.
Data S1. Details of luminescence dating procedures.