Active deformation and Plio-Pleistocene fluvial reorganization of the western Kura fold–thrust belt, Georgia: implications for the evolution of the Greater Caucasus Mountains

Abstract

Since Plio-Pleistocene time, southward migration of shortening in the eastern part of the Greater Caucasus into the Kura foreland basin has progressively formed the Kura fold–thrust belt and Alazani piggyback basin, which separates the Kura fold–thrust belt from the Greater Caucasus. Previous work argued for an eastward propagation of the Kura fold–thrust belt, but this hypothesis was based on coarse geological maps and speculative ages for units within the Kura fold–thrust belt. Here we investigate the initiation of deformation within the Gombori range in the western Kura fold–thrust belt and evaluate this eastward propagation hypothesis. Sediments exposed in the Gombori range have a Greater Caucasus source, despite the modern drainage network in the NE Gombori range, which is dominated by NE-flowing rivers. Palaeocurrent analyses of the oldest and youngest syntectonic units indicate a switch happened between ~2.7 Ma and 1 Ma from dominantly SW-directed flow to palaeocurrents more similar to the modern drainage network. A single successful ²⁶Al–¹⁰Be burial date indicates the youngest syntectonic sediments are 1.0 ± 1.0 Ma, which, while not a precise age, is consistent with original mapping suggesting these sediments are of Akchagylian–Apsheronian (2.7–0.88 Ma) age. These results, along with recent updated dating of thrust initiation in the eastern Kura fold–thrust belt, suggest that deformation within the Kura fold–thrust belt initiated synchronously or nearly synchronously along-strike. We additionally use topographic analyses to show that the Gombori range continues to be a zone of active deformation.

1. Introduction and motivation

The Caucasus system represents the northern margin of the Arabia–Eurasia collision zone, and from north to south includes the following tectonic units: East European Craton, Scythian Platform, the Greater Caucasus (GC), the Rioni (southwest), Kartli (central) and Kura (southeast) foreland basins, and the Lesser Caucasus Mountains (Cowgill et al. 2016). The tectonic boundary between the Arabian and Eurasian plates in the Caucasus region is a complex zone of compressional tectonics represented dominantly by thrust and reverse faulting (Onur et al. 2019). The Kura fold–thrust belt (KFTB) is located between the GC and Lesser Caucasus Mountains and represents a major structural system within this region, accommodating shortening between these two orogenic belts (e.g. Forte et al. 2010, 2013). Closure of the GC back-arc basin in late Miocene time and the transition from subduction to collision in Pliocene time resulted in a fast exhumation phase of the GC (Avdeev & Niemi, 2011; Vincent et al. 2020); however, the exact timing of collision along-strike between the northern and southern margins of the GC relict back-arc basin remains controversial (e.g. Cowgill et al. 2016; Vincent et al. 2016). Since Plio-Pleistocene time, much of the shortening in the eastern half of the GC has propagated southwards, into the Kura foreland basin, and formed the KFTB (Fig. 1). Since initiation of deformation within the KFTB, it has accommodated approximately half of the total Arabia–Eurasia convergence at the longitude of the eastern GC (~48° E) (Forte et al. 2013). Geodetic measurements indicate that there is an along-strike, eastward increasing velocity gradient between the Greater and Lesser Caucasus, with an approximate convergence gradient that increases from ~3 mm yr⁻¹ to upwards of 10 mm yr⁻¹ along the length of the KFTB (Reilinger et al. 2006; Forte et al. 2014). By analysing large, twentieth century earthquakes in eastern Turkey and the Caucasus along with expected Arabia–Eurasia motion, Jackson & McKenzie (1988) and Jackson (1992) hypothesized that the Caucasus must be deforming mostly aseismically, either by creep on faults or by folding. It might be expected that shortening, especially by folding, of thick, possibly overpressurized, sediments, should occur without generating major earthquakes, even if folding were to occur above buried (blind) thrust or reverse faults (Jackson, 1992). Nevertheless, from the eastern domain of the KFTB in Azerbaijan, there are strong indications that the KFTB is actively deforming (Forte et al. 2010, 2013; Mosar et al. 2010) and thus the potential seismic hazard within the fold–thrust belt may be underestimated. There are several Mw 5–5.4 earthquake events within the KFTB area in the Complete Catalogue of Instrumental Seismicity for Georgia (period of 1900–2017) (Godoladze et al. in prep.). The earthquake data indicate a S-dipping low-angle thrust under the Gombori segment of the KFTB, which is consistent with geological observations throughout the KFTB (Forte et al. 2010, 2013; Adamia et al. 2010, 2011). The strike of the fault plane of a M 5.4 event (27 November 1997) was approximately E–W (Tan & Taymaz, 2006), also consistent with the structural geometries within the KFTB (Fig. 2). However, detailed palaeoseismic studies are absent in the region, leaving significant uncertainties with regards to the seismic hazard.

Fig. 1. Location and topography of the KFTB.

Fig. 2. Earthquake events of the KFTB from the Complete Catalogue of Instrumental Seismicity for Georgia (Onur et al. 2019); fault plane solution by Tan & Taymaz (2006) indicates a compressional fault mechanism.

Previous work on the KFTB noted that there is more elevated topography (measured with respect to the adjacent basins), cross-strike width and older structures exposed in the western part of the belt. Forte et al. (2010) argued this pattern could be caused by an eastward decrease in total shortening, timing of initiation or a combination thereof. According to an analysis of growth strata in seismic profiles and oil well data from the Kura foreland fold–thrust belt by Alania et al. (2017), the formation of the Kakheti range (located in the western KFTB, here referred to as the Gombori range), took place in Pliocene time. Over 300 km east along-strike, the initiation of deformation, based on the age of transition between pre- and syntectonic strata, within the eastern segment of the belt was originally estimated to be between 1.8 Ma and 1.5 Ma by Forte et al. (2013), though more recent and higher resolution dating of this same eastern KFTB stratigraphy suggests the pre- to syntectonic transition, and thus, deformation may have initiated closer to 2.2–2.0 Ma (e.g. Lazarev et al. 2019). Even with this new older age of initiation for deformation within the eastern KFTB, this is still consistent with the idea first proposed by Forte et al. (2010) that deformation started in the western KFTB and propagated eastwards, but it depends on the exact timing of initiation in the western KFTB, which at present is poorly constrained.

