Aeolian dynamics at the northern edge of Deliblato (Banat) Sand Sea, Vojvodina, Serbia, at the time of the last deglaciation

Abstract

The Deliblato (Banat) Sand Sea, which is one of the largest areas of аeolian sand in Europe, is located near the Iron Gate, which marks the crossing of the Danube River through the biggest gorge of this river. Here, Danubian alluvium has served as the sand source for the Banat Sand Sea, which was formed primarily through southeasterly (Košava) winds. Utilizing a multi-proxy approach, the objective of this study is to gain a better understanding of the environmental dynamics of the Banat Sand Sea. To achieve this goal, we conducted an analysis of an archive representing an approximately 20-m-thick dune formation on the northern edge of this dune field. Using optically stimulated luminescence (OSL) dating, we calculated aeolian sedimentation rates and dune ages. Sand was deposited here approximately between 17 ka and 13 ka. Magnetic susceptibility, grain size, and colorimetric analyses were interpreted in terms of local paleoenvironmental conditions. Calculated sedimentation rates (SR) indicate intensive aeolian deposition during the study period that range from 483 cm/ka to 502 cm/ka. We compared our data with regional and other European archives, as well as with climatic variations recorded in the Greenland ice core North Greenland Ice Core Project (NGRIP).

Introduction

Improving our understanding of loess and dune-sand dynamics is an essential step in better explaining the interactions among atmospheric systems, sedimentation processes, and paleoclimate (Marx et al., 2018; Albani and Mahowald, 2019). As highlighted in the work of Sipos et al. (2022), the high sensitivity of aeolian activity makes the mobilization of sand a valuable indicator of climatic and environmental changes. The Carpathian (Pannonian) Basin boasts a notable abundance of aeolian landscapes, including loess plateaus and dune fields (Marković et al., 2015; Sipos et al., 2016, 2022; Fenn et al., 2022). Dune fields primarily originated from fluvial deposits of the Danube, the Tisza, and other rivers during the dry and cold climatic phases of the Pleistocene (Borsy, 1990; Smalley et al., 2009). Extensive aeolian activity during the end of the Pleistocene and the beginning of the Holocene predisposed the present-day morphology of the highest sand accumulation in this region, the Banat (Deliblato) Sand Sea (Marković-Marjanović, 1950; Bukurov, 1954; Sipos et al., 2016, 2022; Lehmkuhl et al., 2018).

Erosion and accumulation produced by the Danube River and its tributaries contributed to the overall aeolian process in Carpathian Basin. The steepest slope of lower Danube River course is represented by a narrow gorge known as the Iron Gate (Mihailović, 2021). Before entering this large gorge, the Danube accumulated a vast amount of aeolian sediments during the Pleistocene (Bukurov, 1954; Smalley and Leach, 1978; Buggle et al., 2008; Újvári et al., 2008; Smalley et al., 2009; Obreht et al., 2015; Fenn et al., 2022). During the Holocene and end of the Pleistocene, this material was locally redistributed, mainly under the influence of a southeasterly (Košava) wind system (Gavrilov et al., 2018; Sipos et al., 2022), helping form an extensive sand dune field (Figs. 1 and 2) with the same orientation (SE–NW) as the Košava wind (Fig. 2). Today, the southeasterly Košava wind dominates the climate of the South Banat, Serbia (Tošić et al., 2018). The Košava wind has played a major role in the aeolian geomorphology, both today and as far back as the last glacial maximum, based on current meteorological data, geomorphological evidence, and numerical modeling (Gavrilov et al., 2018). Sipos et al. (2016, 2022) mapped about 1,300 individual dunes in this field, most of which align in the same direction as the prevailing Košava wind (Gavrilov et al., 2018; Fig. 2).

Figure 1. Map of the Quaternary sediments in the Carpathian Basin, with the location of the Banat (Deliblato) Sand Sea highlighted (modified from Lehmkuhl et al., 2018).

Figure 2. Geomorphological sketch of the Banat (Deliblato) Sand Sea and surrounding area on the left (modified from Sipos et al., 2022). The location of the investigated dune is indicated with a red star; cities with meteorological stations are indicated with black dots. Large gray arrow indicates southeasterly (Košava) winds. Wind roses were obtained from Gavrilov et al., 2018. DEM = digital elevation model.

Aeolian sediment is an excellent archive for paleodust dynamics all over the world (Marković et al., 2008, 2011, 2015; Marx et al., 2018; Schaetzl et al., 2018). Research in the Carpathian Basin has been a major asset in this effort, mainly because of the region's loess and loess-like deposits, as well as its aeolian sand/sandy loess, and other loess derivatives (Lehmkuhl et al., 2018). While loess sediments in the southeastern part of the Carpathian Basin have been intensively studied (Varga et al., 2013; Marković et al., 2015), the sandy aeolian sediments have been mostly overlooked. Most of the past research on this sand sea has been largely descriptive (Cholnoky, 1910; Bulla, 1938; Marković-Marijanović, 1950; Zeremski, 1972; Bukurov, 1984; Menković, 2013). The majority of these aeolian sands occur in the South Banat subregion (Fig. 1), where they form the Deliblato Sand Sea (Deliblatska Peščara) (Marković-Marjanović, 1950; Menković, 2013). This vast, sandy terrain evinced longstanding and intensive aeolian processes, with the last dune redeposition occurring 19–11 ka (Sipos et al., 2022).

