Human-environmental interactions and seismic activity in a Late Bronze to Early Iron Age settlement center in the southeastern Caucasus

Long-term human-environmental interactions in naturally fragile drylands are a focus of geomorphological and geoarchaeological research. Furthermore, many dryland societies were also affected by seismic activity. The semi-arid Shiraki Plain in the tectonically active southeastern Caucasus is currently covered by steppe and largely devoid of settlements. However, numerous Late Bronze to Early Iron Age city-type settlements suggest early state formation between ca. 3.2-2.5 ka that abruptly ended after that time. A paleolake was postulated for the lowest plain, and nearby pollen records suggest forest clearcutting of the upper altitudes under a more humid climate during the Late Bronze/Early Iron Ages. Furthermore, also an impact of earthquakes on regional Early Iron Age settlements was suggested. However, regional paleoenvironmental changes and paleoseismicity were not systematically studied so far. We combined geomorphological, sedimentological, chronological and paleoecological data with hydrological modelling to reconstruct regional Holocene paleoenvironmental changes, to identify natural and human causes and to study possible seismic events during the Late Bronze/Early Iron Ages. Our results show a balanced to negative Early to Mid-Holocene water balance probably caused by forested upper slopes. Hence, no lake but a pellic Vertisol developed in the lowest plain. Following, Late Bronze/Early Iron Age forest clear-cutting caused lake formation and the deposition of lacustrine sediments derived from soil erosion. Subsequently, regional aridification caused slow lake desiccation. Remains of freshwater fishes indicate that the lake potentially offered valuable ecosystem services for regional prehistoric societies even during the desiccation period. Finally, colluvial coverage of the lake sediments during the last centuries could have been linked with hydrological extremes during the Little Ice Age. Our study demonstrates that the Holocene hydrological balance of the Shiraki Plain was and is situated near a major hydrological threshold, making the landscape very sensitive to small-scale human or natural influences with severe consequences for local societies. Furthermore, seismites in the studied sediments do not indicate an influence of earthquakes on the main and late phases of Late Bronze/Early Iron Age settlement. Altogether, our study underlines the high value of multi-disciplinary approaches to investigate human-environmental interactions and paleoseismicity in drylands on millennial to centennial time scales.

1 Introduction

The dynamics of human-environmental interactions on centennial to millennial timescales that are characterized by thresholds are a focus of geomorphological and geoarchaeological research (Dotterweich, 2008; Harden et al., 2014; Barton et al., 2016; Menges et al., 2019; Dong et al., 2020; Liu et al., 2021). Especially naturally fragile drylands (Fletcher et al., 2013) are ideal to investigate such long-term interactions: On the one hand, human societies in drylands are highly affected by even small-scale environmental variations (Schmidt et al., 2011; Balbo et al., 2016), and on the other hand, drylands sensitively react towards small-scale human disturbances (Neff et al., 2008; Suchodoletz et al., 2010; O’Henry et al., 2017; Rosen et al., 2019). Furthermore, given that many drylands are located in tectonically active regions, next to environmental and societal factors also destructive seismic events impacted the local prehistoric societies (Nur and Cline, 2000; Force, 2017; Lazar et al., 2020).

Located between Mesopotamia and the Eurasian steppes, the southern Caucasus (Georgia, Armenia, Azerbaijan) occupies a distinctive place in Old World archaeology (Smith, 2005; Manning et al., 2018). Given its location in the active continental Arabia-Eurasia collisional zone this region is characterized by strong seismic activity (Ismail-Zadeh et al., 2020), forming a natural hazard also for former regional societies (Varazanashvili et al., 2011). Here and in the neighboring regions of eastern Turkey and northwestern Iran, the Late Bronze to Early Iron Ages between about 3.5 and 2.8 ka were characterized by large fortified permanent settlements and an increasing sociopolitical hierarchy and complexity, leading to the development of the region’s first territorial polities with complex bureaucracies (Smith, 2005; Sagona, 2018; Erb-Satullo et al., 2019; Herrmann and Hammer, 2019). One important settlement center during that time was the Shiraki Plain in the semi-humid to semi-arid southeastern Caucasian lowlands that is largely devoid of settlements today (Figure 1). Here, several Late Bronze to Early Iron Age city-type fortified settlements of the so-called Lchashen-Tsitelgori tradition (Sagona, 2018) were found during the last years, suggesting early state formation since ca. 3.2 ka (Furtwängler et al., 1998; Maisuradze and Mindiashivili, 1999; Winter, 1999; Pitskhelauri et al., 2016; Bukhrashvili et al., 2019) (Figure 2). However, this phase of intensive settlement, maybe also extending over the Middle Iron Age (Manning et al., 2018), abruptly ended around 2.5 ka, and only minor human activity seemed to have continued afterwards (Arnhold et al., 2020). Based on destroyed wall structures and secondary buildings that were built on top of a former building layer in the intensively studied Didnauri settlement (Figure 2), it was suggested that the region was affected by strong seismic activity during the beginning of the Early Iron Ages (Pitskhelauri, 2018). The Shiraki Plain and its surroundings currently lack surficial water resources and mostly show a typical steppe vegetation dominated by grassland. However, a paleolake was suggested for the lowest part of the plain (Maruashvili, 1971), and pollen records from the nearby Udabno region (Figure 1) with similar climatic and environmental conditions and a similar Late Bronze/Early Iron Age settlement history (Kunze, 2017) suggest intensive forestation of the upper altitudes until the Late Bronze Age under a more humid climate compared with today (Kvavadze and Todria, 1992). Therefore, besides substantial variations of human activity this also suggests different environmental and hydrological conditions in this region during the past. However, unlike the well-studied humid western Caucasian lowlands and regions >1,000 m a.s.l. in the Lesser Caucasus mountain range (Connor and Sagona, 2007a; Connor et al., 2007b; Messager et al., 2013; Joannin et al., 2014; Connor et al., 2018), existing studies from the lower-altitude semi-humid to semi-arid south-eastern Caucasus region are limited to scattered pollen records with generally rather poor chronostratigraphic control (Gogichaishvili, 1984; Connor and Kvavadze, 2008), or to fluvial-geomorphological (Ollivier et al., 2015, 2016; Suchodoletz et al., 2015, 2018), paleoecological (Bliedtner et al., 2018, 2020a; Suchodoletz et al., 2020) or paleoclimatic (Bliedtner et al., 2020b) investigations of fluvial sediments. However, due to their discontinuous character fluvial sediment archives generally show a relatively poor chronological resolution (Lewin and Macklin, 2003). Therefore, Holocene human-environmental interactions in the semi-humid to semi-arid south-eastern Caucasus region are currently not well understood. Furthermore, beyond the historical period the possible impact of seismic events on former regional societies has not been studied so far (Varazanashvili et al., 2011).

FIGURE 1

