Control of inherited accreted lithospheric heterogeneity on the architecture and the low, long-term subsidence rate of intracratonic basins

– Intracratonic basins tend to subside much longer than the timescale predicted by thermal relaxation of the lithosphere. Many hypotheses have been suggested to explain their longevity, yet few have been tested using quantitative thermo-mechanical numerical models, which capture the dynamic of the lithosphere. Lithospheric-scale geodynamic modelling preserving the tectono-stratigraphic architecture of these basins is challenging because they display only few kilometres of subsidence over 1000 of km during time periods exceeding 250 Myr. Here we present simulations that are designed to examine the relative role of thermal anomaly, tectonics and heterogeneity of the lithosphere on the dynamics of intracratonic basins. Our results demonstrate that initial heterogeneity of accretionary continental lithosphere explains long-term subsidence and the arches-basins architecture of Saharan type intracratonic basins at ﬁ rst order. The simulations show that initially heterogeneous lithospheres inherited from accretion are strong enough to resist local isostatic re-equilibration for very long period of time. Indeed, the lateral density variations store potential gravitational energy that is then slowly dissipated by differential erosion and slow vertical movements. For relatively well-accepted coef ﬁ cient of erosion of 10 (cid:1) 6 m 2 /s, the subsidence last longer than 250 Myr. Extensional tectonic forcing and thermal anomalies both result in an effective strength drop of the lithosphere, which allows a temporal acceleration of local isostatic re-equilibration. Periodic changes in far ﬁ eld tectonic forcing from extension to compression complicate the tectono-stratigraphic architecture (intra-basin arches, sub-basins) introducing stratigraphic unconformities between different neighbouring basins such as the ones observed in North Africa.


Introduction
Intracratonic basins also called "cratonic basins", "interior cratonic basins" or "intracontinental sags" host most of freshwater aquifers, minerals resources and hydrocarbon reserves of the world (Allen and Allen, 2013). They have a widespread geographic distribution (see Fig. 6 from Heine et al., 2008) and have in common several features thoroughly reviewed by Allen and Armitage (2011). Here we just summarize them to better define the objects of this study. As their name states, intracratonic basins are located in the interior of continents, far from any active margins (stretched or convergent) upon stable continental lithosphere areas. They are usually large (> 150.00 km 2 in area) circular, elliptical to oval-shaped in plan and saucer-shaped in cross section. Despite their small stretching factors (Armitage and Allen, 2010; Allen and Allen, 2013), they can accumulate large amount of sediments. Their dynamics is characterized by long lasting sedimentation (> 250 Myr) with sublinear to gently exponential shape subsidence curves that corresponds to subsidence rates as low as 5 to 50 m/Myr (Fig. 1). These basins are filled with continental to shallow-water sedimentation indicating low topographic relief. In many cases, their structural framework can be characterized by the association of arches s.l and basins of different kilometric wavelengths reactivated through time (de Brito Neves et al., 1984;Quinlan, 1987;Seyfert, 1987;Perron et al., 2018).
After presenting the main characteristics of accretionary lithosphere and arch and basin structures that typically form on them, focusing on the Saharan Platform, part 3 outlines our working hypothesis before describing and justifying the adopted modelling scheme. Part 4 and 5 details the results of the simulation in term of subsidence curves and basin architecture.
Finally, a ternary classification of intracratonic basins is proposed based on their tectono-stratigraphic architecture. We place the different basins of the Saharan Platform example that inspired this modelling study (Perron et al., 2018;Figs. 2 and 3) into this new classification to discuss what the tectonic architecture of basins tells us about external and internal forcing that are responsible for their long subsidence histories.
The Saharan Platform presents probably one of the rare well-documented example of intracratonic basins in the world where, thanks to recent flexural uplift (Rougier et al., 2013), both the nature of basement and the sedimentary architecture are directly outcropping (Figs. 2 and 3) and can be correlated with seismic and well data. In this area, the terranes have been accreted together during the Eburnean (∼ 2 Ga) and the early Pan-African (∼ 850 Ma) orogeny; (Bertrand and Caby, 1978;Black et al., 1994) long before the Cambrian (540 Ma), which corresponds to the onset of their subsidence (Allen and Armitage, 2011). Geological/geophysical observations of the Hoggar Massif (Bouzid et al., 2008;Brahimi et al., 2018;Perron et al., 2018), where terranes structures are exhumed and outcropps (Fig. 2) indicate the existence of mainly sub-vertical sutures and shear zones of lithospheric scale between terranes of different nature (e.g. Liégeois, 2019). These studies also give an approximative idea of the various dimensions of these terranes (Bouzid et al., 2008;Liégeois et al., 2005;Brahimi et al., 2018).
Focusing on the Mouydir and the Ahnet Basins, Figures 2 and 3 permits to illustrate the concepts of basins and arches through simple syncline-shaped basins and a more tectonised complex-shaped basins. In both cases arches correspond to condense sedimentary series located upon "old" Archean and Paleoproterozoic age terranes while the depocenters rest on The specific zonation of the terranes age with the arches-basins architecture is observed (see also Perron et al., 2018). younger Proterozoic domains (Perron et al., 2018). This typical tectono-stratigraphic architecture in basins and arches (also referenced as paleo-highs) was also described in this area by Eschard et al. (2010).
In both cases, the paleocurrent directions are globally oriented NNW which is more or less parallel to the major lineaments (Beuf et al., 1971). However, local variations of these directions can be observed during Cambro-Ordovician, Caledonian tectonics pulses (e.g. Beuf et al., 1968; Fig. 2A). They can be punctually oriented orthogonally to the main directions. These variations evidence that both a remote sediment supply from upstream and a punctual local erosion/ deposition processes during uplift of arches (i.e. Archean terranes) are active.
During the Paleozoic, the structures of these basins are mainly controlled by N-S sinistral or dextral high dip (> 60°) normal faults (i.e. transtension to transpression) forming horst and graben network associated with forced folds (Fig. 3) weakly inverted and/or reactivated through time (see tectonic history and stresses orientation in Zazoun, 2001;Haddoum et al., 2001;Perron et al., 2018). This structural framework of the sedimentary cover is preferentially nucleated on the basement structures (Perron et al., 2018;Figs. 2 and 3) which are characterized by lithospheric shear zones with predominantly higher dip (Bouzid et al., 2008;Brahimi et al., 2018). Yet, the tectono-thermal intensity in the Ahnet Basin was more significant than elsewhere on the Saharan platform because it was relatively close to the deformation front of the Hercynian orogeny and the West African Craton suture ( Fig. 2A; Boote et al., 1998;Logan and Duddy, 1998;Haddoum et al., 2001;Zazoun, 2001;Coward and Ries, 2003;Akkouche, 2007;Craig et al., 2008). The global syncline-shaped tectono-stratigraphic architecture (1st order pattern) best illustrated in the Mouydir Basin (syncline-shaped basin type; Fig. 3A) can be complexified by the presence of interbasin/intrabasin secondary arch structures (2nd order pattern) observed in the Ahnet Basin (complex-shaped basin type; Fig. 3B). The thickness reaches 1.7 to 7.1 km in the Ahnet Basin which is in average higher than the Mouydir type basins (Beuf et al., 1971;Conrad, 1984;Wendt et al., 2006Wendt et al., , 2009Zieliński, 2012).

