The Lithospheric Structure of the Saharan Metacraton From 3‐D Integrated Geophysical‐Petrological Modeling

We modeled crustal and lithospheric thickness variation as well as the variations in temperature, composition, S wave seismic velocity, and density of the lithosphere beneath the Saharan Metacraton (SMC) applying an interdisciplinary 3‐D modeling. Regardless of the limited data set, we aimed at consistent imaging of the SMC lithospheric structure by combining independent data sets to better understand the evolution of the metacraton. We considered that the SMC was once an intact Archean‐Paleoproterozoic craton but was metacratonized during the Neoproterozoic due to partial loss of its subcontinental lithospheric mantle (SCLM) during collisional processes along its margin. This has permitted the preservation of three cratonic remnants (Murzuq, Al‐Kufrah, and Chad) within the metacraton. These cratonic remnants are overlain by Paleozoic‐Mesozoic sedimentary basins (Murzuq, Al‐Kufrah, and Chad), which are separated by topographic swells associated with the Hoggar Swell, Tibesti Massif, and Darfur Dome Cenozoic volcanism. The three cratonic remnants are underlain by a relatively thicker lithosphere compared to the surrounding SMC, with the thickest located beneath Al‐Kufrah reaching 200 km. Also, the SCLM beneath Al‐Kufrah cratonic remnant is significantly colder and denser. Modeling of the lithosphere beneath the Chad and Murzuq Basins yielded a complex density and temperature distribution pattern, with lower values than beneath the Tibesti Massif. Further, our modeling indicated a uniform and moderately depleted mantle composition beneath the SMC. The presence of a relatively thinner lithosphere beneath the noncratonic regions of the SMC is attributed with several tectonic events, including partial SCLM delamination during the Neoproterozoic, Mesozoic‐Cenozoic rifting, and Cenozoic volcanism.


Introduction
Cratons are parts of Earth's continental lithosphere that have not been affected by major tectonic or magmatic events for a long time, possibly since the time of their formation (Black & Liégeois, 1993;Griffin et al., 2009). The stability of the cratons has been largely attributed to the presence of thick, cold, depleted, and anhydrous subcontinental lithospheric mantle (SCLM) that is isolated from the surrounding convecting mantle. However, cratons can lose their long-lasting stability through geodynamic processes leading to either their total destruction, referred to by Yang et al. (2008) as decratonization, or their partial destruction, referred to by Abdelsalam et al. (2002) as metacratonization. There are now several examples of regions that may have experienced decratonization including the North China craton (e.g., Kusky et al., 2007aKusky et al., , 2007bKusky et al., , 2014Santosh, 2010;Zhai & Santosh, 2011;Zhang et al., 2011Zhang et al., , 2014Zhao & Zhai, 2013), the Wyoming craton (Dave & Li, 2016;Föster et al., 2006;Kusky et al., 2014), the Southwest Kaapvaal craton (Eriksson et al., 2009;Kobussen et al., 2008), and the West Amazon craton (Kusky et al., 2014). This process has also been proposed for the Saharan Metacraton (SMC) (Abdelsalam et al., 2002). The term metacraton (Figure 1) is defined as a craton that has been remobilized during tectonic and/or thermal events but can still be recognized through its rheological, geochronological, and isotopic characteristics (Abdelsalam et al., 2002). Later, Liégeois et al. (2013) suggested that the metacratonization processes could occur either at the craton's boundaries or within the craton's interior during lithospheric plate collisions or development of lithospheric-scale shear zones. Abdelsalam et al. (2011) used low-resolution seismic tomography from the model of Grand (2002) to interpret the upper mantle structures of the SMC. The upper 100 km of the metacraton's SCLM has an S wave velocity structure similar to the West African craton, the Congo craton (COC), and the Kalahari craton. However, they also observed that the S wave velocity of the metacraton's SCLM below 100 km depth is much lower than that of the other cratons. Abdelsalam et al. (2011) proposed that such upper mantle structure was due to partial SCLM delamination or convective removal of the cratonic keel, implying that the entire craton  Ganade et al., 2016;Gray et al., 2008;Liégeois et al., 2013). (b) Sketch map showing the cratons, metacraton, and orogenic belts of Africa (modified after Meert & Lieberman, 2008). (c) Tectonic map of the Saharan Metacraton modified from the Tectonic Map of Africa (Milesi et al., 2010). The boundaries of the cratonic remnants are modified from Liégeois et al. (2013). SMC, Saharan Metacraton; MC, Murzuq cratonic remnants; AC, Al-Kufrah cratonic remnants; CC, Chad cratonic remnants; COC, Congo Craton; HS, Hoggar Swell; TM, Tibesti Massif; DD, Darfur Dome; CVL, Cameroon Volcanic Line; K, Keraf-KabusSekerr Suture; J, Jebel Mara; E, Easter Hoggar; U, Uweinat; B, Bayuda; R, Raghane Shear Zone. has been metacratonized. However, Liégeois et al. (2013) suggested that the upper mantle structure of the SMC is heterogeneous and that metacratonization might be unevenly distributed, leaving some remnants of the original cratonic lithosphere unaltered. These include the Murzuq, the Al-Kufrah, and the Chad cratonic remnants (Figure 1).
The understanding of processes affecting cratons is fundamentally important for advancing our understanding of the rates of continental growth and destruction, as well as the current and past volume of the continental lithosphere (Stern & Scholl, 2010). This is certainly an important aspect for the motivation to tackle a relatively new type of integrated, multiparameter modeling. Therefore, in this study, we present a 3-D model, which can show the distribution of lithospheric density and temperature in the SMC. The modeling of petrological and geophysical models is based on the LitMod approach (Afonso et al., 2008;Fullea et al., 2009). The methodology behind the software will be described, in more details, in section 3.2. The interpretation of these model results allows a first step in understanding the geodynamic evolution and the thermal state of Earth's mantle beneath the SMC in order to constrain the present-day thermal and compositional 3-D structure of the lithosphere/uppermost mantle. At the same time, however, this model also shows the inherent problems with which modeling is confronted when the database in the study area is unusually sparse. For investigations in similar remote areas on the Earth with similar problems, we also want to show how a combination of satellite-based measurements and surface-bound data can help to achieve meaningful conclusions. In section 2, the present tectonic provinces in the study area will be briefly discussed and the relevant literature will be referred to. In section 3, the methods by which the data (section 4) are processed and incorporated into the models (section 5) will be explicitly presented. In section 6, the uncertainty of our lithospheric modeling will discussed briefly, and the modeling results and the input data will be mutually compared then interpreted. Some relevant conclusions will be made in section 7.

