Relative Timing of Uplift along the Zagros Mountain Front Flexure (Kurdistan Region of Iraq): Constrained by Geomorphic Indices and Landscape Evolution Modelling

The Mountain Front Flexure marks a dominant topographic step in the frontal part of the Zagros Fold-Thrust Belt. It is characterized by numerous active anticlines atop of a basement fault. So far, little is known about the relative activity of the anticlines, about their evolution, and about how crustal deformation migrates over time. We assessed the relative landscape 10 maturity of three along-strike anticlines (from SE to NW: Harir, Perat, and Akre) located on the hanging wall of the Mountain Front Flexure in the Kurdistan Region of Iraq to identify the most active structures and to get insights into the evolution of the fold-thrust belt. Landscape maturity was evaluated using geomorphic indices such as hypsometric curves, hypsometric integral, surface roughness, and surface index. Subsequently, numerical landscape evolution models were run to estimate the relative time difference between the onset of growth of the anticlines, using the present-day topography of the Harir Anticline as a base 15 model. A stream power equation was used to introduce fluvial erosion, and a hillslope diffusion equation was applied to account for colluvial sediment transport. For different time steps of model evolution, we calculated the geomorphic indices generated from the base model. While Akre Anticline shows deeply incised valleys and advanced erosion, Harir and Perat anticlines have relatively smoother surfaces and are supposedly younger than the Akre Anticline. The landscape maturity level decreases from NW to SE. A comparison of the geomorphic indices of the model output to those of the present-day Akre 20 Anticline topography revealed that it would take the Harir Anticline about 80-100 kyr and 160-200 kyr to reach the maturity level of the Perat and Akre anticlines, respectively, assuming erosion under constant conditions and constant rock uplift rates along the three anticlines. Since the factors controlling geomorphology (lithology, structural setting and climate) are similar for all three anticlines, and under the assumption of constant growth and erosion conditions, we infer that uplift of the Akre Anticline started 160-200 kyr before that of the Harir Anticline, with the Perat Anticline showing an intermediate age. A NW25 ward propagation of the Harir Anticline itself implies that the uplift has been independent within different segments. Our method of estimating the relative age difference can be applied to many other anticlines in the Mountain Front Flexure region to construct a model of temporal evolution of this belt.