Key observations and questions
From the observations listed above, arches and basins structures can be formed above heterogeneous lithosphere separated by vertical mega shear zones. As many modelling studies implies that strength contrast between blocks is often more likely to be reactivated than fault zones (Ranalli, 2000;Buiter and Pfiffner, 2003;Le Pourhiet et al., 2004;Heron et al., 2016;Lafosse et al., 2016;von Tscharner et al., 2016) whether this inheritance recorded in the sedimentation of the basin is related to the weakness of the shear zones or to the rheological contrast between Archean and the Proterozoic lithospheres remains an open question.
Aside from the compositional and mechanical heterogeneities, we cannot a priori disqualify the potential effect of regional thermal events that could affect old accretionary lithosphere, yet, it is legitimate to assume that the geometry predates any destabilizing thermal event related to subsidence.
Nevertheless, while slow exponentially decaying subsidence is often attributed to thermal relaxation, other diffusive processes such as erosion and sedimentation would result in similar signature but with different rates. It is therefore important to quantitatively test the effect of slowly decaying thermal anomalies versus surface processes on subsidence pattern for different thermo-rheological structures.
Finally, the contrast observed between simple synclineshaped basins (e.g. Mouydir Basin) and complex-shaped basins (e.g. Ahnet Basin) seems to be related to far field tectonics solicitation of the system. The sediment routing in the basin highlights that this local tectonic activity is related to change in sources of sediments. Whether these phenomena are first or second order players in the long-term subsidence of basin and arch structures is unknown. While it is always possible to focus on one explanation using oversimplified models, here we have chosen to test the relative importance of all these parameters and their non-linear feedbacks using thermo-mechanical simulations.

General model set up
In order to test the influence of structural inheritance, thermal perturbations, far field tectonics and sediment supply on the subsidence of intracontinental basins, we use thermomechanical simulations. The model domain is 300 km deep and 1600 km long. The mesh is refined towards the surface to enable modelling large-scale dynamics over 250 Myr with a resolution of 500 m at the surface in reasonable computing time. This surface resolution allows us to visualize and constrain the tectono-stratigraphic architecture of the synthetic basins. Appendix A details how conservations of momentum and energy are solved numerically together with surface process evolution.
All the models use a free upper boundary surface subject to erosion-sedimentation allowing the development of sedimentary basins. These surface processes are modelled using Culling's (1965) law with a constant diffusivity (ke) of 1.10 À6 m 2 · s À1 which account for local erosion/deposition processes. All the models have the same thermal boundary conditions. Temperatures are fixed at the top and base of the model to 0°C and 1400°C respectively, and a null heat flux is assumed on the model lateral boundaries (Fig. 4C). The 300 km model domain thickness permits to limit artificially plate growth by conduction to 275 km depth without modelling the whole mantle convection. This boundary condition is similar to the finite plate approximation that is used for oceanic plate geotherms (Parsons and Sclater, 1977).

Structural inheritance
In order to measure the effect of structural inheritance on subsidence, we compare models with laterally homogeneous composition (P and A) with a model M that is largely inspired from the geodynamic setting of the Saharan platform (Perron et al., 2018) where both the age (correlated to the density of the lithosphere) and the geometry of different terranes control the architecture of the basins and arches (Figs. 2 and 3). Geophysical observations in the area (Bouzid et al., 2008;Liégeois et al., 2005;Brahimi et al., 2018), have also contributed to define the dimensions of model M (i.e. size of terranes and shear zones).
This model consists of three 200 km wide Proterozoic terranes separated by two Archean terranes of 100 km in width sandwiched in between two 400 km wide Archean cratons (Fig. 4A). In order to measure the relative effect of weak shear zones versus contrasting rheological profile on inheritance, all the models, A P and M includes weak 2 km wide vertical shear zones affected with a friction of 0.01 and a cohesion of 10 MPa.
The Archean terranes have a 40 km thick crust (20 km quartz dominated upper crust þ20 km anorthite dominated lower crust) and a lighter mantle lithosphere that reflects their high magnesium number. The Proterozoic terranes have a 30 km thick crust (15 km quartz dominated upper crust þ15 km pyroxene dominated lower crust) and their mantles have the same reference density as the asthenosphere. The rheological parameters of the Archean and Proterozoic lithosphere are compiled in Table 1. Notice that these parameters are not constrained by xenoliths data from the Saharan Platform due to their lack. However, they stay coherent with global examples of Archean and Proterozoic lithospheres (e.g. Djomani et al., 2001;Artemieva and Mooney, 2002;Artemieva, 2009). The strength of the lithospheric mantle is limited to 450 MPa (i.e. stress/viscosity limiter) to mimic dislocation glide and Peierls creep (which is not implemented here). The value is chosen according to Watremez et al. (2013) who calibrated the maximum strength of the Arabian mantle using the topography of the Gulf of Aden.