FIGURE 1. The Caucasus region with the distribution of Late Bronze/Early Iron Age cultures (Sagona, 2018). The Shiraki Plain is shown with a red rectangle, and the Udabno region with similar ecologic and climatic conditions as well as a similar prehistoric settlement history with a red star. The topography is based on an SRTM DTM (http://www2.jpl.nasa.gov/srtm) [GE = Georgia, RU = Russia, AZ = Azerbaijan, AM = Armenia, TR = Turkey, IRN = Iran].

FIGURE 2

FIGURE 2. The Shiraki Plain with Late Bronze/Early Iron Age archaeological sites (Varazashvili and Pitskhelauri, 2011) and the core drillings and DGPS transects of this study. The area shown in Figure 3 is indicated with a yellow rectangle, and Lake Kochebi west of the plain with a blue rectangle [GE = Georgia, AZ = Azerbaijan].

During this multi-disciplinary study that was carried out on the semi-arid Shiraki Plain we combined geomorphological, sedimentological, chronological and paleoecological data as well as hydrological modelling to: (i) Reconstruct regional paleoenvironmental changes during the Holocene, (ii) identify possible natural versus human causes, and (iii) study possible seismic events during the Late Bronze/Early Iron Ages. Our study will help to better understand human-environmental interactions in the fragile drylands of the southeastern Caucasus especially during the Late Bronze to Early Iron Age period when early state formation occurred in this region.

2 Study area

The endorheic Shiraki Plain is located in eastern Georgia between ca. 41°18′ and 41°27′N, and 46°11′ and 46°30′E (Figure 2) at altitudes between about 550 and 650 m a.s.l., covering a catchment area of about 300 km². To the west, south, and north the plain is surrounded by low mountain ridges with highest altitudes of up to 970 m a.s.l. The plain forms part of the Kura Fold-and-Thrust-Belt, consisting of a series of south-vergent faults and thrusts that are composed of deformed Plio-/Pleistocene flysch to molasse sediments. These are formed by conglomerates, sands, loams, clays and sandstones (Gamkrelidze, 2003; Forte et al., 2010).

The regional climate is continental semi-arid. Mean temperatures range between ca. -4°C in January and ca. 23°C in July, and mean annual precipitation is around 490 mm (Furtwängler et al., 1998). The surface of the Shiraki Plain is covered by Vertisols (Matchavariani, 2019). The current xerophytic vegetation is characterized by steppe species such as Allium atroviolaceum, Muscari tenuiflorum and Puccinellia distans, and large parts of the plain are used for intensive agriculture of wheat and sunflowers (Furtwängler et al., 1998).

Permanent surface water is currently largely missing in the plain and on the surrounding mountain ridges, but seasonal creeks descend from the slopes. Due to this lack of permanent surface water resources the area is largely devoid of settlements today, and only in the northern part some larger villages are located. In contrast, west of the Shiraki Plain a seasonal salt lake (Lake Kochebi) is located at an altitude of ca. 780 m a.s.l. (Figure 2).

3 Methods

3.1 Field work

3.1.1 Stratigraphical mapping

To detect and map the extension of the formerly postulated paleolake (Maruashvili, 1971, 478f.), we applied vibracore drillings in the deepest parts of the plain using an Atlas Copco Cobra Pro hammer with a 60 mm diameter open corer: In a first step, 11 drillings were irregularly distributed over the topographically deepest part of the plain to verify and roughly outline the distribution of the possible paleolake sediments. In a second step, 10 further drillings were performed along a transect from the approximate center of the paleolake towards the northwest until the lacustrine sediments disappeared. Whereas the initial distance between the transect drillings in the approximate center of the paleolake was 200 m, to properly map the largest spatial extension of lacustrine sediments we reduced the drillings to distances of 100 m between drillings T6 and T8 at the northwestern transect margin (Figure 3). Supported by Ad-hoc (2005) and the Munsell soil color chart, the sediments were mapped according to basic field properties (approximate grain size, wet sediment color, structure, stone content, sedimentary/pedogenic horizons, carbonate content, hydroxymorphic features, occurrence of gypsum crystals). The grain size distribution was classified according to Jahn et al. (2006). To obtain a detailed stratigraphy of the lacustrine and their over- and underlying sediments, in the approximate center of the paleolake an artificial ca. 3 m deep trench was dug into NNE - SSW direction using a shovel excavator (Figure 3). Similar to the drillings, basic field properties of the sediments were described according to Jahn et al. (2006) and Ad-hoc (2005).

FIGURE 3

FIGURE 3. The lowest part of the Shiraki Plain with Late Bronze/Early Iron Age archaeological sites (Varazashvili and Pitskhelauri, 2011), core drillings, the studied trench and the maximal extension of the endorheic palaeolake based on the results of the occurrence of lake sediments that was extrapolated to the lowest part of the Shiraki Plain. However, given that lake sediments must have been deposited below a water column of unknown depth, the real maximal extension of the lake must have been somewhat larger than shown here.

3.1.2 Delineating the paleolake extension

First, we first built up a well-resolved digital terrain model (DTM) of the central Shiraki Plain by performing differential GPS (D-GPS) measurements with a resolution of some centimeters during September and November 2019 using a Stonex device operating with the Georgian CORS system (http://geocors.napr.gov.ge/SBC/spider-business-center). To cover a large area, these measurements were conducted along transects from a slowly moving vehicle. Later, these point measurements were interpolated in Arc Map GIS 10.4.1 using the “Topo To Raster” tool, and a raster DTM was constructed that was integrated into the existing SRTM DTM. Subsequently, we properly mapped the highest altitude of lake sediments that were detected in the drilling transect, and then extrapolated this altitude to the complete central part of the Shiraki Plain using the newly creating well-resolved DTM.

3.2 Laboratory analyses

3.2.1 Sedimentological analyses

Fourteen sediment samples were taken from the main sedimentological units in the trench for sediment analyses (Figure 4). Prior to the analyses the samples were air-dried, and material >2 mm was removed by sieving.

FIGURE 4

FIGURE 4. The stratigraphy of the trench with the locations of numerical dating, sediment, paleoecological and micromorphological samples. The luminescence results are given for feldspar (FS), quartz (QZ), fine grain (FG) and coarse grain (CG), and the preferred ages are underlined. Please note that the uppermost part of facies SH-E in the left-hand figure was not wetted, so that it appears lighter than the lower part of this material. The two trench walls are rectangular to each other. A version of this figure without any labellings and drawings can be found in Supplementary Figure S1.

Grain size was measured using 10 g of bulk sample material. After removing organic matter with 35% H₂O₂, the samples were dispersed in 10 ml 0.4 N sodium pyrophosphate solution (Na₄P₂O₇), and subsequently treated in an ultrasonic bath for ca. 45 min. The sand content (63–2000 µm) was determined by wet sieving, and the clay and silt content by measuring the material <63 µm with X-ray granulometry (XRG) using the SediGraph III 5120 (Micromeritics).

Calcium carbonate contents were determined according to Scheibler (Schlichting et al., 1995): Depending on the pre-tested approximate carbonate content, 1–10 g of material were filled into an Eijkelkamp Calcimeter apparatus, and 10% HCl was added continuously until the reaction ceased. Carbonate-bound C was calculated based on the CO₂ volume that was produced during the reaction.

Element distributions were measured on 32 mm-pellets that were produced by mixing and pressing 8 g of the ground material with 2 g CEREOX Licowax prior to measurement with a Spectro Xepos X-ray fluorescence spectrometer.

Measurements of mass-specific magnetic susceptibility (χ) were performed using a Bartington MS3 magnetic susceptibility meter equipped with a MS2B dual frequency sensor. After softly grounding and densely packing the material into plastic boxes, volumetric magnetic susceptibility was measured with a frequency of 0.465 kHz (κ_LF). Normalizing κ_LF with the sample mass yielded mass-specific magnetic susceptibility χ.

Total carbon and nitrogen were determined using a Vario EL cube elemental analyser on material that was ground in a vibratory mill for ca. 10 min (required grain size <30 µm). Total organic carbon (TOC) was calculated by subtracting inorganic carbon (calculated from carbonate-bound C) from total carbon.

Black carbon (BC) content and its aromaticity were determined by measuring benzene polycarboxylic acids (BPCA) contents and patterns, respectively, according to Glaser et al. (1998) with the first step modified by Brodowski et al. (2005). This method comprises metal removal with 4 M trifluoroacetic acid (Brodowski et al., 2005), followed by nitric acid digestion (170°C, 8 h under pressure), cation exchange chromatography, derivatization (trimethylsilylation) and gas chromatography with flame ionization detection (Glaser et al., 1998). The sum of the individual BPCAs was converted to the black carbon content by multiplication with 2.27 (conversion factor for charcoal; Glaser et al., 1998). The challenge of the samples under study was their high gypsum content, which interfered with BPCA analysis probably due to complexation with dissolved Ca during the procedure. The best compromise to solve this problem was to use smaller amounts of material for the analysis (100 mg instead of 500 mg).

3.2.2 Numerical dating

Numerical dating was carried out in the CEZ radiocarbon laboratory Mannheim (Germany).

3.2.2.1 Radiocarbon dating

For radiocarbon dating we used one bulk sediment sample (location see Figure 4). The sample was prepared with the Acid-Base-Acid method (Steinhof et al., 2017), and the bulk organic matter after removing NaOH-soluble organic matter was dated. The age was calibrated using the software SwissCal applying the Intcal20 calibration curve (Reimer et al., 2020).

3.2.2.2 Luminescence dating

For luminescence dating three samples were taken (locations see Figure 4) during night, and directly packed into light-proof plastic bags. To sample only material that was not influenced by sunlight, we removed the outer 20 cm of the material right before sampling. Sample preparation under subdued red light included sieving recovering material between 90 and 200 µm and <90 μm, following treatment with 10 and 30% HCl to remove carbonate, and with 10 and 37% H₂O₂ for about 14 days to remove organic matter. Only for sample MAL-10551 (55 cm) the coarse grain fraction 90–200 µm could be obtained in a datable amount. However, given a still very low amount of material we did not apply any density separation or etching to this fraction, resulting in a polymineral coarse grain fraction. Furthermore, for all samples the polymineral fine grain fraction 4–11 µm was separated from the sieve fraction <90 µm according to Stoke’s Law in Atterberg settling tubes. The luminescence measurements were carried out on a Risø TL-DA-20 luminescence reader equipped with a⁹⁰Sr/⁹⁰Y β-source (0.056 Gy/s for coarse grain, 0.074 Gy/s for fine grain). For the polymineral coarse (90–200 µm) and fine grain (4–11 µm) measurements applying the post-IR IR protocol of Buylaert et al. (2012) with a measurement temperature of 290 °C (pIRIR₂₉₀), we used a BG3/BG39 filter combination. We also aimed to compare feldspar and quartz fine grain measurement with respect to bleaching, since quartz bleaches faster than the feldspar pIRIR₂₉₀ signal (Kars et al., 2014). However, since also the fine grain fraction did not consist of enough material to etch away the feldspar, we measured the fine grain quartz signal by applying the single aliquot regeneration protocol of Murray and Wintle (2000) using an U340 filter that largely removes the feldspar emission. The a-values of polymineral and quartz fine grain were measured with an ²⁴¹Am α-source (0.115 Gy/s), and the a-value of 0.08 ± 0.03 for feldspar coarse grain was taken as a mean value from various literature about this kind of material (Wallinga et al., 2001; Rees-Jones and Tite, 2007; Kreutzer et al., 2018). Anomalous fading of the polymineral samples was tested by comparing the luminescence signals of not irradiated discs with those of irradiated discs directly after irradiation, and after seven and after 30 days following irradiation.

In parallel with the luminescence samples, ca. 500 g of material were taken from the same positions during daylight for dose rate determination. Given that samples MAL-10551 and MAL-10552 were located near stratigraphic borders, for both samples additional dose rate samples from facies SH-D were taken from nearby positions assumed to contribute ca. 30% of the gamma dose rate (see Table 1). After drying the samples at 105 °C for 24 h, dose rates were measured with low-level γ-spectrometry using a Canberra GCW4023 device. Given deep sediment desiccation for an unknown period after lake desiccation as could be recognized when opening the trench, a water content of 15 ± 10% was assumed to encompass the water contents of the complete burial periods. The cosmic dose rate was calculated according to Prescott and Hutton (1988), assuming a density of the overlying sediments of 1.6 g/cm³. The parameters used for age calculation and the luminescence ages are shown in Table 1.

TABLE 1

TABLE 1. Measurement parameters and results of luminescence dating from the trench.

3.2.3 Micromorphological analyses

One oriented block was taken for micromorphological analyses. After air-drying and impregnating with Oldopal P 80-21, the hardened block was cut and sliced into an upper (u) and a lower (l) 70 * 50 mm thin section (locations see Figure 4). Micromorphological description and analysis of the thin sections has been done under an Olympus BX51 petrographic microscope, following the guidelines of Stoops (2021).

3.2.4 Paleoecological analyses

Four sediment samples of about 1 kg were taken for paleoecological analyses (locations see Figure 4), and sieved with distilled water with mesh widths of 0.71 and 0.21 mm, respectively. Subsequently, organic material such as charcoal, bones, teeth or seeds was analyzed using an Amscope binocular (magnification = 40).

3.3 Hydrological modelling

Based on a DTM (SRTM C-band; Farr et al., 2007) we delineated the catchment boundaries using ArcGIS 10.3, defining the lowest point of the basin manually. Subsequently, we processed the catchment DTM using the freq function of the R raster-package of Hijmans (2020) to calculate the geomorphologically possible maximal lake size and volume (All calculations are available at Gihub: https://github.com/MLFischer/Shiraki-Plain-Water-Balance, 03/22).

Then, we used three Holocene time periods to calculate the water balance of the catchment. Paleoenvironmental and paleoclimatic data were taken from a pollen reconstruction that covers the period around the Late Bronze/Early Iron Ages in the Udabno region ca. 80 km to the west with very similar ecological and climate conditions (Kvavadze and Todria, 1992): (a) The Early/Mid Holocene with a mainly forest-covered basin above and fertile meadow grassland below the timberline at ca. 700 m a.s.l. For this period a similar annual temperature as today and a maximal precipitation amount as for period (b) were assumed. (b) The Late Bronze/Early Iron Ages with a clear-cut basin covered by fertile meadow grassland, a lake with a size of at least 5 km² (taken from this study), a reconstructed higher annual temperature of ca. 1.5 °C and higher precipitation of ca. 200–300 mm/a compared with today, and (c) a modern-day climate scenario with a drier grassland-covered basin and no lake. For the modern period, for comparison we also calculated the water balance for the altitude of current salt Lake Kochebi (located ca. 230 m higher than the lowest part of the Shiraki Plain) (location see Figure 2).

Based on major (paleo-)environmental properties as summarized in Supplementary Table S1, for each time period and land cover we used an established parametrization approach (e.g. Brutsaert, 1982; Bougealt, 1991; Blodgett et al., 1997; Bergner et al., 2003) to calculate the actual evapotranspiration rate (ETa). Climate variables such as air temperature, relative humidity, surface conditions (surface drag coefficient, albedo, emissivity and soil moisture availability) were parametrized following Bougealt, 1991, and using the paleoenvironmental reconstructions of Kvavadze and Todria (1992). Cloud coverage was achieved using a global cloud coverage dataset (Wilson and Jetz, 2016), and the cloud parameter was set according to Budyko (1974). Insolation was set according to Laskar et al. (2004). Modern-day climate variables were set using global gridded climate data (New et al., 2002). Temperature differences to the modern-day conditions were applied either using the paleoclimate reconstruction of Kvavadze and Todria (1992), or by using adiabatic average temperature gradients. Based on the precipitation and evapotranspiration rates we calculated the annual water flux volumes regarding the different land cover proportions within the basin.

4 Results

4.1 Field work

4.1.1 Stratigraphical mapping

To facilitate the stratigraphical description, we start here with the detailed description of the trench in the lowest part of the Shiraki Plain to establish the different facies, and then continue with the description of the drilling cores.

4.1.1.1 Trench

The facies described from bottom to top are shown in Figure 4:

- Facies SH-A (300 - 144 cm): The base of the trench is formed by brown (10YR4/3) carbonate-rich heavy clay with a polyedric structure. This material contains in-situ carbonate, gypsum, and iron oxide precipitations. Given the very fine grain sizes, we suggest that these deposits could represent limnic sediments of a paleolake. Between ca. 250 and 144 cm facies SH-A and SH-B are interfingered with each other, and the vertical cones in the interfingering zone showed widths of ca. 20 cm. We suggest that these deformation structures were linked with seismic shocks.

- Facies SH-B (144 - 135 cm): This facies consists of very dark gray (10YR3/1) heavy clay with a polyedric structure, and contains less gypsum precipitations and much less carbonate compared with facies SH-A. In some parts also reddish iron oxide stains were found. Given its high clay content and shiny slickensides on the aggregate surfaces, this material was interpreted as the Ahb horizon of a buried pellic Vertisol that had developed on the surface of facies SH-A and was overprinted by light-coloured gypsum precipitations. Given that it is not a surficial but buried soil, possibly formerly existing vertical cracks should not be recognizable any more.

- Facies SH-C (135 - 113 cm): Similar to facies SH-B this facies consists of black (7.5YR2.5/1) heavy clay with a fine polyedric structure. However, in difference to the former slight horizontal layering was recognizable. The matrix of this facies is mostly free of carbonate, and contains some mycelia-like gypsum precipitations and brownish iron oxide stains. This facies was interpreted as layered lacustrine sediments mostly consisting of eroded Vertisol material from the surrounding slopes (lake facies I).

- Facies SH-D (113 - 70 cm): This facies consists of dark reddish gray (5YR4/2), (dark) grayish brown (10YR4/2) and brown (7.5YR5/2) up to 1 cm thick evaporative carbonate and gypsum layers alternating with clastic layers of heavy clay that are broken by desiccation cracks. Numerous in-situ gypsum precipitations cause a pseudosand structure of the greyish evaporative layers, whereof the proportion increases upwards. This facies was interpreted as lacustrine sediments of a successively drying lake (lake facies II).

- Facies SH-E (70 - 0 cm): This facies unconformably overlies facies SH-D, and consists of carbonate-rich black (7.5YR2.5/1) heavy clay with a polyedric structure and shiny slickensides on the aggregate surfaces. This material contains carbonate and gypsum precipitations. It was interpreted as a vertisol-derived colluvial layer. In the WNW trench wall, at ca. 40 cm depth a ca. 10 cm thick layer of material of facies SH-D was intercalated. However, in the neighbouring NNE trench wall it could be recognized that this material forms part of a series of small inclined load structures that were detached from the upper part of facies SH-D and squeezed into the material of facies SH-E. We suggest that these inclined load structures were linked with seismic shocks. The upper part of facies SH-E was strongly disturbed by recent human activity (recent Ap horizon).

4.1.1.2 Drillings

The lowest parts of all drillings were formed by facies SH-A. However, with increasing distance from the former lake center such as in drillings T7 and T8 (Figure 3), the grain size of this facies coarsened and changed from heavy clay to silty clay. The intermixture of facies SH-A and SH-B (Figure 4) was found in all cores (Figure 5). Unfortunately, given the similar sediment properties of facies SH-B and SH-C these could not well be differentiated from each other in the cores. However, black-coloured material potentially originating from one of both facies was mostly found in the center of the plain, but not in cores TPCW, TPSS and T6 - T8. Overlying facies SH-D was detected in all cores but cores T7B and T8. Although facies SH-B was eroded in core T7B, some intermixing of this facies and facies SH-E was observed here. Facies SH-E was found in all cores. The core stratigraphies are listed in Supplementary Tables S2, 3, and the stratigraphy of the drilling transect is shown in Figure 5.

FIGURE 5

FIGURE 5. Stratigraphies of the drilling cores of the transect in the central Shiraki Plain (cores T1 - T8) with derived former maximal lake level that was based on the occurrence of lake sediments. However, given that lake sediments must have been deposited below a water column of unknown depth, the real maximal lake level must have been somewhat higher than shown here.

4.1.2 Delineating the maximal paleolake extension

To build up the well-resolved DTM of the central Shiraki Plain we measured 10,700 points with D-GPS. These points showed altitudes between 554 and 642 m a.s.l. The larger study area was covered with an irregular grid of measurement transects that followed existing earth roads, being denser in the central part of the plain where the study area is located (Figure 2). The maximal paleolake extension was derived by taking the altitude of core T7 that forms the northwesternmost drilling core of the transect towards the former northwestern lake margin where lake sediments (facies SH-D) were observed (555.34 m a.s.l.; Figure 5), and subsequently following this contour line in the newly created well-resolved DTM. This resulted in a maximal lake extent of ca. 6.2 km².

4.2 Laboratory analyses

4.2.1 Sedimentological analyses

The results of the sedimentological analyses are shown in Figure 6, and the data can be found in Supplementary Table S4.

FIGURE 6

FIGURE 6. Analytical results of the trench.

The grain size is largely dominated by clay, with irregularly varying values between 73 and 98% found throughout the sequence. The sand fraction strongly fluctuates between 0.5 and 21% within the sequence, showing the lowest values in facies SH-E. However, most of the sand was formed by in-situ precipitated gypsum crystals (pseudosand).

Calcium carbonate contents varied between 14 and 20% in facies SH-A. Facies SH-B showed a sharp drop to a value of 4%, and in facies SH-C carbonate contents dropped even more down to values between 0.4 and 0.2%. In facies SH-D carbonate values strongly increased to values between 20 and 11%, and similar values around 16% were also found in facies SH-E.

To analyze the inorganic element distributions, we first correlated the elements Si, Ti, Al, Fe, Rb, K, Zr, Ti, Mn, P, Ca, Mg and S in a correlation matrix. These elements are commonly included in element ratios that are used to trace grain sizes (Croudace et al., 2006), possible volcanic provenances (Zielhofer et al., 2017), paleo-redox conditions (Koinig et al., 2003) or human impact (Holliday and Gartner, 2007), or form the main components of silicates, gypsum and carbonates (Salminen, 2005). The results showed that the highest negative correlations exist between the siliciclastic elements Ti, Rb, Zr, Si and Al as well as Fe on the one hand, and Ca and S forming the main components of calcium carbonate and gypsum on the other hand (Supplementary Table S5). Therefore, we used the ratio Si/S to trace changes between siliciclastic layers and layers characterized by a high content of (evaporative) gypsum. This ratio shows high values between 117 and 40 in facies SH-A, an intermediate value of 20 in facies SH-B, slightly lower values between 13 and 10 in facies SH-C, an upwards decreasing trend from 6 to 1 in facies SH-D, and intermediate values between 20 and 26 in facies SH-E.

Relatively high values of mass-specific magnetic susceptibility (χ) between 0.31 and 0.18 * 10^-6 m³/kg were measured in facies SH-A, and a slightly lower value of 0.16 * 10^-6 m³/kg in facies SH-B. Uniform intermediate values of 0.11* 10^-6 m³/kg were found in facies SH-C, and in facies SH-D low values with an upwards decreasing trend from 0.08 to 0.03 * 10^-6 m³/kg were observed. Facies SH-E showed higher values around 0.25 * 10^-6 m³/kg.

Facies SH-A showed relatively low values of total organic carbon (TOC) between 0.3 and 0.5%. In contrast, facies SH-B showed a higher value of 1.0% that only slightly dropped to 0.9% in facies SH-C. Facies SH-D showed lower values between 0.4 and 0.7% that strongly increased to values between 3.1 and 3.3% in facies SH-E. Low TOC/N ratios between 3 and 6 in facies SH-A increase to values around 9 in facies SH-B and SH-C, subsequently decreasing to values between 3 and 7 in facies SH-D. In facies SH-E this ratio increases again to values around 11.

Black carbon (BC) contents ranged between 0.8 and 6.1 g kg⁻¹ in all facies but SH-D, where no black carbon could be detected. Black carbon contribution to TOC ranged between 19 and 37% in all but facies SH-D. The patterns of the individual benzenepolycarboxylic acids showed an equal contribution of benzenetetracarboxylic acids (B4CA), benzenepentacarboxylic acid (B5CA) and benzenehexacarboxylic acid (B6CA), demonstrating a similar high degree of condensation (Supplementary Table S6).

4.2.2 Numerical dating

4.2.2.1 Radiocarbon dating

The results of radiocarbon dating can be found in Table 2, and the age is shown in its stratigraphic position in Figure 4.

TABLE 2

TABLE 2. Results of radiocarbon dating from the trench.

4.2.2.2 Luminescence dating

Luminescence dating yielded different ages for different grain sizes and minerals, ranging between 11.6 ± 0.7 ka and 0.32 ± 0.09 ka. However, the ages of the three samples are in the stratigraphic order. The results of luminescence dating can be found in Table 1, and the ages are shown in their stratigraphic positions in Figure 4.

4.2.3 Micromorphological analyses

Lower sample (l) and approximately the lower half of upper sample (u) were taken from facies SH-C, and showed a blocky structure that was highly separated by fissures. This structure was moderately bioturbated, and most pores were occupied by lenticular gypsum crystals of different sizes that showed frequent inclusions. The birefringence fabric of the micromass is porostriated, i.e. shows stress features along the fissures what reflects some swelling-shrinking processes linked to argilloturbation (Figure 7A). Due to abundant very small organic particles the micromass is speckled. Upwards, organic material increases as layered amorphous intercalations (Figure 7B). Pedofeatures not related to gypsum precipitation include redoximorphic phenomena such as nodules of Fe- and Mn- oxi-hydroxides as well as clusters of small spheres of Fe-oxides after pyrite framboids (Figure 7C). The latter are often located next to organic residues, and some of those spheres still keep some unaltered pyrite inside (Figure 7D). The redoximorphic pattern can be classified as stage F or G of Vepraskas et al. (2018), meaning extreme, permanent reduction conditions.

FIGURE 7

FIGURE 7. Lower sample (l): (A) Porostriation along fissures (arrows) resulting from stress (argilloturbation) due to swelling-shrinking. Some gypsum lenses occur at the top as loose discontinuous infillings. (B) Upper part of the reduced material, with layered sorting of particle sizes and subhorizontal organic intercalations. (C) Hypocoating of Fe-oxihydroxides along a pore (dark line) that has been broken by swelling-shrinking (argilloturbation) and is seen now as infilling of fragments together with gypsum lenses. (D) Hypocoatings of Fe-oxihydroxides around pores, and loose infilling of iron spheres, pseudomorphs after pyrite framboids. Some of the pseudomorphs show luster under incident light, which indicates that pyrite oxidation has not been complete. Upper sample (u): (E) Dark inclusions in the groundmass corresponding to particulate organic matter. (F) Loose discontinuous infilling of Fe nodules, pseudomorphs after pyrite framboids developed on organic materials. (G) Geodic nodules of sparite. (H) Loose continuous infilling of lenticular gypsum in a biopore (note the root section) together with an impregnative hypocoating of Fe-oxyhydroxides.

The upper part of sample (u) was taken from facies SH-C. This sample is composed of a 0.5 cm thick layer of lenticular gypsum crystals that underlies clayey material with a lighter color compared with the clayey material of sample (l), but that also contains microparticles of organic matter (Figure 7E). Fe-oxide pseudomorphs after pyrite are equally present (Figure 7F), and the redoximorphic pattern corresponds to stage B of Vepraskas et al. (2018) indicating moderate reduction conditions. As the most striking pedofeatures of this sample we find micrite and geodic sparitic nodules of probably biogenic origin that could be derived from algae or other organisms living in shallow water conditions (Freytet and Verrecchia, 2002) (Figure 7G). In particular, infillings of calcitic excrements in channels are difficult to explain, unless they were made by fauna mixing materials with different calcite contents (Figure 7G).

Scans of the thin sections of both samples can be found in Supplementary Figure S2.

4.2.4 Paleoecological analyses

Different types of organic material were detected in the samples, i.e. charcoal, plant or wood tissues, plant roots, probably fish bones that could not be assigned to certain species, and in samples I and II fish teeth of the family Cyprinidae (sample locations see Figure 4).

The results are listed in Supplementary Table S7, and a photo of fish teeth and probably fish bones from sample II is shown in Supplementary Figure S3.

4.3 Hydrological modelling

The minimal Holocene lake size of ca. 6.2 km² that was calculated during this study should have hosted a water volume of ca. 0.006 km³. This water volume is much lower compared with the value of ca. 0.75 km³ for a lake with a size of about 70 km² that would have formed at the altitude of the theoretical modern-day outflow of the plain at ca. 577.5 m a.s.l. However, according to our data that water level was never reached during the Holocene, i.e. the Shiraki Plain was always endorheic during that period.

For the Early/Mid Holocene we assumed a similar temperature as today and calculated an ETa of 481–565 mm a⁻¹ on fertile meadow grassland and of up to 1,092–1,648 mm a⁻¹ over forests, and for the Late Bronze/Early Iron Ages with a reconstructed higher temperature of ca. 1.5°C (Kvavadze and Todria, 1992) an ETa of 517–607 mm a⁻¹ over fertile meadow grassland and of 822 mm a⁻¹ over the paleolake. For the modern-day climate we calculated an ETa over drier grassland of 429–508 mm a⁻¹, and for the conditions at the ca. 230 m higher altitude of current salt lake Kochebi this results in a slightly lower current drier grassland ETa of 400–475 mm a⁻¹.

For the Early/Mid Holocene with an assumed maximal precipitation amount of 700–800 mm a⁻¹ (0.21–0.24 km³ a⁻¹) and the timberline located at 700 m a.s.l we calculated an ETa volume over forests of 0.05–0.08 km³ a⁻¹, and for the fertile meadow grassland below the timberline an ETa volume of 0.12–0.14 km³ a⁻¹. This would result in a balanced water budget of -0.02–0.06 km³ a⁻¹ (Figure 8A). For the Late Bronze/Early Iron Ages with a reconstructed precipitation amount of 700–800 mm a⁻¹ (0.21–0.24 km³ a⁻¹; Kvavadze and Todria, 1992), for a lake size of 5–10 km² we calculated an ETa volume of 0.004–0.008 km³ a⁻¹, and a fertile meadow grassland ETa (covering the rest of the basin) of 0.16–0.18 km³ a⁻¹. This results in a positive water budget between 0.02 and 0.08 km³ a⁻¹ (Figure 8B). For the modern-day time with a precipitation amount of about 500 mm a⁻¹ (0.15 km³ a⁻¹) and a drier grassland evapotranspiration of 0.13–0.15 km³ a⁻¹, this results in an overall balanced water budget of -0.002–0.022 km³ a⁻¹, and for the climatic conditions at the ca. 230 m higher altitude of current salt lake Kochebi the model shows a slightly positive water budget of 0.007–0.029 km³ a⁻¹ (Figure 8C).

FIGURE 8

FIGURE 8. Water balance model of the Shiraki Plain for three different scenarios during the Holocene based on different annual precipitation and evapotranspiration values. The resulting annual water balances are shown in the lower parts in blue: (A) Early/Middle Holocene. (B) Late Bronze/Early Iron Age. (C) Today. For the latter period, for comparison the current water balance for the altitude of recent salt lake Kochebi west of the Shiraki Plain is shown (location see Figure 2).

5 Discussion

5.1 Holocene landscape evolution of the Shiraki Plain

Based on our stratigraphical, chronological, sedimentological, micromorphological and paleoecological data we reconstructed the Holocene landscape evolution of the Shiraki Plain:

i) During the Early/Mid-Holocene no lake existed in the central part of the Shiraki Plain (Figure 9A). Instead, the clay-rich material parent material (facies SH-A) was overprinted by Vertisol formation (facies SH-B) what is also confirmed by our analytical data: Depending on the local conditions Vertisols show strongly varying TOC contents that also encompass the measured value of ca. 1%, and similar TOC/N ratios around or slightly higher than 10 were also reported from other Vertisols (Virmani et al., 1982; de la Rosa et al., 2008). BC was hardly studied in Vertisols so far. However, according to Reisser et al. (2016), next to the contents of organic carbon and clay one factor that strongly determines the soil BC dynamics is the regional climate. Accordingly, our measured proportion of BC of TOC of about 30% and the equal proportions of the benzenepolycarboxylic acids B4CA, B5Ca and B6Ca (Figure 6) are typical for Chernozems (Rodionov et al., 2010) that are found in regions with climate conditions that are similar as those of the Shiraki Plain (Strouhalová et al., 2019). No numerical ages are available for the start of Vertisol formation. However, the bulk sediment radiocarbon age of 5.3 - 5.0 cal ka BP from facies SH-B, averaging the age of organic matter in the soil matrix (Wang and Amundson, 1996), demonstrates that Vertisol formation must have lasted at least until that time. Looking at the luminescence ages, fine-grained feldspar pIRIR₂₉₀ dating gave a minimum age of 11.6 ± 0.7 ka (no fading correction applied), and fine-grained quartz an age of 7.1 ± 0.3 ka. Soil turbation processes cause the input of bleached mineral grains into the soil matrix also after sediment deposition (Bateman et al., 2003). Hence, in difference to the radiocarbon ages dating the soil organic matter, the luminescence ages should reflect argilloturbation linked with Vertisol formation. The large difference between both luminescence ages can possibly be explained with the much better bleachability of quartz compared with the pIRIR₂₉₀ feldspar luminescence signal (Kars et al., 2014). The cone-like intercalation of facies SH-A and SH-B indicates a strong seismic event when the Vertisol (SH-B) was already well developed, i.e. this event must have occurred during the Mid-to Late Holocene (Figure 9B).

ii) A permanent lake must subsequently have covered the pellic Vertisol in the central part of the Shiraki Plain. The first lake phase is reflected by horizontally layered facies SH-C (Figure 9C). Besides the layering, its (permanent) lacustrine character is demonstrated by micromorphology that showed pyrite formation and redoximorphic stages F or G (Figures 7C,D), indicating permanent reduction conditions (Vepraskas et al., 2018). These must have been caused by longer-lasting permanent covering of the sediments by lake water. The sedimentological properties of facies SH-C with a TOC content around 0.9%, the grayish to black color, TOC/N ratio around 9, equal proportions of the benzenepolycarboxylic acids B4CA, B5Ca and B6Ca and proportions of BC of TOC of ca. 33 to 22% are similar with those of the underlying Vertisol horizon (facies SH-B; Figure 6). This indicates the origin of facies SH-C from eroded Vertisols around the lake. Lower values of magnetic susceptibility compared with facies SH-B could possibly be explained with partial destruction of the magnetic signal under anoxic conditions (Hanesch and Scholger, 2005) (Figure 6). The origin of facies SH-C from eroded Vertisols is also supported by the observation that the Vertisol horizon (facies SH-B) was not detected between transect cores T6 and T8 located towards the former lake margin: The flat topography and rapid covering by lake water should have protected the Vertisol against soil erosion in the central plain, whereas it was eroded from the not water-covered and steeper parts of the plain beyond the paleolake shore (Figure 5). Fine-grained feldspar pIRIR₂₉₀ luminescence dating gave a minimum age of 7.4 ± 0.5 ka (not corrected for anomalous fading) for facies SH-C, and fine-grained quartz an age of 4.3 ± 0.3 ka. Given the better bleachability of quartz compared with the pIRIR₂₉₀ signal (Kars et al., 2014), the quartz age probably best approximates the depositional age of the lake sediments. Archaeological finds show intensive settlement of the Shiraki region during the Late Bronze to Early Iron Ages since ca. 3.2 ka (Furtwängler et al., 1998, 1999; Maisuradze and Mindiashivili, 1999; Pitskhelauri et al., 2016; Bukhrashvili et al., 2019). Given that no large-scale human settlement is known from the Shiraki Plain prior to that time, the Vertisols must have been eroded by agricultural activity of that culture. Hence, on the one hand the chronological difference between the luminescence-dated start of Vertisol erosion reflected by the deposition of black-colored lacustrine facies SH-C, and the start of intensive human activity some centuries later, can possibly be explained by so far incomplete numerical dating of the archeological finds. Thus, further numerical dating of archaeological sites could result in older ages. On the other hand, this difference could also be explained by insufficient bleaching of the fine-grained vertisol-derived material during its colluvial transport from the slopes towards the central lake basin (Fuchs and Lang, 2008).

iii) After an unknown period of time the lake must have dried out. This process is indicated by facies SH-D formed by the alternation of clastic clayey with evaporitic carbonate and gypsum layers (Figure 9D). Accordingly, unlike lower lake facies SH-C with <1% carbonate facies SH-D contains ca. 10–20% carbonate. Likewise, increasing sand contents up to >21% can be attributed to gypsum crystals (pseudosand) that were not destroyed during preparation of the grain size samples (Figure 6). In facies SH-D the micromorphological analyses show only moderate reduction conditions as indicated by stage B of the redoximorphic pattern (Vepraskas et al., 2018). Furthermore, shallower lake conditions compared with facies SH-C are also indicated by micrite and geodic sparitic nodules of probably biogenic origin likely derived from algae or other organisms living in shallow water conditions (Freytet and Verrecchia, 2002), and regular desiccation of the lake bed is indicated by numerous desiccation cracks within the clayey layers. The upwards increasing proportion of evaporitic carbonate and gypsum compared with clastic clayey material probably reflects an intensified desiccation trend with time. This trend is also indicated by upwards decreasing Si/S ratios and decreasing values of mass-specific susceptibility, reflecting the relative increase of S-containing gypsum and parallel decrease of siliciclastic material towards the top, respectively (Figure 6). Fish teeth of the family Cyprinidae (including species such as carps and minnows) were found in upper facies SH-D. These species are not seasonal, and require fresh or at least brackish water conditions (Kottelat and Freyhof, 2007). Hence, also during the final period of lake desiccation the lake must have had longer stable phases with fresh or brackish water, allowing the establishment of stable fish populations that could potentially have been used by the local populations. This also demonstrates that the desiccation process must have continued over a longer period. Going from the central plain towards its margin, facies SH-C was detected up to transect core T7 (Figure 5). Based on its altitude we derived a maximal lake extent of ca. 6.2 km² and maximal depth of >1.5 m (Figures 3, 5). However, given that deposition of facies SH-C must have been linked with a water column of unknown depth, these calculated maximal values should actually be minimum values. There are no numerical ages for the end of lake desiccation, however, the luminescence sample taken from overlying facies SH-E gave a fine-grained feldspar pIRIR₂₉₀ age of 2.4 ± 0.2 ka, a fine-grained quartz age of 0.5 ± 0.05 ka and a coarse-grained feldspar pIRIR₂₉₀ age of 0.3 ± 0.1 ka. Both feldspar pIRIR₂₉₀ ages did not show anomalous fading. The older fine-grained quartz age compared with the coarse-grained feldspar pIRIR₂₉₀ age probably indicates that slower bleaching of the latter signal compared with quartz (Kars et al., 2014) must have been counterbalanced by the better bleachability of coarse compared with fine grains (Olley et al., 1998). However, it is also possible that due to argilloturbation, causing the input of younger material, these youngest ages are to some degree underestimated (Bateman et al., 2003). Hence, the youngest coarse-grained feldspar pIRIR₂₉₀ age of 0.3 ± 0.1 ka should represent a terminus ante quem for lake desiccation. Accordingly, also historic maps from the 18th century AD (= ca. 0.3 ka) do not show a lake or swamp in the central Shiraki Plain (Vakhushti Bagrationi Institute of Geography Tbilisi, 1997).

