Reconstruction of river valley evolution before and after the emplacement of the giant Seymareh rock avalanche (Zagros Mts., Iran)

The Seymareh landslide, detached ∼10 ka from the northeastern flank of the Kabir-kuh fold (Zagros Mts., Iran), is recognized worldwide as the largest rock slope failure (44 Gm³) ever recorded on the exposed Earth surface. Detailed studies have been performed that have described the landslide mechanism and different scenarios have been proposed for explaining the induced landscape changes. The purpose of this study is to provide still missing time constraints on the evolution of the Seymareh River valley, before and after the emplacement of the Seymareh landslide, to highlight the role of geomorphic processes both as predisposing factors and as response to the landslide debris emplacement.

We used optically stimulated luminescence (OSL) to date lacustrine and fluvial terrace sediments, whose plano-altimetric distribution has been correlated to the detectable knickpoints along the Seymareh River longitudinal profile, allowing the reconstruction of the evolutionary model of the fluvial valley. We infer that the knickpoint migration along the main river and the erosion wave propagation upstream through the whole drainage network caused the stress release and the ultimate failure of the rock mass involved in the landslide. We estimated that the stress release activated a mass rock creep (MRC) process with gravity-driven deformation processes occurring over an elapsed time-to-failure value on the order of 10² kyr. We estimated also that the Seymareh damming lake persisted for ∼3500 years before starting to empty ∼6.6 ka due to lake overflow. A sedimentation rate of 10 mm yr⁻¹ was estimated for the lacustrine deposits, which increased up to 17 mm yr⁻¹ during the early stage of lake emptying due to the increased sediment yield from the lake tributaries. We calculated an erosion rate of 1.8 cm yr⁻¹ since the initiation of dam breaching by the Seymareh River, which propagated through the drainage system up to the landslide source area.

The evolutionary model of the Seymareh River valley can provide the necessary constraints for future stress–strain numerical modeling of the landslide slope to reproduce the MRC and demonstrate the possible role of seismic triggering in prematurely terminating the creep-controlled time-to-failure pathway for such an extremely large case study.

1 Introduction

Tectonically active landscapes are very dynamic systems in which threshold conditions on hillslopes are often reached. Therefore, there are considerable implications for natural hazards related to seismicity and to the geomorphic coupling between hillslopes and rivers, with both fluvial control on hillslopes and landslides affecting the fluvial network. In response to rock uplift, relief and hillslope angles increase linearly in time mainly due to fluvial erosion processes in landscapes affected by low to moderate tectonic forcing (Montgomery and Brandon, 2002; Binnie et al., 2007; Larsen and Montgomery, 2012). Nonetheless, such a linear increase in relief and hillslope angles is limited by the reaching of the threshold slope conditions associated with the hillslope material strength (Schmidt and Montgomery, 1995), until the latter is exceeded by gravitational stress giving rise to bedrock landslides. This leads to a nonlinear increase in erosion rates in landscapes affected by long-lasting or high-rate tectonic forcing, where the increase in the rate of channel incision is accommodated by an increased frequency of slope failure rather than by slope steepening. These slope failures may represent the terminal phase of mass rock creep processes (MRC; Chigira, 1992) contributing to the development of gravity-driven deformations due to time-dependent rheology. These kinds of deformations can evolve into sudden and complete slope collapse when the increased strain rate leads to an accelerating deformation (Saito, 1969; Petley and Allison, 1997). The strain limit conditions are reached when the stationary creep stage evolves into the accelerating creep stage (Saito, 1969; Petley and Allison, 1997). Failure is typically associated with mass strength reduction as an effect of rock mass damage occurring over time during the ongoing gravity-driven deformations (Eberhardt et al., 2004; Stead et al., 2006). Catastrophic failures occurring in the accelerating MRC stage can evolve rapidly into rock avalanches due to the fragmentation of rock masses during their transport (Hungr et al., 2001; Davies and McSaveney, 2009).

This work is focused on the evolution of the Seymareh River valley before and after the valley was dammed by the landslide that is recognized worldwide as the largest subaerial landslide ever observed, the Seymareh rock avalanche (Roberts and Evans, 2013) that occurred in the northwestern Zagros Mts. (Iran). Different interpretations have been proposed so far by the scientific community to explain the generation of such an exceptional event and different scenarios have been hypothesized for explaining the induced changes in landscape. Harrison and Falcon (1937, 1938) provided much of the present knowledge on the rock avalanche, including the geology and structure of the source area, the general geomorphology and the basic geometry of the landslide. Oberlander (1965) included a short appendix on the landslide in his study of Zagros streams and discussed its origin in relation to the activity of the Seymareh River. Later, in the 1960s, Watson and Wright (1969) characterized the geomorphology and stratigraphy of the debris, discussed the origin of the initial rockslide and examined the debris avalanche emplacement mechanisms. Roberts (2008) and Roberts and Evans (2013) provided a detailed model of how the geological and tectonic evolution of the Kabir-kuh fold predisposed the slope to such a large-scale failure, including formation of structural or kinematic and rheological control, and inferred a seismic trigger. Specifically, Roberts and Evans (2013) obtained a ¹⁴C age of 8710 years BP from a charcoal-rich layer approximately 15 m above the base of the lacustrine sequence. Based on the interpretation of three separate radiocarbon ages provided additionally by Griffiths et al. (2001), an estimated radiocarbon bracket age of the Seymareh event was suggested between 9800–8710 ¹⁴C years BP. Yamani et al. (2012) provided some general details on the evolution of the dammed lake drainage, describing a sequence of entrenched lacustrine terraces upstream of the landslide dam. Finally, Shoaei (2014) reviewed the possible mechanisms of failure and interpreted the post-failure geomorphic features through analyzing the processes responsible for the formation and erosion of the landslide dams of the Seymareh, Jaidar and Balmak (called also Chah Javal) lakes by using available annual sedimentation data and field measurements of the deposits in these lakes.

