Constraining a model of punctuated river incision for Quaternary strath terrace formation

In the small fraction of Earth's surface with the highest erosion rates such as the Alps and Himalayas quantifying rates of incision, rock uplift and inferring climatic controls on the landscape can be relatively straightforward once the ages of river terraces cut in bedrock (strath terraces) are constrained. However, in many mid to lower relief settings that are more typical of mountain belts worldwide, periods of net river incision and riverbed lowering are relatively short (punctuated), interrupted by long periods of sediment aggradation or stasis. We define a conceptual model of punctuated river incision and strath terrace formation for the calculation of incision and rock uplift rates, and recommend strategies for geochronological sampling and interpretation. An approach using OSL dating of terrace gravels allows us to constrain a detailed ~150 kyr history of punctuated river incision and strath terrace formation spanning two stratigraphic landform levels in the High Atlas Mountains (NW Af- rica). Extensive preservation and exposure of strath-top gravels, within a post-orogenic setting unaffected by eustatic influences, enables the derivation of rates of base-level fall, integrated over periods of strath-top aggradation and incision, that are consistent with independently constrained regional rock uplift rates. Combining a punctuated river incision model with our well-constrained terrace formation history allows us to demonstrate how assumptions concerning Quaternary river incision and aggradation can lead to the problematic Sadler Effect, an apparent dependence of incision rates on measured time interval. Subsequently, we demonstrate that an approach to reinterpreting previously published data using the punctuated incision model, even when combined with limited terrace age data, results in more consistent conclusions about rates of river incision, rock uplift and base-level lowering across the mountain belt. Our recommendations for sampling strategies to constrain rock uplift rates require samples to be taken just above the strath surface, and in addition towards the top of the deposit for river incision rates. In a setting with punctuated river incision and strath terrace formation, both rock uplift and incision rates require burial dates, as exclusive use of abandonment ages will not yield constraints on accurate rates of rock uplift or incision. Furthermore, we find that only with multiple along-stream locations and multiple burial dates in each terrace deposit, could a reliable climatic signal be extracted: this signal would not have shown up in terrace abandonment ages such as those derived from cosmogenic exposure dates. The demonstrated effects of assumptions about strath terrace formation, and the recommended approaches for sampling and interpretation, have implications for those attempting to constrain palaeoclimatic, tectonic, and geomorphic histories from strath terrace records in regions exhibiting punctuated river incision. This model, even when applied to limited terrace age data, has solved for the Sadler Effect and reconciled rates between studies across the mountain belt. Our punctuated incision model allows for strath formation that is not necessarily cyclical or directly correlated with climatic stages. The 150 kyr Quaternary terrace age-model constrained here shows that in the High Atlas, strath terrace abandonment and river incision are asynchronous between sites and not convincingly correlated with climatic stages, but aggradation experiences a pulse during the interglacial maximum at MIS 5e. Whilst a history of punctuated river incision complicates the interpretation of strath terrace records for constraints on rates of incision and rock uplift, combining the punctuated river incision model and our recommendations on sampling and interpretation of strath terrace age data allows for more consistent and accurate insights into the climatic control on erosion and landscape evolution, Quaternary river incision and rock uplift rates. of