In this study, we present the first systematic, multi-proxy study of the Deliblato Sand Sea, employing detailed OSL dating, coupled with magnetic susceptibility and colorimetric analyses, on a dune near the Devojački Bunar locality. The main goal of this study is to reconstruct the age of this dune, to couple these data to other multi-proxy data from the region, and then compare the derived depositional processes here with similar aeolian formations in northern Europe.

Study area

The Deliblato Sand Sea, with an area of 800 km², is the largest dune field in the Carpathian Basin and contains some of the largest sand dunes in Europe. This area has been called the “Desert of Europe” or the “European Sahara” (Cholnoky, 1902). The oval-shaped dune field is located immediately north of the Danube River, from which it extends farther southeast to northwest. The northern part of the Banat Sand Sea is covered by loess sediments, explaining why many of the dunes northwest of Kovačica village are not as topographically visible as the dunes to the southeast. The dune field that lacks a loess mantle occupies about 350 km² (Sipos et al., 2022; Fig. 2).

On a historic map of the Timisoara Banat region from 1778, the Deliblato Sand Sea area is represented as a large area of actively moving, aeolian sand. To stabilize this sand movement, the Austro-Hungarian Empire authorities began a reforestation campaign in the early nineteenth century by planting black locust (Robinia pseudoacacia) and other types of trees on the dune field. At that time, the Deliblato Sand Sea terrain covered an estimated 40,660 ha. The stabilization of sand that resulted from this reforestation effort was reported, starting in 1818 (Milenković et al., 2018).

The study site is located in a sand quarry near the town of Banatski Karlovac (Fig. 2). Here, the sand dune stands about 20 m above the surrounding terrain. The site is located in the NW corner of the Deliblato Sand Sea, near the border of dunes and loess-covered dune deposits. This suggests that the formation of the dune could have been affected by sand arriving by saltation from the Košava wind (Gavrilov et al., 2018; Sipos et al., 2022).

Materials and methods

Sampling and description of the research profile

In March 2018, we sampled a fresh pit face of the dune on the northern edge of the Deliblato Sand Sea, near Banatski Karlovac at Devojački Bunar (Fig. 2). Two cores were drilled to facilitate additional sampling. The first core was through the crest of the dune, to a depth of about 8 m, and the second core was at the bottom of the dune, within and through the quarry bottom, covering an additional 6 m. Two sampled vertical subsections were reconciled using spirit level and following the profile stratigraphy. Given that the thickness of the dune in the quarry was about 16 m, the total thickness of the analyzed sequences was approximately 22 m (Fig. 3).

Figure 3. (a) Sketch of the open sand quarry, with position and depth cores and (b) photo of the investigated open pit face (photo by M.B. Gavrilov, 2018).

We collected eight samples from the dune for luminescence dating (Fig. 4). Three samples were recovered from the open pit face by hammering steel tubes into the face of the cleaned exposure. Samples from the cores were taken using a percussion drill with a window sampler. Samples for dose-rate and water-content assessment were collected from the surroundings of the OSL samples.

Figure 4. Lithostratigraphy of the sand dune, also showing stratigraphic locations of the OSL samples. The overlap between the investigated profiles is based on presence of a paleosol at an approximate relative depth of 13 m.

Without overlaps, 270 samples for magnetic susceptibility, grain-size, and color analyses were then collected at intervals of 5–10 cm. Charcoal, biomass, and secondary carbonates were avoided during sampling.

Magnetic susceptibility measurements

Magnetic susceptibility (MS) is a physical quantity that expresses the degree of magnetism of material (Thompson and Oldfield, 1986), or the magnetic properties of the material when exposed to an external magnetic field. Thus MS data represent the amount of magnetic material in the sample (Heller and Evans, 1995; Maher, 2011).

The samples were dried at 40°C for an hour in LAPER laboratory, Faculty of Sciences, University of Novi Sad, Serbia. After this, the material was crushed with mortar and pestle, until it was homogenous.

The low field MS of each sample was measured using a Bartington MS2 instrument at the Laboratory of Physical Geography, Faculty of Sciences, University of Novi Sad, Serbia. The samples were placed in non-magnetic plastic boxes, after which they were carefully compressed with a non-magnetic plastic press. Before closing the box, cotton wool was placed to prevent any movement of the material during the measurement.

Grain-size analyses

