Low-temperature thermochronology as a control on vertical movements for semi-quantitative Source-to-Sink analysis: A case study for the Permian to Neogene of Morocco and surroundings Authors

Continental passive margins and their hinterlands in the Atlantic realm have been the locus of many Low Temperature Thermochronology (LTT) and time-Temperature (t-T) modelling studies that evidence pre-, syn-, and post-rift episodic km-scale exhumation and burial episodes. In this study we integrate data from over thirty published LTT and t-T modelling studies from Morocco and its surroundings using a 3-step workflow to obtain 1) exhumation/burial rates, 2) erosion rates, and 3) paleoreconstructions of source-to-sink domains, between the Permian and the Present. Our synthesis of available t-T modelling results predicts high exhumation rates in the Anti-Atlas (0.1 km/Myr) during the Early to Middle Jurassic, and in the High Atlas (0.1 km/Myr) and Rif (up to 0.5 km/Myr) during the Neogene. These rates are comparable to values typical of rift flank, domal or structural uplifts settings. During the other investigated periods, exhumation rates in the Meseta, High-Atlas, Anti-Atlas and Reguibat shield are around 0.04±0.02 km/Myr. Interpolation of the exhumation rates at the regional scale allow calculation of the volume of rocks eroded. Estimates of erosion rates are between 0.2x103 and 7.5x103 km3 (in the Meseta and the Reguibat Shield, respectively). Ten erosional (quantitative, from interpolation results) and depositional (qualitative, from data synthesis) “source-and-sink” maps have been constructed, with emphasis on the Jurassic and Cretaceous periods. The maps integrate the extent of exhumed domains, using information from geological maps, lithofacies and biostratigraphic data from new geological fieldwork and well data from onshore and offshore basins. The results illustrate changes in the source-to-sink systems and allow for a better understanding of the Central Atlantic margin hinterlands evolution.


