Active Deformation and Relief Evolution in the Western Lurestan Region of the Zagros Mountain Belt: New Insights From Tectonic Geomorphology Analysis and Finite Element Modeling

To unravel how and where coseismic and interseismic deformation impacts the spatial and temporal patterns of rock uplift of the Lurestan sector of the Zagros Mountains, we performed an investigation of the large‐scale features of topography and river network coupled with 2‐D finite element modeling. Geomorphological analysis and constraints from parameters such as elevation, local relief, normalized channel steepness index (ksn), river longitudinal profiles, and transformed river profiles (chi plots) were used to unravel the time‐space distribution of vertical motions. Whereas the much longer timescale over which topography grows and/or rivers respond to tectonic or climatic perturbations with respect to even multiple seismic cycles, the outputs of the finite element model yield fundamental information on the source of the late part of the spatiotemporal evolution of surface uplift recorded by the geomorphology. Model outputs shed new light into the processes controlling relief evolution in an actively growing mountain belt underlain by a major blind thrust. The outputs illustrate how coseismic slip controls localized uplift of a prominent topographic feature—the Mountain Front Flexure—located above the main upper crustal ramp of the principal basement thrust fault of the region, while continuous displacement along the deeper, aseismic portion of the same basement fault controls generalized uplift of the whole crustal block located farther to the NE, in the interior of the orogen.


Introduction
The Zagros mountain belt ( Figure 1) represents one of the most active collisional orogens on Earth (Agard et al., 2005;Alavi, 1994;Berberian & King, 1981;Dercourt et al., 1986;Koshnaw et al., 2017;Stampfli & Borel, 2002;Stocklin, 1968). Yet, active tectonics studies and particularly the investigation of large-scale surface motions in the region are challenging due to both (1) the varying erodibility of outcropping rocks, which affect the distribution of topographic high and lows and (2) the almost complete absence of alluvial deposits in the erosion-dominated mountain belt, which complicates reconstructions of the spatial distribution of the highs and lows in the past. Those two issues hinder the time-space reconstruction of vertical movements. The relief of the Zagros mountain belt is strongly controlled by variable resistance of the outcropping rocks to erosion (e.g., Oberlander, 1968Oberlander, , 1985Ramsey et al., 2008;Tucker & Slingerland, 1996).
Features such as anticlinal domes and elliptical hogbacks encircling breached anticlines, both occurring in hard (e.g., Oligocene Asmari Fm.) carbonate rocks, are common in areas where erosion has stripped the stratigraphically higher, weak rocks including Miocene shales and evaporites (e.g., Oberlander, 1965;Ramsey et al., 2008;Tucker & Slingerland, 1996). The Zagros sector known as the Simply Folded Belt is characterized by the widespread occurrence of anticlinal ridges formed in hard rocks ("whaleback anticlines" sensu Ramsey et al., 2008, formed in the Asmari carbonates) onto which erosion landforms such as fossil fluvial paths are preserved. This pattern of outcrops has allowed extracting much information on the lateral development and linkage of individual folds, and inferring information on relative fold age and fold uplift rate (e.g., Bretis et al., 2011;Burberry et al., 2007Burberry et al., , 2010Collignon et al., 2016;Ramsey et al., 2008;Zebari et al., 2019). However, the noise represented by the rugged topography sculpted into alternating soft and hard rocks complicates the identification of geomorphological signatures of both differential and large-scale motions. Within such a setting, a study of the factors that modulate mountainous topography benefits from combined qualitative and quantitative constraints provided by the analysis of topography and river networks, both of which preserve features created in response to local surface changes and/or regional long-term external processes (Dehbozorgi et al., 2010;Forte et al., 2014;Snyder et al., 2000;Whipple, 2004;Whipple & Tucker, 1999;Whittaker, 2012;Willett et al., 2001;Wobus et al., 2006). Such tectonic geomorphology analyses are best integrated with numerical modeling in order to provide useful insights into the tectonic processes controlling active deformation. The value of numerical modeling is particularly important in regions affected by crustal thrust earthquakes that do not rupture the topographic surface (Tavani, Parente, Puzone, et al., 2018;Vajedian et al., 2018), potentially resulting in complex surface deformation patterns (e.g., Ellis & Densmore, 2006;Picotti & Pazzaglia, 2008). In this work, a 2-D elastic finite element model (FEM) of a crustal geological section is discretized in order to simulate interseismic stress and strain accumulation and to obtain information on the vertical surface displacement associated with both interseismic and coseismic stages in a region characterized by thrust-related seismicity. We use the large earthquake of 12 November 2017 (Mw = 7.3), which occurred at the junction between the Lurestan arc (salient) and the Kirkuk embayment (recess) of the Zagros fold and thrust belt (Figure 2), to constrain coseismic slip in the model. This seismic event shows a clear connection with the main thrust fault of the region (Tavani, Parente, Puzone, et al., 2018), which we term here the Main Frontal Thrust (MFT). This thrust underlies at depth the Mountain Front Flexure (MFF), a major morphotectonic feature that, according to Tavani et al. (2020), developed since circa 3 Ma by basement-involved thrusting as deformation migrated downward at deeper crustal levels. Late-stage, crustal ramp-dominated thrusting (Butler & Mazzoli, 2006) would have led to the development of this prominent geomorphological boundary between the high Zagros mountains and the low foothills to the SW (Berberian, 1995;Emami et al., 2010;Falcon, 1961;Sepehr & Cosgrove, 2004; Figure 2).
Modeling of the MFT activity is carried out using a recent, cutting edge FEM methodology that has been successfully applied and validated by various studies in different regions of the world (Candela et al., 2015;Carminati & Vadacca, 2010;Megna et al., 2005Megna et al., , 2008Vigny et al., 2009;Zhao et al., 2004;Zhu & Zhang, 2013). This method analyses the interseismic deformation characterizing zones of plate convergence as a result of the elastic strain affecting a brittle layer resting above a ductile half-space (Savage, 1983). Within

10.1029/2020TC006402
Tectonics this framework, the brittle layer can be assumed as encompassing the entire thickness of the model characterized by stick-slip behavior. Using the same numerical method, the coseismic deformation associated with a specific seismic event can also be modeled. Constraints to the model include available GPS data and geological information derived from both published studies and our own fieldwork in the study area. The outcomes of finite element modeling are then compared with information on the space-time distribution of surface vertical motions inferred from the tectonic geomorphology analysis. In this way, we use a topographic analysis in combination with lithologic boundary conditions to study the locus of recent uplift and then combine it with an FEM to explore the crustal processes (specifically interseismic and coseismic deformation) that could explain the pattern of uplift. The results of this work allow us to obtain a comprehensive model of the evolution of topography in response to the vertical component of surface displacement along a section (M-M′ in Figure 2) across the western Lurestan region of the Zagros mountain chain. Besides increasing our understanding of the somewhat elusive active tectonic behavior of the Zagros fold and thrust belt, our results provide new, general insights into the processes controlling relief evolution in areas affected by large, crustal thrust earthquakes.