Grain-size analyses were performed on a Fritsch Analysette 22 MicroTec laser diffraction grain-size unit at the University of Szeged, Department of Geoinformatics, Physical and Environmental Geography. This device is equipped with both a green laser (λ = 532 nm, p = 7 mW) and an infrared laser (λ = 940 nm, p = 9 mW), enabling it to measure particles within the range of 0.08–2000 μm. The procedure for preparing the samples adhered to the methods outlined in two prior studies (Kun et al., 2013; Serban et al., 2015). Notably, to prevent the modifying effects of a dispersant on the measurements, no chemical dispersion was employed. Instead, a lengthier ultrasonic pretreatment lasting 180 seconds (f = 36 kHz, p = 60 W) was applied (Moine et al., 2008). The generation of grain-size distribution curves involved the application of Mie optical theory to the laser diffraction data. The specific parameters used in this process were a refractive index of 1.52 and an absorption index of 0.1 for the dispersed sample, while a refractive index of 1.33 was used for water (Sipos et al., 2022).

Color analysis

The color of each air-dried sample was measured using a Konica Minolta chromameter 400/410 (Lukić et al., 2014, 2023). The values obtained were given different color coordinates, such as the CIE chromaticity coordinates: Y, x, y; CIE (Commission Internationale de l'Eclairage, 1978) coordinates: L*, a*, b* values. The positive sign in front of +b* value indicates a yellow color, and in front of +a* value represents a red color. The values with negative signs (−b* and −a*) imply that the colors are bluer and greener, respectively. Gray colors impart a value of zero. Indicators a* and b* can have values from −120 to +120. The L* value is a measure of lightness of the sediment, with values ranging from 0, for absolute black, to 100, which represents white (Commission Internationale de l'Eclairage, 1978).

OSL sample preparation

The OSL samples were measured under low-intensity amber light conditions at the Laboratory for Optically Stimulated Luminescence, University of Szeged, Hungary. The sediment for equivalent dose measurements was wet sieved using 63, 90, 150, 180, and 250 μm sieves. The obtained 90–150 μm material (the most abundant material) was first treated for 60 minutes with 10% HCl and for 60 minutes with 10% H₂O₂ in order to remove carbonates and organic matter, respectively. Subsequently the samples underwent a heavy liquid (LST Fastfloat; 2.62 g/cm³) separation to isolate a quartz-rich fraction. The quartz grains (<2.62 g/cm³ material) were next treated with 38% HF for 45 minutes to remove the alpha-irradiated outer coating and any remaining feldspar contamination. The etching was followed by a 10% HCl rinse for 40 minutes to remove fluorides. After each chemical treatment, samples were washed three times in deionized water. Pure quartz extracts were mounted on stainless steel discs using a large 6 mm mask with silicone spray.

OSL measurements and dose rate determinations

OSL measurements were performed using continuous-wave optically stimulated luminescence (CW-OSL) on a Risø TL/OSL reader, model DA-20 (Bøtter-Jensen et al., 2010). The reader was equipped with blue LEDs (470 nm, ~80 mW/cm²), infrared (IR) LEDs (870 nm, ~135 mW/cm²) and a ⁹⁰Sr/⁹⁰Y beta source calibrated for discs using Risø quartz calibration (Hansen et al., 2015). The quartz OSL signals were collected through a 7.5 mm Schott U-340 (UV) glass filter. For the equivalent dose (D_e) determination, a standard, single-aliquot regenerative (SAR) protocol was used (Murray and Wintle, 2000). In order to define the suitable preheat temperature, a combined preheat plateau and dose recovery test was conducted on 15 aliquots (3 aliquots per temperature combination) of sample 1555. First, aliquots were bleached two times using the blue LEDs for 250 seconds at room temperature. Between the two bleaching treatments, a pause of 10 kiloseconds (ks) was employed, after which the aliquots were irradiated with a beta dose of ~20 Gy. Subsequently, preheat temperatures ranging from 180–300°C (20°C increments) were induced while the test dose preheat (cut-heat) was kept at 160°C. The dose recovery ratio showed acceptable values (within 10% of unity) up to a preheat temperature of 260°C, after which the recovered doses underestimated the given doses (Fig. 5). Recuperation and dose recovery recycling ratios were close to unity for all temperatures. A small plateau was observed at 180–200°C.

Figure 5. Results of the preheat plateau–dose recovery test for sample 1555. Three aliquots were measured per temperature point. The cut-heat was kept at 160°C for all preheat temperatures. A dashed line is drawn at an ideal recycling ratio of 1.0. The errors are given as 1 standard error (1σ).

Based on the results of the preheat plateau–dose recovery test (Table S1), a first preheat temperature of 200°C for 10 seconds and cut-heat of 160°C was chosen for all subsequent measurements. To test the suitability of the chosen preheat conditions and the SAR protocol, we performed a dose recovery test on sample 1555 (18 aliquots). Natural aliquots were first bleached twice using blue LEDs at room temperature for 100 seconds (with a 1 ks pause between the two bleaching treatments). Subsequently, a laboratory β-dose was administered to the bleached aliquots, equivalent to the natural D_e (20 Gy) and measured using the chosen SAR protocol (200°C preheat for 10 s; 160°C cut-heat) (Murray and Wintle, 2003). The calculated dose recovery ratio was 0.98 ± 0.02 while the recycling ratio and recuperation displayed values of 1.07 ± 0.02 and 0.15 ± 0.16, respectively. These results suggest that the chosen SAR protocol is able to successfully recover a known laboratory dose prior to any laboratory heat treatment and is suitable for dating samples from the Deliblato sands.