iv) Finally, during the last centuries the dried paleolake was covered by a colluvium (facies SH-E) (Figure 9E), although on the current soil map this material is shown as part of the regional surficial Vertisols (Matchavariani 2019). Besides its black color also the TOC/N-ratio of ca. 10-11 is similar with the values for the buried Vertisol (facies SH-B), supporting an origin of this material from eroded surrounding Vertisols. Significantly higher TOC contents of >3% and higher values of the magnetic susceptibility compared with facies SH-B can possibly be explained with the longer time of Vertisol formation and current agriculture, causing continued input of organic material and magnetic enhancement (Jordanova and Jordanova, 1999). However, the lower proportion of BC of TOC of about 16% should indicate a dilution of BC by other kinds of organic material. Given the large absence of human settlements following the Late Bronze/Early Iron Ages (Arnhold et al., 2020), definite causes for the deposition of the colluvial layer cannot be given here. However, a strong phase of fluvial activity was detected in the regional rivers for a later phase of the Little Ice Age between 0.5 and 0.4 ka (Suchodoletz et al., 2018). This phase must have been linked with at least seasonal stronger discharge caused by more intensive precipitation compared with today. Although the total precipitation amount was obviously not enough to cause renewed longer-lasting lake formation in the Shiraki Plain, we suggest that due to preceding anthropogenic deforestation extreme precipitation events during that time could have caused surface erosion and colluvial deposition in the semi-arid landscape of the central Shiraki Plain also without direct human influence. Following its deposition, colluvial facies SH-E was interfingered with underlying facies SH-D by a strong seismic event in the form of inclined load structures (Figure 9F).