Additional evidence of along-strike variation in structural history is interpretable from the topography and comparisons between the modern drainage network and the palaeo-drainage network of the KFTB as reconstructed from alluvial stratigraphy. Specifically, in the eastern KFTB, south-flowing rivers sourced from the GC still cross the KFTB, but west of where the Alazani river enters the KFTB, no south-flowing river from the GC crosses the KFTB (Fig. 1), an additional observation used by Forte et al. (2010) to argue for potential west-to-east propagation of the KFTB. Based on these broad patterns in the drainage network, Forte et al. (2010) speculated that prior to the development of the western KFTB and during some portion of the deposition of pre- and syntectonic alluvial sediments, now exposed within the western KFTB, some GC-sourced rivers did make it to, or through, the KFTB. Such drainage reorganizations during the progressive growth of fold belts is observed in both other natural examples (e.g. Lawton et al. 1994; Burbank et al. 1996; Delcaillau et al. 1998, 2006; Keller et al. 1999; Delcaillau, 2001; Davis et al. 2005; Bretis et al. 2011) and experiments (e.g. Champel, 2002; Douglass & Schmeeckle, 2007). This study tests the hypotheses that (1) a drainage basin reorganization within the western KFTB occurred as speculated by Forte et al. (2010), but which has not been systematically documented in the region in any prior work, and (2) the KFTB propagated from west to east (Forte et al. 2010). To evaluate the hypotheses, within the Gombori range, we analysed palaeocurrent directions preserved within syntectonic strata, ²⁶Al–¹⁰Be burial dating of a single successful sample and quantitative geomorphological analyses of the Gombori range. Ultimately, we find that there is good evidence for drainage reorganization within the western KFTB and Gombori range, but integration of our new data with all available results suggests we are currently unable to distinguish between an eastward-propagating KFTB and one that initiates nearly synchronously along-strike.

2. Stratigraphic background

The Gombori range, which is the highest relief part of the KFTB and defines the NW edge of the belt, is built by deformed lower and upper Cretaceous, Eocene and Oligocene, Miocene, Plio-Pleistocene and Quaternary sedimentary rocks (Figs 3, 4). Here we focus exclusively on the Plio-Pleistocene sediments of the Gombori range, as this portion of the stratigraphy is the most relevant for establishing the neotectonic history of the western KFTB.

Fig. 3. Stratigraphy of the Gombori range compiled after Buleishvili (1974), Zedginidze et al. (1971), Kereselidze (1950), Sidorenko & Gamkrelidze (1964) and Buachidze et al. (1950). Thicknesses are approximate and likely vary along-strike within the Gombori range.

Fig. 4. Simplified lithology (complied according to Soviet-era maps and ground revisions), sampling sites, palaeocurrent directions and selected catchments. Palaeocurrents measured in the Al₁ facies are dominantly SW-directed, while measurements in Al₃ indicate no dominant flow direction but are generally consistent with the present-day Alazani flow direction.

Previous work has described the Plio-Pleistocene sediments of the Gombori range as a part of the Akchagylian–Apsheronian regional stages, and they are collectively described as the Alazani series. Within the Caspian Sea region and its sub-basins, the Akchagylian regional stage corresponds to the late Pliocene epoch (Jones & Simmons, 1996; Krijgsman et al. 2019). The Akchagylian represents a series of large transgressions, which temporarily re-established marine connections between the Caspian Sea and world ocean (Jones & Simmons, 1996; Forte & Cowgill, 2013; Van Baak et al. 2019). The Akchagylian sediments are broadly considered as having been deposited in a marine environment (Jones & Simmons, 1996), but there are continental facies of the Akchagylian stage within the eastern (Forte et al. 2015a), central and western KFTB as well (Sidorenko & Gamkrelidze, 1964; Chkhikvadze et al. 2000; Alania et al. 2017). The Apsheronian stage, which overlies the Akchagylian, is essentially regressive in character and corresponds to the lower and middle Pleistocene (Jones & Simmons, 1996; Krijgsman et al. 2019). It generally represents shallow marine and continental deposits, but, within the Gombori range, Apsheronian sediments are considered part of the Alazani series, which has previously been interpreted as having been deposited in a terrestrial environment (Sidorenko & Gamkrelidze, 1964). The maximum thickness of the Alazani series in the NE slope of the Gombori range is ~1800 m (Fig. 3) (Sidorenko & Gamkrelidze, 1964) between catchments 7 (Kisiskhevi river) and 12 (Papriskhevi river) (Buachidze et al. 1950; Buleishvili, 1974) and thins to ~1400 m along the SW slope of the Gombori range (Sidorenko & Gamkrelidze, 1964) (Fig. 4).

Three facies, Al₁, Al₂ and Al₃, have been previously defined within the Alazani series. There is an angular unconformity at the base of the Alazani series between it and the older Neogene, Palaeogene and Cretaceous sediments. Angular unconformities are also present between all of the Alazani series facies. Our field measurements show that the Al₁ facies has higher dip angles (50–60°), Al₂ has moderate 20–30° dip angles, and the youngest Al₃ facies has the shallowest dips of 5–15° (Fig. 4), broadly suggestive that these strata are syntectonic, i.e. they are growth strata.