Introduction
The long-term evolution of mountain landscapes is controlled by the competition between rock uplift and bedrock river incision. Strath terraces are benches cut into underlying bedrock by river erosion, which can be topped with a body of river-transported sediment. Within mountains, strath terraces form when a river incises, abandoning the former floodplain and channel belt (Leopold et al., 1964), and thus record a river's downcutting to its modern position in a river valley bottom. Successive periods of fluvial sediment aggradation, river incision and floodplain-channel abandonment can form inset river terrace staircases commonly preserved along river valley margins around the world (Starkel, 2003). In these settings, the development of numerical dating of Quaternary river sediments and rock surfaces (Rixhon et al., 2017) has enabled strath terrace staircases to be commonly used to derive histories of landscapes evolving under changes in climate (Bridgland and Westaway, 2014;Schanz et al., 2018), and to quantify river incision and rock uplift rates (e.g. Lavé and Avouac, 2000;Fuller et al., 2009;Burbank and Anderson, 2011;Pazzaglia, 2013;Zondervan et al., 2020a). In the small fraction of Earth's surface with the highest erosion rates such as the Alps and Himalayas (Burbank et al., 1996;Pan et al., 2003), such quantifications are straightforward once the age of strath terraces is constrained. Such high erosion rate settings are characterised by the highest relief and histories dominated by river incision with brief periods of strath terrace formation, dictated by a cyclical rhythm of Quaternary climate (e.g. Pan et al., 2003;Amos et al., 2007). However, in many mid to lower relief settings river incision is interrupted by long periods of sediment aggradation or stasis (Wegmann and Pazzaglia, 2002;Collins et al., 2016;Schanz and Montgomery, 2016;Foster et al., 2017). The potentially complex responses of such landscapes to climatic changes (Schumm, 1979;Bull, 1991;Wegmann and Pazzaglia, 2002) and the limits of geochronology and preservation when constraining such river incision histories make strath terrace records more challenging to interpret.
Early theories on river incision and strath terrace formation suggested that some settings experience episodic rather than progressive steady incision, where periods of net river incision and riverbed lowering are relatively short (punctuated), interrupted by long periods of sediment aggradation or stasis of the riverbed ( Fig. 1; Bull, 1991;Schumm, 1979). Such punctuated river incision histories have been found using radiocarbon dated strath terraces over a short Holocene timescale (Wegmann and Pazzaglia, 2002;Collins et al., 2016;Schanz and Montgomery, 2016). These studies demonstrate that interpretation of incision rates from such terraces risk the derivation of incorrect or ambiguous incision rate histories. The risk of deriving inaccurate rates from strath terraces follows a wider trend observed in field records of river incision over timescales of several magnitudes approaching 10 7 yr, where erosion rates show an apparent increase towards the present as a bias of time-intervals between measurements ( Fig. 1; Gardner et al., 1987;Mills, 2000;Finnegan et al., 2014). This phenomenon is generally referred to as the "Sadler effect" (Sadler, 1981) named after the same effect found in the sedimentary record (Sadler, 1981), and is a major challenge in the interpretation of records of river incision to derive rates of rock uplift or climatic histories. Nevertheless, recommendations for constraining accurate rates from strath terraces in settings with punctuated river incision histories have been mixed.
For example, Wegmann and Pazzaglia (2002) recommend integrating rate calculations over Quaternary terraces representing at least one glacial-interglacial cycle of incision and aggradation, a similar approach to using the lowest strath terrace instead of the modern riverbed as the reference frame (Gallen et al., 2015). In contrast, Schanz and Montgomery (2016) point out that using the method proposed by Wegmann and Pazzaglia (2002) and Gallen et al. (2015) without taking into account the full range of terrace aggradation would still result in erroneous uplift rates. However, Schanz and Montgomery (2016) did not find an approach to derive accurate rates. None of these collective studies were able to constrain a full history of terrace aggradation and incision over a glacial-interglacial timescale, limited by preservation and the short time range of radiocarbon dating. In addition to the challenge of deriving rates of incision and rock uplift from strath terraces formed via punctuated river incision, studies have questioned the role of climate in Quaternary strath formation: cosmogenic nuclides show that strath terraces can be left un-incised over hundreds of thousands of years (Foster et al., 2017), and a global review demonstrates that strath terraces do not necessarily form during any particular climatic stage (Schanz et al., 2018). These studies suggest that strath terrace formation is not always directly related to Quaternary climate cycles. However, the limited geochronological control on strath-top aggradation histories means that potential signals of climate-driven aggradation within such strath-top fills may remain unexplored. Furthermore, whilst a lack of climatic control on the formation of Quaternary strath terrace staircases in settings with punctuated river incision can have important implications for river incision and rock uplift rate calculations, to fully understand such implications requires a conceptual model relating these, supported by data of Quaternary strath-top aggradation, incision and rock uplift. Overall, the aforementioned studies demonstrate that using strath terraces to constrain the climatic control on landscape evolution and to quantify rates of incision A) The gradual incision model follows the assumption of a normative steady-state incision, which can be briefly interrupted by cyclic climatic changes. B) In the punctuated incision model (modified from Bull, 1991), prolonged periods of meta-stable equilibrium dominate with short intervals of net aggradation and incision. C) metastable equilibrium can be understood as a non-optimised, but locally stable position, where threshold exceedance can tip the system into the next stable state. or rock uplift is riddled with pitfalls in settings with punctuated river incision. These studies also demonstrate the challenge of constraining a Quaternary aggradational record of strath terraces in such settings, as cosmogenic nuclides mainly date the abandonment of terrace surfaces, and radiocarbon dating cannot extend back a full 100 kyr glacialinterglacial maximum cycle (Rixhon et al., 2017). Whilst burial dating of Quaternary strath-top sediment using optically stimulated luminescence (OSL) can cover hundreds of thousands of years (Rixhon et al., 2017), this technique presents challenges in dating coarse grained gravel deposits representing bedload aggradation in strath terraces (e.g. Stokes et al., 2017).
Consequently, whilst absolute dating techniques now allow the age of Quaternary strath terrace deposits and surfaces to be determined; how rock uplift, river incision and climatic controls on landscapes are interpreted from such constraints in settings with punctuated river incision needs to be better understood. In addition, the challenges of dating strath terrace gravels over and across multiple full 100 kyr glacial interglacial cycles, and the derivation of full aggradational histories remain key problems. Addressing these challenges requires a dataset of aggradation on strath terraces in a setting with extensive preservation and exposure of strath-top gravels over a full eccentricity cycle (~100 kyr), unaffected by eustatic sea-level changes, and an independently constrained rock uplift rate. Such an opportunity exists in the intracontinental south-central High Atlas Mountains, where post-orogenic regional uplift rates since 5 Ma are isolated from eustatic sea-level change, and strath terrace deposits are extensively preserved and exposed. Our approach to dating strath terrace gravels using OSL enables us to derive a full aggradational history on these straths over ~150 kyr. This dataset allows us to constrain a model of punctuated river incision and strath terrace formation over Quaternary timescales. In this study, we advance the interpretation of strath terraces in settings with punctuated river incision by: (i) defining a conceptual model of punctuated river incision and strath formation for the interpretation of numerical dates; (ii) constraining an age model of punctuated river incision and its relationship to climate and rates of incision and rock uplift using our OSL approach in the High Atlas, and (iii) recommending strategies for strath terrace sampling and interpretation of geochronology, even in settings with limited age constraints or sediment preservation. Our constraints on a model of punctuated river incision for Quaternary strath terrace formation have implications for tectonic geologists, geomorphologists and sedimentologists attempting to derive palaeoclimatic, tectonic, and geomorphic histories from strath terrace records in regions exhibiting punctuated river incision.