Introduction
Continental passive margins and their hinterlands ( fig. 1) are the locus of a significant amount of Low-Temperature Thermochronology (LTT) studies that evidence pre-, syn-, and post-rift episodic kmscale upward (i.e. exhumation) and downward (i.e. burial) movements (e.g. Gallagher et al., 1994;Green et al., 2018). Such movements are an important component of source-to-sink systems (e.g. Helland-Hansen et al., 2016;Bhattacharya et al., 2016). A number of exhumation episodes are recorded on the margins of the Atlantic Ocean. Pre-rift movements are well documented (e.g. Juez-Larre and Andriessen, 2006;Ruiz et al., 2011;Jelinek et al., 2014;Japsen et al., 2016a). Syn-rift exhumation is observed along the Atlantic rift flanks (e.g. Oukassou et al., 2013;Jelinek et al., 2014;Wildman et al., 2015;Japsen et al., 2016a and several hundreds of kilometres into the plate interior (e.g. Leprêtre et al., 2013). Syn-rift burial episodes in the unstretched continental crust have also been documented, although in fewer places (e.g. Juez-Larre and Andriessen, 2006;Ghorbal et al., 2008;Gouiza et al., 2019). Lastly, postrift km-scale exhumation/burial have been documented along the North (e.g. Japsen et al., 2006;Japsen et al., 2016a;2016b), Central (e.g. Frizon de Lamotte et al., 2009;Bertotti and Gouiza, 2012;Amidon et al., 2016) Equatorial (e.g. Hayford et al., 2008;Ye, 2016;Wildman et al., 2018) and South (e.g. Jelinek et al., 2014;Wildman et al., 2015) Atlantic margins. Vast regions along Atlantic rifted continental margins expose pre-rift rocks, which can be examined using LTT and time-Temperature (t-T) modelling techniques to provide unique understanding of the thermal history. Low-temperature thermochronology ages record the cooling of rock samples, which can be caused by either thermal relaxation and/or exhumation (also called denudation when erosional in nature; e.g. Pagel et al., 2014;Malusà and Fitzgerald, 2019b). This is especially valuable for geologically ill-constrained areas characterised by no or little sedimentary cover (e.g. Gallagher et al., 1998;Ghorbal et al., 2008;Japsen et al., 2009;Teixell et al., 2009;Japsen et al., 2012b;Jelinek et al., 2014). Results of these techniques can be used as proxies to reconstruct exhumation/burial events (e.g. Teixell et al., 2009) and provide a key constraint on sediment source (e.g. Helland-Hansen et al., 2016). Bertotti and Gouiza (2012) proposed that unpredicted km-scale exhumation episodes are concurrent with excessive downward movements in the subsiding domain along the Moroccan Atlantic margin. This relationship was further tested and documented in different segments of the same margin . Exhumation and burial episodes occur at the same time in regions characterised by non-stretched lithosphere, demonstrating that processes other than rifting are at work, or that the effects of the rifting and drifting extend beyond the rifted margin and their direct flanks (e.g. Gouiza et al., 2017;. Several authors have tested aspects of these km-scale movements with numerical models (e.g. Leroy et al., 2008;Gouiza, 2011;Cloetingh and Burov, 2011;Yamato et al., 2013). However, to better constrain the models, quantification of exhumation/burial events over geological time and at the scale of the margin and beyond, is required (e.g. Ye et al., 2017;Wildman et al., 2019). In this study we integrate data from over thirty LTT and t-T modelling studies published in Morocco and its surroundings (see supplementary file 1), including detrital LTT. The majority of the derived LTT ages span the period between the Variscan and Atlas orogenies (ca. 300-40 Ma;Charton, 2018) and the Atlasic inversion event (e.g. Frizon De Lamotte et al., 2009;ca. 40-0Ma). These cooling ages were interpreted to record vertical movements relative to Earth's surface (e.g. Ghorbal et al., 2008;Sebti et al., 2009;Leprêtre et al., 2013;Gouiza et al., 2017;Charton et al., 2018), although currently the mechanism that generates the vertical movements within the passive margin hinterland remains enigmatic (e.g. Ghorbal et al., 2008). The onset of Pangaea break-up, marking continental drift and initiation of the proto-Central Atlantic, occurred at the beginning of the Jurassic, between 190 and 170 Ma (e.g. Davison, 2005;Labails et al., 2010;Sibuet et al., 2012). The development of the Moroccan passive margin during the Jurassic and Early Cretaceous witnessed the accumulation of neritic and deeper marine sediments offshore, and continental to paralic deposits in the coastal basins, the High/Middle Atlas, and in the Meseta basins. The Peri-Atlantic Alkaline Pulse (PAAP) is recorded by plutons and flood basalts in the conjugate margins of the Central and South Atlantic oceans, between 125 and 80 Ma (Matton and Jebrak, 2009;Montero et al., 2016). Hollard et al., 1985). CAMP: Central Atlantic Magmatic Province (Triassic); PAAP: Peri-Atlantic Alkaline Province (Cretaceous). B) Tectonic chart of the Atlantic realm and Morocco, see references therein (numerical ages after ICS 2016/12). LIP: Large Igneous Province.

Figure 2 | A) Geological map of Morocco and its surroundings (after
The opening of the South Atlantic at ca. 135 Ma (Collier et al., 2017) led to African-European plate convergence and sustained far-field intraplate stresses. The related Atlas Orogeny peaked in the Eocene, and is characterised by exhumation and deformation of the Middle and High Atlas ( fig. 1A; Guiraud, 1998. Finally, Cenozoic volcanism (Missenard and Cadoux, 2011) and surface uplift, observed along two axes, from the Canary Islands to the Siroua massif (Anti-Atlas) and from the Siroua Massif to the Rif belt (figs. 1 and 2), have been interpreted as associated to a mantle anomaly Teixell et al., 2005;Missenard, 2006). Although, such interpretation has been questioned by more recent work (e.g. Ellero et al., 2020). Lanari et al. (2020) suggest that the crustal deformation controlled the exhumation pattern and uplift on a short scale, with coexisting mantle dynamics that control both large wavelength topography and short wavelength crustal deformation.
3. From time-Temperature models to Source-to-Sink maps: workflow, intermediate results, and uncertainties 3.1. Key steps of the workflow Step 1) From time-temperature (t-T) models to exhumation/burial rates (part 3.2); S1: sample 1.
Step 3) From exhumation maps to source-and-sink maps (part 3.4).
The workflow, illustrated in figure 3, is composed of the following three steps: Step 1.
Selection of t-T models based on review of geological constraints and other conditions (table 1); 1.3.
Temperature-to-depth conversion using varying and constant geothermal gradient from present-day analogues (figs. 6 and 7), assuming that cooling means erosional exhumation and heating means burial; 1.4.
Divide investigated time into periods according to main trends in t-T models and stratigraphic boundaries (fig. 8); 1.5.
Obtain slopes for t-Depth models based on start and end of previously defined periods, which assumes that a strong exhumation followed by burial in one period is similar to an overall milder exhumation (as it cancels out peaks within each period); 1.6.
Plot rates and ranges, based on the depth conversion using varying and constant geothermal gradients, respectively (fig. 8).
Interpolation inputs for exhumation maps: i1. Compilation of stratigraphic columns to use as constraint in the interpolation ( fig. 9), assuming that sediment preserved in a basin means burial for the related; i2. Localisation of input data with a GIS program. Digitalization of active faults for defined periods from paleogeographies for instance, which may then be used as interpolation limits. This assumes that all faults have similar effect on the exhumation rates; 2.2.
Interpolation of data using the "Inverse Distance Weighted" (IDW) algorithm in the Surfer program ( fig. 10). Using this method means that we assume a linearly weighted variation of exhumation/burial rates between points; 2.3.
Volumetric calculations are performed using the Surfer program, for values higher than 0km/Myr, which yields the amount of material which has been removed per million year (erosion rate in km 3 /Myr; fig. 11 and table 2). Clipping polygons are used to estimate the contribution of each considered area. The range is given by running the interpolation and the volume calculations using rates obtained with constant geothermal gradients.
SandS mapping inputs: i3. Surficial geology from the geological map ( fig. 2A; Hollard et al., 1985) with associated depositional environment when described; i4. Well data (and location) with associated or interpreted depositional environments from literature or company reports ( fig. 12); i5. Outcrop and fossil data with associated or interpreted depositional environments from the literature ( The detail of the method of each step of this workflow and methods for calculation of exhumation/burial rates and generation of exhumation maps is summarized in the following sections, which include a brief description of the intermediate results, and for each of the steps, the assumptions and limitation.

3.2.
Step 1: From time-Temperature (t-T) models to exhumation/burial rates 3.2.1. t-T modelling database In Morocco and its surroundings, over thirty LTT published studies have been conducted in the last 30 years (e.g. Mansour, 1991;Malusà et al., 2007;Ghorbal et al., 2008;Ruiz et al., 2011;Lanari et al., 2020; see supplementary file 1). This provides an extensive database of over 1000 LTT ages including; (U-Th)/He ages from apatite crystals (AHe), Apatite Fission Track (AFT) ages; (U-Th)/He ages from zircon crystals (ZHe), and Zircon Fission Track (ZFT) ages. Exposed Precambrian crystalline basement rocks, Palaeozoic metapsamites within an otherwise marine metapelite dominated column, Meso-Cenozoic clastic sediments, and dykes/sills of all ages were the most sampled lithologies. LTT studies often use the calculated ages, fission track density and length as inputs for t-T inverse modelling (e.g. Pagel, 2014;Malusà and Fitzgerald, 2019a). This inverse modelling technique allows testing of several t-T paths by guiding model realisations with user-defined constraints. This method provides a comprehensive assessment of the thermal history of the analysed sample. Twenty published studies had performed and published t-T modelling, resulting in 141 t-T models (figs. 4 and 5; table 1). The programs that were used in these studies for the inverse modelling of LTT data are: HeFTy (Ketcham, 2005), AFT Solve (Ketcham et al., 2000), and QTQt (Gallagher, 2012; also, see Vermeesch and Tian (2014) for a comparison of HeFTy and QTQt codes), which provide comparable thermal histories. The program outputs are labelled as 'acceptable', 'good', 'best-fit', or 'weighted average' (referring to a goodness of fit) paths for HeFTy/AFTSolve and 'probability' range, 'maximum likelihood' or 'expected' paths for QTQt. Weighted average/expected curves were digitized when available, or if weighted average was not available, the best-fit/maximum likelihood curves, using WebPlotDigitizer (Ankit Rohatgi; https://automeris.io/WebPlotDigitizer/). All digitized curves for each study are available in the supplementary file 2.  table 1). Note that the location of t-T models from Logan and Duddy (1998) in the Reganne basin was approximated and Leprêtre (2015) conducted one t-T modelling in the Western Anti-Atlas and all others in the Reguibat Shield. Table 1 | Selection of t-T models for temperature-to-depth conversion and for exhumation and subsidence rate calculations. Conditions: The t-T modelling results had to i) start before 20 Ma, ii) be based on HeFTy, AFT solve or QTQt results and display a favoured curve, as opposed to envelopes only, iii) be from the most recent realizations, if different models using the same LTT data exist, iv) be based on single rock sample as opposed to vertical profiles and other combined modelling, or from shallowest sample from borehole data, and v) be compatible with valid geological constraints as discussed in the text. * In the case of Sehrt (2014) and Lepretre (2015) theses, more recent publications of Lepretre et al., 2017 andSehrt et al., 2017; reused the same t-T models and are thus not listed here.
Figure 5 | Digitized t-T modelling weighted averages for considered regions (A-I), best-fit curves, and acceptable envelopes (grey area) for HeFTy results or expected curves, maximum likelihood curves and limits of 2σ confidence level (grey area) for QTQt results (see references in table 1); Central M.: Central Massif; AA: Anti-Atlas; RS: Reguibat Shield B.: Basin

Geological Constraints
In the study area, two types of geological constraints are used in considered t-T modelling studies: 1) sediments of known stratigraphy that overlay the sampled basement, with or without an erosion gap, provide a constraint to surface temperatures (ca. 10-40°C) at the time of deposition. In the case of detrital LTT, the stratigraphic age of the sampled sediments is used as a time constraint for surface temperatures; 2) when the age and temperatures of emplacement of igneous body or metamorphism of the sample are known, a corresponding radiometric age constraint is applied to the subsurface temperatures.
In the Rif belt, 40 Ar/ 39 Ar and K-Ar radiometric dating of Pliocene and Miocene sediments were used as constraints (Romagny et al., 2014;Azdimousa et al., 2013). In the Meseta, the age of Variscan granite emplacement and the Permian, Triassic, and Cenomanian sedimentary record in the basins surrounding the Variscan massifs were utilised in the t-T modelling (Saddiqi et al., 2009, Ghorbal et al., 2008Barbero et al., 2011). In the Variscan High Atlas (Massif Ancien), Triassic and poorly dated Lower Cretaceous sediments overlying Precambrian basement rocks provided the surfacetemperature age constraint (Ghorbal, 2009;Balestrieri et al., 2008;Barbero et al., 2007;Domènech et al., 2016). In the central High Atlas, emplacement timing of Jurassic intrusives was used. (Barbero et al., 2007). In the Tarfaya basin, Sehrt et al., (2017b) used the Aptian and Albian stratigraphic age of sampled sediments in the t-T modelling. In the Reguibat Shield, a geological constraint at surface temperatures was defined for the Early/Middle Cretaceous, from sediments (poorly dated) exposed in the Tarfaya and Tindouf basins (Leprêtre et al., 2013;2017). In the eastern Reguibat Shield, the poorly dated Upper Cretaceous sediments of the Reggane Basin were used to guide t-T models (Leprêtre et al., 2017). A review of these studies suggests that the use of geological constraints was applied in a fairly consistent manner in all published studies, conforming to the most-likely geological history (e.g. Michard et al., 2008). Radiometric dating of Variscan metamorphism (Ruiz et al., 2011, Malusà et al., 2007 and emplacement of the CAMP dykes (see Gouiza et al., 2017) were used to constrain the time when the sampled rocks were at great temperature/depths (i.e. higher than LTT retention temperatures). Evidence that the presently exhumed Anti-Atlas basement rocks were close to surface temperatures in the past came from observation of the relationship with overlying sedimentary packages. This included the contact with Triassic sediments (Ghorbal, 2009), Middle Jurassic sediments  in the north, poorly dated terrestrial Infra-Cenomanian sediments in the western, central, and eastern Anti-Atlas (e.g. Ruiz et al., 2011;Oukassou et al., 2013) and Cenomanian fluvial sediments overlying the Variscan basement on the eastern Anti-Atlas (Gouiza et al., 2017). However, comparison between t-T models of the Anti-Atlas indicates some major discrepancies (e.g. Oukassou et al., 2013vs Gouiza et al., 2017fig. 5C). Most of these discrepancies can be attributed to the use of different t-T modelling constraints. An example is the use of an Early Cretaceous modelling geological constraint (Early Cretaceous sediment overlying Basement) in the Anti-Atlas (Ghorbal, 2009;Ruiz et al., 2011;Oukassou et al., 2013;Sehrt et al., 2017a). Recent biostratigraphic work has re-dated much of the terrestrial to shallow marine sediment in the Tarfaya area (mapped as Lower Cretaceous) as being conclusively Bathonian (Middle Jurassic) or older Jurassic /Triassic, with implications for the modelling Arantegui et al., 2019). Finally, in the western Anti-Atlas, recent work by Arantegui (2018) presents new biostratigraphic data that conclusively dates the local redbeds (mapped as Lower Cretaceous) as Bathonian (Middle Jurassic) or older. Extensive paleontological work has also shown that the clastics overlying basement in parts of the Anti-Atlas are Cenomanian in age (Late Cretaceous; e.g. Benyoucef et al., 2015). In the central Anti-Atlas, no recent biostratigraphic study on the undifferentiated Mesozoic clastic sediments has been conducted, but the Cenomanian limestones positioned higher in the stratigraphic column can be used as a time constraint (e.g. Fetah et al., 1990). The underlying clastics could hence be Cenomanian in age, similarly to the sediments in eastern Anti-Atlas. Overall, this invalidates the use of an Early Cretaceous constraint in time-temperature modelling that was previously used to force basement rock surface temperatures in the Anti-Atlas. Moreover, it is also unclear whether Cretaceous sedimentary cover was continuous over all the Anti-Atlas or not (Malusà et al., 2007). In this work, we have used the abovementioned new modelling constraints for the Anti-Atlas, and disregard published models using surface temperatures for the Early Cretaceous as a geological constraint (e.g. Ruiz et al., 2011;Oukassou et al., 2013;Sehrt, 2014).

