Abstract

The Greater Caucasus Mountains are a young (∼5 m.y. old) orogen within the Arabia-Eurasia collision zone that contains the highest peaks in Europe and has an unusual topographic form for a doubly vergent orogen. In the east-central part (45°E–49°E), the range is nearly symmetric in terms of prowedge and retrowedge widths and the drainage divide is much closer to the southern margin of the range (prowedge side) than it is to the northern margin (retrowedge side). Moreover, the divide does not coincide with the topographic crest, but rather the crest is both shifted northward by as much as 40 km and traversed by several large north-flowing rivers. Both the topographic crest and drainage divide appear to coincide with zones of active rock uplift, because they are characterized by bands of high local relief and normalized channel steepness values (>300). This uplift pattern could result from a synchronous initiation of the two uplift zones or propagation of deformation either northward or southward. The two propagating scenarios differ fundamentally in their predictions for the relative ages of topographic features; northward propagation predicts that the topographic crest is younger than the drainage divide, and the southward scenario predicts the converse. Because available geologic and topographic data are consistent with both propagation directions, we use a landscape evolution model to test all three scenarios. Model results indicate that the current topography and drainage network is best explained by a northward propagation of deformation from the south flank into the interior of the east-central Greater Caucasus. Such propagation implies recent out-of-sequence deformation within the Greater Caucasus due to reactivation or development of new structures within the core of the orogen. It remains unclear if such deformation is a transient response to an accretion cycle or stems from a fundamental change in the structural architecture of the orogen.

INTRODUCTION

The location of the main drainage divide in an active orogen represents a fundamental control on its structural and topographic evolution, but is also a dynamic feature, strongly predicated on the imposed tectonic and climatic forcing (Fig. 1; e.g., Willett, 1999; Willett et al., 2001). Divides are important because they often separate areas with different climatic conditions, effective base levels, or structural regimes (e.g., Bonnet, 2009; Hoth et al., 2006; Roe et al., 2003; Willett et al., 1993); this can produce significant differences in the topographic form, structural architecture, sediment flux, or exhumation rates on either side of an orogen (e.g., Bonnet and Crave, 2003; Hoth et al., 2006; Montgomery et al., 2001; Whipple, 2009; Willett, 1999). Drainage divides are rarely static (e.g., Willett et al., 2014); within an active doubly vergent orogen, the location of the main drainage divide is thought to represent a dynamic balance between (1) horizontal advection of topography that drives divides toward the retro side of an orogen (Fig. 1B; e.g., Hovius and Stark, 2006; Stolar et al., 2006; Willett et al., 2001), (2) orographic effects that localize precipitation on one side of a divide, driving it toward the dry side as catchments on the wet side expand (Fig. 1C; e.g., Bonnet, 2009; Hoth et al., 2006; Roe et al., 2003), and (3) local variations that affect erosion rate, such as differences in rock type, structural activity, or frequency of mass-wasting events that drive the divide away from the side with faster headward erosion (e.g., Kühni and Pfiffner, 2001; Pelletier, 2004; Stark, 2010; Willett et al., 2014).

In most active orogens, the effects of horizontal advection and topographic enhancement of precipitation are thought to dominate, leading to relatively predictable divide locations (e.g., Stolar et al., 2007; Willett, 1999). In the absence of strong orographic effects, models of doubly vergent orogens predict that the location of the drainage divide will mirror the asymmetric mass distribution that results from flow of material within the orogen, leading to a prowedge that is wider than the retrowedge (e.g., Willett, 1999; Willett et al., 1993, 2001). Localization of precipitation on one side of an orogen can either increase or reduce such asymmetry, depending on whether precipitation is localized on the prowedge or retrowedge, respectively, and the extent to which spatial patterns in deformation change as a result of increased precipitation (e.g., Hoth et al., 2006; Willett, 1999). In all of these cases, the drainage divide is predicted to be coincident with the topographic crest of the range, and thus still represents the division between the two distinctive structural regions, the prowedge and retrowedge (Fig. 1; e.g., Sanders et al., 1999).

The east-central Greater Caucasus provides an interesting exception to simple versions of standard wedge theory that predict that the topographic crest will coincide with the drainage divide and that both will be displaced toward the retrowedge side of an orogen with minimal cross-wedge precipitation gradients (Fig. 1). Here we define the east-central Greater Caucasus as the section of the range between 45°E and 49°E, characterized by both the Dagestan thrust belt on the north (retrowedge) side of the range (Sobornov, 1994) and the Kura fold-thrust belt on the south (prowedge) side (Forte et al., 2010, 2013), which is separated from the main range by the piggyback Alazani Basin, rendering it topographically distinct from the main range (Fig. 2C). This portion of the range is doubly vergent (Forte et al., 2014a; Mosar et al., 2010) and is the same region referred to as the “doubly vergent eastern zone” in Forte et al. (2014a). North-directed subduction and/or underthrusting of Kura Basin lithosphere beneath the Greater Caucasus is indicated by global positioning system (GPS) data (Reilinger et al., 2006; Vernant and Chery, 2006), seismicity (Mellors et al., 2012; Mumladze et al., 2015), and tomography (Skolbeltsyn et al., 2014), making the prowedge and retrowedge sides clear. We consider the topographic crest of the range as the division between the prowedge and retrowedge (see Forte et al., 2014a). Although the north flank is the retrowedge in this part of the range (Forte et al., 2014a; Mosar et al., 2010), the drainage divide is displaced southward toward, and within, the prowedge side of the orogen, and it does not coincide with the topographic crest, which is 20–40 km to the north and near the center of the range (Figs. 2 and 3; Forte et al., 2014a). In addition, the topographic asymmetry in the east-central Greater Caucasus departs from predictions of a doubly vergent wedge model, with the northern retrowedge having a shallower average surface slope than the main range within the southern prowedge (e.g., Forte et al., 2014a). However it is important to note that much of the width of the southern prowedge is represented by the Kura fold-thrust belt, which has a low topographic slope of ∼0.5° and in places is characterized by out-of-sequence deformation, leading to the hypothesis (Forte et al., 2013) that it may be in part subcritical. The prowedge and retrowedge in this portion of the orogen are nearly equal in width (∼130 km), when defined as the distance between the topographic limits of the surface expression of latest Quaternary active deformation and the topographic crest (Forte et al., 2014a). Thus, the topography across the east-central Greater Caucasus is decidedly different from that predicted by standard models of doubly vergent orogenic wedges (Figs. 1 and 3; e.g., Sanders et al., 1999; Willett et al., 1993).