The test dose for the subsequent D_e measurements was set at 10% of the assumed natural value. Recycling and recuperation tests were included in each measurement (Murray and Wintle, 2003), a maximum deviation of less than 10% and less than 5% was considered acceptable for the recycling ratios and recuperation, respectively. To check purity of the quartz aliquots, OSL IR depletion tests were performed at the end of every measurement using IR LEDs stimulation at a temperature of 50°C for 100 seconds (Duller, 2003). For the quartz D_e calculation, the signal from the first 0.16 s of stimulation minus a late background (35–40 s) was used. The luminescence data were analyzed using Risø Analyst software, version 4.57 (Duller, 2015). The annual dose rates and final ages were calculated using the DRAC dose rate and age calculator, version 1.2 (Durcan et al., 2015). The dose response curves for all the samples were best fitted with a single saturation exponential function.

Dose rate samples were dried at 50°C, then crushed and placed in 450 mL Marinelli beakers and stored for more than 3 weeks to allow ²²²Rn to build up to equilibrium with ²²⁶Ra. Next, the concentrations of ²³⁸U, ²³²Th, and ⁴⁰K were measured for more than 24 hours using a high-purity CanberraXtRa Coaxial Ge detector (Murray et al., 1987). The alpha and beta attenuation factors were estimated according to Brennan et al. (1991) and Brennan (2003), respectively. Dose rate conversion factors were obtained using the calculations of Liritzis et al. (2013). Detailed method description can be found in Vandenberghe (2004). The uncertainties associated with etch depth attenuation were estimated according to Brennan (2003). Long-term sample water contents were estimated by using in-situ water contents measured from borehole samples. The cosmic ray contribution to the total dose rates rate was calculated according to Prescott and Hutton (1994). Data are presented in Table 1.

Table 1. Summary of sample codes, depth information, in situ water content (w), radionuclide concentrations, total dose rates (D*), weighted mean D_e values, OSL ages, and total errors in ka for Deliblato sands samples. Also shown are the number of measured aliquots (n) and the number of accepted aliquots (n_a). OD represents the calculated overdispersion. Error terms are given as 1 standard error (1σ) of the mean.

Age-depth modeling, and sedimentation rates

The OSL data allow us to construct an age-depth model for the sand dune at the study site, which can then be used to determine sedimentation rates (SRs) for the interval of its formation. Numerous methods have been utilized for developing continuous age-depth models from luminescence ages (e.g., Haslett and Parnell 2008; Blaauw 2010; Zeeden et al., 2018; Lougheed and Obrochta, 2019), the majority of which use distinctive approaches. In this study, we utilized the Bacon R age-depth modeling software version 2.5.8 (Blaauw and Christen, 2011). In contrast to most age-depth modeling methods for the reconstruction of SRs, Bacon uses Bayesian statistics by combining numerical age dates with prior information (e.g., accumulation rates, its versatility, and its memory change with depth). SRs are controlled using a gamma autoregressive semiparametric model with any defined number of subdivisions along the sand dune profile (Blaauw and Christen, 2011). Prior information on the SRs can result in lower uncertainties and more accurate modeling. The section thickness can, to a certain extent, define the flexibility of the age-depth model, which is attained by Markov Chain Monte Carlo (MCMC) iterations (Blaauw and Christen, 2011). For the investigated site, the depth interval resolution was set to 5 cm and the section thickness to 20 cm (Blaauw and Christen, 2011).

Results

Stratigraphy

Figure 4 shows in detail the lithostratigraphic description of the investigated sand dune profile and cores, indicating the positions of the luminescence samples. The lowest layer of the section (below 1800 cm) is a homogeneous layer of well-sorted fine-grained sand to coarse silt. Large nodules of cemented, secondary carbonate occur at 2060 cm. Between 1800 cm and 1420 cm, the sand varies from fine grained and silty (loess-like) to medium grained without silt, with some charcoal fragments at 1750 cm. This layer is overlain by laminated, silty fine-grained sand with carbonate accumulations. The remains of plant roots are located in the thin layer of loess (between 1420 cm and 1350 cm), perhaps associated with the clearly differentiated overlying paleosol, which occurs between 1350 cm and 1260 cm. The soil has high MS values, in excess of 30×10⁻⁸ m³/kg (Fig. 6). Charcoal fragments were recovered from the dune sand at 1750 cm, 1260 cm, and 1200 cm.

Figure 6. Magnetic susceptibility values (×10⁻⁸ m³/kg) and grain size variations (%), as they correlate to lithostratigraphy and pedostratigrahy. Chronostratigraphic units: I. The upper unit ranges 0–900 cm.; II. The middle unit with paleosol ranges 900–1350 cm.; III. The lower unit ranges 1350–2200 cm. See Figure 4 for lithostratigraphic key.