Conceptual framework
In mountain belts with high rates of erosion, river incision is often assumed to either be in steady-state equilibrium, or to be readjusting to one while temporarily perturbed, over timescales of 10 4 -10 7 yr (Whipple, 2001;Burbank and Anderson, 2011). River incision tending towards a normative steady-state equilibrium, with periodic climatic perturbations forming terraces, is a simple, established model that has been used in many studies (e.g., Lavé and Avouac, 2000;Hancock and Anderson, 2002;Fuller et al., 2009;Pazzaglia, 2013;Zondervan et al., 2020a). Here, we call this the 'gradual incision' model ( Fig. 1A). An extreme version of this model is one where vertical incision is constant, progressive and uninterrupted, where terraces form by the sinuous movement of the channel across and into the valley floor (Limaye and Lamb, 2016). If a set of conditions hold, correlation between terrace abandonment ages will yield rates that approximate to both actual river incision rates and the long-term integrated base-level lowering rates, which often reflect rock uplift rates (Fig. 1A). These conditions include: (i) incision approximates steady-state gradual downcutting, where periods of non-incision during terrace aggradation are brief; or where terrace formation follows a consistent periodicity often taken to be coincident with Milankovitch-forced climatic cycles (e.g. Gallen et al., 2015), and (ii) where the thickness of deposits on top of terraces are approximately consistent between levels (Fig. 1A). If these conditions hold, cosmogenic exposure dates of strath-top surfaces, a technique that ranges across thousands to millions of years (Rixhon et al., 2017), enables extensive constraints to be placed on such landscapes.
When none of the conditions listed above, that make the interpretation of strath terrace geochronology straightforward, are met, such a setting approximates to a model of punctuated river incision and strath formation that we present here (Fig. 1B). As such, this model is an endmember on the opposite side of the spectrum from the gradual incision model (Fig. 1A). The model of punctuated river incision and strath formation defined here is a combination and modification of complex response theory (Bull, 1991) and the dynamic equilibrium model (Schumm, 1979). Instead of steady-state incision interrupted by brief periods of strath formation, strath occupation and aggradation during periods of metastable equilibrium dominate the riverbed history ( Fig. 1, Wegmann and Pazzaglia, 2002;Collins et al., 2016;Schanz and Montgomery, 2016). Rather than a climatically-forced periodic cyclicity of strath formation, local thresholds needed to exceed metastable equilibrium vary in space and time (Schumm, 1979;Phillips, 2006;Fryirs and Brierley, 2010) and results in an incision history that is not faithfully cyclical (Fig. 2, e.g., Foster et al., 2017). Instead of similar extents of strath surfaces and deposits, these vary widely between stratigraphic levels ( Fig. 1). Thus, in settings described by the metastable equilibrium model of punctuated river incision, strath surface abandonment ages will not capture true river incision rates, and river incision rates are decoupled from the long-term base-level lowering rate (Fig. 1B).
Using the OSL chronology derived in this study we will demonstrate how the metastable punctuated river incision model for Quaternary strath terrace formation allows pitfalls in the interpretation of geochronology to be prevented. We will show how this model can help derive consistent rates and interpretations of strath terraces across a mountain belt, even when age data is limited. Consequently, future workers will be able to use this model and the approaches detailed here to plan their geochronological sampling and interpretation in settings where river incision histories might be punctuated. We will further discuss the use of this model in the discussion, using the constraints established below.