In the case of the drainage divide position, this uncharacteristic asymmetry could theoretically be driven by orographic precipitation (e.g., Bonnet, 2009; Hoth et al., 2006; Roe et al., 2003) or differentials in base-level history on either side of the orogen. However, within the east-central Greater Caucasus, modern climatological data suggest no significant difference in mean annual precipitation on either side of the drainage divide (e.g., Forte et al., 2014a). Instead, the major precipitation gradient is along the strike of the range, driven by eastward transport of moisture-laden air from the Black Sea by the Westerlies (Borisov, 1965; Hijmans et al., 2005; Lydolph, 1977). Detailed paleoclimatological data are not available for the late Cenozoic for either the northern or southern side of the east-central Greater Caucasus, but reconstructions for the Lesser Caucasus to the south and the northwestern Greater Caucasus both show a trend toward more arid conditions ca. 2 Ma (e.g., Gabunia et al., 2000; Joannin et al., 2010; Kovda et al., 2008; Kvavadze and Vekua, 1993; Messager et al., 2010). At larger scales, reconstruction of regional paleoclimate during the late Cenozoic throughout the western Mediterranean and Central Asia suggest that during interglacials, the current climatic regime is a reasonable approximation for the past (e.g., Dodonov and Baiguzina, 1995; Youn et al., 2014). In contrast, during glacial periods, the northern part of the Westerlies, within the latitudinal range of the Caucasus, were partially disrupted by cold, dry air masses flowing southward from the Siberian high pressure system and off large continental ice sheets (e.g., Dodonov and Baiguzina, 1995), presenting the possibility for a more arid northern margin for the Greater Caucasus during glacial periods. Thus, neither the interglacial nor glacial period climatic regime of the east-central Greater Caucasus is consistent with a simple orographic explanation for this atypical topographic profile. Likewise, even though the Caspian has undergone extreme (>1 km) variations in water level (e.g., Forte and Cowgill, 2013; Popov et al., 2010; Zubakov, 2001), it is unlikely that the asymmetry results from greater base-level fall and headward erosion on the north side because base level for both sides of the east-central Greater Caucasus is set by the Caspian Sea. At present, the most reasonable explanation for the locations of the drainage divide and topographic crest in the east-central Greater Caucasus is that their positions are controlled by the geometry of active faults beneath the range, similar to the proposed structural control on the drainage network in the Swiss Alps (Kühni and Pfiffner, 2001). Spatially clustered zones of elevated normalized channel steepness (e.g., Kirby and Whipple, 2012; Wobus et al., 2006) and local relief within the east-central Greater Caucasus led us (Forte et al., 2014a) to conclude that the divide and crest each coincide with active loci of uplift in the interior of the range (Figs. 2 and 3).

Here we investigate the topographic and tectonic evolution of the east-central Greater Caucasus by integrating structural, stratigraphic and topographic observations with simple landscape evolution models (LEMs). The main goal is to determine if the two zones of active rock uplift inferred by us (Forte et al., 2014a) formed by southward or northward propagation. We show that current geologic observations are inconclusive because they are compatible with both scenarios, but that simple LEMs are most consistent with northward propagation. Thus, it appears that the topography of the east-central Greater Caucasus records a recent northward propagation of deformation from the southern flank to within the interior of the range. The timing of this northward propagation of deformation appears to be broadly concurrent with that of southward propagation of thrusting into the Kura Basin (e.g., Forte et al., 2010, 2013; Mosar et al., 2010), and we infer that the two zones of active rock uplift are above structures that are kinematically linked at depth to foreland deformation.

BACKGROUND

East-Central Greater Caucasus Tectonics and Topography

Between the Black and Caspian Seas, the Greater Caucasus Range represents the modern locus of northeast-southwest shortening within the central portion of the Arabia-Eurasia collision zone (Fig. 2B; Allen et al., 2004; Jackson, 1992; Reilinger et al., 2006). The crustal architecture and kinematics of the Greater Caucasus vary significantly along strike; the east-central Greater Caucasus is well described as a doubly vergent orogen, while the orogen both west of 45°E and east of 49°E is one sided and south directed (Forte et al., 2014a). Consideration of the east-central Greater Caucasus in the context of a doubly vergent orogenic wedge model with the southern and northern sides of the range being the prowedge and retrowedge, respectively, is consistent with modeling of GPS data (Reilinger et al., 2006; Vernant and Chery, 2006), relocation of deep earthquakes (Mellors et al., 2012), Rayleigh wave tomography (Skolbeltsyn et al., 2014), and synthesis of local earthquake catalogs (Mumladze et al., 2015), all which suggest north-directed subduction and/or significant underthrusting of Kura Basin lithosphere beneath the southern margin of the east-central Greater Caucasus to depths of at least 160 km (for a summary, see Mumladze et al., 2015). Convergence rates between the Kura Basin and stable Eurasia increase eastward from ∼8 mm/yr at 45°E to ∼14 mm/yr near the Caspian Sea coast (Kadirov et al., 2012; Reilinger et al., 2006), with most convergence equally partitioned into shortening along the southern and northern margins of the orogen within the east-central Greater Caucasus (Forte et al., 2014a).

Along the southern front of the main range, the locations, geometries, and activities of major thrusts within the east-central Greater Caucasus are not well constrained. A north-dipping thrust, typically referred to as the Main Caucasus thrust and abbreviated MCT, is often cited as the primary structure within the southern Greater Caucasus as a whole (e.g., Khain, 1975; Philip et al., 1989); however there is disagreement regarding its exact location and significance within the east-central Greater Caucasus. In this region, some consider the Main Caucasus thrust to be located near the topographic range front, and in part buried beneath sediments of the Alazani Basin (e.g., Forte et al., 2014a; Kadirov et al., 2012; Saintot et al., 2006), while others consider the Main Caucasus thrust to be farther north, with its surface trace in the low-relief area between the topographic crest and the drainage divide (Fig. 4; e.g., Adamia et al., 2011; Mosar et al., 2010). For our purposes, we consider the Main Caucasus thrust to be the range-front thrust system, which places Jurassic and Cretaceous flysch and carbonates over Jurassic and Cretaceous volcaniclastic rocks (Fig. 4). Quaternary shortening along the southern margin has been localized within the Kura fold-thrust belt, a series of predominantly south-vergent folds and thrusts deforming Pliocene–Pleistocene sediments, since at least ca. 2 Ma in the eastern portion of the belt (Fig. 2C; e.g., Forte et al., 2010, 2013). Geomorphic observations along the southeastern Greater Caucasus range front (Forte et al., 2010, 2014a; Mosar et al., 2010) and similarities between average rates of shortening from balanced cross sections within the Kura fold-thrust belt with modern geodetic rates of convergence (Forte et al., 2010, 2013) suggest that once deformation propagated into the southern foreland, active surface slip of the north-dipping thrusts within the main range ceased or slowed considerably.

Along the northern margin of the range, the Dagestan thrust belt is the principal active structural system (Sobornov, 1994). At the surface this system is characterized by north-dipping, south-directed thrusts (e.g., Dotduyev, 1986), inferred from analysis of seismic and well-log data to be backthrusts that sole into a roof thrust above a triangle zone that defines the leading edge of a north-vergent thrust system at depth (Fig. 2C; Sobornov, 1994, 1996). However, numerous regional maps also show the Dagestan belt as a relatively simple, north-vergent thrust system (e.g., Philip et al., 1989; Ruppel and McNutt, 1990). Unconformities and variation in stratal thicknesses suggest that the Dagestan thrust system initiated prior to deposition of Pliocene–Quaternary sediments, although the initiation age is not well constrained (Sobornov, 1996).