Temperature-to-depth conversion
In order to quantify the volume of material that was subsequently eroded above presently exposed rocks and shallowest samples from borehole data, a temperature-to-depth conversion is applied on selected t-T results (figs. 6 and 7). This conversion is often used in t-T model to estimate exhumation rates and was already applied in some of the previous studies considered in this work (e.g. Gouiza et al., 2017). One of the assumptions to apply such a conversion to t-T modelling results is that cooling must be due to exhumation, not thermal relaxation (see Malusà and Fitzgerald, 2019b for a comprehensive review). Exhumation, as defined in England and Molnar (1990), is a decrease of the distance between a considered rock and the Earth's surface. It does not mean 'uplift', as exhumation may occur by the erosion and denudation of a rock overburden, without involving motion of the rock with respect to the geoid or the undeformed reference lithosphere. The conclusion reached in nearly all t-T studies, in Morocco and surroundings, is that samples with cooling ages younger than their stratigraphic ages recorded burial and exhumation. When supported by the identification of local or regional detrital sediments of ages similar to the modelled cooling (i.e. the eroded product), it supports the use of identified cooling events as periods of exhumation. Exceptions were applied to samples that may have exhibited thermal relaxation, such as in the High Atlas (post-rift sampled Jurassic dykes; Barbero et al., 2007), in the Canary Islands (Cenozoic postvolcanism; Wipf et al., 2010) and in the Siroua massif (Cenozoic post-magmatism/volcanism thermal relaxation; Ghorbal, 2009). The interpretation of heating due to burial requires the presence of truncated sediments of coeval age in the rims of the sampled basement massif, or remnants of formerly emplaced thrust sheets that could have loaded the basement (e.g. Cerrina Feroni et al., 2010), or alternatively the absence of evidence for heating by volcanism/magmatism. Five conditions were applied to the selection of representative and valid t-T curves, resulting in a total of 68 t-T curves (~ half of initial dataset were disregarded; detailed in table 1). Published modelling results had to: i) start before 20 Ma (ca. the duration of youngest considered period), ii) be based on HeFTy, AFT solve or QTQt results and display a best fit/weighted curve, as opposed to envelopes only, iii) be the most recent realizations, if different models using the same LTT data exist, iv) be based on single rock samples as opposed to vertical profiles and other combined modelling, or from shallowest sample from borehole data, (justified by the fact that we use points for the spatial interpolation), and v) be compatible with geological constraints discussed above. To achieve the temperature-to-depth conversion ( fig. 7), different paleo-geothermal gradients were used based on the location of the selected t-T curves (fig. 6) and keeping surface temperature constant at 20°C. The applied geothermal gradients are based on several studies, including studies of similar geodynamic settings that serve as analogues for Morocco. Luth and Willingshofer (2008; see references therein) obtained geothermal gradients of 23-35°C/km for the Alps. Based on these values, we considered an average geotherm of 29°C/km for the Variscan orogeny. The geothermal gradient in the rift zone of the present day East African Rift system is ca. 40°C/km (van der Beek et al., 1998) and between 25 and 32°C/km in the Rio Gande Rift (Bridwell, 1976). We applied a geotherm of 34°C/km as representative for the High Atlas Rift zone. The flanks of the East African Rift systems today display geothermal gradients between 25 and 30°C/km (van der Beek et al., 1998) and using this as an analogue, in Morocco 27°C/km was applied for the Central Atlantic/High Atlas rift flanks. For the post-rift Moroccan passive margin c. 100Ma after the continental break-up (mature passive margin) and for the intra-continental domain of the Reguibat Shield, we applied a geotherm of 24°C/km (Zarhloule, 2004). Zarhloule (2004) reports present-day values from Moroccan Passive margin of 20 to 35°C/km and these values were applied to their respective areas. Note that as the selected geotherms are based on distinct geodynamic context, they take into account the thermal relaxation that follows the rift-related heating phase. Applying a varying paleo-geothermal gradient, that is also specific for each area, is required to 1) match the changing thermal conditions of the upper crust based on known geodynamic events (e.g. High Atlas rifting; e.g. Gouiza et al., 2017), and 2) define a preferred model to obtain values of exhumation/burial rates and erosion rates with their associated ranges ( fig. 3), as opposed to only ranges using constant geotherms (see following paragraph). To calculate the range of depth-converted results, two constant paleo-geothermal gradients of 20 and 40°C/km were used, together with both a low case (10°C) and high case (30°C) surface temperature. Both gradients are within the range of present-day values of geothermal gradients in Morocco (Zarhloule, 2004), and are considered as "end-member" values. The low and high constant geotherms, applied to the digitized temperature envelopes, expectedly yield the maximum and minimum depths throughout the considered period of time, respectively ( fig. 7). Both sets of low and high constant geotherms and surface temperatures are later used to calculate the range of exhumation/burial rates ( fig. 3). figure  5, assuming geothermal gradients evolving as illustrated in figure 6 and assuming a surface temperature of 20°C. The thicker lines display vertical movement rate calculations; non-selected thinner curves, in the backgrounds, were also converted to depth but were not used in the later study; conditions for result selection are detailed in table 1. The upper and lower limits of the t-Depth envelope (thick dashed lines) are calculated with geothermal gradients of 40 and 20°C/km, respectively.

Defined time periods and burial/exhumation rates
Exhumation and burial rates (km/Myr; fig. 8) were calculated from the depth-converted t-T curves ( fig. 7). Seven periods of time were defined (periods a to g), resulting in up to seven exhumation/burial rate calculations for each curve: Permian (a; 299-252Ma), Triassic (b; 252-201Ma), Early to Middle Jurassic (c; 201-163Ma), Late Jurassic to Early Cretaceous (d; 163-125Ma), middlelate Cretaceous (e; 125-66Ma), Palaeogene (f; 66-23Ma), and Neogene (g; 23-0Ma). The selection of time periods is based on the stratigraphy, and as much as possible on the timing of exhumation and burial events as recorded by t-Depth models. The calculated rates range from -0.09 to 0 (burial) and from 0 to 0.49 (exhumation) km/Myr. The rates calculated from t-Depth models are based on the start and end the above-defined periods. This implies that, for instance, a strong exhumation followed by burial has a similar overall rate of exhumation in any given time period to an overall milder exhumation and vice versa.  figure  7 using the variable geotherms shown in figure 6. The seven defined periods a to g span between 300 and 0 Ma. The combined ranges are extracted from the results of rate calculations using constant geotherms ( fig. 3). North and South are only marked for period g because Morocco has been significantly rotated since the Permian (e.g. Scotese, 2012). Note that the x-axis is not time. Figure 8 illustrates that there were four periods of active and widespread denudation in the study area: Permian, Early-Middle Jurassic, Late Jurassic-Early Cretaceous, and Neogene (periods a, c, d, and g, respectively), in line with results of previous work. During the Permian (period a), sampled basement in the Meseta (e.g. Ghorbal et al., 2008) and the Anti-Atlas (e.g. Oukassou et al., 2013) were exhumed up to 0.12 km/Myr, while the Reguibat Shield (e.g. Leprêtre et al., 2015) showed a slower rate of exhumation of ca. -0.01 km/Myr. During the Triassic (period b), exhumation in the Meseta and the Anti-Atlas slows down (0.01 to 0.05 km/Myr). The High Atlas and most of the Meseta and Reguibat samples were being buried (down to -0.08 km/Myr). In the Early to Middle Jurassic (period c), the basement of the Anti-Atlas (Gouiza et al., 2017) was mildly to strongly exhumed (rate increases up to 0.16 km/Myr). For this region, we observe an acceleration of the exhumation rates from the Triassic to the Jurassic, characterised by the highest rates recorded in this study, with the exception of Neogene exhumation rates. Concomitantly, the regions surrounding the Anti-Atlas were mostly subsiding in the north (down to -0.09 km/Myr) and mildly exhuming in the south (up to 0.06 km/Myr). The Late Jurassic-Early Cretaceous period (d) is marked by stability or limited burial of the Anti-Atlas (down to -0.05 km/Myr), whereas the sampled basement of the Meseta (e.g. Saddiqi et al., 2009), the Reguibat shield (e.g. Leprêtre, 2015), and the High Atlas to some extent (e.g. Ghorbal, 2009), were exhuming (0 to 0.09 km/Myr). During the middle-late Cretaceous (period e), the Meseta, the High Atlas and the Reguibat Shield regions become stable (weak exhumation and burial) with vertical motion rates relative to Earth's surface between ca. 0.03 and -0.03 km/Myr. Exhumation is renewed along the eastern and western rims of the Anti-Atlas during the Palaeogene (period f; 0.01 to 0.03 km/Myr), while other areas remain characterised by low rates, similar to those of the middlelate Cretaceous. Finally, the Neogene period (g) is characterised by widespread exhumation with a clear orogenic signal; which increased from a typical value below 0.10 km/Myr, to 0.20 km/Myr in the High Atlas, and 0.49 km/Myr in the Rif belt, after conversion of Ghorbal (2009) and Romagny et al. (2014) t-T models, respectively.