Despite the number of scenarios proposed so far, quantitative constraints on the river valley evolution before and after the occurrence of this giant landslide are still missing, which could help in quantitatively modeling the trigger scenario of this end-member case study of massive rock slope failures. In this regard, recently, some works (Bozzano et al., 2012, 2016; Della Seta et al., 2017; Martino et al., 2017) have focused on the role of landscape evolution rates on the development of MRC processes. Among these, Bozzano et al. (2016) demonstrated that erosion rates play a key role in the development of MRC processes within the rock masses and, consequently, in their possible evolution into massive rock slope failures, even without invoking transient external forcing (e.g., earthquakes). Furthermore, evidence of giant rock slides due to MRC that caused valley river damming have been recently documented by Zhao et al. (2019) in the Sichuan River basin; these rockslides are related to the effects of river knickpoint propagation and inner gorge formation and serve as a further confirmation of the combined role played by the fluvial dynamics and the geological structural setting.

This work is aimed at filling the knowledge gap in landscape evolution associated with the Seymareh landslide. Detailed geomorphological mapping, correlation and dating of certain geomorphic markers (Burbank and Anderson, 2012) represented by pre- and post-failure fluvial and lacustrine terraces have been performed upstream and downstream of the landslide dam. We provide new time constraints on the Seymareh River valley evolution and outline the role of the geomorphic processes both as predisposing factors for MRC processes and as response to this giant gravitational instability.

2 Regional geological and geomorphological framework

The Seymareh landslide detached from the northeastern flank of the Kabir-kuh fold, the largest and highest anticline in the Pusht-e Kuh Arc, in the northwestern part of Iran (Vergés et al., 2011). The Zagros mountain range is part of the Alpine-Himalayan orogenic system that originates from the Late Cretaceous–Cenozoic convergence between Africa and Arabia–Eurasia (Talbot and Alavi, 1996; Stampfli and Borel, 2002; Golonka, 2004; McQuarrie, 2004; Mouthereau et al., 2012). The Zagros orogen was traditionally classified by distinctive lithological units and structural styles into four NW-trending tectono-metamorphic and magmatic belts (Fig. 1). These are bounded by defects on a regional scale such as the Main Zagros Thrust (MZT), High Zagros Fault (HZF) and Mountain Front Fault (MFF) (Agard et al., 2005, and references therein). These tectonic units, from the inner to the outer sectors of the belt, are the Urumieh Dokhtar volcanic arch, the Sanandaj-Sirjan Zone, the Imbricate Zone, the Zagros (Simply) Folded Belt and the continental Mesopotamian Foreland (Fig. 1).

Figure 1 Simplified structural map of the Zagros mountain range with location of the MZT, Main Zagros Thrust; HZF, High Zagros Fault; MFF, Mountain Front Fault; BR, Bala Rud fault zone (modified from Casciello et al., 2009) . © OpenStreetMap contributors 2019. Distributed under a Creative Commons BY-SA License.

Seismicity is distributed in a 200–300 km wide area of the Zagros mountain range (Hatzfeld et al., 2010; Paul et al., 2010; Rajabi et al., 2011) with a sharp cut along the Main Zagros Reverse Fault in the NE (e.g., Yamini-Fard et al., 2006), with recurrent earthquakes of M_w 5–6 and exceptional earthquakes of higher magnitude, i.e., up to M_w 6–8 (see the Supplement). The Seymareh landslide occurred in a very seismically active area, and the recurrence rate of nearby strongly felt earthquakes considerably higher than the rate of slope steepening led Roberts and Evans (2013) to hypothesize that seismic forcing may have played a primary role in triggering the landslide.

Figure 2 Geological map, stratigraphic column and cross sections of the study area.

The outcropping formations in the Kabir-kuh anticline date to a time interval ranging from the Late Cretaceous to the Early Miocene and are characterized by different lithological and rheological properties (Vergés et al., 2011). Since the geo-structural setting of the fold flanks represented a crucial predisposing factor for the catastrophic massive rock slope failure (Roberts and Evans, 2013), we referred to the most detailed stratigraphic column proposed by James and Wynd (1965), Alavi (2004) and to the detailed mapping of the Kabir-kuh fold conducted by the National Iranian Oil Company (Setudehnia and Perry, 1967; Takin et al., 1970; Macleod, 1970). Specifically, the investigated area includes the middle and low reaches of the Seymareh River starting approximately 60 km upstream of the Seymareh landslide down to the SE termination of the Kabir-kuh fold. In Fig. 2, the geological map of the study area, the stratigraphic column and two geological cross sections related to different structural sectors are reported.

It is noteworthy that, in the Kabir-kuh anticline, the Pabdeh Formation is composed of three rheologically contrasting members, which crop out in the Seymareh landslide scar area: (i) the lower Pabdeh member (150 m thick), which is dominated by marls and shales; (ii) the Taleh Zang member (50 m thick), consisting of platform limestone; and (iii) the upper Pabdeh member (150 m thick), composed mainly of calcareous marl. The Asmari Formation creates a carapace originally covering the top of the Kabir-kuh fold, while in the synclinal valleys between the Kabir-kuh fold and the adjacent folds, the Asmari Formation is overlapped by a Miocene–Pliocene succession (Homke et al., 2004). Referring to the Changuleh syncline studied by Homke et al. (2004), the foreland stratigraphy includes the following: (i) the Gachsaran Formation (Early Miocene–12.3 Ma, thickness approximately 400 m), composed of salt, anhydrite, marl and gypsum; (ii) the Agha Jari Formation (12.3–3 Ma, thickness approximately 1400 m); and (iii) the Bakhtiari Formation (3 Ma–Early Pleistocene, thickness approximately 900 m). The Agha Jari Formation consists of sandstones and conglomerates, linked to the evolution from deltaic to fluvial transitional environments (Elyasi et al., 2014), and the Bakhtiari formation consists of conglomerates characterized by coarse and mud-supported grains, sandstones, shales and silts and marks the onset of syn-orogenic fluvial environment conditions (Shafiei and Dusseault, 2008).