While active shortening in the Greater Caucasus is accommodated predominantly along the margins of the orogen, the topography of the east-central part of the range suggests that there are at least two range-parallel zones of active uplift within the interior of the mountain belt: a southern zone 30–40 km north of the northern margin of the Kura fold-thrust belt and coincident with the drainage divide, and a northern zone roughly along the topographic crest (Figs. 3 and 4). We (Forte et al., 2014a) found two prominent zones of elevated normalized channel steepness index values (ksn) that coincide with range-parallel bands of elevated local relief (Figs. 3 and 4), which is known to correlate well with channel steepness (e.g., DiBiase et al., 2010; Kirby and Whipple, 2012; Wobus et al., 2006). Accounting for variations in lithology or climate, ksn values generally scale with uplift and erosion rates (Kirby and Whipple, 2012; Wobus et al., 2006). The two zones of high ksn in the east-central Greater Caucasus likely mark zones of elevated rock uplift and erosion rates, because they do not correlate to significant changes in lithology or precipitation (Figs. 3 and 4; Forte et al., 2014a, fig. 7D therein).

Mechanisms for Divide Migration in the Caucasus

In the east-central Greater Caucasus, the geometry and location of the divide and crest appear to be controlled by zones of active uplift, as indicated by the separation between the divide and crest, their correspondence with regions of elevated ksn and local relief, and low ksn and local relief in the region between the divide and crest (Fig. 3). Thus, changes in the location of active structures within or beneath the range should have controlled both the organization of the drainage network and its temporal evolution. We envision two possible end-member relationships between the topographic crest, drainage divide, and zones of uplift, differing in the relative timing of uplift and propagation direction of deformation (Fig. 6). There is also the null hypothesis that there has been no significant change in the drainage network or location of active structures within the interior of the Greater Caucasus, essentially that the two zones of uplift initiated simultaneously.

In the first end member, deformation and uplift propagate northward relative to the upper plate, so that the topographic crest and its associated zone of uplift are younger than the divide. It is important to clarify that this scenario does not require a completely stationary drainage divide. As discussed in Forte et al. (2014a), the juxtaposition of high and low ksn values across the current drainage divide (Figs. 3 and 4) suggests that this feature is unstable, with higher erosion rates on the southern side of the divide, thus implying northward migration over time. However, such migration should be slow and more akin to scarp retreat than an abrupt shift in the divide position and topology of the drainage network stemming from formation of a new structure. Northward propagation could be accomplished through several nonunique sets of structures. One option is the formation of a new structure at depth beneath the range, such as an out-of-sequence duplex, possibly related to propagation of deformation into the Kura foreland basin. Formation of out-of-sequence hinterland duplexes during foreland propagation is observed in other thrust systems, and is interpreted as a deformational response to return the wedge to critical taper (e.g., Gutscher et al., 1998; Konstantinovskaya and Malavieille, 2005; Mitra and Sussman, 1997). An alternative is northward expansion of the uplift zone, possibly due to widening of the orogen or a change in the nature of the underthrusting lithosphere, e.g., a progressive transition from subduction to collision with a gradual increase in the amount of material accreted into the orogen.

In the second end member, the northern zone of uplift predated the southern zone and deformation migrated south relative to the upper plate. In this scenario, the topographic crest represents something like a paleo–drainage divide, with later uplift along the southern flank of the range decapitating formerly south-flowing rivers and reversing their flow direction between the divide and crest by integration into the north-flowing river systems. Such decapitation would produce a dramatic decrease in the length and drainage area of south-flowing rivers. From a structural perspective, southward propagation of deformation could suggest in-sequence thrust propagation and the initiation of a relatively discrete zone of uplift near the southern range front, but would depend on the age of this feature relative to the Kura fold-thrust belt. This structure could also represent reactivation or underplating in relation to a frontal accretion cycle (e.g., Gutscher et al., 1998; Hoth et al., 2007). We expect that this structure is in the form of a duplex or other blind structure, based on neotectonic observations indicating that there is no surface-breaking structure at the range front along the northern margin of the Alazani Basin (Forte et al., 2010).

POSSIBLE GEOLOGIC EVIDENCE OF SOUTHWARD DIVIDE MIGRATION

Results of landscape evolution modeling and additional geologic observations (discussed in the following) appear to refute southward propagation of uplift. However, an argument for such migration can be made based on the extents of modern and paleo–alluvial fans on the south flank of the east-central Greater Caucasus (Fig. 7). As we show here, fan length has decreased by a factor of 2–3 since ca. 2 Ma. Assuming other dominant variables were constant (e.g., subsidence rate), this reduction in fan length could indicate a reduction in drainage area for rivers draining the southern flank of the range, as predicted by southward propagation of uplift, formation of a new drainage divide, and decapitation of south-flowing rivers (Fig. 6C).