The Tectonic Settings of the SMC
The SMC is bounded in the east, west, and south by lithospheric-scale suture zones resulting from Neoproterozoic collisional events of the metacraton with the surrounding terranes (Abdelsalam et al., 2002). To the east, a N-trending arc-continental suture separates the SMC from the Arabian-Nubian Shield (see Figure 1a), which represents the northern part of the East African Orogen (Abdelsalam & Dawoud, 1991;Abdelsalam & Stern, 1996;Stern et al., 1994). Similarly, a N-trending shear zone in the west is interpreted as a suture defining the collision zone between the metacraton and the Tuareg Shield (see Figure 1a), which represents the northern part of the Trans-Saharan Orogen (Henry et al., 2009;Liégeois et al., 1994Liégeois et al., , 2003. In the south, as shown in Figure 1a, the SMC is separated from the COC by the E-trending Oubanguides Orogen, which is interpreted as imbricated Neoproterozoic and Archean-Paleoproterozoic thrust sheets, tectonically emplaced southward onto the COC (Pin & Poidevin, 1987;Toteu et al., 2006).
Significant portions of the Precambrian rocks of the SMC are covered either by Phanerozoic sedimentary rocks or unconsolidated sediment, though outcrops of Precambrian rocks can be found in several places in the region, including the Tibesti Massif (TM), the Uweinat (U), the Bayuda (BD), and the Darfur Dome (DD) besides small scattered outcrops (Abdelsalam et al., 2002) (Figure 1c).
Within the SMC, there are three broad sedimentary basins of Paleozoic-Mesozoic age: the Chad, the Murzuq, and the Al-Kufrah. The Al-Kufrah has been identified as an intracontinental sag basin (Heine et al., 2008), and its evolution has been attributed to the lithospheric cooling and thickening of an initially thin lithosphere (Holt et al., 2010). Also, the SMC is surrounded by Mesozoic continental rifts related to the breakup of the Gondwana, and those rifts are grouped at the West and the Central African Rift Systems (see Figure 1c) (Fairhead, 1979;Fairhead et al., 1988).
In addition, widespread Cenozoic volcanism has more recently affected the SMC and is interpreted to be associated with prior continental rifting (Thorpe & Smith, 1974), Neoproterozoic Pan-African crustal reworking (Ashwal & Burke, 1989), or hotspots beneath the HS, the TM, the DD, and the Cameroon Volcanic Line (CVL) (Figure 1c). The volcanic areas of the HS, the TM, and the DD are associated with topographic swells of uplifted Precambrian crystalline basement (see Figures 1b and 1c). It has been suggested that Cenozoic volcanism is fed by the Afar plume (Ebinger & Sleep, 1998), unconnected plumes (Burke, 1996;Wilson & Guiraud, 1992), or by upwelling of the asthenosphere associated with the Africa-Europe collision (Bailey, 1992).

Methods
An important part of our study focuses on modeling the regional-scale thermo-chemical structure of the lithospheric mantle in the Sahara Metacraton. Therefore, a reliable crustal thickness needs to first be defined.

Moho Depth From Gravity Inversion
To develop the Moho depth of the SMC, we follow the methodology of Uieda and Barbosa (2017), who developed a nonlinear inversion algorithm that integrates both gravity and seismic data. The algorithm takes into account Earth's curvature by using tesseroids for forward modeling. An inversion of the gravity field for the calculation of the Moho interface is known to be ambiguous. However, it has been already pointed out, in section 1, that the study area severely lacks surface-based data, including gravimetric data. Nowadays, as a result of modern satellite missions, high-precision gravity data are globally available, practically covering all parts on Earth. These satellite data were used to determine the crustal-mantle boundary in North Africa in the absence of appropriate data. This inversion technique is based on the modified Bott's method (Silva et al., 2014) with a Tikhonov regularization to stabilize the computed solutions in a two-step procedure. As a result of such an inversion, two types of information are available for further processing, namely, the Moho reference depth and the density contrast between the crust and the upper mantle. Since the mean depth of the Moho boundary and its density contrast are poorly known for the SMC, we examined a wide range of combinations. The Moho reference depth ranges from 20 to 40 km with a 2.5 km incremental step, and the density contrast values range between 200 and 500 kg·m −3 with 50 kg·m −3 intervals. In the next step within the inversion process, the Moho depth results are validated against the Moho depth from seismological estimates. Finally, the combination that produces the lowest mean square error (MSE) is selected to give the best fit with the constraining data. For more details on this gravity inversion methodology, the reader is referred to Uieda and Barbosa (2017). The inversion parameters specifically selected for the SMC are described in section 4.2.1.

Integrated Modeling by LitMod3D and Model Assumptions
The LitMod3D (LIThospheric MODeling in a 3-D geometry) was introduced and explained by Afonso et al. (2008) and Fullea et al. (2009). Here we briefly summarize the main aspects relevant to the modeling purposes of this study.
LitMod3D was developed to perform integrated geophysical-petrological forward modeling of the lithosphere and the sublithospheric mantle down to the top of the mantle transition zone (MTZ) at 410 km depth. Any delimitation of model LAB in the forward modeling process is based on both temperature and compositional distributions. Therefore, the lithospheric mantle is defined • thermally, as the portion of the mantle characterized by a conductive geotherm; • compositionally, as the portion of the mantle characterized by a different (normally more depleted) composition with respect to the fertile primary composition in the sub lithosphere, that is, primitive upper mantle composition (PUM) (McDonough & Sun, 1995).
The lithospheric geotherm is computed under the assumption of a steady-state heat transfer condition in the lithospheric mantle, considering a P-T-dependent mantle thermal conductivity (Afonso et al., 2008;Fullea et al., 2009). As boundary conditions, we impose no lateral heat flow across the vertical limits of the model and fix temperature at the surface (15°C) and at the base of the defined lithosphere (1330°C).
Stable mineral assemblages in the mantle are calculated using a Gibbs free energy minimization using Perple_X, as described by Connolly (2005). The composition is defined within the major oxide system CFMAS (CaO-FeO-MgO-Al 2 O 3 -SiO 2 ). The LitMod3D approach allows for a self-consistent calculation of phase equilibria (identity and amount of mineral phases stable at a certain pressure and temperature) and physical properties. Model densities and seismic velocities are determined according to the elastic moduli and density of each end-member mineral as described by Connolly and Kerrick (2002) and Afonso et al. (2008).
After briefly describing the requirements for the modeling and the software, we can now focus on the database in the study area.

Data Sets and Basic Information for Model Building
Data sets and basic information for geometrical and compositional model building of the lithosphere of the SMC can be split up into three groups as summarized in Table 1. 1. Geophysical observables and fields (elevation, geoid height, gravity, and gravity gradients); 2. Crustal thickness information; lithospheric thickness and geometry; and 3. Petrological information on the SCLM compositions.

Geophysical Fields Observables 4.1.1. Elevation
We used the 1 arc-min spatial resolution ETOPO1 global elevation model shown in Figure 2a (Amante & Eakins, 2009). The SMC has a moderate low-lying topography that is interrupted by a number of shorter wavelength swells and uplifted regions, such as the DD (reaching an elevation of ∼1,200 m), the HS (reaching an elevation of ∼1,000 m), and the TM (reaching an elevation of ∼2,000 m). In places, there is relief lower than ∼500 m, and this has been interpreted to be caused by large-scale crustal flexure (Globig et al., 2016;Pérez-Gussinyé et al., 2009) associated with the development of sedimentary basins such as the Chad Basin ( Figure 2a).

Potential Field Data
The geoidal heights were derived from the fifth release of Gravity Observation Combination model (GOCO5S) (Mayer-Güerr, 2015). The GOCO5S model was developed up to a spherical harmonic degree/order (d/o) 280, exploiting the full GOCE (Gravity field and steady-state Ocean Circulation Explorer) mission data with a spatial resolution of about 80 km. In order to remove the contributions of the signal from below the lithosphere, the spherical harmonics up to a d/o 10 were subtracted. The truncation of N<10 is commonly done to eliminate long-wavelength components from the signal, which are  (1995) Note. In this paper, CGS units are used to represent potential field data; these are still frequently used. However, it has been internationally agreed that SI units should be used since 1972. The conversions yield 1 mGal =10 −5 m/s 2 and 1 E =10 −9 s −2 (Eötvös).
associated with sublithospheric sources ). The obtained geoid anomaly, which reflects the distribution of density and the geometry of masses within the lithosphere (Bowin, 2000), is shown in Figure 2b. Despite the low amplitude, the geoid anomalies are highly correlated with topography. Therefore, positive geoidal heights (1-10 m) correlates with the higher altitudes of the DD, the TM, and the HS, while negative geoidal heights (−10 to −1 m) correspond to topographic low values. These low values are associated with the Murzuq Basin overlying the Murzuq cratonic remnant, the Al-Kufrah Basin overlying the Al-Kufrah cratonic remnant, and the Chad Basin overlying the Chad cratonic remnant. The Free-air gravity ( Figure 2c) was derived from the GOCO5S satellite data at 50 km ellipsoidal height and was used to calculate station complete Bouguer anomalies ( Figure 2d) by a mass correction with a density contrast of 2,670 kg·m −3 for onshore and 1,040 kg·m −3 for offshore areas. The height of 50 km was chosen because such height is a good compromise regarding the resolution of lithosphere structures in the gravity field while keeping the noise amplification at an acceptable level by the downward continuation (Mansi et al., 2018;Sebera et al., 2014;Sobh et al., 2019b). All topographic corrections of GOCE gravity data were computed applying the tesseroids software package of Uieda et al. (2016) and the 1 arc-min spatial resolution ETOPO1 model (Amante & Eakins, 2009) (Figure 2).
A ring of negative anomalies of approximately −50 to −80 mGal (1 mGal=1×10 −5 m·s −2 ) ( Figure 2d) is striking, running from the south (COC) first in a northerly direction into the DD area (−95 mGal) and later in a northwestern to western direction via TM to HS (−85 and −80 mGal, respectively). Girdler (1975) recognizes the so-called North African great negative Bouguer anomaly lineament in this pattern. For him and Fairhead (1979) and Fairhead et al. (1988), the existence of the negative Bouguer anomalies is associated with domal uplift of the Cenozoic volcanism and processes in the East African Rift System.
This does not necessarily mean that the North African gravity anomaly is caused by the same processes that cause the anomalies in the East African Rift System through low mantle densities. However, with respect to the three investigated lithosphere structures, the Chad area (CC) with −20 mGal lies completely outside the mentioned ring of negative anomalies. The Murzuq area (MZ) with the basin of the same name is practically divided in the center by the ring of negative Bouguer anomalies (70 mGal), while the Al-Kufrah (AC) region is only touched by the negative anomalies in the very south. These observations will be considered when interpreting the modeling results in section 6. With regard to a synoptic interpretation of the gradients shown in Figure 3, it should be noted that the image of the gradients is inconsistent with regard to the detection of significant density differences in the SMC and the three cratonic remnants. Most likely, the components V zz , which show a Y-shaped pattern of negative gradients, and V yy with rather positive values suggest that the causing density distribution may indicate the existence of anomalous lithospheric density regions.