Uncertainties and limitations of step 1
In addition to the ranges, other uncertainties need to be recognised to fully understand the limitations of the first step of the workflow. Typically, LTT ages in Morocco have ~10% error for both (U-Th)/He and FT systems (Charton, 2018). The t-T modelling shows a substantial temporal and thermal uncertainty associated to all models (e.g. good and acceptable envelopes for HeFTy) and this uncertainty increases with time ( fig. 5). Furthermore, nearly all t-T modelling studies reviewed in this work (table 1) consider post-Variscan cooling events as resulting from erosional exhumation. As reviewed in Malusà and Fitzgerald (2019b), cooling does not exclusively mean erosional exhumation, but may be linked to thermal relaxation after a regional heating event or tectonic exhumation (e.g. footwall of a normal fault), depending on the geological context. The interpretation that cooling exclusively records erosional exhumation (assumed by previous workers and in this study) has a direct control on the related t-T path(s) used to obtain exhumation rates. If this assumption is incorrect, results would need to be modified, which would in turn results in different exhumation maps. In this study, the t-T modelling does not consider the erodibility of the overburden. While it has been documented as having a potentially strong impact on the interpretation of erosion histories (Flowers and Ehlers, 2018), because of the lack of published paleo-geology maps, this study could not consider rock erodibility. This is an area for future refinement of the models. Paleo-geothermal gradients and the assumption that the surface temperature remains constant will also have an impact on the error range. In this modelling, assumptions made regarding the thermal structure of the crust for each region, the gradients and surface temperature are invariable. However, surface temperature will vary temporally, locally, and regionally, while the geotherms may be influenced locally by independent parameters (e.g. crustal heat production, volcanism, large salt province; e.g. Mareschal and Jaupart, 2004). Another assumption has been that paleo-geotherms are not dynamically modified as a function of erosion and deposition, although it is known that both processes can impact the upper crust thermal structure. Erosion may lead to an increase of the paleogeotherm while sedimentation may lead to its decrease (Ehlers, 2005;Gouiza et al., 2017). The use of different analogues for geothermal gradients may also have an impact on the workflow. For instance, in the Mediterranean Sea a paleo-geothermal gradient of >80°C/km was constrained based on t-T modelling for the rift zone of the Tethysian rifting (Malusà et al., 2016). Thus, the applied value of 34°C/km during rifting in this work might be seen as low. However, assessing the sensitivity of this study, we estimate that using such high gradient in the syn-rift would not affect greatly the obtained exhumation/burial rates, as considered t-T realisations for the rifting period in the High Atlas show a flat path fairly close to surface temperature ( fig. 5B). However, using gradients higher than the 40°C/km (used in range calculation, in the rift flanks during the Triassic for example, would result in significantly lower exhumation rates.

3.3.
Step 2: From exhumation/burial rates to exhumation maps 3.3.1. Inputs for interpolation of exhumation/burial rates In order to estimate the volume of material that has been removed through time in the exhumed areas, seven 'exhumation maps ' (fig. 3) have been constructed, using the calculated exhumation rates recorded by t-T modelling and regional stratigraphy (figs. 9, 10, and 11). This step qualitatively considers burial rates and the onshore and offshore basins published stratigraphic columns ( fig. 9) to define areas undergoing burial. The dataset, composed of exhumation/burial rates, is characterised by dense-data (basement massifs in Meseta/Anti-Atlas) and sparse-data areas (e.g. Reguibat Shield). Areas of burial were included to add geological meaning to the exhumation maps ( fig. 9). For each preserved sedimentary basin, one or several points were created depending on the extent of the area, and positioned either in the centre of small basins (e.g. Souss Basin) or at the edges of larger basins (e.g. Tindouf Basin). If sediments were deposited during one of the selected periods and are still preserved in a basin, we attribute a negative rate to all the synthetic points of this basin. When sediments are not recorded in a basin, because they were not deposited or not preserved, we attribute to synthetic points a rate of 0 km/Myr. Synthetic points away from the present-day coastline, to account for subsidence in the rifted margin, have been added along the Continental-Ocean Boundary (COB; Miles et al., 2012). Prior to the Jurassic, there was no oceanic crust, and therefore no COB was present. Nevertheless, for the Permian and Triassic interpolations (periods a and b), we consider the same COB location in order to add synthetic rates in the Variscan chain and Central Atlantic rift zone, respectively. For the Permian only, points along the COB line are attributed exhumation rates of 0.1 km/Myr. This is to account for the tectonic collapse and associated denudation of the Variscan chain, which is documented in Morocco as occurring somewhere between the Carboniferous and the Triassic (e.g. Michard et al., 2008). Similar exhumation rates during post-orogenic collapse have been used in other published models and documented between 0.15 and 0.7 km/Myr (e.g. Clift et al., 2004;Mazzoli et al., 2010;Casini et al., 2015). In this study we applied a lower exhumation rate of 0.1 km/Myr, comparable to the highest rate we calculated for the Permian (period a). From the Triassic onwards (periods b to g), we attributed a burial rate to the COB of -0.01 km/Myr (smaller than the colour bar limit for stable domain, defined between -0.005 and 0.005 km/Myr). This is equivalent to adding a number of synthetic points with burial rates in the slope and/or basinal domains. These values are not necessarily realistic, but are acceptable as volumes in the buried domain are not the focus of this study. . The seven selected periods (a to g) are shown on the left. *Late Cretaceous in the Tindouf basin is present in the eastern and western parts, but not in its central part (Hollard et al., 1985).   . 9) and exhumation/burial rates ( fig. 8). Three domains are defined on the exhumation maps: a subsiding domain with rates ≤ -0.011 km/Myr, a stable domain characterised by rates between -0.01 and 0.01, and an exhuming domain with rates ≥ 0.011 km/Myr. Note that the western boundary is the Continent-Ocean Boundary (COB). The insert grades differences in areas, between the Permian and the Present-Day, and is based on the geological history of Morocco .

Interpolation method
For the interpolation of the exhumation/burial rates, an Inverse Distance Weighting (IDW) algorithm was used. This algorithm does not extrapolate the rates above and below input values and functions with an interpolation grid that allows for the addition of active faults. The interpolation grid extends between 0 to -18°W and 20 to 36°N (with a 10x10 km spacing) and ends at the COB in the west. For the Triassic, Palaeogene and Neogene periods (b, f, and g, respectively), the Atlas fault system is defined according to Frizon de Lamotte et al., (2004). Here, we assume that exhumation rates change following a normal distribution between data points ( fig. 3; step 2.2) and the interpolation method (IDW) calculates each cell value using a linearly weighted combination of surrounding input points (unless interrupted by faults). The descriptions of the exhumation maps are presented alongside the source-to-sink maps (section 4).

Erosion rates
Volume calculations of eroded material, hereafter called "erosion rates", are performed using the Surfer software between the interpolated surface and the plane of null denudation rate (0 km/Myr). Volumes are separately calculated for the following areas (based on polygons of the present-day extent of the basement exposure): Meseta, High Atlas, Anti-Atlas, and Reguibat Shield ( fig. 11; table  2). The erosion rates are an estimation of the volume of removed material, per million years, during the seven considered erosion periods (see similar approaches in other areas: Guillocheau et al., 2012;Grimaud et al., 2018; expressed in km 3 /Myr). They are the calculated upper limit for the total detrital sediment flux, which is especially important for carbonate successions possibly undergoing dissolution and in the case of bypassing the shelf and slope to more distal sinks. The minimum and maximum values ( fig. 11) are obtained by running the same interpolation and subsequent volume calculation using exhumation rates calculated with the high and low constant geothermal gradients. The calculated erosion rates for the seven exhumation maps range from c. 3x10 3 to 34x10 3 km 3 /Myr for the Triassic and Neogene, respectively ( fig. 11; table 2). For regions of specific interest, these rates are between ca. 600 and 7,500 km 3 /Myr for the Reguibat Shield, ca. 600 and 5,200 km 3 /Myr for the Anti-Atlas, and ca. 200 and 4,800 km 3 /Myr for the Meseta and High Atlas massifs ( fig. 11).

Uncertainties and limitations of step 2
The non-restored base maps used for the exhumation maps contain an element of uncertainty, as some areas may have changed in shape and/or surface (the High Atlas rift before the Cenozoic shortening for instance). The use of null exhumation/burial rates for synthetic points where no stratigraphy is recorded also introduces a potential error, as the lack of sediments may be linked to either deposition/erosion (i.e. burial rates should have been used during the assumed time of deposition) or non-deposition (i.e. null exhumation rates). A final uncertainty is the use of published stratigraphy columns as reference, which may be incomplete and/or incorrect. Potential errors associated with the interpolation algorithm and volume calculations are difficult to quantify. Such errors include the choice of interpolation algorithm and selected parameters, mathematical approximations, using all types of fault in a similar fashion, and the placement of the synthetic points for the stratigraphic columns. In conclusion, it has not proved possible to estimate an error for the calculated erosion rates in this workflow, only ranges using arbitrary high and low geothermal gradients may be presented at this stage.