Modern alluvial fans in the area have an average length of ∼14 km and a maximum of 25 km (Figs. 7A, 7B). Previous work documented an empirical relationship between the area of an alluvial fan and the contributing drainage area within the upstream catchment, and suggested that the fan area is approximately half the catchment area (Fig. 7C; Dade and Verdeyen, 2007; Whipple and Trayler, 1996). Using a topographic contour map derived from the SRTM (Shuttle Radar Topography Mission, http://srtm.csi.cgiar.org/, ∼90 m pixel resolution) digital elevation model, we determined areas of 22 modern alluvial fans within the Alazani Basin and their associated contributing catchment areas in the southeastern Greater Caucasus (Fig. 7A). These fan and catchment areas closely follow the general empirical relationship found in diverse settings (Dade and Verdeyen, 2007), indicating that the fan depositional and catchment accumulation areas are well equilibrated (Fig. 7C).

Foreland-basin deposits now exposed in the Kura fold-thrust belt suggest Pliocene–Pleistocene alluvial fans extended 2–3 times farther from the Greater Caucasus into the Kura Basin than today (Forte et al., 2014b). In Forte et al. (2014b) exposures of Pliocene alluvial fan deposits in the Sarica and Vashlovani sections were described (Fig. 7A), and dated to ca. 2.5 Ma based on a population of detrital zircons from the Sarica paleofan. The distance of these sections from the range front provides minimum paleofan lengths (i.e., distance from the range front to the fan toe) of ∼20–36 km (Fig. 7A). These values are minima because they do not account for younger shortening within the Kura fold-thrust belt, which initiated ca. 2 Ma in the Göy area, ∼80 km to the east of Sarica (Fig. 7A; Forte et al., 2013). To reconstruct Pliocene fan lengths, we estimate ∼12 km of total postdepositional shortening between Sarica and the range front, based on shortening estimates from a balanced cross section near the Sarica area (western cross section in Forte et al., 2010). Thus, restored lengths of the paleofans were at least ∼32–48 km (Fig. 7B); however, this assumes that the topographic range front of the Greater Caucasus has not moved significantly since that time.

While paleofans approximately twice the size of modern are compatible with southward divide migration, such motion is not required because the relation between fan and catchment area is fundamentally modulated by various other factors, including sediment accumulation rate, rock uplift rate, basin geometry, drainage outlet and fan spacing, and properties of the sediment (e.g., Whipple and Trayler, 1996). As a preliminary test, we are interested whether the observed post–2.5 Ma reduction in fan length can be reasonably attributed to a major reduction in drainage area of the scale required for southward propagation and drainage beheading, keeping all other factors constant. Addressing this question requires estimating the axial lengths of the paleocatchments. To estimate this value, we first calculate a least squares regression between fan length and area for modern alluvial fans in the southeastern Greater Caucasus (Fig. 7B). We then make a simplifying assumption that the modern relationship was true for the paleofans; this allows us to estimate paleofan areas from the palinspastic reconstructions of paleofan length (Fig. 7B). We then use the reconstructed paleoalluvial fan areas and the general relationship between fan and catchment areas to estimate paleocatchment area (Fig. 7C; Dade and Verdeyen, 2007). To estimate the position of the paleo–drainage divide, we use the estimated paleocatchment area and an empirical relationship between catchment length and catchment area to estimate paleocatchment length, the mean of which is ∼40 km longer than modern catchments (Fig. 7D). We use a least-squares regression through lengths and areas of modern catchments in the southeastern Greater Caucasus that we extrapolate to the values spanned by the paleocatchments. However, in this case, the calculated linear relation between catchment length and area is not particularly robust (R2 = 0.6379). Such a fit is not surprising, given the variability of catchment shapes, but it provides an approximation of paleocatchment length. The 40 km reduction in mean catchment length over time is similar to the distance between the modern drainage divide and the topographic crest of the range (Figs. 7A, 7D). Thus, these data are compatible with the idea that the topographic crest may represent a paleo–drainage divide and that the reduction in drainage area was driven by a southward propagation of uplift that beheaded drainages in the southeastern Greater Caucasus. However, it is important to note that changes in precipitation, base level, local structural history, or accommodation space could also explain this decrease in alluvial fan area, thus a southward propagation of uplift is not a unique solution. Assuming constant amounts of accommodation space or uplift rates within catchments since 2.5 Ma during a period in which the foreland underwent large-magnitude base-level changes (e.g., Forte and Cowgill, 2013) and active propagation of the thrust front in the foreland basin (e.g., Forte et al., 2010, 2013) provides a reasonable level of doubt against the southward-propagation scenario. We use landscape evolution modeling to attempt to address this uncertainty.

LANDSCAPE EVOLUTION MODELING

To test the feasibility of both northward and southward propagation of uplift as mechanisms for producing the divergence between the crest and divide, we evaluated multiple variations of both scenarios using a landscape evolution model (LEM). The premise behind the use of a LEM is that the distinctive features of the current cross-sectional and plan-form topography of the east-central Greater Caucasus (Figs. 3 and 5) represent a primary constraint that either scenario must satisfy. Thus, it may be possible to distinguish between the two mechanisms by evaluating the extent to which each can faithfully generate the observed topographic patterns. As shown here, the model results suggest that northward propagation of uplift is the most feasible mechanism for generating the observed topography of the east-central Greater Caucasus.

Modeling Methodology and Parameterization

For the LEM, we use the FastScape model (Braun and Willett, 2013). This LEM is relatively simple in terms of the surface processes it simulates, compared to other popular LEMs (e.g., channel-hillslope integrated landscape development model, CHILD; Tucker et al., 2001). However, FastScape is well suited to the present problem because it allows simulation of landscapes of the spatial (thousands of square kilometers) and temporal (millions of years) scale of the Greater Caucasus in model runs that take less than an hour when utilizing the multithread capabilities of the code, thus permitting efficient testing of multiple scenarios. Also of particular importance is the capacity of FastScape to accommodate complicated uplift patterns that vary in time and space, although deformation is limited to vertical motion with no advection.

The east-central Greater Caucasus strike 115°–295° in a clockwise-positive azimuthal coordinate system in which north is 000°. All models are set up such that the positive x and y directions of the model domain coincide with the across-strike (025°) and along-strike (295°) directions in the east-central Greater Caucasus (see Supplemental Text1). For the majority of reported runs, the model domain is 260 km (x) by 500 km (y), which generally represents the dimensions of the east-central Greater Caucasus, i.e., the orogen is ∼260 km wide from the southern to northern margins of the Kura and Dagestan fold-thrust belts, respectively (Figs. 2, 4, and 7). The length of the model (y) is approximately half the length of the entire Greater Caucasus, but we focus our observations and discussion on the central 100–200 km of the model results to avoid edge effects. All models have a spatial resolution of 100 m, similar to the SRTM topographic data, and a time step of 100 k.y. Streams are forced to drain out of the northern (x = 260 km) and southern (x = 0 km) edges of the model domain, to replicate the Greater Caucasus. Base level (0 m), the constant of erosional efficiency (1.0E-5), and mean annual precipitation (1 m/yr) are held constant for all model runs. In addition, we do not consider flexural responses to erosion in any of the presented model runs.

We modeled 13 basic uplift scenarios to evaluate sensitivity to model parameters listed in Table 1, such as the width of the uplift zones and whether one or both of the uplift zones propagates laterally, along the strike of the orogen. Assigned values for variables, uplift functions, and final topographic outputs for all models are provided in greater detail in the Supplemental Text (see footnote 1). All scenarios have two primary peaks in uplift, which for most runs are located such that in cross section they occur in the same general areas as the peaks in uplift inferred from the high ksn values in the east-central Greater Caucasus (Table 1; Fig. 5; Supplemental Text [see footnote 1]). Models Dagestan-1 and Dagestan-2 are modified versions of the others designed to determine if the observed topography can be produced by an early single peak of uplift within the core of the range, followed by initiation of structures near the area corresponding to the Dagestan fold-thrust belt on the northern margin of the range. Models Stationary-1 and Stationary-2 explicitly test whether the topography of the east-central Greater Caucasus could be produced without either northward or southward propagation of deformation.

Other parameters that differ between the 13 base models are the geometries of the uplift-rate profiles, the relative timing of the two uplift zones, whether the uplifts initiate synchronously along strike or propagate laterally in the y direction (along strike), and the difference in maximum rate of the two zones of uplift (Table 1; Supplemental Text [see footnote 1]). Absolute magnitudes of uplift are largely unconstrained within the east-central Greater Caucasus; there are no published thermochronologic data for this portion of the range. The nearest available thermochronologic data are from the far eastern tip of the range (Avdeev, 2011; Fig. 3B), within a south-directed, singly vergent portion of the belt (Forte et al., 2014a). A pair of samples in the northern half of the range (N1 and N2, Fig. 3B) record average cooling rates of ∼10 °C/m.y. since 5 Ma, and a sample near the southern range front (V1, Fig. 3B) yields an average cooling rate of ∼15 °C/m.y. since 5 Ma (Avdeev, 2011). The geothermal gradient in this region is also not well constrained, but choosing standard values of between 20 and 30 °C/km suggests ranges of average uplift rates of 0.3–0.5 and 0.5–0.75 mm/yr since 5 Ma for the northern and southern samples, respectively. It is important to note that the geomorphic and geologic contexts of these samples are decidedly different than the area of interest for this study (Fig. 3B). These samples are from a portion of the range where the topographic crest and drainage divide are either coincident, or nearly so, and the orogenic width and elevation are reduced with respect to the east-central Greater Caucasus (Fig. 3). In addition, for the southern sample (V1), the structural geometry of this portion of the southern Greater Caucasus range front may differ substantially from areas of the range front farther west, as argued in Forte et al. (2013) based on a change in the topographic slope of the range and structural styles of the Kura fold-thrust belt. Thus, while it is not appropriate to simply translate these uplift rates westward, they constrain the magnitudes of uplift within our modeling study. We test models where the maximum magnitude of both uplift peaks is 0.4 mm/yr (South-1, South-2, and Stationary-2), the southern uplift is 0.3 mm/yr and the northern uplift is 0.4 mm/yr (Stationary-1, North-1, North-2, North-3, North-4, North-5, and South-3), the southern uplift is 0.4 mm/yr and the northern uplift is 0.3 mm/yr (Dagestan-1 and Dagestan-2), and one model (South-4) where the northern uplift is 0.4 mm/yr and the southern uplift is 0.8 mm/yr, slightly greater than the maximum average uplift rate from sample V1 (see Table 1).

Early tests of models with a flat initial topography produced extremely linear drainages, the forms of which were heavily influenced by the imposed uplift functions (Supplemental Figure S12). To produce more realistic drainage networks and to ensure that results were not biased by the initial drainage network form, all of the models were run using a common initial topography (Supplemental Text [see footnote 1] and Supplemental Figure S23). This initial topography was generated by letting a model of the same dimensions (x = 260 km, y = 500 km) evolve under a constant uplift rate of 0.01 mm/yr for 100 m.y. (Supplemental Figure S2 [see footnote 3]). Drainages in this model were forced to exit along the southern edge (x = 0). A one-sided initial model geometry was chosen to avoid having a preexisting drainage divide within the model that might influence later results. With the low uplift rate, the topography of this initial model was subdued, with a maximum elevation of ∼17 m above base level. It is important to note that this initial model topography is not meant to represent anything geologically specific to the Greater Caucasus, but only to provide the ensuing models with a preexisting dendritic channel network. In our analysis and discussion of model results, we consider only the dendritic models, but comparison of a representative set of base models (those without the dendritic seeds) and the presented models (those with the dendritic seeds) suggests that while the drainage networks of the presented models appear more realistic, the general topographic forms of the two sets of results are not significantly different (Supplemental Figure S1 [see footnote 2]).

In constructing and running the models, we also experimented with different erosion laws, because FastScape can be run using a unit stream power law with a linear slope dependence (n = 1), general unit-stream power law (n ≠ 1), or the ξ – q model of Davy and Lague (2009), where ξ is the ratio between sediment river load and the depositional flux and q is river discharge. The latter can approximate transport-limited behavior by adjusting the alpha parameter within the FastScape model. In general, these different erosion laws did not dramatically influence the results, so the presented models are those run with a general unit-stream power law with values of m = 0.4 (exponent on drainage area in stream power law) and n = 0.8. We selected values of m and n such that the ratio of m/n = 0.5, consistent with empirical results indicating that m/n ∼ 0.5 is required to the match of topography of most ranges (e.g., Whipple and Tucker, 1999). We used a value of 0.8 for the slope exponent n as this was between values of n in which erosion rate is proportional to shear stress (n = 2/3) and unit stream power (n = 1).

For simplicity, we discuss two exemplar models in detail, models South-3 and North-5, but we present the final outputs of all models in the Supplemental Text (see footnote 1) (Supplemental Figure S34). Both models were run for 10 m.y. and uplift rates for the first uplift peak increase incrementally over 2 m.y. before reaching final constant values, to partially simulate the initial slow exhumation phase observed in the Greater Caucasus (Avdeev, 2011; Avdeev and Niemi, 2011). However, as discussed in more detail in the following, because we do not have adequate constraints for absolute uplift magnitudes, erosional efficiency or precipitation, we do not consider times in any of the models to be realistic with respect to timing of events within the Greater Caucasus. We therefore limit our interpretations and discussion to the relative sequences of events. For both models, we use FastScape to produce digital elevation models for the final time step, which we then use to calculate normalized channel steepness values by using the TopoToolbox-2 (Schwanghart and Scherler, 2014). In this analysis, we use the same parameters of the east-central Greater Caucasus (reference concavity of 0.45) as in Forte et al. (2014a), but with a larger accumulation area to limit the number of streams.

North-5 tests a northward propagation of deformation. An initial southern zone of uplift (U1) initiates with its peak near the location of the drainage divide (Fig. 8A). This uplift starts at time 0 and increases in a stepwise fashion to one-half, three-quarters, and then full magnitude over 2 m.y., while also propagating in the east-west direction (y) as described in the Supplemental Text (see footnote 1). The northern uplift zone (U2) initiates 2.5 m.y. into the run and propagates in the east-west direction (y) for 2.75 m.y. until it reaches the edges of the model domain. The model then runs for another 4.75 m.y. The northern (U2) and southern (U1) uplifts have maximum values of 0.4 mm/yr and 0.3 mm/yr, respectively.

South-3 tests a southward-propagation scenario. The first uplift (U1) is broad and centered in the approximate location of the topographic crest (Figs. 5 and 7B). This uplift starts at time 0 and increases in a stepwise fashion to one-half, three-quarters, and then full magnitude over 2 m.y., while also propagating in the east-west direction (y) as described in the Supplemental Text (see footnote 1). At the end of 2 m.y., the uplift rate function for the northern peak reaches its maximum and remains constant for the rest of the run. At 2.5 m.y. into the run, the second uplift (U2) initiates to south of U1 and propagates in the east-west direction (y) for 2.75 m.y. until it reaches the edges of the model domain. The uplift function then remains unchanged for another 4.75 m.y., until the end of the model run. The northern (U1) and southern (U2) uplifts have maximum values of 0.4 mm/yr and 0.3 mm/yr, respectively.

Modeling Results

Models with northward propagation most successfully produce the separation of the topographic crest from the drainage divide seen in the east-central Greater Caucasus (Figs. 8 and 9). More recent uplift in the northern zone successfully creates high peaks in the interfluves while maintaining transverse drainages. In contrast, scenarios with younger uplift in the south typically fail to decapitate drainages and force the drainage reversals required to generate a jump in divide position. In addition, lateral propagation of a young northern uplift is needed to establish longitudinal tributaries like those observed in the Sulak and Samur drainages (Fig. 2), though significantly curved drainages like the Samur only occur at the very edges of the model domain in most model results. Models Stationary-1 and Stationary-2, which test a scenario with no propagation of deformation, fail to produce a separation between the drainage divide and topographic crest. After 10 m.y., the run time for the majority of other models (Table 1), these two models had produced a continuous elevated region between the two uplift zones, so they were allowed to run for an additional 10 m.y. For Stationary-1, where the southern uplift zone is 0.3 mm/yr and the northern is 0.4 mm/yr, at the end of the 20 m.y. run the divide and crest are merged and located near the northern zone of uplift (Supplemental Figure S3 [see footnote 4]). The final topography of Stationary-2, which has equal magnitudes for the two uplift zones (0.4 mm/yr), also has a coincident divide and crest, but with the location of this topographic feature varying between the northern and southern zones of uplift (Supplemental Figure S2 [see footnote 3]).

The two models Dagestan-1 and Dagestan-2 designed to test the influence of the Dagestan fold-thrust belt also fail to produce topography analogous to that of the east-central Greater Caucasus (Supplemental Text [see footnote 1]). The only southward-propagation models that partially resemble the east-central Greater Caucasus are ones in which the magnitude of the southern uplift zone is at least twice that of the northern zone of uplift, e.g., South-4 (Supplemental File5). In this specific case, the later development of a new, southern zone of uplift causes a drainage reversal, but also causes the new drainage divide to quickly become the topographic crest, with peak elevations significantly higher than those above the northern uplift. Such topography is inconsistent with that observed in the east-central Greater Caucasus.

The exemplar northward propagation model, North-5, produces patterns of normalized channel steepness, topography, and longitudinal profiles that are generally similar to those seen in the east-central Greater Caucasus (Figs. 8 and 9). It is important to note that the absolute values of ksn and topography in the east-central Greater Caucasus and the model outputs are different, thus we only compare the spatial patterns between high and low ksn values and topographic form (Fig. 8). A difference between absolute values of ksn is not unexpected, as the models are not tuned to match the Greater Caucasus because the information required to do so is either unknown (in the case of the erosional efficiency) or lacks the spatial and temporal certainty required for accurate implementation in the model (in the case of the uplift rate and precipitation).

Comparing the synthetic topography in North-5 and the real range, both have high ksn values south of the drainage divide and north of the topographic crest, with lower ksn values in the intervening region. A swath topographic profile oriented across strike through the North-5 model results (Fig. 9B) further illustrates the similarity between the topography produced in this model and the topographic form of the east-central Greater Caucasus (Fig. 5). In addition, longitudinal profiles of selected drainages (Fig. 9B) appear similar in form to south- and north-flowing drainages within the east-central Greater Caucasus (Fig. 5). In contrast, the exemplar southward propagation model, South-3, has high ksn values on either side of the drainage divide, which is also the topographic crest, and a third zone of high ksn values in the vicinity of the southern uplift. The topography of the South-3 results in both map (Fig. 8E) and cross section (Fig. 9B) view differ significantly from the observations within the east-central Greater Caucasus, including a significantly diminished topographic expression of the southern peak and more symmetrical river longitudinal profiles between north and south flowing drainages.

Two other models, North-3 and North-4, have similar uplift geometries as in North-5, but the northern uplift initiates later in the model run, starting 5.5 m.y. after the start as opposed to 2.5 m.y. (Table 1). In addition, the geometry of the southern uplift was wider in North-3 and North-4 (Supplemental Text [see footnote 1]). These two models also successfully produce a separation between the drainage divide and crest, but the geometry of the northern drainage network did not mirror the northern east-central Greater Caucasus as well as North-5, i.e., longitudinal rivers like the Samur did not form (Supplemental Text [see footnote 1]). This further indicates that relatively discrete uplift peaks with a lower uplift zone between them are required to faithfully reproduce the east-central Greater Caucasus topography. It may also indicate that for drainages like Samur and its tributaries to form, initiation of this northern uplift peak after the southern is required, but still relatively early within the history of the range. However, the time scales of the models are not well calibrated to the actual Greater Caucasus.

We also attempt to use the LEM results to distinguish between two tectonic models for the northward-propagation scenario (Fig. 10). The progressive widening of the orogen could produce an expansion of the southern zone uplift, forming a continuous zone of uplift roughly between the drainage divide and the topographic crest (e.g., North-1 or North-2, Fig. 10; Supplemental Text [see footnote 1]), whereas the formation of a hinterland duplex or reactivation of structures might produce a more discrete zone of uplift (e.g., North-3, North-4, North-5, Dagestan-1, Dagestan-2; Supplemental Text [see footnote 1]). Models with a continuous zone of uplift fail to provide a separation between the topographic crest and the drainage divide, and generally produce topography inconsistent with that of the east-central Greater Caucasus. Therefore, the scenario modeled in North-5 best reproduces the modern topography (Figs. 8 and 10).

DISCUSSION

Reconciling Conflicting Results

Results of the landscape evolution modeling suggest that northward propagation of uplift within the interior of the east-central Greater Caucasus is most consistent with the current topography of the range (Figs. 8 and 9). However, apparent reductions in size of paleoalluvial fans based on exposures within the Kura fold-thrust belt stratigraphy are consistent with decapitation of south-flowing drainages and a reduction in drainage area that could be driven by a southward propagation of uplift (Fig. 7). We consider the results of the landscape evolution modeling, and thus the northward propagation of uplift, to be more robust because several alternate solutions can explain the reduction in alluvial fan size along the southeastern Greater Caucasus.

While the approximate twofold decrease in alluvial fan length since 2.5 Ma could result from a structurally controlled reduction in drainage area, it could also be explained by an increase in subsidence rates during growth of the Greater Caucasus, an increase in aridity leading to decreased sediment flux, or a dramatic rise in base-level (e.g., Clevis et al., 2003; Dade and Verdeyen, 2007; Whipple and Trayler, 1996). Natural examples and numerical models of alluvial fans indicate that the amount of accommodation space exerts a primary control on alluvial fan size (e.g., Clevis et al., 2003; Dade and Verdeyen, 2007; Whipple and Trayler, 1996). The amount of available accommodation space (independent of sediment flux into a basin) will be a combination of subsidence rate and regional base level. In general, there is an inverse correlation between alluvial fan size and the amount of accommodation space, because an increase in accommodation causes the same volume of material delivered to a basin to be distributed over a smaller area, thereby reducing fan size. In the case of the southeastern Greater Caucasus, an increase in average subsidence rate within the foreland basin likely occurred coincident with, or shortly after, the onset of rapid exhumation of the Greater Caucasus and presumed growth of topography. An important additional consideration is the formation of the Kura fold-thrust belt as both a barrier to sediment dispersal from the Greater Caucasus and a driver of local subsidence in the Alazani Basin. Timing of fold-thrust belt formation in this region is not well constrained, but the initiation ages of 2 and 1.5 Ma in the Göy area (Forte et al., 2013) suggest that preliminary development of the Kura fold-thrust belt could also influence the extent of alluvial fans in the southeastern Greater Caucasus. The initiation of the Kura fold-thrust belt may have also led to some amount of range front retreat, further contributing to an apparent reduction in alluvial fan length. The southeastern Greater Caucasus range front on the north side of the Alazani Basin is extremely embayed; this is interpreted to indicate cessation or significant slowing of activity on the surface trace of the Main Caucasus thrust (Fig. 6A; Forte et al., 2010). When the surface trace of the Main Caucasus thrust was active, the range front was likely coincident with the trace of this thrust, which throughout much of the east-central Greater Caucasus is 5–10 km south of the modern range front and buried under sediments of the Alazani Basin (Fig. 2C).

Decreases in precipitation can also decrease alluvial fan size during transient response of a landscape, assuming that rates of subsidence or creation of accommodation space remain constant (e.g., Clevis et al., 2003; Whipple and Trayler, 1996). Within the Caucasus, there was a regional shift to more arid conditions ca. 2 Ma (e.g., Gabunia et al., 2000; Joannin et al., 2010; Kovda et al., 2008; Kvavadze and Vekua, 1993; Messager et al., 2010); however, the scale of this aridification is unclear. Also, in general, climatic variations probably do not exert enough control to effect such a significant reduction in fan size once a landscape has returned to steady state (e.g., Clevis et al., 2003; Dade and Verdeyen, 2007; Whipple and Trayler, 1996), thus the importance of climate in reducing fan size in this case is uncertain, but still potentially relevant.

Finally, an increase in regional base level could drive this apparent reduction in alluvial fan size by effectively increasing the accommodation space closer to the range front. In Sarica and Vashlovani, the deposits immediately overlying this alluvial fan complex are interpreted to have been deposited in a lacustrine setting that correlates with the Akchagyl regional stage (Forte et al., 2014b), which represents an extreme sea-level highstand in the Caspian that submerged much of the Kura Basin (e.g., Forte and Cowgill, 2013; Popov et al., 2006).

Interpreting the LEM Results

The LEM results provide good evidence that the current topography of the east-central Greater Caucasus results from northward propagation of uplift. The simplest interpretation of this finding is consistent with previous work that suggests it is difficult to modify an entrenched, antecedent drainage network once it is established (e.g., Castelltort and Simpson, 2006; Kühni and Pfiffner, 2001; Oberlander, 1985). In addition, the modeling suggests that a simple northward expansion of the zone of uplift does not reproduce the relief structure or drainage patterns of the region between the topographic crest and drainage divide; we therefore favor a model with two discrete zones of uplift (Fig. 10).

However, the caveats and assumptions of the modeling approach must be considered. One potentially limiting factor is the lack of data to constrain input variables. While the models were designed to represent the east-central Greater Caucasus, they are not direct simulations because we have no measurements of actual uplift rates or field calibration for variables used in the erosion law (i.e., erosional coefficient, exponents m and n for drainage area and slope terms). Such uncertainties should have secondary impact on formation of the primary topographic patterns, such as the separation between the drainage divide and topographic crest, although they are likely essential for understanding the rates and time scales over which these features form. A second restriction is the lack of advection in FastScape, because horizontal motion and resulting advection of topography are important in simulations of both individual folds (e.g., Hardy, 1995, 1997; Waltham and Hardy, 1995) and orogen-scale deformation (e.g., Braun and Yamato, 2010; Willett et al., 2001). Although we attempted to approximate transport-limited erosion by performing some runs with the FastScape implementation of the Davy and Lague (2009) erosion model, no sediment is deposited between the two uplift zones in runs using this erosion law, which we speculate is due to elevated uplift rates within this zone. A true transport-limited erosion law might not exhibit the same behavior and could be important in allowing models with southward-propagation to produce the observed topography. Transport-limited erosion laws specifically include sediment flux, and consider the role of sediment in both protecting the river bed from incision and enhancing incision through abrasion or saltation; such laws may better capture the behavior of bedrock rivers under certain conditions (e.g., Gasparini et al., 2007; Sklar and Dietrich, 2006, 2008). Transport-limited erosion in south-flowing rivers between the two zones of uplift could help facilitate drainage reversal. In this case, younger initiation of the southern uplift could cause deposition and a reduction in gradient in channel segments between the two zones of uplift, thereby facilitating headward erosion by north-flowing drainages and eventual capture and drainage reversal. However, if this occurred in the east-central Greater Caucasus we would expect remnants of alluvial deposits in the low-relief area between the drainage divide and topographic crest. No such deposits are shown on geologic maps of this region, although the maps are large-scale overviews (e.g., Nalivkin, 1976).

Implications for East-Central Greater Caucasus Tectonics

The initiation of an uplift zone coincident with the topographic crest after the formation of the more southern uplift zone along the drainage divide implies out-of-sequence deformation within the east-central Greater Caucasus. This younger zone of uplift roughly corresponds with the area highlighted by Mosar et al. (2010) as a zone of increased uplift within the core of the east-central Greater Caucasus; however, we do not consider this zone to be bounded by active thrust as interpreted by Mosar et al. (2010). The gradual, decreasing gradients in ksn on either side of the high ksn and local relief areas near the topographic crest are not consistent with predictions for a discrete, surface breaking fault (e.g., Kirby and Whipple, 2012). This strongly suggests that neither the southern range front fault, e.g., our Main Caucasus thrust, nor faults within the interior of the range near the topographic crest or drainage divide, e.g., the Main Caucasus thrust of Mosar et al. (2010), actively break the surface. Instead, we hypothesize that motion on the structures associated with both zones of uplift are blind and transfer slip into the foreland fold-thrust belts along the margins of the range (Fig. 2C). Thus, these zones of uplift appear to illuminate the subsurface locations of geometric complexities (e.g., ramps, duplexes, or blind splay faults) along the main thrusts or shear zones that underlie the east-central Greater Caucasus. However, the topographic observations and LEM results do not provide a clear indication of the subsurface geometries or vergence directions of the structures producing the uplift. In the case of the southern uplift zone coincident with the drainage divide (as noted in Forte et al., 2014a), the prominent northward increasing ksn values along the southern range front (Figs. 3 and 4) are most consistent with a similar northward-increasing gradient in uplift rate (e.g., Gasparini and Whipple, 2014). This gradient in uplift rate could be consistent with growth of a duplex underneath the southern range front and would most likely be a south-vergent set of structures, given the structural vergence within the Kura fold-thrust belt and structures along strike within the southeastern Greater Caucasus. If either a duplex or ramp is causing the southern uplift zone, the evidence for the current inactivity of the range-front thrust systems suggests that these structures may be in the process of being deformed by movement on this blind structure (Fig. 11). The geometry of the northern uplift is less clear and could be above either (1) an additional south-vergent duplex or ramp within the north-dipping thrust system that appears to link the Main Caucasus thrust and the Kura fold-thrust belt, or (2) a north-vergent structure within with south-dipping décollement of the Dagestan thrust belt. From the available data, we can also not eliminate the possibility of a high-angle structure within the interior of the range, similar to that depicted by Mosar et al. (2010), but if such a structure exists, the ksn data suggest that it does not break the surface.

At this point we favor a south-vergent geometry for the northern uplift, as geologic data attest to recent expansion of the south-directed Kura fold-thrust belt (Forte et al., 2010, 2013, 2014a; Mosar et al., 2010), and the Dagestan thrust belt along the northern margin of the range appears less active in comparison. Although we suggested (Forte et al., 2014a) that ∼50% of the geodetically measured convergence across this portion of the Greater Caucasus is accommodated on the north flank of the range, we also found relatively low ksn values for rivers in this area (Figs. 3 and 4); this may suggest minimal activity in the Dagestan thrust belt, but would also be consistent with movement on shallow-dipping thrusts, consistent with geometries illustrated by Sobornov (1996; Fig. 11). Mainstem channels have high ksn where they cross the Dagestan thrust belt, but these are generally correlated with either lithologically supported knickpoints (e.g., knickzones occurring at transition between overlying carbonate and underlying clastics) or anthropogenic features (e.g., Chirkey Dam, Fig. 3; Forte et al., 2014a). More work focused on the neotectonics of the Dagestan thrust belt is needed to resolve this question.

Reactivation of structures within the hinterland of an orogen can also occur during an accretion cycle, which is when a discrete thrust sheet is added to the orogen, causing the deformation front to expand (e.g., Gutscher et al., 1998; Hoth et al., 2007). The formation of the Kura fold-thrust belt may reflect such an accretion cycle. To estimate the time scale of a full accretion cycle we use the following expression from Hoth et al. (2007): 
graphic
in which accretion time scale (t) is a function of the décollement depth (D) and convergence velocity (v). Using a detachment depth of 10 km (Fig. 11) and a convergence velocity for this portion of the orogen of between 6 and 8 mm/yr (Forte et al., 2014a; Reilinger et al., 2006) suggests that the time scale of a full accretion cycle is 5–6 m.y. This time scale is similar to the age of the Kura fold-thrust belt (ca. 3–2 Ma; Forte et al., 2010, 2013), suggesting that active uplift within the interior of the range could result from formation of the Kura fold-thrust belt during an accretion cycle. Analogue models indicate that such hinterland activity can be expressed as either reactivation of older structures (e.g., Hoth et al., 2007) or formation of new hinterland duplexes during underthrusting (Gutscher et al., 1998). The latter case is more consistent with observations noted here that indicate a lack of evidence for surface-breaking structures in the interior of the range or along the northern margin of the Alazani Basin, in the east-central Greater Caucasus.

An alternative option, given the location of the northern uplift zone near the center of the range and in the vicinity of deep (>50 km) earthquakes interpreted as being related to subduction and/or underthrusting of Kura Basin lithosphere (e.g., Forte et al., 2014a; Mumladze et al., 2015), is that the formation of this uplifting zone is related to changes in the subsurface geometry of the subducting and/or underthrusting lithosphere, possibly as increasingly thick, and more buoyant lithosphere is underthrust (Fig. 11). Distinguishing whether the activity of these hinterland structures is related to a structural change within the orogen or simply part of an accretion cycle requires better documentation of both the geometries and histories of these zones and the responsible structures. Regardless of the exact cause of these zones, active uplift in the core of the range is consistent with results from low-temperature thermochronometry, indicating that the youngest thermochronometric ages, and thus the most rapid exhumation, occur in the center of the range (e.g., Avdeev, 2011; Avdeev and Niemi, 2011).

Drainage Networks and Active Deformation

The LEM results presented here are consistent with previous work on the interaction between developing orogens and rivers that suggest it is very difficult to modify antecedent drainage networks (e.g., Castelltort and Simpson, 2006; Kühni and Pfiffner, 2001). Similar to the findings presented by Kühni and Pfiffner (2001), our results suggest that the drainage network of orogens (e.g., divide positions) is strongly influenced by the location of active structures during early orogenesis. Considering this example in the context of the expectations of models of doubly vergent orogenic wedges, the east-central Greater Caucasus demonstrate that local variations in structural geometry can strongly influence the topographic form of a doubly vergent orogen.

We suggested (Forte et al., 2014a) that the Greater Caucasus may represent an orogen in the early stages of transition from single to double vergence. The form of the drainage network may be similarly transitional, and evolving toward a topology that is more consistent with a doubly vergent orogen. As highlighted in Forte et al. (2014a), the juxtaposition between high and low ksn values to the south and north of the drainage divide, respectively, is not stable (e.g., Willett et al., 2014). The drainage divide will only maintain its position as long as the structure beneath it is active. If activity diminished on this structure, we speculate that the drainage divide would migrate northward, and assuming no further change in structural geometries, will eventually be more closely associated with the topographic crest of the range. The rate of this scarp retreat would depend on the rate of structurally driven uplift, the erosive power of the rivers draining the southeastern Greater Caucasus, and the erodibility of the different rock types in these catchments. This would also potentially bring the overall topographic profile of the east-central Greater Caucasus closer to the predictions of a doubly vergent wedge model, with a wider, more shallowly sloping prowedge in the south and, narrower, steeper retrowedge in the north and a coincident drainage divide and topographic crest. However, the reasons for the similar widths and uncharacteristic topographic form of the prowedge and retrowedge within the modern Greater Caucasus remain an active area of research.

CONCLUSION

Coupling topographic, structural, and stratigraphic observations from the east-central Greater Caucasus with results of landscape evolution modeling suggests recent structural reorganization or reactivation within the interior of this portion of the orogen. While the loci of active shortening structures have progressively stepped southward, largely abandoning the Main Caucasus thrust and other north-dipping thrusts within the main range to form the Kura fold-thrust belt, results presented here suggest that active uplift within the interior of the range also stepped northward. LEMs that best reproduce the relief and drainage structure of the east-central Greater Caucasus involve two discrete zones of uplift, with a southern zone near the modern-day range front initiating first, followed by the initiation of a northern zone of uplift, near the modern-day topographic crest. The available data do not allow us to define subsurface structural geometries beneath these zones of elevated uplift rate, but the simplest explanation is formation of south-vergent duplexes or developments of ramps at depth that feed slip to active shortening structures in the Kura fold-thrust belt. While we consider the LEM results robust, they require field verification, especially considering the stratigraphic evidence from the Kura fold-thrust belt that could indicate a southward propagation of uplift. More generally, our results also demonstrate that local structural geometry can dramatically influence the topographic evolution of an orogen, specifically related to maintaining the location of drainage divides.

We thank Jean Braun for providing the source code for the FastScape landscape evolution model and Byron Adams for discussions regarding modeling methodology. We also thank editor Raymond Russo and reviewer Silvan Hoth for helpful comments that improved this manuscript. This material is partially based upon work supported by National Science Foundation grant EAR-0810285 (Cowgill) and a School of Earth and Space Exploration Postdoctoral Fellowship at Arizona State University (Forte).

1Supplemental Text. Additional details about the methods used for the paper. Please visit http://dx.doi.org/10.1130/GES01121.S1 or the full-text article on www.gsapubs.org to view the Supplemental Text.
2 Supplemental Figure S1. Comparing the final time-steps of models North-5 and South-3 using a flat initial topography (left), versus a low-relief landscape with a south-flowing dendritic channel network, shown in Supplemental Figure S2 (right). For the majority of the main text and supplement, we only discuss the models using the initial dendritic channel network topography. Please visit http://dx.doi.org/10.1130/GES01121.S2 or the full-text article on www.gsapubs.org to view Supplemental Figure S1.
3 Supplemental Figure S2. Final topography after running a model for 100 million years of dimensions x = 260 km and y = 500 km that drains to the south (x = 0) with an uplift rate of 0.01 mm/yr. This final topography was used as the initial topography for all of the presented models. Please visit http://dx.doi.org/10.1130/GES01121.S3 or the full-text article on www.gsapubs.org to view Supplemental Figure S2.
4 Supplemental Figure S3. Schematic representation of the uplift rates imposed for the various model scenarios tested, viewed in north-south orientation with respect to the Greater Caucasus. Black lines indicate initial geometries, grey lines indicate portion of the geometry that is different after initiation of second uplift zone. Concentric circle symbol indicates lateral propagation was implemented and dotted lines indicate the partial uplift rates implemented during these lateral propagation periods. Relative differences between the maximum value of the two uplift rates, U1 and U2 are given. Geometry of uplift rate gradient is labeled as either, “c” for a rate that decreases with the cube of distance, “s” for a rate that decreases with the square of distance, or “l” for a linear decrease with distance as described in the supplemental text. Below each uplift function is the output topography of the final timestep of the model. Final panel shows schematic of time-steps in the lateral propagation function. For this panel, the view is 90 degrees from the other diagrams, approximately east-west in terms of the Greater Caucasus. Please visit http://dx.doi.org/10.1130/GES01121.S4 or the full-text article on www.gsapubs.org to view Supplemental Figure S3.
5Supplemental File. Zipped file containing FastScape model input files for models discussed in the text and the Supplemental Text. Please visit http://dx.doi.org/10.1130/GES01121.S5 or the full-text article on www.gsapubs.org to view the Supplemental File.