Basic Information for Model Building
As in many parts of the world with low infrastructure and isolated geographical location, seismic field experiments are rather rare, and thus their interpretation is also rather difficult. Global geophysical studies are most likely to be used to investigate the structure and evolution of these regions. This also applies to the areas where the SMC is postulated. Because of the very satisfying data situation in many parts of the globe after conducting the international gravity missions (CHAMP-GRACE and GOCE), our research also relies on gravimetric data, which can be combined and jointly interpreted with the help of the available seismic constraints, petrological, and geological information. However, in the available models of the lithosphere, large differences are found, which originate from the different spatial resolution of data and/or the modeling technique (Globig et al., 2016;Laske et al., 2013;Reguzzoni et al., 2013;Tugume et al., 2013), spatial constraints, and the evaluation of their quality.

Gravity-Inverted Moho
As already stated in the introduction to this section, we concentrated on the determination of a constrained Moho interface from the GOCO5s satellite gravity model (see sections 3.1 and 4.1.2 for more information) and boundary conditions from seismics. The combined gravity model GOCO05s at 50 km above the ellipsoid (Figure 2) was used. The representation of the field at 50 km is chosen, as it offers a higher level of detail in the signal than at satellite altitude and keeps the noise amplification at an 10.1029/2019JB018747

Journal of Geophysical Research: Solid Earth
acceptable level (Sebera et al., 2014). The inversion methodology (refer to section 3.1) is based on work by Uieda and Barbosa (2017) and has recently been discussed in Sobh et al. (2019a).   (Bouman et al., 2016) at a height of 225 km above ellipsoid. In the single maps, subindices of the gravity potential V denote its second derivatives of the Earth's potential (i.e., the gravity gradients). All gradient components were rotated into an Earth-related coordinate system, which is suitable for forward modeling: X refers to the East, Y to North, and Z points radially to the Earth-related center of the coordinate system in the Earth's center.