3.4.
Step 3: From exhumation maps to Source-and-Sink maps 3.4.1. Source-to-sink maps: principle and inputs The availability of robust sedimentary, stratigraphic, geochronological and provenance studies, and LTT and numerical (e.g. Landscape Evolution Modelling) analyses allows integration to improve source-to-sink models (e.g. Helland-Hansen et al., 2016;Bhattacharya et al., 2016). Source-to-sink studies have limitations, depending on the spatial and temporal resolutions of each data, or simply on the existence and quality of the sedimentary record. Because the nature of data used in this contribution does not constrain or depend on the transport of sediments ('to'), the presented maps are not labelled as 'Source-to-Sink' but instead as Source-and-Sink (SandS) maps. A recent effort at generating paleogeography maps by combining AFT data with more classical datasets has been carried out in southern NW Africa (Ye et al., 2017), at the sub-continental scale. The work resulted in several paleo-reconstructions around the equatorial Atlantic African margin since the early Mesozoic. This work allowed the inclusion of exhuming areas, as based on LTT and t-T results, onto what can also be described as qualitative SandS maps. Overall, Ye et al. (2017) and our contribution are similar in terms of deliverables, and therefore ours maps are, by design, greatly inspired from the former. As the initial control the the SandS maps, the geological map of Morocco at 1:1,000,000 (Hollard et al., 1985) was digitized, or for neighbouring regions, the UNESCO geological map of NW Africa compiled in 1990 (fig 2A). This generally provided surface geology control points, with the exception of the Permian, for which there are few outcrops. Additional outcrop and fossil data was incorporated from studies (table 3). Several "modifications" to the stratigraphy were added to our composite geological map ( fig. 2B), based on new fieldwork or new published data, especially around the Anti-Atlas (see details part 3.2.2) and in the Central High-Atlas. In the latter, the so-called "Couches Rouge" terrestrial redbeds that have historically (and on the 1:1,000000 map) been attributed a Middle and Late Jurassic age (Hollard et al., 1985) are re-dated to the Middle and Late Jurassic and middle Cretaceous (Barremian to Aptian) (Charrière and Haddoumi, 2016). The well database ( fig. 12) is composed of DSDP and IODP well reports, confidential oil exploration well reports and completion logs accessed from the "Office National des Hydrocarbures et des Mines" (ONHYM). Detailed well data from published work (notably Michard et al., 2008), and limited well data such as total depth (TD), formation at TD, or stratigraphy from published studies or company reports were also used. Large scale reviews of the Moroccan geological history have been carried out in several studies (e.g. Ranke et al., 1982;Le Roy et al., 1997;Nemčok et al., 2005;Sibuet et al., 2012;Ye et al., 2017) and the results have been collated and integrated into our SandS maps. Several types of paleoreconstructions are used as the basis for depositional environments (sub-continental scale, e.g. Ye et al., 2017; to sub-basin scale, e.g. Luber, 2017) and also for the interpretation of tectonic regime (e.g. Ait Brahim et al., 2002). Depositional environment reconstructions for the Phanerozoic at the scale of North Africa from Guiraud et al., (2005) were used to broadly constrain all SandS maps. Paleogeography maps of C.R. Scotese (Paleomap project, 2013) are shown in the inserts below each SandS maps, but were not used in the reconstructions. Sediment provenance analysis conducted in Morocco and surroundings are scarce. Five recent works investigated the provenance with detrital zircon U-Pb (Pratt et al., 2015;Marzoli et al., 2017;Domènech et al., 2018;Azdimousa et al., 2019;Perez et al., 2019); one study used traced elements and radiogenic Nd-Sr isotopes (Ali et al., 2014), one using neodymium isotopes and major/trace elements from DSDP wells (Mourlot et al., 2018), and one produced detrital LTT ages (Sehrt, 2014).
Several studies have also documented paleo-current directions in fluvial systems (e.g. Brown, 1980;Baudon et al., 2009, Mader et al., 2017. Provenance and paleo-current data are used to constrain "source-to-sink" arrows on the presented maps. Finally, we use the exhumation maps presented in step 2 to constrain the source domains (areas undergoing denudation), while modifying their contours, based on the control points described above. The SandS maps illustrate the source, transitional and sink domains from the Permian to Neogene, within a simplified structural framework.

Building the Source-and-Sink maps
The SandS maps are defined according to the International Chronostratigraphic Chart (IUGS, 2020; https://stratigraphy.org/icschart/ChronostratChart2020-03.pdf): Permian (300 to 252 Ma; corresponding to period a from step 2; fig. 13 fig. 22). Exhumation maps from step 2 for periods c, d, and e have therefore been used for two SandS maps each. Four types of depositional environments are described in this study: terrestrial, transitional, shallow marine and marine. For most of the map control points, the depositional environments were already interpreted (table 3). In other instances, we interpreted the depositional environment based on lithology and/or fossil data. The "transitional" environment primarily suggests a coastal situation, but may account for areas characterised by fluctuation(s) between shallow marine and terrestrial environments. Source areas with low exhumation rates (between 0 and 0.01km/Myr) may momentarily act as sedimentary bypass zones, as exemplified onshore in the Neogene map ( fig. 21). Conversely, "terrestrial" domains may have contained source areas. The construction of each map was completed by superimposing all the above-mentioned input layers. Source areas and faults were placed first, then data pertaining to the sink domains. Each feature (outcrop, well, fossil…) was then manually associated with a small coloured halo corresponding to the dominant depositional environment. The extent of the source areas was modified at this stage, modification highlighted by the grey hatched lines, and finally the map completed by connecting depositional environments together. A descriptions of each SandS map and observed patterns are presented in Section 4.

Limitations of the source-and-sink maps
Data quality, data density and temporal resolution are highly variable across the area covered by the SandS maps. This leads to variable robustness of the presented maps, which complicates the comparisons from one to another. Stratigraphic age is the primary uncertainty, as several Phanerozoic layers in Morocco and the surrounding areas have undifferentiated ages (e.g. Hollard et al., 1985). These intervals may contain some dated marine sequences or magmatic intrusions, but these are generally local and extrapolation still has to be made to similar but non-constrained facies. The Permian map is characterised by having the most uncertain dominant depositional environments, as it is constrained with very limited data, while the Neogene is a data-rich time interval and the corresponding map is of higher resolution. A variable aspect is the time-window duration that each map covers. The temporal resolution of the Permian and Triassic maps is low, ca. 50 Myr, compared to the short Middle Jurassic interval (ca. 11 Myr). It is also important to note that the Triassic map is mostly composed of Late Triassic data, as Early and Middle Triassic sediments are rarely documented in the stratigraphy of Morocco. The most accurate dating in the Triassic is obtained at the end of the sequence, where CAMP dykes and sills are radiometrically dated. Ye et al., (2017) presented a technique whereby the onlapping relationships between the Mesozoic sediments and their Precambrian/Palaeozoic basement is used as a proxy to reconstruct the original extent of the basin (after Sloss and Scherer, 1975), which in turn provides the areal extent of the exhuming/emerged zones. We do not follow this approach, as we lack the data on the slope values. These angles are often small and hard to measure in the field. Little published data exist on this subject and it is sometimes unknown if the sediments are truncated, onlapping, or simply completely missing. In this work, we assume that the interpolation results (using data from both source and sink domains) are a fair approximation of the extent of the source and depositional areas. Moreover, our SandS maps use the present-day unrestored geography and geology as base maps. In many areas this is acceptable ( fig. 10 insert), however, there may be more substantial implications (and induced uncertainty) in areas such as the Atlas and Rif belts, for which Cenozoic upheaval was mostly caused by N-S shortening (e.g. Michard et al., 2008) associated with major strike-slip motion (Ellero et al., 2020).

Pre-rift: Permian and Triassic
During the Permian ( fig. 13), erosion of the remnants of the Variscan chain (e.g. Lorenz, 1988;Hmich et al., 2006), mostly effected the Meseta and Western Anti-Atlas. Known Permian deposition occurred in the Eastern Meseta, Doukkala and Argana Valley basins. Subsiding domains include the Central and Eastern Reguibat Shield . We estimate the volume of produced sediments to ca. 1.2x10 6 km 3 , of which little is preserved today. Our working hypothesis is that Permian sediments were re-worked during the Triassic. Another possibility is that detrital material may have travelled outside of the study area to the south, using the Reguibat Shield as a transient sink.
From the Permian to the Early Triassic (periods a to b; figs. 13 and 14), the eroded area, which initially included the entire Variscan Orogenic belt, became restricted to the Meseta and the Anti-Atlas. The eastern Reguibat Shield is also affected by erosion in the Triassic. The occurrence of transitional depositional environments is another important change, resulting from Tethysian marine incursions as far as the Tarfaya basin (e.g. Ranke et al., 1982;Scotese, 2012;Leleu et al., 2016). The erosion that occurred in the Anti-Atlas during the Triassic is supported by sedimentary provenance analyses and paleo-currents conducted in the Massif Ancien (Brown, 1980;Baudon et al., 2009;Domènech et al., 2018), that evidence a drainage systems organised perpendicular to the Anti-Atlas trends.  Table  3. Dominant dep. env.: Dominant depositional environment. Source-to-sink*: Simplified source-to-sink systems evidenced with provenance study or paleo-currents. Structural/tectonics: directional result of compression/extension and shortening/stretching in local and regional studies, respectively. Well data: full (white) points mean that sediments of that age were preserved; empty (transparent) points illustrate that sediments were not deposited or not preserved. WAC: Western African Craton. Hatched grey lines highlight modification made to the interpolation results from step 2. Black squares are highlighting sediments from the surficial geology layer if isolated and barely visible at their original scale.

Syn-rift: Triassic and Early Jurassic
During the Triassic (period b; fig. 14), the northern Meseta (estimated volume of eroded material: ca. >40,000 km 3 ), the Anti-Atlas (ca. 30,000 km 3 ), and the Reguibat Shield (ca. 30,000 km 3 ) were being eroded. A large portion of the Meseta shows considerable burial (e.g. Ghorbal et al., 2008), comparable to the Triassic burial documented in the Central Atlantic and High Atlas rift zones (e.g. Gouiza et al., 2010;Moragas et al., 2016, respectively). The Central Atlantic and Atlas rift zones (i.e. present-day rifted margin and High Atlas mountains, respectively) were subsiding, except perhaps for a part of the Massif Ancien, which acted as a major drainage divide as suggested by Domenech et al. (2018;fig. 14). This massif is believed to have formed a positive structural relief sourcing Triassic sediments to the Argana and Oukaimeden valley (paleocurrents; e.g. Mader et al., 2017;Baudon et al., 2012). Red clastics overlying Precambrian basement in the northern central Anti-Atlas are mapped as Triassic or older, where basaltic flows overlying them yielded ages of 205.9±7.9 and 207.8±6.5 Ma (Fiechtner et al., 1992). This implies that, while the core of the Anti-Atlas continued to exhume (c. 0.04 km/Myr), its northern margin was stable or subsiding (including the Souss and Ouarzazate basins). From the Late Triassic to Early Jurassic (periods b to c; figs. 14 and 15), marine transgressions affected the Atlas and Central Atlantic rift (cf. blue arrows on figure 14). The only significant change in source areas is an exhumation event occurring in the western Reguibat Shield (up to 0.06 km/Myr) starting from the Early Jurassic fig. 15). A recent provenance study by Marzoli et al. (2017), using detrital zircon U-Pb ages, shows that sediments intercalated with the CAMP basalts in the western High Atlas have been sourced from the Meseta domain. For at least part of the Early Jurassic, there was no sea connection between the Western and Central High Atlas. This can be seen in the Early Jurassic record, where facies associations show continental and paralic deposits in the western part of the Central High Atlas, and a general deepening towards the east (Krencker et al., 2020). Sediment supply and accumulation over the Massif Ancien was likely large enough to prevent a marine connection (Krencker et al., 2020).

Early Post-rift: Early Jurassic-early Early Cretaceous
The Early and Middle Jurassic (period c; figs. 15 and 16) are marked by enhanced erosion in the Anti-Atlas and Reguibat Shield, and to some extent in the Meseta (see volumes of eroded material in table 3). Period c has the most active erosion rate for the Anti-Atlas and it is likely that the Anti-Atlas formed a topographic swell, as obtained exhumation rates are higher in the central part of the belt (up to 0.16 km/Myr; e.g. Gouiza et al., 2017;fig. 15). Middle Jurassic redbeds are recorded in the onshore basins north and west of the Anti-Atlas (Tarfaya, Agadir-Essaouira, Central High Atlas, Sidi-Ifni Margin, and Souss basins; see table 3 and references therein), while in the basins south and east of the Anti-Atlas, no Jurassic sediments are recognised (e.g. Hollard et al., 1986;Michard et al., 2008). This supports an exhuming Anti-Atlas and Reguibat Shield, linked by an exhuming or stable Tindouf area. Although the Jurassic is dominated by carbonates, there are significant periods of discrete siliciclastic influx. The Early Jurassic is dominated by carbonate sedimentation but has up to 40 % of fine grained siliciclastics (10-20% on average), most abundant in the Central High Atlas and in the Essaouira-Agadir Basin during most of the Early and Middle Jurassic (Duval-Arnould, 2019). Coarse grained clastic sediments are also deposited in the Central and Western High Atlas during the Toarcian. Middle Jurassic rocks are composed of mudstones to coarse clastics (alluvial plain) and mixed carbonatesiliciclastics in the Western and Central High Atlas, respectively (Malaval, 2016;Joussiaume, 2016). These observations point towards the presence of an active source of clastic sediments in the vicinity of these basins during the Early to Middle Jurassic, corroborating the results of t-T studies conducted in the Anti-Atlas (e.g. Ruiz et al., 2011;Oukassou et al., 2013;Gouiza et al., 2017). The Early and Middle Jurassic (period c; figs. 15 and 16) are fairly similar in terms of depositional environments. Exhumation patterns are identical on both SandS maps as they originate from the same exhumation map, presented in figure 10C. Provenance studies carried out on Middle Jurassic sediments show evidence of sedimentary transport from the Meseta to the Middle Atlas (based on detrital zircon geochronology; Pratt et al., 2015), and from the Anti-Atlas to the Essaouira-Agadir Basin (based on paleocurrents; Stets, 1992). Mobilisation of Triassic salt occurred in the High Atlas and Essaouira-Agadir salt basins from the Early until the Middle Jurassic (e.g. Moragas et al., 2018), presumably where the sedimentary overburden was thickest (figs. 15 and 16). Salt mobilisation in the High Atlas basin stopped prior to the Late Jurassic, but continued in the coastal basin until the Late Cretaceous (Moragas et al., 2018). Exhumation/burial rates are negative in the Anti-Atlas (i.e. indicating subsidence). In the western part of the belt an observed positive rate (fig. 8D) could be due to a t-T modelling inconsistency or it could record remnant relief from the previous period. A provenance study of Lower Cretaceous sediments from the north Tarfaya Basin (Ali et al., 2014) showed that they were sourced from the Reguibat Shield only. It suggests that the positive vertical movement rate calculated for the western Anti-Atlas for period d should be discarded, and that the t-T modelling results for that specific sample are inconsistent with the rest of the t-T modelling (Gouiza et al., 2017). We therefore considered no active source area in the western Anti-Atlas during the Late Jurassic in the corresponding SandS map ( fig. 17). From the Middle to Late Jurassic, our results show a shift in location of sediment production, from the Anti-Atlas to the Meseta. The Late Jurassic (period d; fig. 17) shows a major shift in the sediment source areas: the Anti-Atlas is no longer an active source, while the Meseta region was undergoing erosion (up to 0.07 km/Myr; e.g. Ghorbal et al., 2008). A high-resolution clay mineralogy study, carried out in folded Jurassic sediments of the Essaouira-Agadir Basin (Ouajhain et al., 2011), shows a clear shift in either the sediment source lithology or area, between the Middle and Late Jurassic, passing from a chlorite-to an illite-dominated assemblage. It is possible that Middle Jurassic erosion of the Anti-Atlas cut down to the Precambrian, hence reaching metamorphosed Palaeozoic series, a source of chlorite minerals. The Meseta is also a documented source area during the Late Jurassic, as suggested by paleo-currents measured in the western High Atlas (Stets, 1992). The transition between the Jurassic and the Cretaceous (period d; fig. 17 and 18) was fairly quiescent. The coastline shifted towards the north in the Middle Atlas/Rif areas, and towards the west in the Tarfaya basin. The latter change was accompanied by the onset of large Early Cretaceous deltaic systems (Tan-Tan and Boujdour deltas; fig. 18). The entire Reguibat Shield had been an active source of sediments since the Early Jurassic and remained such during the Early Cretaceous. The prograding delta systems suggests an acceleration of exhumation or a change in source lithology in the Reguibat Shield in the earliest Cretaceous. An acceleration is not observed in the t-T modelling for that region ( fig. 10D), therefore it is possible that erosion in the Reguibat Shield reached the granitic basement at the end of the Jurassic, with near complete removal of the metapelitic Palaeozoic cover prior to the Cretaceous. Exhumation of the Meseta/Atlas massifs (rates of up to 0.1 km/Myr; ca. 50,000 km 3 /Myr) from Late Jurassic to Early Cretaceous (period d) was first described in Ghorbal et al., (2008). The preserved onshore basins in the Meseta do not record Upper Jurassic sediments, except in the coastal Doukkala basin ( fig. 18). It suggests an exposed surface larger than that of the presently outcropping basement. In the Late Jurassic and Early Cretaceous, it is not known whether the Tindouf basin was sourcing sediments to the Tarfaya Basin deltas, acting as a transient sink, or simply stable, as the only two wells that were sampled are located at the transition with the Anti-Atlas in the NW of the basin. The t-T models from the shallowest samples of these wells show mild exhumation rates according to our calculations (c. 0.04km/Myr). A sedimentary provenance study conducted in the north Tarfaya Basin showed that the Lower Cretaceous and Upper Cretaceous sediments were sourced from the Reguibat Shield, and from both the Reguibat Shield and the Anti-Atlas, respectively (Ali et al., 2014; fig. 18 and 19). During period d, the Reguibat Shield witnessed substantial erosion (ca. 3,300 km 3 /Myr), and it is interpreted to be the source of the Lower Cretaceous sediments deposited in the Boujdour and Tan-Tan deltas (with potentially a contribution from part of the south Tindouf basin), with a volume of eroded material of about 125,000 km 3 .

Late Post-rift and Atlas orogeny: middle Cretaceous to Neogene
The Cretaceous (periods d and e; figs. 18, 19, and 20) is characterised by another shift in source areas. The Anti-Atlas basement, which was subsiding during the Early Cretaceous, was exhumed again from the middle Cretaceous onwards (Aptian-Turonian; up to 0.04 km/Myr; Gouiza et al., 2017;fig. 18). Meanwhile, the southern massifs of the Meseta underwent burial during the middle Cretaceous after a prolonged exhumation episode (e.g. Ghorbal et al., 2008). Widespread coarse detrital sediments are described from the Early Cretaceous (e.g. Davison, 2005;Frizon de Lamotte et al., 2009, Luber, 2017. The transition from shallow marine to continental depositional environments occurred between the early and the middle Cretaceous in the Essaouira-Agadir basin. The fluvial deposits record a general paleocurrent direction towards the west during the Latest Barremian to Early Aptian (Luber, 2017). In the Late Aptian and Albian times, the area was drowned once again with the establishment of shallow marine conditions (Luber, 2017). In the Tarfaya Basin, Aptian-Albian (middle Cretaceous) fluvial units exposed along coastal outcrops display a mean paleo-current direction towards the northwest (Arantegui, 2018). Sediment provenance analyses suggest that only the Reguibat Shield (including the Mauritanides) was exporting clastic sediments to the north Tarfaya basin during the Early Cretaceous, while the western(?) Anti-Atlas partly supplied clastics to the coastal basin from Late Cretaceous onwards (Ali et al., 2014). Pratt et al., (2015) collected Albian sediments deposited in the Rif basin, and traced the provenance to two sources: the Meseta and another unknown area. Finally, Mourlot et al. (2018) undertook a large-scale provenance study targeting Albian to Maastrichtian sediments from DSDP wells (from north to south offshore Morocco sites 370, 416, 415, and 369 were sampled) along the NW African coast. They used the major and trace element concentrations from these marine sediments, as well as Nd isotopes, to constrain continental crust signatures. Their results show that the sediments from the DSDP wells north of Agadir have signatures from the Meseta, the Anti-Atlas, and the Reguibat Shield during the middle Cretaceous. DSDP Well 369, in the offshore Tarfaya basin, has a continuous record of sediments with a Reguibat Shield signature throughout the Albian and Late Cretaceous (Mourlot et al., 2018). Azdimousa et al. (2019) analysed detrital zircon for U-Pb ages and documented similar potential source areas for Aptian-Albian, Palaeogene (Eocene), and Neogene (early Miocene) sediments deposited in the north-western Rif. These potential source areas, namely the Meseta, the Triassic sediments from the High Atlas, or the northern West African Craton (Azdimousa et al., 2019), would require several local, or regional to sub-continental drainage systems (to export eroded material towards the Rif basin), respectively. In our maps, we have displayed the most local required system, i.e. between the exhuming Meseta massifs and the northern tip of the Rif (figs. 19,21,and 22). Although the source of detrital material deposited in the Rif may indeed be the Reguibat Shield, sediment fairways between the two massifs were probably disconnected by the Anti-Atlas acting as a geographical barrier, deflecting such routing towards the present-day northeast. Cenomanian-Turonian outcrops onlap Palaeozoic basement (western Anti-Atlas; Abioui et al., 2019) and are characterised by having similar facies with a widespread occurrence around the Anti-Atlas ( fig. 19) (Eastern Anti-Atlas; Guimera et al., 2011). Sauropods and crocodilian ichnofacies from the Guir Hamada (Benyoucef et al., 2015) showed that an emerged land was nearby, which was likely the Anti-Atlas (Gouiza et al., 2017;Charton, 2018;Abioui et al., 2019). In the late Early to Late Cretaceous (period e; fig. 20), subsiding domains are dominant in the study area. This period is characterised by a rise in global sea level (Cenomanian-Turonian transgression; e.g. Piqué et al., 2006), which transgressed the interior of Morocco and Algeria (e.g. Upper Cretaceous deposits in the Guir Hamada; e.g. Benyoucef et al., 2015). It appears to have partially submerged the Reguibat shield, the Tindouf basin (except its central part), and the borders of the Anti-Atlas. A Cenomanian-Turonian carbonate platform with low detrital influx prevailed in the Central High Atlas, while on the Atlantic and Tethysian margins, Turonian organic-rich black muds were deposited, in fault bounded basins where relatively deeper marine environment prevailed (Wang, 2018). Erosion was first localized in the centre of the Anti-Atlas (ca. 600 km 3 /Myr), then extended to the eastern and western regions. During the middle Cretaceous, the central Anti-Atlas became an active source area and exhumation of the Meseta slowed down (from up to 0.06 down to 0.02 km/Myr). By the end of the Cretaceous, the entire Anti-Atlas s.s. was sourcing sediments to surrounding basins, and most of the Meseta and High Atlas domains were subsiding and progressively drowned. The Palaeogene and Neogene (periods f and g; fig. 21 and 22) are characterised by the Atlas orogeny, which records high exhumation rates in the High Atlas and Rif (up to 0.20 and 0.49 km/Myr, respectively) and to a lesser extent in the Anti-Atlas and Reguibat Shield (up to 0.05 and 0.10 km/Myr, respectively). We estimate the volume of eroded material from the study area during the Palaeogene and Neogene to ca. 0.3x10 6 and 0.8x10 6 km 3 , respectively. During the Palaeogene (fig. 21), a large portion of the study area was emerged. Epicontinental basins developed around the exhuming massifs of the Meseta and the High Atlas, and shallow marine setting persisted in the Tarfaya Basin. The Neogene (fig. 22) shows only minor differences with the present-day situation, with nearly all of Morocco emergent. Important sediment source areas are the Meseta, the High Atlas, the Anti-Atlas, and the Reguibat Shield. Some shallow marine sinks developed in the North Tarfaya, southern Settat and Gharb basins, and along the Mediterranean coast in the Rif domain.

Implications for Moroccan source-to-sink systems
The investigation of large-scale denudation patterns has highlighted three regions that act as important sediment sources: the Reguibat Shield, the Anti-Atlas, and the High Atlas/Meseta. The Reguibat Shield is marked by burial from the Permian to the Triassic, followed by exhumation from the Jurassic onwards (exhumation rates: 0.01-0.06 km/Myr; cumulative denudation of ca. 970,000 km 3 ). We interpret that the Reguibat shield was the main source of sediments for the Cretaceous deltas, offshore Tarfaya basin. In the Anti-Atlas, basement rocks were deeply buried in the Permian and then exhumed between the Triassic and the Middle Jurassic (0.01-0.16 km/Myr; cumulative denudation of ca. 330,000 km 3 ). Renewed burial during the Late Jurassic/Early Cretaceous was followed by a final exhumation from the Late Cretaceous onwards (0.01-0.05 km/Myr; ca. 130,000 km 3 ). Presently outcropping Variscan age rocks in the High Atlas and Meseta were near the surface during the Permian/Late Triassic, buried until the Middle Jurassic, exhumed in the Late Jurassic/Early Cretaceous (0.01-0.09km/Myr; cumulative denudation of ca. 50,000 km 3 ), buried again during the Late Cretaceous and finally exhumed during the Cenozoic (0.01-0.20 km/Myr; ca. 160,000 km 3 ). The Triassic-Middle Jurassic burial event was synchronous to Atlasic rifting, and temperature-to-depth conversion indicates rapid burial from close to the surface down to 4 km. In turn suggesting that subsidence associated with the Central Atlantic and/or the High Atlas rift zone(s) extended over nearly the entire Meseta. Our results show that sediment source areas in Morocco since the Permian have varied in location and size and were not always active simultaneously ( fig. 23). The present-day extent of the Precambrian and Variscan basement correlates fairly well to the areas that have experienced exhumation since the Permian. A few exceptions exist (e.g. centre Tindouf Basin, Hauts Plateau, etc…; fig. 23), which can be explained by the fact that their sediment cover prevented the basement from being sampled and analysed by LTT or they could represent uncertainty due to the way the exhumation maps were designed ( fig. 10). We conclude that the Anti-Atlas, the High Atlas and the Meseta massifs have all been important sediment sources to the surroundings sedimentary basins through the investigated periods. The results also highlight the Reguibat Shield as a vast and long-lived source areas ( fig. 23), and may therefore be the most important supplier of clastics material to the onshore and offshore Meso-Cenozoic basins in Morocco and its surrounding; either directly as exemplified for the Tarfaya Basin (Ali et al., 2014) or through several sediment reworking cycles. The contrast in exhumation histories between the Anti-Atlas and the Atlas/Meseta regions, characterised by several shifts of detrital sediment sources, can be explained in terms of movements of different fault bounded terrains separated by a long-lived major lithospheric shear zones (Ellero et al., 2020). The SandS maps illustrate the increasing occurrence of terrigenous sediment in the southern regions, with the presence of the Reguibat Shield as a major detrital material source (feeding the Tarfaya-Dakhla and Taoudeni Basins). This interpretation supports the paleo-reconstructions of Ye et al., (2017) that also identified the Reguibat Shield as an important sedimentary source. Important fluvial systems also existed between the massif and the Taoudeni Basin (southernmost basin on figure 1) and the Mauritania-Senegal Basin. These source-to-sink systems likely existed since the Triassic, and were active in the Berriasian (Early Cretaceous), Cenozoic, and to a lesser extent during the Aptian-Albian and Late Cretaceous (Ye et al., 2017).
The coastal Essaouira-Agadir and Tarfaya basins have accumulated a near complete succession of deposits since the Permo-Triassic, as well as the present-day offshore passive margin domain (Tari and Jabour, 2013;fig. 9). Although major erosional unconformities are documented in the Doukkala, Essaouira-Agadir and Tarfaya Basins as well as offshore in the continental margin (Le Roy et al., 1997;Abou Ali et al., 2005;Gouiza, 2011;El Jorfi et al., 2015), the relatively continuous subsidence in these regions suggests that the widespread erosional exhumations were more focused around the documented source areas.  Hollard et al., 1985). * The transparency for each layer depends on the time duration of related period (e.g. transparency for Palaeocene source areas is 43%).

Mass balance: onshore sources to offshore sinks
The total estimated eroded volume rates (table 2) can be compared to published sedimentation rates from Helm (2009; fig. 24), as exemplified in Tinker et al. (2008) in South Africa. The sedimentation rates in the offshore domain and the coastal basins were obtained using interpreted stratigraphic thickness integrated from nine interpreted seismic profiles perpendicular to the coast (Helm, 2009), and extended by extrapolation and/or well control to the basin (using DSDP wells). Detrital sedimentation rates from the Triassic to the Neogene were then estimated and adjusted to the proportion of clastic sediments recorded in the control wells (Helm, 2009).
Figure 24 | Comparison of the total erosion rates to sedimentation rates in Moroccan offshore and coastal basins (after Helm, 2009). Note that it may be expected that the results should not match perfectly, as a proportion of the eroded material could be removed as dissolved carbonate, reworked after deposition, exported towards Algeria, Mauritania and Mali, or as fines detrital material that could have been transported over long distance beyond the study area of Helm (2009) and ours.
Results from the erosion rate modelling are interpreted to approximate the detrital sedimentation flux (i.e. disregarding routing of sediment in onshore basins, climate, lithology of the source area or suspended load, and the detrital component exiting the system to the distal offshore). In this case study, the comparison shows that the siliciclastic volumes deposited in the offshore/coastal basins and the eroded volume of material from the interpolation grid are, for the most part, of the same order of magnitude ( fig. 24). Deposition and erosion rates are very similar, with relative increases/and decreases matching during the Triassic and Late Cretaceous. The exception is during the Jurassic period, when, according to our results, around 12,000 km 3 /Myr of material were eroded, while the offshore seems to record much lower volumes, about 1000 to 2000 km 3 /Myr of detrital sediments (after Helm, 2009). Such discrepancies, with influx higher than recorded sedimentation rates, could be explained by numerous reasons: a) erosion in the Reguibat Shield, Anti-Atlas, Tindouf, and Meseta may have removed finegrained and/or carbonate-rich meta-sediments from the upper part of the Palaeozoic section. Hence, little coarse detrital material would be recorded in the Atlantic basins, were platform limestones are dominant. b) fluvial system pathways could have re-routed sediments to the east, towards the Meso-Cenozoic Algerian and/or Tunisian basins for instance c) the geothermal gradient could have been higher during rifting, which would result in much lower exhumation rates; d) the offshore study of Helm (2009) did not account for compaction, meaning that sedimentation rates plotted against the erosion rates are minimum estimates. e) the limitation of only using nine sections (seismic lines) as reference along ~2000 km of coast means that depocenters may have been missed or overestimated; f) the algorithm extrapolated outside of data point coverage (e.g. centre of Tindouf Basin); g) longshore and/or contourite systems may have transported the temporarily stored material on the Moroccan shelf beyond the study area, or bypassed to the deeper basins. h) the volume of sediments that has been eroded from the shelf/slope, as recognised in offshore unconformities (e.g. Tari and Jabour, 2013), cannot be quantified (e.g. Tinker et al., 2008). For times when the sedimentation rates are significantly higher than the erosion rates (Early Cretaceous and Palaeogene), it is possible that: a) a substantial sedimentary source may be unidentified if it lacks an LTT/t-T study (see previous part and fig. 23); b) selected t-T models are underestimating the cooling in the related time periods in one or several know sedimentary sources; c) the volume of detrital sediments calculated using offshore wells (Helm, 2009) may not be representative when upscaled to the entire Atlantic margin.

