A Detailed Paleoclimate Proxy Record for the Middle Danube Basin Over the Last 430 kyr: A Rock Magnetic and Colorimetric Study of the Zemun Loess-Paleosol Sequence

In mid-latitude Eurasia, loess-paleosol sequences (LPS) provide the most widespread sedimentary records of Quaternary paleoenvironmental evolution. In the Middle Danube Basin (MDB), these archives cover at least the last million years of climate history, and occasionally contain archeological findings. The studied Zemun LPS is located on the right bank of the Danube in Northern Serbia. The site was declared as a protected site, based on Paleolithic artifacts found on the riverbank and stemming from unknown stratigraphic levels of the loess cliffs exposed along the Danube. The present study aims to provide a stratigraphic, paleoenvironmental, and temporal context for the Zemun LPS by means of environmental magnetic and colorimetric methods. Our investigations result in a chronostratigraphic scheme allowing direct comparison with other well-established reference records in the MDB and elsewhere. Two potential tephra layers tentatively assigned to the so-called L2 and Bag tephras, which are both widespread in the MDB and beyond were investigated for their bulk magnetic properties. The resulting integrated age model suggests that the Zemun LPS records a detailed history of a quasi-continuous accumulation of mineral dust from Marine Oxygen Isotope Stage (MIS) 11–5a (c. 430–60 ka). The outcome of our integrative approach indicates a continuous aridification over the last four interglacial/glacial cycles and we discuss potential changes in seasonality over time.

In mid-latitude Eurasia, loess-paleosol sequences (LPS) provide the most widespread sedimentary records of Quaternary paleoenvironmental evolution. In the Middle Danube Basin (MDB), these archives cover at least the last million years of climate history, and occasionally contain archeological findings. The studied Zemun LPS is located on the right bank of the Danube in Northern Serbia. The site was declared as a protected site, based on Paleolithic artifacts found on the riverbank and stemming from unknown stratigraphic levels of the loess cliffs exposed along the Danube. The present study aims to provide a stratigraphic, paleoenvironmental, and temporal context for the Zemun LPS by means of environmental magnetic and colorimetric methods. Our investigations result in a chronostratigraphic scheme allowing direct comparison with other well-established reference records in the MDB and elsewhere. Two potential tephra layers tentatively assigned to the so-called L2 and Bag tephras, which are both widespread in the MDB and beyond were investigated for their bulk magnetic properties. The resulting integrated age model suggests that the Zemun LPS records a detailed history of a quasi-continuous accumulation of mineral dust from Marine Oxygen Isotope Stage (MIS) 11-5a (c. 430-60 ka). The outcome of our integrative approach indicates a continuous aridification over the last four interglacial/glacial cycles and we discuss potential changes in seasonality over time.