Introduction
The Zagros Fold-Thrust Belt is an active orogen that resulted from the collision between the Arabian and Eurasian plates and contains the deformed portions of the NE part of the former Arabian passive margin ( Fig. 1; Berberian, 1995;Mouthereau et al., 2012). Many aspects of the structural configuration and the evolution of the Zagros Fold-Thrust Belt are by now satisfactorily constrained, but the detailed spatial and temporal distribution of deformation across the belt is not yet well 5 understood. This concerns especially the NW part of the belt in the Kurdistan Region of Iraq (KRI) due to a lack of comprehensive studies and for geopolitical reasons that make access to the field challenging. The style, timing, and relative activity of front thrusts, deformation propagation, and along-strike variations have not been sufficiently studied. It is not wellknown which structures are currently the most active ones either.
One of the morphologically most conspicuous structural elements of the Zagros Fold-Thrust Belt is the Mountain Front Flexure 10 (MFF), which separates the High Folded Zone and the Foothill Zone (known in Iran as the Zagros Simply Folded Belt and Zagros Foredeep, respectively;Figs. 1 and 2;Berberian, 1995;Jassim and Goff, 2006;McQuarrie, 2004;Mouthereau et al., 2012;Vergés et al., 2011). In most parts of the Zagros, the MFF marks a pronounced topographic step, separating folds with high amplitudes, narrow wavelengths, and higher topography in the High Folded Zone from folds with relatively low amplitudes, long wavelengths, and lower topography in the Foothill Zone (Fig. 2). The MFF is characterized by numerous 15 active anticlines atop of fault strands emerging from a basement fault. It was suggested that the onset of the MFF activity in the NW Zagros was about 5±1 Ma based on low temperature thermochronology (Koshnaw et al., 2017). The timing of this activity is expected to differ along-strike the belt and, hence, the initiation of uplift of the anticlines on the hanging wall of the MFF is the key to understand this temporal and spatial evolution. In the neighbouring Iranian part, the MFF was a relatively long-lived structure active from 8.1 to 7.2 Ma to about the Pliocene-Pleistocene boundary. After that, only the 20 southwesternmost anticline remained active in front of the MFF. This was inferred from progressive unconformities and magnetostratigraphy (Hessami et al., 2001(Hessami et al., , 2006Homke et al., 2004).
In active orogens, the main factor that contributes to building up topography is ongoing convergence (Bishop, 2007;Burbank and Anderson, 2012;Whittaker, 2012). Recent advancements in the availability of high-resolution digital elevation models (DEMs) and GIS software allowed to quantitatively analyse the landscape (Bishop, 2007;Tarolli, 2014). Tectonic 25 geomorphology approaches and landscape maturity studies have been used extensively and proven to be efficient in studying the relative tectonic activity of different areas in contractional settings (Cheng et al., 2012;Mahmood and Gloaguen, 2012;Ramsey et al., 2008;Regard et al., 2009). Nevertheless, the NW part of the Zagros lacks modern studies on tectonic geomorphology -with few exceptions. Bretis et al. (2011) detected sets of wind gaps (i.e. segments of river valleys abandoned due to lateral and vertical fold growth) in the High Folded Belt, NE of the MFF, suggesting that larger folds grew by linkage 30 of smaller, shorter folds. Zebari and Burberry (2015) performed detailed analyses of various geomorphic indices for numerous anticlines in the High Folded Zone, concluding that the combination of clearly asymmetric drainage patterns and the mountain front sinuosity index (Bull, 2007;Keller et al., 1999) is a valuable tool for identifying putatively active fault-related folds. Obaid and Allen (2017) studied the landscape maturity of various anticlines within the Zagros Foothill Zone and constrained the order of deformation of these anticlines by proposing an out-of-sequence propagation of underlying faults into the foreland.
They proposed that the Zagros Deformation Front was among the earliest faults that have been reactivated within the Foothill Zone.
In an active orogen such as the Zagros, a better understanding of the temporal and spatial distribution of deformation due to 5 ongoing tectonics can be achieved with landscape modelling. In the last two decades, numerical models have been extensively used to study landscape evolution (Chen et al., 2014;Tucker and Hancock, 2010) and several software packages were specifically developed for this purpose (e.g. Hancock and Willgoose, 2002;Hobley et al., 2017;Refice et al., 2012;Salles and Hardiman, 2016;Tucker et al., 2001). Most of these models include algorithms for bedrock fluvial incision and hillslope creep as input parameters. Several studies have constrained the landscape evolution with the involvement of the corresponding 10 tectonics and structures elsewhere (Collignon et al., 2016;Cowie et al., 2006;Miller et al., 2007;Refice et al., 2012).
In this study, we assessed variations in the landscape maturity of three anticlines (from SE to NW, the Harir, Perat and Akre anticlines) located on the hanging wall of the MFF by quantitatively analysing landscape indices (hypsometric curve, hypsometric integral, surface roughness, and surface index) in order to distinguish more mature segments from less mature ones, and to reconstruct the relative variation of uplift time and/or rate along these anticlines. We then computed the difference 15 in the onset of uplift between more mature anticlines and less mature ones using a landscape evolution model. The presentday topography of the least mature anticline served as an input model for computing the time that it takes this anticline to reach the same state as the most mature one. Also, three structural cross-sections were constructed across the three anticlines to delineate their structural style and to link it with their landscape maturity.

Geological Setting 20
The Zagros Fold-Thrust Belt is the result of the collision between the Arabian and Eurasian plates (Fig. 1;Berberian, 1995;Mouthereau et al., 2012). Continental collision started in the Early Miocene following the progressive subduction of Neo-Tethyan oceanic lithosphere underneath Eurasia (Agard et al., 2011;Csontos et al., 2012;Koshnaw et al., 2017;Mouthereau et al., 2012). The Zagros Fold-Thrust Belt extends for about 2000 km from the Strait of Hormuz in southern Iran to the KRI and further into SE Turkey. Since the onset of collision, the deformation front has propagated 250-350 km southwestward, 25 involving the northeastern margin of the Mesopotamian foreland basin and the Persian Gulf into a largely NW-SE-trending foreland fold-thrust belt (Mouthereau, 2011;Mouthereau et al., 2007). The shortening across different sectors of the Zagros Fold-Thrust Belt is estimated to range between 10% and 32% (Blanc et al., 2003;McQuarrie, 2004;Molinaro et al., 2005;Mouthereau et al., 2007;Vergés et al., 2011). GPS-derived horizontal velocities between Arabia and Eurasia show presentday convergence rates between 19 and 23 mm/yr (McClusky et al., 2003). It is suggested that deformation partitioning occurs 30 between the external and internal portions of the Iranian part of the Zagros Fold-Thrust Belt. While the internal Zagros Fold-Thrust Belt currently accommodates 3-4 mm/yr of right-lateral displacement along the Main Recent Fault ( Fig. 1; Reilinger et al., 2006;Vernant et al., 2004), the external part accommodates 7-10 mm/yr of shortening by thrusting and folding (Hessami et al., 2006;Vernant et al., 2004), 2-4 mm/yr of which is taken up by the MFF in the Fars Arc (Oveisi et al., 2009). However, no such estimates are available for the Iraqi segment of the Zagros Mountains. It is hence not known how much of the total Arabia-Eurasia plate convergence is being accommodated across the Iraqi part of the Zagros Fold-Thrust Belt.
The NW segment of the Zagros Fold-Thrust Belt in the KRI is subdivided into several NE-trending morphotectonic zones. 5 These zones from NE to SW are: (i) Zagros Suture, (ii) Imbricated Zone, (iii) High Folded Zone and (iv) Foothill Zone (Figs. 1 and 2; Jassim and Goff, 2006). These zones are bounded by major faults in the area. The faults include Main Zagros Thrust separating the Zagros Suture from the Imbricated Zone, High Zagros Fault that separates the Imbricate Zone from the High Folded Zone, and the Mountain Front Flexure that separates the High Folded Zone from the Foothill Zone (Figs. 1 and 2; Berberian, 1995;Jassim and Goff, 2006). 10 The deformed sedimentary succession is composed of 8 -12 km thick Paleozoic to Cenozoic strata that rest on the Precambrian crystalline basement (Aqrawi et al., 2010;Jassim and Goff, 2006). The thick sedimentary cover consists of various competent and incompetent rock successions separated by detachment horizons. The infra-Cambrian Hormuz salt, which acts as a basal detachment in much of the Southern and Central Zagros Mountains in Iran, pinches out towards northwest (Hinsch and Bretis, 2015;Kent, 2010). Other intermediate detachment horizons influence the structural style of Central Zagros in Iran (e.g., 15 Sherkati et al., 2006;Sepehr et al., 2006), but their behaviour is uncertain in NW Zagros due to limitations in outcrops and insufficient seismic profiles southwest of the Main Zagros Thrust. Some proposed detachment levels include Ordovician and Silurian shales (Aqrawi et al., 2010), Triassic-Jurassic anhydrites (Aqrawi et al., 2010;Hinsch and Bretis, 2015;Zebari, 2013;Zebari and Burberry, 2015), and Lower Miocene anhydrite (Aqrawi et al., 2010;Csontos et al., 2012;Jassim and Goff, 2006;Kent, 2010;Zebari and Burberry, 2015). 20 The exposed geological units within the High Folded Zone are limited to c. 5 km thick Upper Triassic to Recent rocks (Fig. 2;Jassim and Goff, 2006;Law et al., 2014). Most anticlines are made up of Cretaceous carbonate rocks, while Upper Triassic-Lower Cretaceous strata are only exposed in the core of some anticlines. The Tertiary clastic rocks are preserved within the adjacent synclines. Within the studied structures, the Upper Jurassic-Lower Cretaceous Chia Gara and Lower Cretaceous Sarmord formations only crop out in Bekhme and Zinta gorges and consist of medium to thick bedded marly limestone,  Unconsolidated Quaternary sediments in the study area consist of slope deposits, residual soil, alluvial fan deposits, and river terraces.
There is no agreement concerning the overall structural style of the NW Zagros in KRI. Several authors (Al-Qayim et al., 2012;Ameen, 1991;Fouad, 2014;Jassim and Goff, 2006;Numan, 1997) suggested that the Iraqi part of the Zagros Fold-Thrust Belt reveals a combination of both thin-and thick-skinned deformation. Partly relying on reflection seismic data, it was also suggested that contraction has been localized on inherited passive-margin normal faults in the basement, which were inverted during the late stage of deformation since c. 5 Ma (Abdulnaby et al., 2014;Burberry, 2015;Koshnaw et al., 2017).
The structural relief across the MFF (Fig. 2) is likely linked to blind thrusts in the basement (Al-Qayim et al., 2012;Ameen, 1991Ameen, , 1992Fouad, 2014;Koshnaw et al., 2017;Numan, 1997). The same linkage between structural relief and a regional 5 basement blind thrust is also documented in the Iranian Zagros (Blanc et al., 2003;Emami et al., 2010;Leturmy et al., 2010;Sherkati et al., 2006). Alternatively, Hinsch and Bretis (2015) argued that the structural relief in the hanging wall of the MFF is related to an underlying duplex structure that is linked to a stepped detachment horizon rooting in a Lower Paleozoic detachment in the internal parts of the orogen. The relief has been attributed to the accumulation of the Hormuz salt in the Iranian Zagros (McQuarrie, 2004). Even though the MFF is believed to be a major blind thrust in the basement (Berberian, 10 1995), it is usually mapped along the southwestern limb of the last high anticline where the Pila Spi limestones or the Bekhme and Aqra limestones crop out (Fouad, 2014;Jassim and Goff, 2006;Numan, 1997). Given that landforms in the vicinity of the MFF indicate ongoing tectonic deformation, we suspect that these blind faults might be active at present. Unfortunately, however, instrumental seismicity in the entire region is too diffusely distributed to be attributed to any particular faults (Jassim and Goff, 2006). 15 Structurally, this segment of Zagros Fold-Thrust Belt is dominated by NW-SE trending fault-related folds, the trend of folds changes to nearly E-W to the west of the Greater Zab River (Fig 2). The folds are usually S-verging and the related faults emerge to the surface within both Imbricated Zone and High Folded Zone, while they remain blind within the Foothill Zone (Fouad, 2014;Hinsch and Bretis, 2015).

Data and Methods 20
We calculated and analysed landscape indices from DEMs for the studied anticlines and built a landscape evolution model that simulates progressive uplift and erosion of the landscape. We also constructed structural cross-section across these anticlines based on literature data and our own field observations.

Geomorphic Indices
The present-day relief in the study area resulted from a competition between rock uplift triggered by horizontal contraction 25 and erosion destroying it. Parameters controlling these competing processes are the rate of tectonic accretion, rock erodibility and climate (Bishop, 2007;Burbank and Anderson, 2012).
In order to quantitatively analyse the landscape for the Harir, Perat and Akre anticlines (Figs. 2 and 4), we calculated hypsometric curves and determined three geomorphic indices: (i) hypsometric integral, (ii) surface roughness, and (iii) surface index. These are considered proxies for the relative maturity of a particular landscape. The hypsometric curve and integral 30 refer to the distribution of surface area of a landscape with respect to the elevation (Strahler, 1952). The surface roughness value is mainly sensitive to incision (Andreani et al., 2014;Andreani and Gloaguen, 2016;Pike and Wilson, 1971); the surface index is a measure for the amount of erosion. When referring to the results obtained by using this set of geomorphic indices, we colloquially refer to them as "landscape maturity" parameters.

Hypsometric Curve
The hypsometric curve for a basin is the frequency distribution of elevation of the watershed area below a given height 5 (Strahler, 1952). Convex-shaped hypsometric curves represent relatively youthful stages of the basin, s-shaped and concave curves refer to more mature and old stages (Strahler, 1952). Hypsometric curves are usually calculated for a specific drainage basin. In this study, we calculated the weighted mean of the hypsometric curves for basins with areas > 0.25 km 2 within each anticlinal ridge, weighted by the basin area within the anticline. We restricted our analyses to those basins where Upper Cretaceous carbonates are exposed (Fig. 4). This allowed us to make realistic comparisons between the three anticlines, 10 neglecting the differences in rock erodibility that arise when varying lithologies are included. Wind gaps and water gaps as well as the plunging crests of the anticlines were also excluded from the calculation.

Hypsometric Integral
The hypsometric integral (HI) is the ratio of area under the hypsometric curve (Strahler, 1952). It is used to highlight the erosional stage of a landscape with high values corresponding to less mature landscapes and low values indicating advanced 15 stages of erosion. The hypsometric integral is computed for a certain area by the following equation (Pike and Wilson, 1971): where h mean , h min and h max are the mean, minimum and maximum elevations [m] of the examined area.

Surface Roughness
The surface roughness (SR) measures how much an area deviates from being totally flat. It differentiates flat planar surfaces 20 with values close to 1 from irregular surfaces with higher values. It increases with the increase in incision by streams. The surface roughness is calculated using the following equation (Grohmann, 2004): where TS and FS are the areas [m 2 ] of the actual topographic surface and the corresponding projection of that surface onto a planar surface, respectively. 7

Surface Index
The surface index (SI; Andreani et al., 2014) combines elevations, hypsometric integral and surface roughness to map simultaneously preserved and eroded portions of an elevated landscape. It is calculated using equation 3 (Andreani and Gloaguen, 2016): where N HI , N h and N SR are the normalized hypsometric integral, elevations, and surface roughness values, respectively.
Elevated and poorly incised landscapes with high hypsometric integral and low surface roughness show positive surface index values. Highly dissected landscapes with a high surface roughness yield negative surface index values. This means that the surface index is also sensitive to elevation.

Digital Elevation Models 10
The geomorphic indices for this study were calculated from the 12 m resolution TanDEM-X DEM (Krieger et. al., 2007) obtained from the German Aerospace Center (DLR) and the 30 m resolution SRTM1 DEM (NASA JPL, 2013); these two inputs were used since different DEM inputs give slightly different geomorphic results (Andreani et al., 2014;Koukouvelas et al., 2018;Obaid and Allen, 2017). Geomorphic indices were calculated using both the TanDEM-X and the SRTM1 data.
However, the TanDEM-X data revealed numerous artefacts and voids, which made calculations unstable and results unreliable 15 (also see the comparison in the supplementary material). All results of the geomorphic indices and subsequent calculations presented in the following sections were calculated from a 100 * 100 cell (3 * 3 km) moving window on the 30 m resolution SRTM1 data. A larger moving window makes the obtained measurements smoother and vice versa. The size of the moving window must be chosen based on the scale of the target; here we targeted anticlines with wavelengths varying from 5 to 8 km.
A 3 km moving window covered almost an entire limb of an anticline. We also tested the method proposed by Pérez-Peña et 20 al. (2009) in order to account for the neighbouring cells in the calculation of the geomorphic indices. Rather than using a moving window, this approach uses a spatial autocorrelation of neighbouring cells and maps clusters of high and low values of indices using Moran's I statistic (Moran, 1950) and Gi * statistics (Ord and Getis, 1995). We have tested the same method here by calculating the HI for a 500 * 500 m grid of the SRTM data. We applied a hot spot analysis using Gi * statistics with a distance of 1.5 km to define neighbour cells. Then, we resampled the HI map calculated from a 100 * 100 cell (3 * 3 km) 25 moving window to 500 * 500 m grid from SRTM data. The calculations were performed using the focal and zonal toolsets in ESRI ArcGIS 10.4 software. In addition, the SRTM DEMs with 30 m resolution were used to extract topographic profiles, drainage networks, watersheds, stream slopes, and upstream drainage areas wherever required.

Modelling Landscape Evolution
We built a landscape evolution model to quantify the time difference in between the maturity level of the Akre and Harir anticlines by comparing the geomorphic indices of the evolved landscape with those of both anticlines based upon the opensource Landlab toolkit (Hobley et al., 2017; http://landlab.github.io).
We used two components in our model: one simulating erosion due to fluvial action and another simulating sediment transport 5 along slopes due to hillslope diffusion processes. Chen et al. (2014) showed that consideration of only these two components is sufficient for many landscapes but cannot model fluvial sedimentation. However, from field observations and from satellite imagery, we know that no significant fluvial sedimentation takes place on the slopes of the analysed anticlines. On slopes of anticline flanks, the detachment-limited erosion due to the fluvial system tends to be the dominant process (Howard, 1994).
To detect changes in the landscape due to fluvial erosion through time, we applied the commonly accepted idea that the rate 10 of stream incision is directly proportional to the hydraulic shear stress of a stream (Braun and Willett, 2013). Consequently, we used the stream power incision law (Sklar and Dietrich, 1998;Whipple and Tucker, 1999): where ∂z/∂t is the erosion rate [myr -1 ]; K is an erodibility coefficient [yr -1 m (1-2m) ] that encompasses the influence of climate, lithology, and sediment transport processes; A is the upstream drainage area [m 2 ] and is typically taken as a proxy for discharge 15 (Wobus et al., 2006); S = ∂z/∂x is the local channel slope [m/m]; z is the elevation [m]; and m and n are the area and slope exponents, respectively. The stream power incision law (Eq. 4) is derived since the upstream drainage area A scales with channel discharge and channel width. The magnitude of the sediment flux in the channel is assumed to equal unity in the standard detachment-limited stream power model (Perron, 2017;Whipple, 2002). In the model, an incision threshold (C) was included, below which no incision occurs (Hobley et al., 2017). 20 To account for the provision of sediment due to hillslope diffusion processes from slopes outside the river system, we used the hillslope diffusion equation (Culling, 1963;Tucker and Bras, 1998): where Kd is the diffusivity coefficient [m 2 yr -1 ], z is the elevation [m], and ∇ 2 is the Laplace operator, i.e. the divergence of the gradient. 25 Finally, the overall evolution of the landscape in different time steps was calculated as the uplift rate subtracted by the changes due to both fluvial erosion and the hillslope diffusion (Temme et al., 2017): where U is the uplift rate [myr -1 ].
A DEM raster grid of the present-day Harir Anticline and the surrounding basins (Fig. 5a) served as model input. The advantage of using Harir Anticline was that the evolved drainage network overprinted the pre-existing one. The boundary conditions were set as closed on all sides except in pre-existing outlets in the input grid. The basins surrounding Harir Anticline were also included in the input grid to minimize the effect of the boundary conditions on the Harir Anticline itself. In the input raster grid, a flow route of each cell was connected with neighbouring cells both diagonally and orthogonally. This means that each 5 cell had the possibility to be linked with eight surrounding cells across its sides and corners (Hobley et al., 2017;Tucker et al., 2016).
Concerning the parameter used in the model, the value of m/n, n, and K were found following the methodology described by Harel et al. (2016), Mudd et al. (2014), and Perron and Royden (2013), and by comparison with data from Harel et al. (2016).
The value of m/n was found by plotting the elevation against X (elevation-X plot) for streams in the input grid (Fig. 5a), where 10 X is found following the equation described by (Perron and Royden, 2013): where A 0 is the reference drainage area [m 2 ] of 166160 m 2 and x is the horizontal upstream distance [m]. In this approach, we ascribed values for m/n range from zero to one, and X was calculated for each time from Eq. 7. The value of m/n with maximum regression (R 2 ) value in the elevation-X plot was taken as the best-fitting value, which was 0.41 in our case for the present-15 day Harir Anticline's drainages (Fig. 5b). This value of m/n is located within the theoretically predicted values of m/n, which ranges from 0.3 to 0.7, based on the stream power incision model (Kwang and Parker, 2017;Temme et al., 2017;Whipple and Tucker, 1999).
In the model, n = 1.7 and K = 3.0E-6 yr -1 m -0.4 were used; these values were estimated as mean of K and n in Harel et al. (2016) for those areas that are comparable with our study area in aspect of lithology, climate, and precipitation. The value of m was 20 0.7. We used an incision threshold of C = 1.0E-5myr -1 , which is widely adopted for erosion of an upland landscape (Hobley et al., 2017). A present-day annual mean precipitation of c. 0.7 myr -1 was used through the time due to the lack of nearby paleoclimate data with good quality. The average of the modelled precipitation anomaly data for Lake Van in SE Turkey (200 km to the NNW of the studied anticlines) is close to zero (Stockhecke et al. 2016). The current elevation of the Bekhme and Aqra formations in the crest of Harir Anticline is about 1500 m above sea level. Above that, 2070 m of Upper Cretaceous-25 Miocene units (Law et al., 2014) and 300 m of Upper Miocene Lower Bakhtiari were exhumed before exposure of the Bekhme and Aqra formations. If we consider that the Lower Bakhtiari have been deposited close to sea level before onset of the MFF c. 5 Ma, there would be 3872 m of rock uplift at a rate of ~0.0007 myr -1 , which was used in the model. Since soil (regolith) is rare and very thin when present on the slopes, a low diffusivity coefficient of Kd = 0.001 m 2 yr -1 was used (Fernandes and Dietrich, 1997). 30 There are minor variations in lithology between the three anticlines (they consist of a thick pile of dominantly Cretaceous carbonate) and no variation in climate can be expected on such a relatively local scale. Therefore, no significant variances are expected in the used parameters. Lastly, the parameters were calibrated by comparing the nature of the evolved landscape to other anticlines within the High Folded Zone that are cored by Cretaceous carbonates and more mature than the Harir Anticline to evaluate how realistic the evolved landscape is.

Landscape Maturity 5
The three studied anticlines are composed of the raised Cretaceous carbonates in their crests, whereas the Tertiary clastic rocks have been denuded, but conserved in the adjacent synclines. The three anticlines are dissected by rivers that form water gaps across them. Bekhme and Zinta gorges cut the Perat and Akre anticlines, respectively. We also observed wind gaps, such as those in the NW end of Harir Anticline (Zebari and Burberry, 2015). Therefore, neither the location of these water and wind gaps nor the plunging tips of anticlines have been considered in interpreting the geomorphic indices as proxies for relative 10 landscape maturity.
The anticlines reach up to c. 1500 m asl. The minimum altitude is c. 700 in the adjacent synclines and c. 400 m in the Greater Zab river course. The hypsometric curves for the three anticlines are presented in Fig. 6. Harir Anticline's curve is more convex, and its shape is close to the youthful stage of Strahler's diagram (Ohmori, 1993;Strahler, 1952) with 78% of the area above the mean elevation, while Akre Anticline is less convex and close to a mature stage with only 50% of the area above 15 the mean elevation. Perat Anticline's values are located in between and closer to the Harir Anticline curve with 64% of its area above the mean elevation.
The three calculated geomorphic indices (HI, SR, and SI) seem to be substantially influenced by the local structure, and wind and water gaps (Fig. 7). In our calculation, the results of geomorphic indices change with changing the size of the moving window and the resolution of the input data (see supplementary material). This was also detected by Andreani et al. (2014) and Obaid and Allen (2017). Andreani et al. (2014) found that the DEM resolution does not affect the hypsometric integral, but it affects the surface 10 roughness, while the size of the moving window affects both hypsometric integral and surface roughness. The results become smoother with an increasing size of the moving window. Here, we found that it is reasonable to use a 100 * 100 cell (3 * 3 km) moving window, which covers approximately one limb of the anticline with 6-7 km width. It therefore highlights the desirable signal. The cluster map for the HI Gi* statistics was calculated following the approach by Pérez-Peña et al. (2009) for the 500 * 500 m grid. We obtained results similar to the HI map calculated from a 100 * 100 cell (3 * 3 km) moving 15 window and resampled to 500 * 500 m grid in terms of highlighting the cluster of high and low HI values (Fig. S19 a and b in supplementary material). This comparison proves that our method is equally applicable and valid. We therefore ran all analyses based on the 100 * 100 cell moving window as described above.
According to the hypsometric curves and the geomorphic indices, we found that there is a measurable difference in landscape maturity between the three anticlines. We classified our anticlines as relatively mature (Akre Anticline), moderately mature 20 (Perat Anticline), and less mature (Harir Anticline). The difference in the maturity level must be due to a difference in one or more of the factors tectonics, climate, or rock erodibility. No variation in the climate is expected for the scale of the studied area, therefore its impact on the landscape maturity can be neglected. The three anticlines show essentially the same lithology (Figs. 2 and 3), i.e. similar lateral rock type and erodibility. Thus, the only factors that may vary along the anticlines are uplift rate or onset time of the uplift. This can be interpreted with one of the following scenarios: either the anticlines started to uplift 25 in the order 1) Akre, 2) Perat and 3) Harir from west to east, or all of them started at the same time but with different rates. In the latter case, the uplift rate would have been highest at Akre and lowest at Harir, exposing Akre to erosion earlier than Harir.