The lower Al₁ is represented by well-consolidated conglomerates and cobbles with 0.2–1.5 m thick lenses of silts and clays. The lowest boundary of the Al₁ facies is marked by a bluish colour conglomerate (Fig. 5). The longest axis of cobbles within this conglomeratic interval averages between 10 and 15 cm in length. Sandstone, black shale, limestone and marl clasts are the dominant rock types of the cobbles and conglomerates within the Al₁ facies. Some of these rock types here (e.g. black shale) are typical of the GC and suggest that these sediments are sourced broadly from the north (Buachidze et al. 1950; Sidorenko & Gamkrelidze, 1964), but detailed provenance analyses of these sediments have yet to be performed. The thickness of the Al₁ facies is ~700 m. The Al₁ layers broadly define the Gombori range as an anticlinorium, with Al₁ layers dominantly N-dipping at ~50–60° along the NE slope and SE-dipping at ~20–45° along the southern slope (Fig. 4). The surface elevation of exposures of the lower boundary of the Al₁ facies within the Gombori range is at 481 m, but as it is challenging to distinguish between the lower Al₁ and upper Al₃ facies in the field, a clear upper limit for the Al₁ facies is not yet estimated, but it may reach up to 1991 m.

Fig. 5. Base of the Al₁ series from catchment 7; view to the NW showing steeply NE-dipping conglomeratic (a) and silt beds (b).

The overlying Al₂ facies is mostly dominated by silt and clay, but small amounts of cobbles and conglomerates are also present. As in the Al₁ facies, the sediments in this facies are also suggestive of a GC source (Jurassic black shale) (Buachidze et al. 1952). The maximum thickness of this facies is ~500 m at catchment 6 and gradually decreases to 50 m to a southeastern direction (Buachidze et al. 1952). In the northern slope of the Gombori range, the Al₂ facies dips to the NE, but at shallower angles with respect to the underlying Al₁ facies, with average dips in the Al₂ facies rocks being ~20–30°. This facies contains a thin layer of volcanic ash (Fig. 6). The surface elevation of the Al₂ facies exposures within the Gombori range varies between 448 and 1569 m.

Fig. 6. NE-dipping volcanic ash layer exposed in catchment 12, facies Al₂.

The upper Al₃ facies is dominantly composed of conglomerates with minor interbeds of silts and clays. According to some reports (e.g. Kereselidze, 1950), another volcanic ash layer of 0.4 m thickness is traceable within the silt layer of catchments 6 and 7, but we did not observe this ash layer in the field. The thickness of this facies is between 150 and 250 m. Layers within the Al₃ facies exposed along the NE edge of the Gombori range dip shallowly to the NE at 5–15°. The surface elevation of exposures of the Al₃ facies within the Gombori range varies between 390 and 1210 m. There are also isolated packages of conglomerate higher in the Gombori range that are unconformable with the underlying, older stratigraphy and are likely exposures of the Alazani series. These exposures may be associated with the Al₃ facies, but could also be associated with the Al₁ facies. We attempted to date one of these isolated packages and thus identify to which facies it belonged (see Section 3.c.2), but we were unfortunately unsuccessful.

3. Methods

Plio-Pleistocene sediments are well exposed on the northern slopes of the Gombori range; moreover, three different facies are well expressed in the outcrops of local catchments. Therefore, the 12 largest area catchments (>10 km²) draining the northern slope were selected as primary study areas (see Figs 4, 7).

Fig. 7. TRMM 3B42 pixel extents (black) and catchments of the study area (red) and the identifying numbers for those catchments referenced in the text.

3.a. Palaeocurrent analyses

Modern rivers draining the NE slopes of the Gombori range flow NE and drain into the Alazani basin, but archival data from geological reports (Buachidze et al. 1950, 1952) suggest that the alluvial sediments of the Alazani series (Al₁ and Al₂) contain rock types typical of the GC, suggesting a southward flow of rivers during the deposition of at least some of the sediments.

Alluvial channels are very sensitive to active tectonics and adjust to vertical deformation or base-level change by channel modification (Merritts et al. 1994). Research on fluvial terraces (such as abandoned floodplains) using gravel or pebble imbrication, is one of the reliable indicators of palaeocurrent in coarse-grained deposits and can shed light on the tectonic evolution of the site (Miao et al. 2008). The direction of imbrication of oblate clasts in a conglomerate can be used to indicate the direction of the flow that deposited the gravel (Nichols, 2009).

Based on the quality of exposure and access to these exposures of Alazani series sediments in the walls of canyons along the main stem rivers of catchments 7 and 11, we selected these two catchments for palaeocurrent analyses. A total of 265 clasts were measured from four sites of the Al₁ and Al₃ facies of both catchments (see Table 1). The analysis was not performed for the Al₂ facies because it is dominated by silt. In this study, we measured the orientation of the clast imbrication with a Brunton compass and performed unfolding and further processing using Stereonet 10 software (Allmendinger et al. 2011). We performed this palaeocurrent analysis to specifically test whether there was evidence of flow reversal and/or drainage reorganization during the deposition of the potentially syntectonic Alazani series sediments and whether any change was diachronous between these two catchments.

Table 1. Von Mises distribution results for the palaeocurrent measurements

3.b. Tectonic geomorphology

Topography reflects the balance between rock uplift, driven by tectonics, and erosional and depositional processes modulated by climate and lithology. With careful consideration of potential climatic and lithological complications, quantitative geomorphological analyses can constrain relative differences in rates of rock uplift, and thus improve our understanding of tectonics (e.g. Kirby & Whipple, 2001, 2012; Wobus et al. 2006; Dibiase et al. 2010; Whittaker, 2012; Whittaker & Boulton, 2012; Rossi et al. 2017; Gallen & Wegmann, 2017). Importantly, in the absence of other data, e.g. dense geodetic networks and/or long-term and complete seismic and palaeoseismic records, tectonic geomorphology can also be useful in highlighting areas of active tectonics and potential seismic hazard (e.g. Kirby et al. 2003).

