Disentangling Sedimentary Pathways for the Pleniglacial Lower Danube Loess Based on Geochemical Signatures

Introduction

Since the earliest days of loess research, the origin and source of aeolian sediments such as loess have been a matter of debate. Richthofen (1882), a pioneer in loess research, named non-vegetated deserts and desiccated lake beds as potential source areas for loess in China. Throughout the decades, the definition of loess as well as its sources and origin have been subjects of debate in numerous studies (see Smalley et al., 2001; Marković et al., 2016; Sprafke and Obreht, 2016 and references therein). However, most studies agree that the formation of dust (i.e., silt-sized particles) essentially involves a series of sedimentary processes prior to aeolian transport and loess formation (Smalley and Smalley, 1983; Pécsi, 1990; Pye, 1995; Muhs and Bettis, 2003; Badura et al., 2013; Sprafke and Obreht, 2016). The knowledge of the sedimentary pathways, thus the transport chains from the initial production of silt-sized particles to its final deposition and loess formation, is crucial to understand past environmental and climatic systems in varying geographical settings (Muhs, 2018). The reconstruction of atmospheric circulation patterns and climatic gradients is crucial to understand past environmental constraints on the spread of early modern humans (Obreht et al., 2017; Staubwasser et al., 2018), especially in areas which are believed to be important trajectories for their migration into Europe (Anghelinu and Niţă, 2014; Chu, 2018; Fitzsimmons et al., 2020).

The steps involved in dust production can be generalised and summarised in the pathways of desert and glacial loess genesis (Tsoar and Pye, 1987; Wright, 2001; Muhs, 2013; Muhs et al., 2014; Lancaster, 2020 and references therein). The main differences of these two pathways are the processes of dust formation, thus the production of silt-sized particles. Silt-sized particles of desert loess are mainly generated by physical weathering in high altitude areas and desert basins. This mode of genesis is constrained to the margins of drylands in low‐ to mid‐latitudes as well as extreme continental regions (Wright, 2001; Lancaster, 2020). For glacial loess, however, grinding by mountain glaciers and continental ice sheets (CIS) as well as frost shattering of periglacial mountainous areas are the main dust formation processes (Smalley and Smalley, 1983; Smalley et al., 2014). Furthermore, various authors highlighted the general significance of rivers as agents of reworking and storage of silt-sized particles, respectively (Smalley and Leach, 1978; Smalley et al., 2009; Badura et al., 2013; Lehmkuhl et al., 2018a, 2018b). The distinction between dust formation by CIS and processes in the high mountain areas was recently postulated by Li et al. (2020), who differentiated between the continental glacier-river transport (CR) mode and the mountain provenance-river transport (MR) mode of loess genesis.

The two glacial modes of loess genesis can be assumed as the main sedimentary pathways of European loess (Figure 1; Li et al., 2020). Loess deposits in northern Central Europe, roughly extending over southern England, France, Germany, Poland, Belarus, and Ukraine, are believed to partially originate from glacially ground material derived from the Fennoscandian and British-Irish ice sheets (Antoine et al., 2009; Skurzyński et al., 2019, 2020; Stevens et al., 2020). In Central and Southeastern Europe, silt particles forming loess deposits were subjected to a multi-step transport, essentially involving river transport of silt particles generated by intensive physical weathering and/or glacial activity in the adjacent mountain ranges (Smalley and Leach, 1978; Buggle et al., 2008; Smalley et al., 2009; Újvári et al., 2012).

FIGURE 1. (A) Distribution of aeolian sediments (modified according to Haase et al. (2007). Continental glacier provenance-river transport (CR) and mountain provenance-river transport (MR) modes of loess genesis (modified according to Li et al., 2020) for Europe are displayed. Danube (Da.), Dniester (Dn.) and Dnieper (Dp.) rivers are shown for orientation. Extent of the dry continental shelfs during the last glacial maximum (LGM) is modified according to Willmes (2015). (B) Location of the Lower Danube Basin in southeast Europe. Black triangles show investigated and reference loess-palaeosol sequences (LPS): 1: Urluia; 2: Rasova; 3: Vlasca; 4: Sageata; 5: Balta Alba Kurgan; 6: Langenlois; 7: Ságvár; 8: Bodrogkeresztúr; 9: Orlovat; 10: Petrovaradin; 11: Irig; 12: Kurortne. Elevations >300 m above sea level (a.s.l.) are shown as the upper topographic boundary for loess and archaeological open-air findings in Southeast Europe (Hauck et al., 2018). (C): Detailed view of the Lower Danube Basin with investigated LPS.

However, global models of loess genesis simplify transport processes and atmospheric circulation patterns, which hampers their (supra-) regional application. Regions can be influenced by several dust sources, which correspond to both modes of loess genesis. In Eastern Europe for example, the boundary between the two modes is roughly parallel to the Eastern Carpathians (Smalley et al., 2009; Li et al., 2020; Pańczyk et al., 2020; Skurzyński et al., 2020). Hence, the Lower Danube Basin (LDB) is located near the boundary between these two modes and thus it represents a transitional depositional zone, where a mixing of several potential source areas occurs (Buggle et al., 2008).

Therefore, the aim of this study is to evaluate potential source areas for the loess of the LDB and to refine our understanding of the modes of loess genesis in the region. We hypothesise five possible regional sedimentary pathways that may account for loess formation in the LDB and asses their importance alongside results of high-resolution geochemical analyses of five loess-palaeosol sequences (LPS) distributed throughout the basin. First, geochemical compositions are compared in order to evaluate signatures of each loess-palaeosol sequence that may reflect commonalities and differences in terminal dust sources. The continuous sampling in equidistant intervals allows us to reliably assess shifts in sediment provenance over time. Second, the elemental compositions and respective ratios are compared to loess records from reference regions upstream the Danube in Lower Austria, Central and Northeastern Hungary, and the Vojvodina region in Northern Serbia as well as loess deposits from the shore of the Black Sea in Ukraine. After establishing commonalities and differences of the different loess regions, the results are compared with findings from published provenance studies from the reference regions.

Regional Setting and Study Sites

The Lower Danube Basin

The LDB is the last depositional segment of the Danube River catchment, spanning between the Iron Gates and the Danube Delta (Figures 1, 2). The Lower Danube drainage basin sensu stricto covers large parts of the Transylvanian Basin, most of the Romanian Carpathians as well as the north facing slopes of the Balkan Mountains (Figure 2). The LDB sensuJipa (2014), who defined it as a loess sedimentary basin, only covers the lowlands of the Danube Plain (Romanian/Wallachian Plain and Bulgarian Plain, hereafter combined as Danube Plain) as well as the plateau of the Dobrogea. The basin is bounded by the Carpathians to the north and west, the Balkans to the west and south, and opens toward the Black Sea in the east. The Moldavian Plateau adjoins in the northeast (Figure 2). The basin is mainly covered by detrital Late Neogene to Pleistocene/Holocene deposits of predominantly marine, fluvial, and aeolian origin (Krézsek and Olariu, 2021). The aeolian sediments consist of loess, sandy loess, and Pleistocene to Holocene sand dunes, particularly in the western part of the basin, and essentially loess and loess derivates in the central and eastern part (Lehmkuhl et al., 2021).

FIGURE 2. (A) Perspective view of the LDB, looking from the southeast. Extents are indicated in Figure 1. Shaded areas indicate landscape units within the LDB: I: Wallachian Plain; II: Bulgarian Plain; III: Dobrogea; IV: Moldavian Plateau; V: Subcarpathians; VI: Pryedbalkans; VII: Transylvanian Basin; VIII: Carpathians; IX: Balkans; X: Danube Delta. Hypothesised sedimentary pathways are shown as black arrows (DA: Danube Alluvium pathway; CB: Carpathian Bending pathway; BS: Black Sea pathway; GF: Glaciofluvial pathway). Note that the scale bar is only accurate in the front part of the figure, due to the perspective view (B) Stratigraphic sketches for the investigated sections with correlation scheme.

The loess deposits in the LDB are among the thickest in Europe and laterally continuous over large areas (Haase et al., 2007; Fitzsimmons et al., 2012; Jipa, 2014; Marković et al., 2015). The landscape of the basin is diverse and allows for its division into seven subunits. The Carpathian foothills are intersected by a multitude of Danube tributaries such as Jiu, Olt, Argeş, Ialomiţa, and Buzău (Figure 1C). The Bulgarian Plain is characterized by a hilly landscape with short, steep tributaries discharging into the Danube that incises into mainly Mesozoic sediments underlying the Late Neogene to Quaternary mantles of the plain (Figure 2; Krézsek and Olariu, 2021). The Wallachian Plain is dominated by a flat topography, gently inclined toward the south. The left-side tributaries of the Danube, draining the Carpathians and their foothills, are sediment laden and show braided river characteristics. Therefore, they intersect the widespread Late Neogene to Quaternary deposits in broad, alluvial valleys, leaving vertical bluffs which often expose thick loess sequences. The topography of the Dobrogea plateau shows many ephemeral valleys, draining the area towards the Danube and the Black Sea. The Moldavian Plateau, which adjoins to the northeast, is dissected by the broad valleys of the Siret and the Prut, large tributaries of the Danube which are draining the Eastern Carpathian flysch. Additionally, the plateau is deeply dissected by a dense dendritic network of tributary systems.

Study Sites

The investigated loess sections are located along a south-north transect of the eastern part of the LDB in different geomorphological settings (Figure 2).

1) The Urluia (URL) LPS is located in an abandoned limestone quarry on the limestone plateau of the Dobrogea (Fitzsimmons et al., 2013; Fitzsimmons and Hambach, 2014; Obreht et al., 2017). The outcrop preserves five loess-palaeosol couplets reaching to the Middle Pleistocene. In this study, we focus on the uppermost 8.5 m, spanning from Marine Isotope Stage (MIS) 3 to 1. The Campanian Ignimbrite (CI) tephra, which is an important isochron in the area (Fitzsimmons et al., 2013; Veres et al., 2013), can be found in a stratigraphic depth of approximately 6 m at the studied profile. At the sampling site, tephra particles are syn-depositionally mixed into the loess and are exposed in an approx. 10 cm thick characteristic bed (Obreht et al., 2017). The tephra layer is dated in its source region of the Phlegrean Fields in Italy to an age of roughly 39–40 ka (De Vivo et al., 2001; Giaccio et al., 2017). The URL LPS was additionally dated by luminescence dating. The results were incorporated into a correlative age model (Obreht et al., 2017). This age model dates the uppermost 8.5 m of the section to a time span reaching back to ca. 55 ka.