Landscape Evolution Model
The aged landscape from the model run (Figs. 8 and 9) is the result of fluvial erosion and hillslope diffusion on the one hand, and uplift due to folding on the other hand. In the landscape modelling, various simulations with different parameters and time 30 spans were performed (supplementary material S24-S27). Harir Anticline was used as an input model and the landscape evolution model was run for a time span of 10 kyr up to 100 kyr and then it was run for a time span of 20 kyr. The evolving drainage system overprints the pre-existing one in the input and gradually becomes more deeply incised from the anticline flanks carving toward its core (Fig. 8). Harir Anticline is a box-shaped anticline with a wide and flat crest area. With ongoing incision towards the core of the anticline, this plain crest narrowed gradually and finally became a sharp ridge that divided the drainage basins of the SW flank from those of the NE.
We compared the hypsometric curves of the model outputs to the present-day curves of the anticlines. The hypsometric curves were calculated first as total weighted mean for basins within the anticline and later calculated for the entire anticlinal ridges 5 ( Fig. 9). Statistically (with minimum RMS), the hypsometric curve of Harir Anticline was closest to the present-day Perat Anticline after 100 kyr and 80 kyr of erosion when using total weighted mean and entire anticlinal ridge, respectively, in calculation of the hypsometric curves. The output curve after 160 kyr and 200 kyr matched best with present-day Akre Anticline, when using total weighted mean and entire anticlinal ridge, respectively, in calculation of the hypsometric curves.
We conclude that it will take Harir Anticline roughly about 80 to 100 kyr to reach the maturity level of Perat Anticline and 10 about 160 to 200 kyr to reach the level of Akre Anticline based on our model and comparison of the hypsometric curves if the uplift rates of the three anticlines were the same. The other possibility is that the anticlines started to grow at the same time but with different uplift rates. In this case, it is not possible to find the difference in uplift rates via our landscape modelling.
Since the factors that control geomorphology (lithology, structural setting, and climate) were similar for all three anticlines, and under the assumption of constant growth and erosional conditions, we infer that uplift of Akre and Perat anticlines started 15 respectively 160-200 kyr and 80-100 kyr before Harir started to grow if their uplift rates were the same.

