Late Pleistocene‐Holocene Slip Rates in the Northwestern Zagros Mountains (Kurdistan Region of Iraq) Derived From Luminescence Dating of River Terraces and Structural Modeling

A significant amount of the ongoing shortening between the Eurasian and Arabian plates is accommodated within the Zagros Fold‐Thrust Belt. However, the spatial and temporal distribution of active shortening within the belt, especially in its NW part, is not yet well constrained. We determined depositional ages of uplifted river terraces crossing the belt along the Greater Zab River using luminescence dating. Kinematic modeling of the fault‐related fold belt was then used to calculate long‐term slip rates during the Late Pleistocene to Holocene. Our results provide new insight into the rates of active faulting and folding in the area. The Zagros Mountain Front Fault accommodates about 1.46 ± 0.60 mm a−1 of slip, while a more external basement fault further to the SW accommodates less than 0.41 ± 0.16 mm a−1. Horizontal slip rates related to detachment folding of two anticlines within the Zagros Foothills are 0.40 ± 0.10 and 1.24 ± 0.36 mm a−1. Basement thrusting and thickening of the crust are restricted to the NE part of the Zagros belt. This is also reflected in the regional topography and in the distribution of uplifted terraces. In the southwestern part, the deformation is limited mainly to folding and thrusting of the sedimentary cover above a Triassic basal detachment. In the NE, deformation is associated with slip on basement thrusts. Our study sheds light on the distribution of shortening in the Zagros Mountains and helps to understand the regional tectonic system. Our results may be the foundation for a better seismic hazard assessment of the entire area.

SW of the Turkish-Iranian Plateau (Allen et al., 2013;Berberian, 1995;Nissen et al., 2011;Talebian & Jackson, 2004). Large earthquakes in the Zagros appear to occur within the upper part of the basement and at the base of the sedimentary cover at depths of 8-15 km (Nissen et al., 2011;Talebian & Jackson, 2004). The instrumental seismicity within the northwestern part of the Zagros Belt in Kurdistan Region of Iraq (hereinafter NW Zagros) is lower and more disperse compared to the Iranian Zagros ( Figure 1; Dziewonski et al., 1981;Ekström et al., 2012). Also, geodetic data in the NW Zagros are sparse (Khorrami et al., 2019;Reilinger et al., 2006). Therefore, spatial and temporal constraints of the present-day shortening distribution  (Zebari et al., 2019). Rivers and their names are in blue. The overview map shows tectonic outlines of the Arabia-Eurasia collision, Arabian Plate motion with respect to fixed Eurasia after Reilinger et al. (2006), and earthquake focal mechanisms with Mw >4.6 from 1976 to 2019 using the Global CMT database (Dziewonski et al., 1981;Ekström et al., 2012) are far from being well understood, and no data on the distribution of active shortening are available for the NW Zagros.
Several rivers drain the Zagros towards the Persian Gulf. Fluvial terraces along the banks of the Mand River have been successfully used to study active tectonics and Late Pleistocene slip rates in the Fars Arc (Collignon et al., 2019;Oveisi et al., 2009). Within the NW Zagros, Quaternary river terraces are abundant along the banks of the Tigris River and its tributaries. Here, we have used the river terraces of the Greater Zab River, which crosscuts several tectonic units of the NW Zagros, to estimate long-term fault slip rates in this part of the belt for the first time. First, we have mapped a system of staircase terraces along the banks of the river and put them into relation with the fault-related folds traversed by the river. We subsequently sampled sediment material from numerous of these terraces. These samples were dated using a combination of infrared stimulated luminescence (IRSL) and post-infrared stimulated luminescence (pIR) dating to estimate their relative vertical uplift rates (i.e., without an a priori assumption as to whether river incision occurred with or without surface uplift). The relative terrace uplift rates were then integrated with kinematic models for all traversed fault-related structures along the studied transect, including modeling of particle trajectories during progressive slip. This approach led us to conclude that terrace uplift is largely induced by surface uplift triggered by fault-related folding. We subsequently converted uplift rates into horizontal slip rates on the faults within the corresponding fault-related structures. Our results provide long-term slip rate estimates for the Late Pleistocene to Holocene. We were also able to detect a regional uplift rate in the area unrelated to fault-related folding above shallow detachments, but rather caused by crustal thickening and/ or basement thrusting.