2) The Rasova (RAS) sequence is exposed in a small valley on the right bank of the Dunarea Veche branch of the Danube bordering the Dobrogea to the west. The section is 7.2 m thick and overlies weathered Cretaceous limestone (Zeeden et al., 2018). The CI tephra is found at a depth of around 6 m and was confirmed by geochemical evidence and by luminescence dating (Anechitei-Deacu et al., 2014). The tephra is preserved as a thick bed of pure volcanic ash, with glass shards also incorporated into the overlying loess. Furthermore, the chronology of RAS relies on luminescence dating and a correlative age model (Zeeden et al., 2018) linked to the NGRIP δ¹⁸O record (North Greenland Ice Core Project members, 2004). The resulting age model dates the base of the sequence to ca. 47 ka.

3) The Vlasca (VLA) LPS is exposed at the left bank of the Braţul Borcea branch of the Danube River bordering the Wallachian Plain toward the east (Obreht et al., 2017). The loess bluff has a length of more than 2 km and a height of approx. 30–40 m above the current Danube level. The investigated part of the outcrop, covering MIS 3 to 1, is 12.5 m thick and the CI tephra is situated at a depth of 9.1 m. The tephra is preserved as a lateral continuous bed, several centimetres in thickness. The section was dated by a correlative approach to URL and Greenland ice core data (Obreht et al., 2017), which led to an age model with a maximum age of ca. 44 ka in a depth of 10.72 m. Furthermore, the underlying 1.8 m of the here investigated section were correlated to the published proxy data from URL covering the entire MIS 3 (Supplementary Material 3).

4) The section Sageata (SAG) is situated on the left bank of the Buzău River and is 5.5 m thick. The LPS consists of sandy loess sediments and overlies river pond loams, which developed in the floodplain of the Buzău River. A more than 50 cm thick layer of pure ash, assigned to the CI tephra, is intercalated in the sequence. The section was not directly dated yet but correlated to the other LPS records using the CI tephra and the Holocene topsoil as tie points (Supplementary Material 3).

5) The Balta Alba Kurgan (BAK) site is situated 1.5 km west of the alkaline Balta Alba lake in the northern Danube Plain. The section is exposed in a road cut opening the edge of a loess plateau. Within the upper 2 m of the section, remains of a Neolithic to Bronze Age tumulus can be found. Below the tumulus, the sequence consists of loess with several intercalated interstadial palaeosols and weakly pedogenically overprinted horizons. The CI tephra lies at a depth of 8.06 m and is preserved as a horizontal bed of patchy tephra lenses. The investigated section is 10.9 m thick. The chronological framework for this section is based on an integrative age model, using luminescence dating, radiocarbon dating as well as (palaeo-) magnetic correlation (Scheidt et al., 2021). The resulting model shows a maximum age of ca. 53 ka at a depth of 9.8 m.

All investigated LPS were correlated using proxy data and their respective age models revealing a set of records covering the timeframe from MIS 3 to MIS 1.

Possible Sedimentary Pathways of the Lower Danube Loess

The loess deposits of the LDB are located in the transitional zone between the East European loess zone and the Caspian-Black Sea loess zone (Stephens et al., 2002). The provenance of the Danube alluvium (DA) in the Lower Danube Basin is diverse, including sediments derived from the upper reaches of the river in Central Europe but also more locally derived material carried by its main tributaries within the basin. The DA is regarded as the main source area for loess deposits within the LDB (Buggle et al., 2008; Jipa, 2014). The Carpathian Bending (CB), narrowly defined here as the flysch area, drained by the Buzău and Ialomiţa rivers consists mainly of Cretaceous and Neogene flysch deposits, as well as unconsolidated but uplifted Pliocene to Quaternary marine to brackish sediments evolving into lacustrine deposits and locally interfingered with fluvial and alluvial fan deposits (e.g, Krézsek and Olariu, 2021 and references therein). The sediment yields of these rivers are among the largest in Europe and their carrying capacities are strongly reduced when entering the basin. Their river flow regime changes from an intra-mountainous, hence erosive, to a depositional fluvial system. These rivers built up large alluvial fan systems and even nowadays show braided river characteristics (Râdoane et al., 2003; Floroiu, 2011). These sedimentary cascades are also influenced by complicated local tectonics (Matenco et al., 2016; Krézsek and Olariu, 2021). The strong tectonic subsidence in the Lower Danube area, that comprises several km thick Quaternary marine to brackish, lacustrine, and fluvial deposits (Matenco and Andriessen, 2013; Matenco et al., 2016; Krézsek and Olariu, 2021) coupled with uplifting in the Carpathian Bending area (Necea et al., 2013; Krézsek and Olariu, 2021) considerably influence the riverine morphocharacteristics and sediment yield. Thus, the exposed large valley plains act as deflation areas within and in the foreland of the Carpathian Bending area (Jipa, 2014). The third pathway, which similarly to the DA and CB pathways corresponding to the MR mode (Li et al., 2020), is the drainage of the Eastern Carpathians (EC). This mode is mainly fed by eroded Cretaceous and Neogene flysch sediments, as well as widespread unconsolidated fluvial, aeolian, and shallow marine deposits of Late Miocene to Quaternary age covering the Eastern European Platform. These areas are drained by the Danube River system with its tributaries Siret and Prut, as well as the Dniester River that discharges into the Black Sea. This sediment pathway was identified as an important dust source for loess deposits in Poland and (southwestern) Ukraine (Nawrocki et al., 2018; Pańczyk et al., 2020). The glaciofluvial outwash plains of Central and Northern Ukraine pose a major dust source for Ukrainian loess-palaeosol sequences (Nawrocki et al., 2006, 2019; Buggle et al., 2008; Pańczyk et al., 2020). This pathway, hereafter called glaciofluvial (GF) pathway, can be related to the CR mode (Li et al., 2020). Another potential sediment pathway is the dust entrainment from exposed sediments on the shelf of the Black Sea (BS) during low sea level stands. This pathway cannot be directly linked to any loess genesis mode sensuLi et al. (2020), since the catchment of the Black Sea between southern Ukraine and Bulgaria covers riverine input from several mode domains, including the DA. For example, the sediment yield of the Danube corresponds to the MR mode, whereas the sediment yield of the Bug, Dnieper, and Don is derived from the CR mode. Since the BS pathway is assumed to be an important dust source area for loess deposits in Bulgaria (Avramov et al., 2006; Jordanova et al., 2007), it is considered as a potential pathway for the LDB loess.

Sampling and Laboratory Analyses

Methodological Background of Loess Provenance Studies