The next unit (between 1100 cm and 900 cm) is a homogeneous, continuous body of sandy silt (loess-like). Charred plant remains were observed at a depth of 720 cm. Below these remains, at about 850 cm, we also collected charcoal fragments. Between about 900 cm and 180 cm, strata of fine-grained sand layers alternate with layers of medium-grained sand. A zone of secondary carbonates has developed at the top of this layer, between 190 cm and 180 cm. In the dune sand between 430 cm and 420 cm, traces of former grass vegetation are present, while at 480 cm, we observed a thin layer of dune sand rich in muscovite. Fine-textured dune sand occurs below from about 180 cm to 160 cm. Sediments become gradually much siltier towards the top, where an approximately 1-m-thick, poorly sorted, sandy loess was deposited. A 50-cm-thick arenosol is present at the dune crest.

Magnetic susceptibility

Pedogenesis tends to cause higher MS values in pedogenetic horizons than observed in sand and loess layers (Thompson and Oldfield, 1986; Heller and Evans, 1995; Evans and Heller, 2001; Maher, 2011). At this site, MS values range between 12×10⁻⁸ m³/kg and 81×10⁻⁸ m³/kg, with an average value of 28.8×10⁻⁸ m³/kg (Figs. 6 and 7). The lowest MS values occur between 1350 cm and 1750 cm. The highest values were recorded at 890 cm (81×10⁻⁸ m³/kg) and 1930 cm (67×10⁻⁸ m³/kg). The upper MS peak coincides with a layer of charcoal. There are also peaks at depths around 700 cm and 120 cm, but it is not clear what influenced this increase in MS values. The high MS values between 1210 cm and 1345 cm coincide with the weakly developed fossil steppic soil (Figs. 6 and 7).

Figure 7. Magnetic susceptibility values (×10⁻⁸m³kg⁻¹) and sediment colors (L*, a* and b*), as they correlate to lithostratigraphy and pedostratigraphy. Chronostratigraphic units are the same as shown on Figure 6; see Figure 4 for lithostratigraphic key.

Grain size

Variations in the grain size (GS) distribution in comparison to the pedostratigraphy of the investigated dune are presented in Figure 6. The grain-size distribution indicates the variability in the fine-grained (< 2 μm and 2–62 μm) fraction and coarse-grained material (> 62 μm) contents. Generally, the < 2 μm and 2–62 μm fractions correspond well to MS fluctuations, as well as pedostratigraphy. The higher contribution of two finer-grained fractions is related to pedogenic layers and loessic strata. However, GS variations show more cyclic variations than magnetic and color proxies.

We divided the site into three stratigraphic units, each with similarities in MS, GS, and color analyses (marked with roman numbers on Figs. 6 and 7). The lowermost 4 m of the analyzed profile has finer GS composition, characterized by slightly higher contributions of 2–62 μm than the > 62 μm fractions (unit III). Upper part of unit III is represented dominantly by > 62 μm fractions. Lower part of unit II is represented with paleosol, which is clearly visible in the GS data, with significant increases in < 2 μm and 2–62 μm fractions. Interestingly, this section is still coarser than the lowermost 4 m. Unit I is represented by fluctuating GS values and increases in the > 62 μm fractions. The uppermost part of unit I contains the soil, as well as some loess and sand layers, and is represented in GS values.

Color data

Color, which is an important sediment property, may indicate changes in processes during the depositional period (Lukić et al., 2014). Figure 7 shows the changes in the values of L*, a*, and b* across the sediment column. The data can be separated into three stratigraphical units of the dune, all of which correlate well with GS and MS values (Figs. 6 and 7). The modern soil has low values of L* and b*, and high values of a*, indicating, as is typical, that pedogenesis has led to darker, browner colors. Below, the loess has lighter colors, such that by about 200 cm, colors are similar to those of stratified dune sand. This stack of stratified sand extends down to approximately 1250 cm and is interpreted to represent very dynamic environmental conditions of aeolian sand deposition. The paleosol that occurs between 1250 cm and 1350 cm shows up well in the color data; the values of L* and b* are lower, while a* is higher than in the overlying sediment. Interestingly, paleosol developed in about one meter of loess, as indicated by the color data. Finally, the basal sandy sequence is again characterized by significant variations in color, as in the uppermost sandy sequence.

OSL results and ages

Most of the measured aliquots show low OSL IR depletion ratios and recycling ratios (close to unity), suggesting that the correction for sensitivity changes was performed accurately throughout the SAR cycles. The response to a zero dose (recuperation) was less than 1% from the corrected natural signal, where 90% of the measured aliquots passed the rejection criteria of the SAR protocol. The laboratory dose response curves for a representative aliquot of sample 1565 are presented in Figure 8. The natural and regenerated decay curves are almost indistinguishable from each other from the Risø calibration quartz.