The reported cross sections intersect the synclinal valley of the Seymareh River. The dip angle of the northeastern flank of the syncline considerably decreases from NW to SE from 45^∘ (section A–A') to 18-20^∘ (section B–B'). Therefore, along the section A–A' the Cretaceous-Paleogene bedrock (from the Sarvak Formation to the Asmari Formation) offers a greater accommodation volume to the continental and epicontinental formations (Gachsaran Formation and Agha Jari Formation) as the synclinal axis is located at a lower elevation than in the B–B' section.

The Zagros Range globally provides one of the most spectacular examples of landscape evolution in response to active tectonics (Bourne and Twidale, 2011) because its drainage network clearly adapted to the growth of the thrust-fold structures (Ramsey et al., 2008) and to the erodibility of the outcropping formations (Oberlander, 1985). Oberlander (1968) suggested that the drainage network in the NW Zagros was superimposed from structurally conformable younger horizons. In his model, the breaching of hard geological units of the antiformal ridges follows a phase of river cutting and expansion of the fold axial basins through the softer overlying units. In the Kabir-kuh fold, the transverse cutting of the Asmari limestone, and the exposure of the underlying more erodible Pabdeh-Gurpi marls, leads to the formation of a low-relief landscape with synformal ridges on which the new through-going drainage system can be developed. In Oberlander's hypothesis, it is the Pabdeh and Gurpi marls that facilitate the creation of a low-relief landscape across the anticline crests and are therefore integral to the story of drainage superimposition.

Tucker and Slingerland (1996) computed a numerical landscape evolution model, calibrated on the Kabir-kuh fold, to understand how the growth and propagation of the folds, the different lithologies and the drainage network can influence the sediment flux from a tectonically active belt towards the foreland basin. The authors calibrated the landscape evolution model with the current topography of the range, obtaining time constraints for landscape evolution modeling. According to the Oberlander model, Fig. 3 shows four main steps that describe the landscape evolution of the Kabir-kuh fold with the timing provided in the model by Tucker and Slingerland (1996).

• Step 1 – Approximately 4.3 Ma, in response to the initial stages of fold growth, an orthoclinal drainage develops parallel to the main structures. The tributaries flowing along the flanks of the folds transport debris, which is deposited in the synclines. In the Kabir-kuh fold the carbonate core is still buried by the Miocene cover units.

• Step 2 – Approximately 3.8 Ma, as soon as the deformation front migrates towards the SW, new folds raise with a progressive adjustment of the drainage to these morpho-structures. The previously deposited sediments are remobilized and transported towards the depocenter of the syncline basins and partly outside; the syn-orogenic deposits are strongly eroded along the crests of the anticlines, thus exposing the underlying formations. This causes a topography characterized by resistant hogbacks that flank the inner cores.

• Step 3 – Approximately 2.4 Ma, with the ongoing deformation, the drainage develops in a “trellis” pattern. The river erosion affects the erodible units located stratigraphically between the limestone of the Asmari Formation and the inner core of the fold. At the end of this step the Miocene cover is completely removed from the ridges and the river erosion also affects the marls and evaporites of the syn-orogenic formations in the valleys, exposing the underlying limestone of the Asmari Formation.

• Step 4 – Approximately 1.6 Ma, due to the continuous uplift and exhumation of younger more external folds, the sediment accumulation becomes negligible and the Asmari limestone is strongly eroded giving rise to syncline ridges. The following Quaternary landscape evolution is then likely driven by the evolution of the drainage network and is also influenced by climatic factors and by the slope-to-channel dynamics.

The model by Tucker and Slingerland (1996) is the unique numerical model existing on the Kabir-kuh fold and this motivates our choice of using it as a reference for the medium- to long-term evolution of the Seymareh River valley. The Seymareh River valley is arranged parallel to the Kabir-kuh fold and its evolution was inevitably influenced by the exceptional landslide event that temporarily dammed it, causing the formation of the three-lake system which included the Seymareh, Jaidar and Balmak lakes (Fig. 4). The valley evolution before and after the event is well recorded by Quaternary landforms preserved along the valley. Yamani et al. (2012) focused on the post-failure evolution of the valley describing four levels of terraces upstream of the landslide dams as a sequence of lacustrine terraces. Shoaei (2014), in addition to evaluating the longevity of the Seymareh landslide dams, identified the merging of Seymareh River with the Kashgan River left tributary as the reason for strong river incision at the base of the northeastern flank of the Kabir-kuh fold and as a possible causal factor for the Seymareh landslide collapse.

Figure 3 Evolution of the drainage network in the Zagros chain sector of the Kabir-kuh fold, according to the model by Oberlander (1985) and with the timing provided in the landscape evolution model by Tucker and Slingerland (1996). See the text for explanation of the four steps (modified from Oberlander, 1985).

Figure 4 Overview and zoomed-in view (in the red box) of the Seymareh landslide, Zagros Fold–Thrust Belt, NW Iran (modified from © Google Earth).