The range of methodological approaches used in loess provenance research is broad. Conceptual research assessed dust sources based on the geographical and geomorphological situation of loess deposits (Smalley and Leach, 1978; Leger, 1990; Smalley et al., 2009; Badura et al., 2013; Lehmkuhl et al., 2016, 2018b) or prevailing atmospheric regimes favoring dust transport and deposition (Antoine et al., 2001, 2009; Vandenberghe et al., 2006). Spatial variability in grain size distributions can be used as a proxy for the distance from the source area (Frazee et al., 1970; Jipa, 2014; Schaetzl et al., 2018). More complex statistical analyses of grain size data, e.g. based on end member models, have been used to unmix grain size distributions and to constrain transport processes, pathways, as well as potential source areas (Prins and Vriend, 2007; Bokhorst et al., 2011; Nottebaum et al., 2015). The imbrication of grains in sediments, recorded by the anisotropy of magnetic susceptibility (AMS) allows the reconstruction of palaeowind directions and thus dust transport pathways (Lagroix and Banerjee, 2002; Nawrocki et al., 2006, 2018; Bradák, 2009; Bradák et al., 2020; Zeeden and Hambach, 2021). Heavy mineral assemblage data can be used to fingerprint loess deposits and compare them to potential source sediments. However, using heavy mineral distribution assemblages as a provenance proxy may bear the drawback of mixing multiple distinct sources (Mikulčić Pavlaković et al., 2011; Römer et al., 2016; Meng et al., 2019; Wolf et al., 2019). Recently, the advent of single-grain provenance techniques such as detrital zircon U-Pb dating allowed disentangling multiple dust sources and potentially sign specific proto-source terranes (Stevens et al., 2010; Újvári et al., 2012; Fenn et al., 2018; Pańczyk et al., 2020). However, the aeolian system is subjected to well-known gravitative distance sorting processes, potentially leading to an over- or underrepresentation of proximal sources, depending on the sample location (Schaetzl et al., 2018).

Loess geochemical provenance proxies circumvent misinterpretations of proximal and distal dust sources by including sedimentary components representative for all grain size classes and mineralogical components (Buggle et al., 2008; Skurzyński et al., 2020). Despite the aforementioned mixing of sources, distinct geochemical differences between loess sequences provide information about the detrital background of the sediment and therefore allow the interpretation of dust sources (Muhs, 2018). Usually, the geochemical compositions of loess deposits are compared to those of potential source areas to determine geochemical commonalities (Buggle et al., 2008; Obreht et al., 2016). Despite methodological limitations (Buggle et al., 2008; Schatz et al., 2015a), geochemical signatures provide a tool which is neither time-consuming nor expensive. Its application has enabled for example to distinguish dust sources of the Chinese Loess Plateau and loess deposits in southern China (Hao et al., 2010) and to exclude local sources for aeolian silt mantles on the Californian Channel Islands (Muhs et al., 2007, 2008). Campodonico et al. (2019) confirmed the dominance of Andean detritus as source material for Pampean loess deposits of Central Argentina, based on bulk geochemical evidence.

Due to its complex formation, loess represents a geological and lithological mixture of its source areas. The different mineralogical components of loess deposits are prone to a multitude of methodological issues and sensible to, e.g., instrumentation or pretreatment of samples (Goff et al., 2020). Differences between different geochemical analysis techniques, such as inductively coupled plasma mass spectrometry (ICP-MS) and x-ray fluorescence analyses (XRF) can account to up to 10% of the measured contents (Miyazaki et al., 2016). To exclude such biases (Demir et al., 2006), we only use XRF data from samples pretreated and measured at the Physical Geography Laboratory of the RWTH Aachen University (Germany).

Sampling

Sampling of all sections was conducted in a continuous sampling line. Before sampling, several decimetres of sediment were removed from the vertical face of the loess wall to avoid contamination from freshly weathered exposed loess or downslope relocated material. Samples were taken contiguously in 2 cm increments for URL, RAS, VLA, and BAK and in 10 cm increments for SAG. The reference sections discussed in our regional comparison were sampled in 3–10 cm resolution (Table with sampling resolution and number of samples is provided in Supplementary Table S1). Sampling was conducted with cleaned tools and samples were packed into sterile plastic bags.

Geochemical Analyses

All samples were analysed by means of an energy dispersive x-ray fluorescence (EDPXRF, hereafter XRF) analysis using a Spectro Xepos. The samples were sieved to the silt fraction (<63 µm) and dried at 105 °C for 12 h. A subsample of 8 g for each sample was mixed with 2 g Fluxana Cereos wax and homogenised in a shaker. The admixture was pressed to a pellet with a pressure of 19.2 MPa for 120 s. The measurements were conducted by means of a pre-calibrated method. Each sample was measured in duplicate and the pellets were rotated by 90° between measurements to avoid matrix effects. The results of every sample were averaged. Conspicuous samples with striking differences between both measurements were remeasured in duplicate to avoid analytical artifacts. Major elements are given in oxide form (%; Na₂O, Al₂O₃, SiO₂, K₂O, TiO₂, Fe₂O₃), whereas trace elements are expressed in mg/kg (V, Cr, Ni, Rb, Sr, Y, Zr, Ba, Th). A similar analytical procedure was applied for all other studies used in this comparison (Obreht et al., 2015; Bösken et al., 2018, 2019; Marković et al., 2020; Tecsa et al., 2020).

The results were visualised in bi- and depth plots to investigate differences and similarities in space and time. To visualise the relations between soluble and insoluble elements and therefore interpret post-depositional alterations of sediments due to weathering, the Al₂O₃-CaO*+Na₂O-K₂O (molar proportions) ternary (A–CN–K) diagram (Nesbitt and Young, 1984) was plotted for the investigated LPS to display commonalities and differences between the sites. CaO* reflects the content of Ca in silicates. The ratio of the total Ca content (CaO) and Na₂O in silicates is assumed to never exceed a value of 1. Therefore, when CaO content is higher than the Na₂O, the siliceous Ca content (CaO*) is set to the content of Na₂O (McLennan, 1993; Buggle et al., 2011).

Grain Size Analyses

For grain size analyses, the samples were dried at 35 °C, sieved to the fine earth fraction (<2 mm), and two subsets of each sample (between 0.1 and 0.3 g) were pretreated with 0.7 ml H₂O₂ (30%) at 70 °C for 12 h. This process was repeated until a bleaching of the material was visible (Allen and Thornley, 2004), but not longer than three days. To keep the particles dispersed during the analysis, samples were treated with 1.25 ml Na₄P₂O₇ (0.1 mol*L⁻¹) in an overhead shaker for 12 h. The grain size was determined with a Beckman Coulter LS 13320 laser diffractometer using the Mie-theory (Fluid RI: 1.33, Sample RI: 1.55, Imaginary RI: 0.1; Özer et al., 2010; Nottebaum et al., 2015; Schulte et al., 2016). The results of grain size analyses are visualised as heatmaps. This visualisation allows interpretations of shifts in grain size distributions with depth (Schulte et al., 2016; Schulte and Lehmkuhl, 2018). From each loess sequence, an interval of 30 cm apparently unaffected by post-depositional processes such as weathering or pedogenesis was chosen. These intervals were determined using the grain size distributions and magnetic susceptibility. Over each of these intervals, the average of the mean grain size was calculated as representative for undisturbed loess of the respective section.

Environmental Magnetic Analyses and Correlative Age Model

Samples for environmental magnetic analyses were dried, ground using a porcelain mortar, and packed into plastic boxes. The samples were fixed with cotton wool to avoid movement of particles during transport and measurement. The volumetric magnetic susceptibility was measured at two frequencies (300 and 3,000 Hz) in a static field of 300 mA/m, using a Magnon International VSFM at Bayreuth University (Germany). The results were corrected for drift and the weak diamagnetism of the sample boxes and normalised for density. Environmental magnetic data was used for correlation of the investigated LPS to create a tentative regional age model (Supplementary Material 3, Figure S1).

A correlative age model (Supplementary Material 3, Figure S1) was established based on published age models for the investigated LPS URL and VLA (Obreht et al., 2017), RAS (Zeeden et al., 2018), and BAK (Scheidt et al., 2021). SAG was correlated based on tephrochronological evidence and trends in magnetic susceptibility. The LPS were correlated on the MIS scale spanning from MIS 3 to 1. To ensure comparability to the reference regions, only (parts of) sections covering the same or comparable time frames were used (Marković et al., 2007, 2020; Obreht et al., 2015; Bösken et al., 2018, 2019; Tecsa et al., 2020; Pötter et al., 2021).

Results

Geochemistry

Within the geochemical data set, the tephra layers show the largest fluctuations in their geochemical composition compared to the surrounding sediments, especially for RAS and SAG (see below). However, to ensure the full integrity of the presented records, the data of the samples affected by the tephra are presented together with other samples, but not included into the discussion. The same applies for samples where volcanic glass shards are incorporated into the sediment in considerable amounts, influencing the geochemical signature noticeably. In the data visualisations, extreme outliers caused by tephra particles are displayed together and additionally marked.