INTRODUCTION
Over the past decades, many loess-paleosol sequences (LPSs) have been investigated predominantly in the extensive northern hemisphere loess belt (Marković et al., 2015;Schaetzl et al., 2018;Lehmkuhl et al., 2021). Unlike lake and other terrestrial records, which are rather sparsely distributed, the spatial continuity of LPSs makes them valuable archives of past environmental change in the prevailing climatic past regimes that sustained loess formation and preservation (e.g., Basarin et al., 2014;Marković et al., 2015). Glacial and interglacial cycles, usually associated with shifts in humidity and temperature resulting from longterm variations in orbital parameters of eccentricity, obliquity, and precession driving ice volume (Imbrie and Imbrie, 1980;Heslop et al., 2000;Lisiecki and Raymo, 2005;Sun et al., 2006;Abe-Ouchi et al., 2013), are recorded through the alternation of loess and paleosol horizons, respectively. The paleoclimatic relevance of the quasi-cyclic alternation of loess and paleosols in LPSs of the Chinese Loess Plateau (CLP) was demonstrated in the 1980s through the correlation of magnetostratigraphically dated susceptibility records to the marine oxygen isotope records Liu, 1982, 1986), and similarly in the 1990s based on long-term variability in grain size distribution (e.g., Ding et al., 1994). Most of the early work relied on analyses of rock-magnetic parameters, often limited to magnetic susceptibility (Ding et al., 1993) and its frequency dependence. However, quantitative analyses of loess color spectra are also valuable indicators of shifts in mineralogical assemblages, as well as in organic matter content (Ding et al., 2002b;Lukić et al., 2014). Studies of LPSs from the CLP (Ji et al., 2002) and the Danube Basin (Lukić et al., 2014;Obreht et al., 2016) demonstrated the strength of combining magnetic and colorimetric parameters. While many multi-proxy studies have been performed on archives covering the last glacial cycle, investigations of European LPSs spanning multiple glacialinterglacial cycles are fewer (e.g., Jordanova et al., 2007;Necula et al., 2013;Basarin et al., 2014;Marković et al., 2015;Zeeden et al., 2016;Sümegi et al., 2018;Antoine et al., 2019;Obreht et al., 2019). In this study, we extend colorimetric data to about 430 kyr with the Zemun LPS record from the Middle Danube Basin (MDB), spanning from Marine Isotope Stage (MIS) 11 to MIS 5a.
Besides a dominating amount of quartz, feldspar, phyllosilicate and carbonate grains, comprising the average composition of the upper continental crust, Eurasian loess consists of heavy minerals and measurable relevant ferromagnetic (s.l.) particles such as a broad variety of iron oxides (Maher, 2016). Rock magnetic investigations allow differentiating magnetic particles formed insitu via pedogenic processes from the initial detrital content of magnetic particles in wind-blown loess. Relevant iron oxides for deciphering between loess and paleosol-units, which are readily detectable and relative concentrations quantifiable by roomtemperature magnetic investigations, consists of magnetite, maghemite and hematite (Heller and Evans, 1995). Pedogenesis takes place under relatively warm and humid conditions (Maher, 2016), involving various abiotic and bio-mediated chemical reactions not yet completely understood (e.g., Torrent et al., 2007; sections 2.5.2 and 2.6 in Lagroix et al., 2016) and resulting in the neo-formation of ultra-fine magnetic minerals (Maher and Taylor, 1988;Dearing et al., 1996;Maher, 1998;Torrent et al., 2007;Hu et al., 2013). Pedogenetically neoformed magnetite/maghemite have magnetic grains sizes ranging from unstable single domain [superparamagnetic (SP)] to stable single domain (SD) sizes, where the SP/SD particle size threshold at room temperature is about 30 nm (e.g., Peters and Dekkers, 2003). SP particles display a frequency dependence of susceptibility, which can be used to quantify their relative concentration and help discriminate between paleosols (high SP concentration) and loess (low to SP concentration). Weathering of loess can also occur under glacial conditions (Maher, 2011), but it is generally less intense due to increased aeolian sediment accumulation rates and generally drier conditions (Kohfeld and Harrison, 2003). Such weak glacial pedogenic alterations can also be detected by rock magnetic investigations (e.g., Taylor et al., 2014).
In most Eurasian LPSs, magnetic susceptibility values of paleosols are enhanced with respect to loess units (e.g., Maher, 1998Maher, , 2016Marković et al., 2009Marković et al., , 2015 from which a first order (chrono)stratigraphy can be established. For example, Marković et al. (2015) proposed a composite loess stratigraphy nomenclature scheme for the Danube loess belt region to facilitate pan-Eurasian LPS comparisons based on trends in magnetic susceptibility and other chronostratigraphic constrains. The stratigraphic system is based on the "S" (for soil) and "L" (for loess) labeling without any specific regional prefixes as is well established for the Chinese loess stratotype sections (Kukla and An, 1989). In the present study, the proposed composite stratigraphy is used to compare magnetic susceptibility variations of the Zemun LPS, evolving under non-monsoonal controlled (paleo-) climate with that of a monsoonal-dominated loesspaleosol sequence from the CLP (Luochuan, e.g., Hao et al., 2012).
The Zemun site hosts a variety of Middle and Late Paleolithic, and Neolithic artifacts found on the river bank nearby the LPS outcrop and upstream to the northwest (Šarić, 2008). The stratigraphic position within the Zemun LPS of these artifacts is, to date, unknown. However, providing a chronological framework for the protected loess site and establishing a range of environmental conditions these settlements may have been subjected to is important and will be beneficial to future archeological work.
Establishing a reliable chronology for LPSs remains challenging when absolute dating techniques such as luminescence reach their dating limits. Correlative age models are a mean to partly overcome this challenge. For LPS, correlative age models rely on the key process of pedogenic alteration, which induces variations in physical parameters and consequently records the environmental impact of interglacial/glacial cycles. In this study, magnetic susceptibility measurements coupled to diffuse reflectance spectrometry (DRS) analyses provide insights into past changes in environmental conditions. Rock-magnetic parameters provide insight into the source of variations in magnetic susceptibility, and together trace changes in the strength of pedogenesis along the Zemun LPS. Lastly, the developed correlative age model provides a mean to compare DRS and magnetism-based indicators of pedogenesis with global climate stacks (here the LR04 stack; Lisiecki and Raymo, 2005).
Volcanic ashes are often identified in European loess units, but in most cases, these tephras are mixed within aeolian loess or occur as crypto-tephras invisible to the naked eye. As most widespread tephras stem from highly explosive silicic magmas they are rather dominated by volcanic glass instead of by minerals or rock fragments (e.g., Lowe et al., 2017). Therefore, tephra layers are expected to have magnetic grain sizes similar to rapidly cooled magmatic rocks, characterized by dominantly SD and pseudo-single domain (PSD) particles, representing "frozen" magma at the time of eruption (e.g., Till et al., 2011). Both the magnetic grain size and composition of magnetic minerals contained in the tephra are also expected to contrast with that of the loess in which they are embedded. In the Zemun LPS, two tephra layers were identified through distinguishing magnetic properties similar to previous reports (Marković et al., 2015. Despite the lack of geochemical and mineralogical constraints, they are tentatively assigned to the so-called L2 and Bag tephra, both widespread regional marker-horizons.

PROFILE SETTING AND STRATIGRAPHY
The Zemun LPS (N 44.9246 • , E 20.3197 • ) is located in the southern part of the Middle Danube Basin and is exposed by the escarpment of a landslide scar developed in the high cliffs on the right bank of the Danube at the northwestern border of the Belgrade conurbation (Šarić, 2008). The profile belongs to the Srem loess plateau, which is delimited to the East and Northeast by the Danube and to the South by the Sava (Figure 1A). The Zemun LPS reveals several loess-paleosol couplets, reflecting warm/humid (interglacial) and cool/dry (glacial) past environmental conditions. The Zemun LPS shows similar litho-and pedo-stratigraphic features as neighboring profiles at Batajnica (Marković et al., 2009) and within the Titel loess plateau (Marković et al., 2015).
Parallel to sampling, the lithology of the sequence was described in the field focusing on loess structure, pedogenic features, and sediment color. Additionally, sediment samples were investigated under the binocular microscope in order to estimate mineralogical compositions and degree of roundness of grains. In loess deposits, aeolian origin is assumed when single grains predominantly show little rounding due to direct source to sink transport, whereas well-rounded grains reflect complex sediment recycling, long term transport and/or fluvial processes. Soil units were described following the World reference base for soil resources (IUSS Working Group WRB, 2015 1 ), even though we are fully aware that this system is not designed to classify buried soils also referred to as paleosols. Additionally, we follow the concept of accretionary soil formation characteristic of aeolian landscapes of western Eurasian dry steppe regions Jordanova and Jordanova, 2020;Lehmkuhl et al., 2021). These observations combined with results of laboratory analyses (colorimetric and magnetic parameters) form the base for a genetic interpretation of sediment formation.
Sampling for sedimentological analysis started at the transition from a basal brownish-yellowish loess unit to the lowermost exposed interglacial pedocomplex (S4) characterized by humus infiltrations in ancient root channels, carbonate concretions (up to 5 cm in diameter) and hydromorphic features, such iron coatings and iron/manganese oxide patches ( Figure 1C). The basal red-brownish fossil pedocomplex (24.35-22.30 m relative profile depth) is c. 2.05 m thick and can be classified as a Cambisol typical of forest-steppe environments. A gradual transition is observed to the overlying loess unit (L4, 22.25-21.90 m). From 21.85 to 21.70 m, a 0.20 m thick bed of loosely cemented loess can be recognized. At 21.65 m depth, a drastic change in sediment toward sandier material is observed and referred hereafter as marshy floodplain sediment (MFS) (see section "Integrated Stratigraphy"). The MFS persists from 21.65 to 16.50 m with alternation of sandy iron-stained and sandy-silty beds with dark humus rich partly clayey beds dominating.
The subsequent polygenetic light reddish-brown paleosol (S3, 16.45-18.80 m) can be classified as a Phaeozem typical of steppe environments with mild and relatively humid winters. It is built from the bottom to the top by a lower pedogenetically overprinted unit. This is followed by a thin, lightly colored horizon with carbonate nodules and numerous krotovinas likely representing a pedogenetically overprinted loess unit overlain by a pedogenetic horizon with again many krotovinas. A gradual transition to the overlying c. 2.65 m thick yellowish-grayish loess (L3, 13.75-11.10 m) is observed.
The contact between L3 and the pedocomplex above is massively bioturbated and interspersed with carbonate concretions of sub-centimeter to decimeter size. This interglacialtype pedocomplex is c. 3.50 m thick and comprises in its lower part a brownish-blackish to dark grayish c. 2.00 m thick paleosol unit (S2, 11.05-8.70 m), which is classified as a polygenetic steppe soil of Kastenozem to Chernozem type typical of cool mid-latitude steppe environments. The lower strongly developed darker horizon has many carbonate pseudomycelia. Krotovinas are scattered throughout the middle lighter colored horizon and the upper weakly developed interstadial-type paleosol horizon (L2SS1,. S2 and L2SS1 are separated a thin (∼25 cm) pale yellowish loess unit (L2LL1, 8.65-8.45 m), which is bioturbated in its uppermost part transitioning into L2SS1 interstadial paleosol. The overlying pale-yellowish to partly yellowish-grayish loess unit (L2, 7.30-3.60 m) is c. 4.00 m thick and shows numerous fine lamina consisting of sandy loess beds in its middle part. At c. 5.8 m, a 2-3 cm reddish to brownish loose loess layer was observed.
The contact between L2 and the uppermost interglacial pedocomplex (S1, 3.55-0.60 m) is also characterized by massive carbonate concretions and krotovinas. The S1 pedocomplex is c. 2.95 m thick, composed of a lower dark brownish-grayish paleosol overlain by lighter paleosols marking variable degrees of pedogenesis during continued dust deposition. The entire complex can be characterized as a typical Chernozem soil with decreasing intensity of pedogenesis with time. The uppermost horizon of the S1 pedocomplex is rather weakly developed and contains numerous krotovinas. The uppermost layer of the studied Zemun LPS is composed of pale porous sandy loess (L1, 0.60-0.00 m), which is loosely cemented with a few thin laminated beds of sandy-silt.
With the exception of the MFS sand dominated unit between 21.65 and 16.50 m, the entire sequence consists of coarse primary eolian silt with alternating contributions of fine sand in the pure loess and medium to fine silt in the paleosol units. Similar stratigraphic features can be observed in LPSs close to the Zemun LPS, e.g., from Batajnica (Marković et al., 2009), and the Titel-Stari Slankamen composite profile (Marković et al., 2011(Marković et al., , 2015.

Sampling
Sampling was performed from abseiling the wall (subsections V, A) and in dug trenches (subsections B, C, D, E, P 1-3). Subsections V and A cover the uppermost c. 15 m and are c. 2 m apart, whereas subsections B, C, D, E, and P 1-3 are located between c. 6 and 10 m upstream to the northwest ( Figure 1D). All subsections overlap and provide the basis for the construction of the composite profile. After thorough cleaning of the profile, bulk samples were taken with 5 cm spacing for a total of 479 sampling depths. For laboratory analyses described in following section "Room Temperature Susceptibility Measurements", bulk samples were dried at 30 • C for 48 h, homogenized and compressed into non-magnetic 6.4 cm 3 plastic boxes (n = 479) and ∼1.1 cm 3 gelatin capsules (n = 262). Masses of bulk samples within boxes and capsules were measured and recorded.

Room Temperature Susceptibility Measurements
Magnetic susceptibility measurements were carried out in 2017 at the Environmental and Palaeomagnetic laboratory at the University of Bayreuth. Measurements of low-field mass normalized magnetic susceptibility (χ), performed on all 479 box-samples, were obtained with a Magnon kappa bridge (Magnon, Dassel, Germany-VFSM) operating with a 320 A/m applied field at low frequency (300 Hz, χ lf ) and high frequency (3,000 Hz, χ hf ). The absolute frequency dependence of magnetic susceptibility χ is determined from: and the percent increase (χ fd ) with respect to χ lf from: (Mullins and Tite, 1973;Dearing et al., 1996;Eyre, 1997).

Hysteresis Measurements
Hysteresis measurements were performed with a Princeton Measurements Corporation Model 3900 Vibrating Sample Magnetometer (VSM) at the Institut de Physique du Globe de Paris (Paleomagnetism Research Group). A total of 262 samples were investigated across sub-sections of the Zemun profile covering the interval from 12.05 m depth (L3) to the top of the profile continuously (n = 240), an interval from 22.25 to 21.50 m continuously (n = 16) and 6 other discontinuous intervals within the sandy MFS (Figure 1). For 6 of 10 samples stemming from the MFS sample preparation required modification to prevent grain movement favored by the higher sand content during hysteresis, FORC and IRM acquisition measurements on the VSM. The preparation consisted in piercing a hole at the base of the capsule to allow air to be released during the filling process. A thin layer of diamagnetic cotton was inserted to prevent single grains to fall through the air-release hole. Between 120 and 170 mg of the coarse grain sediment was added to the capsule and the netweight measured. Grains were immobilized by adding neoprene gel glue, verified to be diamagnetic (c. −3 × 10 −8 m 3 kg −1 ). Complete saturation with glue was confirmed by leaking glue through the air-release hole. The lid was also filled with neoprene glue before sealing the capsules with the temperature resistant tape (see Supplementary Figure 1). Samples were dried overnight.
Prior to hysteresis measurements, the χ of each sedimentfilled capsule was determined with an AGICO KLY-3 kappabridge operating in a 300 A/m applied field at a frequency of 875 Hz. These low-field susceptibility measurements were used to (1) test whether the gelatine capsule subsamples are representative of the larger box samples, and (2) calculate the ferrimagnetic component of susceptibility (χ ferri ) by calculating the difference of the KLY3 bulk low-field susceptibility and the calculated highfield susceptibility (χ hifi ) from the slope of the high-field linear segment of the hysteresis loop: Three experiments were conducted with the VSM. (1) Hysteresis loops were measured in maximum applied fields of ± 1.5 T over 100 ms measurement times at each applied field step. The applied field increment was set to 5 mT. Experiments began at positive maximum field following a 1 s pause. The following parameters were derived from each hysteresis loop: coercivity (Hc), saturation magnetization (Ms), saturation remanent magnetization (Mrs) and high-field magnetic susceptibility (χ hifi ), calculated from the linear high-field slope of the magnetization above 1.05 T. The ratio between χ ferri and Ms can be used to track variations in relative concentration of superparamagnetic (SP) particles assuming a constant magnetic mineral assemblage. (2) The coercivity of remanence (Hcr) was determined from backfield direct current demagnetization of the forward maximum field (1.5 T) isothermal remanent magnetization (IRM). 70 backfield steps were logarithmically spaced over the 0-500 mT range. (3) Backfield direct current demagnetization of the forward maximum field (1.5 T) IRM was also measured using linear 100 mT backfield steps from 0 to 1.5 T. From these data, relative contributions to the total IRM 1.5T , such as the S-ratio are determined, as well as absolute contributions to the total IRM 1.5T within various coercivity windows, such as HIRM, are quantified. The S-ratio was calculated after King and Channell (1991) as follows: resulting in values possibly ranging from −1 to +1, where values decrease from 1 as the relative contribution of hard magnetic minerals (defined here as having a coercivity of remanence greater than 300 mT) to the total IRM 1.5T increases. The hard isothermal remanent magnetization (HIRM) was calculated from the difference in magnetization between the 0.3 T and 1.5 T backfield steps following a forward field magnetization in 1.5 T (e.g., Taylor and Lagroix, 2015;Liu et al., 2016): The ratios of Hcr/Hc and Mrs/Ms (indicative of a samples mean magnetic grain size) were plotted on a modified Day-Dunlop plot (Day et al., 1977;Dunlop, 2002a,b), tracing differences in loess/paleosol samples and potentially tephrabearing samples, the latter expected to contain higher PSD and SD content than loess.

First Order Reversal Curve (FORC) Investigations
FORC investigations were carried out for six samples using the VSM. Specifically, FORCs were acquired on two samples suspected to contain tephra material (ZV 117, 5.80 m profile depth, containing the L2 tephra) and ZP1 001 (21.70 m, containing the Bag tephra, and reflecting the highest magnetic susceptibility signal out of 4 samples (ZP1 001-ZP1 004). To identify differences in magnetic grain size and possible differences in magnetic interactions, FORCs were also carried out on samples representative of loess units bracketing the suspected L2 and Bag tephra occurrences ZV 112 (5.55 m), ZV 126 (6.25 m), ZE 047 (21.65 m), and ZP1 016 (22.45 m). ZE 047 marks the base of the MFS and is hydromorphically overprinted but it is the sample immediately above the suspected Bag tephra sample analyzed (ZP1 001). Samples were exposed to a saturating field of 500 mT, Hu (min, max) were selected after initial test-measurements to ± 80 mT, and Hc (min, max) to 0 and 100 mT. Averaging time was set to 300 ms, reducing measurement noise. For each diagram, 150 FORCs were acquired at a field increment of ca. 1.9 mT. For weak samples originating from loess units, 2-9 series of FORC measurements were measured and averaged. The number of FORCs averaged for each sample is reported as n-values in Figure 6. FORC data analyses were carried out using FORCinel version 3.06 (Harrison and Feinberg, 2008). The main pre-processing consisted in drift and high-field slope corrections, first point removal, and lower-branch subtraction. FORCs of loess samples, for which multiple FORC were measured, were averaged in FORCinel. FORC diagrams were improved by color rescaling selecting a rectangle area around the central ridge and using the "Autoscale" option. Smoothing was conducted according to the VARIFORC approach (Egli, 2013) with the following smoothing factors: vertical ridge Sc0 = 4, central ridge Sb0 = 3, horizontal smoothing Sc1 = 7, vertical smoothing Sb1 = 7, horizontal and vertical lambdas 0.1, output grid = 1 and central ridge offset = 0. Horizontal profiles were extracted using FORCinextras.

High-Resolution IRM Acquisition
Complementary to FORC measurements, high-resolution stepwise isothermal remanent magnetizations (IRM) were acquired on the same samples using the VSM. Samples were demagnetized before IRM acquisition. The initial field was set to 15 µT, the final field to 1.5 T, using 300 logarithmically spaced measurement steps.
Tentative un-mixing of the IRM acquisition curves was carried out with aid of the MAX UnMix RShiny web-application (Maxbauer et al., 2016). Smoothing factors were set between 0.5 and 0.6, depending on measurement noise. Before using MAX UnMix, first derivatives of the raw IRM acquisition data were calculated, distributions of suspected tephra and bracketing loess samples compared to determine roughly the number and characteristics of geologically realistic components. These components were used in the fitting panel of MAX Unmix, and then refined in the optimization panel by minimizing the residual sum of squares (RSS). Final determination of the components characteristics were obtained in the Error Analysis panel, providing mean values and uncertainties based on 100 Monte-Carlo random re-samplings of 95% of the data.

Colorimetric Measurements
Initial colorimetric measurements were carried out at the Environmental and Palaeomagnetic laboratory at the University of Bayreuth. Colorimetric analyses have been done in the past decades by visual identification of Munsell color charts (e.g., Tsatskin et al., 1998;Günster et al., 2001;Machalett et al., 2006), but its importance for paleoclimatic reconstruction remained little explored for a long time. These color charts provide alphanumerical indices for each sample by comparing the sampling material with a single-color rectangle. This procedure has some disadvantages (Post et al., 2015) such as subjective evaluation by a researcher and dependence on specific moisture and lighting. Spectrophotometers eliminate these uncertainties (Sun et al., 2011) and provide objective numerical values of luminance (L * ), redness (a * ) and blueness (b * ). L * ranges from 0 (black) to 100 (white), whereas a * ranges from negative (green) to positive (red) and b * from negative (blue) and positive (yellow). The redness a * is commonly interpreted as reflecting weathering intensity (Yang and Ding, 2003). These three parameters span a 3-dimensional color sphere and derived Lab-values can be transformed to RGB colors to make their real colors visible (as conducted in Figure 3). Furthermore, spectrophotometers provide backscattered intensities of wavelengths each 10 nm from the visible light spectrum (here 400 to 700 nm only). These backscattered spectra are capable of tracing distinct minerals out of a heterogeneous bulk-sample mineral assemblage. An additional advantage is the fast (∼2 s per sample) and nondestructive nature of the measurement.
Diffuse Reflective Spectroscopy (DRS) measurements were carried out in 2020 at the Leibniz Institute for Applied Geophysics (LIAG, Grubenhagen) with a Konica Minolta CM 700d spectrophotometer, in 10 • observer angle and using the D65 norm-light calibration. An oculus width was obtained with a 0.8 cm open adapter. Backscattered reflectance spectra were converted into 1st derivative values with a newly developed R-script. 1st derivative values were used to calculate the hematite/goethite ratio (HGR), following the same approach described in Wu et al. (2018). HGR was calculated from the ratio of backscattered intensities (I) at 565 nm, associated to hematite, and 435 nm, associated to goethite, so that HGR = I 565 nm /I 435 nm (Barranco et al., 1989;Deaton and Balsam, 1991;Debret et al., 2011). The backscattered intensity bands were selected based on peaks in the 1st derivative spectrum (Figure 2). For the transformation of Lab values into RGB colors, a modified R-script (R Core Team, 2020) was used .

Generation of the Age Model
The age model (Figure 9) was constructed by correlating χ variations and the LR04 stack (Lisiecki and Raymo, 2005), reflecting mainly global ice volume by benthic δ 18 O. Tiepoints were selected based on (1) variations in magnetic parameters interpreted to fluctuate predominantly as a function of pedogenic intensity, and (2) similarities in fluctuation and amplitude of magnetic parameters when compared to the LR04 stack and Imbrie and Imbrie ice model (Imbrie and Imbrie, 1980). Identified tephra layers provide additional tie-points. This approach resulted in 13 tie-points (Table 1). Two additional tie-points constrain the age model stemming from the two potential tephra layers identified magnetically and compatible with sedimentological observations in the field. For the age model construction, the "astrochron" package version 0.9 (Meyers, 2014) was used for linear interpolation.

Magnetic Susceptibility Parameters and Stratigraphical Assignment
Based on the record of magnetic susceptibility with depth, a general pattern with low-frequency susceptibility (χ lf ) values being generally low in loess units, and high in paleosols can be observed for the Zemun LPS (Figure 3 and Table 2). Such a pattern is well known from Eurasian LPS and predominantly reflects past environmental changes between glacials and interglacials during the Pleistocene. Following the stratigraphic nomenclature introduced by Marković et al. (2015), results reveal that the sampled profile extends from a basal S4 pedo-complex, assigned to MIS 11 to the oldest part of the L1 directly overlaying the last interglacial-early glacial pedo-complex S1. Beside distinct elevated magnetic susceptibilities in interglacial paleosols, the highest χ lf value (157 × 10 −8 m 3 kg −1 ) is observed as a sharp peak at 5.80 m depth within the L2 unit. This peak coincides with the field description of a 2-3 cm reddish to brownish loose loess layer and is stratigraphically and magnetically consistent with the L2 tephra layer described in other LPS regionally (e.g., Laag et al., 2018;Antoine et al., 2019). Mean χ lf values of paleosols S3 to S1 decrease slightly from ∼115 to 100 to 90 × 10 −8 m 3 kg −1 , FIGURE 3 | Stratigraphic view of magnetic susceptibility parameters χ lf , χ, χ fd [%], and selected colorimetric parameters L* and a*. The variability of magnetic susceptibility parameters with depth reflects the pacing of interglacial (S4-S1) and glacial (L4-L1) cycles. Potential tephra layers (red dotted lines) are referred to as the L2 tephra and to as the Bag tephra. Note that the scale for L* is inversed. Lab color values are transformed to RGB colors for each sample and plotted in the background.  Figure 3). Mean values of χ also decrease from S3 to S1. The S4 pedocomplex displays mean values of χ lf , χ and χ fd % that are much lower than S3-S1. Finally the MFS layer has moderately high χ lf values, greater than S4 but less than S3-S1, which fluctuate over short depth intervals. χ and χ fd % values of the MFS are low comparable to loess units. The magnetic enhancement of the paleosol is readily observable on the χ vs. χ lf bi-plot shown in Figure 5A. Depth intervals within the MFS and containing tephra material plot above the main trend line. Two populations are distinctively identified along this trend with one predominantly represented by samples showing no magnetic enhancement characterized by nearly all samples within loess units, and a second reflecting pedogenesis (higher χ and χ lf ) are all stemming from paleosol units. The strongest enhancement is expressed in paleosol S3, while S1 and S4 pedocomplexes spread over a wide range of χ and χ lf . A background low-frequency magnetic susceptibility of 24.6 × 10 −8 m 3 kg −1 is obtained by fitting a linear regression line through the data excluding the MFS and tephra intervals. This background value is higher than what is observed in Central Asia and China (17 ± 2 × 10 −8 m 3 kg −1 Forster et al., 1994) or at Semlac in Romania (10 × 10 −8 m 3 kg −1 Zeeden et al., 2016) but similar to LPS north of Prague in the Czech Republic (22 × 10 −8 m 3 kg −1 for Sedlec and 23 × 10 −8 m 3 kg −1 for Zemechy Forster et al., 1996).

Hysteresis Parameters
Hysteresis loop derived parameters are summarized in Table 3 and plotted in Figures 4, 5. Combining data sets acquired on  Figure 2) by KLY3 χ lf vs. VFSM χ lf measurements defining a linear fit with an r of 0.99. The commonly used magnetic proxy for pedogenic intensity is χ, which increases with increasing degree of pedogenesis is compared to χ ferri /Ms, Hc, HIRM and S-ratio in Figure 4A through D, respectively. The χ ferri /Ms ratio, like χ, track relative changes in SP particle concentration, assuming the mineral assemblage contributing is monomineralic or that proportions between mineral components are constant. The data acquired clearly shows that the magnetic mineral assemblage is not monomineralic but that the soft ferrimagnetic component dominates both Ms (Figure 5C) and χ ferri (Figure 5G). The magnetic mineral assemblage of the MFS and tephra samples do differ in terms of composition, concentration and magnetic grain size with respect to the loess and paleosol units, thus falling off the main trend line correlating (r = 0.9) χ ferri /Ms and χ. At higher pedogenic intensities, χ ferri /Ms increases less rapidly than χ suggesting an increased contribution to Ms of a mineral, likely hematite, with a lower saturation magnetization than magnetite and maghemite. Figure 5I shows HIRM increasing as χ increases in paleosols corroborating the previous observation. In Figure 4B, χ and Hc are displayed, showing for all samples (tephra layers excluded) a negative exponential behavior. Loess and MFS samples reflect the highest coercivities with low amounts of pedogenetically formed SP particles. The tephra samples reflect a diametral behavior, indicating the presence of SP particles but at similar SP concentrations the tephra samples have a higher bulk coercivity due to either a higher proportion of SD and PSD grains and/or of high coercivity minerals. Figure 4C confirms the latter and results presented in section "FORC Analysis" supports the former. HIRM and χ correlate positively with a weaker coefficient (r = 0.69) than for soft remanencebearing ferrimagnetic (IRM −100 mT ) component ( Figure 5F, r = 0.92), indicating that high-coercivity minerals present are also in paleosols and their concentration increases with increasing pedogenic intensity ( Figure 4C). Interestingly, samples from the MFS, L2SS1 and samples stemming from the L4 loess unit also define a linear trend but with a much higher slope where increases in HIRM are associated with only small increases of χ. Figure 4D correlates χ and the S-ratio. S-ratio values equal to 1 indicate the presence of only soft magnetic minerals like maghemite and magnetite, whereas reduced values underline the presence of high-coercivity minerals like hematite and goethite. The relative proportion of hard-magnetic minerals contributing to the IRM 1.5T is greater in Zemun loess, especially L2, than in paleosols units. The highest S-ratio values are observed for samples stemming from the MFS and the tephra samples. Excluding MFS and tephra samples, the S-ratio correlates positively with χ with r = 0.79 ( Figure 4D) indicating that the relative proportion of soft-magnetic minerals contributing to IRM 1.5T increases with increasing degree of pedogenesis. At the same time, absolute concentrations of both soft remanencebearing ( Figure 5F) and hard ( Figure 4C) magnetic mineral components increase individually with increasing pedogenesis.
Finally, hysteresis data presented on a modified Day-Dunlop plot (Supplementary Figure 3) reveals that all loess and paleosol samples fall within the PSD range as mean domain size, but samples stemming from the MFS as well as the tephra samples indicate a higher amount of SD particles. FORC diagrams and central ridge profiles of the bracketing loess for each subset interval, are very similar (Figures 6A,C,D,F). Samples ZV112 and ZV126 display a FORC diagram globally expected for PSD assemblages, in accordance with hysteresis parameters (Hcr/Hc = 3.18 and 3.29, Mrs/Ms = 0.14 and 0.13). The horizontal profiles are characterized by one main coercivity peak around 10 mT, and a second peak close to Hc = 0, along with a broad distribution up to about 80-100 mT. The PSD character of the FORC diagrams of samples ZE047 and ZP1016 is less pronounced and also reflected by the main hysteresis parameters (Hcr/Hc = 2.62 and 2.72, Mrs/Ms = 0.17 and 0.17). For these samples, the horizontal profiles are also characterized by a peak around 10 mT, but lack the very low-coercivity peak, probably due to a lesser content in grains either close to the superparamagnetic threshold or in the multidomain range. Sample ZV117 (Hcr/Hc = 2.46, Mrs/Ms = 0.17) displays a FORC diagram with a main distribution more elongated along the central ridge axis and shifted to higher Hc values, making the diagram appear more SD-like than those of the bracketing loess samples (Figure 6B). A similar observation can be made from sample ZP1001 (Hcr/Hc = 2.66, Mrs/Ms = 0.16), although it is not as prominent as for sample ZV117 (compare Figures 6E,B). This is also seen in the horizontal profiles, where the main distribution peak are broadened and shifted to slightly higher values by 5-10 mT. This indicates the input of additional and non-loessic material.

High Resolution IRM Acquisition
Detailed IRM acquisitions can provide alternative insights on the coercivity distributions because focusing on remanent magnetizations instead of direct in-field measurements, which is the case for FORCs. These are therefore differently affected by magnetic particle sizes (for instance, superparamagnetic particles do not carry a remanent magnetization). IRM acquisition curves were obtained for the same six samples that underwent FORC  Table 3. For correlation coefficients and regression fits, samples from the MFS, and the tephra layers were excluded. as shown in Figures 7A,F. All six samples fail to saturated by 1.5T, indicating the presence of high-coercivity minerals. Raw coercivity distributions of samples ZV112 and ZV126 (Figure 7G) are nearly identical, as expected from both hysteresis parameters and FORC diagrams. Samples ZE047 and ZP1016 are also very similar ( Figure 7G). The coercivity distribution of all four bracketing loess samples is characterized by two main peaks at about 20-30 mT (1.3-1.5 in log scale) and 60-100 mT (1.8-2.0 in log scale) (Figures 7C,E,H,J). The distribution of sample ZV117 (Figure 7D) differs significantly from that of samples ZV112 and ZV126 and is dominated by a peak around 45 mT (1.65 in log scale). For sample ZP1001, the distinction with its bracketing loess samples is less marked (Figure 7I) Colorimetric/DRS Analysis L * values are generally higher in light loess units meaning brighter colors due to higher amounts of quartz and carbonates than in paleosols, where clays and organic matter decrease the luminance (Figure 3 and Table 4). The highest L * value (72.5) is found in the L2SS1 and the lowest L * value (50.8) in the S3 pedocomplex. Luminance values follow stratigraphic units from the S3 pedocomplex to the L1 loess at the top of the profile. Below S3, the correlation of L * with the stratigraphy is more complex. The pedocomplex S4 is weakly expressed in the Zemun LPS characterized by a broad range of L * values. S4 can be considered as brighter than the other pedocomplexes with brightest values in the lower part at the base of the profile, possibly as a result of carbonate precipitation.
Similar to luminance (L * ) values, redness values (a * ) also reflect alternations between loess and paleosols (Figure 3). However, a clear separation between loess and paleosols is not as obvious for all paleosols from a * data. In the L1 loess unit, a * values vary between 3.2 and 4.2. Compared to the mean a * from the S1 pedocomplex, redness is reduced in L1 (3.8, in S1 4.4). The mean a * of L3 is similar to L2 (4.2) and relatively low with respect to S1. Higher a * values are reached in the S3 pedocomplex (mean = 5.7), but the highest a * values are reached in the S4 pedocomplex with a mean of 5.8 and a maximum of 6.6. Over the entire profile, a decreasing trend from S4 to S1 samples can be observed with even lower a * values for the weakly developed interstadial L2SS1 paleosol. Blueness (b * ) values, which are inverse to yellow, are generally low in paleosols and high in loess units. The lowest b * values are found in the pedocomplex S1 and the highest b * values in the S4 pedocomplex. However, the pedocomplex S2 (compared to the other pedocomplexes) does not show these well-expressed low values of b * . In summary, a * , a commonly used pedogenic indicator suggest the highest intensity of pedogenesis is found in S4 which contradicts interpretations based on magnetic proxies and field observations. DRS derived goethite-and hematite concentrations (1st derivatives of intensities at 435 and 565 nm; I 435 and I 565 ) and the ratio of hematite to goethite (HGR) mainly follow glacial/interglacial cyclicity. Both I 435 nm and I 565 nm values are reduced in paleosol complexes and are high in loess units with the exception of S4, which has values characteristic of the loess units (Table 4). HGR values are maximum in paleosols and like magnetic susceptibility parameters decrease from S3 to S1. Luminance (L * ) and redness (a * ) correlate with r ∼ −0.3, while L * and blueness (b * ) correlate positively with r∼0.78. HGR correlates with L * with r∼-0.82, with a * r∼0.76, and b * with r∼-0.5.
Additionally, selected colorimetric parameters are tested against magnetic indicators of relative concentration of SPparticles ( χ and χ ferri /Ms) and of high-coercive minerals (Hcr and S-ratio) (Figure 8). In Figures 8A-C     correlation of χ with HGR (r = 0.87) reflects a relatively higher amount of hematite than goethite in paleosol samples. This is supported by correlations based on magnetic indicators for hard magnetic minerals like hematite and goethite (high Hcr and low S-ratio values) for loess samples reflecting positive correlations for hematite and goethite and a negative correlation with HGR. This underlines the presence of hematite and goethite in loess samples and a relatively higher presence of hematite compared to goethite in paleosol samples. Overall, the variation of HGR and I 565 is high in paleosols and low in loess units. I 435 varies distinctively less in loess and paleosol samples. The luminance L * reflects toward younger samples a continuous evolution of brighter paleosols, where the S3 is the darkest (Figure 3). a * decreases as well, even more than L * , for all paleosols from S4 to S1 (Figure 3). With respect to the hematite/goethite ratio (HGR), the decreasing trend is most obvious. While there is a general increase of the colorimetric properties with increasing χ, these are clearly non-linear.

Age Model
The age model (Figure 9) was constructed based on 15 tiepoints listed in Table 1 and shown as dotted lines in Figure 10.
Thirteen tie-points are based on synchronizing variations in χ and benthic δ 18 O values from the LR04 stack (Lisiecki and Raymo, 2005). In total, fifteen tie-points were used including the two identified tephra layers (L2 and Bag tephras). Two additional tie-points correspond to tephra layers demonstrated to be present based on magnetic properties presented in sections "Magnetic Susceptibility Parameters and Stratigraphical Assignment" through "High Resolution IRM Acquisition", especially evidenced by FORC analysis and IRM unmixing, and compatible with field-based sedimentological observations described in section "Profile Setting and Stratigraphy".
The L2 tephra is a widespread tephra occurrence within southeastern European loess records (Laag et al., 2018;Avram et al., 2020). Based on its stratigraphic position in several other well-dated profiles, like the neighboring Batajnica LPS (Avram et al., 2020), and the Harletz LPS in Bulgaria (Antoine et al., 2019), the L2 tephra can potentially be correlated to the Vico-Ignimbrite B eruption, dated to 160.6 ± 2.0 ka . However, uncertainty remains with the age assignment, because (1) no geochemical evidence is to our knowledge available to trace back the L2 tephra to a specific eruption and (2) the archives of Lake Ohrid  and Fucino Basin (Leicher et al., 2016;Giaccio et al., 2017Giaccio et al., , 2019 exhibit several prominent tephra layers during MIS 6, which may be assigned to the L2 tephra with ages spanning 150-168 ka. For the present study, we anchored the L2 tephra at 160.6 ka (Table 1) following the hypothesis of a correlation with the Vico Ignimbrite B eruption and the age provided by Mannella et al. (2019).
The Bag tephra, demonstrated to be contained in the Zemun LPS in 4 samples between 21.85 and 21.70 m depth, has been FIGURE 9 | Age model based on the correlation of absolute frequency dependent magnetic susceptibility values ( χ) with LR04 (Lisiecki and Raymo, 2005) and Imbrie and Imbrie Ice model (Imbrie and Imbrie, 1980). (A) χ as function of stratigraphy. (B) Extract of the LR04 stack (last 430 ka). (C) Imbrie and Imbrie ice model (last 430 ka). (D) The age-tuned χ record based on the 15 tie-points, indicated by dashed lines (see Table 1). associated elsewhere to an eruption originating from the Alban Hill Volcanic region and dated to 350-360 ka (Pouclet et al., 1999;Marković et al., 2015). A corresponding tephra layer can be found in the Fucino Basin record (TF-85), which is dated to 367 ± 1.6 ka (Marra et al., 2009(Marra et al., , 2019Giaccio et al., 2012) and thought to originate from the Colli Albani Villa Senni eruption (Marković et al., 2015). We tentatively use this age estimate for the tie-point represented by the Bag tephra at Zemun (Table 1 and Figure 9).
Attempts to extrapolated ages for samples depths above and below the youngest and oldest tie-points lead to an age of 50 ka corresponding to MIS 3 interglacial, which is not contained in the Zemun LPS and therefore too young for the top of the profile. Based on magnetic susceptibility variations, the loess at the top of the profile was likely deposited during the early MIS 4 (∼60 ka).
The uncertainty in age extrapolation at the base of the profile is highlighted by the gray shading in Figure 9.

Integrated Stratigraphy
The Zemun LPS comprises four interglacial paleosol complexes, S4 (Cambisol), S3 (forest-steppe soil), and S2 and S1, both considered as steppe soils (Marković et al., 2011;Obreht et al., 2016). Our correlative age model (Figure 9) suggests that the Zemun loess profile covers the time interval from MIS 11 to MIS 5/4, which is additionally constrained by the identification of two widespread tephra layers of known stratigraphic position within FIGURE 10 | Magnetic susceptibility parameters χ lf , χ, χ fd [%], and the hematite/goethite ratio (HGR) as function of age. Note that the axes for L* and b* are inverted. The gray dotted lines indicate the tie-point-ages from Table 1. other Danube loess sequences (Marković et al., 2015. However, the selected correlative tie-points (Table 1) include the tephra layers, but do not account for the colorimetric data. The general stratigraphy of loess and intercalated paleosols at the Zemun LPS is in good agreement with independently dated records in the southern MDB, based on luminescence dating and paleomagnetism (Marković et al., 2009(Marković et al., , 2011(Marković et al., , 2015Basarin et al., 2014;Song et al., 2018;Avram et al., 2020). It is noteworthy, that uncertainties of the direct correlation of χ to the LR04 stack exist , mainly because of uncertainty in the exact position of tie points and correlation targets. The stacked LR04 record provides ages for the last 1 Ma with 4 kyr uncertainty (Lisiecki and Raymo, 2005).
The fluvial sediment facies (MFS unit) below S3 is widespread in the southern MDB and is interpreted as an intercalation of marshy floodplain sediments (MFS) consisting of an alternation of sandy iron-stained and sandy-silty beds with dark humus rich partly clayey beds (e.g., Marković-Marjanović, 1970;Gaudenyi et al., 2015). The MFS unit represents overbank deposits of a river system flowing through a steppe/forest steppe landscape laterally interfingering with aeolian loess. They probably represent relatively short-term intervals of aquatic deposits formed after prominent flooding events. The loess unit L4 in which the MFS are intercalated measures ∼ 1.1 m at the neighboring section Batajnica (Marković et al., 2009) and is generally relatively thin in the entire MDB (Marković et al., 2015) and across Eurasia (Song et al., 2018).
Magnetic susceptibility parameters of the Zemun LPS convincingly track the environmental impact of interglacials and glacials during the Middle to Late Pleistocene (Figure 9). All four paleosols (S4-S1) are characterized by elevated values of χ lf , χ, and χ fd indicative of higher concentration of magnetic minerals and specifically superparamagnetic particles formed insitu during pedogenesis (Maher and Taylor, 1988). Based on χ lf , χ and χ fd parameters, the intensity of pedogenesis increases from S4 (MIS 11) to S1 (MIS 5), to S2 (MIS 7) to S3 (MIS 9).
Colorimetric parameters are used to quantify hematite and goethite (e.g., Scheinost, 1998). The assumption that the HGR reflects changes in goethite remains questionable. Studies from Liu et al. (2008) and Jiang et al. (2018) show that goethite does not play a crucial role for LPS and additionally the process of the neo-formation of goethite is different from processes related to pedogenesis. In contrary, Bilardello et al. (2020) demonstrate in their experimental simulations that a simultaneous neoformation of goethite, hematite and magnetite is possible, and they add further evidence that biogeochemical conditions are a keycontroller of the proportions of the neo-formed minerals. Our good correlation of magnetic parameters (Figure 8) like the S-ratio and remanence of coercivity (indicating the presence of hard magnetic minerals) with 1st derivative intensity band values for goethite (I 435 nm ) indicate the presence of goethite. However, the good correlation of HGR with the pedogenesis indicators (χ ferri /Ms and χ) indicates an overall higher presence of hematite over goethite in the paleosols of the Zemun LPS. This said, the magnetic experiments conducted in the present study do not readily permit the identification of goethite, which is facilitated in low-temperature magnetic experiments (e.g., Lagroix and Guyodo, 2017). Following that, we consider for the Zemun LPS, that the HGR can reflect changes in goethite and/or hematite. Our quantitative colorimetric data show variations that are related to the different soil types. S3, the forest-steppe soil and the most mature paleosol based on magnetic susceptibility parameters and χ ferri /M s has the lowest luminance, which would result from high amounts of mineralized organic matter. With decreasing precipitation, the predominant soils are chernozems, which are represented by S1 and S2. Cambisols are brighter, which is reflected by S4 and its increased luminance (L * ).

Paleoclimatic Implications for the Zemun LPS
The stratigraphic interpretations of both rock magnetic and colorimetric data in combination with the correlative age model provide insights into the environmental evolution of the southern Middle Danube Basin for the last 430 kyr.
This would suggest that for S4 the intensity of pedogenesis was weaker compared to S3-S1, which is not supported by the general pedogenetic trend in the MDB and Western Eurasia (e.g., Marković et al., 2015, and references therein). Colorimetric data of redness (a * ), an indicator for weathering intensity (Yang and Ding, 2003), leads to the opposite conclusion. The a * parameter is highest in S4 paleosol (Table 4) suggesting it has the highest intensity of pedogenesis, which is in line with the general paleoenvironmental evolution in the Eurasian loess belt (e.g., Buggle et al., 2013Buggle et al., , 2014 and our pedologic interpretation. Both I 435 nm and I 565 nm 1st derivative intensities are highest in S4 (Figure 8) shows that both high coercivity goethite and hematite are important mineral components in S4. Both are magnetically weak, with mass specific magnetic susceptibilities 3-4 orders of magnitude less than magnetite and maghemite (Hunt et al., 1995). While the absolute concentration of SP particles in S4 is half that of S1 and one third of S3, χ fd varies significantly less being 8.2% for S4, 8.9% for S1 and 9.8% for S3 ( Table 2). The lower χ lf and χ for S4 may simply result from the mineralogy of the neo-formed pedogenic iron mineral, paleoclimate and local biogeochemical conditions favored the formation of hematite and goethite or at least their preservation.
Progressively reduced precipitation and temperature characterize the general climatic trend over the last 430 kyr at Zemun when considering the colorimetric data, gradually leading to more arid climate conditions. A fundamental factor for reduced weathering intensity, indicated by decreasing a * values, may be the general aridification as proposed by Marković et al. (2009) and Buggle et al. (2013) for the Middle and Lower Danube Basins. Buggle et al. (2013) proposed the progressive uplift of the Carpathians, Dinarides, and Eastern Alps during the Pleistocene as main cause of this of the aridification trend.
Since, at Zemun, magnetic susceptibility is dominated by ferrimagnetic minerals, maghemite and magnetite (Figure 5B), and colorimetric analyses by hematite, an integrative view of the results from both methods is needed and beneficial. Regarding the formation of magnetite in LPS, the preferred climatic conditions are seasonally alternating wet (reducing) and dry (oxidizing) conditions controlling pedogenetic processes subsequently leading to higher susceptibilities (Maher, 1998). Differences in relative amounts of maghemite and hematite are probably caused by distinct characteristics and seasonality of climatic periods (Maher, 2011). A climate characterized by warm and dry summers in combination with mild and humid winters leads to relatively higher amounts of hematite (Balsam et al., 2004;Liu et al., 2008;Jiang et al., 2018). The interpretation of our results are in good agreement with interpretations from other studies focusing on the Lower and Middle Danube Basin (Buggle et al., 2013. With regard to the paleosols, a decreasing HGR is associated with a lower seasonality of precipitation and/or a reduction of temperature (Balsam et al., 2004). Based on our results, we assign HGR as an index for increased heat and dryness during estival periods. With respect to the susceptibility parameters, in interglacial periods the mean precipitation over the year cannot be that different. A similar trend from MIS 9 to MIS 5 is reported by Obreht et al. (2016) from the Stalać LPS, who interpret this as changing influence of different climate systems over western Eurasia.

Comparison of the Zemun LPS With Titel-Stari Slankamen and Luochuan
Thick loess-paleosols sequences are spread over the northern hemisphere, located in different geomorphological settings which are influenced by diverse prevailing climate regimes (see Schaetzl et al., 2018;Lehmkuhl et al., 2021). To detect similarities and differences in magnetic and colorimetric properties conserved in loess-paleosol sequences, we compare the Zemun LPS with two well-known reference profiles. The Titel-Stari Slankamen LPS (Serbia) covers the last million years and the interval from MIS 11 to present is archived in c. 38 m.
The Luochuan record from the CLP (Hao et al., 2012) covers approximately the last 1.1 million years and the interval MIS 11 to present in c. 30 m of loess/paleosol couplets. For comparison, the low-frequency magnetic susceptibility data (the only data set available for all three LPSs) are shown on an age scale (Figure 11). Generally, the Zemun LPS's long-term χ lf variability reflecting interglacial/glacial alternations are in excellent agreement with Titel-Stari Slankamen and Luochuan, even though the latter is located in a monsoon prevailing climate. Indeed, it is well known that different prevailing climate regimes and primary sediment sources affect the χ lf values (Maher, 2016;Schaetzl et al., 2018), which are nearly 60% higher in Luochuan than in the both Serbian LPSs. More than doubled χ values in Luochuan compared to Titel-Stari Slankamen (Song et al., 2018) and Zemun indicate that the neo-formation of magnetite/maghemite SP particles during pedogenesis was enhanced at Luochuan, potentially due to more moisture during interglacials. The S4 paleosol at Titel-Stari Slankamen and Zemun shows-compared to the S3 paleosol -less prominent values for χ lf and χ (Figure 11). In all records, the S3 paleosol appears in the χ lf data as a double peak. The timing of onset and demise of the corresponding interglacial MIS 9 (S3) is in better agreement between Zemun and Luochuan than between Zemun and Titel-Stari Slankamen. Variations in χ lf for paleosol S2 (MIS 7) are similar for both Zemun and Titel-Stari Slankamen, but the χ lf values are 2-3 times higher at Zemun. The S2 pedocomplex at Luochuan displays a clear three folded paleosol based on χ lf . Song et al. (2018) correlated the two older paleosols of MIS 7 with the S2 and L2SS1 horizons at Titel-Stari Slankamen. However, comparing Zemun and Titel-Stari Slankamen, a correlation of the S2 pedocomplex with the older two paleosols of Luochuan and a further correlation of the Zemun L2SS1 with the youngest MIS 7 paleosol seems more plausible. Similar to Titel-Stari Slankamen, the S1 of Zemun shows two pedogenetic horizons, a lower strongly expressed paleosol, and a younger, weakly expressed paleosol. This feature is also observable in Titel-Stari Slankamen but differently expressed in Luochuan, where the pedogenic intensity remains fairly constant throughout the S1 pedocomplex. Age differences for the onsets and demises of the paleosols might result from the applied dating techniques. For Zemun, the correlative age model is based on fluctuation in χ and two tephra marker horizons. For Titel-Stari Slankamen, composite record ages were provided by tuning with aid of an astronomical target, correlating peaks in χ lf with June perihelia . Ages for the Luochuan LPS were determined by direct correlations of the fluctuations in χ lf and isotopic data from deep-sea sediments (Ding et al., 2002a;Sun et al., 2006).

Tephra Layers as Widespread Marker Horizons Using Magnetic Proxies to Highlight the Presence of Tephra Layers in Loess
In tephrochronology (Lowe, 2011) ash layers of known age and origin can be used as reliable dating tools and correlation targets. This approach plays an increasingly important role in Quaternary stratigraphy (see Abbott et al., 2020 and references therein) as for example, the accuracy of luminescence numerical dating, the choice of method of dating the emplacement time of loess deposits is still impacted by many limiting factors such as large offsets between the different methods currently applied on the same sample (Avram et al., 2020), or even between different grain-size aliquots (e.g., Timar-Gabor et al., 2011;Veres et al., 2018). The most prominent tephra layer preserved in Eastern European loess and archeological deposits is the Campanian Ignimbrite/Y5 tephra Giaccio et al., 2017). Dated to c. 40 ka BP, it provides an exceptional tie-point for linking records within MIS 3 and in testing the accuracy of luminescence and radiocarbon dating for loess records (Constantin et al., 2012;Fitzsimmons et al., 2013;Anechitei-Deacu et al., 2014;Obreht et al., 2016Obreht et al., , 2017Scheidt et al., 2021). However, several more tephra layers have been identified in south-eastern European loess profiles (see Marković et al., 2015Marković et al., , 2018 which, if well assessed, not only chronologically but also by their bulk sediment geochemical (see Pötter et al., 2021) and magnetic properties (this work) may serve as important anchor points for more secure lateral stratigraphic correlations . At present, glass shard geochemical data are not available for the two tephra layers identified at Zemun (and elsewhere in the wider region) because all attempts of geochemical fingerprinting have failed due to the strong weathering and alteration of the volcanic products.
Indeed, our integrated approach combining high-resolution FORC analysis and IRM unmixing provide magnetic evidence for the presence of volcanogenic material not only for the L2 tephra, but also for the Bag tephra (Figures 6, 7). Both ZV 117 (L2 tephra) and ZP1 001 (Bag tephra) samples are characterized by a higher relative content of SD and PSD grains than the bracketing loess as well as a distinctive coercivity distribution not present in the bracketing loess. As most volcanic ashes identified in European loess units are mixed with loess or occur as cryptotephra, our approach may provide an excellent screening potential for assessing the tephrostratigraphic potential of loess records beyond the classical glass-shard geochemical analyses.
FIGURE 11 | Comparison of the χ lf record from the LPS at Zemun, Titel-Stari Slankamen, and Luochuan on their individual age scales. Note the very similar chronostratigraphic depiction of the potential tephra layers at Zemun and Titel-Stari Slankamen. Differences in onsets and demises of different paleosol units may result from different age models and applied indirect (correlative) dating. Details are discussed in the text.

CONCLUSION AND SUMMARY
In this study, we present a new high-resolution loess-paleosol paleoclimate record from the Middle Danube Basin based on the application of environmental magnetic and colorimetric techniques for the last four-interglacial-to-glacial cycles spanning from MIS 11 to MIS 4. Data stemming from the S2 pedocomplex reveals different stages of environmental change during an interglacial, reflected by different coercivities, magnetic grain sizes, and different magnetic mineralogies. The combination of rock-magnetic and colorimetric data for the Zemun loess profile confirms the proposed aridification trend by Marković et al. (2009), Buggle et al. (2013), and Obreht et al. (2016). Furthermore, our HGR dataset indicates a general decreasing summer heat or dryness over the last 430 kyr in the Middle Danube Basin. Our study shows tracking only the soft ferrimagnetic (magnetite/maghemite) content may not provide a clear interpretation of changes in pedogenic intensity between interglacials, and that colorimetric data, providing insight into hematite and goethite content (high coercivity minerals) in addition are important for reconstructing summer conditions in addition to mean annual soil moisture. For loess samples stemming from the younger L3, L2, and L1 units, both hysteresis data and colorimetric data provide evidence for higher relative concentrations of high-coercivity magnetic minerals (like hematite and goethite) than in paleosols. The hematite/goethite ratio reflects a relatively higher amount of hematite compared to goethite in paleosol layers. Additionally, the investigation and comparison of FORC and highly resolved IRM acquisitions of potentially tephra bearing and pure loess samples provides further evidence of the deposition of volcanic ashes. Even when volcanic glass shards are too weathered for geochemical investigations, the remaining magnetic signal (increased PSD and SD particles in combination with reduced magnetic interactions) can be clearly identified and used as individual anchor points for future comparison and stratigraphy in the MDB and beyond.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CL, UH, and SM designed the study. CL measured all magnetic and parts of the colorimetric parameters, did the data analysis and interpretation with close advice from CZ, FL, YG, DV, and UH, and programmed all analytical R-scripts for this study. MJ conducted the sampling. CL, UH, CZ, FL, DV, and YG wrote the manuscript with close advice from all co-authors. All authors contributed significantly to the article.