Geological Background
The NW-SE striking Zagros mountain belt formed due to the convergence between the Arabian and Eurasian plates since the Late Cretaceous (Stampfli & Borel, 2002;Stocklin, 1968;Talbot & Alavi, 1996). Following the subduction of Neo-Tethyan oceanic lithosphere beneath Eurasia, continent-continent collision occurred during the Early Miocene (Agard et al., 2011;Csontos et al., 2012;Koshnaw et al., 2017;Mouthereau et al., 2012;Vergés et al., 2011). GPS measurements show that the northward relative motion of the Arabian Plate, oblique with respect to the NW-SE trend of the Zagros belt, is still active today and occurs at~2 cm/year with respect to fixed Eurasia (Vernant et al., 2004). It is generally believed that a major fault (Main Recent Fault, MRF), parallel to the strike of the thrust belt, accommodated the entire strike-slip component of the Zagros orogeny during Cenozoic continent-continent collision (Allen et al., 2004;Authemayou et al., 2006;Berberian, 1995;Blanc et al., 2003;Talebian & Jackson, 2002). The Zagros Fold and Thrust belt is bounded to the SW by the MFF and to the NE by the MRF, which separates the Arabian Plate from the Sanandaj-Sirjan Zone (Iran block; Berberian, 1995). The MFT represents the main thrust in the outer part of the orogen, where it exerts a major control on relief evolution (Berberian, 1995;Emami et al., 2010;Koshnaw et al., 2017;Tavani, Parente, Puzone, et al., 2018). It is marked by intense seismic activity occurring at depths between 10 and 20 km (Engdahl et al., 2006), including the Mw = 7.3 earthquake of 12 November 2017 (Tavani, Parente, Puzone, et al., 2018). The Simply Folded Belt is bounded to the NE by the High Zagros Fault (HZF). NE of the HZF, the Imbricate Zone includes the telescoped distal part of the original continental margin of the Arabian Plate (e.g., Tavani, Parente, Vitale, et al., 2018). The Simply Folded Belt and Imbricate Zone are formed by the typical succession of the Arabian margin (Upper Triassic to Quaternary), largely folded and faulted due to the continental collision, resting on top of crystalline basement (Casciello et al., 2009;Colman-Saad, 1978;Gavillot et al., 2010;Jassim & Goff, 2006;Law et al., 2014;Rudkiewicz et al., 2007;Sepehr & Cosgrove, 2004;Tavani, Parente, Vitale, et al., 2018). The Sanandaj-Sirjan Zone is formed mainly by metamorphic rocks, by Jurassic to Early Eocene calc-alkaline magmatic rocks, and by the products of Middle Eocene gabbroic plutonism (Alavi, 1994;Baharifar et al., 2004;Leterrier, 1985).

Seismicity
To produce a model of stress and strain accumulation during the interseismic stage and a simulation of the coseismic behavior of the MFT in western Lurestan, large-to moderate-magnitude earthquakes that nucleated in the study area are taken into account. Earthquake data for western Lurestan are available from the Global Centroid-Moment-Tensor (CMT) and U.S. Geological Survey (USGS) catalogues (USGS.gov| Science for a Changing World, 2020) starting from 1967. For the study area, we obtained earthquake data from the USGS catalogue (USGS.gov|Science for a Changing World, 2020). This data set includes 112 events with Mw ≥ 4.5 and 9 events with Mw ≥ 5.5 ( Figure 3). Focal mechanisms outline two distinct types of fault plane solutions, including thrust and strike-slip earthquakes. In particular, earthquakes nos. 2, 3, 4, 5, and 9 in Figure 3 are of thrust type with a depth ranging from 19 km (main seismic event, no. 4 in Figure 3) to 6 km (no. 2 in Figure 3). Tavani

Tectonic Geomorphology Analysis
Information on rock uplift along the investigated transect of the Zagros mountain belt was inferred from the qualitative and quantitative analyses of the features of the landscape and drainage network carried out using a Digital Elevation Model (the 30 m resolution ASTER GDEM), satellite images (Google Earth, 2019) and orthophotos, and a GIS-based analysis of the Digital Elevation Model (DEM). The qualitative analysis of the topography was focused on a segment of the transect centered on the epicenter of the 2017 earthquake ( Figure 4) and was aimed at collecting information that could discriminate between contributions exerted by variable erodibility of outcropping rocks versus surface motions in the formation of the landscape. This procedure included an analysis of the active and relic drainage network, carried out through, e.g., identification and mapping of abandoned fluvial paths (i.e., beheaded valleys and wind gaps) and points of capture, and identification of erosional and depositional landforms suggestive of a multistage development of topographic relief.
Quantitative constraints on the evolution of relief were obtained from the analyses of the spatial distribution of elevation and local relief and from river longitudinal profiles. The elevation parameter depends on both resistance to weathering of outcropping rocks and rock uplift. However, the maximum elevation is considered as primarily influenced by resistance to erosion of outcropping rocks, while the mean elevation is considered as more closely representing a response to surface uplift (England & Molnar, 1990). On the other hand, local relief is considered a robust indicator of rock uplift, particularly in regions underlain by rocks with rather homogeneous resistance to erosion (DiBiase et al., 2010). The elevation and local relief were analyzed both along a profile and in map view. A swath profile was constructed following the methodology of Pérez-Peña et al. (2017) along a 320 km long and 40 km wide transect centered on the trace of the cross section M-M′. Maps showing the spatial distribution of the mean elevation and local relief parameters were constructed for the area framing the analyzed transect using a moving window of 5 × 5 km.
The drainage network was analyzed by the construction of longitudinal profiles of rivers, the spatial distribution of the normalized channel steepness index (k sn ) and the construction of transformed river long profiles or chi plots. River long profiles are sensitive to both small-and large-scale tectonic perturbations. Deviations from the theoretical graded, concave-upward profile (Hack, 1957), that is, rectilinear or convex upward profiles, characterize river long profiles of areas of spatially variable uplift, spatially variable erodibility, and/or regions with a change in uplift rate at some point within the timeframe of the river response time (e.g., Attal et al., 2011;Duvall et al., 2004;Kirby & Whipple, 2001;Pazzaglia et al., 1998;Whittaker et al., 2008). River long profile analysis was applied to 22 rivers that include the trunk of the Diyala River that crosses the investigated area, five of its main tributaries from the SE, and rivers pertaining to their correlative hydrographic basins.
An investigation of the large-scale features of the river network was carried out by applying a slope/area analysis to all rivers dissecting the transect and plotting the results on a swath profile. The slope/area analysis we applied to the river network, using Topotoolbox (Schwanghart & Kuhn, 2010;Schwanghart & Scherler, 2014) relates the channel slope with the drainage area (Flint, 1974;Hack, 1957;Kirby & Whipple, 2001;Snyder et al., 2000;Whipple & Tucker, 1999). This analysis suggests that, at the reach scale, the slope of bedrock rivers is inversely proportional to the drainage area, following the hyperbolic equation: where S is the river long profile slope, k s is the steepness index, A is the drainage area, and Θ is the concavity index. For our analysis we adopted a smoothing window of 500 m and a reference drainage area A 0 = 1 km 2 .
To derive a normalized steepness index (k sn ), we used a reference concavity value of 0.41, which was obtained from a chi plot analysis to derive the best fit value. The k sn index is considered as a combined indicator of (i) active uplift, (ii) enhanced incision associated with knickpoints, and (iii) bedrock erodibility (e.g., DiBiase et al., 2010). The spatial distribution of the k sn index was also synthetized in a curve showing the mean k sn values centered along the trace of the swath profile. To obtain this curve, we interpolated the k sn values to convert them from a vector format to a raster format and thus derive the k sn map. We then applied the swath profile method to the k sn map and derived the curve of the mean values.
The features of 21 rivers were also investigated by transforming river longitudinal profiles into chi plots. The transformation removes the effect of the downstream increase in drainage area (which creates curvature), thus allowing for a meaningful comparison of river profiles at different spatial scales and with different uplift rates and erodibility values, besides visually enhancing knickpoints and transient signals (e.g., . The transformation of the river long profiles into this nondimensionalized form predicts that bedrock rivers dissecting a uniform substratum and equilibrated with uplift will have a linear chi plot, and the uplift rate will be reflected in the slope of the transformed profiles, which equates the steepness index . To obtain chi-transformed profiles, Equation 1 can be rewritten as follows: where U is the rock uplift rate, K is an erodibility coefficient, A is the drainage area, and m and n are constants. Separating variables in Equation 2, with the assumption of invariable U and K and integrating them, produces  Table 1).
where z(x) is the elevation of an observation point along the river long profile, z(x b ) is the elevation of the local base level, A(x) is the drainage area at the observation point z(x), and A 0 is a reference drainage area. Like for the k sn analysis, we adopted a smoothing window of 500 m and A 0 = 1 km 2 .
Variations from the linear shape of a chi plot may be due either to variable erodibility or temporally or spatially variable uplift . To verify whether a river is close to steady-state conditions, it is crucial to recognize the best fit m/n ratio at the drainage basin scale, whereas to compare chi plots among rivers in different drainage basins it is common to use the average m/n value among all of the analyzed population . Therefore, we first analyzed the single river long profiles and evaluated for each of them the m/n exponent of Equation 2 using the Matlab tool Topotoolbox (Schwanghart & Kuhn, 2010;Schwanghart & Scherler, 2014), determining the best fitting value for the m/n ratio. To compare chi plots from the 21 investigated rivers, we then calculated the mean m/n ratio and transformed the chi plots using the obtained value.