Geological Background
The Zagros Fold-Thrust Belt extends for c. 2000 km from the Hormuz Strait in southern Iran NW-ward to the Kurdistan Region of Iraq and SE Turkey (Agard et al., 2011;Berberian, 1995;Mouthereau et al., 2012). It results from the convergence and subsequent continent-continent collision between the Arabian and Eurasian plates during Oligocene-Early Miocene following the obduction of the ophiolite sequences on Arabia's margin in the Late Cretaceous and the subduction of the Neotethys oceanic crust beneath the Eurasia Plate (Agard et al., 2011;Khadivi et al., 2012;Koshnaw et al., 2018;Mouthereau et al., 2012). Since the onset of collision, the deformation front has propagated 250-350 km SW-ward during the Neogene, involving the northeastern margin of the Arabian Plate into a largely NW-SE-trending foreland fold-thrust belt (Agard et al., 2011;Mouthereau, 2011;Mouthereau et al., 2007).
The Zagros Belt is subdivided into several NW-trending, belt-parallel morphotectonic zones that are bounded by major faults (Figure 1). In the NW Zagros these are from NE to SW: (a) the Zagros Suture Zone, which consists of the Tethyan allochthonous nappes of Mesozoic-Cenozoic ages; (b) the Zagros Imbricated Zone, which consists of thrusted anticlines and overriding imbricated structures; (c) the High Folded Zone that is dominated by closely spaced anticlines of high amplitude with Paleogene or Mesozoic carbonates exposed in their cores; and (d) the Foothill Zone with low amplitude anticlines that have a Neogene core and wide synclines of thick Miocene-Quaternary molasse (Figure 1; Jassim & Goff, 2006 berian, 1995;Jassim & Goff, 2006).
Structurally, this segment of the Zagros Fold-Thrust Belt is dominated by NW-SE trending fault-related folds. West of the Greater Zab River, the fold trend changes to nearly E-W (Figure 2a; Csontos et al., 2012;Le Garzic et al., 2019;Tozer et al., 2019;Zebari et al., 2020). The folds are usually S-verging and the related faults emerge at the surface within both the Imbricated Zone and the High Folded Zone, while they remain blind within the Foothill Zone.
In the NW Zagros, the sedimentary succession reaches a thickness of 8-13 km ( Figure 2b). These deposits consist of various competent and incompetent rock successions separated by detachment horizons above the crystalline basement of the Arabian Plate (Aqrawi et al., 2010;Jassim & Goff, 2006). The Infra-Cambrian Hormuz salt is unlikely to be present in the NW Zagros (Talbot & Alavi, 1996). Cambrian to Ordovician units consist mainly of epicontinental clastic deposits followed by Devonian thin transitional sequences of marine sedimentation and Carboniferous-Permian shales and carbonates. The Mesozoic sequence starts with the continental passive margin setting of Arabia in the Neotethys. Thick platform carbonates and evaporites with some open marine shale intervals make up most of these units. During the Cretaceous, platform and basinal carbonates became more dominant, and radiolarian cherts of the Qulqula formation were  Zebari et al., 2019;Zebari & Burberry, 2015). For location see Figure 1. (b) Stratigraphic column in the study area with scaled formation thicknesses (Zebari et al., 2020). Dashed red line denotes thicknesses used for constructing cross-sections, with a thickness of (c) 12 km in total. deposited in what later became the Suture Zone. In the Late Cretaceous, a foredeep basin formed in front of the obducting ophiolites and nappes of the Suture Zone. This basin was characterized by carbonates that changed upward to siliciclastic sediments. Paleogene units are clastic sediments that record obduction, with intervals of reef and lagoon carbonates. The High Folded Zone was mainly raised during the Oligocene. A structural step separated these units in the High Folded Zone from the Foothill Zone, which still produced shallow marine carbonates. The Zagros foreland basin has developed during the Miocene with the deposition of marine carbonates, evaporites, and clastic sediments that changed to continental fluvial and alluvial successions (Aqrawi et al., 2010;Jassim & Goff, 2006).
Mechanically, the Paleozoic succession consists of competent units of carbonates and sandstones with thin shale intervals. The Paleozoic units reach a thickness of more than 5 km and are linked with the crystalline basement (Le Garzic et al., 2019;Tozer et al., 2019;Zebari et al., 2020). The Triassic is about 1.8 km thick and the shales and evaporites of the Lower Triassic are decoupled from the Paleozoic units. This level is considered as the basal detachment in the NW Zagros (Le Garzic et al., 2019;Tozer et al., 2019;Zebari et al., 2020). Mechanically weak layers of evaporites and shales continue within the Upper Triassic. Assigning Triassic as the basal detachment is consistent with the wavelengths of folds (Le Garzic et al., 2019) and with an area balancing approach of calculating the depth of the detachment (Zebari & Burberry, 2015). Several incompetent units within the Jurassic and Lower Cretaceous can act as minor detachment levels and show some changes in bedding dip (Le Garzic et al., 2019;Sapin et al., 2017). However, they do not control the overall geometry of the constrained structural cross-sections (Le Garzic et al., 2019;Tozer et al., 2019;Zebari et al., 2020). The Lower-Upper Cretaceous units consists of thick competent carbonates. The Upper Cretaceous-Paleogene succession consists of incompetent units with carbonate intervals showing minor disharmonic folding. Miocene Lower Fars evaporites are considered as a weak layer. These units do not influence the style of folding of the younger units (Csontos et al., 2012;Le Garzic et al., 2019;Tozer et al., 2019;Zebari & Burberry, 2015).

Regional Cross-Section
We constructed a balanced structural cross section across the belt and along the Greater Zab River to study the deformation geometry of the NW Zagros and the style of basement thrusting (Figure 3). Field measurements were the main input for this section along with an interpreted seismic line (Csontos et al., 2012)