The geochemical compositions of the investigated LPS exhibit broad variations (Supplementary Table S2). Mean values for the Na₂O content vary between 0.92% for URL and 1.34% in SAG. The latter shows the highest values within the CI tephra layer (max.: 4.84%). Mean Al₂O₃ contents spread between 10.28% for URL and 12.01% for SAG. BAK shows the highest variability of SiO₂ contents, spanning from 39.84 to 66.72%, with an average of 50.41%. The other sections, with exception of SAG (mean: 52.78%) exhibit lower contents (URL: 43.46%; RAS: 45.61%; VLA: 48.73%). The mean content of K₂O of the southern sites of URL (1.68%), RAS (1.75%), and VLA (1.74%) are generally lower than those of the northern sites of BAK (1.89%) and SAG (2.71%). The high values in SAG again occur in the thick CI tephra bed. TiO₂ values for SAG again show the broadest spread (min.: 0.48%; max.: 0.84%), whereas the mean value (0.68%) is the lowest of the investigated LPS (URL: 0.78%; RAS: 0.83%; VLA: 0.79%; BAK: 0.74%). The range of mean Fe₂O₃ contents of 0.16% for all sections is narrow (VLA: 4.05%; RAS: 4.21%). The trace element compositions as well as the data for the reference section can be found in Supplementary Table S2.

The bi-plot visualisations of element contents and ratios show distinct patterns when comparing the investigated LPS. The ratios of major elements (Figures 3, 4) overlap between the sites. This is especially the case in the Al₂O₃/K₂O ratio (Figure 3). In this ratio as well as the Al₂O₃/SiO₂ ratio, URL and RAS show very similar trends. Especially URL is characterised by the lowest values for both elements, and the point clouds of both sites are roughly parallel to the diagonal of the plot. VLA shows higher SiO₂ contents compared to URL and RAS, especially in samples with lower Al₂O₃ values. The overlap of these three sites is generally quite large, whereas BAK and SAG are characterised by higher SiO₂, K₂O, and Al₂O₃ contents. The TiO₂ values are depleted for these two sites. For SAG, the CI tephra layer is visible in all major elements with either extremely elevated or depleted values. The TiO₂/Ni plot, adapted from Újvári et al. (2008), shows different trends within the five sites. BAK, SAG, and URL show elongated point clouds with a trend parallel to a positive correlation between Ti and Ni. VLA, however, shows a more circular point cloud with a trend roughly perpendicular to BAK and URL. The pattern of RAS shows both, a circular and a linear component in its pattern.