To evaluate the extent to which tectonic activity within the western end of the KFTB may still be localized in the Gombori range, we selected the 12 largest catchments (14–108 km²) along the northern slope of the Gombori range and calculated several morphometric parameters (including maximum local relief, mean slope and channel steepness index) using the Topographic Analysis Kit (TAK) (Forte & Whipple, 2019), TopoToolbox (Schwanghart & Scherler, 2014), QGIS and a digital elevation model (DEM) acquired through the ALOS AW3D30. The DEM is produced by the Japan Aerospace Exploration Agency (JAXA) and has a horizontal resolution of ~30 m (available from https://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm). The AW3D30 DEM dataset was generated based on the 0.15-arcsec AW3D DEM dataset. Two resampling methods were applied to obtain one pixel value on the AW3D30 from 7 by 7 pixels on the AW3D. The first one used the averaging method (Ave), which is simply calculated as an average value from ~49 pixels except for masked-out values. Another is the medium method (Med), which selects a medium height value, i.e. 25th height, from 49 pixels. If it shows a masked value, the same value is kept in the AW3D30. Both the Ave and Med datasets are contained in an individual AW3D30 dataset, which can be downloaded free of charge. The AW3D30 Ave DEM has a vertical accuracy of 5 m (RMSE) (Tadono et al. 2016) using the EGM96 vertical reference frame (JAXA, 2017). In this study, the average dataset is used.

We attempted to limit our analyses to areas that were bedrock streams, as many of the metrics were designed for application to bedrock rivers. Thus, we avoided the lower portions of catchments, as these portions of the rivers are likely more alluvial in character and, additionally, are in zones subject to intense agricultural activities and other human modifications.

Topography in actively deforming regions is not only influenced by rates of uplift; thus, care must be taken to ensure that tectonic interpretations are not unduly influenced by compounding factors, such as spatially variable precipitation (e.g. Kirby & Whipple, 2012). In regions with strong orographic forcing of precipitation, the influence of variations in discharge on the relationship between channel profile form and erosion rate can be diagnosed with careful analysis (e.g. Bookhagen & Strecker, 2012).

To check how variable precipitation is between catchments, we used satellite data from the Tropical Rainfall Measurement Mission (TRMM) 3B42 V7 collected from 1998 to 2017. The TRMM dataset contains daily rainfall information recorded in 30 km size pixels. TRMM-derived rainfall data is well tested in tectonic geomorphological studies in the Caucasus (Forte et al. 2016), Andes (Bookhagen & Strecker, 2008) and Himalayas (Bookhagen & Burbank, 2006). Figure 7 shows that all 12 catchments are covered by five TRMM pixels.

Similar to climatic influences on topography, contrasts in the relative erodibility of different lithologies can produce patterns in topography that may be confused with tectonic signals (e.g. Mitchell & Yanites, 2019). Lithological contacts and catchment-dominant rock types were identified to include in tectonic geomorphological analyses. To check the correlations between lithological units (including rock properties) and topographic indices, we compiled several Soviet-era geological maps with new field observations and mapping. For each catchment, we calculated the dominant rock types (according to surface area) to correlate this data to other tectonic geomorphological proxies.

For our quantitative topographic analyses, we calculated the normalized channel steepness index (k _sn), catchment-averaged normalized channel steepness index, catchment relief, catchment-averaged hillslope gradient (S _avg), catchment-averaged local relief calculated using a 1 km radius circle and drainage area for all selected catchments.

3.b.1. Channel steepness index

The normalized channel steepness index is an important topographic metric (e.g. Dibiase et al. 2010). Despite incomplete understanding of the varied processes contributing to fluvial erosion, the stream profile method has proven an invaluable qualitative tool for neotectonic investigations. When controlling for differences in precipitation and lithology, empirical observations and simple models of fluvial erosion suggest a positive correlation between channel gradient and rock uplift rate (e.g. Wobus et al. 2006), and thus the normalized channel steepness index can be used in active ranges to illustrate relative rates of rock uplift (e.g. Dibiase et al. 2010), which in our case may illustrate the neotectonic activity of the Gombori range.

Typical river longitudinal profiles, for both bedrock and alluvial rivers, are concave and can be described by an empirical power law relationship between slope and area:

$${k_{sn}} = {\rm{S}}{{\rm{A}}^{\rm{\theta }}}$$

where k _sn is the normalized channel steepness index, S is slope, A is the upstream contributing drainage area and θ is the channel concavity index (Flint, 1974). Numerous studies indicate that most channels have uniform concavity regardless of the (spatially constant) uplift rate (Snyder et al. 2000; Whipple, 2004), because the concavity index (θ) is relatively insensitive to differences in rock uplift rate, climate or substrate lithology at steady-state (provided such differences are uniform along the length of the channel), while the steepness index (k _sn) varies with these factors; therefore, the steepness index is a useful metric for tectonic geomorphological studies (Kirby & Whipple, 2012).

To normalize channel steepness indices, we used a reference concavity (θ_ref) of 0.5, because, in practice, it is found that values of θ_ref between 0.4 and 0.5 work well for most mountain rivers (Kirby & Whipple, 2012). Normalization of the channel steepness index allows for the comparison of river profile morphology between streams and watersheds of different drainage areas.

3.b.2. Local relief

Local relief is the difference between minimum and maximum elevations within a specified distance and is strongly correlated with erosion rate (Ahnert, 1970; Montgomery & Brandon, 2002; Kirby et al. 2003; Dibiase et al. 2010), which is well correlated to rock uplift rate (e.g. Kirby & Whipple, 2001; Lague, 2014). We used a 1 km radius circle to generate local relief.

3.c. Cosmogenic nuclide burial age dating

The ages of the Alazani series sediments are particularly important as the age of these syntectonic sediments could help constrain the age of initiation of this portion of the KFTB. Because the Alazani series sediments lack abundant ash horizons and are mostly too coarse grained for magnetostratigraphy or the preservation of microfauna useful for biostratigraphic correlation, we attempted to constrain the age of these sediments through the use of cosmogenic nuclide burial age dating. Terrestrial cosmogenic nuclides (TCNs), such as ¹⁰Be and ²⁶Al, are produced by the interaction of secondary particles, produced in the Earth’s atmosphere during interaction with cosmic rays, with Earth materials (e.g. see review by Gosse & Philips, 2001). The accumulation of TCNs in Earth materials is a function of depth, the duration of exposure, the erosion rate of the surface and the production rate of the isotope in question, which is a function of latitude and elevation. Importantly, production of TCNs goes to zero below an attenuation depth such that virtually all production occurs in the first 1–3 m of the Earth’s surface. Measuring the cosmogenic nuclide abundances in sediment eroded from upland catchments and then deposited in adjacent basins records both a palaeo-erosion rate and a time since burial (e.g. Granger et al. 1997; Granger & Muzikar, 2001; Granger, 2006). As there are two unknowns, i.e. the exposure history prior to burial and the time of burial, it is necessary to measure the concentration of two separate TCNs with different half-lives, which in this study are ¹⁰Be and ²⁶Al. It is assumed that both burial of these sediments and shielding from any further production of TCNs post-burial (i.e. burial below several metres) occurs rapidly and that the sediments in question remain shielded until nearly the time of collection (e.g. see review by Granger et al. 1997).

3.c.1. Sample collection and preparation

We targeted three samples for cosmogenic nuclide burial age dating to estimate the depositional age of the youngest facies (Al₃) and of an unidentified isolated package that could belong to Al₃ or Al₁ (sample ID: GOMSS01). Owing to the poor concentration of quartz, we sampled the Al₃ facies at two different locations (sample IDs: GOMSS02 and GOMSS03) (see Table 2 and Fig. 4).

Table 2. Burial age sampling site information

GOMSS01: The site is within either facies Al₁ or Al₃ (i.e. we attempted to date this package to help constrain the age of sediments unconformably exposed at the crest of the Gombori range), within catchment 7, 1.5 km southeast from the highest peak of the range called Tsivi (1991 m). The sample was collected from the bottom of a 1.0 m deep pit that we dug. The upper 0.2 m of this pit was soil and the rest was conglomerate. The sampling site was stable and undisturbed.

GOMSS02: The sample was taken from the lowest edge of an outcrop exposed along the Turdo river (catchment 6) from Al₃. The sampling spot was already eroded ~1.5–2.0 m back by the active channel, and we excavated an additional 0.4 m back into the vertical face. Sample GOMSS02 was taken from 0.5 m above the floodplain and 14 m below the surface of the canyon wall (horizontal depth dug: 0.4 m; dip: 5°; dip direction: 2°). The sampling site was stable and undisturbed.

GOMSS03: The sample was taken 500 m upstream from GOMSS02 within the same outcrop belt along the Turdo river (catchment 6), from 1.88 m above the floodplain from the Al₃ facies. The outcrop was eroded back ~0.9 m by the active channel, and we excavated an additional 0.4 m into the vertical face. The sampling location was 66 m below the surface of the canyon wall (dip: 10°; dip direction: 5°).

Of the three samples collected, two (GOMSS01 and GOMSS03) yielded sufficient quartz for dating. The isolation and purification of quartz, dissolution, column chemistry and precipitation of Be and Al oxides was performed in the cosmogenic isotope laboratory at Arizona State University. Isolation of quartz in these samples required modification of standard methods (e.g. Kohl & Nishiizumi, 1992), because of significant fractions of fine-grained, quartz-rich lithic material that dissolved at similar rates in HF and HNO₃ leaches as the quartz being targeted for analysis. Thus, after the initial step of cleaning in Aqua Regia, instead of proceeding directly to leaching in HF and HNO₃, we first used the hot phosphoric acid (HPA) method (Mifsud et al. 2013) to remove feldspars and break up these lithic clasts. After HPA, samples were leached with HF and HNO₃ as in the standard procedure (e.g. Kohl & Nishiizumi, 1992). After cleaning and during dissolution, samples were spiked with commercial ¹⁰Be carrier. We then extracted ¹⁰Be and ²⁶Al through column chromatography (Ditchburn & Whitehead, 1994), and nuclide ratios were measured via accelerator mass spectrometry at the Purdue Rare Isotope (PRIME) Laboratory at Purdue University. We measured native Al concentrations for the two samples using a Thermo iCAP6300 inductively coupled plasma optical emission spectrometer (ICP-OES) at Arizona State University’s Goldwater Environmental Laboratory.

3.c.2. Modelling burial age dates

For the two burial age samples that yielded sufficient quartz, we used CosmoCalc v3.0, a Microsoft Excel add-in, for cosmogenic nuclide calculations (Vermeesch, 2007). We used the default settings for calibration sites for ¹⁰Be and ²⁶Al production and production mechanisms within CosmoCalc v3.0 and report the results of using the Burial-Exposure function within CosmoCalc’s Age/Erosion rate calculator, though we also tested the Burial-Erosion function, which produced similar estimations of burial age. CosmoCalc provides two different numerical methods for fitting burial dates, the Metropolis and Newton’s method. We tested both methods and found that the Metropolis method, which is more complicated, produced variable burial ages, i.e. running the calculation multiple times yielded different results, but that given the magnitude of the uncertainty, this variability in burial ages was small and the error ranges for the simpler Newton’s method were extremely large. Importantly, for most runs, the reported burial ages using the Newton’s and Metropolis methods were similar, and the error ranges reported from the Metropolis method were largely consistent between runs. We elected to report values from the Metropolis method, as these likely reflect a more reasonable range of uncertainties on the burial ages (e.g. Vermeesch, 2007). To account for the variability in reported burial age from multiple runs of the Metropolis method, we report the average of the result of ten runs.

To determine a burial age, a production scaling factor must be assumed for the area that originally contributed the sediment that was eventually eroded, transported, deposited and then buried. While the exact parameters included in different scaling schemes vary, in general, latitude and elevation will be the most important factors controlling the production rate (e.g. Gosse & Philips, 2001). Because the source of sediment for the Alazani series sediments is not well constrained, we tested four different scaling schemes assuming different source areas. Specifically, we tested a ‘local’ sourcing using a latitude and mean elevation appropriate for a representative catchment in the northern Gombori range, and then three different sources from the GC with representative latitudes and mean elevations for a river draining the higher portions of the central GC (e.g. the modern Aragvi river), one draining an intermediate set of elevations (e.g. the modern Iori river) and one draining lower elevations coming directly from the small catchments that drain into the Alazani valley from the central and eastern GC (for detailed parameters see online Supplementary Material Table S2). For the calculation of scaling factors, we used the CosmoCalc implementation of the Desilets et al. (2006) scheme. The calculated burial ages are reported in online Supplementary Material Table S1 for GOMSS03; calculations were not performed for GOMSS01 as an age is not interpretable for this sample as it plots in the region above the constant exposure line, outside the range of physically possible results.

4. Results

4.a. Palaeocurrent analyses

Palaeocurrent analyses of outcrops of Al₁ in two catchments indicate that Al₁ sediments were deposited by a river flowing in a SW direction through the modern Gombori range, counter to the modern drainage direction and consistent with rivers sourced from the GC. For the younger, stratigraphically higher Al₃, palaeocurrents no longer indicate a single, dominant flow direction but are generally consistent with the present-day Alazani eastward flow direction (see Fig. 4 and Table 1). At the same time, results from catchment 7 indicate that during the deposition of Al₃, north-flowing rivers were present.

4.b. Tectonic geomorphology

Quantitative tectonic geomorphological analyses show higher channel steepness indices (Fig. 9a, b) from the western catchments. Maximum local relief is also higher in the western catchments (Figs 8, 9c), which is consistent with the observation in many landscapes that mean normalized channel steepness and local relief are often linearly related (Dibiase et al. 2010). A simple interpretation of these two indices would suggest that the western part of the Gombori range is uplifting faster than its eastern segment, as channel steepness and local relief are often positively correlated with rock uplift (e.g. Kirby & Whipple, 2012).

Fig. 8. (Top) Topography and local relief maps of catchments. (Bottom) Catchment-averaged and stream K _sn values. See text for details of these calculations.

Fig. 9. (a) Swath profile of topography, K _sn values, (b, c, d) along-swath geomorphological indices and (e) rainfall data. Each of the points corresponds to a catchment labelled above. Standard errors are represented by bars and labels.

As noted above, tectonic geomorphological proxies could be influenced by rainfall and lithology. Indeed there are strong correlations between rainfall and each of catchment-mean elevation, local relief and mean k _sn (r² = 0.84, r² = 0.68 and r² = 0.89, respectively) (see correlation matrix in online Supplementary Material Table S3). This likely reflects expected orographic enhancement of rainfall such that areas of high relief, channel steepness and mean elevation driven by high rock uplift rates are associated with high rates of precipitation. Importantly, a climatic control on topography would imply reduced relief and channel steepness in areas of enhanced precipitation, which we do not observe. Thus, interpreting topography as reflecting rock uplift-rate patterns alone is a conservative assumption. We also evaluated whether lithology significantly influenced our tectonic geomorphological indexes, but correlations between dominant rock types and geomorphological proxies are low, as the correlation coefficients between mean k _sn and K (Cretaceous rocks) and Ak–Ap (Akchagylian–Apsheronian sediments) are 0.42 and −0.46. The higher slopes of conglomerate-dominated catchments could be explained by the tendency of the conglomerate deposits to be exposed as cliffs.

4.c. Burial age dates

Online Supplementary Material Table S1 summarizes the analytical results. Unfortunately, one of our samples, GOMSS01, yielded a ²⁶Al/¹⁰Be ratio that even within the uncertainty bounds plots entirely above the constant exposure line of the standard erosion island plot, in the so-called ‘forbidden zone’ (Fig. 10). Data that plots in this region is physically impossible as the ²⁶Al/¹⁰Be ratio cannot exceed the ratio of the production rates of the two isotopes because ²⁶Al decays faster than ¹⁰Be. This suggests that there was a methodological error during processing; thus, a burial age is not interpretable from this sample. The other sample, GOMSS03 from Al₃, did yield an interpretable age, but because of relatively high concentrations of native Al and low concentrations of ²⁶Al, the analytical precision of this measurement is quite low, yielding a burial age of ~1.0 Ma, with lower and upper bounds of 0.005 Ma and 2.5 Ma, respectively (for complete results see online Supplementary Material Table S1). While imprecise, given that there are no published geochronological ages for the Alazani series, or more broadly for any of the sediments in this region of the KFTB, this age is still meaningful as it confirms that these sediments are most likely Apsheronian in age. Because of the relatively low ¹⁰Be concentration, and thus the relatively high implied palaeo-erosion rates, the uncertainty in the source area for the sediment and associated uncertainty in applicable production scheme does not significantly influence the interpreted age for sample GOMSS03, but does have implications for the implied palaeo-erosion rate (Fig. 10). The minimum and maximum scaling for sample GOMSS03 would imply palaeo-erosion rates within the source area of between ~20 cm ka⁻¹ and ~35 cm ka⁻¹ (or 0.2–0.35 mm a⁻¹), respectively.

Fig. 10. Erosion island plot for Gombori range samples. Variability in production rate scaling for the two samples, GOMSS01 and GOMSS03, are reflected in the pairs of points. Sample GOMSS01 plots in the forbidden zone and is thus uninterpretable. Sample GOMSS03 has mean ages of ~1 Ma regardless of the exact scaling relationships used. The relatively high uncertainties on the ages reflect high native Al concentrations. Burial isochrons are reported in Ma and bounds for estimated palaeo-erosion rates in cm ka⁻¹. Plots produced using CosmoCalc (Vermeesch, 2007).

5. Discussion

5.a. Initiation and development of the western Kura fold–thrust belt

The results of our palaeocurrent analyses suggest that a major drainage reorganization and flow reversal of rivers within the western KFTB started during or after the deposition of the Al₁ facies within the Alazani series and finished during or after deposition of the Al₃ facies. As a result, portions of rivers flowing from the GC in a SW direction changed their course to SE, flowing along the piggyback Alazani basin, and new rivers started flowing from the newly emerged Gombori range in a NE direction, into the Alazani basin. We attribute this drainage reorganization to initiation, or intensification of uplift of the western KFTB at this longitude during the time period spanning the deposition of the Alazani series (Fig. 11).

Fig. 11. Fluvial system evolution diagram for the western KFTB. (a) During the deposition of Alazani Suite 1 (Al₁), rivers draining from the Greater Caucasus were still able to flow directly south across what is now the KFTB. (b) Alazani Suite 2 (Al₂) represents deposition in a lacustrine setting, which could relate to damming of rivers by growth of the KFTB, or could be related to broader, basin-wide changes in base level. (c) By the time of deposition of Alazani Suite 3 (Al₃), the river network in the northwestern KFTB had developed into something similar to the modern situation, with rivers draining northward out of the Gombori range and with a well-defined axial drainage occupying the Alazani basin.

The timing of drainage reorganization in the western KFTB is an important constraint on the structural and topographic evolution of this portion of the KFTB and thus helps constrain the along-strike evolution of the KFTB overall. The sediments of the Alazani series were previously mapped as being a part of the Akchagylian–Apsheronian stages. The reported age for the base of the Akchagylian is variable between publications and regions (e.g. Krijgsman et al. 2019), but it has been constrained to be ~2.7 Ma near the Azerbaijan Caspian Sea coast based on ⁴⁰Ar–³⁹Ar dating of an ash horizon (Van Baak et al. 2019). It is suggested that the base of the Akchagylian may be time transgressive, and in a section ~150 km to the east of the Gombori range it has been constrained to be ~2.5 Ma based on the maximum depositional age from detrital zircons in the strata below the Akchagylian (Forte et al. 2015a). The boundary between the Akchagylian and Apsheronian stages is similarly variable, but in the vicinity of the KFTB, the Apsheronian has been dated to extend from 2.2 Ma to 0.88 Ma (e.g. Krijgsman et al. 2019).

According to this information, we can make an attempt to estimate the ages and reconstruct the depositional environment and tectonic context of the Alazani series. Deposition of Al₁ sediments started not earlier than ~2.7–2.5 Ma years ago by the streams flowing from the GC to the southwest through the location of the modern Gombori range area into the Kura basin. We hypothesize that during deposition of Al₁, uplift of the Gombori range initiated and potentially dammed the formerly south-flowing rivers, which could explain the finer, more lacustrine sediments in Al₂, though given the uncertainty in the exact age of the Al₂ facies and the broad context of the Akchagylian stage as a transgressive event, it is not possible to rule out a more regional explanation for the lacustrine character of the Al₂ facies. Regardless, by the time of deposition of Al₃, sufficient deformation and uplift had accrued in the Gombori range to effect a significant drainage reorganization and the development of (1) a set of north-flowing rivers on the Gombori range (see Fig. 11) and (2) an axial valley, i.e. the Alazani valley, between the Gombori range and the GC. We interpret the lack of a dominant palaeocurrent direction in these Al₃ facies sediments to reflect possible deposition within this axial valley, which today is dominated by a set of meandering fluvial systems. This would imply that the northeastern extent of the Gombori range has expanded since the deposition of Al₃, i.e. at the time of deposition the palaeocurrent sites were not within the deformed part of the Gombori range, but have subsequently been incorporated into the range. Comparison between the interpreted palaeo-drainage network and the modern drainage network suggests that uplift in the Gombori range was sufficiently rapid such that river(s) could not maintain antecedent gorges (Humphrey & Konrad, 2000) like they currently do in the eastern KFTB (see Forte et al. 2010). Considering the burial age sample from Al₃ in the context of the palaeocurrent results implies that, during deposition of this facies, the palaeo-erosion rate in the source region for these sediments in the Gombori range was 0.35–0.2 mm yr⁻¹ (e.g. Fig. 10). This provides at least some constraint on the relative magnitude of uplift rates within the Gombori range, though the extent to which this palaeo-erosion rate is relatable to an uplift rate of the Gombori range depends on whether the assumption of steady-state was valid at that time, i.e. did the average erosion rate equal rock uplift rate, which is challenging to assess retrospectively.

The lack of precise age control for the Alazani series sediments, and that our one successful burial age date only provides a constraint for the time by which a drainage reorganization was in progress, introduces uncertainty in terms of when deformation initiated in the western KFTB. However, if we assume that (1) the age of the base of the Al₁ strata is between 2.7 Ma and 2.5 Ma (the maximum permissible age of the Akchagylian stage in this region), and (2) Al₁ reflects deposition before significant development of the western KFTB and that the age of the Al₃ strata is ~1 Ma (from our burial age date of sample GOMSS03), and (3) deposition of Al₃ reflects a time by which the drainage reorganization had made significant changes in flow direction, this brackets the initiation age of the western KFTB to between 2.7 Ma and 1 Ma. Comparison of this range of possible initiation ages with those observed in the far eastern end of the KFTB, which based on new age constraints (e.g. Lazarev et al. 2019) likely initiated at ~2.2–2.0 Ma, suggests that if there was eastward along-strike propagation of the KFTB as suggested by Forte et al. (2010), it took no more than 0.5–1 Ma. Given the lingering uncertainty in the initiation age of the western KFTB and the newly revised, older age of initiation in the eastern KFTB, it is equally viable that there was no significant propagation along-strike. This uncertainty highlights the need for additional work to establish the ages of the Alazani series stratigraphy in the western KFTB and identify additional areas where the timing of initiation of the KFTB can be assessed along-strike.

5.b. Implications for regional tectonics

Coarse spatial resolution GPS-derived crustal motion velocity data suggest an eastward horizontal velocity increase along-strike within the KFTB (see Reilinger et al. 2006). However, our tectonic geomorphological analyses suggest that the rates of uplift along-strike within the Gombori range are not well correlated with GPS horizontal velocities (with respect to Eurasia). In detail, our results indicate that the western Gombori range may be experiencing more rapid uplift, leading to its generally higher elevation, normalized channel steepness and local relief. If our assumption of almost simultaneous deformation initiation along-strike within the KFTB is correct, there could be several explanations for this apparent disconnect between an eastward increase in GPS velocity with an eastward decrease in local relief within the Gombori range: (1) the along-strike decrease in relief reflects structural complexity, with larger portions of the total convergence being taken up by additional structures to the southeast of the Gombori range; (2) an along-strike change in the ratio of shortening accommodated either currently or through time between the KFTB and the interior of the GC; (3) an along-strike change in structural geometry between steeper to shallower dipping structures from west to east within the KFTB that would result in an eastward decrease in the relationship between incremental total shortening and vertical rock uplift; (4) a first-order control from lithology such that once there is sufficient exhumation to expose older, more-resistant units in the core of folds, this led to an increase in relief compared to adjacent areas, which expose younger, less-resistant units, even if those areas are experiencing greater rates of rock uplift; (5) the modern GPS velocity field is not representative of the long-term, i.e. several million-year, rate of convergence in the region, a suggestion that has been made more broadly for the GC as a whole (Forte et al. 2016); or (6) the Gombori range itself reflects an eastward-propagating set of structures.

At present, we do not have the data to uniquely select among these hypotheses. Option 1 would be consistent with coarse resolution syntheses of structures and estimation of the activity of those structures presented in Forte et al. (2010), but without quantitative assessments of the amounts of total shortening accommodated by structures southeast of the Gombori range (or in the Gombori range itself), this is hard to validate. Similarly, option 2 would be consistent with an eastward along-strike decrease in range front sinuosity for the frontal GC, used as a proxy for time since the GC range front fault was active at the surface, as noted by Forte et al. (2010), but generally not consistent with other observations within the eastern GC of no clear differences along-strike in terms of the tectonic geomorphology of this portion of the range (e.g. Forte et al. 2014, 2015b). Option 3 is not broadly consistent with the observed bedding orientations within the Gombori range as, at least within the Alazani series, there does not appear to be any clear change in the orientation of units along-strike, e.g. Al₁ facies sediments uniformly dip 50–60° along the exposed portion of the Alazani series. For option 4, our analyses of the topography did not indicate that lithology exerts a strong control, but importantly our analyses did not extend beyond the Gombori range. Fully evaluating options 5 or 6 requires detailed estimations of total shortening and timing of initiation along-strike within the KFTB and the Gombori range; however, it is worth noting that as discussed in the previous section, our results along with an updated chronology for stage boundaries have narrowed the range of time over which the KFTB would need to propagate eastward along-strike, leaving open the possibility that a fundamental disconnect between GPS rates and long-term geological rates is viable. Ultimately, this work further highlights the necessity for detailed estimates of the amounts of total shortening and ages of deformation initiation throughout the KFTB.

Previous work from the eastern (Forte et al. 2013) and central (Alania et al. 2017) KFTB concluded that the Kura foreland is an active fold–thrust belt. Our study revealed that the western portion of this belt has experienced large-scale tectonic movements and drainage reorganization that are still in progress. GPS data from the western neighbouring region showed that Tbilisi and the northern boundary of the Lesser Caucasus is a zone of active convergence (Sokhadze et al. 2018), and the sparse GPS network that spans the area south of the Gombori range and GC indicates a horizontal velocity gradient, i.e. convergence, between the Gombori and the GC (Onur et al. 2019). All these data lead us to assume that the western KFTB is actively deforming, and it should be considered during seismic hazard assessment of the region.

6. Conclusions

Our synthesis of the tectonic geomorphology, absolute age dating of syntectonic Plio-Pleistocene sediments within the Kura fold–thrust belt and palaeocurrent analyses within those same sediments reveal a Plio-Pleistocene drainage reorganization event within the northwestern corner of the southeastern foreland of the GC Mountains, which appears linked to initiation and development of the KFTB. If the timing of this drainage reorganization event, constrained to have occurred between ~2.7 Ma and 1 Ma, is representative of initiation of this western-most segment of the KFTB, then along-strike propagation of the fold–thrust belt along its ~300 km length took no more than ~1 million years and opens the possibility of no significant along-strike diachroneity in fold–thrust belt initiation within the KFTB. However, given the lingering uncertainty in the age estimations for the stratigraphy in the western KFTB, the idea of an eastward-propagating KFTB as originally proposed by Forte et al. (2010) cannot be totally excluded, though it appears less likely.

This result has potential implications for the tectonics of the KFTB and GC, as nearly synchronous initiation of the KFTB along its length would further support disconnects between the modern GPS velocity field and the structural and topographic history of the region (e.g. Forte et al. 2016). It, however, remains unclear how to reconcile nearly synchronous initiation of the KFTB along-strike with apparent geological and geomorphological evidence that the western KFTB has experienced more total exhumation and shortening than the eastern KFTB (e.g. Forte et al. 2010). This further highlights the necessity of more quantitative estimates of total shortening and timing of initiation throughout the KFTB to enable testing of tectonic models and firmly establishing the extent to which there are along-strike gradients in shortening or timing of initiation.

Quantitative tectonic geomorphological analyses of the Gombori range indicate that the western KFTB is still a zone of active deformation, especially its NW segment. This is consistent with recently published preliminary GPS velocity data (Onur et al. 2019), suggestive of an across-strike velocity gradient between the western KFTB and GC Mountains.