FIGURE 9

FIGURE 9. Landscape evolution of the Shiraki Plain: (A) Early/Middle Holocene. (B) Middle/Late Holocene. (C) Late Bronze/Early Iron Ages. (D) Later part of Late Bronze/Early Iron Ages and/or following period. (E) Last centuries–I. (F) Last centuries–II.

5.2 Changes of the hydrological balance between natural and human controls

Hydrological modelling helped to retrace the observed changes of the regional hydrological balance of the Shiraki Plain, and to identify possible natural and human drivers:

(i) From the neighbouring Udabno region showing comparable ecological and climatic conditions, we know that prior to the Late Bronze/Early Iron Ages the slopes above 700 m a.s.l. have been intensively forested, and that species-rich meadows were found below 700 m a.s.l. (Kvavadze and Todria, 1992). The detailed regional climate evolution of that period is not known, so that we applied the maximal precipitation amount of 700–800 mm that was reconstructed for the Late Bronze/Early Iron Ages and the current annual temperature. Given that at least parts of that period in the southern Caucasus were warmer compared with today (Connor and Sagona, 2007a) and drier compared with the Late Bronze/Early Iron Age (Connor and Kvavadze, 2008; Bliedtner et al., 2020b), the calculated water balances for that time should represent maximum values. Using those parameters, we obtained a balanced instead of a positive water budget that would be required for lake formation. Hence, the observed existence of a terrestrial Vertisol and not a lake in the central Shiraki Plain for this time could be supported by hydrological modelling.