Figure 8. Representative sensitivity-corrected dose-response curves constructed for sample 1565 using one aliquot. The sensitivity-corrected natural signal is depicted as an open square and the dashed line indicates the equivalent dose. The response to a zero dose is depicted as an open circle. The recycling ratio point is represented as an open triangle. The inset shows a typical decay curve of natural CW-OSL signal (open circle) in comparison to a regenerated signal (open diamond) induced by a beta dose approximately equal with the equivalent dose and a calibration quartz signal (closed circle).

The calculated OSL ages and dosimetric data, presented in Table 1, provide age control for the sediments at the Banatski Karlovac sand quarry site. The calculated D_e values vary from 21 ± 0.03 Gy for sample 1552 to 36 ± 0.04 Gy for sample 1570, while the measured dose rates ranged from 1.36 ± 0.07 Gy/ka (sample 1559) to 2.25 ± 0.07 Gy/ka (sample 1569), with no apparent depth trend.

The OSL dates may be subdivided into three main groups. Sediment from the first group of samples (1552, 1555, 1557) was deposited at approximately 14–13 ka. These samples represent the last vestiges of sediment accumulation at the study site. The second group includes samples 1559 and 1563, which were deposited at approximately 15–14 ka. The third group of samples, from the bottom of the section (1565, 1569, 1570), is chronologically distinct from the overlying sediment. The basal sediment appears to have been deposited at approximately 17–16 ka, suggesting that the Banatski Karlovac sand dune started forming during the latter part of Marine Isotope Stage 2.

Timescale and sedimentation rates

The modeling results of the Bacon R software are represented in Figure 9. The results are based on more than 62 million Monte Carlo iterations and yield mean 95% confidence age ranges of 1570 years, with 100% of the ages overlapping with the age model. The data displayed only one age inversion of the mean ages (outside the uncertainty window), which is noticeable within 1σ, however, with a much better fit within a 2σ uncertainty within the profile. For this reason, the age-depth models were generated with only a few resampling attempts and yielded stratigraphically consistent data. The Bacon age-depth model illustrates low sensitivity to the OSL ages, which was followed by a linear age-depth function.

Figure 9. Age modeling results and OSL ages from the Banatski Karlovac site using the Bacon model of Blaauw and Christen (2011). The original data and uncertainties are plotted as closed circles with 1-sigma (1σ) error bars. The gray shaded area shows the minimum and maximum age limits with 95% certainty, while the dashed line represents the mean age model.

The calculated SRs based on the OSL ages are shown in Figure 10 and Table S2. The data were smoothed, for a more realistic representation, using Sigmaplot software version 11.0 (sampling proportion = 0.100; polynomial degree = 1). The SR values range from 484.2 ± 1 cm/ka to 501.7 ± 2 cm/ka with a median and mean of 497.6 cm/ka and 495.3 cm/ka, respectively. Three peaks in sedimentation rate occur at 14 ka, 15 ka, and 16 ka with average values of approximately 500 cm/ka and approximately 502 cm/ka, respectively. Sudden decreases in sedimentation rates occurred at ca. 14.5 ka and 15.5 ka, where lower values of approximately 484 cm/ka and approximately 489 cm/ka were determined (Fig. 10).

Figure 10. Correlations between (a) oxygen isotope values and (b) dust variations in the NGRIP2 ice core (Rasmussen et al., 2014), with (c) magnetic susceptibility data and (d) sedimentation rates for the investigated dune. Dust shown as particles (n) per mL (n/mL). GIS1 = Greenland Interstadial 1; GS2 = Greenland Stadial 2

Discussion

Environmental dynamics

Based on the profile description and presented luminescence chronology, we identified three chronostratigraphic units that represent the main environmental phases that occurred during dune formation (Figs. 4, 6, and 7). Unit I, the upper unit, ranges from 0–900 cm and was deposited between 14–13 ka. Unit II, the middle unit with paleosol, spans from 900–1350 cm and was deposited between 15–14 ka. Unit III, the lower unit, extends from 1350–2200 cm and was deposited around 17–16 ka.

However, based on our data, we were unable to reconstruct the complex geomorphological relationships among aeolian sand and loess deposition, intervals of soil formation, and prevailing wind direction(s) at the northern limit of the huge Deliblato Sand Sea (Sipos et al., 2022). The chronologies for the sediments here, obtained via OSL dates, provide the timeframe for the environmental dynamics at this site, using sedimentation rate data (cm/ka).

Relationships to aeolian deposits in the Carpathian Basin