Geometry of the Studied Anticlines
The structural cross-sections for the three anticlines ( Fig. 10) constructed from field data and literature (Syan, 2014) shows that the anticlines are box-shaped with broad crests. They are asymmetrical verging toward the SW with a nearly vertical or overturned forelimb. The three anticlines are thrust-related, in accordance with published studies of the area (Csontos et al. 20 2012), and have a thrust in their forelimb. Perat Anticline additionally exhibits a back thrust in its NE limb. The shortening across the three anticlines was calculated using line-length balancing. We found 26%, 28% and 29% shortening for Harir, Perat and Akre anticlines, respectively, and conclude that there is no significant variation. A difference in the fold amplitude between the three anticlines can be discerned in the cross-sections. The amplitude of Perat Anticline is higher than that of both Harir and Akre anticlines. Another difference concerns the thickness of the Upper Cretaceous -Middle Eocene clastic 25 succession, which dwindles toward the Akre Anticline. The whole Miocene succession, however, is the thickest in the Akre Anticline. Both the Upper Cretaceous -Middle Eocene and Miocene successions consist of highly erodible clastic rocks, and the thicker Miocene succession in the Akre Anticline counterbalances the thinner Upper Cretaceous -Middle Eocene succession. Therefore, it is not expected that these variations in the structural geometry and stratigraphic thickness have a great impact on the variation of the landscape maturity in the three anticlines and the landscape model.