However, none of the previous studies on the Quaternary evolution of the Seymareh River valley provided absolute dating of geomorphic markers (mainly fluvial terraces), preserved upstream as well as downstream of the landslide dam, or provided robust and quantitative constraints on the pre-failure valley evolution as a possible geomorphological factor for failure occurrence.

3 Methods

The geomorphological study of the area was carried out first through the analysis and interpretation of remote sensing data, such as aerial photos (National Cartographic Center of Iran, aerial photo, scale: 1 : 20 000, acquired on 24 August 1955), Google Earth satellite optical images (2018 Landsat Imagery) and vector topographic maps (National Cartographic Center of Iran, topographic map of Kuhdasht, scale: 1 : 25 000), which led to the first detection of possible geomorphic markers within the Seymareh River valley. Vector topographic data also allowed the construction of a 10 m digital elevation model (DEM) for terrain analyses and led to the projection of the possible geomorphic markers along the river longitudinal profile (Wilson and Gallant, 2000; Burbank and Anderson, 2012). The DEM was obtained by the ArcGIS 10^® software package, starting from vector topographic data (contour lines, hydrography and point elevation) and using the ANUDEM interpolation algorithm (Hutchinson et al., 2011, and references therein). To automatically extract the hydrographic network from the DEM and then to project the geomorphic markers along the longitudinal river profiles, some of the ArcGIS^® 10 tools of the “Hydrology toolbox” were used (Flow Direction, Flow Accumulation, Reclassify, Stream Order and Stream to Feature), setting the flow accumulation threshold according to that proposed for the fluvial domain (10⁻¹ km²) by Montgomery and Foufoula-Georgiu (1993). The longitudinal profile was therefore transformed into a route along which the elevation of the top surfaces of geomorphic markers identified in the area were projected through the linear referencing tools (Create Route and Locate Features along Route).

A geological and geomorphological field survey was then carried out with the aims of mapping the most significant active and relict landforms for the Quaternary evolution of the Seymareh River valley and of sampling the corresponding deposits in order to date them with the OSL method (optically stimulated luminescence; Murray and Olley, 2002; Wintle and Murray, 2006; and references therein).

OSL sampling is a very delicate and quite complex technique. In fact, it is absolutely necessary to prevent the sample from being exposed to light because the luminescence signal could be reduced or even reset. In choosing the most suitable site to sample, of course, levels were identified with original sedimentary structures, avoiding bioturbations and post-depositional alterations. Once the site for sampling was identified, it was important to carefully clean off the slope and prepare, according to the consistency or cementation of the material, the equipment necessary for taking the sample, without it being exposed to light. Furthermore, to minimize the effects of cosmic radiation and to thereby avoid the risk of rejuvenated ages, the samples were taken at least one meter below the topographic surface (or below eventual erosional surfaces identified within the deposits). All of the samples, mainly characterized by fine-grained loose sediments (size <2 mm) were taken directly by using a hammer to insert a metal tube horizontally into a vertical face, which must be isolated from light and humidity immediately after collection. To maximize the uniformity of the natural radioactivity of the burial period, the tube was inserted into zones of homogeneous sediment at least 30 cm wide and thick. From the same level where it was sampled, an additional 500–800 g of sediment was extracted to evaluate natural radioactivity (if the annual dose rate measurement is not performed in situ), for the mineralogical and granulometric analyses, as well as to determine the moisture content.

The OSL dating was performed by the LABER OSL Laboratory, in Waterville, Ohio, USA. Quartz was extracted for equivalent dose (D_e) measurements. In the OSL laboratory, the sample was treated first with 10 % HCl and 30 % H₂O₂ to remove organic materials and carbonates, respectively. After grain-size separation, the fraction of 90–125 µm size was relatively abundant, so this fraction was chosen for D_e determination. The grains were treated with HF acid (40 %) for approximately 40 min to remove the alpha-dosed surface, followed by 10 % HCl acid to remove fluoride precipitates. Luminescence measurements were performed using an automated Risø TL/OSL-20 reader. Stimulation was carried out by a blue LED ($\mathit{\lambda }=\mathrm{470}\pm \mathrm{20}$ nm) stimulation source for 40 s at 130 ^∘C. Irradiation was carried out using a 90Sr/90Y beta source built into the reader. The OSL signal was detected by a 9235QA photomultiplier tube through a U-340 filter with 7.5 mm thickness. For D_e determination, the single-aliquot regenerative-dose (SAR) protocol (Murray and Olley, 2002; Wintle and Murray, 2006) was adopted. The preheating temperature was chosen to be 260 ^∘C for 10 s and the cut heat was 220 ^∘C for 10 s. The final D_e is the average of the D_e of all aliquots, and the final D_e error is the standard error of the D_e distribution. For each sample, at least 12 aliquots were measured for D_e determination. The D_e was measured using SAR on quartz, and the aliquots that passed criteria checks were used for final D_e calculation.

Recycling ratios were between 0.90–1.1, and recuperation was relatively small. The cosmic ray dose rate was estimated for each sample as a function of depth, altitude and geomagnetic latitude. The concentration of U, Th and K was measured by neutral activation analysis (NAA). The elemental concentrations were then converted into the annual dose rate, taking into account the water content (lab measured) effect. The final OSL age is then D_e/Dose rate.

4 Results

The best geomorphic markers preserved in the study area are represented by a lacustrine terrace and two suites of fluvial terraces. The lacustrine terrace and the fluvial terrace suite upstream testify to the evolution after the landslide emplacement. The fluvial terrace suite downstream testifies to the evolution before the landslide emplacement. Fluvial terraces consist of gravel, sand, silt and clay, while the lacustrine deposits are mainly composed of silt (detailed information on each sampled deposit is provided in the Supplement). Conglomerates outcropping immediately upstream and downstream of the landslide pertain to inactive alluvial fans connected to a relict position of the valley floor, likely of a paleo-Seymareh river.