Thermal perturbations
The "cold" lithosphere corresponds to the steady state solution of the heat equation for the boundary conditions imposed on our 300 km thick modelling domain and the 1300°C isotherm often referred as thermal lithosphere asthenosphere boundary is located at 270 km depth. It is similar to the solution of a finite plate cooling model with radiogenic heat production (Turcotte and Schubert, 2014), assuming 300 km plate thickness as the maximum thickness for continental lithosphere and a 1400°C basal temperature.
Nevertheless, although these basins are currently located on stable shields, we cannot disregard that the system may have been out of thermal equilibrium at the onset of subsidence. We therefore design models with homogenously "hot" geotherm that are coherent with the model of Holt et al. (2010) and with a lithosphere having undergone an orogenic cycle (Beuf et al., 1971;Guiraud et al., 2005) and models with localised "thermal anomaly" that are intended to explore the effect of more localised thermal subsidence related to deep thermal events which could correspond to igneous activity observed for example on the Saharan platform, (e.g. Liégeois et al., 1991;Moreau et al., 1994;Derder et al., 2016;Perron et al., 2018). The "hot" and "temperature anomaly" models are designed with non-steady state initial conditions using Burov and Diament (1995) analytical solution for finite plate cooling including radiogenic heat production decreasing from a surface value of 3.10 À9 W kg À1 with a characteristic exponential decay of 10 km. This analytical solution accepts two parameters: 1) a plate thickness and 2) a thermal age which are usually interpreted as due to the last tectono-metamorphic event. The heat production imposed in the sediments is 1.10 À9 W kg À1 , a non-null value motivated by the singular presence of (organic-rich) "hot shales" identified in the North Africa (Lüning et al., 2000(Lüning et al., , 2003. In the "hot" lithosphere models, we use set the 1300°C isotherm at 120 km and use a thermal age of 400 Myr. For the "thermal anomaly" models, the thermal age follows a Gaussian distribution in x-axis (abscissa) from 50 Myr at the centre of the model to 400 Myr on the borders with a standard deviation of 600 km (Fig. 4C). In both cases the initial geotherm is prolongated down to 300 km where temperature reaches 1400°C using a constant linear gradient (∼ 0.5°C/km). During the simulations, these thermally young lithospheres can be destabilized convectively (usually in the initial stages) but more generally, they cool down over time by conduction until they reach a steady state temperature that corresponds to the "cold" lithosphere initial geotherm.

Sediment routing and far field tectonics
In some models, we also add a source term in the Culling model that corresponds to remote sediment supply (equivalent to regional source) by drainage network that runs perpendicular to the cross-section of the model (see Jourdon et al., 2018 for implementation and discussion).
effects of far field orogenic cycle. A sinusoidal variation with time period of 40 Myr peak to peak has been chosen as representative between two shortening events (Fig. 4B). This regular cyclicality of far field tectonic loading is a rough approximation/simplification for the different pulses of extensional and compressional tectonics that have affected the Saharan platform (Ziegler et al., 1995;Coward and Ries, 2003;Perron et al., 2018). The very small velocity vx applied on the boundary, i.e. ± 0.5 mm/Myr (1.5 e À11 m.s À1 ) at peak, ensures a minimal amount of shortening and stretching per cycle (10 km over 1600 km). In order to compensate for stretching and shortening, a small velocity (6vx/16, i. e. ± 0.2 mm/yr at peak) is applied at the base of the model to ensure the volume of the modelling domain remains constant. This value together with the size of the model domain corresponds to background strain rate of 10 À17 s À1 that is considered as rigid plate on the world strain-rate map (Kreemer et al., 2014).

Accretionary vs homogeneous lithosphere
In this part, after briefly describing the results of the simulation with homogeneous Archean (model A) and Proterozoic (model P) lithosphere structured by an imposed tectonic heritage (4.1), we will concentrate on characterizing how the initial geotherm (4.2), far field tectonic forcing (4.3) and interplay between the two (4.4) are recorded in accretionary type of lithosphere (model M). The input parameters for each model (A, P and M) are referenced in Table 2.

Limits of homogeneous lithosphere
The main purpose of the experiments presented in Figure 5 is to show that homogeneous Archean or Proterozoic lithospheres affected by weak vertical shear zones do not produce intracratonic basins when submitted to thermal anomaly (A1 and P1, Fig. 5) nor far field tectonics (A2 and P2 Fig. 5).
On the one hand, thermal anomaly alone (P1 and A1 in Fig. 5) results in a complete lack of deep sedimentary basins after 250 Myr. The subsidence rate displays an exponential decay that is characteristic of thermal subsidence. In both simulations, the subsidence actually ceases after 150 Myr (see W1-P1 and W1-A1 in Fig. 6). On the second hand, the two homogeneous models submitted to far field tectonic loading display small basins that are controlled by the presence of the imposed pre-existing faults. In both cases, the faults, despite their weakness, are not reactivated in a strict sense. Instead, they help initiating new dip slip faults. During the extension phases Archean and Proterozoic terranes have similar thermorheological structure, i.e. crust and mantle are coupled and brittle, and four narrow basins are formed. Yet basins are wider in the Archean than in the Proterozoic terranes because they root slightly deeper according to the rheology of the lithospheric mantle. During compression however, the thick Archean crust displays a decollement at the mid-crust that is not present in the Proterozoic lithosphere (see yield-strength envelope, Fig. 4E). This induces a different mechanical response. The Proterozoic lithosphere deforms preferentially on the inherited weak zones splitting the extensional basin in two sub-basins (P2, Fig. 5) while the Archean upper crust pops up along the dip-slip faults (A2, Fig. 5). In these two runs, the maximum of strain (and of basin thickness) is concentrated along the second shear zone from the limits of the models.
The subsidence curves show a linear decreasing trend with alternation of up and down deviations of amplitude of 110 m for P2 and 400 m for A2 (Fig. 6). The Figure 7 aims to analyse the behaviour of accretionary lithosphere (model M) in response to initial geotherm and tectonic loading. It clearly evidences that all the models displayed capture the first order feature of low rate intracontinental basins unlike the experiments shown in Figure 5. A subtle difference between homogeneous models and heterogeneous ones is the location of active faults when loaded by tectonics. All models with Table 2. List of parameters inputs for each model. Note that duration of all models is 250 Myrs. TA: Thermal anomaly initial geotherm; HL: "Hot" initial geotherm; "Cold" initial geotherm. M5', M7', and M7 models are presented in Supplementary Data.

Models
Lithosphere accretionary lithosphere and without remote sediment supply (Fig. 7) display larger basins on the sides of the model domain than in the centre. In contrast, models with homogeneous lithosphere (A2 and P2 in particular) display larger fault displacement in the centre of the domain. The reason for this contrasting behaviour is the known positive feedback between sedimentation and fault activity (e.g. Beaumont et al., 1992;Burov and Poliakov, 2001;Jourdon et al., 2018). The presence in the model M domain of large stable and buoyant Archean cratonic areas located on either side of the heterogeneity provides the lateral basins with an extra source of sediments that does not exist in model P2 and A2 (i.e. Proterozoic and Archean homogeneous lithospheres with tectonics and only local erosion/deposition). A more obvious difference between the homogeneous lithosphere models displayed in Figure 5 and the heterogeneous ones displayed in Figure 7 is of course the presence of intracratonic basins separated by arches of Archean lithosphere. Consequently, the rest of this contribution will focus on model M that is a mechanically heterogeneous accretionary lithosphere.

Impact of initial geotherm
The first line of Figure 7 represents a "hot" lithosphere which can cool down with time (M2 and M3 in Fig. 7), the second one a thermally equilibrated "cold" (M2' and M3' in Fig. 7), and the last one is a "hot" lithosphere affected by an initial thermal perturbation (M1 and M4 in Fig. 7) that could be related to igneous activity.
The experiment M1 (i.e. heterogeneous lithosphere with a thermal anomaly and local erosion/deposition) displays three bowl-shaped basins created upon the 3 Proterozoic terranes, in our configuration i.e. two peripheral basins and one central basin, separated by basement inter-basin arches upon Archean terranes (Fig. 7). The peripheral basins are 200 km wide and 1.25 km deep. The central basins are thinner (i.e. 0.8 km) and narrower (i.e. 140 km) than peripheral ones.
During the simulation, due to the relief creation, the uplifted Archean terranes get eroded and sediments deposited upon Proterozoic ones. No sediment is deposited on Archean terranes (i.e. arches). The basins form topographic lows, which indicate that sedimentation rate does not compensate for the creation of accommodation space. The basins are starved. The thickness of the different sedimentary layers increases towards the centre of the basins and decreases progressively approaching the arches, forming growth strata. Truncations (M1 in Fig. 7) show that the strata are successively eroded before the next deposits. Consequently, the width of the basins remains stable through time.
For each model with the "cold" lithosphere inputs (model M2' with tectonics and M3' without tectonics in Fig. 7), the major trends of the architecture described previously remain identical. Nevertheless, the thickness of the different basins is smaller than in models with an initially "hot" lithosphere (M1 and M3 in Fig. 7). The maximum of thickness of model M3' (i. e. "cold" heterogeneous lithosphere with only local erosion/ deposition) reaches 1.1 km, which is 300 m less than the model M3 (i.e. "hot" heterogeneous lithosphere with only local erosion/deposition). The analysis of subsidence curves between model M3' and M3 shows more or less a same trend with a shift of 300 m at the end of the 250 Myr (W1-M3' and W1-M3 in Fig. 8). We have shown that thermal anomaly below homogeneous lithosphere (model P and A) cannot explain the long-lived subsidence of intracratonic basins (1st order pattern) and whatever the characteristics of the lithosphere. In less than 150 Myr the equilibrium is reached (see dotted lines W1-A1 and W1-P1 in Fig. 8). On the contrary, with a heterogeneous lithosphere formed by the accretion of different terrains of different densities, the morphology of the curve (see W1-M1 in Fig. 8) indicates that the subsidence, in that case related to buried loads, is still active after 250 Myr. Besides, this part shows that there is no relevant control of initial geotherm on the 1st order subsidence pattern of intracratonic basins (M1, M3 and M3' in Fig. 8).

Impact of far field stresses (tectonics)
The purpose of this second part is to analyse the behaviour of the three types of lithosphere in response to far field tectonic periodic loading by comparing simulation P2, A2 (i.e homogeneous lithospheres; Fig. 5) and M2 (i.e. heterogeneous lithosphere; Fig. 7). The accretionary lithosphere model M2 (i.e. "hot" heterogeneous lithosphere with tectonics) display the formation of arches and basins, which are very similar at first order to the structural pattern obtained with simulation M1 (i.e. heterogeneous lithosphere with thermal anomaly; Fig. 7). At second order, some dissimilarities are however to be noted. The bottom of the saucer-shaped basins of M2 is flat, with angular shape and with some weak undulations that do not happen in the purely thermally-driven model (M4 in Fig. 7). Also, the main basins are formed in peripheral positions rather than in the centre of the heterogeneous corridor. This second part demonstrates that for similar tectonic loading, homogenous lithosphere with faults remain stable, and only allows the formation of small basins at the apex of the shear zones (i.e. the imposed structural heritage). Their subsidence curves display a linear trend with deviations (W1-A2 and W1-P2 in Fig. 6). Heterogeneous lithosphere on the contrary forms basin and arch structures associated to inherited lithospheric heterogeneities in buoyancy and rheology that are very similar to the results obtained with a thermal anomaly. The far field stresses trigger period of acceleration, deceleration and inversion of the subsidence (W1-M2 in Fig. 8) that were identified on the data displayed in Figure 1 (2nd order  signal). It also complicates the architecture of intracratonic basins.

Interplay between tectonic and thermal anomalies
Having shown that basins and arches only form and last for long in case of heterogeneous lithosphere, we now want to evaluate the relative effect of thermal anomaly and far field tectonics on the location and rate of subsidence. This third part therefore aims (1) to dissociate the role of the thermal anomaly from the role of heterogeneous lithospheric column, (2) understand the interplay between tectonic and thermal anomaly using two extra experiments M3 and M4. The M3 model does not include any thermal anomaly nor tectonic forcing, the M4 model at the opposite includes a thermal anomaly and a tectonic forcing. Both models are displayed in Figure 7 together with M1 and M2.
After 250 Myr, M3 simulations show globally the same features as M1 (Fig. 7). The peripheral basins are 200 km wide and 1.5 km deep. The central basins are thinner (i.e. 1 km) and narrower (i.e. 160 km) than peripheral ones. This peculiarity can be explained by the larger surface of erosion of the two Archean blocks at the ends of the models (i.e. the source of sediments is more important) directly feeding the peripheral basins (i.e. sedimentation rate varying according the different basins) and by the thermal doming of the central basin due to the initial thermal anomaly. The starving of sediments in the central basin could be seen as an effect of the 2D model approach. The only difference between these two models is indeed a negative vertical shift of 125 m of subsidence curve of M3 as compared with M1 (W1-M3 and W1-M1 in Fig. 8) that we interpret as initial thermal doming in M1.
After 250 Myr, M2 and M4, the two models subjected to tectonic forcing display more differences than M1 and M3 (Fig. 7). While M2 displays the same overall distribution of depocentres as M1 and M3 and only differs by the flat angular base of the basins, M4 displays more complex distribution of depocentres. Central and boundary sub-basins are indeed separated by inter-basin arches, inter-basin arches secondary arch and intra-basin secondary arches. This specific structural framework has been described in Perron et al. (2018). Moreover, the maximum of deformation is localized in the peripheral basins for M2 and in the central basin, above the initial thermal anomaly, for M4 (Fig. 7).
The secondary arches and basins are controlled by steeply dipping conjugated normal faults (synthetic and antithetic), forming graben structures located from either side of terranes boundaries (and shear zones; M4 in Fig. 7).  Figure 4B for the boundary conditions. These structures are repeatedly activated and re-activated during extensions and inversed during compressions due to far field tectonics (i.e. sinusoidal boundary conditions). The comparison between M2 (i.e. "hot" heterogeneous lithosphere with tectonics) and M4 (i.e. heterogeneous lithosphere with thermal anomaly and tectonics; Fig. 7) demonstrates that the thermal anomaly favours the formation of new faults in the early stage of the simulation. Fault softening allows continuous activity after the thermal anomaly dissipation. The presence of thermal anomaly is recorded by higher strain intensity in the basins that are located above them. Thermal weakening of the lithosphere indeed reduces the thickness of the effective brittle layer and the integrated strength of faults is reduced as their surface diminishes. The subsidence curves are so impacted and show a linear decreasing trend with alternation of up and down deviations of amplitude of 280 m for M2 and 960 m for M4 (W1-M2 and W1-M4 in Fig. 8). The deformation as well as the amplitude of deviations are much more significant above the central thermal anomaly.
The comparison of these four tests (Fig. 7) clearly indicates that slow isostatic rebalancing by differential erosion between different accreted terranes with heterogenic rheological properties (Archean and Proterozoic) can be considered as a driving force for the creation of accommodation. Thermal or tectonic forcing are not necessary conditions for the creation of basins and their preservation through time. However, while thermal forcing alone does not induce very large changes in the distribution and shape of the basins, tectonic forcing is sensitive to the presence of thermal anomalies.
5 Architecture of basins in accretionary lithosphere 5.1 Covering the arches: Impact of remote sediment supply One may note that none of the simulations present sediments covering the arches, as it is the case in the sub-Saharan platform (Fig. 3; Perron et al., 2018). As surface processes permit the local isostatic re-equilibration and controls it (M3 in Fig. 7), we expect that variations in remote sediment supply implemented as a source term might also affect the subsidence of the basins and arches. Thus, we After 250 Myr, we observe the same configuration than the last simulations (M1 and M3 in Fig. 7). The peripheral basins are characterized by a thickness of 3 km and the central basin by a thickness of 2.75 km. Contrary to all previous models, we observe now the presence of sediments on arches up to about 1.5 km thick. The width of the peripheral basins is about 350 km and 300 km for the central one (from the edge to the centre of the arches). The curves display an exponential decay trend and almost reach equilibrium after 250 Myr (Fig. 9C). They show a differential subsidence between peripheral (W1), central basins (W3) and arches (W2). The average rate of subsidence is 12 m/Myr in peripheral basins (W1), 11 m/ Myr in the central basin (W3) and 6 m/Myr on arches (W2).
During the simulation, the Proterozoic terranes and the Archean terranes are differentially subsiding one relative to each other. The addition of remote sediment supply increases the temperature of the basins. The temperatures of the basins in simulations of Figure 7 are < 50°C and > 100°C in Figure 9A (see also Fig. 10). We observe an unexpected rise up of the isotherms under Proterozoic terranes (i.e. basins) and a go down under Archean terranes (i.e. arches). It is caused by the slow burial of the radiogenic heat production of the basement that follows the relative uplift of the Archean terranes regarding the Proterozoic terranes (Fig. 10). This phenomena is discussed by Sandiford and McLaren (2002).
The analysis of these two simulations M3 (i.e. "hot" heterogeneous lithosphere with only local erosion/deposition) and M5 (i.e. same as M3 with remote sediment supply in addition) shows the importance of remote sediment supply rate as a forcing factor on the dynamics and filling of basins (Fig. 10). This parameter also reveals the limits of 2D modelling. Adding more sediments than what is eroded allows to rapidly reach the isostatic compensation by increasing the subsidence rate. It also brings sediment on arches and enlarges the width of the basins. The stratigraphic structures (e.g. onlaps and erosional truncations) upon the arches allow us to understand the complexity of the sedimentary history of the basin.

More complex models
We have now circumscribed the first-order trend controlling the low long-lived subsidence rate and the architecture (arches and basins) of the intracratonic basins. Besides, the second-order trend featured by deviations with periods of acceleration, deceleration and inversion of the low subsidence rate can be explained by far field stresses alternating compressional/extensional pulses (i.e. changes in tectonic processes occurring at the adjacent plate margins). We now want to compare the results at smaller scale by comparing the internal architecture of M5 (Fig. 9), a simple model driven by buoyancy and erosion with remote sediment supply, to M6 (Fig. 11), the same model submitted to both thermal anomaly and tectonic forcing.
In M5, the association of arches and basins (Fig. 9A) is evidenced by divergent onlaps (i.e. growth strata), truncations and reduced thicknesses when approaching the arches (i.e. Archean terranes). The stratigraphic succession is featured by many hiatuses that can be followed on the model at the local scale (one basin) and the regional scale (three basins). The unconformities are particularly well marked on arches where some entire stratigraphic units are missing or are amalgamated. Besides, some stratigraphic units are present in central basin and not in peripheral basins. The differential movement of the basement caused by the relative uplift of the Archean terranes ("lighter") regarding the Proterozoic terranes ("heavier") is at the origin of these observed sedimentary structures (truncations, hiatus and divergent onlaps) near the arches (Figs. 9A and 9B).
Simulation M6 (i.e. heterogeneous lithosphere with thermal anomaly, remote sediment supply and tectonics; Fig. 11) displays the same sub-basins and sub-arches than M4 (i.e. heterogeneous lithosphere with thermal anomaly and tectonics), with more sediments due to remote sediment supply. The left peripheral basin and the central basin are both characterized by a thickness of 4.8 km while the right peripheral one is slightly less thick (4.25 km). We observe the presence of sediments on arches about 2.2 km thick. The maximum of thickness is reached in the central boundary basin which is nearly 5 km thick. The initial thermal perturbation favours the accumulation of more strain on the normal faults that control the central basin. This faster rate of frictional softening is responsible for asymmetrical shape to the central basin (Huismans and Beaumont, 2002) that is not as marked in simulations performed without initial thermal perturbation.  Figure 4B for the boundary conditions. The basins display divergent onlaps (i.e. growth strata), truncations and reduction in thickness when approaching the different arches (i.e. inter-basins or intra-basins on Fig. 11A). The stratigraphic succession features many unconformities. Some entire stratigraphic units are missing in the sub-basins, intra-basin arches and inter-basin secondary boundary arches while present in others (Fig. 11B). The minimum of thickness and the maximum of amalgamated erosional surfaces are detected on the inter-basin principal arches (W5 in Fig. 11). In the central boundary basins (W2 in Fig. 11A), unconformities are observed in the depocentre (maximum thickness recorded in the model) where a continuous conformable stratigraphy would have been expected. It is worth to notice these differences between the basins in terms of truncations and missing series even though the boundary conditions are the same.

Basins evolution: key to deciphering past geodynamics
The far field stresses associated with thermal perturbation parameters bring specificities on the tectono-stratigraphic architecture of the basins (Fig. 11). The arch and basin structural first-order pattern (Fig. 9) is remodelled by the formation of grabens near terrane boundaries during extension, they are inverted during compression. It is defined by subbasins, intra-basin arches and inter-basin secondary boundary arches a characteristic identified in the Saharan intracratonic basins (Perron et al., 2018).
In our case, the lithospheric heterogeneities associated with newly created faults on weakness zones by far field stresses control the compartmentalization and the tectonic kinematics. This individualization of the different structural units and the disparate propagation of the deformation through them explain the diachronic features in the subsidence signal (i.e. acceleration, deceleration and inversion) and the stratigraphic succession architecture between neighbouring basins. For instance, we have highlighted that the layers present in the footwall can be eroded in the hanging wall where the maximum of thickness is usually expected.
First of all, the analysis of 1D well burial history shows that, the initial rate of subsidence in the centre of the basins is greater in models with tectonics (W1-M6, 4 and 6 in Fig. 11C) than in models without far field tectonics (W1-M5 and 3 in Fig.  9C). Nevertheless, a clear tectonic signal is only recorded in the sedimentary architecture of the central basin at the onset of the models with tectonics. The peripheral basins do subside faster, but they do not display large temporal oscillations with 40 Myr cycle before 80 Myr. We infer that this delay reflects the reduced strength of the lithosphere at the apex of the thermal anomaly at the onset of the model. When the thermal signal disappears, after 80 Myr, tectonic deformation tends to distributes itself more equally across the three basins.
This suggests that variations of the sedimentary record of tectonic oscillation in subsidence rate is a good indicator for lateral variations in strength of the lithosphere and that whether these variations are stable or not in time can be interpreted as local thermal (non-stationary) or chemical (stationary) weakening of the lithosphere. In other words, the analysis of the subsidence curves (W1-M6 to W4-M6 in Fig. 11C) shows how the far field stresses and thermal anomaly reduce the strength of the lithosphere. Indeed, models without these forcing parameters have basins with less subsidence (e.g. W1-M5 in Fig. 9C vs W4-M6 in Fig. 11C). showing increase of basin depths and the deposits on the arches in case on higher sediment supply due to remote sediment supply. One may note also that the temperature in the basin and in the crust is affected by this higher sedimentation rate, it is due to heat production of the basement. The isotherms rise up below the basins and down below the arches.  Figure 9 with the set-up of inter-basins boundary secondary arches and intra-basin secondary arches. The strain is concentrated in the central basin with creation of grabens above shear zones (i.e. limits of terranes). (B) Stratigraphy and hiatus repartition between the wells W1-M6 to W5-M6. (C) Subsidence curves of wells W1-M6 to W6-M6 display an exponential decay with deviations of different amplitude depending on their localization (i.e. near maximum strain zones or not). Oscillations related to tectonic loading might be in phase, in antiphase or out of phase between arches and basins or between neighbouring basins but also depending on tectonic context. Comparison with Figure 1 subsidence bibliographic data (dashed lines). See Figure 4B for the boundary conditions. Looking in more details at the sedimentary record of the far field tectonic sinusoidal loading, it is clear that the subsidence curve response is (1) not always sinusoidal and (2) different whether the wells are located on interbasin principal arches (W5-M6, Fig. 11A), intra-basin arches (W3-M6, Fig. 11A First of all, tectonic loading is stronger in the well located in central boundary basins (W2-M6, Fig. 11C) than in any other well. At first order, these basins subside rapidly during extension and uplift in lesser amount during compression. In the details the phases of subsidence last longer than the phases of uplift. Indeed, subsidence starts during the slowing down of far field compression and last to the very end of the extension cycle. It is easier to understand how the system behaves by studying the effect of one tectonic cycle (see M6 movie in Supplementary Materials).
The peripheral basins display very short periods of uplift, which corresponds to (1) maximum subsidence rate in the central boundary basins, (2) a marked increase in subsidence rate in the central basin, (3) maximum uplift of arches and (4) onset of subsidence in the central boundary basins. This short period of time corresponds to the period during which extension rate increases at the boundary. During that time period, the system behaves like a rift bordered by the external normal faults of the central boundary basins and where the arches and the peripheral basins behave like uplifting rift shoulders. As extension decelerates at the boundary, the central boundary basins (W2-M6) continue to subside until extension ceased but the outer part of the system relaxes as shown by the subsiding trend of the arches (W5-M6) and the peripheral basins (W1-M6). The central basin (W4-M6) continues to subside at smaller rate than the central boundary basin (W2-M6), which highlights that the conjugate normal faults are still active.
At the onset of compression, the peripheral basin (W1-M6) subsidence accelerates while the central basin (W4-M6) and the central boundary basin (W2-M6) mark a rapid uplift. This corresponds to a phase of tectonic inversion of the principal boundary faults. At the peak of compression, central basin and the central boundary basins start to subside together with the peripheral basins marking the end of tectonic activity on faults for the tectonic cycle. During that phase, the system in buckling down as a whole.
In summary, the principal faults that bounds the central basins are active through all the extensional phases, regardless the rate, but the activity of the conjugate faults starts only towards the peak of extension rate. During compression, principal boundary faults are only active during the phase of shortening acceleration. After the peak of compression, the system is locked and responds by downward buckling of the whole lithosphere. This asymmetric behaviour between extension and compression phase is well explained by the fact that the lithosphere is stronger in compression than in extension (Brace and Kohlstedt, 1980). These delays in inversion of the fault system versus global buckling may explain why during the single tectonic event, both extensional or compressional structures can be locally identified in the different sub-basins (see Fig. 17 in Perron et al., 2018).

new genetic classification of arch and basin system
While most previous models have considered mainly largescale geodynamic processes occurring in the lithosphere and did not account for the peculiar architecture and intrinsic characteristics of the sediments filling these basins, our study shows that the local differential strength between the terranes permits to build the complexity recorded in the stratigraphic record.
Based on the simulations, we propose a ternary classification of intracontinental basins that relates basin architectures to three forcing parameters: remote sediment supply, far field tectonics and local thermal perturbation (Fig. 13). A best match with geological data (tectono-stratigraphic architecture and subsidence curves) of Saharan intracratonic basins with the results of models can be proposed and plot on that ternary classification (Fig. 13).
We can observe a particularly good fit of the Mouydir Basin (syncline-shaped basin type; Figs. 2 and 3A and see well W21 in Fig. 12) with the geometry of the peripheral or central basins in the numerical model M7 coupling lithosphere heterogeneity, local erosion/deposition processes, tectonics Fig. 12. Comparison of total subsidence curves between numerical models and geological data extracted from the Mouydir Basin, Algeria (9: Well W21; Perron et al., 2018) and the Ahnet Basin, Algeria (5: well W1, Kracha, 2011). The Algerian wells are localized in Figure 2A. The curves are presented in Figure 1. and remote sediment supply (see plot in Fig. 13). Conversely, the data from Ahnet Basin (complex-shaped basin type; Figs. 2 and 3B and see well W1 in Fig. 12) better fits with the model M6 which assumes lithospheric heterogeneities, local erosion/ deposition processes, tectonics, remote sediment supply and thermal anomaly (see plot on Fig. 13). Modeling infers that tectonic forcing can explain the differences in architecture between the two neighboring basins of Ahnet and Mouydir. According to this same procedure of classification, we can plot in Figure 13 other intracratonic basins of the Saharan Platform (e.g. Tamesna and Tim Mersoï Basins).
Upon tectonic forcing, localized vertical movements accommodate basins creation not only by flexure but also by buckling (in compression) or normal faulting (in extension). During extensional events, localized displacements occur along normal faults which do not correspond to vertical weakness zones in the initial geometry. During compressive events, after a short phase of partial inversion of the normal faults, the system locks and weaker ribbons of lithosphere are lifted up by buckling as suggested by many authors (Lambeck, 1983;Cloetingh et al., 1985;Cloetingh, 1986Cloetingh, , 1988Klein and Hsui, 1987;Xie and Heller, 2009). However, buckling instabilities described in these studies are short lived (see also Cloetingh and Burov, 2011), and our proposition to alternate them with period of extension as the advantage of capturing the full complexity of the stratigraphic record. The presence of local and short-lived thermal anomaly in the lithosphere is found to favour faults activity. Here we show that far field stress forcing allows explaining the second order trend characterized by the subsidence deviations pattern and complexification of the structural framework (2nd order pattern). The variations in subsidence rate related to tectonic in our models (Figs. 8 and 11C) are coherent with the data range provided by the literature (Fig. 1).
6.2 Gravitational potential energy: a driver for longterm subsidence?
Several studies have attempted to explain the subsidence at low rate for a long duration (1st order subsidence pattern) of intracratonic basin using classical models (McKenzie, 1978). Alternative models involve cooling following a compressional phase (McKenzie and Priestley, 2016) or mantle delamination (Avigad and Gvirtzman, 2009). The last category involve cooling and thickening of the lithosphere (Holt et al., 2010(Holt et al., , 2015 sometimes accelerated by phase transitions in the mantle lithosphere (Naimark and Ismail-Zadeh, 1995;Baird et al., 1995;Hamdani et al., 1991;Kaus et al., 2005;Eaton and Darbyshire, 2010;Gac et al., 2013). Here, we find that thermal anomalies alone do not produce the first order characteristics of basins and arches framework but that contrast in rheologies and densities of Archean and Proterozoic terranes (all M models) are actually necessary and sufficient to drive slow long-term subsidence observed in intracratonic basins (1st order subsidence pattern). The values of subsidence rate in our models (Fig. 8) are coherent with the range of data provided by the literature (Fig. 1). Fig. 13. Ternary classification of low rate intracratonic basins. The diagram relates typical tectono-stratigraphic basin architecture (flat synclineshaped, syncline-shaped or complex-shaped basin) to internal and external forcing parameters. Example of plotted basins according to their architecture, geometry and subsidence history (see Fig. 12): 1: Ahnet Basin, North Africa (Perron et al., 2018); 2: Mouydir Basin, North Africa (Perron et al., 2018); 3: Tim Mersoï Basin, North Africa (Perron, 2019); 4: Tamesna Basin, North Africa (Lessard, 1961;Perron, 2019). Notice that here the "hot" geotherm ( Fig. 4C) is by default if it is not specified. Basic local erosion/deposition processes are active in each model.
The main driver of long-term subsidence is therefore stored gravitational potential energy inherited from variation in densities between heterogeneous lithospheric columns. Using the initial variation in thickness and density it is possible to determine an isostatic potential for subsidence (Supplementary Data 1). It is interesting to find out that after 250 Myr none of our simulations ever reach the 5 km sedimentary thickness predicted analytically. With large additional remote sedimentary supply, the models reach regional isostatic equilibrium for less than 3 km of sediments in the basins (M5 in Fig. 9) unless tectonic forcing disrupts the strength of the lithosphere (M6 in Fig. 11). This is consistent with the preservation of gravity anomalies related to change in basement lithologies that is observed today in intracratonic area where adjacent columns of Proterozoic and Archean lithospheres are not yet locally isostatically compensated (Gwavava et al., 1996;Perron et al., 2018). These isostatically uncompensated ancient mass excess can also be related to ancient rift zone (DeRito et al., 1983) or dense body in the lower crust (Haxby et al., 1976;van der Pluijm, 1990, 1999) acting as buried loads. The strength of the lithosphere resists to local isostatic readjustment causing downward surface flexure of the lithosphere that is proportional to the buried load. Our simulations show that erosion and sedimentation dissipate this potential energy by redistributing and slowly dissipating the gravitational potential energy stored by these buried loads. It is beyond the scope of this study to discuss the general hypothesis that models should always start at isostatic equilibrium but out studies show that even with very weak vertical shear zones, in initially "hot" lithosphere, the strength of the lithosphere is sufficient to freeze gravitational potential for uplift and release this energy over long time period. The preservation of subsidence over long time scale in our model is the result of local isostatic re-adjustment to erosion/denudation of uplifted reliefs and deposition in depressions (Avouac and Burov, 1996). This process is time dependant because internal loading is modified through time by erosion. As erosion and sedimentation are diffusive processes just like thermal relaxation, the computed subsidence curves can be described approximately by exponential decay. Our model implies that the decay observed in nature is more sensitive to the rate of basin infilling rather than erosion as high sedimentation rates indeed correspond to faster decay (Fig. 9C).

Conclusions
Through a 2D thermo-mechanical modelling, we have applied internal (initial geotherm, far field stresses) and external forcing factors (surface erosion/deposition and constant sedimentation rate) to homogenous lithospheres and accretionary lithospheres (Archean and Proterozoic). From the analysis of the simulations, we can state that: The presence of a thermal anomaly alone is not sufficient to create long-lived basins. Even with erosion sedimentation processes, thermal subsidence is largely reduced after 50 Myr and completely ceased after 150 Myr.
Arches and Basins can emerge from the amplification of the geometry of the terranes through differential erosion sedimentation of Archean/Proterozoic terranes (columns) with different rheologies/densities. The sedimentation rates control the duration of subsidence, typically over 250 Myr in intra-continental context where there are no mountain ranges to provide large sediment supply.
Local sediment supply alone can't cover the arches while remote sediment supply is necessary to both increase the thickness of basin and rise the temperature (supported by the slow burial of radiogenic heat from the basement).
Far field stresses lead to more asymmetric basins and result in the formation intra-basins arches and inter-basin boundary secondary arches delimited by fault-related sub-basins (grabens). They can explain dissimilarities of sedimentary fillings between neighbouring basins as well as the presence of unconformities in the deeper part of the basins.
The effects of tectonics are amplified when a deep-seated thermal anomaly weakens the lithosphere.
The sedimentary piles record complex hiatus, truncation and onlap terminations that differ between the various basins even if the boundary conditions in term of far field stresses are the same. Such a complexity has been often noticed by many geologists when trying to do well correlations in the Saharan Platform (e.g. Perron et al., 2018), for instance.
Erosion/deposition alone cannot dissipate the isostatic potential (gravitational potential energy) because of the strength of the lithosphere (regional isostasy) resists to complete re-equilibration. However, far field stresses or/and thermal events can temporally drop the lithospheric strength and allow subsidence driven by isostatic potential.
Taken as keys to interpret real dataset, we believe that the simulations presented here are simple but realistic enough to constitute a step forward in tectono-stratigraphic trap prediction and heat flux analysis in intracratonic basins.