(ii) In the Udabno region the forests were largely clear-cut during the Late Bronze/Early Iron Ages, and annual temperatures were up to 1.5 °C and annual precipitation 200–300 mm higher compared with today (Kvavadze and Todria, 1992). Using those parameters also for the Shiraki Plain that was settled by the same cultures, we obtained a positive water budget that even remained positive when the reconstructed paleolake surface with a higher ETa was taken into account. Thus, observed lake formation, indicating a significant change of the regional water budget during that time, could be retraced by hydrological modelling. Since we did not change the annual precipitation amount compared with the period before, and even used a higher annual temperature leading to higher ETa values for hydrological modelling, the change towards a positive water budget must have been caused by the large-scale human-induced regional forest clearcutting that reduced regional

(iii) Using the current climate and the drier grassland vegetation of the Shiraki Plain for hydrological modelling, we again obtained a balanced water budget. This could explain the observed lake desiccation. Since grassland vegetation in the whole catchment has probably persisted since clear-cutting during the Late Bronze/Early Iron Ages, that change of the water budget must have been linked with the regional aridification trend that was observed since that time (Kvavadze and Todria, 1992; Bliedtner et al., 2020b). Interestingly, when using the present environmental conditions of the ca. 230 m higher altitude where the current salt lake Kochebi is located, we obtained a slightly positive water balance what could possibly explain the existence of that lake. Although next to climatic parameters lake formation is also influenced by hydrogeological factors (winter, 1999), however, two other similar salt lakes in the Gareja region are located at about 760 m a.s.l. (Jikurebi lake) and at about 820 m a.s.l. (Kopatadze lake), respectively. Hence, the existence of those lakes seems to confirm that the current water budget of the Shiraki Plain is very close to positive values, i.e. that also smaller-scale changes of the controlling factors could currently cause a shift towards a positive water budget allowing lake formation.

However, limitations of such kind of hydrological modelling remain, as unquantified processes may affect the water balance as well. These include e.g. cloud cover changes, whereof no paleo information is available. Further unquantified processes could be seasonality changes, rapid snow melt and/or intensified run-off within the basin. Additionally, specific values for a class-specific ETa may differ interannually, as the vegetation reacts in a dynamic way to the available moisture. Therefore, our ETa estimates represent the best-possible long-time averages using this method.

Altogether, our hydrological modelling suggests that the Holocene hydrological balance of the Shiraki Plain was and is situated near a major hydrological threshold (White, 2019). That means that also relatively small-scale human or natural changes of the regional water balance caused a change from balanced/negative to positive values or vice versa, with significant geoecological changes due to lake formation or desiccation in the central Shiraki Plain. These changes should also strongly have influenced the living conditions of former local societies by supply or removal of water and fish resources.

5.3 Paleoseismic activity

The Shiraki Plain forms part of the NW-SE trending Kura Fold-and-Thrust-Belt which accommodated ∼25 km of shortening (Forte et al., 2010). While thrust front advance is generally towards the SW defining the southern margin of the Shiraki Plain, Forte et al. (2010) postulated the presence of a backthrust at its northern margin. According to these authors, the neotectonic expression of the latter structure in the field suggests its relatively young age. It follows that SW-directed thrust belt dynamics, including backthrusting towards the NE, controlled the formation of the endorheiic Shiraki Plain. Historic as well as recent seismic activity (Varazanashvili et al., 2011; Ismail-Zadeh et al., 2020; Kldiashvili, 2021) suggests that these processes are still ongoing. The latter authors documented three large earthquakes in the Alaverdi region located ca. 100 km to the northwest during the last 500 years. These occurred in 1530 (reconstructed magnitude 5.7), between 1667 and 1668, and in 1742 AD (reconstructed magnitude 7).

Seismites are deformational structures attributable to seismic shocks, and classification schemes including tentative links to earthquake magnitudes were proposed by Rodriguez-Pascua et al. (2000) and Moretti and Sabato (2007). During the last years such structures in lacustrine sediments received growing interest as archives for paleoseismic activity (Archer et al., 2019). However, many of such studies have focused on lakes in high-relief regions that are also susceptible to slope failures not related to earthquakes, rendering the seismicity interpretation more difficult. In contrast, the stratigraphic sequence of the Shiraki region records deposition in a relatively low relief region and should therefore provide an ideal archive for past seismicity. The geometries observed at the base of the trench walls (Figure 4) indicate an interfingering of facies SH-A and SH-B in the form of vertical cones, and are interpreted as load structures (Moretti and Sabato, 2007). Given that the internal layering is difficult to recognize, these structures may represent pillow structures according to Rodriguez-Pascua et al. (2000) which these authors tentatively linked with an up to magnitude 7 earthquake. As shown in Figure 4, these cones are unconformably overlain by lake facies SH-C and SH-D. Within the uppermost part of SH-D, a series of inclined load structures, verging towards WNW, are observed (Figure 4). A similar situation is also interpreted from cores TPNW and TPE, since material of facies SH-D was found within facies SH-E here (Supplementary Table S2; location of the cores see Figure 3). The absence of internal layering hinders a further definition of these soft-sediment deformation features. However, we note that the load structures between SH-A and SH-B are vertical, whereas a sense of shear with top to the WNW-direction must have been present during the formation of the inclined load structures in SH-D and SH-E (Supplementary Figure S1). The strike of these inclined load structures is at a right angle to the main tectonic transport direction of the Greater Caucasus with top to the SW-direction and available fault plane solutions (Forte et al., 2010; Ismail-Zadeh et al., 2020). According to the latter authors, all active faults in the Shiraki Plain trend NW-SE and no perpendicular accommodation structures have been mapped yet (Supplementary Figure S4).

The above observations indicate the presence of two distinct seismites. These are separated from each other by at least one angular unconformity that is represented by lake facies I (SH-C) and the lower part of lake facies II (SH-D), and thus record two individual earthquakes. Based on the available chronological data the first earthquake must have occurred when the Vertisol, represented by facies SH-B, was already well developed but vertisol-derived lake facies I (SH-C) was still not deposited. Hence, the first event must have occurred during the Middle to Late Holocene prior to the intensive phase of Late Bronze/Early Iron Age settlement with large-scale soil erosion. The second earthquake must have occurred when the colluvial layer (SH-E) already covered the paleolake sediments (facies SH-D), i.e. during the last centuries. However, given the available dating we cannot safely attribute the latter seismic event to a certain historic earthquake mentioned above (Varazanashvili et al., 2011; Kldiashvili, 2021). Altogether, based on our available data the two recorded strong seismic events should not have occurred during the main and late phases of the Late Bronze/Early Iron Age settlement period. Therefore our data do not support the hypothesis of Pitskhelauri (2018), who suggests strong seismic activity in the Shiraki Plain during the start of the Early Iron Age.