Table 1 OSL ages obtained for the geomorphic markers recognized in the Seymareh River valley. Detailed information and photos of sampling sites are available in the supplementary material.

^* Minimum ages from samples saturated due to their low concentration of quartz grains.

In the upper reach of the Seymareh River valley, the geomorphic markers include inactive, terraced conglomeratic alluvial fans (Cg_m), a terraced lacustrine deposit and a suite of four orders of fill terraces (named from Qt1_m to Qt4_m). The geometry of the terraced conglomerates (section A-A' in Fig. 2) can be associated with alluvial fans generated on the flanks of a former synclinal valley by streams likely forming the tributaries of a paleo-Seymareh river whose path was to the SW of the present one. The fill terraces are entrenched in the terraced lacustrine deposit of Seymareh Lake upstream of the landslide, in the area where Harrison and Falcon (1938), Roberts and Evans (2013) and Shoaei (2014) hypothesized this natural damming lake could be located (Figs. 5 and 6). Prograding lacustrine fan deltas formed by tributaries of Seymareh Lake have been recognized (Fig. 5). In this sector, we successfully dated four samples (SEY4, SEY5, SEY6 and SEY8; Table 1 and the Supplement); in particular, the lacustrine deposit at two different stratigraphic levels, 560 and 590 m a.s.l., which provided OSL ages of 10.4±0.90 ka (sample SEY8) and 7.37±0.73 ka (sample SEY4), respectively. The OSL age of 17.9±1.50 ka (SEY6) obtained for an alluvial deposit at the base of the lacustrine deposits predated by ca. 7 kyr the emplacement of the Seymareh landslide, according also to the time constraints provided by Roberts and Evans (2013).

Figure 5 Geomorphic markers upstream of the Seymareh landslide, represented by a suite of four orders of alluvial terraces entrenched in the lacustrine deposits (Lac) of Seymareh Lake upstream of the landslide, in the areas where Harrison and Falcon (1938), Roberts and Evans (2013), and Shoaei (2014) hypothesized the natural damming lake could be extended. (a) Overall view of the suite of terraces; (b) example of fluvial terrace deposit; (c) example of lacustrine deposit; (d) evidence of a prograding lacustrine fan delta formed by one of the right tributaries of Seymareh Lake during its early emptying phase.

Figure 6 Map of the lacustrine and alluvial terrace suite and of the most significant landforms for the valley slope evolution upstream of the Seymareh landslide.

A suite of two strath terraces and a flood plain formed onto the landslide debris along the Seymareh River gorge have been identified (Fig. 7a and b), which are important markers of the evolution of the natural dam because they formed after its cut likely due to an overflow of the damming lake. Here, we successfully dated one sample taken on a strath terrace (SEY9; Table 1 and the Supplement), which provided an age of 6.59±0.49 ka as time constrain for the initial stage of lake emptying.

Figure 7 Geomorphic markers upstream of the Seymareh landslide. (a) Strath terraces and a flood plain developed over the landslide debris, which are important markers of the evolution of the natural dam since they testify to the overcoming of the damming lake and the overflooding of the river onto the landslide debris, respectively; (b) detail of the strath terrace deposit sampled for OSL dating; (c) bedrock exposure below the landslide debris along the modern river profile passing through the relict landslide dam; (d) the suite of fluvial terraces downstream of the Seymareh landslide; the Qt1 level is poorly preserved and not visible in this photo.

Figure 8 Map of the alluvial terrace suite and of the most significant landforms for the valley slope evolution downstream of the Seymareh landslide.

In the lower reach of the Seymareh River valley, the downstream geomorphic markers include inactive, terraced conglomeratic alluvial fans (Cg_l) and a suite of four orders of fill terraces (named from Qt1_l to Qt4_l) downstream of the Seymareh landslide (Figs. 7c and 8). Here, we successfully dated three samples from the fill terrace deposits (SEY3, SEY10, SEY11; Table 1 and the Supplement). The ages obtained provide useful time constraints on the main depositional events during the pre-failure valley evolution. Minimum ages of 373±34 and 312±45 ka have been obtained for samples SEY3 and SEY10, respectively, since these samples were saturated due to their low concentration of quartz grains, and SEY11 was dated at 60±5 ka.

The above-described geomorphic markers of the Seymareh River valley have been mapped and reported in morpho-stratigraphic profiles. The most significant landforms for the valley slope evolution are presented with a detail for the post-failure fluvial and lacustrine terrace suites upstream of the landslide dam (Fig. 6) and the pre-failure fluvial terrace suite downstream of the landslide dam (Fig. 8), respectively.

Figure 9 Projection of the pre-failure geomorphic markers along the longitudinal profile of Seymareh River. They are named according to the legend of Figs. 6 and 8 (see the text for explanation). The obtained OSL ages are indicated.

Figure 10 Projection along the longitudinal profile of Seymareh River of the post-failure geomorphic markers They are named according to the legend of Fig. 6 (see the text for explanation). The obtained OSL ages are indicated.

Figures 9 and 10 report the longitudinal profile of Seymareh River, along which, in addition to the geomorphic markers, the following were projected also: the benchmarks of the basal contact of the Quaternary deposits on the bedrock, the projection of points corresponding to the top of the Seymareh landslide debris, the upstream and downstream limits of the landslide, the location of the OSL sampling, and the projection of the outcrop of the Bakhtiari Formation (Fig. 9), which is rarely preserved and marks the initial alluvial infill of the Seymareh River valley. Figure 9 shows the height distribution of the pre-failure geomorphic markers. The benchmarks along Seymareh River indicate a mostly bedrock channel, and the longitudinal profile is characterized by two knickpoints located upstream of the Seymareh landslide and downstream of the lowest suite of alluvial terraces (as indicated by the black arrows in Figs. 9 and 10). The geomorphic markers downstream and upstream do not belong to the same suite of terraces as their projections along Seymareh River do not have any topographic correlation to each other (Figs. 9 and 10). The tops of all the fluvial terraces downstream of the Seymareh landslide are located lower in height than the most important knickpoint located immediately upstream and sculpted in the bedrock. Figure 10 shows the height distribution of the post-failure geomorphic markers. The markers are represented by (i) the horizontal lacustrine terrace formed by the incision of the deposits pertaining to the Seymareh Lake, formed as a consequence of the landslide damming; (ii) the two levels of the strath terraces and a flood plain formed on the landslide debris during the initial stages of dam cutting and emptying of the lake; and (iii) the four fill terrace levels formed after the emptying of the Seymareh Lake.

Figure 11 Scar area of the Seymareh landslide. (a) Front view of the scar area with the location of sites where evidence of buckling has been recognized (Source © Google Earth: Image © 2018 CNES/Airbus, Image © 2018 DigitalGlobe, Image © Landsat/Copernicus); (b) ductile buckling deformations within the layers of the upper Pabdeh member which cannot be referred to as parasitic structural folding due to their localization just below the Seymareh landslide sliding surface and to their reduced persistency (both lateral and vertical) within the rock mass; (c) brittle buckling deformation of the Taleh Zang member along the scar area.

The geomorphological field survey, supported by a remote survey based on optical satellite and aerial images, also allowed us to demonstrate the evidence of gravity-induced features of the downslope-dipping strata, along the scar of the Seymareh landslide. As shown in Fig. 11, evidence of MRC driving towards stress concentration and failure has been recognized in gravity-induced folding within the thin-layered Pabdeh Formation just below the sliding surface of the Seymareh landslide (i.e., that cannot be ascribed to parasitic structural folding). Furthermore, impressive buckling of the downslope-dipping strata, which crop on the sliding surface of the Seymareh landslide, have been interpreted as a release of concentrated stresses due to the post-failure rebound caused by the collapsed rock mass.

5 Discussion

5.1 Constraints on pre-failure valley evolution

The longitudinal profile of the Seymareh River and the geomorphic markers preserved mainly downstream of the landslide dam provided new constraints on the pre-failure valley evolution. The major knickpoint located immediately upstream of the Seymareh landslide is the most interesting to be analyzed in relation to the landslide event (Fig. 9). Its shape in the longitudinal profile clearly lets us identify it as a “slope-break knickpoint” (Kirby and Whipple, 2012; Boulton et al., 2014), thus developed as a knickpoint retreating in response to a persistent perturbation to the fluvial system (Tucker and Whipple, 2002), as frequently observed in tectonically active regions. The location of this knickpoint upstream of the Seymareh landslide and the exposure of the basal contact of the landslide at the bottom of the Seymareh River gorge (Fig. 7c) suggest that this shape of the longitudinal profile was already developed before the failure, meaning that the erosion wave which generated the knickpoint affected the Seymareh landslide slope foot before the failure occurrence.

The poorly preserved well-cemented alluvial fan conglomeratic deposits outcropping upstream of the landslide lie on the Miocene Agha Jari Formation at a higher elevation than the outcrops of the Bakhtiari Formation. Their remnants are aligned in correspondence with the axis of a relict synclinal valley, likely corresponding to a very early stage (possibly Pliocene) of the Seymareh valley evolution. On the other hand, the conglomerate deposits outcropping downstream of the landslide (Cg_l) are closer in height to the major knickpoint, thus suggesting that they were in equilibrium with a local base level corresponding to the early propagation of the major knickpoint. Furthermore, they must be younger than the Bakhtiari Formation, which is preserved at higher elevation.

The alluvial terraces preserved downstream of the Seymareh landslide likely mark the valley evolutionary stages during the major knickpoint retreat (Demoulin et al., 2017). Along the longitudinal river profile, the uppermost outcrops of each level of this terrace suite were swept away by the landslide, which unfortunately prevents estimation of the rates of knickpoint retreat. Nonetheless, according to what was observed by Bridgland et al. (2017) about river terrace development in the NE Mediterranean region, the sedimentation phases should correspond to cold periods. In particular, Bridgland et al. (2012) observed – in the valleys of the Tigris and Ceyhan in Turkey, the Kebir in Syria and the transborder rivers Orontes and Euphrates – a regular terrace formation in synchrony with 100 ka climatic cycles that can be correlated with MIS 12, 10, 8, 6 and 4–2 (marine isotope stages). Therefore, the minimum ages obtained for the SEY3 and SEY10 samples could be reasonably extended to 478 ka (MIS 12) and 374 ka (MIS 10), respectively, and the OSL age of the SEY11 fits well with the Last Glacial Period.

5.2 Constraints on post-failure valley evolution

The geomorphic markers preserved upstream of the landslide dam provided new constraints on the geomorphic response of the Seymareh valley to the 44 Gm³ natural dam (Roberts and Evans, 2013). Such a response was first the formation of three lakes (Seymareh, Jaidar and Balmak; Fig. 4) whose persistence and evolution is well recorded by the deposits outcropping in the valley. In this regard, the estimation of a sedimentation rate of 10 mm yr⁻¹ in the Seymareh Lake was obtained using the OSL ages of 10.4±0.90 and 7.37±0.73 ka for the lacustrine deposit sampled at 560 and 590 m a.s.l., respectively.

The strath terrace sculpted on the landslide deposit and dated at 6.59±0.49 ka is the oldest time-constrained marker of the natural dam breaching due to overflow, which caused the lake to empty. The lake overflow was likely caused by the gradual filling of the reservoir with lacustrine deposits, which progressively reduced the dam infiltration section. Many studies treated the effects of infiltration on the longevity of large landslide dams (Ischuk, 2011; Schuster and Alford, 2004; and references therein). Internal drainage through the Seymareh landslide dam – or at least parts of it – is extremely likely given the coarse nature of the debris (Watson and Wright, 1969; Roberts and Evans, 2013). Nevertheless, the possible role of groundwater seepage within the pervious natural dam in balancing the Seymareh River discharge and delaying the dam overflow remains a questionable topic to be approached and solved in future studies. Despite their interpretation as progressively younger lacustrine deposits by Yamani et al. (2012), the four terrace levels entrenched in the terraced lacustrine deposit show a longitudinal downstream gradient, which, along with their sedimentological characters, identify them as fluvial fill terraces. Furthermore, the OSL age obtained for the lacustrine deposit at the base of the Qt2_m terrace (sample SEY8) is 10.4±0.90 ka, demonstrating that the suite of alluvial terraces is entrenched into the same (and unique) lacustrine deposit. The OSL age of 4.49±0.48 ka obtained for the Qt1_m terrace (sample SEY5) provides time constraints on the emptying phase of Seymareh Lake. Such time constraints are fine-tuned by the age of the strath terrace formed on the landslide debris, which corroborate the initial stage of lake emptying at 6.59±0.49 ka (SEY9). As indicated by the age of the Qt1_m terrace (of 4.49±0.48 ka), the Seymareh Lake likely persisted up to ∼5 ka, much longer than the 935 years estimated by Shoaei (2014). Since the top of the lacustrine deposit lies at 630 m a.s.l., an increased sedimentation rate of ∼17 mm yr⁻¹ can be inferred for the late stage of the lake evolution, which is in agreement with an increased sediment yield from tributaries during the early stages of lake emptying (Fig. 5d).

The oldest strath terrace sculpted onto the landslide deposit and dated at 6.59±0.49 ka is just a few meters below its top. Therefore, it can be reasonably used to calculate the erosion rate affecting the landslide deposit after the overflow. The ratio between the thickness of the eroded sediment below the strath surface (∼120 m) and the time elapsed since the beginning of the process (∼6.59 kyr) allows for estimation of an erosion rate of 1.8 cm yr⁻¹ for Seymareh River along the gorge. The cut of the landslide dam induced a new change in the fluvial base level, bringing the slope-to-valley floor system into disequilibrium. For this reason, a dense drainage system was set on the scar area, which, due to the high erodibility and low permeability of the less competent Pabdeh-Gurpi Formation immediately below the sliding surface on the Kabir-kuh ridge NE slope, has generated the badlands mapped in Fig. 6.

5.3 Evolutionary model of the Seymareh River valley

The landscape evolution of the Seymareh River valley before and after the failure occurrence can be summarized in the following six phases:

1. Setting of a paleo-Seymareh river into a synclinal valley, likely developed in the Pliocene, to the west of the present position of the Seymareh River and deposition of fan deposits (Cg_m) (Fig. 12a).

2. Development of the valley with local base level correlated to the Seymareh longitudinal profile segment upstream of the major knickpoint along the Seymareh River and coeval to the deposition of the Bakhtiari Formation (Late Pliocene–Early Pleistocene) (Fig. 12b).

3. Emplacement of the downstream fan deposits corresponding to the Cg_l conglomerates (Early Pleistocene) and generation of the four orders of Middle–Late Pleistocene alluvial terraces (Qt1_l–Qt4_l) preserved downstream of the landslide and formed during the progressive migration of the major knickpoint, which is presently located upstream of the landslide (Fig. 12c).

4. Seymareh landslide event (∼10 ka), according to the ¹⁴C ages by Roberts and Evans (2013) and to the OSL ages provided in this work for the lacustrine deposits (Lac) (Fig. 12d).

5. Formation and permanence of the Seymareh Lake (∼ 10–6.6 ka), according to the ¹⁴C estimated ages by Roberts and Evans (2013) and to the OSL ages provided in this study for the lacustrine deposits (Lac) (Fig. 12e). The progressive infilling of the lake reservoir progressively reduced the infiltration section on the upstream side of the landslide dam. The presence of a minor emissary on the downstream side of the landslide debris cannot be excluded.

6. Overflow of the lake and cut of the natural dam with formation of the first strath terrace (6.59±0.49 ka), followed by a second strath terrace and a flood plain during the emptying of the lake, which upstream is associated with the sedimentation of a fluvio-lacustrine sequence at the top of the lacustrine sediments (Fig. 12f).

7. Complete emptying of the lake and generation of the suite of fill terraces entrenched in the deposits of Seymareh Lake (4.5 ka–present) (Fig. 12g).

Figure 12 Evolutionary model of the Seymareh River valley. See text for explanation. Traces and legend of geological cross sections are reported in Fig. 2.

5.4 Implications of the evolutionary model for future back-analysis of the Seymareh landslide

According to the multimodeling approach proposed by Martino et al. (2017), Quaternary landscape evolution modeling of slope-to-valley floor systems plays a key role as a tool for chronological constraints on the creep evolution of entire slopes (Bozzano et al., 2016; Della Seta et al., 2017).

The geomorphic processes developed before the failure of the Seymareh landslide likely acted as predisposing factors for MRC processes as the rock mass successively collapsed. Kinematic freedom, both at the top and on the fold flank, was created by the incising network of streams that dissect the Asmari Formation carbonate caprock following the major joint set in the Asmari Formation already described in Roberts and Evans (2013, and references therein). In particular, the headward erosion of streams towards the anticline's structural high, described by Oberlander (1968), caused the expansion of the fold axial basins through the softer units, determining the upslope kinematic freedom. In the timing proposed by Tucker and Slingerland (1996), the latter was reached at approximately 1.6 Ma. Stress release at the slope base was definitely produced by the Middle–Late Pleistocene upstream migration of the knickpoint along the Seymareh River longitudinal profile. Unfortunately, since the emplacement of the landslide swept away the uppermost outcrops of the alluvial terraces formed in response to the upstream knickpoint migration, the rate of knickpoint migration cannot be inferred. Nonetheless, an elapsed time-to-failure value on the order of 10² kyr, since the kinematic freedom at the slope base was reached, can be reasonably estimated by the age of the oldest terrace in the lower reach of the river minus the age of the landslide occurrence.

It is noteworthy that the stratigraphy of the source rock mass, also described in detail by Roberts and Evans (2013), accounts for different rheological behaviors, which could have induced differential strain rates within the slope leading to failure according to a MRC process (Fig. 11). To date, a quantification of such rheological properties is lacking for the lithological units of the Kabir-kuh fold. Nonetheless, some inferences can be done according to previous studies on rock masses affected by MRC (Apuani et al., 2007; Bozzano et al., 2012; Bretschneider et al., 2013; Della Seta et al., 2017). More particularly, the time-dependent viscoplastic behavior, more typical of clayey and marly deposits, which have lower viscosity values, can justify time-dependent (creep) strains which could have generated a high-stress concentration within the higher viscosity level over time (i.e., mostly characterized by elastoplastic rheology), inducing their cracking and leading to failure. In fact, a stiffness contrast exists between the upper member of the Pabdeh Formation and the overlying Asmari Formation. The attitude of the strata is moderately dipping downslope (15–20^∘), leading to a lower vertical thickness of (and consequently a reduced lateral confining effect by) the epicontinental and continental units, such as the Gachsaran and Agha Jari Formation. Moreover, the low dip angle of the strata reduces also the vertical thickness of the Asmari Formation caprock (Roberts, 2008) and consequently the incision necessary for kinematic freedom of the landslide, relative to a more steeply dipping sequence, which would require a greater depth of incision.

Therefore, the results of this work have implications for a future back-analysis through stress–strain numerical modeling of the Seymareh landslide because they can be used to constrain the elapsed time since MRC initiation and ultimate failure conditions. Such a perspective is to be regarded as a key challenge for dimensioning such an end-member event in regards to the spatiotemporal distribution as well as for evaluating the possible role of impulsive triggering actions (i.e., strong to very strong earthquakes) in anticipating the time-to-failure value of the slope.

6 Conclusion

In a multimodeling approach to the study of MRC processes affecting slopes at a large spatiotemporal scale, the performed geomorphic analysis allowed us to constrain the evolution of the Seymareh River valley in the northwestern Zagros Mts., before and after the failure of the largest landslide ever recorded on the exposed Earth surface. The identification and OSL dating of different suites of lacustrine, alluvial and strath terraces constrained the major pre- and post-failure evolutionary steps of the river valley system in time.

The oldest geomorphic markers in the Seymareh River valley are represented by relict conglomerates preserved upstream of the landslide, which demonstrate the early (possibly Pliocene) position of a paleo-Seymareh river flowing into a synclinal valley close to the northeastern flank of the Kabir-kuh fold.

Drainage evolution associated with the growth of the Kabir-kuh fold was characterized by the deep incision of the stream network, which allowed the kinematic release of the rock mass involved in the giant Seymareh landslide. Such a stream incision was accompanied by the retreat of a major “slope-break knickpoint” along the Seymareh longitudinal profile, time constrained by the age of a suite of river fill terraces. According to the age of pre-failure terraces, in the Middle–Late Pleistocene, the erosion wave reached the portion of the Kabir-kuh fold that was affected by the Seymareh landslide at ∼10 ka. According to the timing of the landscape evolution model proposed by Tucker and Slingerland (1996), the upper slope underwent kinematic release at about 1.6 Ma. Therefore, the collapse was prepared by MRC processes acting over a time window of 10² kyr.

The geomorphic response to the landslide dam consisted of the formation of three lakes, among which Seymareh Lake persisted for ∼3500 years before its emptying phase started ∼6.6 ka due to lake overflow. A sedimentation rate of 10 mm yr⁻¹ was estimated for the lacustrine deposits, which increased up to 17 mm yr⁻¹ during the early stage of lake emptying due to the increased sediment yield from the lake tributaries. Since ∼4.5 ka, a suite of four alluvial terraces upstream of the landslide demonstrates the alternating erosion and deposition phases of the re-established Seymareh River.

An incision rate of 1.8 cm yr⁻¹ was estimated since the beginning of the landslide cut by Seymareh River, and such a strong erosion started propagating up to the landslide source area where badlands developed, eroding the marly Pabdeh-Gurpi Formation.

The results obtained here provide new constraints on the valley evolution in view of future stress–strain numerical modeling of the MRC process that involved the Seymareh landslide slope before its generalized collapse. Such a modeling could also be considered to discuss the possible role of impulsive triggering (earthquakes) in anticipating the time-to-failure value due to the gravity-driven deformational processes.