Landscape Maturity and Modelling
Any relative change in the base level induced by tectonics or climate leads to the change of erosion rates. A relict landscape survives when its uplift is not completely counterbalanced by erosion (Andreani and Gloaguen, 2016;Burbank and Anderson, 2012;Pérez-Peña et al., 2015). The relative relict landscape and its distribution on these three anticlines atop the MFF reveal 5 clues about underlying tectonics since there is no significant variation in climate and lithology.
Within the three studied anticlines, the geomorphic indices effectively highlighted their incision. The surface index which combines both hypsometric integral and surface roughness, sets apart relict landscapes of positive and high values from transient landscapes of negative values that are preferentially incised (Andreani et al., 2014). There is a notable relative deviation in areas where anticlines are crossed by rivers e.g. Bekhme and Zinta gorges, which show a high surface roughness. 10 Also, variations in surface index are found in comparable areas in the three anticlines. Focussing on the crest of the anticlines, we observe that Harir Anticline shows higher values than the two others. The lowest values are found in Akre Anticline. Harir has low incision at elevated surfaces while Akre has more incised uneven landscape, because erosion has worked deeper into the core of the anticline. This can also be inferred from the valley shapes. We observe a tight V-shape valleys in the flanks of Harir and a open V-shape valleys in Akre (Figs. 11a and 11b). We relate this difference to the tectonic uplift and to the effects 15 of longer erosion acting on Akre. The observation can be interpreted with one of the following two premises: either the anticlines started to uplift successively (first Akre, then Perat, and finally Harir), or all of them started at the same time but with different uplift and exhumation rates (Akre the fastest, Harir the slowest). In other words, the Cretaceous carbonates in Harir Anticline were exposed to the erosion later than in Akre Anticline and, consequently, incised less.
The current landscape of these anticlines exposes Cretaceous carbonates of the Qamchuqa, Bekhme and Aqra formations, 20 which became exposed to erosion only after unroofing of the entire Palaeogene -Neogene succession. The Upper Miocene-Pliocene Bakhtiari Group, which is the youngest stratigraphic unit in the area, is affected by folding, as observed from growth strata (Csontos et al., 2012). This has also been observed in the Upper Bakhtiari (Pliocene-Pleistocene) close to the MFF (Koshnaw et al., 2017). In between Bekhme and Aqra and the Upper Bakhtiari formations, 2.37 km of the Upper Cretaceous -Miocene clastic rocks interbedded with thin units of limestone (Law et al., 2014) have been exhumed due to successive rock 25 uplift in the crest of the studied anticlines, triggered by shortening and erosion. They are only preserved in the adjacent synclines. The Cretaceous carbonates themselves have been exposed in the crests of Akre, Perat, and Harir anticlines for c. 0.9 km above the level of the other exhumed units. Based on the thickness, the amount of the exposed Cretaceous carbonate makes c. 28 % of the total exhumed and exposed thickness in the crest of the anticlines. Therefore, with both scenarios (different uplift time or different uplift rate) and with assumption of constant (linear) rock uplift rate through time, the Cretaceous 30 carbonate in Harir Anticline was exposed to erosion later than in the Akre Anticline (Fig. 12a).
The steeper valley flanks in Harir Anticline compared to those of Akre also support higher uplift rates of the Harir Anticline.
Furthermore, the relationship between stream slope and upstream drainage area at any given point of streams in the Harir Anticline is positive (Fig. 12b). This means that the streams have a convex shape and the streams' segments with steeper slopes are still located in the flanks of the anticline. In the Akre Anticline, this relationship is negative (Fig. 12b), which means that the streams have a concave shape and the segments with steeper slopes have migrated toward the core of the anticline. This implies that tectonic activity in the Harir Anticline is younger than in the Akre Anticline. Therefore, the premise of having Harir Anticline starting its uplift later than Akre Anticline is most likely. This is our preferred scenario in the model for the 5 successive tectonic evolution of the study area presented below.
The geomorphic indices have been widely used to assess the landscape maturity and, subsequently, active tectonics (Andreani and Gloaguen, 2016;Mahmood and Gloaguen, 2012;Pérez-Peña et al., 2009). The challenging aspect of landscape maturity modelling is to obtain an absolute quantification of tectonics and the relevant time spans. The same holds true with using a landscape evolution model to estimate the relative time difference between two landscapes, because it is difficult to compare 10 two landscapes in terms of maturity by absolute means. The results of the landscape modelling approach yielded a numerically derived estimate on the relative age difference between the studied anticlines but without absolute growth ages.
In the model, various parameters and two well-known landscape evolution equations for the fluvial erosion and hillslope diffusion were used, but in general it is impossible to mimic nature perfectly. The relative time difference of landscape evolution of these anticlines was estimated from the model assuming that the climate has not changed much during the 15 evolution of the landscape since there was not much variation in the precipitation based on the modelled data (Stockhecke et al. 2016)and for the sake of simplicity, admiting that climatic change has a significant impact on the landscape. In addition, neither rock fall nor karstification were included in the model for simplicity. Field observations suggest that karstification does not have a significant impact on the landscape. Overall, the evolved landscape from the model seems to be plausible in comparison with the other anticlines that surround Harir Anticline, and the landscape models are more mature with respect to 20 the developed topography and to the overall drainage patterns.