6 Conclusion

Using a multi-disciplinary approach, during our study in the semi-arid Shiraki Plain we reconstructed Holocene human-environmental interactions in an important Late Bronze/Early Iron Age settlement center in the southeastern Caucasus. Our data show a balanced to negative water balance during the Early and Middle Holocene, so that no lake could develop in the lowest part of the plain but a terrestrial vertisol had formed instead. It is very likely that similar with the neighbouring Udabno region forestation of the higher altitudes contributed to the balanced to negative water balance. Subsequently Late Bronze/Early Iron age human activity since at latest 3.2 ka obviously changed the regional water balance by forest clearcutting, leading to the formation of a lake with a size of several square kilometres in the lowest part of the plain. However, the lake level never reached the necessary altitude for an outflow so that the lake always remained endorheic. Strong human impact on the landscape during that intensive settlement period is also demonstrated by the deposition of vertisol-derived material on the lake bottom resulting from soil erosion processes. Subsequently the lake dried out during a longer period, probably caused by a regionally observed aridification trend. Given that freshwater fishes could thrive in the lake even during parts of the desiccation period, the lake potentially offered valuable ecosystem services for the regional prehistoric societies of that period. During the last centuries the lake sediments were covered by a vertisol-derived colluvium, and given the absence of large-scale human activity between the Early Iron Ages and the recent period its deposition must have been linked with hydrological extremes during the Little Ice Age. Also seismites were detected in the studied sediments that were probably related with two individual earthquakes. Given that one event must have occurred prior to the intensive phase of Late Bronze/Early Iron Age settlement, and the second after large-scale colluvial deposition during the last centuries, strong earthquakes did obviously not affect the main and late phases of Late Bronze/Early Iron Age settlement in the Shiraki Plain.

Our study demonstrates that the Holocene hydrological balance of the Shiraki Plain was and is situated near a major hydrological threshold. Consequently, also relatively small-scale human or natural influences on this semi-arid landscape can cause significant geoecological changes that strongly affected the living conditions of local societies due to lake formation or desiccation in the central part of the plain. Hence, the results of our study underline the high fragility of drylands towards also small-scale external perturbations, and demonstrate the high value of multi-disciplinary approaches to investigate their long-term evolution on millennial and centennial time scales.

Data availability statement

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

Author contributions

HvS—study design, stratigraphical field work and analysis, data compilation, writing; GK—D-GPS and GIS analyses and mapping, writing; TK—paleoecological analyses, writing; MLF—hydrological modelling, writing; RMP—micromorphological analyses, writing; AK—stratigraphical field work and analysis, writing; BS—laboratory analyses, writing; BG—laboratory analyses, writing; SL—numerical dating, writing; SH—evaluation of paleoseismicity, writing; AS—stratigraphical field work and laboratory analyses; LN—stratigraphical field work, writing; ML—stratigraphical field work, writing; MA—archeological field work and data compilation; LL—archeological field work and data compilation, writing; ME—study design, data compilation, writing.

Funding

This study was financially supported by the Georgian National Science Foundation (SRNSF) (project number FR-18-22377).

Acknowledgments

We thank Ulrich Göres (Dresden) for his support during field work. Furthermore, we are indebted to Katja Pöhlmann (Leipzig) and Heike Maennicke (Halle) for their help during the laboratory analyses. The authors acknowledge support from the Open Access Publishing Fund of Leipzig University supported by the German Research Foundation within the program Open Access Publication Funding.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

Ad-hoc, A. G. B. (2005). Bodenkundliche Kartieranleitung. Hannover: Federal Institute for Geosciences and Natural Resources in collaboration with the State Geological Services. 5th ed.

Google Scholar

Archer, C., Noble, P., Rosen, M. R., Sagnotti, L., Florindo, F., Mensing, S., et al. (2019). Lakes as paleoseismic records in a seismically-active, low-relief area (Rieti Basin, central Italy). Quat. Sci. Rev. 211, 186–207. doi:10.1016/j.quascirev.2019.03.004

CrossRef Full Text | Google Scholar

Arnhold, S., Bukhrashvili, P., Fassbinder, J., Tskvitinidze, Z., Abele, J., and Davitashvili, S. (2020). Untersuchungen in Samreklo 2019. Mittl. Dtsch. Orient-Gesellschaft Berl. 152, 111–123.

Google Scholar

Balbo, A. L., Gómez-Baggethun, E., Salpeteur, M., Puy, A., Biagetti, S., and Scheffran, J. (2016). Resilience of small-scale societies: A view from drylands. Ecol. Soc. 21, 53. doi:10.5751/ES-08327-210253

CrossRef Full Text | Google Scholar

Barton, M. C., Ullah, I. I. T., Bergin, S. M., Sarjoughian, H. S., Mayer, G. R., Bernabeu-Auban, J. E., et al. (2016). Experimental socioecology: Integrative science for Anthropocene landscape dynamics. Anthropocene 13, 34–45. doi:10.1016/j.ancene.2015.12.004

CrossRef Full Text | Google Scholar

Bateman, M. D., Frederick, C., Jaiswal, M. J., and Singhvi, A. K. (2003). Investigations into the potential effects of pedoturbation on luminescence dating. Quat. Sci. Rev. 22, 1169–1176. doi:10.1016/s0277-3791(03)00019-2

CrossRef Full Text | Google Scholar

Bergner, A. G., Trauth, M. H., and Bookhagen, B. (2003). Paleoprecipitation estimates for the Lake Naivasha basin (Kenya) during the last 175 ky using a lake-balance model. Glob. Planet. Change 36, 117–136. doi:10.1016/s0921-8181(02)00178-9

CrossRef Full Text | Google Scholar

Bliedtner, M., Suchodoletz, H. v., Schäfer, I., Welte, C., Salazar, G., Szidat, S., et al. (2020a). Age and origin of leaf wax n-alkanes in fluvial sediment-paleosol sequences, and implications for paleoenvironmental reconstructions. Hydrol. Earth Syst. Sc. 24, 2105–2120. doi:10.5194/hess-24-2105-2020

CrossRef Full Text | Google Scholar

Bliedtner, M., Zech, R., Kühn, P., Schneider, B., Zielhofer, C., and Suchodoletz, H. v. (2018). The potential of leaf wax biomarkers from fluvial soil-sediment sequences for paleovegetation reconstructions - upper Alazani River, central southern Greater Caucasus (Georgia). Quat. Sci. Rev. 196, 62–79. doi:10.1016/j.quascirev.2018.07.029

CrossRef Full Text | Google Scholar

Bliedtner, M., Zech, R., Schäfer, I., and Suchodoletz, H. v. (2020b). A first Holocene leaf wax isotope-based paleoclimate record from the semi-humid to semi-arid south-eastern Caucasian lowlands. J. Quat. Sci. 35, 625–633. doi:10.1002/jqs.3210

CrossRef Full Text | Google Scholar

Blodgett, T. A., Lenters, J. D., and Isacks, B. L. (1997). Constraints on the origin of paleolake expansions in the central Andes. Earth Interact. 1, 1–28. doi:10.1175/1087-3562(1997)001<0001:cotoop>2.0.co;2

CrossRef Full Text | Google Scholar

Bougealt, P. (1991). “Parameterization of land surface processes in numerical weather prediction,” in Land surface evaporation: Measurement and parameterization. Editors T. J. Schmugge, and J. C. Andre (New York: Springer), 55–92.

Google Scholar

Brodowski, S., Rodionov, A., Haumaier, L., Glaser, B., and Amelung, W. (2005). Revised black carbon assessment using benzene polycarboxylic acids. Org. Geochem. 36, 1299–1310. doi:10.1016/j.orggeochem.2005.03.011

CrossRef Full Text | Google Scholar

Brutsaert, W. H. (1982). Evaporation into the atmosphere: Theory, history, applications. Dordrecht: D. Reidel Publishing.

Google Scholar

Budyko, M. I. (1974). Climate and life. New York: International Geophysics Series 18, Academic Press.

Google Scholar

Bukhrashvili, P., Blocher, F., Tskvitinidze, Z., and Davitashvili, S. (2019). Ausgrabungen in Nazarlebi, Kachetien (Georgien) 2017 und 2018. Mittl. Dtsch. Orient-Gesellschaft Berl. 151, 271–294.

Google Scholar

Buylaert, J.-P., Jain, M., Murray, A. S., Thomsen, K. J., Thiel, C., and Sohbati, R. A. (2012). A robust feldspar luminescence dating method for Middle and Late Pleistocene sediments. Boreas 41, 435–451. doi:10.1111/j.1502-3885.2012.00248.x

CrossRef Full Text | Google Scholar

Connor, S. E., Colombaroli, D., Confortini, F., Gobet, E., Ilyashuk, B. P., Ilyashuk, E. A., et al. (2018). Long-term population dynamics: Theory and reality in a peatland ecosystem. J. Ecol. 106, 333–346. doi:10.1111/1365-2745.12865

CrossRef Full Text | Google Scholar

Connor, S. E., and Kvavadze, E. V. (2008). Modelling late Quaternary changes in plant distribution, vegetation and climate using pollen data from Georgia, Caucasus. J. Biogeogr. 36, 529–545. doi:10.1111/j.1365-2699.2008.02019.x

CrossRef Full Text | Google Scholar

Connor, S. E., and Sagona, A., (2007a). “Environment and society in the late prehistory of southern Georgia, Caucasus,” in Les Cultures du Caucase (Vie-IIIe millénaires avant notre ère). Leurs relations avec le Proche-Orient. Editors B. Lyonnet (Paris: CNRS. Éditions Recherche sur les Civilisations), 21–36.

Google Scholar

Connor, S. E., Thomas, I., and Kvavadze, E. (2007b). A 5600-yr history of changing vegetation, sea levels and human impacts from the Black Sea coast of Georgia. Holocene 17, 25–36. doi:10.1177/0959683607073270

CrossRef Full Text | Google Scholar

Croudace, I. W., Rindby, A., and Rothwell, G. (2006). “Itrax: Description and evaluation of a new multi-function X-ray core scanner,” in New techniques in sediment core analysis. Editors R. G. Rothwell (London: Geological Society London, Special Publications) 267, 51–63.

CrossRef Full Text | Google Scholar

de la Rosa, J. M., Knicker, H., López-Capel, E., Manning, D. A. C., González-Perez, J. A., and González-Vila, F. J. (2008). Direct detection of Black Carbon in soils by Py-GC/MS, carbon-13 NMR spectroscopy and thermogravimetric techniques. Soil Sci. Soc. Am. J. 72, 258–267. doi:10.2136/sssaj2007.0031

CrossRef Full Text | Google Scholar