FIGURE 3. Bi-plots of major element composition for all investigated sections and reference regions. Symbols represent the here investigated loess-palaeosol sequences (Urluia (URL): green; Rasova (RAS): blue; Vlasca (VLA): black; Sageata (SAG): gray; Balta Alba Kurgan (BAK): red. Extreme outliers from the CI tephra layers and mixed loess-tephra samples are indicated by the letter T. The dashed lines show the geochemical composition of the reference regions: Lower Austria = orange; Lake Balaton (Central Hungary) = purple; Bükk Mountains (Northeast Hungary) = blue; Vojvodina = green; Black Sea Shore (Ukraine) = yellow (for location of the reference sections, see Figure 1). Note that extreme outliers for the reference regions are excluded.

FIGURE 4. Bi-plots of major element ratios for all investigated sections and reference regions. Legend according to Figure 3. Note that extreme outliers for the reference regions are excluded.

The Na₂O/Al₂O₃ vs. K₂O/Al₂O₃ diagram shows a large overlap between URL, RAS, VLA, and most of the BAK record (Figure 4). The rest of BAK as well as the majority of SAG are shifted toward higher K₂O/Al₂O₃ and lower Na₂O/Al₂O₃ ratios. In RAS and SAG, the CI tephra layer is clearly visible in strongly elevated values. The distinction between the investigated LPS is more evident in the other major element ratios. Especially the SiO₂/TiO₂ vs. K₂O/Al₂O₃ diagram shows a clear differentiation between most sites. RAS and URL again are overlapping, whereas VLA is characterised by higher SiO₂/TiO₂ ratios. BAK shows a narrow range, with similar K₂O/Al₂O₃ and elevated SiO₂/TiO₂ ratios. The range of the SAG data, even when the tephra samples are not included, is quite high for both ratios. Whereas the SiO₂/TiO₂ is comparable to BAK, with higher maximum values, the K₂O/Al₂O₃ ratio shows a spread of 0.17 to more than 0.25. The Al₂O₃/TiO₂ and to lesser extent the Fe₂O₃/TiO₂ ratio show a similar pattern.

Trace element compositions vary strongly between the five LPS (Figure 5). Whereas URL and VLA overlap in large parts, RAS shows elevated Zr/Ni and Zr/V ratios. The samples from the tephra layer in SAG have, similarly to RAS, extremely high values in these ratios. The loess layers are characterised by depleted values. The Cr/V vs. Y/Ni diagram shows a trend subparallel to the x-axis. Here, SAG exhibits the lowest Y/Ni values, but the CI tephra samples again show very high values for VLA, BAK, and SAG compared to the other LPS. With the samples covering the CI tephra excluded, the Cr/V ratios for these three LPS vary between 0.5 and 1.5. RAS and URL show a large spread, from values around 1 to more than 4. The Ba/Sr vs. Rb/Sr diagram depicts an almost linear trend between these two ratios for most of the sample. BAK is slightly shifted toward lower Rb/Sr ratios. SAG, which shows the largest spread in both ratios, also has an almost horizontal trend.

FIGURE 5. Bi-plots of trace element composition for all investigated sections and reference regions. Legend according to Figure 3. Note that extreme outliers for the reference regions are excluded.

The A–CN–K diagram shows point clouds subparallel to the CN joint for all sections (Figure 6). Whereas this is the case for all samples from URL, VLA, and BAK, RAS and SAG have a strong vertical component in their pattern. The point clouds of the respective sections have relative offsets to each other. VLA shows the largest spread on the Al₂O₃ joint, whereas RAS is shifted toward the CN joint, and BAK and especially SAG are shifted toward the K₂O apex. The point clouds of VLA and BAK have very similar starting points trending toward the composition of the upper continental crust (UCC).

FIGURE 6. A–CN–K plot for the investigated sections (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). Compositions of minerals as well as grain size, weathering and sorting trends are modified according to Buggle et al. (2008). UCC = Upper Continental Crust. Extreme outliers from the CI tephra layers and mixed loess-tephra samples are indicated by the letter T.

Grain Size Analyses

The results of the granulometric analyses for the investigated LPS are depicted in Figure 7. All five sites are dominated by strong modes of coarse silt (20–63 µm). However, the five sites show distinguished patterns of grain size distributions. The base of the investigated section at URL is characterised by a distinct mode of coarse silt, with a slightly fining tendency up-section (Figure 7A). Between 7 and 6 m below surface (b.s.), the contents of finer particles, particularly fine silt, increase. Around the CI tephra layer (approx. 6 m b.s.), a thin but distinct peak in finer particles occurs. The mode of coarse silt is less well expressed but still visible. Above the tephra, the sediment becomes coarser until 5 m depth and the mode of coarse silt as well as fine sand contents increase. Between 4 and 3 m b.s., two coarser layers with increased amounts of sand occur. Above 3 m, a decreasing trend in grain size can be seen toward the top. Generally, the grain size distributions are well sorted.

FIGURE 7. Grain size distributions for the investigated sections, sorted from south to north (A: Urluia; B: Rasova; C: Vlasca; D: Sageata, E: Balta Alba Kurgan). F: Grain size distributions for unweathered loess from the investigated sections [URL: Urluia (3–3.3 m); RAS: Rasova (2–2.3 m); VLA: Vlasca (5.1–5.4 m); SAG: Sageata (2–2.3 m); BAK: Balta Alba Kurgan (3.5–3.8 m)].

RAS generally shows a coarser distribution than URL, as the main mode is between coarse silt and fine sand (Figure 7B). The basal part is characterised by a strong mode of these two grain size classes. The CI tephra is clearly visible as a strong distinct peak of coarse grain sizes, with a strong mode in fine sand and increased middle sand contents. Above the tephra, the strong mode continues until a depth of around 5.5 m. Within the 50 cm above, two distinct coarse peaks are visible. Above that, the grain size distribution is well sorted, with a fining trend toward the top. Despite this trend, no major variations are visible.

The mode of coarse silt is also visible in VLA, where the basal part shows high percentages of this class (Figure 7C). Between 11.5 m b.s and the CI tephra at a depth of 9.1 m, increased proportions of medium and fine silt occur. The coarse silt mode dominates with increased grain size between 7 and 6 m b.s.. Above that, the distributions remain rather consistent until a depth of around 2 m, where the fine silt and especially the clay contents increase toward the top of the profile. The loess packages with a prominent mode in coarse silt are very well sorted.

The grain size distribution of SAG is very complex (Figure 7D). The base of the section, which overlies alluvial loams, shows no clear sorting. Within the thick tephra bed, a prominent mode within the fine sand fraction is visible. This continues above the ash layer, until a depth of around 3.1 m. Between three and 2.5 m b.s., the sediment fines up and shows a mode between coarse silt and fine sand. Above, the sediment is coarser and less well sorted. Between 1.7 and 1 m, the sediment is relatively well sorted and the mode ranges in the coarse silt fraction. Above 1 m b.s., middle clay contents increase.

The mode of the grain size distribution of BAK is from the base to approx. 7 m b.s rather constant, ranging in the coarse silt grain size (Figure 7E). Within this package, however, there are four layers of increased clay content as well as a coarser interval. Above that, a slight coarsening of the sediments texture is noticeable, with a more pronounced coarse silt mode. The sediment is generally very well sorted. Until 4 m b.s., four coarser layers are intercalated. Between 4 and 3 m, a coarsening trend is visible, whereas between 3 and 2 m, the clay content increases. The uppermost 2 m of the profile exhibit a mode of coarse silt to fine sand, with varying clay contents. However, since this is the anthropogenic infill of the kurgan, it is not included in the provenance discussion.

The mean grain size distributions for unaltered loess intervals are shown in Figure 7F. The distributions of URL and BAK show almost identical unimodal courses, with a prominent peak in the coarse silt fraction. Both sites also exhibits shoulders in coarse clay, as well as fine sand. The peak in VLA is shifted slightly toward coarser grain sizes, whereas the clay and fine silt contents are lower than in URL and BAK. The mode for RAS lies on the boundary between coarse silt and fine sand, and the content of fine sand is higher than at the other sites. The most complex distribution is found in SAG, where the mode is less expressed and not as narrow as for the other sites. The content of fine sand is considerably higher than, e.g., in BAK. SAG also shows the highest contents of clay, especially in the coarse and middle clay fraction.

Discussion

The Geochemical Imprint of the Campanian Ignimbrite Tephra

Tephra layers alter the chemical composition of the sediment in which they are incorporated (e.g., Obreht et al., 2016). The CI tephra, which is used here as a chronostratigraphic marker horizon (Constantin et al., 2012; Fitzsimmons et al., 2013; Veres et al., 2013; Obreht et al., 2017), leaves a clear geochemical signature compared to the embedding loess. Due to its trachytic characteristics the CI tephra can be identified using geochemical proxies of bulk sediment samples even when present as cryptotephra (see Obreht et al., 2016). The influence of the CI tephra on the geochemical composition of the embedding loess varies strongly between the investigated sections. In BAK and URL, where the CI tephra is only preserved in patchy lenses or mixed with loess, respectively, the geochemical composition of the embedding loess samples seems less affected (e.g., Figures 8, 9). Contrary, for VLA the Zr/V and Zr/Ni ratio for example are increased since the ash layer is preserved as a laterally continuous stratum of several cm thickness (Figure 10, Supplementary Figure S3). The strongest and most prominent peaks in Zr/V can be seen in RAS (30 cm thick tephra) and SAG, where the tephra bed consists of pure pyroclastic material. Especially in SAG, where the CI tephra is more than 50 cm thick, parts of the layer consist of pure ash beds, which strongly influence the geochemical signal of the analysed samples. Since volcanic glass shards are usually incorporated into subsequently deposited dust (Obreht et al., 2017), the geochemical effects of the tephra can be found even some cm to dm above the actual layer forming a gradual transition to the overlying loess. Geochemical proxy data, especially those based on diagnostic elements, can be easily used to detect volcanic material in loess sequences. In the here investigated LPS, the tephra layers show typical characteristics imprinted in the geochemical data. The Al₂O₃/TiO₂ ratio in RAS, e.g., shows a rapid increase at the lower boundary of the tephra and a gentle decline toward the above-lying loess, indicating incorporation of volcanic glass shards into the sediment. In other localities, such as BAK where the tephra is preserved as a patchy layer, the geochemical influences up-profile of the tephra are reduced (Figure 10). However, these large variations originating from the volcanic ash, both in samples of pure ash beds as well as varying amounts of glass shards in the sediment, cannot be interpreted as true loess provenance signals but must be kept in mind when discussing transport pathways of contemporaneous dust deposition.

FIGURE 8. Depth variations of the SiO₂/TiO₂ ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

FIGURE 9. Depth variations of the K₂O/Al₂O₃ ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

FIGURE 10. Depth variations of the Zr/V ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

Geochemical Signature and Mineralogical Composition of the Lower Danube Loess

Some mineralogical approximations can be made based on loess geochemical data (Muhs, 2018). Generally high SiO₂ contents show the predominance of quartz in the mineralogical composition of all five LPS. The varying contents, both between and within the LPS, are mainly governed by differences in the mineralogical compositions and grain size distributions (Figure 8, Supplementary Figure S2). The clear positive relationship between Al₂O₃ and SiO₂ (Figure 3) indicates the presence of aluminosilicates such as mica or feldspars in all LPS investigated here. This is supported by the positive relation between Al₂O₃ and K₂O, which suggests that the aluminum content is mainly bound to phyllosilicates or K-feldspars. The presence of Ti-bearing phyllosilicates is also supported by the weak positive correlation between Ti and Ni (Újvári et al., 2008), since Ni tends to be enriched in these minerals. This is, to varying degrees, visible in all sections (Figure 3).

The A–CN–K ternary diagram developed by Nesbitt and Young (1984) has been widely applied to LPS from varying environments to investigate influences of chemical weathering and mineralogical compositions on the geochemical composition of sediments (e.g., Yang et al., 2004; Schatz et al., 2015b; Vlaminck et al., 2018; Skurzyński et al., 2019). All investigated LPS show a trend of plagioclase weathering (Figure 6). This trend indicates a rather stable ratio of K-feldspar and albite for most of the samples. The (sub-) parallel point clouds of especially URL is characteristic for loess deposits of the Dobrogea, such as the Mircea Voda (Buggle et al., 2008) or the Costineşti LPS (Tugulan et al., 2016). The shifts of BAK and especially SAG toward the right edge of the plot indicates either increased feldspar weathering or different mineralogical composition of source areas compared to the other sites (Vlaminck et al., 2018). In the latter case, it can be interpreted to reflect more K-feldspar rich parent material (Újvári et al., 2008). The abundance of K-bearing minerals, such as K-feldspars, muscovite, or illite is indicated for all sections by the shift toward the K-apex in the A–CN–K ternary diagram (Figure 6) and a strong correlation between the Rb/Sr and Ba/Sr ratio (Figure 5; Újvári et al., 2008). High Ba/Sr ratios also show the presence of plagioclase, as indicated by the A–CN–K plot. However, the shifts of SAG and BAK toward higher Rb/Sr and Ba/Sr values, as well as their extreme values, can be interpreted as Sr depletion due to leaching processes during soil formation (Újvári et al., 2008).

The geochemical and mineralogical compositions of loess and palaeosols is not solely determined by the composition of the parent material. In subaerial deposits such as loess, post-depositional processes including weathering or pedogenesis frequently alter the material’s geochemical composition (Vlaminck et al., 2018; Bösken et al., 2019; Skurzyński et al., 2020). As seen in the A–CN–K plot (Figure 6), the majority of the samples show a discernible plagioclase weathering trend. Especially soluble or easily weathered compounds and minerals are prone to depletion by leaching. The Ba/Sr and the Rb/Sr ratios, e.g., are used in loess research as weathering indices, since Sr is easily leached during carbonate dissolution while Ba and Rb are rather immobile elements, which are relatively enriched during leaching (Chen et al., 1999; Buggle et al., 2011). When discussing loess provenance based on geochemistry, palaeosol samples are often excluded to overcome issues of geochemical alteration during pedogenesis (Buggle et al., 2008; Újvári et al., 2008; Schatz et al., 2015a). However, a clear distinction between loess, pedogenically overprinted loess and fully developed palaeosols is not always straightforward in semi-arid environments. Therefore, we included all samples and plotted provenance proxies as depth functions to complement the bi-plot visualisations. In these plots (Figures 8–13), data from the reference regions are displayed as box-whisker plots on top of the graphs. The interquartile ranges are indicated as shaded boxes of the respective color in the depth plots. The effects of solution and weathering can be seen in a variety of element contents and ratios. The Na₂O/Al₂O₃ ratio is depleted in the recent topsoil at URL, VLA, and BAK, and increased and/or constant in the loess intervals (Supplementary Figure S4). Post-depositional alterations by weathering or pedogenesis are also expressed by the (neo-)formation of (aluminous) clay minerals or iron oxides. As the former may lead to a relative enhancement of aluminum oxides (Figure 11, Supplementary Figure S2), the latter leads to an increase in iron contents in the affected sediment (Figure 12). Since the formation if these minerals is usually linked to pedogenesis, the elevated Al and Fe contents can be related to the clay content (Supplementary Figure S2). The Al₂O₃/TiO₂ ratio plotted as a depth function (Figure 11) as well as the direct comparison between the Al content and the grain size (Supplementary Figure S2) reflects this issue and must therefore be carefully discussed. Shifts in these ratios, such as those observed below the CI tephra in VLA, cannot be interpreted as provenance signal as they are of post-depositional nature.

FIGURE 11. Depth variations of the Al₂O₃/TiO₂ ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

FIGURE 12. Depth variations of the Fe₂O₃/TiO₂ ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

Since trace elements such as Ni and V are compatible elements within igneous systems, they can be used as indicators for the provenance of rocks and sediments (Újvári et al., 2014; Obreht et al., 2016). Ratios based on these trace elements reveal striking similarities between URL and VLA (Figures 3, 5). RAS shows large overlaps with the two sites in, e.g., the Th/Ni or Th/V ratios (Figure 5). However, the Zr based provenance ratios vary substantially between those sites. These differences can be explained by the enrichment of zircon minerals in RAS; in detrital sediments Zr is almost entirely bound to zircons. These heavy minerals can be enriched in loess by selective deflation of less dense particles. This so-called hydraulic sorting (Malusà et al., 2016) causes enrichment of zircons which in turn leads to elevated Zr contents in the geochemical screening of such samples. Therefore, the elevated Zr/Ni and Zr/V ratios can be interpreted as an indicator of sorting processes (Figure 5; Újvári et al., 2008).

The hydraulic sorting of minerals by water or air not only affects the trace element geochemistry, such as the Zr content of loess deposits, but is also expressed in other geochemical proxies (Ohta, 2004). The position and the shape of the point cloud of a sample set in the A–CN–K plot can be indicative for such sorting processes. Due to the abundance of Al in aluminous clay minerals, finer, often weathered sediments tend to plot close to the Al₂O₃ apex of the diagram, whereas coarser sediments show a contrary trend (Nesbitt and Young, 1996; Nesbitt et al., 1996). When the grain size differences within a LPS are solely driven by (plagioclase) weathering, point clouds plot (sub-) parallel to the CaO*–Na₂O join (see Vlaminck et al., 2018; Bösken et al., 2019; Skurzyński et al., 2019). However, when hydraulic sorting skews the grain size distribution, point clouds in the A–CN–K plot tend to develop (sub-)vertical shapes, trending toward the Al₂O₃ apex (Ohta, 2004; Buggle et al., 2008). For BAK and URL, the point clouds sub-parallel to the CaO*–Na₂O join indicate no major influence of hydraulic sorting processes. VLA shows a similar pattern, but the point cloud shows a slightly stronger vertical component (Figure 6). The strongest evidence for sorting processes can be found in RAS and SAG, with strong vertical components. The direct vicinity of RAS and SAG to the recent courses of the Danube and the Buzău might be the main reason for these sorting processes. The high intensity of these processes in RAS can be explained by enhanced fluvial activity of small Danube tributaries acting as dust sources, as this section is exposed in the recent dry valley of such a tributary and its base shows a sharp fluvial erosional contact (Zeeden et al., 2018). The pre-aeolian setting is therefore fluvial as it is in SAG. This section lies on top of alluvial loams, which were presumably deposited in river ponds in the floodplain of the palaeo-Buzău. Therefore, the general high fluvial activity strongly influenced the formation of these loess deposits during the Late Pleistocene. The alluvial deposits near the sites of RAS and SAG, and to a lesser extent of VLA, acted as local sources for coarse material deposited in direct vicinity while finer particles were transported over longer distances. The grain size distributions also reflect these processes, since the three sites show generally coarser distributions than URL and BAK (Figure 7).

Generally, the investigated LPS show striking similarities in some of the presented provenance proxies. Especially the sites from the Dobrogea, URL and RAS, show significant common variabilities, indicating geochemically similar parent material as expressed in the Cr/V ratios (Figures 5, 13). This ratio serves as a good indicator for the contributions of detrital material derived from mafic to ultramafic rocks since Cr is almost entirely bound to mafic strata (Zimmermann and Bahlburg, 2003; Mikes et al., 2006; Újvári et al., 2008). For the other sites, the values of this ratio are in a rather narrow spectrum without large variations. However, data from URL and RAS exhibit large scatter of this ratio with maximum values of >4. This feature is in accordance with elevated Cr contents compared to the UCC established in loess sequences from the southern Dobrogea (Tugulan et al., 2016), which indicates the influence of distinct sources for the loess deposits of the area (Figure 5, Supplementary Figure S7).

FIGURE 13. Depth variations of the Cr/V ratios for the investigated loess-palaeosol sequences, sorted from south to north (URL: Urluia; RAS: Rasova; VLA: Vlasca; SAG: Sageata, BAK: Balta Alba Kurgan). The box-whisker plots on top show the composition of the sections from the reference regions: I: Lower Austria; II: Vojvodina; III: Balaton (Central Hungary); IV: Bükk Mountains (Northeast Hungary) and V: Black Sea Shore (Ukraine). The interquartile area (between 25 and 75% quartile) are shown as shaded boxes. The correlative age model as well as the CI tephra are shown for temporal comparison.

VLA shows similarities to URL in several of the presented geochemical provenance proxies. The correspondence with URL is evident for example in the Na₂O/Al₂O₃ vs. K₂O/Al₂O₃ and the Al₂O₃/TiO₂ vs. SiO₂/TiO₂ plot (Figure 4), suggesting similar sourced material contributed to both records. RAS, however, shows some overlap with URL, but also exhibits a distinct pattern, e.g., in the Al₂O₃/TiO₂ ratio, which is depleted compared to the other sites. Besides the general commonalities identified within the studied LPS, significant variations are also discernible. It is noteworthy that the variations in the K₂O/Al₂O₃ vs. SiO₂/TiO₂ diagram tend to plot closer to the data from BAK, with which VLA also has large overlaps. SAG, however, shows some commonalities to BAK, but is distinguishable from the other sites, indicating a substantial contribution from a different dust source area. Therefore, we propose similar dominant source areas for the URL, RAS, and VLA sequences, with a more diverse input from multiple sources for the latter. For SAG, different and especially local sources, such as the proximal input from the Carpathian Bending can be assumed while BAK represents a more mixed provenance signal of different sources. While the use of multiple insightful bulk geochemical proxies has enabled us to sufficiently differentiate loess deposits in the Lower Danube Basin, further reference data is needed to assign the respective dust sources.

Estimations of Sedimentary Pathways of the Lower Danube Loess

As stated previously, at this stage of research we compare our results with geochemical data obtained using the same instrumentation and pretreatment from other loess deposits in Central and Southeastern Europe (i.e., the Middle Danube Basin (MDB) and Ukraine), where previous multi-proxy provenance studies have constrained potential dust sources (Supplementary Table S1). This comparison potentially allows to indirectly assign potential sources to the loess deposits in the Lower Danube Basin, a region yet largely underrepresented in loess provenance studies. The discussed provenance studies were chosen according to their proximity to the reference sections from Austria, Hungary, the Vojvodina, and Ukraine. However, local and regional differences in deflation and accumulation cannot be excluded. Therefore, only a tentative identification of potential source areas can be postulated in addition to the discussed geochemical characterisation. An overview of the discussed and additional provenance (Supplementary Table S3) and palaeowind direction studies (Supplementary Table S4) can be found in the Supplementary Material.

The Danube Alluvium (DA) Sedimentary Pathway

The large overlap of especially URL, RAS, and VLA with the geochemical compositions of some of the reference sites in Lower Austria, Central Hungary, and the Vojvodina, indicates similar geochemical signatures of the dust sources to the respective loess deposits (Figures 8, 9, Supplementary Figure S4). The alluvial sediments of the MDB are recognized as the main source for the loess deposits in the area (Újvári et al., 2012; Obreht et al., 2019). This is in accordance with findings from Buggle et al. (2008), who demonstrated geochemical similarities between the Mircea Voda site in the LDB and the Batajnica/Stari Slankamen composite sequence of the Vojvodina. The authors interpreted these similarities as an indicator for the importance of the Danube alluvium (DA) as the main sediment supplier and dust source in the LDB, corresponding to the MR mode of loess genesis (Li et al., 2020). However, the geochemical composition of the alluvial sediment changes laterally, as major tributaries in the MDB contribute to the sediment load of the Danube (Nenadić et al., 2016; Vandenberghe et al., 2018; Mitrinović et al., 2020). The influence of other local to sub-regional source areas for the loess of the Middle Danube Basin is evident in the range in the geochemical data of the reference sections from the MDB. The importance of this influence increases with distance from the course of the Danube (Basarin et al., 2011; Obreht et al., 2014, 2016, 2019). Accordingly, the loess deposits of Northeast Hungary, which cover the upper to middle reaches of the Tisa catchment (Schatz et al., 2015a; Bösken et al., 2019), show different geochemical signatures as expressed in the K₂O/Al₂O₃ vs. SiO₂/TiO₂ diagram (Figure 4). These patterns are also evident for other LPS of the area, which show considerable geochemical differences to deposits in direct proximity to the Danube (Schatz et al., 2015a). Enhanced contributions from other local sources, such as the Bohemian Massif for Austrian loess (Újvári et al., 2013; Újvári and Klötzli, 2015), exposed local sandstones in Hungary (Újvári et al., 2008, 2012; Sebe, 2013) or from the Deliblatska Peščara sand field in Serbia (Obreht et al., 2015) were detected using various methods. However, for most loess deposits in the Carpathian Basin, the Danube alluvium appears to be the main dust source. The large overlap between the geochemical compositions of especially VLA, RAS, and URL with sections from the Middle Danube Basin indicates that the Danube alluvium in the LDB area is geochemically similar to the upper reaches of the river and acted as the dominant regional source area for dust deflation and ultimately, loess formation. Geochemical differences are either due to local sources within the LDB, but can also be attributed to palaeoclimatic gradients (Buggle et al., 2008).

The Eastern Carpathian (EC) Sedimentary Pathway

The SiO₂/TiO₂ ratio also shows a good agreement between BAK and SAG and the reference section of Kurortne (Black Sea Shore) in southwest Ukraine (Tecsa et al., 2020), which is visible in the K₂O/Al₂O₃ vs. SiO₂/TiO₂ diagram as well (Figures 4, 8). In general, two sediment pathways were reconstructed for Ukrainian loess covers. On the one hand, the loess deposits of Central Ukraine are linked to the glaciofluvial (GF) sedimentary pathway (Buggle et al., 2008). These deposits, inter alia located in the catchment of the Dnieper River (Veres et al., 2018), consist of aeolian dust entrained from the glaciofluvial outwash plains by strong northerly winds (Buggle et al., 2008; Nawrocki et al., 2019). On the other hand, the loess deposits of southern Ukraine rather correspond to the Eastern Carpathian (EC) sedimentary pathway sourced from the sediments of the Dniester River and its tributaries. Northwesterly palaeowinds have been reconstructed for dust transport in this region (Nawrocki et al., 2006, 2018; Pańczyk et al., 2020). Although the GF pathway cannot be ruled out completely for Kurortne, its location on the southern Ukrainian shore of the Black Sea (Figure 1; Tecsa et al., 2020) close to the Dniester estuary suggests the EC pathway to account for loess deposits in this area. Therefore, the striking overlap with especially BAK suggests similar source areas for both sections and agrees with the reconstructed strong and consistent N-NW palaeowinds in the area. This is supported by residual forms of loess dunes, so-called gredas (Rozycki, 1967; Leger, 1990; Mason et al., 1999), which can be found in the entire eastern foreland of the Carpathians, and which are usually oriented NW-SE (Figure 14; Rozycki, 1967; cited in Leger, 1990). Although detrital material derived from the Carpathians has a strong influence on the sediment yield of the Danube, the sediments of tributaries draining the (Outer) Eastern Carpathians can be distinguished by geochemistry (Buggle et al., 2008) and detrital zircon age assemblages (Ducea et al., 2018). Therefore, the differences between BAK and loess deposits further south in the LDB such as VLA and URL can be linked to the dominance of the EC sedimentary pathway for the BAK area.

FIGURE 14. A: potential sedimentary pathways for loess deposits in the Lower Danube Basin (DA: Danube alluvium; CB: Carpathian Bending; BS: Black Sea; EC: Eastern Carpathians; GF: glaciofluvial). Modes of loess genesis according to Li et al. (2020). Residual forms of gredas and longitudinal dunes in the Middle and Lower Danube Basin are adapted from (Rozycki 1967; cited in Leger, 1990). The division of the Carpathians is given to demonstrate different (proto-) source areas such as the flysch terrains of the Outer Eastern Carpathians, the volcanic chains of the Inner Eastern Carpathians or the crystalline and metamorphic bed rocks of the Southern Carpathians. The study sites are shown as triangles. B: surface trends for zircons (pink) and hornblende minerals (green). The isolines represent surface concentrations (in %) and are derived from statistical computation of mineralogical analyses of hundreds of surface samples (Andar and Codarcea, 1981; cited in; Jipa, 2014). The arrows indicate reconstructed transport directions (pink: towards Danube; green: towards Carpathians), based on the heavy mineralogical evidence.

The Carpathian Bending (CB) Sedimentary Pathway

Another potential mode of origin for loess at SAG and BAK is the Carpathian Bending (CB) sedimentary pathway, with primary material being eroded from the bending sector of the Eastern Carpathians and transported by rivers like Buzău or Ialomiţa. This pathway corresponds to the MR mode of loess genesis (Li et al., 2020). There are no direct provenance studies on loess deposits from this region, but the floodplains of the rivers draining the Southern Carpathians to the west of the CB were detected as potential source areas based on heavy mineral composition data (Jipa, 2014). Decreasing zircon concentrations toward the Danube plain can be linked to the decreasing influence of sediment input from the Carpathians where detrital material is enriched in zircons. The opposing trend in hornblende contents was interpreted as the transport direction of dust deflated from the Danube alluvium, as the hornblende minerals are enriched in the deposits of the Danube (Figure 14B; Andar and Codarcea, 1981; cited in Jipa, 2014). Figure 14B shows that BAK and SAG are located in a transitional area, where the influence of the Eastern Carpathians is strong, but also some amount of hornblende can be found within the loess deposits, pointing to the Danube alluvium as a potential source area. The EC pathway is supported by high zirconium concentrations for both sites but especially in SAG (Supplementary Table S2). As stated above, the inferred zircon enrichment is indeed an indicator for sorting processes at this site. However, the fact that according to the A–CN–K plot (Figure 6) not all samples of the SAG section are affected by sorting, the generally high values could be a sign for the input of detrital material from the Carpathians for these loess deposits. The overlap with data from northeast Hungary in proxies such as the Al₂O₃/TiO₂ ratio (Figures 4, 11) or the elevated Zr contents (Figure 5 and Supplementary Table S2) support the influence of detrital material from the Carpathians, as northeast Hungary is strongly controlled by the Tisa and its tributaries (Thamó-Bozsó et al., 2014; Schatz et al., 2015a). Therefore, these areas experience a strong input from the Inner Eastern Carpathians and their extensive volcaniclastic forelands, which are also drained by the upper reaches of some main tributaries of the Danube in the LDB such as the Olt River.

The Black Sea (BS) and Glaciofluvial (GF) Sedimentary Pathways

The main characteristics of the BS and GF sedimentary pathways cannot be illustrated in detail due to the lack of reference data. However, the importance of these pathways has been discussed for loess deposits of the LDB and adjacent areas. Especially the BS pathway was linked to the formation of loess deposits in Bulgaria (Avramov et al., 2006; Jordanova et al., 2007). The complexity of the sedimentological composition of the Black Sea shelf, caused by multiple large rivers discharging into the basin, hinders a clear link to this pathway based on geochemistry. Large amounts of the sediments building the Black Sea shelf were, e.g., deposited by the palaeo-Danube, indicating similar geochemical and mineralogical compositions to the DA pathway (Ducea et al., 2018). The GF pathway was reconstructed for LPS in Poland and Ukraine (Buggle et al., 2008; Skurzyński et al., 2019, 2020; Pańczyk et al., 2020) and was mainly driven by enhanced aeolian activities of katabatic winds influenced by fluctuations of the Eurasian CIS (Nawrocki et al., 2019). Pańczyk et al. (2020) constrained a divide of dust provinces based on detrital zircon ages, which divides the GF from the EC mode in Poland and Ukraine. In combination with the reconstruction of atmospheric circulations for the last glacial (Meyer and Kottmeier, 1989; Obreht et al., 2017), a major input of entrained dust from the GF pathway seems unlikely, especially for the southern part of the LDB.

Spatial and Temporal Provenance Variations of Lower Danube Loess and Palaeoclimatic Implications

The provenance of Lower Danube loess deposits varies spatially. While the sites from the Dobrogea plateau and the Danube plain seem to be linked mainly to sediment input via the Danube alluvium, sites in the northern part of the LDB are influenced by a mixture of potential source areas, including the DA, the EC, and the CB sedimentary pathways.

The hypothesized sedimentary pathways for the Lower Danube loess highlight the importance of proximal dust sources, as previously stated for the Middle Danube Basin (Obreht et al., 2019). However, far traveled dust originating from, e.g., the Sahara Desert can contribute especially to the fine silt fraction of aeolian deposits (Stuut et al., 2009). Although the importance of these distal sources for Southeastern European dust records was discussed for interglacial periods (Longman et al., 2017), the proposed contribution to glacial loess is of minor importance (see Varga et al., 2013).

Changes of dust source areas over time are potentially linked to palaeoclimatic fluctuations driving dust source expansion and/or changes in dust transport pathways (Zhang et al., 2018; Nawrocki et al., 2019; Pańczyk et al., 2020; Skurzyński et al., 2020). For MIS 3, URL, VLA, and BAK show some variation in the zirconium-based provenance proxies (Figure 10, Supplementary Figure S3). Especially in VLA, the Zr/Ni and Zr/V ratios show a decreasing trend from the bottom of the investigated section towards the CI tephra, whereas above it the values remain rather constant. The same trend but to lesser extent is visible in BAK. Since the grain size and especially the median grain size do not show large variations during this period (Figure 7), zircon enrichment due to sorting seems unlikely. The decreasing ratios are rather an indicator for the weakening of the influence of the Eastern Carpathians as a source area during this period. The weakened input of material derived from the Eastern Carpathians can be explained by shifts of atmospheric patterns during this time. During the middle MIS 3, westerlies reached Southeastern Europe on a NW-SE trajectory (see also Staubwasser et al., 2018). During the late MIS 3, the Siberian High pressure system intensified, leading to a shift of westerlies southwards (Obreht et al., 2017). Local maxima in the Zr based ratios in VLA suggest temporally increased input of Carpathian detrital material, indicating short-term fluctuations in source material. In URL, the Cr/V record is clearly divided during MIS 3 (Figure 13). Below the CI tephra, some oscillations occur indicating the input of detrital material potentially originating from more mafic sources. Above the tephra, these oscillations decrease and are within the interquartile range of the reference regions upstream. In RAS, the Cr/V proxy shows large variations throughout the whole profile, with exceptions during mid-MIS 2 and MIS 1, indicating an almost constant influence of mafic sources for these intervals. The pattern during MIS 2 is similar to URL, but the values are generally higher. Accumulation effects can explain the higher values and the large variations through time. Cr rich material was eroded in the Dobrogea and transported via small rivers toward the Danube Plain. The sediment was accumulated near RAS, deflated and deposited as aeolian dust (Zeeden et al., 2018).

The early MIS 2 was generally characterized by stable source areas for the investigated LPS. URL, RAS and VLA are clearly dominated by the DA sedimentary pathway during this period, whereas BAK shows large overlaps with the EC pathway (e.g., Figure 8). In the middle of MIS 2, however, the geochemical evidence shows changes in all sequences, e.g., in the K₂O/Al₂O₃ ratio (Figure 9). URL, RAS, and BAK show a noticeable decline in this proxy, whereas VLA and SAG show increasing values. These changes could be explained by the shift toward prevailing easterlies in large parts of Europe during the LGM (Schaffernicht et al., 2020). This would imply that the Black Sea (BS) sedimentary pathway, or long-range transport of material deflated even further to the east, becomes progressively more important throughout MIS 2, possibly contributing dust into the aeolian system of the LDB more significantly than hitherto considered. This coincides with sea level low stands of the Black Sea during the LGM (Constantinescu et al., 2015). However, since no reference data from sediments of the Black Sea shelf is available for this study, this hypothesis remains tentative. The same applies for long-range dust transport from the Caspian Basin further east, which was postulated for southwestern Ukraine especially for glacial periods (Stephens et al., 2002).

Methodological Limitations and Potential for Further Research

The general limitations of this methodological approach, such as the averaging of input from multiple dust sources, must be considered when discussing geochemical source estimations. The admixture of material from different proto-sources in a large sedimentary system such as the Danube Basin, can complicate the interpretation of distinct dust sources. When using bulk sediment (<63 µm) geochemistry as a provenance tool, signatures of multiple sources may be averaged out and assigning distinct sources might not be feasible (Bird et al., 2015). To identify specific dust (proto-)sources, it is essential that future research complements this method with further analyses, preferably single grain techniques that allow for more conclusive discrimination between potential source sediments of similar geochemical composition (Stevens et al., 2010; Fenn et al., 2018). Nevertheless, the geochemical composition of bulk sediments can be used to either rule out sources or to assess the main variability in loess provenance based on the varying contents of diagnostic elements and identifiable geochemical similarities or differences (Buggle et al., 2008; Újvári et al., 2008), as presented here for a transect of records from the Lower Danube loess.

Summary and Conclusions

Our high-resolution geochemical study on several loess-palaeosol sequences from the Lower Danube Basin is among the first multi-record analyses contributing to the understanding of loess provenance and major sedimentary pathways in different sectors of the basin. Based on the presented data, we subdivide the basin into two major domains of loess geochemical signatures. The southern part of the basin including the southern Wallachian Plain and the Dobrogea with the LPS of Urluia, Rasova, and Vlasca can be distinguished from the northern Wallachian Plain with the LPS of Sageata and Balta Alba Kurgan. The boundary between these two domains appears to be transitional rather than a sharp contact.

The comparison with reference sections from Austria, Hungary, the Vojvodina, and Ukraine reveals that the Urluia, Rasova, and Vlasca records show striking geochemical similarities with loess records from the Middle Danube Basin. Since these similarities can be seen in many provenance proxies, a similar dominant geochemical signature of source areas in the Middle and Lower Danube Basin can be assumed. Furthermore, the geochemical composition of the Balta Alba Kurgan LPS rather resembles the one of Kurortne in Ukraine, pointing to similar source areas for these regions. Lastly, Sageata, situated in the floodplain of the Buzău River, shows rather unique geochemical patterns compared to the other sections, despite some overlaps with Balta Alba Kurgan.

The link of these geochemical findings with published provenance studies from the reference regions allows a tentative identification of potential source areas for the Lower Danube loess. It confirms the previous perception that the alluvial sediments of the Danube and its major tributaries act as a main storage and as source areas for dust deflation in the basin. The comparison to loess deposits in the middle and upper reaches of the Danube supports the Danube alluvium pathway hypothesis, especially for the southern sites of Urluia, Rasova, and Vlasca. Further to the north, other sources gain more importance. For Sageata, detrital material eroded from the Carpathian Bending area and redeposited by sediment-laden rivers such as the proximal Buzău River plays a major role as a local dust source. For Balta Alba Kurgan, the Eastern Carpathian sedimentary pathway is of particular importance. Here, material is deflated from the riverbeds draining the northern slopes of the Carpathians. Therefore, the dominant mode of loess genesis here is the mountain provenance-river transport loess sensuLi et al. (2020). The input of material deflated from the glaciofluvial outwash plains of Ukraine can neither be confirmed nor ruled out. However, we consider this sedimentary pathway, corresponding to the continental glacier provenance-river transport mode sensuLi et al. (2020) to be of minor importance, especially for the southern Lower Danube Basin. Our study confirms the general assumptions of prevailing northwesterly to westerly winds during the Pleniglacial in the area. These findings should be complemented using single grain techniques such as detrital zircon U-Pb ages to illuminate the sedimentary pathways accounting for dust deposition in the Lower Danube Basin in greater detail.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

SP, PS, and UH designed the study. SP compiled the figures and wrote the first draft of the manuscript. SP, DV, PS, and UH defined the methodology. PS supervised the laboratory work and generated the discussed data. SP, DV, JN, UH, PS, and FL among others conducted fieldwork. All authors discussed the data, visualisations, and text.

Funding

The investigations were carried out in the frame of the CRC 806 “Our way to Europe,” subproject B1 “The Eastern Trajectory”: “Last Glacial Palaeogeography and Archaeology of the Eastern Mediterranean and of the Balkan Peninsula,” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 57444011-SFB 806.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Marianne Dohms and her team for the laboratory framework. Alina Blume and Bruno Boemke helped vectorising parts of Figures 1, 14. We thank Dagmar Haase for providing the shapefiles of loess distribution for Figure 1. We thank Jacek Skurzyński for his input on the methodological discussion. DV’ work on the tephrochronology of investigated sequences was supported through a fellowship provided by Alexander von Humboldt Foundation. We thank Annett Junginger for the editorial work and MQ and RS for their helpful comments, which improved the manuscript substantially.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.600010/full#supplementary-material.

References