Study area
The study area is located on the southern flank of the central High Atlas Mountains, where rivers drain into the intracontinental Ouarzazate basin (Fig. 2), unaffected by cyclic eustatic sea-level fluctuations, with the Atlantic coast >800 km downstream from the mountain front (Boulton et al., 2014). Though localised base level fluctuations are likely to have occurred, studies of basin fill sediments and river profiles suggest that such fluctuations were small to negligible over the Plio-Quaternary (Boulton et al., 2014(Boulton et al., , 2019. Long term relative base-level fall in the High Atlas is driven by a low rate of rock uplift: long-term isostatic rock uplift rates across the central High Atlas have been 0.17-0.22 mm yr − 1 due to lithospheric thinning, calculated from 1200 m of uplifted Messinian deposits across the Atlas area from the Middle Atlas in the north to the Anti-Atlas in the south, dated at 7.1-5.3 Ma . Over the Quaternary local tectonic advection is very limited in the High Atlas, especially upstream of the mountain front as exemplified by terrace-treads of the Dades River forming equably spaced and parallel paleo-river long profiles throughout the fold-thrust belt and thrust front (Stokes et al., 2017). In the western Ouarzazate basin, landforms such as river terraces in the incising foreland basin alluvial fan systems record localised fault movement of up to 0.3 mm yr − 1 (Pastor et al., 2012).
Drainage development in the High Atlas is considered to relate primarily to the exhumation of structurally distributed lithologies with different hardness (Zondervan et al., 2020b), controlling where river terraces develop (Stokes et al., 2017). Fluvial incision dominates bedrock erosion, with the influence of Pleistocene glacial processes restricted to elevations above ~3300 m (Hughes et al., 2004). The Quaternary climatic history of the area is recorded in northwest African climate humidity and Green Sahara Periods which follow a signal strongly controlled by precessional cycles (~26 kyr) albeit still paced by eccentricity (~100 kyr) (Tjallingii et al., 2008). Quaternary strath terrace abandonment dating has been previously applied to five levels of strath terraces over the last 200 kyr along the Madri River in the Ouarzazate basin using cosmogenic nuclides ( Fig. 2; Arboleya et al., 2008), and on one terrace level upstream of the mountain front along the Dades River dating back to up to ~100 kyr using OSL of fluvial overbank sands capping terrace fluvial gravels (Fig. 2;Stokes et al., 2017). These studies report rates of incision which are inconsistent up to an order of magnitude from 0.1 to 5.2 mm yr − 1 , reporting an increase of rates towards the present, and correlation with climatic cycles is not straightforward . These studies also demonstrate the challenge of dating terraces in the High Atlas, since cosmogenic nuclide dating of strath surface abandonment harbours significant uncertainties where the erosion rate of strath-top deposits is unknown , and OSL dating in the Dades has targeted the few available pockets of overbank sands covering the top of these deposits (Stokes et al., 2017). In contrast, this study focuses on dating the full aggradational history of Quaternary strath terraces along the Mgoun River, situated upstream of the mountain front between these two previously dated sites (Fig. 2), using an approach to OSL burial dating of gravels. Strath terraces along the Mgoun River are characterised by low relief (less than a metre) strath surfaces in weak red mudstones and siltstones forming valley bottoms, varying from ~ ten meters to hundreds of meters in width (Fig. 3), covered by gravel and sand fill derived from the headwaters which erode primarily through limestones and red sandstones.

Methods and approaches
We mapped and recorded the sedimentary characteristics of strath deposits (Figs. 2, 3). To document the stratigraphy and elevations of strath surfaces and deposits we conducted field surveys using a laser range finder and GNSS, combined with analysis of the 12 m per pixel TanDEM-X supplied by the German Aerospace Centre (https://tand emx-science.dlr.de/). A letter and numbering system names terraces at each site's river reach with T1 the lowest and youngest strath terrace, with subsequent higher straths corresponding to T2 and T3 etc. We refrain from making any apriori inferences about correlations of such levels between sites. Two sites (A & B), positioned downstream and upstream, were selected for geochronology based on the ability to sample from the basal bedrock strath surface to the top of the terrace sediment cover. Both sites have a bedrock strath surface in red mudstones and siltstones. Site A exhibits a strath with two sequences of bedload gravels and overbank sands (Fig. 3A). Four gravel and three sand samples (Fig. 3A) were taken from this sequence. Samples were also taken in the gravels of a lower strath terrace at this site (Fig. S1). In site B, four samples of gravels on top of the strath were taken within a 3 m thick gravel deposit (Fig. 3C). In addition to the two main sites, two gravel samples were taken from the strath terrace at site C (Fig. 3B) upstream of site A and downstream of site B (Fig. 2) to assess the synchronicity of incision and aggradation along the length of the Mgoun River.

OSL methods
We applied conventional multi-grain single aliquot regenerativedose quartz OSL dating along with the routine reliability checks of the protocol Wintle, 2000, 2003) together with innovative methods of dose-rate (Ḋ) determination (Cunningham et al., 2018;Supplementary Material). Ḋ measurements and calculations for graveltopped strath terraces are complicated by the heterogeneity of the matrices and potentially complex water saturation histories. This study demonstrates an approach to both of these challenges, enabling a chronology to be generated throughout the strath-top fill. As gravels markedly do not conform to assumptions of a homogeneous infinite matrix (Aitken, 1985), the Ḋ contributions of groups of grain sizes need to be explicitly considered. The effective ionising ranges of alpha, beta and gamma radiation emitted by grains are ~10 μm, ~2 mm, and ~ 30 cm respectively ( Fig. 4; Aitken, 1985). The alpha component can be safely ignored because of the HF etch that removes the outer rim of sand grains during sample preparation. As most sand grains are in a matrix that is several mm wide in between cobbles and pebbles (Fig. 4), beta radiation comes primarily from surrounding sand grains, whilst gamma radiation integrates the bulk of sand and pebble/cobble material. To account for grain size-dependent Ḋ variations four samples with at least 500 g of Ḋ material were wet-sieved through 500 μm, 1-, 2-, 4-, and 20mm sieves and prepared for gamma spectrometry and beta counting. We then calculate a beta dose rate correction factor which we apply to all other samples with bulk sediment dose rates. By quantifying dose rate contributions from several grain size populations, we estimated the Fig. 3. Valley cross-sections of the three sites organised from upstream to downstream, with sections facing downstream direction. Sampled terraces are highlighted by arrows. Section logs and photos include the luminescence sample locations indicated by white and black arrows. A) Site A has three strath terrace levels, with a wide strath terrace T2 in red mudstones and siltstones bordered by conglomerates. Slope material from the conglomerates lap onto the fluvial gravels of T2. T2 at site A has an intercalation of sand and gravel material sampled from bottom to top. B) Site C has unpaired strath terraces on weak red mudstones and siltstones, with valley walls consisting of strong limestones. C) Site B exhibits a staircase of strath terraces >1 km wide, overlying weak mudstones and siltstones dipping in the direction of the modern river. T2 on the other bank is perched on a mix of mudstones and limestones.
overall dose rate to the measured grain sizes from the gravel deposits ( Supplementary Fig. S6). A total of 11, 6 and 2 OSL ages were obtained from sites A, B, and C respectively (Table S1).
In addition to the grain-size effect on dose rates, estimation of water content is important for Ḋ calculations as moisture content affects the absorption of dose into grains. Gravels of river strath terraces have two main stages of water content during their history of aggradation and abandonment: saturation during perennial river aggradation, and nonsaturation after incision and abandonment. Typically, both "natural or as-found" and saturated water contents are measured in the lab and both values are taken into account for estimating the life-time average water content. The natural moisture content from the field sample, together with the saturated moisture content measured in the laboratory (supplementary materials Table S1) form inputs for Ḋ calculations using an iterative python script (https://github.com/jessezondervan/corr_lum_ dose_rate, Fig. S4). Where abandonment of the strath terrace is constrained by a sample towards the top of the sequence, the lengths of time underlying samples were saturated before terrace abandonment are calculated. Iteration of resulting ages continues until a correction of <2 ka is reached (Table S1, Fig. S4).
Finally, to test the stratigraphic consistency of age data and to restrict the distribution of uncertainties, we analyse the age estimates using Bayesian statistical methods (e.g. Rhodes et al., 2003). Bayesian analysis is a quantitative method for taking into account prior information such as the stratigraphic order of samples. Where samples are relatively closely spaced chronologically, and where age estimates overlap, stratigraphic relationships and Bayesian methods can be used to provide additional resolution. Here, we use the approach described in Rhodes et al. (2003) and in the supplementary material, using the radiocarbon calibration programme OxCal v 4.3 (Bronk Ramsey, 2009) to provide the final chronology (see Supplemental Material Fig. S8, Table S1). The Bayesian mode provided an agreement index of 79.9 % and 73.4 % for site A and B respectively, with a value of 60 % considered acceptable (Rhodes et al., 2003). The final ages corrected for water saturation history and after Bayesian analysis are presented in Table S1 and used in Figs. 5 and 6.

Results
Two distinct trunk river strath terrace levels (T1 and T2) are identified along the length of the Mgoun River in our study area with at least five levels (T1-5) present locally (Fig. 3). Terrace deposit thicknesses range from 2 to 10 m, and overbank sands cap bedload gravels on T2 at site A. Basal strath heights range from ~10 to ~40 m above the modern channel (Fig. 3).
A combination of sedimentology, the relative elevations of terraces and OSL burial dates allows us to constrain periods of aggradation and river incision (Figs. 5, 6). Burial dates of overbanks sands, such as those found along the Dades River (Stokes et al., 2017), record the start of river incision and terrace abandonment when the river channels have a reduced sediment supply or discharge to transport coarse sediment, but where floodwaters can still overtop incised terrace surfaces (e.g. Wang et al., 2014). These sands on top of gravel deposits are laterally continuous for tens to hundreds of meters across strath tops at a constant elevation, which point to widespread terrace abandonment rather than lateral channel migration. This explanation is also observed in the modern Mgoun and Dades River flood hydrology where flood waters overtop valley floor terraces depositing sand and silt. More detail on overbank sands and their relation to terrace incision in the southern High Atlas can be found in Stokes et al. (2017).On the other hand, burial dates of bedload gravels record periods of riverbed aggradation (Stokes et al., 2017;Zondervan, 2021). The Mgoun River terrace stratigraphy and OSL geochronology combine to convey punctuated incision histories at sites A and B (Figs. 5, 6). Over the last 200 ka, site A records the start of strath planation and aggradation of T2 at ~180 ka, which is then incised and abandoned at ~100 ka, followed by the T1 terrace which aggrades from ~70-80 ka and is incised after 50 ka (Figs. 5, 6). Conversely, site B records the start of strath planation and aggradation of T2 at ~130 ka and incision and abandonment at ~60 ka (Figs. 5, 6). Of these, only the planation of T2 of site B at the start of the marine oxygen isotope stage (MIS) 5e interglacial maximum (~130 ka) could be correlated to changes in climate recorded in the marine oxygen isotope and NW African humidity records (Fig. 6). Overall, the periods of river incision and abandonment of strath levels T1-2 across the studied sites do not convincingly and consistently correlate with climatic stages or transitions (Fig. 6).
The gravel burial dates throughout the strath deposits allow not only a constraint on periods of river incision and abandonment, but also the pulses of gravel aggradation within each terrace deposit. In site A, the first accumulation event starts at 170 ka and is followed by a major accumulation event from 124 to 90 ka (Figs. 5, 6). At site B, accumulation of ~1 m of gravel took place around 128 to 115 ka, and another ~1.5 m around 74 and 59 ka (Figs. 5, 6). When the gravel bedload burial ages throughout the strath deposits from sites A, B and site C, upstream of site A and downstream of site B (Fig. 2), are plotted against time (Fig. 6), we observe a synchronous pulse of terrace bedload aggradation during MIS 5e. Thus, whilst the timing of the start of strath planation and aggradation, or periods of river incision and terrace abandonment between the three sites do not yield a systematic relationship to Quaternary climate (see previous paragraph), through careful analysis throughout the strath terrace age-model we do find a climatic signal in the form of pulses of aggradation within prolonged periods of strath-top aggradation.
When the elevation data is combined with the age constraints on aggradation and incision, river incision rates can be constrained and compared to rock uplift rates. The incision rate between T1-T2 aggradation at site A is at least 0.95 mm yr − 1 (Fig. 5a), which is up to an order of magnitude larger than the lowering rate ~ 0.22 mm yr − 1 which is

Fig. 4.
Diagram showing texture of sampled gravels and sample homogeneity of alpha, gamma and beta dose rates. Each radiation type has a different penetration depth (Aitken, 1985). Alpha radiation has the shortest penetration depth and in sand size grain fraction influence from surrounding grains affect only on the outer rim which is dissolved during preparation. Since most sand grains lie within the inter-clast matrix, beta radiation comes from the surrounding sand grains mostly. Gamma radiation penetrates to around 30 cm, and thus comes from both the surrounding sand grains and larger clasts. derived by integrating aggradation and incision over the last ~200 kyr (Fig. 5a). This long-term base-level lowering rate is consistent with High Atlas long-term rock uplift rates of 0.17-0.22 mm yr − 1 . If rates were calculated from the abandonment ages of terraces T2 and T1 this would result in apparent rates of 0.58 and 0.33 mm yr − 1 (Fig. 5b), which are values lower than the punctuated incision rates and higher than the time-averaged base-level lowering rate (Fig. 5a). At site B, rates based on abandonment ages to the modern riverbed reach at least 0.22 mm yr − 1 and time-averaged rates -including aggradational time -reach 0.09 mm yr − 1 (Fig. 5b). However, considering the position of the modern riverbed in the cycle of aggradation, the age data from both sites enable a fit of long-term (10 5 yr) base-level lowering rates in equilibrium with rock uplift rates of 0.17-0.22 mm yr − 1 (Fig. 5b). The time between the T2 and T1 burial ages might contain unpreserved periods of planation or aggradation in addition to vertical incision. Thus, whilst the results in this study constrain a maximum length of time for vertical incision, the duration of incision could be even shorter.

Discussion
The strath terrace stratigraphic age model constrained using OSL burial dates presents periods of up to 80 kyr between strath formation and incision, with these incisional periods lasting 20 kyr or less (Figs. 5,6). Consequently, this age model records a punctuated river incision history. Furthermore, the timing of strath terrace incision varies between locations, and the time between strath incision varies from 30 kyr to up to 80 kyr (Figs. 5, 6). Thus, the terrace chronology here constrains a mode of strath formation that does not faithfully follow periodic cycles exhibited by the climatic history of the region (Fig. 6). Hence, our strath terrace age model constrains a metastable equilibrium model of punctuated river incision and strath formation presented in Fig. 1B. In the following sections we combine the punctuated river incision model with our well-constrained terrace formation history to demonstrate pitfalls and recommended strategies of numerical dating of strath terraces and their interpretation of climatic signals, rock uplift and incision rates. Subsequently, we demonstrate how the punctuated river incision model presented here can aid in the reinterpretation of sites with limited age data to arrive at rates that are consistent across a mountain belt experiencing punctuated river incision.

Determining climatic control on landscape evolution
The key to understanding punctuated incision and its importance when deriving rates is the relationship between river incision and the Quaternary periodicity in climate. In the gradual incision endmember (Fig. 1A), in instances where strath terraces form in response to climate, these form through periodic cycles of incision forced top down by the periodicities of Quaternary climate (e.g., Hancock and Anderson, 2002;Bridgland and Westaway, 2014;Gallen et al., 2015). However, the Quaternary incision history recorded in the terrace sites of this study is asynchronous and does not relay a systematic relation to climatic cycles: neither precession-driven climate periods nor eccentricity-paced interglacial-glacial cycles (Fig. 6). Such a result appears consistent with a study of terraces of the western High plains in Colorado, formed under punctuated river incision, which revealed a strath terrace which was left un-incised over hundreds of thousands of years, suggesting that strath terrace formation is not always directly related to Quaternary climate cycles (Foster et al., 2017). Similarly, Schanz et al.'s (2018) review of strath terrace chronologies showed that strath terraces globally do not necessarily form during any particular climatic stage. Whilst these studies demonstrated that strath terraces may form asynchronously across a mountain belt or between sites globally, we here confirm that such asynchronicity can be expected within sites along the same river too. Thus, one correlation between strath formation and a climatic stage at one site along the river does not guarantee a climatic origin for such a terrace: even if formed through climatic forcing, direct correlations between a climatic stage and timing of one site's terrace-formation might not be causative, especially when responses to climate are asynchronous along a river's length. In addition, we have an OSL burial age model throughout deposits from strath to fill-top, unlike the previous studies which had age control on terrace abandonment and in some cases the start of strath formation.
Whilst the history of incision constrained by the start and end of strath terrace occupation does not reveal a climatic control on the landscape evolution of the High Atlas, a synchronous pulse of riverbed aggradation appears in the burial dates throughout the strath deposits from multiple sites (Fig. 6). This pulse of strath aggradation across the sites is dated to the glacial maximum to interglacial maximum transition occurring at an eccentricity maximum ~130 ka (Fig. 6). This points to the indirect role of climatic control on sediment generation and transport in pacing strath terrace formation in this mountain belt experiencing punctuated river incision. For example, in the High Atlas, it is likely that the effect of increased frequency and magnitude of tropical storm incursions in the southern High Atlas at interglacial maximum times has influenced the coupling of lateral sediment sources to the main stem of the Mgoun River (Zondervan, 2021).
Our results demonstrate that strath terraces formed under punctuated river incision may not reveal a signal of climatic control on landscape evolution unless burial dates are obtained throughout deposits from multiple sites along a river's length. This requires a rigour of geochronology that far exceeds those of rock uplift and incision rate calculations, as explained in the following section. However, the punctuated river incision model presented here (Fig. 1B) frees workers from the urge to box numerical dates and strath formation into climatic stages. This is important for when these terraces are used to constrain rates of rock uplift and river incision.

Calculating incision and rock uplift rates
The history of punctuated incision and long periods of terrace aggradation constrained using stratigraphy and chronology (Fig. 5) illustrates the potential pitfalls in calculating rates of river incision and base-level fall when not considering the position of dates within the cycle of aggradation and incision. At site A, rates calculated between terrace abandonment times were lower than incision rates as calculated between the terrace abandonment and start of terrace aggradation on the next level (Fig. 5). Nor did the rates calculated between terrace abandonment times reflect long term rock uplift rates, derived here from samples dating the start of aggradation on each terrace, taken just above the strath surface, consistent with independently constrained rate of uplift of 0.17-0.22 mm yr − 1 . Such results demonstrate the problem with previous recommendations of rate calculations in setting with punctuated river incision. In the Clearwater River basin of the Olympic Mountains, Wegmann and Pazzaglia (2002) found that over the last 11 kyr, vertical incision was relegated to brief 1 kyr intervals, with the rest of the record dominated by aggradation on strath terraces. They pointed out that incision rates reconstructed using the youngest Holocene terraces were higher than those calculated from Pleistocene terraces, because the latter averaged over periods of glacial-interglacial aggradation and incision. Consequently, Wegmann and Pazzaglia (2002) recommended integrating rate calculations over terraces representing at least one glacial-interglacial cycle. In a similar fashion, Gallen et al. (2015) proposed that rather than the modern riverbed, the lowest strath terrace should instead be used as the reference frame to calculate incision. Site B illustrates the pitfall in using the lowest strath level to constrain a rate of rock uplift (Fig. 5B), as the elevation of the modern riverbed cannot be used as a point to calculate rock uplift rates, as demonstrated by Gallen et al. (2015). However, Gallen et al.'s (2015) approach to calculating rates from multiple terrace levels only works if terrace formation follows a consistent periodicity, which is an assumption that does not always hold in settings with punctuated river incision (see section 6.1). Indeed, at site A rates calculated between terrace abandonment times reflect neither long term rock uplift rates nor river incision rates (Fig. 5A). The limits of Wegmann and Pazzaglia's (2002) and Gallen et al.'s (2015) approaches were previously recognised by Schanz and Montgomery (2016): in the Willapa and Nehalem river basins of the Pacific Northwest, these authors used their radiocarbon record of a 10,000 year break in river incision represented by terrace deposits to point out that even using the method proposed by Wegmann and Pazzaglia (2002) and Gallen et al. (2015) would still result in erroneous uplift rates. Nonetheless, Schanz and Montgomery (2016) were not able to demonstrate strategies to find accurate rates, limited by a lack of preservation and the time range of radiocarbon dating. We can here demonstrate that the limitations of the approaches recommended by Wegmann and Pazzaglia (2002) and Gallen et al. (2015) originate from the lack of differentiation between rates of river incision and rates of rock uplift, and that these rely on the assumptions that terrace formation is consistently periodic, and that the vertical extent of aggradation on each level is similar (Figs. 1, 5). Furthermore, combining the metastable punctuated river incision model with our 150 kyr strath terrace allows us to demonstrate and recommend approaches to deriving accurate rates even when such assumptions break down. These approaches differ depending on which rate is desired: river incision or rock uplift. For rock uplift rate constraints, we advise basal burial dates from just above the strath surface, at a minimum of two strath levels, and to leave the modern riverbed out of the calculation (Fig. 5). On the other hand, if river incision rates are to be constrained, a burial date or surface Fig. 6. A) OSL results as points and summed probability distributions for aggradation on terraces at three sites with their stratigraphic context. Terrace levels, timing of terrace aggradation and incision are labelled. B) The marine oxygen isotope record, based on the LR04 benthic δ 18 O stack (Lisiecki and Raymo, 2005). C) Quaternary climate records of north-western Africa: Green Sahara periods (Larrasoaña et al., 2013), Saharan humidity index from borehole GeoB7920 (Tjallingii et al., 2008), and the δ 18 O record of relative rainfall (Dixit et al., 2020). Synchronous strath aggradation correlated with the interglacial maximum stage 5e is marked. exposure date from the top of the strath deposit and a basal burial date just above the strath surface at the next terrace down are needed (Fig. 5). Whilst an abandonment age of the lowest terrace to the modern riverbed can reveal a river incision rate that is significantly higher than rock uplift, and thus signal a punctuated river incision history (Fig. 7), site B demonstrates how the time between terrace abandonment and the modern river plain may also include a period of aggradation on the modern valley bottom and thus result in a lower apparent rate that may hide the punctuated nature of river incision (Fig. 5). Thus, both reliable rock uplift and river incision rate constraints require the use of burial dating in settings experiencing punctuated river incision. Understanding the difference between constraining rock uplift and river incision rates (Figs. 1, 5) is also important when workers attempt to compare river incision rates to geomorphic measures such as channel width, slope and metrics of sediment character and transport (e.g. Brocard and Van der Beek, 2006;Allen et al., 2013;Zondervan et al., 2020aZondervan et al., , 2020b.

Calculating rates with limited constraints
Whereas the High Atlas provides terraces with well-preserved strathtop deposits, not all settings will be candidates for basal burial and abandonment dates to constrain a punctuated river incision history. Furthermore, not all published terrace studies report the age constraints recommended by us, which can be a challenge when attempting to derive a tectonic or climatic history for a region. However, even with limited age constraints, the punctuated river incision model can be used to interpret terrace ages. The process of interpretation is demonstrated visually (Fig. 7), and requires knowledge of the distance between the top of the terrace deposit and the strath surface and either an independent constraint on rock uplift rates, or an estimate for the river incision rates experienced during periods of punctuated incision, or both (Fig. 7). For example, west of the Mgoun, the interpretation of abandonment ages of terraces of the Madri River in the Ouarzazate foreland basin resulted in incision rates of 0.3-5.2 mm yr − 1 (Fig. 7), seemingly increasing towards the present . However, using the constraints on strath-top sediment thickness and accounting for punctuated incision, long-term rates of base-level fall of 0.3 and 0.5 mm yr − 1 can be calculated for the two profiles (Fig. 7). These rates, derived from sites <3 km apart along stream, include additional structural rock uplift and local relative base-level fall related to slowly moving active blind thrusts in the western Ouarzazate basin which account for 0.1-0.3 mm yr − 1 (Pastor et al., 2012). Similarly, dated abandonment ages of terraces east of the Mgoun, along the Dades River, fit a background base-level fall rate of 0.17 mm yr − 1 (Stokes et al., 2017) consistent with regional rock uplift rates across the High Atlas, when accounting for punctuated incision. Thus, a reinterpretation of previously published age-height data using our model of punctuated incision in the High Atlas provides a consistent understanding of regional low-rate rock uplift and base-level fall. Conversely, the dashed lines in Fig. 7 illustrate how without accounting for punctuated incision, the derived rates of incision from terrace abandonment ages would have produced a Sadler effect, which without due attention could have been interpreted as an increase in tectonic uplift rates. Fig. 7. Age-height data from Arboleya et al. (2008) along the Madri River in the foreland basin west of the Mgoun, plotted as two profiles in the punctuated incision framework. Series a-c guides the interpretation of limited age-height data, illustrating how data that seemingly demonstrates an increased rate of incision to the present can be interpreted to present a punctuated river incision history instead, with a constant rate of base-level fall. (d) In the example of the Madri, reinterpretation based on a punctuated river incision history demonstrated in our detailed sites leads to rates of base-level fall that are consistent with regional and local rates of uplift. Profiles 1 and 2 present sites upstream of an active thrust fault with 0.1-0.3 mm yr − 1 of structural uplift (Fig. 2, Pastor et al., 2012).

Where to use the punctuated river incision model for strath terrace formation
Settings in which punctuated river incision is expected, as opposed to gradual incision, are settings which have rock erodibility or rates of base-level fall which allow for long-term strath bevelling and aggradation, such as low-uplift-rate, wide-valley, low-rock-strength settings . Climatic periodicity within the fluvial record may also be stronger in glaciated catchments, and large downstream river systems (e.g. Bridgland and Westaway, 2014), where the stream power and transport capacity of the trunk stream dominates over any complexities from lateral input from hillslopes or riverbanks (Zondervan, 2021).
We recommend that terrace data is interpreted in the context of other regional constraints: it is more likely to find proof for punctuated river incision within a strath-top deposit, whilst proof that a river has experienced no punctuated incision and only gradual incision is more difficult to find because such a history would not generate or preserve extensive strath-top sediments. Consequently, most records could be interpreted as punctuated incision where there is an absence of sufficient constraints. Thus, factors such as the erodibility of bedrock, valley width and uplift rates which might affect the style of river incision and terrace planation should be considered (e.g. Schanz and Montgomery, 2016), and rates between sites should be compared to find the most consistent incision history.

Conclusions
By highlighting the often-made assumptions in deriving incision rates from strath terraces, we show how interpretations can easily arrive at measured time-dependent and inconsistent incision rates in settings with punctuated river incision. The punctuated river incision model for strath terrace formation over Quaternary timescales, illustrated by our 150 kyr history of strath terrace formation in the High Atlas Mountains, allows for a consistent understanding of Quaternary river incision rates, rates of base-level lowering, and rock uplift across the High Atlas. This model, even when applied to limited terrace age data, has solved for the Sadler Effect and reconciled rates between studies across the mountain belt. Our punctuated incision model allows for strath formation that is not necessarily cyclical or directly correlated with climatic stages. The 150 kyr Quaternary terrace age-model constrained here shows that in the High Atlas, strath terrace abandonment and river incision are asynchronous between sites and not convincingly correlated with climatic stages, but aggradation experiences a pulse during the interglacial maximum at MIS 5e. Whilst a history of punctuated river incision complicates the interpretation of strath terrace records for constraints on rates of incision and rock uplift, combining the punctuated river incision model and our recommendations on sampling and interpretation of strath terrace age data allows for more consistent and accurate insights into the climatic control on erosion and landscape evolution, Quaternary river incision and rock uplift rates.

Funding
The IAS Postgraduate Grant, BSG Postgraduate Research Grant, and the QRA New Research workers Award supplemented J.R.Z.'s University of Plymouth PhD scholarship. We also acknowledge funding in kind from the German Aerospace Center (DLR) in the form of the TanDEM-X data.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data and code are shared in the supplementary files and links therein.