Dong, G., Li, R., Lu, M., Zhang, G., and James, N. (2020). Evolution of human-environmental interactions in China from the late Paleolithic to the Bronze Age. Prog. Phys. Geogr. Earth Environ. 44, 233–250. doi:10.1177/0309133319876802

CrossRef Full Text | Google Scholar

Dotterweich, M. (2008). The history of soil erosion and fluvial deposits in small catchments of central Europe: Deciphering the long-term interaction between humans and the environment - a review. Geomorphology 101, 192–208. doi:10.1016/j.geomorph.2008.05.023

CrossRef Full Text | Google Scholar

Erb-Satullo, N. L., Jachvliani, D., Kalayci, T., Puturidze, M., and Simon, K. (2019). Investigating the spatial organisation of Bronze and iron age fortress complexes in the south Caucasus. Antiquity 93, 412–431. doi:10.15184/aqy.2018.191

CrossRef Full Text | Google Scholar

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., et al. (2007). The shuttle radar topography mission. Rev. Geophys. 45, RG2004–33. doi:10.1029/2005rg000183

CrossRef Full Text | Google Scholar

Fletcher, W., Faust, D., and Zielhofer, C. (2013). Fragile landscape systems. Catena 103, 1–2. doi:10.1016/j.catena.2012.02.006

CrossRef Full Text | Google Scholar

Force, E. R. (2017). Seismic environments of prehistoric settlements in northern Mesopotamia: A review of current knowledge. Bull. Am. Sch. Orient. Res. 378, 55–69. doi:10.5615/bullamerschoorie.378.0055

CrossRef Full Text | Google Scholar

Forte, A., Cowgill, E., Bernardin, T., Kreylos, O., and Hamann, B. (2010). Late Cenozoic deformation of the Kura Fold-thrust belt, southern Greater Caucasus. Geol. Soc. Am. Bull. 122, 465–486. doi:10.1130/b26464.1

CrossRef Full Text | Google Scholar

Freytet, P., and Verrecchia, E. P. (2002). Lacustrine and palustrine carbonate petrography: An overview. J. Paleolimnol. 27, 221–237. doi:10.1023/a:1014263722766

CrossRef Full Text | Google Scholar

Fuchs, M., and Lang, A. (2008). Luminescence dating of hillslope deposits - a review. Geomorphology 109, 17–26. doi:10.1016/j.geomorph.2008.08.025

CrossRef Full Text | Google Scholar

Furtwängler, A. E., Knauß, F., and Motzenbäcker, I. (1998). Archäologische Expedition in Kachetien 1997. Ausgrabungen in Shiraki. Eurasia Antiq. 4, 309–364.

Google Scholar

Furtwängler, A. E., Knauß, F., and Motzenbäcker, I. (1999). Archäologische Expedition in Kachetien 1998. Ausgrabungen in Shiraki. Eurasia Antiq. 4, 233–270.

Google Scholar

Gamkrelidze, I. P. (2003). Geological map of Georgia 1:500,000. Tbilisi: Georgian State Department of Geology and National Oil Company “SAQNAFTOBI”.

Google Scholar

Glaser, B., Haumaier, L., Guggenberger, G., and Zech, W. (1998). Black carbon in soils: The use of benzenecarboxylic acids as specific markers. Org. Geochem. 29, 811–819. doi:10.1016/s0146-6380(98)00194-6

CrossRef Full Text | Google Scholar

Gogichaishvili, L. K. (1984). “Vegetational and climatic history of the Western part of the Kura river basin,” in Palaeoclimates, palaeoenvironments and human communities in the eastern Mediterranean region in later prehistory. Editors e. d. Bintliff, and W. van Zeist (Oxford: BAR International Series), 133, 325–341.

Google Scholar

Hanesch, M., and Scholger, R. (2005). The influence of soil type on the magnetic susceptibility measured throughout soil profiles. Geophys. J. Int. 161, 50–56. doi:10.1111/j.1365-246x.2005.02577.x

CrossRef Full Text | Google Scholar

Harden, C. P., Chin, A., English, M. R., Fu, R., Galvin, K. A., Gerlak, A. K., et al. (2014). Understanding human-landscape interactions in the “Anthropocene”. Environ. Manage. 53, 4–13. doi:10.1007/s00267-013-0082-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrmann, J. T., and Hammer, E. L. (2019). Archaeo-geophysical survey of Bronze and Iron Age fortress landscapes of the south Caucasus. J. Archaeol. Sci. Rep. 24, 663–676. doi:10.1016/j.jasrep.2019.02.019

CrossRef Full Text | Google Scholar

Hijmans, R. J. (2020). Raster: Geographic data analysis and modeling. Available at: https://CRAN.R-project.org/package=raster (Accessed October 11, 2021).

Google Scholar

Holliday, V. T., and Gartner, W. G. (2007). Methods of soil P analysis in archaeology. J. Archaeol. Sci. 34, 301–333. doi:10.1016/j.jas.2006.05.004

CrossRef Full Text | Google Scholar

Ismail-Zadeh, A., Adamia, S., Chabukiani, A., Chelidze, T., Cloethingh, S., Floyd, M., et al. (2020). Geodynamics, seismicity, and seismic hazards of the Caucasus. Earth. Sci. Rev. 207, 103222. doi:10.1016/j.earscirev.2020.103222

CrossRef Full Text | Google Scholar

Jahn, R., Blume, H. P., Asio, V. B., Spaargaren, O., and Schad, P. (2006). Guidelines for soil description. Rome: FAO. fourth edition.

Google Scholar

Joannin, S., Ali, A. A., Ollivier, V., Roiron, P., Peyron, O., Chevaux, S., et al. (2014). Vegetation, fire and climate history of the Lesser Caucasus: A new Holocene record from Zarishat fen (Armenia). J. Quat. Sci. 29, 70–82. doi:10.1002/jqs.2679

CrossRef Full Text | Google Scholar

Jordanova, D., and Jordanova, N. (1999). Magnetic characteristics of different soil types from Bulgaria. Studia Geophys. Geod. 43, 303–318. doi:10.1023/a:1023398728538

CrossRef Full Text | Google Scholar

Kars, R. H., Reimann, T., Ankjaergaard, C., and Wallinga, J. (2014). Bleaching of the post-IR IRSL signal: New insights for feldspar luminescence dating. Boreas 43, 780–791. doi:10.1111/bor.12082

CrossRef Full Text | Google Scholar

Kldiashvili, D. (2021). “Earthquakes in ancient Georgia IV-XVIII centuries,” Manuscript of SRNSFG grant FR 217804 (Tbilisi: Korneli Kekelidze Georgian National Centre of Manuscripts).

Google Scholar

Koinig, K. A., Shotyk, W., Lotter, A. F., Ohlendorf, C., and Sturm, M. (2003). 9000 years of geochemical evolution of lithogenic major and trace elements in the sediment of an alpine lake - the role of climate, vegetation, and landuse history. J. Paleolimnol. 30, 307–320. doi:10.1023/a:1026080712312

CrossRef Full Text | Google Scholar

Kottelat, M., and Freyhof, J. (2007). Handbook of European freshwater fishes. Berlin: Publications Kottelat, Cornol and Freyhof.

Google Scholar

Kreutzer, S., Martin, L., Dubernet, S., and Mercier, N. (2018). The IR-RF alpha-Efficiency of K-feldspar. Radiat. Meas. 120, 148–156. doi:10.1016/j.radmeas.2018.04.019

CrossRef Full Text | Google Scholar

Kunze, R. (2017). Living and working in Late Bronze/Early Iron Age Georgia: The settlements of Udabno in Kakheti (eastern Georgia) and a contribution to metallurgy based on a field survey in the upper Alazani River basin. Stud. Cauc. Archaeol. 3, 54–83.

Google Scholar

Kvavadze, E. V., and Todria, V. T. (1992). Ecologic conditions for the humans of the Bronze and early Iron Age in Udabno Gareji (eastern Georgia) according to paleoclimate data. Tbilisi: Metsniereba. (in Russian).

Google Scholar

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B. (2004). A long-term numerical solution for the insolation quantities of the Earth. Astron. Astrophys. 428, 261–285. doi:10.1051/0004-6361:20041335

CrossRef Full Text | Google Scholar

Lazar, M., Cline, E. H., Nickelsberg, R., Shahack-Gross, R., and Yasur-Landau, A. (2020). Earthquake damage as a catalyst to abandonment of a middle Bronze age settlement: Tel Kabri, Israel. PLoS ONE 15 (9), e0239079. doi:10.1371/journal.pone.0239079

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewin, J., and Macklin, M. G. (2003). Preservation potential for Late Quaternary river alluvium. J. Quat. Sci. 18, 107–120. doi:10.1002/jqs.738

CrossRef Full Text | Google Scholar

Liu, X., Lu, R., Jia, F., Li, X., Li, X., Li, M., et al. (2021). The strategy and environmental significance of Neolithic subsistence in the Mu Us Desert, China. Quat. Int. 574, 68–77. doi:10.1016/j.quaint.2020.12.006

CrossRef Full Text | Google Scholar

Maisuradze, B., and Mindiashvili, G. (1999). Aeroarchaeology of Shiraki. Dziebani 4, 29–36. (in Georgian).

Google Scholar

Manning, S. W., Smith, A. T., Khatchadourian, L., Badalyan, R., Lindsay, I., Green, A., et al. (2018). A new chronological model for the Bronze and Iron Age south Caucasus: Radiocarbon results from project ArAGATS, Armenia. Antiquity 92, 1530–1551. doi:10.15184/aqy.2018.171

CrossRef Full Text | Google Scholar

Maruashvili, L. I. (1971). “Geomorphology of Georgia,” in The relief of the Georgian SSR in the aspects of layers, origin, dynamics and history (Tbilisi: Edition Metsnerieba). (in Russian).

Google Scholar

Matchavariani, L., (2019). The soils of Georgia. Cham: Springer Nature.

Google Scholar

Menges, J., Hovius, N., Andermann, C., Dietze, M., Sowboda, C., Cook, C. L., et al. (2019). Late Holocene landscape collapse of a trans‐Himalayan dryland: Human impact and aridification. Geophys. Res. Lett. 46, 13814–13824. doi:10.1029/2019gl084192

CrossRef Full Text | Google Scholar

Messager, E., Belmecheri, S., Grafenstein, U. von, Nomade, S., Ollivier, V., Voinchet, P., et al. (2013). Late Quaternary record of the vegetation and catchment-related changes from lake Paravani (Javakheti, south Caucasus). Quat. Sci. Rev. 77, 125–140. doi:10.1016/j.quascirev.2013.07.011

CrossRef Full Text | Google Scholar

Moretti, M., and Sabato, L. (2007). Recognition of trigger mechanisms for soft-sediment deformation in the Pleistocene lacustrine deposits of the SantʻArcangelo Basin (Southern Italy): Seismic shock vs. overloading. Sediment. Geol. 196, 31–45. doi:10.1016/j.sedgeo.2006.05.012

CrossRef Full Text | Google Scholar

Murray, A. S., and Wintle, A. G. (2000). Luminescence dating of quartz using an improved single-aliquot regenerative-dose protocol. Radiat. Meas. 32, 57–73. doi:10.1016/s1350-4487(99)00253-x

CrossRef Full Text | Google Scholar