Structural Style and Regional Tectonics
An orogenic bend is depicted in the area where the trend of structures changes across the Greater Zab River from NW-SE at the eastern side of the river to nearly E-W at its western side. The course of the Greater Zab River is suggested to overlie a NE-trending transversal basement fault with right-lateral displacement (Ameen, 1992;Burberry, 2015;Jassim and Goff, 2006;25 Omar, 2005) evidenced by an offset of the High Folded Zone propagation foreland-ward. At the eastern side of the river the deformation has propagated for about 25 km further than on its western side (Figs. 1 and 2). The origin of this fault reaches back to the Late Proterozoic tectonic history of the Arabian Plate. This fault has been reactivated later in subsequent tectonic events (Ameen, 1992;Aqrawi et al., 2010;Burberry, 2015;Jassim and Goff, 2006). This can also be noticed in the thickness of the sedimentary cover, which is thinner to the west of the Greater Zab River (Ameen, 1992;Zebari and Burberry, 2015). 30 This change in thickness is attributed to a series of uplift events and erosional/non-depositional gaps during the Mesozoic (Ameen, 1992;Aqrawi et al., 2010), which in turn may have influenced the foreland-ward propagation of deformation (Marshak and Wilkerson, 1992). The deference in propagation of deformation may also be due to the convergence being accommodated differently across the curved fold-thrust belt (Csontos et al., 2012) from NW-SE in the eastern side where the convergence is accommodated through belt-normal slip and right-lateral strike-slip across the eastern NW-SE trending part, and by belt-normal slip across the western E-W trending segment (Reilinger et al., 2006). Zebari and Burberry (2015) found that anticlines to the east of the Greater Zab River (Harir, Shakrok and Safin anticlines) demonstrate pronounced NW-ward propagation based on their geomorphic criteria, and the start point of the NW-ward 5 propagation of the Harir Anticline is close to its SE end. This implies that progressing uplift in the hanging wall of the MFF was not gradually continuing from the Akre Anticline towards the Perat Anticline and further SE-ward to the Harir Anticline.
The uplift progress is probably rather partitioned into segments along the belt. In addition, other anticlines to the south (Safin Anticline) and to the southwest (Shakrok Anticline) of Harir Anticline are more mature than Harir Anticline itself, based on their hypsometric curves (Fig. 13) and geomorphic indices (supplementary material; S1-S18), implying that the foreland-ward 10 propagation of the deformation was also out-of-sequence in this part of the High Folded Zone. This has also been noted in the Foothill Zone based on thermochronological dating (Koshnaw et al., 2017) and landscape maturity (Obaid and Allen, 2017).
Thus, the most plausible scenario is that deformation in the Harir segment started sometime after that in Akre segment (160-200 kyr according to our landscape evolution modelling). Harir Anticline uplift would also postdate Perat Anticline uplift (80-100 kyr) to the west and the onset of Safin and Shakrok anticlines to the south and southeast, which are not included in the 15 model (Fig. 14). As discussed by Csontos et al., (2012), the fold relay corresponds to the change in strain partitioning and rotation of the horizontal stress direction from NE-SW to N-S in Late Pliocene (Navabpour et al., 2008;Navabpour and Barrier, 2012). During the latest stage of the N-S convergence, a right-lateral shear and superposed folding along the NW-SE trending anticlines (Csontos et al., 2012) can be observed from the relay of the Shakrok, Harir and Perat anticlines (Figs. 2 and 14).
Applying this concept requires a comprehensive paleostress analysis investigation especially within these studied anticlines, 20 which is beyond the scope of this paper.