Finite Element Modeling
The modeled crustal section ( Figure 5) runs across the hypocenter of the 12 November 2017 earthquake ( Figure 2). To investigate both interseismic and coseismic deformation along the section, the methodology of 2-D FEM was used. This method consists of the geometric construction of a model containing the fault setting of interest. Using Marc software (MSC Software Corporation), the prebuilt model was divided into several domains, to which values of Young's modulus, Poisson's ratio, and density were assigned considering an elastic rheology. To resolve the system, the model was divided into an equivalent assemblage of small structures (mesh). As a result, for each unit, a solution was formulated and combined to obtain the solution for the entire system. In our instance, the geometric model was divided into three different homogeneous zones characterized by average elastic parameters and density values included in Table 1. A sensitivity analysis was carried out in order to quantify the impact of variations of density values and Young's modulus. Varying the density from 2,700 to 2,900 kg/m 3 (Basilici et al., 2019(Basilici et al., , 2020Teknik et al., 2019), and Young's modulus from 45 to 55 GPa in the Iran block and from 55 to 65 GPa in the basement (Zhao et al., 2004;Zhu & Zhang, 2013), we obtained uncertainties always smaller than 4%. In particular, the uncertainty varies nonlinearly along the cross section M-M′ for Young's modulus variations within the previously specified ranges, with a mean value of 2% for the central zone of the transect (175-225 km in Figure 5). On the other hand, the uncertainty is almost constant in the case of the density, with a maximum value of 3.7%.
The friction coefficient was set as μ = 0.6 according to Byerlee (1967Byerlee ( , 1978. The model was divided into 14,314 quadrangular elementary cells and 15,298 nodes. Near (i) faults, (ii) homogenous zone boundaries, and (iii) top surface (horizontal in the model), the sides of each single quadrangular cell are about 1 km long, increasing in size away from the contacts.
Numerical modeling included two independent FEM simulations carried out using the same mesh ( Figure 6): A first procedure was used to analyze interseismic stress and strain accumulations, while a second procedure was applied to investigate surface vertical motions associated with the MFT coseismic stage (by modeling the characteristic earthquake).
Modeling of the interseismic stage was divided into two different steps: The first step consisted of setting the boundary conditions of the model, so as to observe the effects of gravity alone. The boundary conditions applied to the model were the following: (i) The surface was free to move in all directions, (ii) the SW and NE boundaries were locked in the horizontal direction and free to move in the vertical direction, and (iii) the base was treated as a Winkler's foundation (Williams & Richardson, 1991). This model base was used to simulate the hydrostatic pressure of the Earth's mantle: Free horizontal movement was allowed and vertical motion was controlled by an elastic spring with stiffness coefficient K equal to K = A/L · E, where A is the base length, L is the thickness of the model, and E is the average of Young's modulus of the rocks included in the model. Moreover, a constant friction coefficient (μ) was attributed to the fault surfaces. The first step allows the achievement of "equilibrium conditions" between gravity, hydrostatic pressure of the Earth's mantle, and compaction of rocks and contacts. The second step, implemented sequentially, consisted of the horizontal movement of the SW boundary toward the NE boundary to simulate observed movements constrained by the GPS stations located at the surface. This procedure allowed us to compute the amount of equivalent Von Mises stress (σ VM ): where σ ij is the stress component and δ ij is the Kronecker delta,and of equivalent (Von Mises) strain (ε eq ): where ε ij is the strain component (Zhuang et al., 2019). This procedure also yields values of surface motion in the vertical direction, accumulated in a given time interval during the interseismic stage.  Basilici et al., 2019Basilici et al., , 2020. Values of Young's modulus and Poisson's ratio were defined according to previous FEM modeling by Zhao et al. (2004) and Zhu and Zhang (2013). Formations and groups belonging to the Arabian sedimentary cover (refer to Figure 4) are from Tavani

Tectonics
Modeling of the coseismic stage shares the first step with the interseismic simulation procedure. However, in this case, the second step consisted of the sudden movement of a part of the fault plane to simulate stick-slip behavior. This behavior was obtained by imposing a slip value consistent with the magnitude of the seismic event to be modeled (characteristic earthquake in our instance). The only GPS data available in our study area are those recorded by the ILAM station (Vernant et al., 2004) for the period 1999-2001. The related velocity projected along our model section results in a value of 2.3 ± 1 mm year −1 toward the NE, considering a fixed Iran block reference frame. Therefore, the motion of the SW boundary during the second step of the interseismic simulation was set up to move northeastward with a horizontal velocity of 2.3 mm year −1 in correspondence with the projected position of the ILAM GPS station in the model ( Figure 5).
Modeling of the interseismic stage ( Figure 6b) involves keeping the seismogenic MFT patch locked (this is the fault patch that is inferred to have slipped during 12 the November 2017 earthquake; Gombert et al., 2019;Nissen et al., 2019), together with the upper crustal fault splays branching out from the upper portion of the main MFT. On the other hand, the deeper (NE) portion of the MFT was let free to move (by stable sliding), thus simulating a creeping fault in the middle to lower crust. The calculation procedure was set up to cover an interseismic period of 1,000 years to reach appreciable values of accumulated stress, strain, and surface Step 2 for the interseismic stage: A horizontal northeastward velocity of 2.3 mm year −1 is introduced in correspondence of the projected position of the Ilam GPS station. The fault portion colored in white is that let free to move by stable sliding, while the rest of the fault (colored in red) is locked. (c) Step 2 for the coseismic stage: The boundary conditions are the same as those for Step 1, but the fault portion colored in white is let free to move instantaneously, while the rest of the fault (colored in red) is locked. Adopted rock units parameters are listed in Table 1. 10.1029/2020TC006402 Tectonics displacement. In contrast, the coseismic scenario ( Figure 6c) considers a portion of the fault located at a depth between 14 and 20 km free to move (by stick-slip), with a total slip of 4 m (as inferred for the 12 November 2017 earthquake; Vajedian et al., 2018).

Geomorphological Constraints on Surface Uplift
The relief of the Zagros mountain belt is strongly controlled by variable resistance to erosion of the outcropping rocks (e.g., Oberlander, 1968Oberlander, , 1985. Features such as anticlinal domes and elliptical hogbacks encircling breached anticlines, both occurring in hard carbonates, are common in areas where erosion has stripped the stratigraphically higher, weak rocks including Miocene shales and evaporites. These features are characterized by widespread exposure of the Oligocene Asmari Formation and, due to greater exhumation, also by Mesozoic carbonates (cropping out in the core of the anticlines; e.g., Oberlander, 1965;Tucker & Slingerland, 1996;Ramsey et al., 2008). The frontal part of the Zagros mountain belt in the investigated area is characterized by a low elevation (within few hundreds of meters above sea level) gradually increasing toward the NE, and a subdued relief, which is associated with smooth hogbacks and cuestas formed in poorly deformed strata of the Neogene foreland basin infill (Figure 7). The hogbacks and cuestas encircle eroded, breached anticlines and synclinal basins filled with alluvial deposits. A rise of the mean elevation up to values around 500 m a.s.l. in the mountain belt foothills roughly coincides with the occurrence of less eroded anticlinal ridges formed in the resistant Asmari carbonates and younger conglomerates. Toward the NE, the mountain front is marked by a sharp increase of the mean elevation, which peaks around 1,250 m at a distance of around 160 km from the SW edge of the swath profile ( Figure 7). Beyond that step, elevation values decrease below 1,000 m in a belt that includes the epicenter of the 2017 earthquake. This belt features a smooth relief, eroded in folded rocks originally underlying the Asmari carbonates-namely, foredeep sediments of the Gurpi Formation-punctuated by anticlinal ridges (e.g., Azgaleh and Mirinjeh anticlinal ridges; Figure 4) formed by carbonates of the Mesozoic Ilam Formation. From the Sheykh Saleh ridge toward NE, a rugged local relief formed in Cretaceous to Jurassic and Triassic carbonates alternating with marls ( Figure 4) characterizes the mountain belt, the elevation of which rises gradually up to~1750 m (Figure 7). A further elevation step located at a distance around 230 km and roughly following the surface trace of the HZF bounds the NE part of the profile, where the mean elevation fluctuates around 2,000 m and the local relief suddenly decreases. The increase of mean elevation beyond the 230 km step is fundamentally controlled by a rise in the minimum elevation, a parameter that, in the 5 × 5 km moving window used for the analysis of topography, may be considered as approximating the elevations of floors of the main valleys (e.g., Valente et al., 2019). Thus, the change of local relief around the 230 km step separates a deeply incised region, to the SW, from an elevated, smooth landscape to the NE.

Tectonics
The region between the two elevation steps identified along the swath profile has been investigated in detail to assess to what extent the distribution of the mean elevation along the profile responds to either the lithological or rock uplift controls, and to infer information useful to the reconstruction of the relief evolution. Although we attempted to minimize the lithological influence on the mean elevation portrayed in the swath profile by averaging across a large (40 km wide) belt, the elevation culmination at 160 km distance is affected by the presence of narrow ridges formed in hard rocks. Indeed, the mountain front is defined by the 1820 m high Mount Bamo ridge, which is one of a series of aligned hog backs formed by the Asmari carbonates that continues toward the south to the Maladizega ridge (Figure 8). Such a hog back alignment separates a low relief sector in the SW, characterized by alluvial reaches and punctuated by substantially flat, nonincised synclinal alluvial basins, from a region where alluvial basins are absent, and only thin veneers of terraced alluvial/colluvial deposits occur along the footslopes, or mark smooth surfaces eroded in weak rocks (Figure 8).
In the elevated area between the Azgaleh anticline and the Sheyik Saleh ridge, alluvial deposits are lacking, with the exception of small alluvial fans located in the piedmont of the Sheyik Saleh ridge and graded to the current river paths. That area is characterized by an overall smooth topography eroded in the deposits of the first foredeep infill (mainly composed of shales with conglomerate and sandstone layers) and the underlying carbonates of the Ilam Formation. Planar surfaces standing at~1,000 m a.s.l., or slightly higher, which are preserved on the interfluves, appear as incised remnants of a low-gradient landscape predating recent

10.1029/2020TC006402
Tectonics deepening of the drainage network ( Figure 8). Other few and sparse relic erosional features are (i) paleovalleys that, to the SW of the Sheyik Saleh ridge, are preserved on limited outcrops of the carbonate Ilam formation, namely, the Vanisar, Pshta, and South Sheyk Saleh (labeled SSS in Figure 8) anticlinal ridges, and (ii) shallow, NE draining beheaded valleys associated with the erosional surfaces.
In the area spanning from the Mount Bamo-Maladizega ridges to the Azgaleh anticline, thin alluvial fans originated from the Azgaleh and Miringeh anticlinal ridges and terraced surfaces eroded in weak rocks of the Gurpi Formation stand several tens of meters (around 30-50 m) above the current local base levels ( Figure 8). The crest of the Azgaleh anticlinal ridge is incised by a series of active and relic transverse fluvial paths. The oldest paleovalleys consist of wind gaps and beheaded valleys located above 1,250 m a.s.l. and are all characterized by a NE oriented paleoflow direction. Conversely, wind gaps and water gaps that incise the lateral termination of the Azgaleh anticlinal ridge at around 800-900 m a.s.l. all feature a SW oriented drainage ( Figure 8). The spatial distribution of those wind and water gaps points to reorganization of the drainage with progressive abandonment of elevated valleys in favor of lower drainages at the termination of the anticline. A similar distribution of wind and water gaps characterizes the Miringeh anticline ridge (Figure 8), where beheaded, SW dipping abandoned valleys occur (Figure 8). These active and relic drainages could result either from the growth of folds along their lengths (Keller et al., 1999), or from the progressive exposure of resistant rocks originally overlain by weaker lithology, both in folded and faulted terrains (e.g., Ascione & Cinque, 1999;Oberlander, 1968Oberlander, , 1985Ramsey et al., 2008). A key indicator of fold growth coeval with the development of drainage sculpted in the hard rocks is the spatial organization of the consequent drainage in the anticlinal ridges, which does not correspond to the present-day topography (Burberry et al., 2010;Keller et al., 1999;Keller & DeVecchio, 2013;Ramsey et al., 2008). However, in the Azgaleh and Miringeh anticlinal ridges, the consequent drainage is fundamentally consistent with the local slope orientations, thus pointing to the progressive exhumation of the hard rocks in the core of the anticlines (limestones of the Ilam Formation) as controlling capture phenomena and shifting fluvial paths toward both the NW and SE along the anticline crests. In the case of the Azgaleh anticlinal ridge, the original orientations of drainage along the short wind gaps are difficult to identify based on topography alone.

Drainage Features
The investigated transect is drained by the Diyala River, a tributary of the Tigris River, which flows transverse to the mountain belt, and tributaries of the Diyala River from the left (southeastern) side. The tributaries of the Diyala River include longitudinal rivers, characterized by narrow and elongated hydrographic basins trending NW-SE, roughly parallel to the main structural trend (hydrographic basins C, D, and E; Figure 9), which are dominant in the NE part of the transect. Mainly transverse rivers, trending oblique/ orthogonal to the strike of the structures in hydrographic basins A and B (Figure 9), are a distinctive feature of the SW part of the mountain belt. The occurrence in the NW portion of the Zagros orogen (including the Lurestan region) of transverse rivers that are not frequent in the SE part of the range (i.e., the Fars region), is related to its relatively wet climate and resulting effective drainage (Obaid & Allen, 2019). The main transverse rivers are characterized by some reaches aligned along the NW-SE oriented hogbacks and anticlinal ridges (hydrographic basin B; Figure 9). The boundary between these two main types of drainages runs along the crest of Mount Bamo and, southeastward, partly following the crest of the Azgaleh anticlinal ridge (Figure 9). This area roughly coincides with the transition from the elevated area that peaks with Mount Bamo, to the NE, to the low frontal part of the mountain belt. Information useful to reconstruct the evolution of the drainage network in the belt that follows the divide between hydrographic basins A and B, which are overall steeper in their upper reaches, and basin C, is inferred from the comparison of the longitudinal profiles of streams pertaining to those basins. This comparison suggests a tendency of the steeper streams of basins A and B to capture the less steep rivers of hydrographic basin C, which flow through the elevated region to the NE of Mount Bamo (Figure 9b; see section 5.3). Overall evidence suggests that piracy phenomena induced by downcutting by the transverse rivers in the hydrographic basins A and B were the main responsible for reorganization of the drainage network and related abandonment of the river paths suspended in the crest of the Azgaleh and Miringeh anticlinal ridges.

River Long Profiles and Derived Parameters
Information on the large-scale features of the river network can be extracted from the map of the normalized channel steepness index shown in Figure

10.1029/2020TC006402
Tectonics plotted against the maps of the local relief and mean elevation, respectively. Values of k sn are based on calculations for bedrock-substrate rivers (Kirby & Whipple, 2001, 2012Snyder et al., 2000;Wobus et al., 2006). In the investigated area, such a condition applies fundamentally to all rivers, except for some of those in the southwestern part of the transect, located in the intermontane basins, which include alluvial reaches. Since the analyzed streams were not individually selected (section 4.1), values from alluvial reaches (i.e., lower reaches) are shown for reference but are not considered in the interpretation (Buscher et al., 2017;Sklar & Dietrich, 1998;Stock & Dietrich, 2003;Whipple & Tucker, 2002). Along the transect, a central belt of high k sn values in the river middle reaches (marked by dashed white lines in Figure 10) can be distinguished from the areas to the SW and NE, respectively. The distribution of k sn values is not correlated with the distribution of precipitation. The latter increases from the lowlands (in the SW) toward the mountain range, without any substantial variation along the northeasternmost part of the transect (e.g., Obaid & Allen, 2019). In the region of high k sn , an overall SW to NE gradient of k sn values can be identified ( Figure 10). In particular, the increase of k sn values in correspondence of the white dashed line at the 160 km step correlates with the increase of mean elevation that characterizes the relatively low-relief region between Mount Bamo and the Sheyk Saleh ridge ( Figure 10). A belt of high k sn values is located to the NE of the Sheyk Saleh ridge (Figure 10b). Within that belt, the higher k sn values are associated with tributaries that dissect the valley flanks, thus possibly affected by rock resistance (e.g., local outcrops of upper Cretaceous rocks, namely, the Ilam carbonates). In the area to the NE of the 230 km topographic step (dashed white line in Figure 10) a general decrease of both k sn values and local relief coupled with high values of mean elevation occurs.
A more in-depth investigation of the drainage network with high k sn values was carried out through the construction of longitudinal profiles and chi plots for the main rivers and tributaries. The analyzed river population dissects a lithologically inhomogeneous bedrock (Figure 9a) of expected variable resistance to erosion. The longitudinal profile of the Diyala River shows that an extended knickzone separates the upstream reach, which dissects the high mean elevation and low local relief region in the northeasternmost part of the transect, from the much lower downstream reach (Figure 9). Profiles of rivers pertaining to hydrographic basins C, D, and E ( Figure 9b) are characterized by overall rectilinear to convex or poorly concave shapes and multiple slope breaks, which in some instances are associated with water gaps. The rivers with irregular profiles of hydrographic basins D, E and C (namely, Rivers 11 and 12) flow through the region of rugged relief to the NE of the Sheyk Saleh ridge (Figure 9a). Steep and irregular profiles characterize also Rivers 13 to 17 of basin C, which for most of their lengths dissect longitudinally the low local relief landscape located between the Azgaleh and Sheyk Saleh ridges. Rivers belonging to hydrographic basins A and B collectively feature concave profiles, although rivers of basin B appear less concave and are characterized by slope breaks and/or rectilinear reaches (Figure 9b).
We explored the features of each single river, with the exception of the Diyala river that is characterized by both bedrock and alluvial sections, with the construction of chi plots. The longitudinal profiles and chi plots for the 21 analyzed rivers are reported in Figure 11. Almost all of the transformed profiles are characterized by nonlinear trends and are punctuated by knickpoints ( Figure 11). Both features may reflect varying erodibility of the incised substratum, which may obscure transient signals related to changes in uplift rate both in space and time (upstream migrating knickpoints; Goren et al., 2014). To assess whether the knickpoints identified through the profile transformations reflect varying erodibility of the underlying bedrock or are an indication of transient river profiles, the locations of the knickpoints were compared with the bedrock lithology ( Figure 11) and knickpoints located at lithological contacts, labeled lithology-controlled knickpoints, were distinguished from nonlithology controlled ones (Figure 9a). Our findings highlighted that several rivers of basins A, B and C (Rivers 1,2,7,8,10,13,15,16, and 17; Figure 11) are affected by net slope changes at outcrops of the Ilam Formation carbonates, while they are insensitive to the main lithological contacts between the upper and lower units of the first foredeep infill. However, a hard sandstone interval included in the upper unit (not mapped in the geological map of Figure 4 and Figure 9a) has shown to affect the river profiles with knickpoints. For rivers of basins D and E, a correlation of some knickpoints with the contacts of the Cretaceous marls with the underlying Jurassic and Triassic rocks (limestones and dolostones), besides those between the Ilam carbonates and the Gurpi Fm., are evidenced ( Figure 11). However, not all of the changes of slope, concavities, or convexities are correlated with changing bedrock lithology.

10.1029/2020TC006402
Tectonics Figure 11. Longitudinal profiles and transformed longitudinal profiles (chi plots) of 21 bedrock rivers that dissect the area shown in the map of Figure 9a (river numbering as in Figure 9a). The chi plots were constructed using a smoothing window of 500 m and a reference drainage area A 0 = 1 km 2 . The colored bars above the long-profile diagrams indicate lithology of the river bedrock.

Tectonics
To explore the nature of the signals that govern the nonlinearity of the chi plots, and to analyze their spatial distribution, we compared the features of the studied rivers following the procedure described in section 4.2.
For the analyzed river population, the best fit m/n value is 0.41 ± 0.26, a value in line with global scale evaluations of mean m/n value for "arid" environments (Harel et al., 2016). The chi plots constructed for each hydrographic basin are shown in Figure 12. For hydrographic basin A, the transformed profiles of fluvial paths located to the SW of the MFF (Rivers 3 to 6 and lower segments of Rivers 1 and 2; Figure 12) are all linear, with a low slope angle. Upstream, the main river and tributary (Rivers 1 and 2; Figure 12) show comparable increase of slope independent from the bedrock lithology ( Figure 11). Their shapes and slopes are similar to those of Rivers 7 and 8 of hydrographic basin B that, as Rivers 1 and 2, cross the MFF, while the chi plot of River 9, which flows for much of its length to the SE of the Mount Bamo culmination, is less steep. The main rivers and tributaries of basins A and B (Rivers 1 and 2, and 7 and 8, respectively) are all characterized by upper knickzones that plot at about the same w(x) position and elevation. In hydrographic basin C, two main groups of rivers may be distinguished. Beyond the knickpoints in their lower courses, which are located at lithological contacts, the transformed profiles of tributaries that flow between the Azgaleh and Sheyk Saleh ridges (Rivers 13 to 17) are rather colinear and less steep than those of rivers flowing to the NE of the Sheyk Saleh ridge. The latter rivers display, in their middle parts, a slope break (River 12) and a convex shape (River 11). As the chi plots for Rivers 11 and 12, the overall composite profiles of rivers of hydrographic basins D and E are all characterized by multiple slope breaks.

Finite Element Model
The values of resulting total strain, equivalent of Von Mises stress, and vertical surface displacements for both the interseismic and coseismic scenarios are shown in Figure 13. For the interseismic stage, a zone of maximum resulting total strain occurs in the hanging wall of the MFT, above the deeper portion of the seismogenic thrust segment. This region of maximum deformation is comprised between the Sheykh Saleh Fault to the NE and the Miringeh Fault to the SW (Figure 13a), being located just NE of the 12 November 2017 earthquake hypocenter. The latter main shock falls within a zone of relatively high strain (in the range of 9.0E−5 to 9.5E−5). It is worth noting that the zone of marked strain accumulation reaches the surface, maintaining similar values to those characterizing the locked thrust segment at depth. The accumulated strain decreases gradually both northeastward and southwestward, defining an~140 km wide perturbed area. The Von Mises stress, besides displaying an expected peak (exceeding a value of 6.2E+6 Pa) at the junction between creeping and locked thrust segments, is characterized by roughly elliptical, concentric regions of high stress elongated in a subhorizontal direction (Figure 13b). Stress accumulation at the surface is not as marked as that occurring at depth (particularly in the 10-20 km range), as stress diffusion appears to follow a horizontal preferential direction. In fact, the perturbed (stressed) region exceeds 175 km in the horizontal direction (i.e., the whole model area shown in Figure 13b).
The model output of vertical surface displacement for the interseismic stage (integrated over a 1,000 year time span) is characterized by a plateau in the SW part of the model, increasing from 150 km northeastward to define a wide bulge between 190 and 240 km and then gently decreasing to the NE (Figure 13c). For the coseismic stage, the modeled hypothesis of a characteristic earthquake of Mw = 7.3 yields curves of vertical surface displacement characterized by two peaks (Figure 13d). The first peak is positive and shows an increment from 0 to 1.2 m at a distance around 162 km, while the second peak is negative and shows a

Discussion
The expression at the surface of blind thrusts underlying mountain ranges is strongly dependent on several factors governing the shaping of topographic relief besides the fault plane geometry. These include climate, resistance to weathering of outcropping rocks and fluvial erosional/depositional dynamics, which is also affected by local base levels that flank the range (e.g., Ellis & Densmore, 2006;Forte et al., 2014;Whipple & Tucker, 1999;Willett et al., 2001). Although resulting from the interplay between all those variables, topography preserves features that may help to unravel the time-space distribution of vertical

10.1029/2020TC006402
Tectonics motions in mountain ranges, particularly if the analysis of topography is combined with an analysis of the drainage network and reconstructions of its development. All of such information provides constraints crucial to the reconstruction of the geometry and long-term temporal evolution of structures at depth (e.g., Eizenhöfer et al., 2019;Ellis & Densmore, 2006). In our study, this approach allows us to constrain and test the modeling of the development of deep-seated structures of well-known geometry. In Figure 14, the main large-scale features of topography and river network along the transect M-M′ are synthesized and compared. In particular, Figure 14 shows the spatial variation of normalized steepness index, which is synthesized in terms of mean k sn values, and transformed longitudinal river profiles, revealing fundamental features on both generalized (i.e., large scale) and differential uplift (e.g., Whipple & Tucker, 1999;Wobus et al., 2006). The spatial pattern of k sn values may be affected by the distribution of precipitation, particularly in arid regions; however, even in arid/subarid environments where the steepness index varies nonlinearly with discharge and erosion rate, k sn is considered a key indicator of rock uplift in active mountain ranges (e.g., Bookhagen & Strecker, 2012;DiBiase et al., 2010;DiBiase & Whipple, 2011).

Varying Erodibility of Outcropping Rocks
Although the control exerted by varying erodibility of outcropping rocks on the topography features of the investigated transect is difficult to be isolated, geomorphological information allows us to perform qualitative evaluations. Consistent with findings from other regions of the Zagros mountain range

10.1029/2020TC006402
Tectonics (e.g., Oberlander, 1965;Tucker & Slingerland, 1996;Ramsey et al., 2008), several lines of evidence indicate that in our study area, the carbonate rocks of the Asmari and Ilam formations have lower erodibilities compared with the other outcropping rock types. In particular, (i) the coincidence of the main topographic highs with outcrops of carbonates of the Asmari Formation and Ilam Formation and (ii) the occurrence of slope breaks in the longitudinal profiles and chi plots of the analyzed rivers at contacts between these carbonates and rocks of the first foredeep infill (Figures 9a and 11), point to the above mentioned limestones as resistant rocks. Furthermore, based on evidence from the chi plots, the upper and lower units of the first foredeep infill appear characterized by overall comparable erodibility. Our analyses have also highlighted that within the rocks of the first foredeep infill, even conglomerate and sandstone intervals a few tens of meters thick appear less erodible than the shales, which constitute most of that unit (section 5.3 and Figure 11).
Based on the local relief and, particularly, stream profiles and chi plots (Figure 11), we infer that the rocks of the first foredeep infill are overall less resistant than the underlying Mesozoic succession, which crops out in the Sheyk Saleh and ridges to the NE of it. In addition, the occurrence of slope breaks at contacts between the Cretaceous marls and underlying Jurassic carbonates and, for tributaries that dissect the valley flanks, high k sn values that can be associated with relatively hard rocks (section 5.3 and Figures 9-11) are suggestive of varying erodibility within the Cretaceous to Triassic rocks. In comparison to the Mesozoic succession, the rocks outcropping in the Imbricate zone appear as less resistant, as it can be inferred from the features of the river profiles ( Figure 11).

Relief Evolution in the Mount Bamo Area and Distribution of Uplift Along the Transect
Along the transect, the elevation rise at 160 km distance corresponds to the roughly N-S trending Mount Bamo-Maladizega alignment of hog backs, representing the remnants of a breached, pronounced anticlinal structure that marked the mountain front. Although being severely eroded, the Mount Bamo-Maladizega ridges mark the net topographic step of the MFF, which appears as the most reliable topographic expression of the MFT. In the area between the Mount Bamo-Maladizega hog backs and Sheyk Saleh ridge, stripping of the carbonate carapace has exposed rocks of the first foredeep infill and, where erosion was greatest, the originally deeper carbonates of the Ilam Formation (Figures 4, 8, and 9). There, the occurrence of rock types of variable resistance to weathering may have affected the topographic expression of rock uplift, making the signature of recent displacements subdued.
In the belt between Mount Bamo and the Sheyk Saleh ridge, aspects of the topographic relief follow the underlying fold structures, this area being characterized by well-preserved anticlinal ridges incised by active and relic valleys. However, both the relic and active drainages do not show any of the plan view forms that are typically associated with the vertical/lateral development of folds (e.g., Keller et al., 1999;Keller & DeVecchio, 2013), which conversely are widespread in the anticlinal ridges in the Fars region and in the frontal part of the Zagros belt in the Kurdistan region of Iraq (e.g., Burberry et al., 2010;Ramsey et al., 2008;Zebari et al., 2019). Such evidence indicates that even the paleovalleys located at highest elevations (above 1,250 m) along the crest of the Azgaleh ridge, which represent the oldest ones in the area mapped in Figure 8, substantially postdate short wavelength folding. This implies that drainage features possibly coeval with the development of individual folds have been erased as the stratigraphically higher rocks were stripped away. Our reconstruction is consistent with folding in the outer Lurestan province starting at circa 7.6 Ma (i.e., in the late Tortonian) and predating the early Pliocene (<5 Ma; Homke et al., 2004;Koshnaw et al., 2017) initiation of uplift, which caused shifting toward the SW of the paleo-Tigris path (Vergés, 2007) and likely governed the establishment of a transverse drainage in the investigated area. However, the location of the paleodivide that is inferred from the NE oriented transverse paleodrainage (Figure 8) suggests the occurrence in the past of a topographic high located to the west/southwest of the divide that separates hydrographic basins A-B and C (Figure 9a), that is, roughly aligned with the Mount Bamo culmination (Figure 8). In fact, the NE oriented paleovalleys that stand above 1,250 m appear as consequent drainages developed in the rear flank of a growing topography that has continued to be uplifted well after breaching and stripping of the Asmari cover. In the area of the Azgaleh anticlinal ridge, the NE flowing drainage has been maintained until, in response to piracy phenomena triggered by deepening of rivers (pertaining to the modern hydrographic basin A and B) that dissect the mountain front and related lowering of the topographic surface, an inversion of drainage was recorded at the lateral terminations of the anticline. Such a reconstruction is supported by the correspondence of the points of capture in the uppermost trunks of 10.1029/2020TC006402

Tectonics
Rivers 1 and 7 ( Figure 8) with knickpoints (Figures 9a and 11). In other words, overall evidence suggests that topographic growth of the roughly N-S oriented belt that includes the Mount Bamo culmination on a scale larger than that of individual short wavelength folds has continued after the development of those folds and subsequent stripping of the hard carbonate carapace. Such a continuing differential uplift of the Mount Bamo ridge, which would have also favored the maintenance of a longitudinal drainage in hydrographic basin C, is consistent with continuing thrust activity at depth.
The MFF limits a belt of high-gradient increase of mean elevation and rugged relief that separates the frontal region of low mean elevation and relief of the mountain range from the hinterland zone, where only a subdued increase of mean elevation occurs and the relief decreases ( Figure 10 and Figure 14). Slow uplift of the frontal region to the SW of the MFF is inferred by low values of mean elevation, local relief and k sn , besides the river chi plots. In fact, the straight and smooth transformed long profiles of fluvial paths limited to the SW of the MFF, which are rather completely colinear and with low slope angle ( Figure 12 and Figure 14a), are indicative of rivers that dissect bedrock of overall homogeneous erodibility equilibrated to slow uplift. To the NE of the 230 km elevation step, the topographic features of the Imbricate zone and, particularly, the Iranian block are suggestive of an elevated smooth landscape overall adjusted to high-standing local base levels. Multiple evidences such as the occurrence of a large knickzone in the Diyala River profile (Figure 9b), the general decrease of k sn values and local relief coupled with high values of both the mean and minimum elevations (which are proxies for the elevation of the valley floors; Figure 14), are collectively indicative of a "deceleration" of vertical incision in that elevated area. Therefore, multiproxy evidence highlights that the smooth, elevated landscape located to the NE of the HZF has not yet been reached by local base level lowering, which affects the rugged relief region to the SW. In contrast, a belt of high k sn values spans over the region of high local relief located between the topographic steps at 160 and 230 km ( Figure 10 and Figure 14a), which appears as a locus of accelerated incision of a more mature landscape. The high k sn belt exceeds the width of high local relief by less than 10 km on both sides. Coupling of high k sn values and high local relief is considered an indication of localized uplift (e.g., DiBiase et al., 2010;Forte et al., 2014) and, specifically, the spatial distribution of high k sn values is used to identify areas subject to active rock uplift and infer the distribution of both uplift and horizontal motions in mountain ranges (e.g., Eizenhöfer et al., 2019, and references therein). In our study area, the belt of coupled high relief and k sn encompasses rock types of inferred variable erodibility ( Figure 14) and as such the region spanning from the MFF to the HZF may be considered as a locus of relief growth in response to concurrent long-term uplift and river incision. Within that region, the local relief and k sn mean values are variable, with higher values being associated with the Mount Bamo-Maladizega alignment and greatest k sn and local relief (exceeding 2000 m; Figure 14a) with the area to the NE of the Sheyk Saleh ridge (Figure 14a). Such a distribution suggests inhomogeneous uplift and, consistent with findings from the detailed geomorphological analysis, overall slower uplift in the belt between the Mount Bamo alignment and the Sheyk Saleh ridge.
More detailed information on river steepness is extracted from the transformed river profiles and, particularly, from chi plots of the rivers of basins A and B that cross the MFF. Despite their different substratum, the main rivers and tributaries of basins A and B (namely Rivers 1, 2, 7, and 8; Figures 11, 12, and 14a) show comparable concave-up slope changes that plot at approximately the same elevation (of 1,100-1,200 m), a feature that may be indicative of migrating knickpoints. Although it is generally difficult to assess if a knickpoint is migrating or stationary (e.g., Demoulin et al., 2017;Goren et al., 2014;, and references therein), the lack of correlation of the considered slope changes with either changes of bedrock lithology ( Figure 11) or known structures possibly accounting for differential vertical motions, suggests that they reflect a common transient signal originated at about the same position and propagated upstream (e.g., . Thus, overall features of the chi plots of Rivers 1, 2, 7, and 8 are compatible with uplift originated at the MFF, although the interpretation of those composite transformed profiles in terms of recent decrease/increase of the uplift rate is difficult. This difficulty arises from the dependence of the transformed profile shape on the slope exponent n of Equation 2 (particularly, whether n is greater or less than 1; , which cannot be evaluated unless the rates of erosion and uplift are independently constrained (e.g., Demoulin et al., 2017;Mudd, 2017, and references therein). Consistent with evidence from k sn mean values, the transformed profiles of rivers located between the Azgaleh and Sheyk Saleh ridges are suggestive of rivers essentially equilibrated to less pronounced uplift relative to the region to the NE (Figure 14a). The latter region is dissected by rivers characterized by more irregular to convex chi plots.

Tectonics
Although the main slope changes can be related with changing bedrock erodibility (Figure 11; while the main break in the slope of River 11 appears as a response to capture phenomena), overall convex shapes of chi plots in the right part of Figure 14a and the noncorrelation of some knickpoints (e.g., the lower ones in the main trunks of basins C and D, namely, Rivers 11 and 19; Figure 11) with bedrock lithology, collectively suggest transient signals.

Activity of the MFT and Orogenic Growth
Combined geomorphological-stratigraphic evidence from the investigated region suggests that the recent uplift has reshaped the topography constructed by earlier deformation. In particular, the widespread exposure of units stratigraphically underlying the Asmari Formation (e.g., rocks of the first foredeep infill and Ilam Formation) in the region where the modern divide between hydrographic basins A-B and C is located (Figure 9) may be interpreted as resulting from topographic growth of that region predating the rise of the Mount Bamo-Maladizega alignment, in which patches of the hard Asmari cover are still preserved. Only the later uplift of the MFF would have affected the landscape of the investigated transect with the drainage evolution synthesized in section 6.2. In contrast, the more mature landscape of the region to the NE of the Sheyk Saleh ridge, where the stratigraphically lowermost (Mesozoic) units outcropping in the entire study area are largely exposed, appears as a region that experienced longer uplift. This reconstruction implies a correlation of mountain range uplift with horizontal advection in the long term. Consistently, the spatial distribution of k sn in relation to width of the inferred belt of focused uplift is also suggestive of horizontal advection accompanying uplift, in agreement with recent modeling results by Eizenhöfer et al. (2019).
Thermochronological data indicate that uplift of the mountain range initiated in the Pliocene (after 5 Ma; Homke et al., 2004;Koshnaw et al., 2017), earlier than the uplift at the MFF, which occurred likely in response to the activation of basement thrusting at the mountain front (and related development of the MFT) later than circa 3 Ma . Superposed on the long-term trend discussed so far is recent uplift that, based on evidence from topography and river network parameters, has been focused into two main zones: the MFF and the area to the NE of the Sheyk Saleh ridge. Within this framework, our results suggest that the vertical component of motion along the basement-involved MFT has played a major role in the recent orogenic growth of the investigated mountain range. The style of long-term growth of the topography associated with the MFT reconstructed so far can be compared with the output of our finite element model ( Figure 14). To carry out the comparison, vertical surface displacements calculated for both interseismic and coseismic stages have been integrated over a time span of 1,000 y. As data on the recurrence interval of the characteristic earthquake are not available, the reasonable hypothesis of a characteristic earthquake occurring every 1,000 ± 500 years on average has been adopted to produce a cumulative vertical surface displacement curve (Figure 14b). It is worth noting that the uncertainty involved in the extrapolation of coseismic displacements over a 1,000 year period is not a major issue in this context, as our focus is on defining the position of areas of major surface deformation rather than on the absolute values of vertical displacement in Figure 14b. Within this framework, a likely contribution of postseismic afterslip taking place during years or even decades following the main shock (e.g., Copley, 2014;Copley & Reynolds, 2014) would also be taken into account by our modeling of surface displacement, assuming slip occurred along the same thrust originating the earthquake. A different issue is that the reference timescale of the finite element model (e.g., 10 2 -10 3 years) differs by several orders of magnitude from that involved in the construction of modern topography (e.g., 10 6 years). Nevertheless, the modeled coseismic and interseismic behaviors of the MFT provide useful insights into the source of the late part of the surface displacement recorded by topographic features and river network, also allowing one to infer whether one or more structures are contributing to such a displacement. Our results show that the local high relief in correspondence of the Mount Bamo-Maladizega ridge, whose SW slope defines the MFF, also occurs in the region of maximum accumulated strain, which reaches shallow depths according to finite element modeling of interseismic deformation (Figure 13a).
A further broad maximum of vertical surface displacement in the model occurs at~215 km, the uplift decreasing then slightly to the NE (Figure 14b). Consistently, indications from the DEM-based drainage analyses highlight that this same region is experiencing active uplift. Whether such an uplift is occurring at a rate comparable to-or, as it could be inferred from the k sn magnitude, even higher than-that of the MFF, is questionable, as the channels with less steep reaches to the SW of the Sheyk Saleh ridge dissect a 10.1029/2020TC006402 Tectonics mainly shaly bedrock that is overall weaker than the substratum of the rivers to the NE. Overall, this entire area of the model defines an uplifted crustal block in the hanging wall of the midcrustal, gently dipping segment of the MFT. Diagrams (c) and (d) in Figure 13 point out how the cumulative vertical surface motion pattern of Figure 14 results from rather distinct contributions provided by (i) continuous regional uplift characterizing interseismic stages (uplift starting at approximately km 150, reaching a maximum at approximately km 210) and (ii) focused deformation in the area of the Mount Bamo-Maladizega ridge (uplift starting at approximately km 130, reaching a maximum at approximately km 160) associated with coseismic displacement. Therefore, our FEM suggests that the Quaternary development of the relief defining the MFF is mainly related to coseismic deformation, while generalized uplift of the orogen segment located more to the hinterland (i.e., to the NE) is associated with stable sliding along the deeper portions of the MFT. Consistent with the results of numerical modeling of topography growing over seismically active blind thrusts (Ellis & Densmore, 2006), the two uplifted areas are separated by a belt characterized by coseismic subsidence associated with the large (M > 7) earthquakes occurring in the region. This belt defines a saddle clearly imaged by both the modern geomorphological setting of the investigated transect and relic landforms indicating the persisting occurrence of an elongated topographic low flanking the SW slope of the Sheyk Saleh ridge.

Conclusions
The landscape of the investigated transect represents the expression of a complex combination of uplift and denudation acting over a lithologically inhomogeneous substrate of variable erodibility. Coupling qualitative and quantitative analyses of the topography allowed us to unravel the major controls on the development of the relief and provided information on the time-space distribution of uplift. Overall information points to the locus of long-term uplift migrating toward the foreland through time. The reconstructed pattern of more recent uplift indicates that differential motions are affecting the entire region to the SW of the HZF. There, the elevated, breached anticline of M. Bamo and the region to the NE of the Sheyk Saleh ridge are being uplifted faster than the saddle between them. Such a scenario is consistent with the output of a finite element model that shows how the pattern of vertical surface displacement is the result of a combination of slip accumulated during large (M > 7) seismic events and continuous displacement along a gently dipping, midcrustal thrust. Finite element modeling of interseismic and coseismic stages allowed us to gain new insights into the relative contribution of each process in the development of the relief: While interseismic deformation produces a generalized uplift of the whole crustal block in the hanging wall of the mid crustal segment of the major thrust, coseismic displacement controls localized uplift of a distinct topographic feature located above the main upper crustal ramp of the same thrust, defining the prominent geomorphological boundary known as the MFF.

Conflict of Interest
There are no financial conflicts of interests for any author.

Data Availability Statement
The data supporting the interpretations and data products related to this paper are available in the Figshare repository (Basilici,