The Carpathian Basin has many kinds of aeolian landscapes, including loess plateaus and dune fields (Lehmkuhl et al., 2018; Obreht et al., 2019; Fig. 1). Its dune fields were primarily shaped by deposition of sediment derived from the Danube, Tisa, Tamiš, and other river valleys during the dry and cold periods of the Pleistocene (Borsy, 1990). In cold/dry times when vegetation was sparse, winds could erode the former alluvial surfaces, alluvial fans, and terraces that had been abandoned due to changes in major river courses (Bukurov et al., 1982; Gábris, 2003; Mezősi, 2016). Various geomorphological and dating methods have been applied to study the dune fields, mostly in Hungary, resulting in an extensive collection of research (Borsy, 1991; Kiss et al., 2009, 2012a, b; Buró et al., 2016). Based on these studies, the primary phase of aeolian activity was dated to the last glacial maximum (LGM), during which dune formations could arise wherever sand was accessible due to the lack of vegetation. Nevertheless, evidence suggests that renewed dune formation took place after the LGM, from 19–11 ka, as well as into the Holocene (Nyári et al., 2007, Shakun and Carlson, 2010; Kiss et al., 2012a, b; Sipos et al., 2016, 2022). This phase of climate-driven aeolian activity persisted for the longest duration in the dry, central and eastern parts of the Carpathian Basin (Mezősi, 2016).

In contrast to the well-studied dune fields in Hungary, similar aeolian deposits in Serbia have received limited investigation. Researchers studying the Deliblato Sand Sea have primarily relied on geomorphic observations, leading to various hypotheses and assumptions about the origin and timing of its formation (Cholnoky, 1910; Marković-Marijanović, 1950; Zeremski, 1972; Menković, 2013). However, recent papers by Gavrliov et al. (2018) and Sipios et al. (2022) offer fresh perspectives that enhance our understanding of dune formation in the Deliblato Sand Sea.

Based on dune morphologies and elevations, the majority of the dune field is dominated by elongated, hairpin-like, parabolic dunes, indicating a relatively low supply of sand during their formation (Kiss et al., 2012b; Mezősi, 2016). However, another set of dunes, primarily of the filled parabolic type, overlies the previous ones, pointing to a later period of much higher sand supply. According to our OSL ages, sand began accumulating at our study site during the cold and arid LGM, approximately 19,000 years ago. Subsequently, a significant phase of additional dune formation continued from the Younger Dryas and into the Boreal period of the Holocene. However, due to increased precipitation after 8000 years ago, vegetation helped stabilize the dune landscape. Consequently, on most of the Deliblato Sands, the dunes are older and more stable than previously assumed (Sipos et al., 2022).

Poleva Cave, in southwestern Romania, is situated in the southern part of the Locvei Mountains, approximately 10 km north of the Danube Gorge and about 50 km from the investigated dune. Stalagmite PP10 from Poleva Cave provides a late Pleistocene and Holocene isotopic record characteristic of this region. The δ¹⁸O signal is considered to reflect primarily regional changes in average temperatures, following a positive relationship (i.e., heavier δ¹⁸O values are interpreted to indicate warmer conditions; Constantin et al., 2007). This record reveals a growth hiatus between approximately 22 ka and 11 ka, which could correspond to a dry climatic period favorable for the intensification of aeolian activity in the investigated area. Additionally, there have been a number of glacial advances in the southern Carpathians in the post-LGM phase (Urdea et al., 2023). However, the geomorphological landscape of the highest southwestern Carpathian massifs indicates that during the deglaciation period, the glaciers remained positioned within the major relief landforms created by the last glacial maximum (LGM) glaciers, indicating smaller amplitude of climatic and environmental dynamics.

Relation to North European aeolian environments

Aeolian sands associated with periglacial conditions extend over a vast belt in the North European continent, from northern France to northern Russia (Koster, 1988; Kasse, 1997, 2002). Patchy deposits also have been reported in SW France (Sitzia et al., 2015; Bertran et al., 2016).

Typically, coversands, deposited as sheets, are characterized by a specific, very finely horizontal–parallel laminated structure (Schwan, 1988, Van Huissteden et al., 2000; Vandenberghe and Wolfe, 2023). As demonstrated by mineralogical analyses, these sheet sands become homogenized by saltation over distances on the order of tens of km (Vandenberghe and Krook, 1981, 1985), resulting in a relatively consistent grain-size composition. In contrast, periglacial dune sands are characterized by their well-expressed dune morphology, some up to a few tens of meters in relief. Their sedimentary structures consist of (low- to high-angle) crossbedding, (sub)horizontal lamination, and occasionally homogeneous beds. Transport was by saltation or in low suspension clouds (De Ploey, 1977) over short distances (tens or hundreds of meters). In contrast to the coversands, these sands have variable granulometric and mineralogical compositions due to their local provenance.

After 17 ka, postdating the last permafrost maximum (LPM), coversands were deposited across a North European coversand belt (e.g., Vandenberghe, 1985, 1991; Kozarski, 1987; Goździk, 1991; Kasse, 2002; Kasse et al., 2007, Zieliński et al., 2011). Coversand deposition was fostered by aridity with reduced vegetation cover, attributed to (apart from drier climatic conditions) the disappearance of permafrost, which allowed for increased infiltration (Kasse, 1997). During the late glacial (14.7–11.9 ka), dune formation was more prominent, although coversand deposition locally continued until 12.7 ka (Kasse, 1997; Kasse et al., 2007). Especially during the Older and Younger Dryas, dunes expanded considerably in north-central Europe (Nowaczyk, 1986; Manikowska 1994, Bohncke et al., 1995; Zeeberg, 1998), and during the climatically dry end of the Younger Dryas (10.5–10.1 ka). Dunes formed especially at this time along valley margins, where sand supplied from (braided) floodplains was captured by the vegetation on the adjacent higher dry areas (Vandenberghe, 1983, 1991; Bohncke et al., 1993).