SW NE
The cross section was constructed using the Move Software (version 2019.1) from Petroleum Experts. We used this cross section to interpret river terrace uplift in terms of slip on faults.
The section extends from the Zagros Suture Zone to the undeformed foreland. It covers the entire Imbricated and High Folded Zones and the NE part of the Foothill Zone (Figures 2a and 3). Here, the Zagros Suture Zone has propagated further southwest than in adjacent areas ( Figure 1). The geometry and the length of the footwall below the Main Zagros Fault are not constrained. The total width of both the Imbricated and High Folded Zones within the transect from the Main Zagros Fault to the Mountain Front Flexure is less than 25 km. Along the considered transect, only one anticline (Bradost) is present within the Imbricated Zone and one anticline (Perat) within the High Folded Zone.
The Bradost Anticline is asymmetrical. It verges towards SW with a gentle back limb and a nearly vertical forelimb. It has a tight interlimb angle. The anticline is cut by an array of reverse faults in its forelimb and by a minor backthrust. These faults nucleated from the Triassic detachment level; their displacement is related to the so-called High Zagros Fault. The displacements on these faults range from 100 to 1,200 m on the main fore-thrust. The Perat Anticline is box-shaped and asymmetrical with a steeper backlimb. It has an open interlimb angle and is associated with a forethrust and a backthrust.  (Berberian, 1995). To the SE of the Greater Zab River, the Mountain Front Fault emerges from the Triassic detachment to the SW of Bina Bawi Anticline ( Figure 2a; Le Garzic et al., 2019) along strike of the Pirmam and Sarta anticlines. The anticlines to the SE of the Greater Zab Rivers show prominent growth towards the NW (Zebari & Burberry, 2015). These fault-related anticlines reveal contrasting maturity, implying they started growing at different times (Zebari et al., 2019). Therefore, we expect that the segment investigated here extends NW-ward along the SW limb of Pirmam and Sarta anticlines with less displacement based on the structural relief of the hanging wall. To the SW of the Perat Anticline, Safin (only it's NW segment is crossed by our section), Sarta, and Girda Baski anticlines are detachment folds above the Triassic detachment level. Generally, the anticlines in the Foothill Zone are widely spaced and have low amplitudes with gentle interlimb angles.
The amount of shortening along the entire transect of Figure 3 is c. 7.8 km. This is less than in the surrounding transects (34, 18.2, 19, and 24.8 km by Le Garzic et al., 2019;Koshnaw et al., 2020;Tozer et al., 2019;Zebari et al., 2020, respectively). Along our transect the shortening is likely a minimum estimate due to underestimations of the amount of displacement on the Main Zagros Fault, which has been displaced more foreland-ward than in the surrounding areas ( Figure 1). The structural step across the Mountain Front Flexure is related to the displacement on a basement thrust. The displacement on this thrust was accumulated within the frontal anticline of the Mountain Front Flexure (Perat Anticline) first, and then within the Foothill Zone by detachment folding (Zebari et al., 2020)

Field Work and Mapping
We conducted field work in a swath of c. 80 km length along the banks of the Greater Zab River (Figures 2  and 3). The terraces were mapped based on our field investigations and by using the TanDEM-X DEM (12 m horizontal resolution; Krieger et al., 2007) as well as satellite images (Esri World Imagery). Additionally, samples of sand layers within the river terraces were collected for luminescence dating, and structural data were collected during the field work along the transect for cross-section construction.

Estimating Fault Slip Rates From Fault-Related Uplift
We used the uplifted and tilted terraces of the Greater Zab River as markers for surface deformation. Then we constructed kinematic models that explain the observed surface deformation by fault-related folding at depth. Kinematic models for fault-related folding are well established (Brandes & Tanner, 2014 and references therein). The slip on faults underlying a folded succession can be constrained from the geometry of surface folding and vice versa (Brandes & Tanner, 2014). Mass balance across the structures is essential for models that describe the kinematics of deformation. Deformation of the geomorphic markers (e.g., river terraces) can be interpreted in terms of these kinematic models, and it can be related to the displacement on the underlying fault (e.g., Lavé & Avouac, 2000;Le Béon et al., 2014). The different kinematics for the deformation of the structures each provide a distinct deformation of the passive markers. If incision of river terraces is associated with tectonic uplift, the horizontal slip can be estimated using fault-related fold kinematics (Lavé & Avouac, 2000). This allows us to interpret the vertical uplift and tilting of terraces of the Greater Zab River in terms of horizontal shortening and displacement on relevant faults.
Both fault-related and detachment folds are present within our cross-section along the Greater Zab River (Section 3). The fault-related folds are characterized by flat-ramp geometries of the underlying fault and were modeled using the fault-parallel flow algorithm ( Figure 4a). This algorithm considers only materials in the hanging wall of a fault by translating particles along fault-parallel trajectories (Ziesch et al., 2014). Provided that the flat segments of the underlying detachment are horizontal, uplift only occurs when material moves over a ramp; tilting of the markers occurs only when they pass through axial surfaces ( Figure 4a). When the geometry of the fault at depth is constrained, a known amount of hanging-wall uplift of any given marker can be converted to slip on the fault itself using simple trigonometry. The displacement on the related fault for a defined period can be calculated from: where θ is the dip angle of the fault, u is the vertical uplift and d is the displacement on the related fault ramp for a defined period. Using fault-parallel flow, the displacement on the fault ramp will be equal to the displacement on the fault flat.
In case of a detachment fold (Figure 4b), the related fault is nearly horizontal, and the fold amplifies by accumulation of ductile layers in its core, leading to limb rotation (Mitra, 2003;Poblet & McClay, 1996). Detachment folds grow by limb rotation, hinge migration, or combinations of both (Mitra, 2003;Poblet & McClay, 1996;Scharer et al., 2006). Uplift and/or tilt of the markers (river terraces) occurs across the fold along with limb rotation and hinge migration (Scharer et al., 2006). The vertical uplift of the markers can be interpreted in terms of horizontal slip on the fault with consideration of mass balance, in which the uplifted area (A u ) in a section view above a specified datum across the fold is equal to the displacement area Fault-parallel flow Detachment fold (A d ) between the datum and the related fault (detachment) in the hinterland. Assuming that the original geometry of the uplifted river terraces was parallel to the present-day river allows calculating the uplifted (deformed) area above the river level (datum) since deposition of the uplifted terraces (Scharer et al., 2006;Thompson Jobe et al., 2017;Viveen et al., 2020). For the detachment folds, the excess area between the dated terrace level and the present-day river level was calculated and converted to the amount of shortening (slip on fault) using the area conservation method. Additionally, we assumed that the material moves parallel to the section (plane strain) and that there is no volume loss during the deformation.
Most studies (Burbank et al., 1996;Lavé & Avouac, 2000;Le Béon et al., 2014;Oveisi et al., 2009;Scharer et al., 2006;Thompson Jobe et al., 2017;Vance et al., 2003;Viveen et al., 2020: Yue et al., 2011 use the tilt of river terraces to reconstruct the pre-deformed state of the terraces, to calculate the deformed area, and to derive slip rates. In our study area, the bedrock underneath the river terraces mainly consists of highly erodible Late Miocene-Pliocene clastics. Furthermore, the river meanders. Therefore, the terraces here are highly dissected, and a certain terrace level cannot easily be traced to obtain the tilt angle. Consequently, we used the uplift rate difference between terraces on the ramp and those on the flats for the calculation of the slip rate in case of the fault-parallel flow kinematic algorithm, and we projected the associated terrace level to the cross-sections in order to calculate the uplift area for the detachment folds. The Move Software (Petroleum Experts) was used for constraining structural cross-sections, for the projection of river terraces, and for the calculation of the uplifted (excess) area.

Luminescence Dating
Luminescence dating determines the last daylight exposure of mineral grains such as quartz and feldspar, that is, the time of sediment deposition (cf., Fattahi, 2009;Preusser et al., 2008;Rhodes, 2011). For the present study, we sampled sand layers within fluvial sediments accessible in natural exposures along the river course, by forcing opaque tubes into the deposits. The tubes were sealed and sent to the laboratory, where they were opened. The samples were further processed under subdued red-light. The fraction 150-200 μm was gained by wet sieving, followed by removing carbonates and organic matter using HCl (10%) and H 2 O 2 (30%), respectively. A K-feldspar and a quartz fraction were isolated using heavy liquids, the latter was etched by HF (40%) for one hour.
First tests revealed that most quartz (usually the preferred mineral for dating) from the study region emits only weak luminescence signals, with no dominant fast component ( Figure 5a). This phenomenon has been reported from several young mountain ranges (e.g., Preusser et al., 2006) and led us to use K-feldspar for dating for all samples but 19K5, where quartz was also used in addition. While K-feldspar exhibits strong signals (Figure 5b), it is associated with some methodological challenges. The Infrared Stimulated Luminescence (IRSL) signal measured at low temperature (50°C) is usually considered unstable and tends to underestimate the real age of deposition, a phenomenon known as anomalous fading. To overcome this problem, a procedure has been developed, where during a second stimulation at higher temperature after initial readout a more stable signal is sampled (Buylaert et al., 2009;Zhang & Li, 2020). This method is known as post-IR IRSL (pIR) dating and different temperatures during the second stimulation step have been suggested. It applies that the higher the temperature, the more stable the signal, but at the cost of signal bleachability (the rate at which the signal is removed during daylight exposure). While partial bleaching can be recognized by the spread of replicate measurements, the relatively slow bleaching rate of the pIR signal might be challenging for samples, where exposure to daylight at deposition has been limited (Colarossi et al., 2015). Hence, using the pIR signal may cause overestimation of the real depositional age.
In this study, Equivalent Dose (D e ) measurements were done using a Freiberg Instruments Lexsyg device (Richter et al., 2013). For quartz Optically Stimulated Luminescence (OSL), we used green LED stimulation (60 s) with a detection filter combination comprising a Hoya U-340 and an AHF-Brightline 365/50 EX-Interference filter. A modified version of the Single Aliquot Regenerative Dose (SAR) protocol (Murray & Wintle, 2000) was used with preheat to 230°C for 10 s prior to all OSL measurements and by including an IR depletion test (Duller, 2003). K-feldspar was measured using a Schott BG-39 together with an AHF-Brightline 414/46 Interference filter. We used a three-step SAR protocol with preheating to 180°C for 10 s prior to optical stimulation (Figure 5c). The signal was recorded during three subsequent 90 s exposures to IR LEDs at 50°C (IRSL), 150°C (pIR-150), and 225°C (pIR-225). The performance of the two protocols was    tested by dose recovery tests, where OSL, IRSL and pIR-150 gave values close to 1.00. However, the pIR-225 protocol significantly underestimates the given dose in this test and also delivered much lower D e values (<50%) than IRSL and pIR-150. The results for this signal are hence omitted from the further discussion. We speculate this is explained by a change in trapping probability during preheating or by the previous optical stimulation (e.g., Qin et al., 2018;Zhang, 2017).
Most of the samples show non-Gaussian D e distributions, partially positively skewed, with overdispersion values ranging from 0.15 to 0.80 (Figures 5d-5f). The spread of values implies the presence of aliquots in which the luminescence was not fully reset prior to deposition signal and calls for the application of statistical approaches to extract the mean D e reflecting the time elapsed since the last daylight exposure (cf., Galbraith & Robert, 2012). For feldspar, sample 19K17 (IRSL) displays a normal D e distribution (Figure 5d) and the Central Age Model (CAM) was applied. The overdispersion of this sample (0.15) was used for the required input value sigma_b when testing the Minimum Age Model (MAM) and the Finite Mixture Model (FMM, using 2-5 components). The most appropriate age model was selected based on statistical tools such as likelihood and Bayesian Information Criterion (Table 1). For quartz, the CAM value was used despite the overdispersion value of 0.26, as the distribution was not skewed. The higher overdispersion is likely explained by a stronger effect of microdosimetry from surrounding grains (Mayya et al., 2006) compared to feldspar, where a substantial part of the dose rate comes from the internal K-content (c. 30% for the samples under consideration).
Samples for dose rate measurements of ca. 200 g were dried and the concentration of dose rate relevant elements (K, Th, U) was determined using high-resolution gamma spectrometry (cf., Preusser & Kasper, 2001). Dose rate and age calculations were done using the ADELE-2017 software (Degering & Degering, 2020). Average sediment moisture of 5%-15% during burial and an internal K-content of 12.5 ± 5% (Huntley & Baril, 1997) have been used in the calculation. Alpha efficiency (a-value) was assumed at 0.07 ± 0.02 for feldspar, while for quartz the outer layer was removed during etching. The cosmic dose rate was calculated after Prescott and Hutton (1994) taking into account longitude, latitude, and sample depth. The dosimetric data are compiled in Table 1.

Calculating Uplift Rates and Fault Slip Rates
We interpret the incision of the Greater Zab River as being caused by surface uplift, which changes the base level. Then the river incises to regain a steady state (Burbank & Anderson, 2012;Kirby & Whipple, 2012). The uplift rate for each sampled terrace was calculated from its elevation above the present-day river level divided by its age. The ages of the river terraces were obtained by luminescence dating using the combination of IRSL and pIR. The elevation of the terraces was measured from the TanDEM-X DEM (Krieger et al., 2007), which has a horizontal resolution of 12 m and a relative vertical accuracy of 2 m (slope ≤ 20%) and 4 m (slope > 20%), respectively. Applying the kinematic models of fault-related folding (Section 4.2) for deriving fault slip from surface uplift onto our regional cross-section (Figure 3), we calculated the slip rates on the relevant faults. The errors in the obtained terrace ages and the relative vertical accuracy of ±2 m were accounted for during the calculation of uncertainties in the amount of the uplift, the uplift area, and, hence, the slip rates of faults. However, there may be other source of errors, such as in depicting subsurface fault geometries and fault dips (Figure 3), which can hardly be accounted for without additional data, for example, from subsurface exploration. Additionally, there may be uncertainties in the calculated amount of river incision related to unsteady erosion (Gallen et al., 2015). This would then also affect the uplift derived from the incision. Gallen et al. (2015) propose to solve this problem by dating the lowest river terrace and using it as a reference instead of the modern river bed. However, in our case the lowest terrace levels were not at the same elevation above the river and they were not dated along the river. Therefore, we could not use them as a datum as Gallen et al. (2015) recommended.

Geomorphology of Greater Zab River Terraces
The Greater Zab River shows abundant terraces that are elevated from a few meters up to c. 300 m above the present-day river level (Figure 6a). Laterally, these terraces extend for more than 5 km from the river course especially within the synclines. Within the anticlines, terraces are restricted to a narrower lateral extent. The terraces extend to maximum 1 km from the river course within the High Folded Zone. During our field campaign, we mapped several levels of terraces (Figures 6a and 6b). The strath terraces are carved into the bedrock and host a thin layer of fluvial sediments, which usually consists of poorly sorted and well-rounded gravel-to sand-sized grains with a finer-grained matrix (Figure 6c). Pockets and lenses of sand are present occasionally (Figures 6c and 6d). To the NE of Perat Anticline, the lowest terrace is overlain by a 6-10 m thick layer of unconsolidated, fine-grained loess (Figure 6e). The terraces are partly very well cemented, especially in higher levels. This made collecting samples for luminescence dating at these high levels impossible. The deposits have been derived from a wide range of geological units within the tributaries of the Greater Zab River. Therefore, the lithology varies extensively. The thickness of the preserved terraces ranges from less than 1 m to more than 20 m.
The vertical offset between the higher successive terraces that are far from the river is usually high and the vertical spacing decreases between the low terraces that are close to the river (Figures 6a and 6b). This can either be related to the lack of preservation of older terrace levels due to erosion, or to higher uplift rates during their formation (Demoulin et al., 2017). In some places the vertical spacing of the lower terraces is so low that it is difficult to tell different terrace generations apart. Bedrock is clearly exposed below the high terraces; therefore, we considered them as strath terraces. In contrast, the lower terraces are aggradational or degradational fill terraces. Along the river, the maximum elevation of the terraces above the river decreases SW-ward and does not exceed 50 m in the southwest (Figures 6a and 7). In the Safin and Sarta anticlines, the terraces are dissected by lateral streams and a single terrace cannot be traced across the anticline. To the southwest on the Girda Baski Anticline and another small anticline that crosses the river and plunges E-ward before reaching the trace of our transect (Figure 2a), single terraces are preserved across the anticlines (Figure 6a) and terrace profiles follow the tilting caused by ongoing folding there (Figure 7).

Age of the River Terraces
For all 15 samples investigated both IRSL and pIR-150 ages are available. These ages range from 6.2 ± 0.4 ka (19K5 IRSL) to 165 ± 12.7 ka (19K21 pIR). For each individual sample, the pIR age is higher than the IRSL age (Figure 8) but the age offset ranges from 12% (19K29) to more than 100% (19K7) and is hence not systematic. Two phenomena can likely explain the observed age offset, either the effect of signal instability (fading) or that of partial bleaching. In the first case, the IRSL ages would underestimate, in the second case the pIR ages would overestimate the real age of deposition. In fact, both phenomena could occur simultaneously, in which case the true age might be somewhere between the IRSL and the pIR ages (if the latter does not underestimate as well). However, the non-systematic difference between IRSL and pIR is interpreted to represent incomplete bleaching of the pIR signal as fading usually produces rather systematic offsets. This  19K27 19K25 19K31 19K23 19K29 19K21 19K15 19K17 19K19 19K13 19K1 19K3 19K5 19K7 19K9 U  19K25 19K31 19K23 19K29 19K21 19K15 19K17 19K19 19K13 19K1 19K3 19K5 19K7 19K9 300 200 100 0 Elevation (m)

SW NE
view is supported by the OSL age of sample 19K5 (6.0 ± 0.5 ka), which is in excellent agreement with the IRSL age of 6.2 ± 0.4 ka; the pIR age of 7.4 ± 0.4 ka is some 20% higher. Hence, we will rely on the IRSL ages in the discussion of uplift rates.

Relative Uplift and Associated Slip on Faults
During the Pliocene-Quaternary, crustal thickening within the Zagros Belt was associated with continuing Arabia-Eurasia collision that led to regional uplift (Mouthereau, 2011). Uplift of terraces within the synclines is less likely to be associated with uplift due to any underlying local faults (Figure 7). It can either be related to the regional uplift from the crustal thickening or to the uplift related to the slip on deeper basement faults. We detect the lowest uplift rates in the Harir Syncline (Figure 9a). Sample 19K15 indicates an uplift rate of 0.1 ± 0.04 mm a −1 . The uplift of this sample can neither be attributed to displacement on any fault-related detachment folds to the SW above the Triassic detachment, nor to the displacement on the Mountain Front Fault as it is located on a flat (Figures 3 and 7). It can, however, be correlated to the regional uplift due to crustal thickening and/or basement thrusting. A constant uplift rate of 0.1 ± 0.04 mm a −1 throughout the Pliocene would explain the present-day topography (400-500 m above sea-level) of the northeastern synclines (e.g., Harir Syncline) of the Foothill Zone close to the Mountain Front Flexure, which do not show evidence of distinctive erosion. Therefore, we interpret a regional background uplift rate of 0.1 ± 0.04 mm a −1 as caused by crustal thickening or by displacement on the basement thrust to the SW of the Mountain Front Fault. However, the displacement on a basement thrust to the SW of the Mountain Front Fault below the Safin and Sarta anticlines is more plausible since this basement thrust has been constrained with a higher displacement along-strike to the SE of our study area (Koshnaw et al., 2020;Le Garzic et al., 2019). Also, uplifted river terraces within synclines to the SW of the Sarta Anticline are rare (Figure 6a). This background uplift rate was subtracted from the uplift rates of all other samples for calculating the horizontal slip on faults responsible for uplifting the river terraces in the affected detachment folds.

Basement Thrusts
The samples that are located within the Hanara Syncline to the NE of the Perat Anticline, that is, in the hanging wall (ramp) of the Mountain Front Fault (sample set Ta; Figures 9b and 9e), and the samples in the Harir Syncline are located on a flat (sample set Ta; Figures 9a and 9e). Due to their age difference, the corresponding terrace levels cannot be categorized and fitted between the Ta and Tb sets. Therefore, we use the average uplift rates here. The terrace set Ta shows an average uplift rate of 0.6 ± 0.17 mm a −1 and a nearly linear age-elevation relationship. The set Tb does not exhibit linear age-elevation trends (Figure 9c). For instance, samples 19K15 and 19K17 have the same ages within errors, but sample 19K17 is located c. 26 m above sample 19K15 (Figures 9c and 9d). This indicates fast aggradation of the terraces within the Harir Syncline around 62 ka and then subsequent incision (Figure 9c). The aggradation of terraces there can be caused by the growth of the Safin and Sarta anticlines downstream, by a climatic force, or by a landslide, which led to a local decrease of channel steepness upstream. Terrace sample 19K15 is now c. 6 m above the river level. This means that the actual uplift rate there is 0.1 ± 0.04 mm a −1 since 62.7 ± 4.1 ka. The Perat Anticline itself forms the Mountain Front Flexure and the frontal anticline of the High Folded Zone. Accordingly, the samples located to the NE of the anticline are in the hanging wall above a NE-dipping ramp segment of the Mountain Front Fault (sample set Ta, Figure 9e), whereas sample set Tb is located above the flat segment of this fault. This means the difference between the Ta and Tb uplift rates is 0.5 ± 0.20 mm a −1 , and this can be attributed to the displacement on the ramp that uplifted terrace set Ta. Using an average dip angle of the underlying fault of c. 20° towards the northeast (Figure 9e) and by applying the fault parallel-flow algorithm, we calculate 1.46 ± 0.60 mm a −1 of slip on the Mountain Front Fault in the last c. 62 ka on average ( Table 2). The higher terraces on the southern limb of the Perat Anticline are more than Figure 8. Comparison of ages determined using different approaches reveals that pIR ages are systematically higher than those determined by infrared stimulated luminescence. This is interpreted to result from incomplete resetting of the pIR signal prior to deposition.  (Figure 6b). This indicates that there is no ongoing limb rotation in the Perat Anticline. If there is any horizontal slip in the anticline at present, this slip must have been caused by motion along the underlying thrusts.
For the basement fault to the SW of the Mountain Front Flexure that is underlying the Sarta and Safin anticlines (Figures 3 and 7), we do not have good constraints on the slip rate. However, if the 0.1 ± 0.04 mm a −1 uplift of the terrace sample 19K15 that is located within the Harir Syncline is associated with displacement on this thrust, then the slip on this thrust must not exceed 0.41 ± 0.16 mm a −1 in order to cause 0.1 ± 0.04 mm a −1 of uplift by applying the fault parallel-flow algorithm and a fault dip of 14º (Table 2). Therefore, we have accounted the 0.41 ± 0.16 mm a −1 of slip on this thrust when calculating the horizontal slip related to the detachment folding for the Safin and Sarta anticlines.

Detachment Folds
The uplift of terrace samples on the detachment folds to the SW of the Mountain Front Flexure was interpreted to result from the ongoing shortening and amplification of these folds, and it was translated in terms of horizontal slip on the detachment fault. The uplift of the dated terrace level caused by the detachment folding has been converted to the uplift area that is considered to be equal to the displaced area on the fault. The anticlines are deformed above the Lower Triassic detachment and this detachment depth was used to covert the deformed area into the displacement (Figure 4b). Here, the assigned terrace levels do not correspond in different anticlines, since there is no age constraint for all of the terrace levels.
Sample 19K25 is located on the crest of Safin Anticline and it is 60.2 ± 3.4 ka old. Its elevated c. 42 m above the present-day river and indicates total uplift rates of 0.7 ± 0.07 mm a −1 . It was taken from the terrace level T4 (Figures 10a and 10b). The uplift due to slip on the basement thrust below the Safin Anticline was deducted and the uplift area below the terrace level T4 was calculated. The detachment level is at a depth of 6,350 m, and by applying the area conservation method, it requires a minimal horizontal shortening (displacement on the associated fault) of c. 24 m. This results in a horizontal slip rate of 0.40 ± 0.10 mm a −1 over the last 60.2 ± 3.4 ka for the Safin Anticline (Table 2), where it is crossed by the Greater Zab River. None of the layers and terraces shows tilting in the crest of the anticline. In contrast, the higher terraces on both limbs of the anticline were tilted by 5-6°. The tilt fits the limb rotation during fold growth ( Figure 11).
Sample 19K27 was collected in the crest of the Sarta Anticline. The sample is 31.2 ± 3.1 ka old and it is currently located c. 42 m above the river within terrace level T4 (Figures 10c and 10d), which means it has been uplifted at a rate of 1.35 ± 0.2 mm a −1 . Using the same approach as for the Safin Anticline and a detachment depth of 6,550 m, this uplift requires c. 38 m of horizontal displacement. The corresponding horizontal slip rate on the detachment is calculated to be 1.24 ± 0.36 mm a −1 over the last 31.2 ± 3.1 ka (Table 2).

Basement Versus Cover Deformation
Within the Foothill Zone of the Kurdish part of the Zagros Mountains, elevated river terraces are present along the bank of the Greater Zab River (Figure 6a) and along other rivers (Tigris, Lesser Zab, and Sirwan). These terraces are located to the NE of the Zagros Deformation Front. Topographically still higher terraces levels are located close to the Mountain Front Flexure. The elevation of the terraces gradually decreases  (Figures 6a and 7). Terraces with elevations of c. 50 m and more above the present-day river course are restricted to the NE parts of the Foothill Zone. Here, they are present only within an elongated zone of 40-50 km width to the SW of the Mountain Front Flexure (Figure 12a). The SW boundary of this zone is located approximately at the 400 m smoothed contour line (dashed pink and white line in Figure 12a). To the SW of this line, the terraces are low and only present where the river cuts through anticlines. These terraces have been uplifted only from the fault related folding above the Triassic detachment. In addition, there is a change in the regional topographic gradient across this boundary line (see the swath profiles in Figures 12b-12d). For example, in the swath profile 2 (SP 2) in Figure 12c, the elevation increases from 100 to 400 m within a distance of 120 km while it increases from 400 to 700 m within a distance of 35 km. This change in topography coincides with a NE-ward increase of the crustal thickness (Jiménez-Munt et al., 2012). Therefore, we interpret that thickening of the crust and thrusting in the basement only occurs NE of the boundary line. Also, there are a few thrust-related earthquakes to the NE of this line. These earthquakes occurred at depths of 8-15 km, which is clearly within the Paleozoic or Precambrian basement (see cross section of Figure 3). This is similar to earthquake depths reported from the Iranian part of the Zagros (Nissen et al., 2011;Talebian & Jackson, 2004).
The folding and thrusting of the sedimentary cover above the basal detachment within the Foothill Zone has started as early as Late Miocene (Koshnaw et al., 2017). The onset of basement thrusting related to the formation of the Mountain Front Flexure was also dated to the Late Miocene (Homke et al., 2004;Koshnaw et al., 2017). The deformation above the basal detachment within the Foothill Zone can be related to the basement thrusting below the Mountain Front Flexure and the more external basement fault towards SW, as portrayed in several regional cross-sections (Le Garzic et al., 2019;Tozer et al., 2019;Zebari et al., 2020). In other words, the displacement on the basement thrusts in the NE of the Foothill Zone is transferred up-section into folding and thrusting above the Triassic detachment in the Foothills Zone. The fact that  deformation in the sedimentary cover above the detachment occurs as far as 140 km to the SW from its kinematically related basement thrust suggests very low friction along the basal detachment.

Impact of the Oblique Convergence
Due to the present-day Arabia-Eurasia oblique convergence, the shortening in the NW Zagros is mostly partitioned between right-lateral slip along the Main Recent Fault in Iran (Reilinger et al., 2006;Vernant et al., 2004) and belt-normal shortening across the folded belt (Hessami et al., 2006;Khorrami et al., 2019;Vernant et al., 2004). However, Late Pliocene-Recent paleostress analyses (Navabpour et al., 2008) and present-day geodetic surveys (Khorrami et al., 2019;Reilinger et al., 2006) denote N-S directed shortening in the NW part of the folded belt. The oblique collision also impacts the right-lateral slip in some N-S trending faults such as the Khanaqin and Kazerun faults that cross the belt (Authemayou et al., 2006). It also has imprints on the structures there, for example, wrenches in the fold axes and formation of en-échelon fold structure (Csontos et al., 2012). The structures investigated in the present study have NW-SE trends oblique to the shortening direction, while the calculated slip rates from the uplift of river terraces here account for the structure-orthogonal component of slip along the cross-sections. Therefore, an oblique slip is expected on these structures. This can be evidenced by focal mechanisms of earthquakes, indicating active oblique thrusting (Figures 1 and 12a; Nissen et al., 2019;Yang et al., 2018), as well as by fault slip data ( Figure 13).
The actual slip rates can be calculated using the following equation:   cos calculated slip rate Actual slip rate where α is the angle between the actual slip direction and the direction of our transects (nearly perpendicular to the trend of structures) that we used for calculation of slip rates. Our own field observations of fault-slip data collected from a reverse fault on the NE limb of Safin Anticline, close to the crest with dip direction/dip angle of 056/60 indicate an oblique reverse slip vector with an average rake angle of 65° (Figure 13). This implies an obliquity of 25° for the actual slip rates, which increases our calculated slip rates by c. 10%. However, it is not certain if this 25° of obliquity represents the present-day slip direction, because the absolute age of the fault slip is not known, and the fault-slip data is only limited to a single structure. Obtaining the most recent slip direction on the oblique faults requires more comprehensive analyses of fault slip data, which is beyond the scope of this paper. In case of slip on horizontal detachment faults, the obliquity of the corresponding fold axis and the direction of the cross-sections ( Figure 10) with respect to the shortening direction should also be taken into account. The anticlines to the SW of the Mountain Front Flexure have an obliquity of 30°-50° to the shortening direction. This suggests that the actual slip rates  (Delvaux & Sperner, 2003) illustrating the angle between the actual slip direction and that of the constructed regional cross-section in Figure 3. Red great circle and slip arrow characterize the average value of the fault-slip data. Shaded gray area defines the dilational quadrants of a right-dihedra solution for the fault-slip data, indicating an oblique reverse mechanism (compare it with the reverse earthquake focal mechanisms of the study area in Figures 1 and 12). Ac tua l sli p C a lc u la te d s li p could be 15%-55% higher than those calculated along the cross-sections, assuming that the slip on these faults occurs parallel to the recent N-S shortening direction. Table 2 shows the calculated fault slip rates for all studied structures. The total slip rate on the two basement thrusts that emerge to the Triassic detachment level is 1.87 ± 0.76 mm a −1 . The total horizontal slip rate accommodated by detachment folding on the Triassic detachment level across the Safin and Sarta anticlines is 1.64 ± 0.46 mm a −1 (Figure 14a). In our structural model, the displacement on the basement thrusts to the NE including the Mountain Front Fault and the external basement fault is accommodated up-section by those detachment folds. Additionally, there are several anticlines to the SW of the study area close to the limit of the Zagros Deformation Front (Figure 14b) that have not been investigated in this study. Therefore, the total horizontal slip rate that triggered folding of the sedimentary cover above the detachment within the Foothill Zone and along a transect in a direction parallel to the transect shown in Figure 3 is expected to be higher than calculated here. It has been documented from morphometric analyses and tilted river terraces that these southwestern anticlines are also accommodating ongoing shortening (Fouad, 2012;Obaid & Allen, 2017).

Spatial Distribution of Fault Slip Rates
Within the two studied detachment folds, the dated terrace level in the Sarta Anticline shows higher uplift rates than that of the Safin anticline, and the horizontal slip rate in the Sarta Anticline is also much higher than the slip rate in the Safin Anticline. Sarta Anticline is located NW of the Bina Bawi and Pirmam anticlines ( Figure 2a) that mark the frontal anticlines of this segment of the Mountain Front Flexure to the SE of the Greater Zab River. The basement thrust underlying this segment extends NW-ward and emerges into the Triassic detachment below the Sarta Anticline.
The obtained fault slip rates here are average values for the Late Pleistocene-Holocene. They are lower than the GPS velocities of 3.7 ± 0.4 mm a −1 of belt perpendicular shortening for this part of the Zagros Belt ( Figure 14b; Reilinger et al., 2006). The cumulated shortening rates from the two detachment folds or those from the two basement faults in this study sum up to less than half of the geodetically obtained shortening rates in the belt. This difference can also be related to the fact that this study does not cover the entire belt and that the actual slip rates due to the oblique convergence have not been determined. Due to the unavailability of comprehensive geodetic studies with high resolution, it is not possible to make a comparison of slip rates within specific structures. Similar studies on the frontal anticlines in the Fars Arc show that these anticlines accommodate much higher slip rates (up to 3-4 mm a −1 ), which account for up to 50% of the geodetically derived shortening rates (Oveisi et al., 2007(Oveisi et al., , 2009

Conclusions
We interpret extensive flights of alluvial terraces around the Greater Zab River in the NW Zagros Fold-Thrust Belt to have formed due to ongoing regional uplift related to fault-related folding during the Late Pleistocene to Holocene. Our results suggest a slip rate of 1.46 ± 0.60 mm a −1 on the Zagros Mountain Front Fault and less than 0.41 ± 0.16 mm a −1 on a more external basement fault to the SW based on the relative uplift rates of terraces in their hanging walls. The horizontal slip rate related to detachment folding of two anticlines within the Foothill Zone is 0.40 ± 0.10 mm a −1 and 1.24 ± 0.36 mm a −1 , respectively. These slip rates were calculated along transects oblique to the present-day N-S convergence direction; therefore, the actual slip rates must be higher than these calculated rates. Within the Foothill Zone, relatively higher river terraces and higher topography are restricted to the NE part of this zone. This, together with evidence from topography, seismic records and the crustal thickness suggests the occurrence of crustal thickening and/or basement thrusting there. While deformation in the SW part of the Foothill zone is limited to folding and thrusting of the sedimentary cover above the basal detachment, deformation in its NE part is associated with reverse faulting on a basement thrust. Our study shows that geomorphic signals of distributed neotectonic deformation recorded by alluvial terraces reveal processes acting at depth when combined with a careful examination of the structural setting, alongside with kinematic fault modeling techniques. Thus, we were able to infer ongoing uplift at rates of only a few millimeters per year. Applying our method to other river terrace systems in the area or elsewhere in the Zagros Mountains may hence be the key to obtain more robust long-term slip rate estimates, which can be compared with GPS-derived, shorter term slip rate estimates. Discrepancies in estimates obtained via these two approaches might thus allow identifying locked fault strands that pose an increased seismic hazard.

Data Availability Statement
The data derived from luminescence dating of river terraces and used in this study can be accessed at PAN-GAEA: https://doi.pangaea.de/10.1594/PANGAEA.933900 (Zebari et al., 2021). The German Academic Exchange Service (DAAD) is acknowledged for providing a scholarship (Research Grants-Doctoral Programs in Germany, 2016/17; 57214224) to M. Zebari to conduct his PhD research in Germany. P. Navabpour and K. Ustaszewski express their gratitude to the German Research Foundation (DFG) project no. 393274947 for providing additional financial support. The German Aerospace Agency (DLR) is thanked for providing TanDEM-X digital elevation models "© DLR 2018." Gamma spectrometric measurements were carried out and analyzed by A. Fülling, Freiburg. The open access costs were covered by the Projekt DEAL (www. projekt-deal.de). Burhan Hostani, Ferix Reshoo, Kaiwan Fatah, Renjbar Harni, and Soran Syan helped us during field work. Thoughtful reviews from Jessica Thompson Jobe and two anonymous reviewers have helped to further improve our work.