Conclusions
The geomorphic indices used in this study allowed us to quantitatively differentiate between variably degraded landforms in the frontal Zagros Mountain of NE Iraq. This area is characterized by active folding due to ongoing convergence between the Eurasian and Arabian plates. Three active thrust-related anticlines that are aligned along-strike the MFF were studied in detail. 25 While the Akre Anticline shows deeply incised valleys indicative of advanced erosion, the Harir and Perat anticlines have relatively smooth surfaces and show younger landscape than Akre. We related this difference to the underlying tectonics. This can be interpreted with one of the following concepts: either anticlinal growth started at different times or all of them started to grow at the same time, but with different surface uplift and exhumation rates.
A comparison of the geomorphic indices values of the model output with those of the present-day Akre Anticline topography 30 revealed that it will take Harir Anticline about 160-200 kyr to reach the maturity level of today's Akre Anticline, and about 80-100 kyr to reach the maturity level of the Perat Anticline assuming constant uplift rates along the three anticlines. Due to similarity in the lithology, structural setting and climate along the three anticlines and by assuming constant growth and erosion conditions, we infer that Akre Anticline started to grow 160-200 kyr before Harir Anticline. The onset of growth of Perat Anticline lies in between that of Harir and Akre anticlines. A NW-ward propagation of Harir Anticline itself implies that the uplift has been independent within different segments rather than having been continuous from the NW to the SE. Our method of estimating relative age differences in variously degraded anticlines can be applied to many other anticlines along the MFF 5 and could eventually develop into a model of the temporal evolution of this fold and thrust belt.

