Malaco temperature reconstructions and numerical simulation of environmental conditions in the southeastern Carpathian Basin during the Last Glacial Maximum

We investigate the glacial climate conditions in the southeastern Carpathian Basin (Vojvodina, Serbia) based on the reconstruction of malacological palaeotemperatures and results from a high‐resolution regional climate simulation. Land snail assemblages from eight loess profiles are used to reconstruct July temperatures during the Last Glacial Maximum (LGM). The malacological reconstructed temperatures are in good agreement with the simulated LGM July temperatures by the Weather Research and Forecast model. Both methods indicate increasing temperatures from the northwestern towards the southeastern parts of the study area. LGM aridity indices calculated based on the regional climate model data suggest more arid conditions in the southeastern parts compared with more humid conditions in the northwestern parts. However, for present‐day conditions, the moisture gradient is reversed, exhibiting more humid (arid) conditions in the southeast (northwest). An explanation for the reversed LGM aridity pattern is provided by an analysis of the prevailing wind directions over the South Banat district (Serbia). The prevailing moist northwesterly winds during summer are not able to compensate for the annual lack of moisture induced by the dry winds from the southeast that are more frequent during the LGM for the other seasons.


Introduction
The rates of temperature and precipitation changes are the most important indicators in determining the effects of past, recent and future climate change. To get an even more complete picture of the climate conditions, a joint analysis of these variables is of great benefit. Parameters defined by a mathematical quotient/ratio between precipitation and temperature values are known as the aridity indices and are used to determine environmental conditions regarding water availability or deficit. According to the American Meteorological Society (AMS, 2000), aridity is the degree to which a climate lacks effective, life-promoting moisture. An overview of the characteristics and classification of aridity indices can be found in Gavrilov et al. (2019a,b), for example.
Recent research on aridity in the southeastern Carpathian Basin, with a special focus on the Province of Vojvodina (Serbia), has been conducted on several occasions. A study by Hrnjak et al. (2014) describes the aridity obtained by the De Martonne (DM) Aridity Index (De Martonne, 1925) and the Pinna Combinative Index (Zambakas, 1992) for the period 1949. Gavrilov et al. (2019a used the Forestry Aridity Index (FAI) for the analysis of the interaction of climate and vegetative processes, especially in forestry for the period from 1949 to 2015. In the third case (Gavrilov et al., 2019b), aridity was investigated again using the DM Aridity Index and the FAI, where the aridity results were compared with the results of tree ring records (Tree Ring Width Index) and the Palmer Drought Severity Index (PDSI, Palmer, 1965) for the period 1927-2017. The results show a good agreement between climate conditions obtained by the two mentioned aridity indices, the PDSI and the reconstructed tree ring climate as independent parameters. Finally, the aridity of the common territory of Hungary and Vojvodina was investigated using the DM Aridity Index for the period from 1931 to 2017, and drought conditions were analysed using the PDSI for the period 1920-2014 . Since the DM Aridity Index is frequently used in studies about aridity in the considered region, it is applied to data from a regional climate simulation to determine the aridity for the Last Glacial Maximum (LGM; Clark et al., 2009) in this study. This enables a direct assessment of the changes in recent and palaeoclimate environments. Additionally, the UNEP (United Nations Environment Program) Aridity Index (UNEP, 1992) is applied to the LGM simulation. This index is suitable to describe aridity based on the data from (regional) climate models (Ludwig et al., 2018) since it requires potential evapotranspiration, a quantity provided by the models directly, but apparently not available by default in observational data sets. The UNEP Aridity Index is closer to the AMS definition of aridity because it is composed of the ratio of two wet parameters (precipitation and potential evapotranspiration, see Equation (1)), and thus will provide a better insight into LGM aridity. To quantify the simulated palaeoclimate environmental conditions, a comparison with proxy data is essential (Braconnot et al., 2012;PAGES Hydro2k Consortium, 2017;Ludwig et al., 2018). Such a comparison allows (i) the capability and uncertainties of atmospheric models to represent glacial conditions to be assessed, (ii) the uncertainty of climate reconstruction based on proxy data to be assessed, and (iii) changes in (large-scale) atmospheric dynamics to be determined to understand (regional) environmental changes (Ludwig et al., 2016).
For the reconstruction of palaeoclimate environments, studies on loess profiles play a crucial role. Loess covers a large part of the mid-latitudes in the northern hemisphere and thus provides a great archive for spatial and temporal reconstructions of past climates and environmental dynamics . However, a significant limitation is a lack of quantitative proxies for the validation with palaeoclimatic reconstructions (e.g. Obreht et al., 2019). In this study, we compare the results of reconstructed July malacological temperatures during the LGM period from land snail assemblages at eight loess sections (Madaras in southern Hungary, Zmajevac in northeastern Croatia, and Crvenka, Zemun, Mošorin, Titel, Irig and Požarevac in northern Serbia) with simulated temperatures and aridity from the high-resolution output of a regional climate model. Additionally, sedimentological and geomorphological evidence is compared with the numerical modelling results for the prevailing winds in the study area.

Study area
The southeastern part of the Carpathian Basin encompasses the confluence area of the Danube, Sava, Drava, Tisa and Morava rivers (Fig. 1). Over 60% of this lowland area is covered by loess and loess-like sediments (Marković et al., 2008;Lehmkuhl et al., 2018). The loess-palaeosol sequences of the investigated region represent the most detailed archive of climatic and environmental fluctuations during the Middle and Late Pleistocene on the European continent (Marković et al., 2009. The recent climate of the investigated region is moderate continental, with cold winters, hot and humid summers, a wide range of extreme temperatures and irregular distribution of monthly rainfall amounts. The diversity of the climate in the study area is also influenced by surface winds, which blow from two prevailing opposite directions; from the northwest when it is cold and humid, and from southeast when it is warm and dry Tošić et al., 2018). The mean annual surface air temperature is about 11°C (Gavrilov et al., 2015(Gavrilov et al., , 2016 and the annual amount of precipitation is approximately 600 mm . The diversity of the climate of Vojvodina and the adjacent areas also makes it a challenging task for climate modelling. Global climate models used in CMIP5/CMIP6 are operated on relatively coarse grids (300-120 km horizontal grid spacing (HGS)); Yukimoto et al., 2019;Sidorenko et al., 2019;Hajima et al., 2020). This coarse resolution may lead to an inadequate representation of the regional to local climate. The utilisation of regional climate models can help to overcome these issues and to provide climate data at a much higher resolution by considering regional features like topography (Ludwig et al. 2019;Giorgi, 2019).
Reconstructions of temperature and humidity provide evidence for increasing mean July palaeotemperatures and decreasing humidity from the northwestern to the southeastern parts of the Carpathian Basin during the LGM period (Krolopp and Sümegi, 1995;Sümegi and Krolopp, 2002). Northern Serbia is presently the warmest and driest part of the Carpathian Basin, and similar climatic gradients might have existed during the LGM. In contrast to the "classical" European periglacial loess zone, the loess of the southeastern Carpathian Basin was deposited beyond the periglacial area. Previous investigations by Ložek (1964Ložek ( , 2001 and Kukla (1975Kukla ( , 1977 suggest that most of the European loess belt formed under a cold and relatively humid continental climate. The absence of cryogenic features and land snail assemblages with only a few cold-resistant species being found, provides strong evidence that the LGM loess deposits in the investigated region accumulated mainly under dry and relatively warm conditions (Marković et al., , 2007(Marković et al., , 2008Hupuczi et al., 2010;Molnár et al., 2010;Sümegi et al., 2016;Gavrilović et al., 2020).

Climate model data
Climate models are based on well-established physical principles and are widely accepted tools to gain a better Figure 1. Geographical overview of the study area and locations of considered loess profiles (see Table 3) for the reconstruction of palaeotemperatures from malacological data.
understanding of the complex interactions between the different components of the climate system (Flato et al., 2013). They are used for a variety of research questions including studies of climate system dynamics, projections of future climate and studies of the evolution of past climates, environment and life, known as palaeoclimate modelling (Haywood et al., 2019). Over time, more complex physical assumptions, e.g. about atmospheric behaviour, have been implemented and appropriate mathematical and sophisticated numerical methods are being continuously further developed. Increasing computing power enables more complexity and a higher spatial and temporal resolution that ultimately provides more realistic climate data sets (e.g. Flato et al., 2013;Gavrilov et al., 2014). A related research direction deals with the application of regional climate modelling. While global climate models usually provide global data sets with a coarse resolution (300-120 km HGS), regional climate models are able to considerably refine the data (50-10 km HGS) over a limited area of interest. To achieve this aim, regional climate models are driven by initial and boundary conditions of the global model simulations in a so-called downscaling approach. The main advantages of regional climate modelling are, e.g. a more realistic representation of the topography and an improved physical representation of atmospheric processes. Recently, regional climate models have been increasingly used to tackle palaeoclimate research questions (Ludwig et al., 2018(Ludwig et al., , 2019(Ludwig et al., , 2020. In particular, they are able to bridge the gap between the coarse resolution of global climate models and the regional to local scales, and may be key to advancing a meaningful joint interpretation of climate model and proxy data. To obtain high-resolution LGM climate data for the study region, we use results from a numerical downscaling experiment with the regional WRF (Weather Research and Forecast) model, version 4.1.2 (Skamarock et al., 2019, and references therein). The WRF model has been forced by initial and lateral atmospheric boundary data from the global PMIP3 (Paleoclimate Model Intercomparison Project phase 3, Braconnot et al., 2012), MPI-ESM-P simulations for the LGM (Jungclaus et al., 2012a) and the pre-industrial period (PI; Jungclaus et al., 2012b). The PI experiment is used as present-day reference data for comparison with the simulated LGM climate conditions. The global LGM MPI-ESM-P simulation provides climate data at a coarse resolution (~200 km HGS) for a 100-year time slice being representative of the climate 21 000 years ago. The simulation takes into account the different boundary conditions like orbital forcing (insolation), trace gas concentrations (CO 2 , CH 4 ), the northern hemisphere ice sheets and a related drop in sea level (−120 m) applicable at that time. These boundary conditions are obtained either by theoretical models (insolation) or by different sources of proxy data and derived reconstructions. A summary of those guidelines has been made available to the modelling community by PMIP3 (https://pmip3.lsce.ipsl.fr/). To account for the LGM boundary conditions in the regional WRF model, several adjustments based on PMIP3 standards have been made to consider changes in trace gas concentrations, orbital forcing, ice sheets and land-sea-mask distribution with respect to the lower sea level. Fig. 2 provides an overview of the adapted LGM surface boundary conditions over Europe in the WRF model. The climate data used in this study includes a 30-year time slice simulation with high resolution (~8.5 km HGS) over a limited area by the WRF model. A more detailed description of the WRF LGM simulation design as well as an evaluation of the model performance in comparison with present-day and proxy data can be found in Ludwig et al. (2020).
We have extracted the relevant climate variables for this study from the 6-hourly output of the WRF simulation. Mean values (monthly, seasonal, annual) were computed for the entire 30-year period covered by the simulation. The relevant variables are mean surface (2 m height) air temperature (T2m), mean maximum surface air temperature (T max 2m), surface skin temperature (Ts), all in°C, the zonal and meridional wind components (u, v) at 10 m height in m/s, the total precipitation amount (Tp) in mm, and potential evapotranspiration (PET) in mm.

UNEP Aridity Index
According to the original definition, the UNEP Aridity Index (AI UNEP ) is used to analyse the availability or deficit of water on a monthly, seasonal and annual basis, as a ratio of the total amount of precipitation (Tp) and potential evapotranspiration (PET) (Gao and Giorgi, 2008;Ludwig et al., 2018;Derdous et al., 2021). The AI UNEP is a non-dimensional quantity, defined as: where all values have the previously defined meaning. The climate classification based on the AI UNEP is given in Table 1.

The De Martonne Aridity Index
The DM Aridity Index has been used to analyse aridity on a monthly, seasonal, and annual basis in many countries, including Greece (Baltas, 2007), Turkey (Deniz et al., 2011), Romania (Croitoru et al., 2013), Iran (Tabari et al., 2014), Spain (Moral et al., 2017), Algeria (Derdous et al., 2020) and Hungary , as well as Serbia Bačević et al., 2017;Radaković et al., 2018;Gavrilov et al., 2019bGavrilov et al., , 2020. The DM Aridity Index (AI DM ) is defined as: where all values have the previously defined meaning. C = 10°C is the DM constant, and AI DM has the dimension mm/°C. In this study, the numerical value of AI DM will be presented in the usual manner without the mentioning of dimension. The climate classification based on the DM Aridity Index is given in Table 2.

Simulated summer temperatures
The comparison between the simulated MPI-ESM and WRF July 2 m temperatures (T2m) for the LGM as well as the comparison with the corresponding simulated PI climate conditions demonstrates the strong benefit of using high-resolution regional climate model data over the study area ( Fig. 3(a)-(d)). The WRF simulation reflects much better the regional conditions, e.g. terrain information, and thus provides more realistic temperature maps over the study area. While the temperature distribution in MPI-ESM over the Carpathian Basin is more or less uniform, the topographic exposures of the mountain ranges (the Dinaric Alps and the Carpathian Mountains) are evident in the WRF simulation ( Fig. 3 in Ludwig et al., 2020). The simulated T2m in Vojvodina for the LGM are higher in the WRF model (mean: 16.7°C) compared with MPI-ESM (14.2°C). Note that T2m for Vojvodina refers to only one grid point in MPI-ESM and 396 grid points in WRF. This enables a more detailed regional view of the temperature distribution based on the WRF simulation. Likewise, the WRF-simulated temperature differences between LGM and PI are not uniformly distributed over the domain in contrast to the differences in MPI-ESM. The WRF model produced more pronounced temperature differences in the north of the study area (−7.5°C colder) than in the southeastern parts (−5°C colder; Fig. 3(d)). For the WRF simulations, the gradient of the temperature differences between the northern and southern parts of the domain is reflected in all temperature variables.

UNEP Aridity Index
The annual UNEP Aridity Index (AI UNEP ) calculated from WRF data exhibits semi-arid conditions (AI UNEP : 0.2-0.5) over  Vojvodina, except for a small area where dry subhumid conditions dominate (AI UNEP : 0.50-0.65) (Fig. 4a). This subarea coincides with the position of Fruška Gora Mountains, which is an important refugium for the restoration of vegetation after the LGM (Marković et al., 2008). This regional diversity supports again the need for high-resolution climate data, to realistically reflect the regional environmental variations.
For the JJA (June-July-August) period, predominately semiarid conditions (AI UNEP : 0.2-0.5) are simulated including the zone of the Fruška Gora Mountains (Fig. 4(b)). In southeastern Vojvodina (South Banat district), a zone of increased aridity is simulated (AI UNEP : 0.05-0.20). It is interesting to note that the humidity in this subregion slightly increased in the recent period in relation to the rest of Vojvodina Gavrilov et al., 2019a,b). The simulated July aridity is similar to that for JJA, but the area under arid conditions (AI UNEP : 0.05-0.20) is extended to the eastern half of the investigated territory (Fig. 4(c)).

De Martonne Aridity Index
The WRF-simulated annual DM Aridity Index (AI DM ) points to very humid conditions (AI DM values of 35-55) in the entire territory of Vojvodina and the surrounding areas ( Fig. 5(a)). Over the Fruška Gora Mountains humidity is again increased (AI DM > 55). The recent DM annual aridity of Vojvodina is in the range of 26.4-27.5 to 29.8-30.9 . Since AI DM points to very humid conditions on an annual basis in comparison with more arid conditions provided by AI UNEP , it is questionable whether this index is useful under the glacial climate conditions during the LGM.
parts of Vojvodina, while drier conditions occur in the southeastern parts (semi-arid, AI DM values of 10-20, Fig. 5(b)). The Fruška Gora Mountains remain slightly more humid (semi-humid; AI DM values of 24-28). Likewise, the western parts of the study area remain more humid, which is in contrast to overall drier conditions as given by AI UNEP . It should be noted that recent summer climate conditions in Vojvodina and adjacent regions are on average slightly more humid with AI DM values reaching 20-24 to 24-28 . Thus, the climate conditions in the LGM during JJA were drier, which is known (Marković et al., , 2007(Marković et al., , 2008, and can be confirmed by the WRF simulation. While AI DM and AI UNEP differ substantially for annual aridity conditions, both indices provide similar aridity patterns for the summer season over the study area, except for the more humid conditions given by AI DM along the far western part.

Malacological palaeotemperatures and simulated LGM temperatures
Mean LGM July malacological palaeotemperatures are calculated for eight locations within or close to the study area (Madaras, Zmajevac, Crvenka, Zemun, Irig, Mošorin, Titel and Požarevac) based on the malacothermometer method of Sümegi (1989Sümegi ( , 1996 and Herelendi et al. (1992). In general, July malacological palaeotemperatures were increasing from the northwestern to the southeastern parts of the Carpathian Basin during the LGM period Sümegi and Krolopp (2002). A similar spatial distribution of the LGM July malaco palaeotemperatures is also reconstructed in this study for the investigated region (Table 3). The highest July malaco palaeotemperatures (19.6 ± 1.2°C) were calculated for the loess site Požarevac in the southeast and the lowest determined July malaco palaeotemperature (15.1 ± 2.3°C) is obtained for the Madaras loess section in the northwest. July palaeotemperatures obtained for Titel (18.8 ± 0.4°C), Irig (18.7 ± 0.4°C), Zemun (18.8 ± 0.7°C) and Mošorin loess sections (17.5 ± 0.3°C) are slightly cooler than for the Požarevac loess sequence. Similar July temperatures were calculated for Zmajevac and Crvenka loess profiles, with temperatures between 15.7 ± 2.3°C and 16.6 ± 2.1°C, respectively. The reconstructed July palaeotemperatures are colder than the current July air temperature in that region, which is between 21.4°C and 21.7°C, according to Radaković et al. (2018), pointing towards a stronger cooling in the northern parts of the study area during the LGM. The reconstructed malacological temperatures indicate July temperatures for LGM that were 2.0°C to 7.0°C colder than present-day conditions. This difference is in the same order of magnitude as simulated by the WRF model for LGM and PI conditions (Fig. 3(d)). Likewise, the gradient of the WRF-simulated temperature differences between LGM and PI is increasing towards the north, which corresponds to the stronger cooling of the reconstructed temperatures at the northern loess sections. However, such a gradient is not reproduced in the MPI-ESM simulations as the study area is only represented by a single grid point ( Fig. 3(b)). This in turn demonstrates again the added value of regional climate modelling for such a model-data comparison. For comparison with the reconstructed malacological temperatures, the WRF-simulated temperatures at the loess profile sites are indicated in Table 3. The simulated mean 2 m temperatures (T2m) underestimate the reconstructed malacological temperatures by approximately 1-2°C, except for the loess site at Madaras, where simulated temperatures are warmer, and the Zmajevac site where the temperatures are similar. However, the geographical temperature gradient between the coldest (northwest) and warmest (southeast) reconstructed loess sites is reproduced correctly by the model simulation. Comparison of the malacological palaeotemperatures with the mean maximum 2 m July temperatures (T max 2m) indicates that malaco temperatures reflect, rather, mean July temperatures, since T max 2m largely overestimates the reconstructed values. Likewise, the simulated surface skin temperatures (Ts) are higher than the reconstructed malaco palaeotemperatures, although they correspond well with the malaco temperature at Irig.

Discussion
The comparison of simulated and reconstructed July temperatures for LGM reveals a general agreement between the two methods. Although the simulated mean July temperatures (at 2 m height) by the WRF model slightly underestimate the reconstructed malacological temperatures at loess profiles in the southeastern locations, the reconstructed temperature gradient (colder in the northwest to warmer in the southeast) is rendered by the model. The mean difference between reconstructed and simulated mean July temperatures (T2m) at the eight considered loess profiles is −1.0°C. A comparison of the reconstructed malaco palaeotemperatures with maximum mean July temperatures and surface skin temperatures exhibits larger differences: T max 2m overestimates the reconstructed temperatures on average by 5.8°C, the overestimation obtained by Ts is 2.1°C. Based on our findings, we conclude that the reconstructed malaco temperatures are most representative for T2m.
An interesting finding in this context is the relationship between reconstructed malaco temperatures and the simulated aridity indices for JJA/July. While the colder malaco temperatures in the northwestern part of the study area are correlated with less arid conditions, the locations of the loess profiles with warmer reconstructed temperatures exhibit more arid conditions. For the reconstructed temperatures at the studied profiles, we thus can assume a correlation with both temperature and aridity, which will be tested in a separate subsequent study. A short note should be given on the two utilised aridity indices. Both indices (AI UNEP and AI DM ) provide comparable results for environmental conditions during JJA and July. However, both indices differ substantially from each other based on annual mean climate data. This can be explained by the different variables included in the calculations. In particular, the consideration of T2m in the calculation of AI DM could be potentially problematic. In comparison with present-day data, the annual mean temperature decrease during the LGM is much stronger than the temperature decrease considering only summer. Thus, the AI UNEP, which is based on precipitation and potential evaporation obtained from the numerical simulations, might be more suitable for characterising aridity/humidity conditions for glacial periods.
Finally, the pattern of the LGM aridity during summer is discussed. The increasing aridity from the northwestern towards the southeastern parts of the Carpathian Basin during the LGM period (Krolopp and Sümegi, 1995;Sümegi and Krolopp, 2002) is also confirmed by our regional climate model simulation. This is in contrast to present-day observations , where the southeastern parts of Vojvodina exhibit slightly more humid conditions, which is also reflected by the regional PI climate model simulation ( Fig. S1(d), (f), Fig. S2(d), (f)). A possible explanation for this can be given by considering the prevailing wind directions in South Banat in more detail. Fig. 6 shows the monthly frequency of the main simulated LGM wind directions (northwesterly and southeasterly) over the South Banat region. The simulated northwesterly winds (Figs. 4 and 5) agree with dune and sand formations in the northern part of the Carpathian Basin that were controlled by strong northerly/ northwesterly winds during glacial times (Sebe et al., 2011), whereas in the southeastern part the sand and loess deposits of the Banat Sands are driven by southeasterly winds . Southeasterly winds predominate throughout the year (see also Fig. S1(a)), except during summer (JJA), where their frequency in the South Banat region is lower than for northwesterly winds. In comparison with the present-day wind direction frequencies, where northwesterly winds dominate on an annual basis (Figs S1(b) and S2(b)), a strong increase in the frequency of southeasterly winds is simulated for LGM conditions, particularly during winter and spring (Fig. 6). We assume that during the LGM the prevailing northwesterly winds in summer (JJA) are not able to compensate for the annual lack of moisture induced by the increased frequency of the dry Košava wind from the southeast dominating the other seasons. Furthermore, the prevailing northerly wind directions in summer can explain the strong north-south gradient of temperature differences in July between the LGM and PI (Fig. 3). Since the northern parts of the study area are closer to the sphere of influence of the Fennoscandian ice sheet, the larger temperature differences between LGM and PI can be associated with the advection of cold air masses by the northerly winds.

Conclusion
The comparison of reconstructed and simulated July temperatures provides new insights into the regional climate of Vojvodina during the LGM (southeastern Carpathian Basin). The analysis of reconstructed July malaco palaeotemperatures and regional climate model data for LGM conditions obtained by the WRF model indicate that the malaco palaeotemperatures are most likely representative for the mean 2 m temperature. Since the reconstructed and simulated July temperatures are in good agreement, we assume that both approaches, the regional climate simulation and the malacothermometer calculations, provide acceptable results regarding the glacial temperature conditions in the considered area. However, when comparing the reconstructed and simulated temperatures, the different sources of uncertainty of the two approaches must also be considered. While the malaco temperatures are obtained by an implicit estimation from proxy data based on the subjective recognition of land snail assemblages, the simulated temperatures are determined by explicit (mathematical) calculations from regional climate model data. It is not possible to determine which of these two approaches provides better results since there are no temperature measurements for the LGM period with which to compare the reconstructed and simulated temperatures. In perspective, both approaches will progress in the future and may provide results closer to the hypothetically accurate/measured temperatures. In particular, the highresolution simulations offer a great potential to advance model-data comparisons in the future, due to further developments of physical parameterisations, numerical methods and computing power. Likewise, the development of new analytical methods and/or the refinement of the existing implicit methods will reduce the uncertainties of environmental reconstructions based on different sources of proxy data. Thus, a joined analysis of both methods opens new perspectives to obtain more reliable reconstructions of past climate states in the future.
Our results show that the application of a high-resolution regional climate model is necessary to reflect the variable regional environmental conditions as expressed by the reconstructed malaco palaeotemperatures. The resolution of the utilised global MPI-ESM model is too coarse to reflect the regional characteristics of the study area (Ludwig et al., 2019). The aridity conditions during summer (JJA/July), determined by two different methods (UNEP, 1992;De Martonne, 1925), are almost identical and show increased aridity from the northwestern to southeastern parts of the study area. This agrees with reconstructions for LGM conditions, e.g. by Krolopp and Sümegi (1995) and Sümegi and Krolopp (2002). Based on annual LGM climate data, the DM Aridity Index exhibits very humid conditions throughout the study area in contrast to semi-arid conditions proposed by the UNEP Aridity Index. Since these very humid conditions are not expressed by other proxy data for the study area (e.g. Marković et al., 2008Marković et al., , 2018, the DM Aridity Index should be considered with caution when applying to glacial climate conditions because it is predominantly defined to describe the aridity where mean air temperatures (T2m) are above −10°C (Equation 1; Gavrilov et al. 2020). The increase in aridity during JJA and July from the northwest toward the southeast is reversed in comparison with the present-day climate. The increased LGM aridity in the southeastern parts of the study area can be explained by more frequent (compared to today) dry winds from the southeastern quadrant during the year.

Supporting information
Additional supporting information can be found in the online version of this article.This article includes online-only Supplemental Data. Figure S1. UNEP Aridity Index based on precipitation and potential evapotranspiration and prevailing surface wind (10 m above ground; streamlines) for LGM (left column) and PI (right column) climate conditions based on regional climate (WRF) simulations. Panels (a, b) annual mean, (c, d) seasonal (JJA) mean, and (e, f) mean July climate and wind conditions. Figure S2. De Martonne Aridity Index (shaded) based on precipitation and temperature, and prevailing surface wind (10 m above ground; streamlines) for LGM (left column) and PI (right column) climate conditions based on regional climate (WRF) simulations. Panels (a, b) annual mean, (c, d) seasonal (JJA) mean, and (e, f) mean July climate and wind conditions.