Lithospheric Thickness
The LAB as the lowest boundary for the model to be created was taken from the large-scale seismic model of Africa (Fishwick & Bastow, 2011); refer to Figure 5c. A glance at this map clearly shows that at this depth, there are no great similarities between the position of the three remnants (MC-AC-CC) and the depth profile of the LAB. CC seems to be located at relatively shallow LAB, while AC corresponds to a pronounced deeper LAB (with different directional trends), and MC is located at different levels. Additionally, we would like to point out that the published LABs show considerable differences in depths. The difference can be as high as 80 km in areas with sparse seismic coverage (Fairhead & Reeves, 1977;Priestley & McKenzie, 2006;Priestley & Tilmann, 2009). In this study, the LAB depths of the Fishwick and Bastow (2011) compilation has been used as an input for the initial lithospheric modeling. However, the LAB, like the Moho interface, will be substantially modified in the course of modeling.

SCLM Composition
The SCLM composition mainly depends on the tectono-thermal age of the overlying crust (Griffin et al., 2003O'Reilly & Griffin, 2006). Generally, the SCLM is highly heterogeneous in composition, ranging from dunites and harzburgites to lherzolites. Such heterogeneities lead to significant density differences at a given P-T condition. The age of the SMC SCLM ranges from Archean to Neoproterozoic. The oldest rocks in the SMC were found in Uweinat, at the Egyptian-Sudanese-Libyan borders. These rocks are Archean aged (André et al., 1991;Bea et al., 2011;Harris et al., 1984;Klerkx & Deutsch, 1977).
Multiple tectonic events affected, heated, and refertilized the lithospheric mantle during the Pan-African orogeny (Kröner & Stern, 2005), and no significant mantle xenolith compositional data are available for the SMC. The complexity of these processes, the very sparse database, and the inherited limitations of the utilized software forced us to make the modeling more general at a larger scale. This is why we used a representative lithospheric mantle composition of Proterozoic age for our starting model (Table 2).
Various representative lithospheric mantle compositions (Archean and Phanerozoic compositions) have been tested in order to evaluate the effects of these changes on the modeled geophysical observables (see Figure S1 in the supporting information for more details). In our 3-D modeling, the composition is defined within the major oxide system CFMAS (CaO, FeO, MgO, Al 2 O 3 , and SiO 2 ) from Griffin et al. (2009) for the SCLM of the SMC. These compositions range from less fertile, that is, more depleted (high  Table 2 summarizes the SCLM compositions considered in this study.

Setup of Model Geometry and Rock Parameters
Our model covers an area of ∼3,000 × 3,000 km with a vertical resolution of 2 km, down to a depth of 410 km and a lateral resolution of 50 km. Given the lack of knowledge about the characteristics of the lithosphere in The continental crust is divided into an upper and lower crust to realistically model the thermal field by introducing differentiated radiogenic heat production rates and thermal conductivity (Puziewicz et al., 2019) ( Table 3).

Results and Interpretation
In the following sections, the results of the final model are explicitly discussed.

Modeling Results
Using the initial parameters as well as the model geometries described in section 4.2, the results of the initial model did not fit observables (in particular, the gravity field) and is far from the isostatic compensation in most regions (see Figures S2a and S2b). It also shows high residuals in the gravity gradients (see Figure S2c). Positive elevation residuals extending north-south, in correlation with the LAB depths used in the initial model (M 0 ), are clearly seen. That means that the initial Moho and LAB depths are not able to correctly represent the deep lithospheric structure.
In the course of modeling, the lithospheric geometry of both the Moho boundary and LAB depth was iteratively modified minimizing the topographic misfit (residual elevation) with the aim of adapting the long wavelength gravity field and isostatic elevation. The adjustment of the model gradients to the measured ones was only possible because the assumption of isostatic compensation was dropped. The geometry of Moho and LAB in the model was modified manually by matching the modeling results with measured gradients. However, it was not possible to minimize the deviations of model topography and of gravity gradients together. Figure 6 shows that for the maximum gradients, deviations are only ±0.2 Eötvös, but the mismatch in the topography remains large. Sensitivity tests have shown that the gravity gradient data are highly sensitive to the density structure of the crust and uppermost mantle (Bouman et al., 2016;Martinec, 2014); this justifies using the model that fits the gravity gradient data better as our preferred model.

Crustal Thickness
The modeled crustal thickness is displayed in Figure 7a. The crustal thickness map shows a thin crust (∼28 km) along the coastal area and Chad Basin. Along the Al-Kufrah Basin, the crustal thickness reaches ∼38-40 km, increasing to 41 km beneath the Murzuq craton and DD. Our results are generally consistent with the Globig et al. (2016) crustal model with some differences. Inside the three cratonic remnants, the few depths obtained from seismic receiver functions (Figure 7a) coincide well with the gravity-inverted depths, where differences approximately range between −5.0 and 7.5 km with a mean value of 0.52 km and an STD of 2.6 km (see the mismatch histogram).

Lithospheric Thickness
In the forward modeling process, the LAB was modified to get a better fit with measured gravity and gravity gradient fields. At the end, the modeled LAB is deviated from the original LAB compiled by Fishwick and Bastow (2011) with discrepancies up to ±50 km in the area of SMC (Figure 7b). In Figure 8b, the differences Δ between seismic depths (Δ Seismic ) from tomography and modeled depths (Δ Model ) (where Δ=Δ Seismic −Δ Model ) are visualized. It has been recognized that the modeled LAB results are deeper than the initial LAB depths (Figure 5c) of Fishwick and Bastow (2011), particularly in the Al-Kufrah cratonic block. In the South, that is, COC, the modeled LAB is shallower than the initial LAB by approximately 30 km. However, in the cratonic remnants CC and MC, the model results

10.1029/2019JB018747
Journal of Geophysical Research: Solid Earth respectively. These regions are characterized by a relatively faster S wave velocity (Fishwick & Bastow, 2011).

Discussion
Different geophysical methods can be used to detect the LAB; the differences between LAB models arise from the different definitions of the LAB and the different implemented methods for its compilation Artemieva (2009);Eaton et al. (2009). Hence, they are sensitive to other properties. For example, from seismic methods, the LAB can be defined as the change in anisotropy or as the boundary between the high S wave velocity lid and the low velocities in the asthenosphere. Correspondingly, various and different definitions of the LAB exist. In our model, the LAB is described by a thermal (1330°C isotherm), rheological (base of rigid layer), and compositional boundary and is more sensitive to elastic and isostatic changes. Moreover, the thermal and compositional structure of the lithosphere and uppermost mantle of the SMC has been assessed using the integrated geophysical-petrological approach, allowing for the self-consistent combination of geophysical data sets (Afonso et al., 2008;Fullea et al., 2009).
In the SMC, the mantle compositions cannot be constrained by direct data due to the lack of xenolith samples. Therefore, we tested representative mantle compositions for the SMC corresponding to the last tectono-thermal age of overlying crust, as suggested by Griffin et al. (2009). Griffin et al. (2009) and Djomani et al. (2001) have classified the main mantle domains based on their age and the main mineralogical composition. Such a classification seems to be a good proxy for modeling the real world in the studied region, where no mantle samples are available.
We have modeled the SMC lithosphere using two different compositions, namely, the Proterozoic (Proton) and Phanerozoic (Tecton) (see Table 2). In the modeling approach, the mantle densities and velocities are calculated as a function of the temperature, pressure, and bulk composition based on the thermodynamic calculations. The results of the Proton composition are clearly favored over the Tecton composition in terms of the match of the observed potential field data and topography. In the next sections, we will discuss vertical sections of the 3-D model of temperature, density, and velocity structure resulting from our modeling approach.

Density Structure
The modeled lithospheric mantle density at different depths. that is, 40, 75, 100, and 200 km, is presented in Figure 9 alongside four cross sections that cover the key elements of the SMC's lithospheric structure. The distribution of lithosphere densities in the area of SMC and in that of the remnants shows that the lithospheric densities in the Al-Kufrah region differ considerably from the other two (MC and CC) at all studied depths. Lower densities than in the surroundings were modeled in the TM and COC region and perhaps still in HB (Figure 9). In the regions of Chad (CC) and Murzuq cratonic fragments, density differences are always higher compared to their surroundings but relatively modest compared to the densities in the Al-Kufrah area (3340−3350 kg·m −3 ). The cross sections A-A′′ and C-C′′ show that the densities in the Murzuq and Chad area are higher than in the TM area; refer also to cross sections C-C′′ and D-D′′. For the lithospheric mantle, density modeling points to low density values beneath the CVL, TM, and DD due to higher temperatures. Especially in the volcanic areas TM and CVL, very low densities were modeled. In the depth range of 200 km, these disappear again and take on the average densities of the surrounding area in the SMC.
Along the Chad and Murzuq cratonic fragments, the density differences are very small at lower depths. From the cross sections A-A′′ and C-C′′, the density values are high beneath the Murzuq cratonic fragment in comparison to the TM, but their density differences are still very small. The same pattern occurs beneath the Chad cratonic fragment, as shown in cross sections C-C′′ and D-D′′.
Previous studies of Proterozoic SCLM, based on xenolith data (Griffin et al., 2003), have estimated an average density value of 3,340 kg·m −3 at a depth of ∼100 km and a temperature of ∼700°C, corresponding to an upper mantle partly depleted in high-density elements (Fe) and mineral phases garnet (Gt).

Temperature Structure
The subsurface temperatures at various horizontal depth slices and as vertical cross sections are displayed in Figure 10. The modeled temperature distribution does not show a uniform identification of the three

10.1029/2019JB018747
Journal of Geophysical Research: Solid Earth cratonic remnants. Strong lateral temperature heterogeneities are observed, at depths ranging between 40 and 100 km. Also, temperature differences of about 500°C ( Figure 10) and ∼150°C at 200 km depth were modeled. A large-scale low-temperature anomaly is modeled beneath the Al-Kufrah block and the northern edge of the COC, that extends down to ∼200 km. Such results are highly consistent with the seismic tomographic results, showing high-velocity values in the upper mantle (Fishwick, 2010;Schaeffer & Lebedev, 2013). Cross sections C-C′′ and D-D′′ ( Figure 10) show the lateral temperature variation of the lithospheric mantle between the MC and TM and between the CC and CVL, respectively. Along the . The larger Sahara Metacraton is outlined by dashed red lines. Archean shields are exposed at the surface within the Congo region, and the proposed boundaries of the Chad (south), Murzuq (northwest), andAl-Kufrah (northeast) cratonic fragments within the larger Sahara Metacraton are outlined by light gray lines Liégeois et al., 2013;Pasyanos & Nyblade, 2007;Raveloson et al., 2015).
Cenozoic volcanic areas of the DD, TM, CVL, and HS, high temperatures are imaged and characterized by upward deflected isotherms near the LAB (Figure 10).
The high-temperature values reflect the presence of a relatively thin lithosphere above the upwelling mantle beneath the Tibesti and CVL. As a consequence, lower densities ( Figure 9) as well as lower seismic velocities ( Figure 11) can be expected. Recent seismic tomography studies inferred a significant mantle velocity reduction of the S wave velocity in regions of Cenozoic volcanism (Emry et al., 2019;Fishwick, 2010;Schaeffer & Lebedev, 2013), which resolved similar features with an anomalous, low-velocity upper mantle.
Recently, a 10 km thinning of the MTZ was observed beneath the Cenozoic volcanic provinces in central Libya (Lemnifi et al., 2019), which also suggests higher-than-average temperatures at MTZ depths. Along the Al-Kufrah cratonic block, the thickness of the MTZ increases, indicating a colder-than-average MTZ. These results proposed low-velocity upwelling are deeply sourced which correlate with the upper mantle temperature anomaly retrieved from our model. Figure 11a shows a slice of the synthetic seismic velocity computed by our final model (M Final ) at a depth of 125 km. By first approximation, our final model can explain the main velocity structure and reproduces the main velocity structure of global and regional seismic tomography models (both in terms of location and magnitude). Our lithospheric modeling results show high seismic shear velocities that are coincident with the previously proposed Al-Kufrah and Murzuq cratonic blocks (Fezaa et al., 2010;Liégeois et al., 2013) but not with the Chad cratonic block. In our results, low-velocity features were also identified in the shallow upper mantle beneath prominent volcanic areas throughout the SMC, including the HS, TM, and DD ( Figure 11). These low-velocity values coincide with both the lithospheric thinning and the high-temperature values along these regions.

Seismic Velocity Structure
The seismic velocity distribution computed by our preferred model (Figure 11a) shows interesting correlations with the published tomography models. More specifically, the prominent high-velocity values along the Al-Kufrah cratonic block and the northern edge of the COC, which have been resolved by the three recent tomographic models at a global scale, are shown in Figure 11b (Schaeffer & Lebedev, 2013) and at a continental scale, displayed in Figures 11c and 11d (Fishwick, 2010). In Emry et al. (2019), the three proposed cratonic blocks are characterized by high seismic shear velocities at middle-upper mantle depth. These high-velocity features are subdued at shallow upper mantle depths (∼100 km), possibly due to lateral smearing from adjacent low-velocity regions beneath the Hoggar, Tibesti, and Darfur Domes. We conclude that our model cannot properly reproduce the actual velocity structure compared to published tomographic models.

Upper Mantle Thermal Anomaly and Cenozoic Volcanism
The origin of the Cenozoic volcanism of the SMC is debatable, and three contrasting models have been proposed: (1) The association is a manifestation of a deep mantle plume (Burke & Wilson, 1976;Sicilia et al., 2008), possibly linked to the postulated Afar Plume (Ebinger & Sleep, 1998;Sebai et al., 2006); (2) alternatively, the swell and volcanism are the results of the reactivation of lithospheric structures due to the intraplate stress generated by the Africa-Europe collision (Liégeois et al., 2005;Pik et al., 2006;Radivojević et al., 2015); or (3) widespread Cenozoic volcanism affected the SMC, and these are thought to be associated with recent continental rifting (Thorpe & Smith, 1974).
Based on our modeling, we can state with some certainty that the Cenozoic volcanism is correlated with the existence of topographic swells, even with the negative Bouguer anomaly in Figure 2d. The evidence for the distribution of temperature is also obvious, as just shown in Figure 7, only the DD area seems less clear. As a consequence of the modeling results mentioned above, lower densities, lower seismic velocities (less clear) are obtained in areas of the swells. These findings are relatively similar to the seismic modeling results by Emry et al. (2019), Fishwick (2010), and Schaeffer and Lebedev (2013).

Uncertainty in the Lithospheric Modeling
Because the implemented model is too computationally expensive to follow a probabilistic approach, providing a precise quantification of the uncertainties in connection with all individual geophysical properties cannot be made in this study. As expressed previously, it should also be noted that the main characteristic of 10.1029/2019JB018747

Journal of Geophysical Research: Solid Earth
both depths of the Moho and the LAB is the nonuniqueness, which is always the case in gravity-based models. Nevertheless, in this study, the modeling strategy carried out is based on the combination of petrology, mineral physics, and geophysical observables; consequently, different geophysical observables are simultaneously fitted, thus reducing the uncertainties associated with the modeling of these observables alone or in pairs. In Figure S1, we discuss the sensitivity of our models to compositional variations within the lithospheric mantle.
The accuracy of the estimated gravity-based Moho, the lithospheric thickness, and the associated properties of the final model depends on several factors, related to the quality of the input gravity data, the applied corrections, and the geometry of the constraining model from seismic data. A summarized explanation of possible uncertainties related to both data and modeling is made in the following: • Uncertainties associated with satellite gravity data should be acceptably small (Bouman et al., 2016). The gravity model GOCO05s (Mayer-Güerr, 2015), which is synthesized at a height of 50 km, also has only small errors according to Sebera et al. (2014), (approximately 5 mGal). Topographic corrections computed with a constant density of 2,670 kg·m −3 in the Bouguer gravity anomalies (has uncertainties of approximately 4 mGal) may lead to some errors in the estimation of the Moho depths; • Another source of uncertainty arises from the calculation of gravity effect for the elimination of the sediments. Furthermore, a new uncertainty is introduced by the interpolation of the sedimentary map, which is available in low resolution compared to the available high-resolution gravity data (approximately ±3 km in term of Moho depths); • Uncertainties in the Moho depth: In the absence of reliable seismic data, gravity constraints on the crustal structure (such as the Moho depth or the sedimentary thickness) are highly uncertain, given that the gravitational effects of different crustal parameters are similar in amplitude and are also similar in terms of their effects on mantle density anomalies (Herceg et al., 2015). Large parts of our study area lie within the typical range of seismological uncertainty in the Moho depth estimation (approximately ±4 km); • Uncertainties associated with the LAB depth: In our modeling, basic assumptions such as the composition and crustal rock parameters would imply uncertainties (approximately ±20 km) due to the lack of reliable local constraints.

Lithospheric Structure and Metacratonization of the Saharan Region
The interdisciplinary modeling carried out here reveals similarities and differences in the structure of the SMC lithosphere and its characterization by physical rock parameters. A review of the density, temperature, and S wave velocity distributions makes this immediately clear. When assessing the modeling results, it is also important to take into account the fact that there is there is very little and sparse local seismic data in the investigated region and that data sets used here are of global nature (GOCE satellite gravity field and its gradients, S wave tomography, and petrological data of the lithosphere).
Before we try to interpret the modeled parameters density, temperature, and S wave velocity jointly in relation to the initial question, we would like to emphasize another difficulty in the joint modeling/ interpretation of lithospheric densities and seismic velocities. Root et al. (2016) have pointed out in a case study that a simultaneous interpretation of densities and seismic velocities in lithosphere modeling is subject to specific uncertainties. These uncertainties have a great influence on the overall interpretation and are therefore subjected to limitations that are often underestimated, which is also expected to occur in our 3-D modeling. A further look at Figure 11 clearly shows that the three available seismic velocity models result in a diverse maps, which has been accounted for as a starting point for our 3-D modeling.
Moreover, Figure 11 clearly shows such three velocity models reveal an inconsistent picture, because they represent different data and methods. So, the comparison between the calculated velocity from our model and the available seismic tomography model should be qualitative.
In the upper mantle range, density and temperature distributions show a rather variable picture. It is in the nature of the two parameters that the distributions shown are very similar and allow the same conclusions with respect to the three cratonic remnants. Thus, the region of Al-Kufrah (AC) is basically different from the Murzuq (MC) and Chad (CC) remnants. At a depth of 40 km, AC has a relatively high density compared to its surroundings and correspondingly lower temperature values. Also, the remnants (MC) and (CC) show differences compared to the ambient values of density and temperature, but these are comparatively low.
The areas of the North African swells (TM) and (CVL) are well represented in the depth range of 40 km (lower densities and higher temperatures, less pronounced is the area around DD). The result is also supported by the topography (see cross sections in Figures 9e-9h): The topography of the swells could be modeled satisfactorily. The surface extensions of the cratonic remnants coincide with broad depressions around relatively elevated regions (topographic swells: DD, TM, and HS in the corresponding figures) with the location of the Paleozoic-Mesozoic sedimentary basins including the Al-Kufrah, Murzuq, and Chad basins. It seems that during the Pan-African orogeny and the Phanerozoic, the three cratonic blocks played an important role as rheologically rigid bodies.
At the depth range of 100 km (up to 125 km in Figure 11), the picture is not much different: Again, the Al-Kufrah region is striking with comparatively higher densities and lower temperatures compared to its surroundings. With respect to seismic velocities (Figures 10b-10d), the most conspicuous features in the tomographic slices are the high-velocity values along the Al-Kufrah cratonic block, This high-velocity anomaly is clearly visible in our model results (Figure 10a). Apart from these observations, the areas of MC and CC, which show weaker negative deviations in density, remain rather modest in the temperature distribution. The picture in the velocities is not really clear. The models (Figure 10a) for MC show slightly higher velocity in the northern part, and CC is not visible. In the distribution of the model parameters, density and temperature at a depth of 200 km, the SMC shows a tendency toward more uniformity. The presence of magmatism along the cratonic block boundaries, where the Cenozoic volcanic field is present, might have triggered the process that modified the lithospheric mantle beneath the SMC, which, in turn, eventually led to any gravitational instability and destruction of the lower part of the lithosphere (Fezaa et al., 2010;Liégeois et al., 2013).
It remains to be noted that the modeling could only partially confirm the existence of cratonic remnants because important information concerning the SMC is still missing. Thus, the development of the cratonic remnants can only be thought about in a more speculative way. Liégeois et al. (2003) proposed that the metacratonization of the interior of the former Saharan craton was caused by collisional events in the Neoproterozoic along the craton's entire margins. In the east, it collided with the Arabian-Nubian Shield.
To the west, it collided with the Tuareg Shield, which in turn collided with the West African craton. In the south, it collided with the COC. Liégeois et al. (2003) suggested that the Saharan craton collided against an unidentified continent in the North that is buried under the Phanerozoic sediments of North Africa. Hence, the geodynamic setting of the Saharan craton during the Neoproterozoic favors significant shortening of the entire craton. This would certainly lead to instability in the form that will be followed by a total or partial SCLM delamination. It is possible that such delamination occurred along ancient collisional zones that stitched the fragments of the Saharan craton together.

Conclusions
Our 3-D lithospheric model of the SMC combines geophysical, geological, and petrological constraints at various scales into one robust and self-consistent model. The modeling approach calculates the temperature, density, and seismic velocity of the upper mantle as a function of different P-T conditions and petrological composition, aiming to examine the potential existence of three preserved cratonic remnants within the SMC.
Our lithospheric model is able to explain the measured geophysical and geodetic data the general trend of the observed seismic velocities, under the assumption that it is uniformly moderately depleted (i.e., typical Proterozoic mantle). The new inverted Moho depth map provides insight into the crust-mantle interface in a region of sparsely distributed ground data and therefore serves as an example for similar studies all around the world.
Furthermore, our preferred model suggest that the lithospheric structure of these cratonic remnants is different from that of the main SMC in terms of lithospheric thickness. The cratonic remnant of Al-Kufrah is modeled with a relatively thick and, consequently, colder and denser SCLM, while the SCLM beneath the remaining SMC is modeled as relatively thin and, consequently, hotter and less dense.
Based on our developed model, we propose that the SCLM of the cratonic remnants has relatively faster S wave velocities compared to the metacraton, which is in general agreement with the seismic tomography 10.1029/2019JB018747 Journal of Geophysical Research: Solid Earth models, both at continental and global scales. This can be interpreted as localized and partial metacratonization of the ancestral craton of the SMC during the Neoproterozoic collisional-related delamination event.
Later changes in the SMC due to Phanerozoic events such as Mesozoic-Cenozoic rifting and Cenozoic volcanic activity are reflected in the topography, Bouguer anomaly, and temperature distribution.
Still, the mechanisms and the extent of the lithospheric evolution remain ambiguous. Joint inversion of seismological and gravity data in a Bayesian manner is an option to further exploit the solution space (e.g., Afonso et al., 2019;Tork Qashqai et al., 2018). Ideally, such analysis should be combined with geodynamic, numerical modeling to explore different possible scenarios of the evolution of the SMC.

Data Availability Statement
GOCE gravity gradient data can be downloaded from (https://eo-virtual-archive1.esa.int/GOCEGradients. html) (Bouman et al., 2016). The Tesseroids software can be downloaded from (https://tesseroids.readthedocs.io/en/latest) (Uieda et al., 2016). Our modeling results can be found in the supporting information and available at https://doi.org/10.5281/zenodo.3653768. All the used data and the modeling software are cited within the manuscript and listed in the references. This study was carried out in the course of the PhD project of the corresponding author. The project has been supported by the German Academic Exchange Service (DAAD). We would like to thank Folker Pappa for his helpful hints for handling LitMod3D. All grid files and maps were created using Generic Mapping Tools (GMT) version 5 (Wessel et al., 2013). The authors are grateful to the Editor, Prof. Martha Savage, and the two anonymous reviewers for their time, constructive comments, and valuable recommendations on early versions of this manuscript, which really helped us a lot to improve the quality of our article.