Acknowledgments
The German Academic Exchange Service (DAAD) is acknowledged for providing a scholarship (Research Grants -Doctoral Programmes in Germany, 2016/17;57214224) to the first author to conduct his PhD research in Germany. The authors express 10 their gratitude to the German Research Foundation (DFG) project no. 393274947 for providing financial research support. The German Aerospace Agency (DLR) and Mr. Thomas Busche in particular are thanked for providing TanDEM-X digital elevation models. We also thank two anonymous peers for their partly very thorough reviews that helped us to further clarify and improve numerous aspects of our work.  ( S a n a n d a j -S i r j a n Z o n e )  A k re A n ti c li n e Pe ri s A n ti cl in e    (300) U. Fars (980) L. Fars (262) Pila Spi (100) Gercus (50) Khurmala (39) Kolosh (277) Tanjero (192) Shiranish (172) Aqra (128) Bekhme (324) Qamchuqa (231) Sarmord (335) Chia Gara (117)

5
A k re A n ti c li n e P er is A n ti cl in e Ch ia Ga ra An tic lin e S h ir in A n t ic li n e P e ra t A n t ic li n e B ij e e l A n t ic li n e S a r t a A n t i c l i n e A k re A n ti c li n e P er is A n ti cl in e Ch ia Ga ra An tic lin e S h ir in A n t ic li n e P e ra t A n t ic li n e B r a d o s t A n t i c l i n e S a n A n t i c l i n e S h a k r o k A n t i c l i n e K o r e k A n t i c l i n e H a n d r e n A n t ic li n e P e l e w a n A n t i c l i n e M a k o o k A n t i c l i n e R a n y a A n t i c l i n e