Exhumation/burial wavelengths
The half wavelengths of the exhumation/burial events can be estimated from the evolution of subsidence and exhumation observed along a cross-section trending perpendicular to the strike of the Jebilets Massif, High Atlas, Anti-Atlas and Reguibat Shield ( fig. 25). The half wavelengths, expressed as the distance between successive null points, are between ca. 50 and 600 km for the exhuming domains, and ca. 40 and 200 km for the subsiding domains. The Tindouf Basin remains fairly unconstrained in terms of vertical movements and therefore these wavelengths should be treated with caution, especially before the Late Cretaceous. The longest half wavelength is observed during the syn-rift stage for the Reguibat Shield to Anti-Atlas area. The half wavelengths shorten towards the north, where subsidence prevails in all illustrated periods. The exhumation/burial events may have driven regression or transgression, with positive feedback with eustatic sea level changes. Proposed mechanisms for the vertical movements are large-scale processes (see Teixell et al., 2009 for a comprehensive review) that act at wavelengths from one to several hundreds of kilometres (e.g. Şengör, 2003;Babault et al., 2008;Frizon de Lamotte et al., 2009). They include rift flank uplift (evidence does not support this in the study area), mantle driven doming, lithospheric flexure, crustal-scale folding and erosional unloading. For the subsidence episodes, sedimentary loading was likely enhanced by crustal thinning (rift zone), thermal cooling (old rift), lithospheric flexure and crustal-scale folding (Teixell et al., 2009).