Recently, van Hateren et al. (2022) introduced a combination of endmember analysis of particle size and shape to identify detailed depositional facies in the considered end-pleniglacial eolian sands. Sands of specific grain size transported by saltation in a humid environment occur especially in the sediments shortly after about 17 ka at our study site. Sand grain size shifted to an aeolian facies, deposited in drier conditions, with horizontal bedding during the final phase of the pleniglacial. Similar dry and wet coversand facies were recognized before by Ruegg (1983), Schwan (1988), Kasse (2002), and Vandenberghe (1991). They clearly differ from the somewhat coarser-grained and well-sorted dry eolian dune sands that were deposited especially during the colder phases of the late glacial. A silty fraction, supplied in suspension, was trapped within the vegetation cover during the Bölling interstadial (van der Hammen and Wijmstra 1971).

For further direct comparison between these two European aeolian sand zones we still need to improve our knowledge of chronological frame and more detailed understanding of sand dune-field formation in the scope of environmental dynamics reconstruction. But we can say that both regions had intensive aeolian activity between the LGM and Early Holocene (Vandenberghe, 1985, 1991; Bohncke et al., 1993; Kasse et al., 2007, Zieliński et al., 2011; Kiss et al., 2012b; Buró et al., 2016; Sipos et al., 2022). Due to different air circulations in these regions, we can conclude that there were favorable conditions for strong wind and less vegetation in this period in central and southeastern Europe (Gavrilov et al., 2018).

Relation to NGRIP2 ice core record

The main sedimentation interval of the dune spans the period between approximately 17 ka and 13 ka (Fig. 10). Simultaneously, on a global scale, several significant climatic fluctuations occurred, such as Heinrich event 1 (18.5–16.8 ka), the Older Dryas (16.8 to ca. 14.7 ka), Bølling (14.7–14.1 ka), and Allerød (14.0–12.9 ka). The Older Dryas and HE1 are equivalent to the Greenland stadial 2 (GS-2) and interstadials Bølling and Allerød correspond with Greenland interstadial 1 (GIS-1) (Björck et al., 1998). Generally, sedimentation rates at the studied site are slightly higher during stadial phase GS-2 than during the interstadial GIS-1. Calculated SR values for the dune are generally high, however, varying from 484 cm/ka and 502 cm/ka and appear to suggest continuously favorable conditions for rapid aeolian sedimentation, at least locally. Generally, the SR data do not correspond well to the main shift in δ¹⁸O and dust concentration records from the NGRIP2 ice core (Ruth et al., 2003; Gkinis et al., 2014; Rasmussen et al., 2014) (Figure 10). Therefore, variations of SR appear to be controlled more so by changes in local environmental dynamics than by global-scale forcings. For example, the drop in SR values at about 15.4 ka coincides with the formation of weakly developed steppic fossil soil at between 1260 cm and 1350 cm. The appearance of this paleosol is also manifested in a slight increase in MS, finer GS composition, and increase of a* color values, simultaneously with a decrease in L and b* color proxies. This environmental event is not evinced in NGRIP2 ice core records (Figs. 6, 7 and 10). The soil was likely caused by changes in local environmental dynamics, as characterized by the dominance of stable grassland vegetation associated with soil formation, leading to reduced inputs of aeolian material.

Conclusions

Our multi-proxy investigation of a 22-m-thick sand dune in the Banat Sand Sea, near Banatski Karlovac in northern Serbia, and at the southeastern margin of the Carpathian Basin, has established this site's high-resolution record of late pleniglacial climatic and environmental dynamics. Our aeolian record provides a unique opportunity for reconstructions of local and regional environmental processes and conditions during the period from approximately 13–17 ka. The dune, which is part of a larger sand sea that grades into sandy loess cover to the north, thus represents an ecotone of the aeolian systems for this region. The dune is part of a larger sand sea that grades into sandy loess cover to the north, thus representing an ecotone of the aeolian systems for this region.

Data on morphology, color, and magnetic susceptibility, plotted against an age model based on eight luminescence dates, provide detailed temporal information on local environmental and aeolian dynamics represented with three environmental phases of dune formation: 13–14 ka, 14–15 ka, and approximately 16–17 ka. Based on obtained OSL ages, sand accumulation in the north part of Banat Sand Sea commenced after the cold and arid last glacial maximum (LGM). These findings are largely consistent with the results reported by Sipos et al. (2022) for other regions within the Banat sand zone. Although SR values are generally high (484 cm/ka to 502 cm/ka), they suggest consistently favorable conditions for rapid aeolian sedimentation in this part of the Carpathian Basin. In summary, this study provides an important link between the aeolian dynamics of the Banat Sand Sea and the previously intensively investigated surrounding loess deposits in this region.