Neff, J. C., Ballantyne, A. P., Farmer, G. L., Mahowald, N. M., Conroy, J. L., Landry, C. C., et al. (2008). Increasing eolian dust deposition in the western United States linked to human activity. Nat. Geosci. 1, 189–195. doi:10.1038/ngeo133

CrossRef Full Text | Google Scholar

New, M., Lister, D., Hulme, M., and Makin, I. (2002). A high-resolution data set of surface climate over global land areas. Clim. Res. 21, 1–25. doi:10.3354/cr021001

CrossRef Full Text | Google Scholar

Nur, A., and Cline, E. H. (2000). Poseidon’s horses: Plate tectonics and earthquake storms in the late Bronze Age Aegean and eastern Mediterranean. J. Archaeol. Sci. 27, 43–63. doi:10.1006/jasc.1999.0431

CrossRef Full Text | Google Scholar

O’Henry, D. O., Cordova, C. E., Portillo, M., Albert, R.-M., deWitt, R., and Emery-Barbier, A. (2017). Blame it on the goats? Desertification in the Near East during the Holocene. Holocene 27, 625–637. doi:10.1177/0959683616670470

CrossRef Full Text | Google Scholar

Olley, J., Caitcheon, G., and Murray, A. (1998). The distribution of apparent dose as determined by optical stimulated luminescence in small aliquots of fluvial quartz: Implications for dating young sediments. Quat. Sci. Rev. 17, 1033–1040. doi:10.1016/s0277-3791(97)00090-5

CrossRef Full Text | Google Scholar

Ollivier, V., Fontugne, M., Lyonnet, B., and Chataigner, C. (2016). Base level changes, river avulsions and Holocene human settlement dynamics in the Caspian Sea area (middle Kura valley, South Caucasus). Quat. Int. 395, 79–94. doi:10.1016/j.quaint.2015.03.017

CrossRef Full Text | Google Scholar

Ollivier, V., Fontugne, M., and Lyonnet, B. (2015). Geomorphic response and ¹⁴C chronology of base-level changes induced by Late Quaternary Caspian Sea mobility (middle Kura Valley, Azerbaijan). Geomorphology 230, 109–124. doi:10.1016/j.geomorph.2014.11.010

CrossRef Full Text | Google Scholar

Pitskhelauri, K. (2018). “Brief report of archaeological works on Didnauri settlement and burial,” in Collection of brief reports from 2017 archaeological expeditions. Editors N. Antidze (Tbilisi: National Agency of Cultural Heritage Protection of Georgia), 24–29.

Google Scholar

Pitskhelauri, K., Elashvili, M., Kirkitadze, G., Navrozashvili, L., Adikashvili, L., Janelidze, Z., et al. (2016). Reconstruction of the paleoenvironment of the Shiraki plain - traces of early state formations in southern Caucasus. Preliminary results. Proc. Georgian Natl. Acad. Sci. 1, 86–100. (in Georgian).

Google Scholar

Prescott, J., and Hutton, J. (1988). Cosmic ray and gamma ray dosimetry for TL and ESR. Int. J. Radiat. Appl. Instrum. Part D. Nucl. Tracks Radiat. Meas. 14, 223–227. doi:10.1016/1359-0189(88)90069-6

CrossRef Full Text | Google Scholar

Rees-Jones, J., and Tite, M. S. (2007). Optical dating results for British archaeological sediments. Archaeometry 39, 177–187. doi:10.1111/j.1475-4754.1997.tb00797.x

CrossRef Full Text | Google Scholar

Reimer, P. J., Austin, W. E. N., Bard, E., Bayliss, A., Blackwell, P. G., Bronk Ramsey, C., et al. (2020). The IntCal20 northern hemisphere radiocarbon age calibration curve (0-55 cal kBP). Radiocarbon 62, 725–757. doi:10.1017/rdc.2020.41

CrossRef Full Text | Google Scholar

Reisser, M., Purves, R. S., Schmidt, M. W. I., and Abiven, S. (2016). Pyrogenic carbon in soils: A literature-based inventory and a global estimation of its content in soil organic carbon and stocks. Front. Earth Sci. 4, 1–14. doi:10.3389/feart.2016.00080

CrossRef Full Text | Google Scholar

Rodionov, A., Amelung, W., Peinemann, N., Haumaier, L., Zhang, X., Kleber, M., et al. (2010). Black carbon in grassland ecosystems of the world. Glob. Biogeochem. Cycles 24, GB3013. doi:10.1029/2009GB003669

CrossRef Full Text | Google Scholar

Rodríguez-Pascua, M. A., Calvo, J. P., de Vicente, G., and Gómez-Gras, D. (2000). Soft-sediment deformation structures interpreted as seismites in lacustrine sediments of the Prebetic Zone, SE Spain, and their potential use as indicators of earthquake magnitudes during the Late Miocene. Sediment. Geol. 135, 117–135. doi:10.1016/s0037-0738(00)00067-1

CrossRef Full Text | Google Scholar

Rosen, A. M., Hart, T. C., Farquhar, J., Schneider, J. S., and Yadmaa, T. (2019). Holocene vegetation cycles, land-use, and human adaptations to desertification in the Gobi Desert of Mongolia. Veg. Hist. Archaeobot. 28, 295–309. doi:10.1007/s00334-018-0710-y

CrossRef Full Text | Google Scholar

Sagona, A. (2018). The archaeology of the Caucasus. Cambridge: Cambridge University Press.

Google Scholar

Salminen, R., (2005). Geochemical Atlas of Europe. Part 1: Background information, methodology and maps. Espoo: Geological Survey of Finland.

Google Scholar

Schlichting, E., Blume, H. P., and Stahr, K. (1995). Bodenkundliches praktikum. Berlin: Blackwell Wissenschaftsverlag.

Google Scholar

Schmidt, A., Quigley, M., Fattahi, M., Azizi, G., Maghsoudi, M., and Fazeli, H. (2011). Holocene settlement shifts and palaeoenvironments on the Central Iranian Plateau: Investigating linked systems. Holocene 21, 583–595. doi:10.1177/0959683610385961

CrossRef Full Text | Google Scholar

Smith, A. T. (2005). Prometheus unbound: Southern Caucasia in prehistory. J. World Prehist. 19, 229–279. doi:10.1007/s10963-006-9005-9

CrossRef Full Text | Google Scholar

Steinhof, A., Altenburg, M., and Machts, H. (2017). Sample preparation at the Jena 14C laboratory. Radiocarbon 59, 815–830. doi:10.1017/rdc.2017.50

CrossRef Full Text | Google Scholar

Stoops, G. (2021). Guidelines for analysis and description of soil and regolith thin sections. Hoboken, NJ: John Wiley & Sons. 2nd ed.

Google Scholar

Strouhalová, B., Ertlen, D., Šefrna, L., Novák, T. J., Virágh, K., and Schwartz, D. (2019). Assessing the vegetation history of European chernozems through qualitative near infrared spectroscopy. Quaternaire 30, 227–241. doi:10.4000/quaternaire.12101

CrossRef Full Text | Google Scholar

Suchodoletz, H. v., Gärtner, A., Zielhofer, C., and Faust, D. (2018). Eemian and post-Eemian fluvial dynamics in the Lesser Caucasus. Quat. Sci. Rev. 191, 189–203. doi:10.1016/j.quascirev.2018.05.012

CrossRef Full Text | Google Scholar

Suchodoletz, H. v., Menz, M., Kühn, P., Sukhishvili, L., and Faust, D. (2015). Fluvial sediments of the Algeti River in southeastern Georgia - an archive of Late Quaternary landscape activity and stability in the Transcaucasian region. Catena 130, 95–107. doi:10.1016/j.catena.2014.06.019

CrossRef Full Text | Google Scholar

Suchodoletz, H. v., Richter, C., Walther, F., Bliedtner, M., Eloshvili, M., Losaberidze, L., et al. (2020). Land snail assemblages in Holocene floodplain research - an example from the southern Caucasus. E G - Quat. Sci. J. 69, 247–260. doi:10.5194/egqsj-69-247-2020

CrossRef Full Text | Google Scholar

Suchodoletz, H. von, Oberhänsli, H., Faust, D., Fuchs, M., Blanchet, C., Goldhammer, T., et al. (2010). The evolution of Saharan dust input on Lanzarote (Canary Islands)—influenced by human activity in the northwest Sahara during the early Holocene? Holocene 20, 169–179. doi:10.1177/0959683609350385

CrossRef Full Text | Google Scholar

Vakhushti Bagrationi Institute of Geography Tbilisi (1997). Atlas of Georgia by Vakhushti Bagrationi. Tbilisi: Tbilisi State University.

Google Scholar

Varazanashvili, O., Tsereteli, N., and Tsereteli, E. (2011). Historical earthquakes in Georgia (up to 1900): Source analysis and catalogue compilation. Tbilisi: Institute of Geophysics I. Javakhishvili Tbilisi State University.

Google Scholar

Varazashvili, V., and Pitskhelauri, K. (2011). Results of archaeological research of Ilia State University at the Iori Upland. Khornabuji: Ilia State University Publishing, 42–105. (in Georgian).

Google Scholar

Vepraskas, M. J., Lindbo, D. L., and Stolt, M. H. (2018). “Redoximorphic features,” in Interpretation of micromorphological features of soils and regolithsEditors e. d. G. Stoops, V. Marcelino, and F. Mees (Amsterdam: Elsevier), 425–445.

CrossRef Full Text | Google Scholar

Virmani, S. M., Sahrawat, K. L., and Burford, J. R. (1982). Physical and chemical properties of vertisols and their management. New Delhi, India: Twelfth International Congress of Soil Science. 8-16 February 1982.

Google Scholar

Wallinga, J., Murray, A. S., Duller, G. A. T., and Törnqvist, T. E. (2001). Testing optically stimulated luminescence dating of sand-sized quartz and feldspar from fluvial deposits. Earth Planet. Sci. Lett. 193, 617–630. doi:10.1016/s0012-821x(01)00526-x

CrossRef Full Text | Google Scholar

Wang, Y., Amundson, R., and Trumbore, S. (1996). Radiocarbon dating of soil organic matter. Quat. Res. 45, 282–288. doi:10.1006/qres.1996.0029

CrossRef Full Text | Google Scholar

White, H. E. (2019). Theses and dissertations (comprehensive). 2176. Available at: https://scholars.wlu.ca/etd/2176 (Accessed March 21, 2022).

Google Scholar

Wilson, A. M., and Jetz, W. (2016). Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions. PLoS Biol. 14 (3), e1002415. doi:10.1371/journal.pbio.1002415

PubMed Abstract | CrossRef Full Text | Google Scholar

Winter, T. C. (1999). Relation of streams, lakes, and wetlands to groundwater flow systems. Hydrogeol. J. 7, 28–45. doi:10.1007/s100400050178

CrossRef Full Text | Google Scholar

Zielhofer, C., Suchodoletz, H. v., Fletcher, W. J., Schneider, B., Dietze, E., Schlegel, M., et al. (2017). Millennial-scale fluctuations in Saharan dust supply across the decline of the African humid period. Quat. Sci. Rev. 171, 119–135. doi:10.1016/j.quascirev.2017.07.010

CrossRef Full Text | Google Scholar