Conclusions
We present a methodology to use t-T modelling results as a proxy to define and quantify exhumation events. We display a series of exhumation maps from the Permian to the present-day in eastern Morocco and its surroundings, from which erosion patterns and erosion rates have been extracted. This allows analysis of the possible mechanism(s) responsible for the observed vertical movements. The presented findings have implications for the evolution of the Central Atlantic rifted margins and for our understanding of the Meso-Cenozoic Moroccan source-to-sink systems. The investigation of denudation patterns at a large scale, has highlighted three sediment source regions for the NW Africa basins: the Reguibat, the Anti-Atlas and the High Atlas/Meseta. During the Permian, terrestrial basins located across the Meseta were sourced with material eroded from the Variscan chain. In the Middle to Late Triassic, widespread rifting allowed more extensive deposition across Morocco. Active sedimentary source areas were the northern Meseta, the western Anti-Atlas and the central Reguibat Shield. Initially, most of Morocco was dominated by continental deposits, with a gradual transgression inland from the Tethys to the east and the proto-Atlantic to the west extending terrestrial/transitional marine environments, across the High Atlas, the Meseta and the Tarfaya basins as well as part of the Reguibat Shield. Throughout the Jurassic, shallow marine and frankly marine environments dominate, with periods of discrete siliciclastic input. Active sedimentary source areas were the Anti-Atlas, the Reguibat Shield, and the Meseta massifs. A substantial shift in detrital source areas was evidenced from the Anti-Atlas to the Meseta/High Atlas at the transition between the Middle and the Late Jurassic (ca. 163 Ma). In the Early Cretaceous, terrestrial environments cover a substantial portion of the study area, especially between 145 and 125 Ma. Another considerable source shift from the Meseta/High Atlas to the Anti-Atlas was evidenced at the transition between the early Early and the middle Cretaceous (ca. 125 Ma). Finally, during the Cenozoic, almost all the presently outcropping basement areas provided sediments to the coastal and foreland basins. The presented workflow enables past large-scale source-to-sink domains to be constrained, which is lacking by design in local LTT studies (e.g. Charton et al., 2018).

Supplementary Files
1. spreadsheet: LTT studies from Morocco, Mauritania, Algeria, and the Canary Islands 2. zip: files of workflow steps 1 and 2

Data availability statement
The well data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to confidential restrictions. All other data is available as supplementary file 2.