A New Lithospheric Density and Magnetic Susceptibility Model of Iran, Starting From High‐Resolution Seismic Tomography

We propose a new model for the crust and upper mantle in Iran by joint inversion of gravity and magnetic fields, constrained with a seismic tomography model. We then calculate shear modulus from the Vs velocities and densities. The crust and mantle tomography model is first converted to a density cube through empirical and petrological velocity‐density relations. The starting susceptibility is assigned to a two‐layer homogeneous model, above a heat flow‐derived Curie depth. Considering the uncertainties in the density‐velocity relations, and the starting layered susceptibility variation, we refine the model by a constrained inversion of the gravity and magnetic fields with a Bayesian approach, producing the final 3D density and susceptibility model. The area is tectonically active with high seismicity and active faulting which are regulated by the crustal density and rigidity variations. Higher rigidity matches lower seismicity and extended deserts and basins, suggesting the control of their development. The Neo‐Tethys suture, extending ∼1,500 km long, as well as the Paleo‐Tethys suture match crustal scale density variations, defining characteristics of the lighter Arabian plate and denser Eurasian crust. The South Caspian Basin is enigmatic, due to focusing on the seismicity along all its borders, but with relatively low average rigidity, which is contrary to what is observed for Iran, where the reduced rigidity correlates with higher seismicity. The 3D density model will be useful for numerical geodynamic models and obtaining geologic inferences from the crustal‐scale units.

The extreme structural complexity and the presence of numerous faults make the area seismically very active.In his quantitative analysis, Raeesi et al. (2017) identified 11 faults in the mature stage of the seismic cycle, located in the northern and eastern parts of Iran.They also report four seismically mature regions in southern Iran and four additional fault segments potentially active even without a big sign of maturity.
One open question regards the reasons for the extreme inhomogeneity of the spatial distribution of seismicity observed in the region, which we study in terms of correlation with the density structure of the crust.The density and the seismic velocity are a means to calculate the elastic rigidity, which in turn is the parameter that together with the bulk modulus controls the strain response to tectonic stress.A further question relates to the directions of faults in Iran which take up well-defined homogeneous orientations in tectonic units or between tectonic and geologic blocks, and which must be guided by rheologic inhomogeneities that we find to correlate with the density structure.A third question we address is whether the important sutures that have been defined in Iran can be associated with crustal or lithospheric scale transitions across the sutures, identifying the properties and extents of the original fragments.The sutures demonstrate the amalgamation of lithospheric fragments detached from Gondwana, which were dislocated from the Gondwana margin through the opening and closure of oceanic basins, as the Paleo and Neo-Tethys.The high spatial resolution of the final density model (25 km), achievable particularly for the upper crust compared to the tomographic seismic model, allows to considerably improve understanding of the mentioned problems and acquire new insights.
To exploit gravity and magnetic data, we use a joint Bayesian inversion algorithm.This approach needs, as starting point, a-priori information about density and magnetic susceptibility distributions.To build this starting model, we use known empirical and experimental relations (Brocher, 2005) and techniques for the conversion from velocity to density within the crust, while we exploit petrological properties for the upper mantle.The starting susceptibility is a homogeneous layer above the Curie depth.The inversion delivers a final model of the Iranian lithosphere, with the distribution of densities and susceptibilities, which suggests valid geologic and geodynamic interpretations.

Geological Setting
The collision zone of Arabia-Eurasia in Iran is part of the Alpine-Himalayan orogenic belt which began forming in the late Eocene (Ballato et al., 2011;McQuarrie & van Hinsbergen, 2013) and intensified in the Miocene (McQuarrie & van Hinsbergen, 2013).The factors controlling the kinematics of the area have been diverse in time, like the extension of the Mediterranean Sea and the rifting of the Red Sea and the Gulf of Aden, or the different stages of opening and closure of the Paleo-Tethys and Neo-Tethys oceans (Irandoust et al., 2022a;Zhu et al., 2021).Presently the tectonics is driven by the northward movement of the Arabian plate relative to the Eurasian plate, leading to high rates of seismicity (Raeesi et al., 2017), including destructive events as the Turkey-Syria earthquake of February 2023 (Barbot et al., 2023), or the 2003 Bam earthquake in South East Central Iran (Nadim et al., 2004).Iran has a complex structure, outcome of the amalgamation of blocks of different ages and formed due to different geologic events, including collision of lithospheric plates, and subduction of continental and oceanic lithosphere (Agard et al., 2011;Ballato et al., 2011).The Central Iranian territory (Figure 1) is boarded, from a structural point of view, at south-west by the Zagros orogen, at south by the Makran oceanic subduction zone, and north by the Alborz and Kopeh-Dagh (KD) orogenic belts, and east by the Sistan Suture zone, which marks the transition to the Afghan Block (Irandoust et al., 2022a).The ongoing convergence between the Eurasian and Arabian plates is estimated from GPS at a present rate of ∼22 mm/yr, and is accommodated by northwards shortening within the Zagros, Central Iran, Alborz, and relative movement of the South Caspian Basin (SCB) and movement along strike-slip faults of different orientations, according to the region (Vernant et al., 2004).Deformation of the Iranian orogenic belts is presumably guided by variations of crustal rigidity and the presence of rigid blocks, a question we address in this work.We identify different geological 10.1029/2023JB027383 3 of 25 provinces and blocks in the study area (Gharibnejad et al., 2023;Irandoust et al., 2022a;Rabiee et al., 2019), as shown in Figure 1.We briefly describe the principal properties of these provinces in the next subsections.

Zagros Orogenic Belt
The Zagros orogenic belt separates the Iranian domains on the Eurasian plate, from those on the Arabian plate.It is part of the much larger Alpine-Himalayan orogenic system and extends from the northwestern Iranian border to southeastern Iran, from the East Anatolian strike-slip fault to the Oman Line (Alavi, 2007;Mouthereau, 2011).It is divided into three different parallel tectonic zones which are, from the northernmost, the Urumieh-Dokhtar Magmatic Arc (UDMA), the Sanandaj-Sirjian Zone (SSZ), and the Zagros Fold-and-Thrust Belt (ZFTB).UDMA is the most internal part of the orogen.It is a magmatic system composed of both intrusive and extrusive rocks in correspondence with the back-arc area of the Neo-Tethys oceanic lithosphere subduction.The oldest rocks are Early Cretaceous, while the youngest mainly include recent lava flows and pyroclastic layers (Alavi, 2007).The SSZ forms the boundary of the Eurasian plate.It is mainly composed of both sedimentary and metamorphic rock facies, with rocks that were transported in a period from the Paleozoic to the Cenozoic by the thrust system (Alavi, 2007;Mouthereau, 2011).In the SSZ, we find the suture line due to the closure of Neo-Tethys ocean, crossing the province from north-west to south-east.Last, the ZFTB is the southern and outermost part of the orogen, and it is underlain by the basement of the northeastern margin of the Arabian plate.It is the part that in the zone accommodates a good percentage of the continental plate movements, with the presence of numerous strike-slip and thrust faults along the belt.

Alborz Belt
The Alborz belt (AB) is a tectonically active orogen that extends in northern Iran from the Lesser Caucasus southwards and wraps around the southern SCB.According to different studies (Stöcklin, 1974;Şengör et al., 1988), the belt formed through the Cimmerian orogeny which affected parts of Cimmeria separated from Gondwana and collided with the Eurasian plate during the Triassic.The evidence of this collision is the Paleo-Tethys suture which separates the belt from the SCB (Guest et al., 2006), while the presence of Paleocenic folds and faults are explained by the accommodation of compressional events that occurred during the closure of the Neo-Tethys along the Zagros suture (Guest et al., 2006).In the suture area, the presence of pre-Jurassic metamorphosed ophiolites and deep-sea sediments, deformed by intense thrusting and folding, and partially covered by Cenozoic strata has been defined (Alavi, 1996).The main faults are thrusts at the northern and southern external domains of the Alborz range, indicating double vergence, southwards/northwards dipping at the northern/southern domains, respectively.Inside the range, vertical left lateral strike-slip faults dominate.The strain displacement along the faults is estimated to vary between 2 and 7 mm/yr according to the location and fault mechanism (Ballato et al., 2011;Djamour et al., 2010;Rashidi et al., 2023).

East Central Iran and Sistan Suture Zone
East Central Iran has been defined as a microcontinent (ECIM), comprised of the Anarak-Jandagh Block (AJB), Yazd Block (YB), Tabas Block (TB), and Lut Block (LB) (Gharibnejad et al., 2023).Toward east the ECIM is separated by the Sistan Suture Zone (SSZ) from the Afghan Block or microcontinent, also termed Helmand Block (HB) (Irandoust et al., 2022a) (see Figure 1).The Sistan suture in east Iran formed as an accretionary complex from the closure of the Neo-Tethyan ocean basin, which extended in north-to-south orientation between the LB and the Afghan microcontinent.The suture extends NS for 700 km and forms the border between Iran, Afghanistan, and Pakistan.The suture is post-Cretaceous and includes melanges of ophiolitic rocks emplaced in schistose rocks.Of relevance for the density modeling is the documented presence of high-density rocks such as gabbros, basalts, and slabs of ultramafic and mafic rocks, and eclogite or amphibolite (Bröcker et al., 2022).The geochronology of the metamorphism and post-metamorphic intrusions, and the age of the protolith gabbro are discussed in Bröcker et al. (2022).Documented ophiolites line the western margin of the LB (Saccani et al., 2010).The center of the LB hosts the Dash-e Lut desert (DLD).Strike-slip faulting with highly destructive earthquakes characterizes the eastern and western margins of the LB, generally with strike-slip faulting.The LB is described based on geologic observations of faulting as a rigid block, limited to the east by the Nehbandan fault system which overprints the Sistan suture.At the western margin, the LB is separated from central Iran by NS or NNW-SSE trending right lateral fault systems (Rashidi et al., 2019).The extremely active seismicity struck with seven earthquakes with magnitude M > 6 in the short time period between 1981 and 2011, which includes the destructive M 6.5 Bam earthquake of 2003 (Rashidi et al., 2019).We shall quantitatively ascertain if the rigidity variation can be shown from the density and seismic models.

Makran Subduction Zone
The Makran subduction zone (MZ) is a 1,000-km-long belt bounding the southern border of the Eurasian plate, beyond the southern end of the Zagros belt.In this area, the subduction of the Oman oceanic crust occurs.Makran has a W-E trend and is delimited on both sides by two N-S transform fault systems: to the west, the Minab Fault (MF), which separates the oceanic subduction from the continental collision between the Arabian and Eurasian plates; to the east, the Chaman fault system (CF) (Figure 2), which separates the Eurasian and Indian plates, being an important fault continuing northeastwards (Abdollahi et al., 2018).We can identify five different geological units into which we can divide the belt, each one separated by major thrust zones.These are from the north to the south, the North Makran, the Inner Makran, the Outer Makran, the Coastal Makran, and the offshore Makran (see Abdollahi et al. (2018)).The main differences are between the North Makran which preserves remnants of the period before the beginning of the subduction, and the central parts of the belts with sedimentary layers resulting from the Eocenic subduction.The North Makran hosts the Jaz Murian Depression (JMD, see Figure 1).We also report the presence of a volcanic belt in the northernmost part of the Makran area with an SW-NE trend of the three main volcanic centers, which are in order the Bazman, Taftan, and Sultan volcanoes (Figure 2).These volcanoes are the youngest representation of arc-related volcanism, still active in Iran (Abdollahi et al., 2018).

Iranian Magmatism
Magmatism related to the recent geologic history as well as to the subduction of the Neo-Tethyan ocean lithosphere is present in Iran (Omidianfar et al., 2023).Mostly, the magmatism is Cenozoic, and we can divide the For the information about the legend relative to the magmatism, we refer to Pollastro et al. (1997).Faults database from Styron and Pagani (2020).
principal outcrops present into three different magmatic belts: the UDMA, the Alborz Magmatic Belt (AMB), and the Sistan Magmatic Belt (SMB).About UDMA, we have already discussed in the previous chapters.It is one of the principal components of the Zagros belt, directed from north-west to south-east, composed of calc-alkaline rocks.For this reason, it is interpreted as an Andean-type magmatic arc, following the Neo-Tethyan subduction (Chiu et al., 2013).SMB is a magmatic area in the east of Iran, crossing LB and Sistan Suture zone.Rocks within the SMB include magmatic intrusive rocks from Jurassic and Late Cretaceous, and Late Cretaceous lavas (Chiu et al., 2013), while Cenozoic lavas and rocks form a massive magmatic sequence.Magmatism is not directly related to the Lut-Sistan collision, but could be the result of Eocene-Oligocene asthenospheric upwelling (Pang et al., 2013).The AMB originates from the collision of Iranian Block with the Eurasian plate, which continues in the present days with the ongoing intra-continental deformation caused by the convergence of the Arabian and Eurasian plates (Castro et al., 2013).It consists mainly of Precambrian to Eocene volcanic successions, which are intruded by Mesozoic to Cenozoic plutons (Castro et al., 2013).To make final considerations, we use the catalog of the magmatic outcrops from USGS (Pollastro et al., 1997), as shown in Figure 2.

Regional Velocity Tomography Model
The 3D shear wave velocity (Vs) model of Kaviani et al. (2020) is the data set used for calculating the starting crustal density model, obtained from the velocity to density conversion based on given empirical and physical equations, as explained below.The velocity model was obtained from group-velocity dispersion curves analyzing continuous seismic ambient noise from 709 broad-band stations, mainly distributed in the Turkish and Iranian territory, and regional earthquakes from 143 permanent broad-band stations of the Saudi National Seismic Network.On the Arabian Peninsula, continuous seismic ambient noise records were not available, so here the earthquake records contributed most to improving ray coverage across the area (Kaviani et al., 2020).
Data sets resulting from ambient noise and earthquake records overlap in the period range of 8-75 s, in the Saudi Arabia ray paths.To avoid problems due to the proximity of the measurements, all the analogous paths have been grouped and considered as single ones (Kaviani et al., 2020).Group-velocities have been then inverted to obtain a depth-dependent shear wave (Vs) model, performing the inversion for periods between 2 and 150 s.For all the other aspects of the tomography and the detailed information regarding the inversion technique used, we refer to Kaviani et al. (2020).From the tomography velocity model, we select our area of study, using a bounding box of 25°N by 40°N in latitude, and of 44°E by 63.4°E in longitude (Figure 3).Spatial resolution used for our final model is 0.2° (∼22 km), while maximum depth reached is 105 km, with a variable vertical resolution that starts from a minimum of 1 km in the first 60 km, then gradually increases to a maximum of 7 km in the deepest part.For details regarding the vertical resolution and the rest of the model, we refer to the "readme" file of the model.

Gravity and Magnetic Data
The gravity field used to perform the Bayesian inversion of the lithospheric density cube is calculated using the model XGM2019e (Zingerle et al., 2020).It is a global gravity field model developed by the Institute of Astronomical and Physical Geodesy of the Technical University of Munich, and represented in spheroidal harmonics up to maximum degree and order (d/o) of 5,399.This high d/o guarantees a nominal final spatial resolution of 2 arc min (∼4 km).Basically, it is a combined model, composed of the satellite model GOCO06s (which collects data from GOCE and GRACE missions) to cover degree and order up to 300, and by two different ground gravity data sets, a 15 arc min (∼28 km) and a 1 arc min (∼2 km) augmentation data set, to resolve the shorter wavelengths.The first, developed by NGA (US National Geospatial-Intelligence Agency), is for land and ocean gravity observation and is limited to d/o 719; the second one is derived from altimetry and topography observations and is used for the high-resolution of the model, resolving the band from d/o 719 up to 5,399.Consequently, in our study, we limit the maximum d/o to 719, corresponding to 27 km half-wavelength resolution (Barthelmes, 2013), coherently with the spatial resolution of our 3D volume of 0.2°.Starting from the XGM2019e model, we compute the Bouguer anomaly exploiting the rock equivalent topography model dV_ELL_RET2014_plusGRS80 (Rexer et al., 2016).The spherical harmonics coefficients of the models have been retrieved from the ICGEM web server (Ince et al., 2019).The gravity field was calculated at the height of 4,000 m above the ellipsoid.The gravitational effect of the deep mantle has been also removed by considering only spherical harmonic degrees larger than 12, according to Tadiello and Braitenberg (2021).Regarding the magnetic field, we use the global model EMAG2v3 (Meyer et al., 2017) as observable, developed by the NOAA National Center for Environmental Information.It was realized with the combination of satellite and magnetic observatory data, and airborne and marine magnetic measurements, with a resolution of 2 arc min (∼4 km), at an altitude of 4 km above the geoid.Differently from previous models, this grid is based only on interpolation of observed data and does not use a-priori information about geologic structure or ocean bottom age.This is because it has been observed that in the other global magnetic anomaly grids, this inclusion forced artificial linear patterns (Meyer et al., 2017).In Figure 4, we can see the final maps of the two data sets.For further information about the two models, we refer to the cited articles.

Sediment Thickness, Curie Depth Models, and Magnetic Susceptibility Distribution
In the inversion methodology we pursue, a division of the starting model into different layers is required, each one with its characteristics and values.Leaving aside for a moment the Moho discontinuity, and focusing inside the crust only, we decide to divide the model into four different layers.These crustal layers are the sedimentary basin and the crystalline crust, divided by the basement depth surface, and the magnetized and non-magnetized crust, separated by the Curie surface.To our best knowledge, an exhaustive and resolute sediment thickness model for the entire Iran is not publicly available to date.We have found localized sediment thickness models for limited areas, but an integration to formulate a comprehensive sediment thickness model was not feasible.We therefore adopt a thickness model defined by seismologic investigations, where a seismologic Vs velocity value of 3.1 km/s was used to define the basement depth variations (Irandoust et al., 2022b).Regarding the depth of the Curie surface, we choose to use the one proposed by Mousavi and Ardestani (2023), who define a 3D surface heat flow model in Iran, derived and calculated by geophysical and petrological observations, and who then extract the Curie surface corresponding to the depth of the nodes on the isotherm at 580°C. Figure 5 shows the maps of the two surfaces used with the relative depths.
With the magnetized crustal layer, it was necessary to suppose a magnetic susceptibility distribution.In detail, in the absence of direct information about the possible magnetic susceptibility of crust and sediments in the area, we consider a constant average value of 0.018 SI for the former (derived from Hemant and Maus (2005)) and of 0.005 SI for the latter Hunt et al. (1995).Modeling the Moho discontinuity, as will be described in the next chapter, we encountered intersections between the crust-mantle transition layer calculated and the available Curie surface.We decided to add a magnetized layer in the mantle, to which we assigned a starting susceptibility value in the a-priori model of 10 −5 SI.

Methods
A multidisciplinary approach is adopted to solve the problem of defining the density structure for Iran and surrounding areas.In this section, we explain all the processing stages starting from the initial tomographic data set and arriving to the new model of the Iranian lithosphere.The first step (step 1) is to fill void pixels in the regional velocity model aimed to obtain a regular 3D grid, with which we define a new map for the depth of the Moho (step 2).The calculated Moho is then used to divide the voxels of the velocity cube which are in the crust, from those that are in the mantle.This distinction is used in the third and fourth steps, namely in the velocity-to-density conversion, which requires different methodologies in crust (step 3) and mantle (step 4).
Step 5 begins with the 8 of 25 addition of mean susceptibility to the magnetized crust and magnetized mantle layers and continues with the Bayesian 3D inversion from which we obtain the final model of the Iranian lithosphere.To conclude, we also use the final inverted model and starting velocity tomography to calculate the shear modulus in the area (Figure 6).

Regional Tomography Interpolation on Regular Grid Through Ordinary Kriging
The first problem faced during the processing is that the Kaviani tomography model (Kaviani et al., 2020) contains cells with no values, which gives the model an irregular shape, difficult to process.To resolve this problem, we operate with a geostatistical approach, interpolating the regional tomography on a regular grid with a resolution of 0.2° applying an Ordinary Kriging technique (Wackernagel, 1995).In detail, a set of bi-dimensional interpolations, one for each height of our model, has been performed.The empirical variograms are modeled with a "Stable" function (Wackernagel, 1995).The result of the interpolation, for one slice, is shown in Figure 7.

Moho Undulation (Gradient Method)
The interpolated velocity cube is analyzed to define the depth of the Moho discontinuity.The chosen method studies the vertical gradient of the seismic velocity model (Kaviani et al., 2020;Tadiello & Braitenberg, 2021), seeking the largest vertical velocity gradient in the velocity interval identified as typical for the region above and below the Moho.First, we conduct a velocity analysis using the global tomography of Simmons et al. (2015),  where the Moho depth has been published on a grid spacing of 1°, to define the mean shear wave velocity at the Moho transition of the area.This analysis constrains the depth of the Moho to be located in a range of Vs velocities from 3.87 to 4.57 km/s.We therefore consider this range of variability as the set of velocities to be further investigated to define our Moho map.Then, for each node of the area, we define the depth with the maximum vertical velocity gradient, calculated in the range proposed.The depth of the node with maximum gradient is associated with the Moho depth.From this analysis, we find an evident similarity between the resulting Moho and the Moho proposed by Kaviani et al. (2020), mainly due to the similarity of the techniques used to derive the depth, both based on the velocity vertical gradient (Figure 8).We prefer using the recalculated Moho depth as it is consistent with our re-interpolated velocity model.

Vs-Velocity to Crustal Density Conversion
The abovementioned Moho depth allows us to separate the velocity cube into two main layers: crust and mantle.This division is necessary for the velocity-to-density conversion in crust and mantle separately, since the two conversions require different methods.Considering the crust, numerous equations have been proposed, most of which convert density from body waves (Vp).Therefore, it is necessary to first convert Vs from our cube into Vp, using the "Brocher's regression fit" (Equation 1) proposed by Brocher (2005), valid for Vs between 0 and 4.5 km/s, for crystalline and sedimentary rocks, estimated from boreholes and laboratory data. (1) For the density conversion from Vp, we consider the polynomial relation (Equation 2) defined by Brocher (2005), based on the "Nafe-Drake curve" theorized by Ludwig et al. (1970), who empirically found the parameters for the conversion valid for Vp between 1.5 and 8.5 km/s.The result is the following relation, which returns density values in g/cm 3 We find that the above equation can be suitably reduced to a fourth-order polynome, without loss of representativity of the experimental laboratory data reported in Brocher (2005).We repeat Brocher's best fit on his data using the Least Square method and we investigate the minimum polynomial order explaining the data by means of statistical inference.In detail, we ask ourselves whether, with a significance level of 95%, the hypothesis that each estimated coefficient of the polynomial is different from zero is satisfied (details about this statistic test can be found in Cazzaniga ( 2008)) (Figure 9).It turns out that, with the above confidence level, the fifth-order polynomial and the fourth-order one are statistically equivalent.We therefore choose to use the fourth-order polynomial, defining a new relation equally valid for the velocity-density conversion.
This last equation (Equation 3) returns density values in kg/m 3 , and is used for all the crustal density conversion steps.

Mantle Modeling (Perple_X)
The modeling of the mantle density requires another approach compared to the crustal modeling, due to the fact that the empirical relation of Brocher (2005) and therefore our Equation 3 are limited to crustal velocities and rocks.We compare the mantle velocities from tomography with synthetic upper mantle velocities and densities calculated with the open-source software Perple_X (Connolly, 2005).Perple_X is a collection of FORTRAN programs used for the calculation of phase equilibrium diagrams minimizing Gibbs free energy to map the phase relations (Connolly, 2005).The necessary starting point is a specified rock compositional model: for the selected zone, we take it from Tunini et al. (2014) assuming a standard mantle chemical composition expressed according to the NCFMAS system (Na 2 O -CaO -FeO -MgO -Al 2 O 3 -SiO 2 ), derived from global studies (Griffin et al., 2009).
In Table 1, the percentage of the elements used to build the model is shown.The software requires a few calculation parameters which we describe here.The thermodynamic model refers to Holland and Powell (1998), which is the default, and results in phase equilibria and elastic properties.Further necessary for the processing is the computation setting without saturated fluids or components, and the use of chemical potentials, activities, or fugacities as non-independent variables.Pressure and temperature are considered as independent of each other.We calculate the phase diagram in a range between 0.1 and 7 GPa for pressure, and from 500°C to 1,300°C of temperature.For the description of the default thermodynamic model and further computational details, we refer to the software documentation.At the end, the result is an equilibrium phase diagram, showing the stability fields for different minerals, depending on pressure and temperature.The software is also able to return several thermodynamic parameters, derived from the calculated phase diagram.From it, we extract synthetic seismic velocity and density values, related to a specific combination of pressure, temperature, and chemical composition.To derive the densities we use in the starting model, we compare the seismic velocities obtained from Perple_X with those in the tomographic model, using pressure values calculated at the given depth.The comparison gives also information on the temperature, since the synthetic velocities are both pressure-dependent and temperature-dependent.We start calculating the pressure at the Moho depth using the crustal densities previously defined,  (Brocher, 2005) with different number of polynomial coefficients.We find the best-fit parameters, for equations from fifth to second order by computing the Fisher test and calculating the acceptance degree of the new equations.
and then continue for greater depths in the mantle using the synthetic densities.At the end, we match mantle velocities calculated by the software with the velocities from the tomography to obtain the calculated mantle densities, relative to the tomographic depths.

Lithosphere Model (Bayesian Inversion)
Since the geometries as well as the density and susceptibility distributions of our 3D volume show large uncertainty due to the different hypotheses and approximations presented above, we exploited potential field data to improve our initial model.As outlined at the beginning of Section 4, the result of the previous steps together with the data reported in Section 3 are used to create an initial three-dimensional geological model, characterized by several layers (namely sediments, magnetized and non-magnetized crust, magnetized and non-magnetized mantle) each one with its own density and magnetic susceptibility distribution.For all these quantities, we also define

Shear Modulus Calculation
The inversion produces the final model of the lithosphere in Iran, with density and susceptibility values.These final cubes, joined to the starting tomography, allow to calculate the elastic shear modulus, and therefore investigate the rheology of the area.The equation relating shear wave velocity and density values, allows the calculation of the shear modulus as follows: We express the shear modulus by the product of the square of the shear wave velocity tomographic model and the final density model (Equation 4).With this simple equation, we are able to obtain a 3D model of the lithosphere rigidity.

Results
We now show the inversion results, starting to quantify the reliability of the final model with the residual maps.In Figure 10, the residuals after the Bayesian inversion are plotted, calculated from the differences between the gravity and magnetic fields used as observations, and the gravity and magnetic fields generated by the final model.Considering the two maps, our final model fits gravity data with a standard deviation of ∼7 mGal and the magnetic data with ∼15 nT.Almost all anomalies are resolved by the model, even if some localized anomalous high-residual patterns remain.These can be explained by the presence of density and susceptibility superficial bodies (in the first 10 km), not defined in the a-priori model and therefore undetectable from our inversion.To verify this condition, we perform some tests starting from the obtained final model trying to invert only the first 10 km, to define what the amount of density and susceptibility variations would be required to explain the Bayesian residuals.This estimate is to be considered in general terms only, as a detailed inversion of these upper structures would require more constraining data.This second simplified inversion solves all the residuals (reducing their standard deviation to less than 1 mGal and 4 nT for the gravity and magnetic field, respectively) identifying the presence of different superficial bodies.These density variations reach up to ∼150-200 kg/m 3 in localized extreme spots.The same condition is valid for the magnetic residual, almost entirely solved by the presence of two isolated magnetized bodies (∼0.04-0.06SI).Maps of the estimated density and susceptibility distributions are reported in Supporting Information S1.
With this kind of inversion, we also obtain a datacube with the uncertainties of the final models.We give below the accuracy calculated for each layers: for the density model, uncertainties have a maximum value of 19 kg/m 3 for the sediments layer, 29 kg/m 3 for the crust layer, and 14 kg/m 3 for the mantle layer, while considering the susceptibility inverted model, uncertainties reach maximum value of 0.0015 SI for the magnetized sediments layer and 0.006 SI for the magnetized crust.Maximum uncertainty values for the density model are concentrated in the South-West Caspian Sea area and in the Amu-Darya Basin (ADB), while in the susceptibility model, uncertainties are more distributed within the Iranian plateau, with higher values in the ECIM and in the Sanandaj-Sirjan Zone (SSZ).From the resulting inverted model, we extract the new interfaces between the modeled geological horizons: the inversion is allowed to modify the input layer in the bounds dictated by the uncertainty of the layer depth, expressed in Section 4.5.In the final model, the position of the sediment basins is quite changed, making some corrections on the depths, which are generally less deep (Figure 11a).The SCB results collectively shallower and linked in a more continuous basin, with reduced depths extending southwards, and connecting eastwards to the West Turkmenistan Basin (WTB) and westwards to the Kura Basin (KB).The ADB became deeper and occupied a larger area of the Turkmenistan Depression (TD), while considering the Persian Gulf, we have an important reduction of the depth along the entire area.The Mesopotamian Foredeep Basin (MFB) passed from 26 km in the original model to a maximum of 15 km, and its depocenter is shifted southwards compared to the original, while in the Makran area (MZ), the basins present in the starting model almost disappear, even if we can identify two small depocenters (Bandar Abbas Basin [BAB] and Hamun-Mashkel Basin [HMB]).The Moho, marking the interface between the lower crust and upper mantle, has systematically shallower depth along the entire Zagros belt, which seems to have a more superficial crustal root, while the rest of the area basically is uplifted, even if the main features were conserved (Figure 11b).Finally, considering the Curie depth surface, in the final model, it is basically flattened, with the deeper nodes being uplifted and the more superficial parts of the model being lowered (Figure 11c).Also, in this case, we can retrieve the same features between the initial and final model, with the entire Zagros belt (ZFTB, SSZ, and UDMA) and part of the ECIM that has a superficial Curie depth surface, which becomes deeper under the SCB and in TD.A shallow Curie depth for the thick crust below Zagros is expected due to a greater radiogenic heat production of the thick crust, compared to a thinner crust.The thick crust has a greater thickness of radiogenic heat-producing rocks, leading to higher crustal temperatures and a steeper temperature gradient with depth, compared to a thinner crust.The consequence is a shallowing of the temperature of isothermal surfaces, and therefore also a shallowing of the Curie depth.The Curie depth is tied to the crustal temperatures, with a shallower Curie depth for higher temperatures.

General Description of the Inverted Density Model
We describe the inverted density model, which defines the 3D density distributions of the sediment, crust, and mantle layers, illustrating density maps at different heights between 10 and 70 km (Figure 12).At 30 km, high-density values in the eastern SCB start to appear, due to the fact that we are approaching the mantle in that area at this depth, and a new localized high-density body (∼2,850 kg/m 3 ) is also found in the Straits of Hormuz (SH), marking the transition between the southern end of the continental collision of the Arabian Plate to the oceanic subduction of the Makran area.The rest of Iran remains quite consistent with the superficial slice, with the low density along the ZFTB and Alborz (AB), and the high-density area in CI and YB.The Dasht-e Lut desert (DLD) high density is still visible, with the exception of the reduced density in the Birijand Block (BB), connecting to the SSZ, and an area protruding westwards to the TB.The Afghan Block (HB) extends the high density down to 30 km depth, defining a separate crust compared to the Sistan suture (SS).At 50 km depth, we are entirely in the mantle, and the density values range between 3,200 and 3,350 kg/m 3 .By qualitatively analyzing the density distribution, we observe a trend reversal along the Zagros belt, where the previously higher density areas along the SSZ and the UDMA become lower density than the ZFTB.A low-density anomaly, as opposed to the upper layer, is also visible in the Strait of Hormuz (SH), while the rest of Iran becomes quite homogeneous.Higher-density areas are located in the southeastern part of the SCB, in northern TD, and in the western Persian Gulf.At 60 km depth, the density is more homogeneous over the entire Iran, with lower density bodies in the LB and in the SH, and an extended higher density area in TD and in the Persian Gulf.These features are reinforced at 70 km depth, with the high density of the Persian Gulf that invades the ZFTB, and the studied area acquires greater inhomogeneity: the Central Iran province (CI) is separated into two opposing blocks, one less dense west of Teheran and one more dense south of the Aladagh-Binalud mountains (ABM) and LB and TB become less dense in comparison to the rest of the East Central Iran microplate (ECIM).The only contrast to the previous mantle slice is in the SH, where we have an inversion compared to the lower densities of the upper layer.

The Inverted Density Model in Relation to the Geologic and Tectonic Structures
We make some considerations regarding the density volume generated and described in the previous section, relating it to the known geological and tectonic aspects in the area.In particular, for this analysis, we focus on the first 20-30 km of the model and consider the average density of the crustal layer.Average crustal density is calculated considering the density column for each node in the grid, using only the density values related to the voxels in the crystalline crust layer.First, we consider the higher-density anomalies in the average crust, and compare them with the position of the magmatic bodies previously presented, extracted from the USGS catalog (Pollastro et al., 1997).
We can see that some trends have very similar shapes, even if they do not correlate perfectly.Regarding the mean crustal density (Figure 13a), the high-density lineaments of the Alborz area (AB) correspond to the presence of Tertiary magmatic bodies, and some high-density patterns are also detectable along the UDMA, although these do not correspond in an optimal way with the map of magmatic outcrops.The situation changes if we consider the mid-crustal slice at 20 km (Figure 13b).Here, we can appreciate how all the magmatic features correspond to high values in density, in particular along the UDMA and Alborz area (AB), where some similarities were already visible, now the two maps are perfectly related.In this 20-km slice, additional high positive anomalies are visible also in the JMR south of KD which correlates with Tertiary and Mesozoic magmatic bodies, and in the Lake Urmia (UL) volcanic area, in north-western Iran, which is affected by Cenozoic magmatism.The only anomaly that persists in both average values and in the 20 km layer is in the SSZ, where the magmatism in the area seems to correlate with low-density values compared to other areas.Turning to more tectonic aspects, considering again the average crustal density, we observe that high-density lineaments which stand out from the surrounding areas, almost overlap with the locations proposed by literature for the Paleo-Tethys and the Neo-Tethys suture  lines (Irandoust et al., 2022a;Wan et al., 2021).We therefore try to draw a density-derived suture, calculating the horizontal gradient of the average crustal density and then following the anomalous lineaments.Our suture lines are proposed in Figure 15, based on the gradients shown in Figure 14.
In particular, our Paleo-Tethys suture follows the KD, Alborz (AB), and then proceeds to the northwest following the Caucasus, while the Neo-Tethys suture extends from the north-western Iranian border along the Zagros belt until the Minab fault (MF, see Figure 2), after then it moves south toward Makran subduction (MZ).The Paleo-Tethys position seems giving a high-density trend following the northern border of Iran, while the Neo-Tethys marks the separation between a high-density area northeast of the suture and a low-density area in the opposite part.We also note from the gradient the position of another suture along the SSZ marking the differences between the lower densities in the Iranian plate and the higher densities in the Afghan plate, which, according to Bröcker et al. (2022), we associate with the Neo-Tethyan suture, that after continuing parallel to the Makran Zone, crosses Iran again, along the Sistan Suture Zone.Focusing on the eastern Iran, differences in density can be identified, reflecting the tectonic lineaments between the LB and the SSZ.In particular, a high-density body is found in the southern part of the LB, corresponding to the DLD surrounded by lower-density zones, which  follow up with the lineaments of the SSZ.In Figure 16, it is clearly visible that the strike-slip faults limiting the western and eastern LB mark the separation between the higher-density nucleus of the LB and the lower-density zones.The lower density body in the Birjand Block (BB) seems to rotate around the central dense body of the Lut Block (LBDB), rotation of the block that is known in literature (Rashidi et al., 2022), and that is also explainable observing the orientation and movement of the strike-slip faults in the area.Our model suggests the presence of a crustal dense and stable body in the LB, which acts as block that guides the crustal movement to the east and west of it in response to the compressive stresses present due to the Makran subduction.The principal compressive stress orientation along the right-lateral strike-slip fault system along the western margin of the LB is estimated to be N15E (Rashidi et al., 2019).

Description and Consideration of the Mean Crustal Susceptibility
We now describe the map with the average crustal susceptibility resulting from the inversion (Figure 17).Like the average crustal density, average crustal susceptibility is calculated by averaging the column of the susceptibilities for each node from the inverted model, only considering the values included in the magnetized crustal layer.The inversion detected a series of magnetized bodies within the crustal layer, in particular in certain areas along the SSZ and the UDMA, and in the western part of the Makran accretionary prism (MZ).Areas with considerable magnetization are also detected in the KD, in the eastern part of the SCB, in the West-Turkmenistan Basin (WTB), and in the East Central Iran microplate (ECIM), in particular in the northern part of the LB.The entire Central Iran seems to have no magnetization as well as the ZFTB.Also, the dense crustal body in the LB, found and described in the previous chapters, appears without a significant magnetization.Compared with the magmatic outcrops catalog also used for the considerations regarding the density model, we partly observe a correlation between magmatic bodies and magnetization.In general, the southern UDMA seems to not have relevant magnetization, like the magmatic belt along the Sistan suture, even if also in these apparently unrelated areas, it is possible to observe some small correlated spots.Considering the northern part of the UDMA, we found an important link between the position of the outcrops and high magnetization, in particular on the border with Azerbaijan, in the Urmia Lake volcanic area, as well as in a central segment of the SSZ.In the west of Makran Zone (MZ), a magnetized area could be associated with the presence of lineaments of ophiolitic outcrops, on the border with the JMD.Also, the Alborz magmatic belt gives some small areas with high magnetic trends in accordance with the presence of outcrops.In summary, most of the magmatic presence is evidenced by the magnetic anomalies in the area, but there are many high anomalies (KD, SCB, and WTB, Eastern Central Iran ECIM) that are not explained by a plain presence of rocks outcropping at the surface.Proceeding to the depth of 20 km (Figure 18b), the low rigidity anomaly moves toward north-east, maintaining the same parallel orientation to the Zagros, but being more aligned with the central axis of the ZFTB.Also, the earthquakes in this case perfectly match the low rigidity area.The rest of Iran becomes more homogeneous, and most of the rigidity anomalies previously presented are reduced, while CI becomes more rigid, in particular south of Teheran, in the Qom area, and south of JMR.To the East of the JMD, located in the forearc of the Makran subduction, the area becomes more ductile, merging northwards into the SSZ.The SCB also becomes a very ductile area, even if no earthquakes match with the zone.At 30 km (Figure 18c), the major part of the anomalies in the previous slice are preserved, and become more marked.The minor band along the Zagros and Makran foreland basin is interrupted by an anomalous very rigid spot in the SH.The SSZ shows increasingly marked differences in rigidity compared to the surrounding areas, with the ductile body from the Birjiand Block (BB) crossing the ECIM to the TB, while a high rigidity zone starts to appear in the YB.High rigidity anomalies are still visible in CI, south of Tehran, and south of JMR.Entering into the mantle, at 50 km (Figure 18d), we observe that the low rigidity anomaly continues its displacement toward north-east, causing a sort of reversion from the crustal trend, since the area along the SSZ and the UDMA that previously was the most rigid, now becomes the only ductile part in the area.The remainder of Iran is quite homogeneous.Considering the entire model, this low rigidity anomaly seems to form a ramp in a northeastern direction extending toward increasing depths.
Considering the average crustal shear modulus, the seismicity in general correlates with the regions of reduced rigidity (Figure 18e), with the low rigidity perfectly matching with the faulting, and with the location of the superficial earthquakes, with depths from 0 to 30 km, and covering the interval from 1905 to 2021.In particular, this correlation is important along the ZFTB, the SS, and the ECIM.The SCB does not show this correlation with the earthquakes location, because it has values of reduced rigidity, but the seismicity follows its border, and divides the southern from the northern Caspian Basin through the seismicity of the Aspheiron ridge (AR).The map of the final average crustal shear modulus is more coherent compared to the initial shear modulus shown in Supporting Information S1.

Discussion
To our knowledge, this work proposes the first well-defined high-resolution 3D density model of the entire Iran, for crust and lithospheric mantle down to a depth of 105 km.The results are based on one of the most recent 3D seismic tomographic models publicly available, and therefore should guarantee the maximum resolution that today is possible to be achieved in the public domain (Kaviani et al., 2020).A more recent seismic tomography was published after we started the inversion processing, but the 3D tomographic model is unavailable (Irandoust et al., 2022a).
The inversion of the gravity field in the Caspian and Alborz has been reported by Motavalli-Anbaran et al. (2013) and along profiles crossing our area in Motavalli-Anbaran et al. (2011).The method uses gravity, geoid, and topography to constrain a forward modeling or a Bayesian inversion which aims to define Moho depth, lithosphere-asthenosphere boundary (LAB), and average crustal density.Thermal modeling through solving the steady state heat transport equation constrained by gravity, geoid, isostatic equilibrium condition, and heat flow values has been made along a north-south profile extending from the Makran to the Turan platform.The density of the mantle lithosphere is temperature-dependent and controls the modeled geoid (Entezar-Saadat et al., 2017).The approach is different from ours, because the reduction in ambiguity of the inversion is achieved by fixing a-priori parameters as the LAB temperature, the assumption of linear variation with depth of lithosphere density, and the assumption of an average density for the crustal column including sediments, mid and lower crust.A further constraint is the local Airy-type isostatic compensation of topography through the underlying density model.Although the methods and assumptions are different, many features are common, like the thickening of the Moho below the SSZ to over 50 km, the shallower Moho around 40 km over the Turan platform, the CI with a shallower Moho compared to the Zagros around 40 km, and the shallow Moho below the SCB.Our average crustal density is consistent with Motavalli-Anbaran et al. ( 2013), apart from some local deviations.
Due to the rock-physical relation connecting seismic velocity and density, it is expected that the higher velocities in the tomographic model will correspond to higher density, an observation that our model allows to quantify, bringing further insight into the Iran lithosphere.The higher spatial resolution of the density and magnetic model compared to the seismic model is given by the higher resolution with which the fields are available, and through which particularly for small-scale and upper-crust features the density model brings a refinement of the tomographic model.This brings to a clearer match of geological structures defined by age, geodynamic evolu tion, and inheritance of rock formation, with the density and susceptibility variations, especially in the upper crustal layers.The Bitlis-Zagros suture (BTS) marks the transition between the Arabian and Eurasian plates, and therefore is expected to produce a major variation in rock properties.The suture marks a density boundary over the entire crustal thickness, with systematically lower density in the ZFTB down to the Moho depth to 40-50 km (up to 200 kg/m 3 density difference), and an inversion of the sign of the difference in the mantle down to 80 km (40-50 kg/m 3 ), a depth after which another inversion occurs, with the Arabian lithosphere that returns less dense than the Iranian lithosphere (30-40 kg/m 3 ) down to the bottom of the model.The suture marks also a difference in the crustal rigidity, with systematically lower rigidity in Zagros to the SW of the suture, and higher rigidity to the NE of it.This rigidity discontinuity shifts toward the UDMA for increasing depths, up to the Moho (Zagros rigidity 17 GPa at 10 km to 27 GPa at 30 km, SSZ and UDMA 38 GPa at 10 km, to 42 GPa at 30 km).Below the Moho, the rigidity matching the suture is limited to a narrow band below the SSZ and UDMA, with values that are systematically higher below the Zagros and the CI.This systematic difference is not visible in the magnetic susceptibility variation.The higher rigidity of CI could be the cause of the lower occurrence rate of significant earthquakes, which are found in the lower rigidity crust.The Zagros has a very high density of earthquakes, up to the higher rigidity margin, but only in the northern half of Zagros, the southern half having epicenters limited to the external (southwestern) parts of the Thrust and Fold belt, whereas the areas closer to the suture have low seismicity.One reason for the absence in this deformable zone could be the presence of the Late Precambrian Hormuz salt layer overlying the basement (Jahani et al., 2009), seen even in satellite images through the folding that salt diapirs generate, and which could favor deformation along detachments.The ECIM, an amalgamation of four blocks, shows density and rigidity variations that cross the blocks, and differentiates rather than the blocks, subunits of the blocks.The SMB, shared between the north SSZ and north LB, is systematically less dense over the entire crustal column, compared to the remainder of the mentioned blocks.Below the Moho the density values are homogeneous, and units composing the East Central microplate do not show up.The rigidity variations are similar, with the SMB of the microplate being less rigid in the crust.The layer below the Moho is differentiated, as the southern YB is less rigid.We find that the higher crustal rigidity controls the development of large basins as opposed to topographic ranges, a rule we find valid for large flat desert areas such as the Turkmenistan Depression, the Central Iran Dasht-e Kavir, the Dasht-e Lut (DLD), and the Hamun Jaz Murian.This observation is valid when considering the rigidity values below 10 km down to the Moho.The Rigidity at 10 km has only broader features, and does not show this correlation with the flat topography development.
Considering the susceptibility variations in the crust resulting from the inversion, the areas with the most intense magnetized bodies are identified in the MZ, in the northwestern part of the UDMA, and in the north of ECIM.
It is difficult to match our results with two previously published susceptibility models, due to the different approaches used to explain the magnetic anomalies with the susceptibilities.One approach is to define the average susceptibility over the crustal thickness through spectral methods (Teknik et al., 2020), which inherently will give a different result compared to ours.In fact, our model constrains the magnetization to be above a thermally defined Curie depth.A shallow Curie depth requires a reduced value of the susceptibility to explain the anomaly, since the magnetized layer has a smaller volume, with respect to a deeper Curie depth.This is a reason for differences in our results with the model of Mousavi and Ebbing (2018), due to some significant differences in the a-priori Curie depth.For example, in the MZ the Curie depth of Mousavi and Ardestani ( 2023) is about 30-35 km, whereas Mousavi and Ebbing (2018) has a deeper value of over 50 km.The much deeper Curie depth leads to a lower value of susceptibility.Another difference between the conceptual setup of our model and the one of Mousavi and Ebbing (2018) is the depths for which the susceptibility is allowed to vary.In our case, the susceptibility is inverted from the surface to the Curie depth.In the final model of Mousavi and Ebbing (2018), the susceptibility is inverted between base sediments and a constant boundary at 20 km depth, with constant low susceptibility value below 20 km down to the Curie depth (0.01 SI).This may be the reason our values are on average smaller compared with those found in Mousavi and Ebbing (2018), because our magnetized layer is generally thicker than 20 km.In detail, the area of the MZ results is more magnetized than the JMD, a trend that appears reversed in the cited models.However, we note that all three models perfectly border the position of the ophiolites present in the MZ.Another similar situation occurs in ECIM, where our model detects an increased magnetization bordering the position of magmatic bodies, anomaly that in Teknik et al. (2020) and Mousavi and Ebbing (2018) is extended in CI, incorporating the JMR.For the remainder, the anomalies identified by the inversion are quite consistent with the above two models, in particular along most of the UDMA and in the Urmia Lake (UL).In general, where the Curie depths are comparable, as in the UDMA, the susceptibility values are also comparable.

Conclusion
The target of this work is to realize a model of the lithosphere in Iran performing a gravity and magnetic joint Bayesian inversion, to study the different density and susceptibility distributions and bind them to different geological and tectonic aspects.We also wanted to calculate, starting from previous results, a shear modulus model to investigate the rheological aspects of the area.To realise the model we start from an available regional tomography model and through a series of processing steps that led to the formation of an a-priori model, including the identification of the different layers and a first transformation from velocity to density using empirical relation, the inversion delivers a reliable and interpretable model.Thanks to this model, we are able to make some tectonic and geological considerations about the lithosphere of Iran.In particular: • the density distribution in the upper crust reflects for most areas the trend of magmatic lineaments; • position of Neo-Tethyan and Paleo-Tethyan sutures can be detected through the density model, as it generates a strong and measurable gradient signal over the entire crustal thickness; • the faulting system in the SSZ and the rotation of the BB are enhanced by the presence of a dense and rigid body within the LB, around which the faulting and orientation of structural lineament occur; • the rigidity model shows a more ductile zone under the Zagros, moving northeastward with increasing depth, in which most shallow earthquakes occur.• the calculated density and susceptibility cube are an innovative data product that are published alongside the paper, and find various applications, including starting rate calculation, if a stress field is given as input.The Curie depth model based on the magnetic inversion is an input for new thermal calculations.For the University of Trieste, the authors acknowledge funding through PhD grant (GM) PON-Green thematics-DM1061 (10.08.2021); further for funding through project MIUR-PRIN2017 and ESA "Quantum Space Gravimetry for monitoring Earth Mass Transport (QSG4EMT)."This work was partially funded by Geomatics Research & Development srl (https://www.g-red.eu/) in the framework of the ESA "eXperimental jOint inveRsioN" project."The authors thank the reviewer Lu Li and an anonymous reviewer for the meticulous reviews that contributed to improve our manuscript.Paolo Ballato is thanked for providing useful scientific references on Iran.

Figure 3 .
Figure3.Map with the area covered by the regional tomography (black area).The green box is the area selected for this work.

Figure 5 .
Figure 5. Depth of base sediment (a) and Curie depth (b) layers used in the realization of the starting model.

Figure 7 .
Figure 7. Example Vs tomography slice (5 km) that illustrates the performance of the Kriging interpolation.On the left, the slice from the starting cube (a); on the right, the result after the interpolation (b).

Figure 6 .
Figure 6.Flowchart of the processing operations.

Figure 8 .
Figure 8. Moho depth proposed by Kaviani et al. (2020) (a) and recalculated from the tomography model through the maximum gradient method (b).
ranges of variability by means of a careful analysis of the available literature on the area(Irandoust et al., 2022b;Kaviani et al., 2020;Mousavi & Ardestani, 2023;Mousavi & Ebbing, 2018) and considering the uncertainties in the starting seismic velocity model and in the conversion to density: density and susceptibility range of variability used have values of ±30 kg/m 3 and ±0.005 SI for the sediments layer, ±50 kg/m 3 and ±0.018 SI for crust layer, and ±20 kg/m 3 and ±10 −5 SI for the mantle layer, while considering the variability of the surface, we adopt values of ±12.4 km for the basement, ±8 km for the Curie, and ±15 km for the Moho depths.We use the joint Bayesian inversion described inSansó and Sampietro (2021) andSampietro et al. (2022) to modify this initial model in such a way that its gravity and magnetic responses fit the observed potential fields, while always preserving the original intrinsic characteristics of the model.In other words, the algorithm, which has already been tested in several synthetic and real case studies such asCapponi et al. (2022),Sampietro et al. (2023), andSampietro and  Capponi (2023), modifies at the same time the initial geometries and the distribution of density and susceptibility values, within the confidence intervals of the a-priori model to fit the observed gravity and magnetic anomaly data.For a more detailed description of the Bayesian inversion algorithms, we refer to Supporting Information S1 and to the above references.

Figure 10 .
Figure 10.Gravity (a) and magnetic (b) residuals before the inversion, and gravity (c) and magnetic (d) residuals after the inversion.The maps represent the differences between the original observation and the forward calculated fields of the starting model (a, b) and the final model (c, d).

Figure 12 .
Figure 12.Slices at different depths representative of the final density model.Geological provinces and blocks are also indicated with their acronyms.For references and acronyms explanation, see Figure 1.
Starting at 10 km, we can see the presence of the sedimentary basins (the ADB on the Turan Platform, the SCB, the MFB, and the BAB), with lower values in Central Iran, and also the sedimentary rocks making up the Zagros orogenic belt, the salt deposits in Southern Zagros, and the JMD.High-density trends are detected along the Sanandaj Sirjan Zone (SSZ) and the UDMA, in the Yazd Block and in the Helmand Block.More localized increased density areas are found northwest of the SCB, in the northwestern part of the Paleo-Tethyan suture.A last increased density is found in the southern part of the LB, in the DLD.At 20 km depth, the pattern changes, and the areas with increased density are less extensive, while the low dense area along the Zagros Fold and Thrust belt (ZFTB) and the East Central Iran microplate (ECIM) are more noticeable.In Central Iran (CI), we find two high-density spots, in the Qom area south of Teheran, and south of the Jaghatai Mountain Range (JMR), separated by a low-density area.Other high-density trends are found at the northern extreme of the Paleo-Tethys suture following the suture southeastwards along the Alborz (AB), and from north-west to south-east boarding the Urumieh-Dokhtar Magmatic Arc (UDMA) and the SSZ, until the JMD.Even the high in the southern LB persists in the deepest layers, while around it the SSZ has a lower density.

Figure 14 .
Figure 14.Density gradient (a) calculated from the final model, used to track the position of the Paleo-and Neo-Tethyan sutures (b).

Figure 13 .
Figure 13.Average crustal density (a) and density slice at 20 km depth (b) with superposed the distribution of the magmatic outcrops.The greater part of the high-density trends follows the position of the magmatic lineaments.The acronyms are defined in Figure 2.

Figure 15 .
Figure 15.Average crustal density and position of the three founded sutures, derived from the final model.

Figure 16 .
Figure16.Focus on the Central East Iran zone, comprising the Lut Block (LB) and the Sistan Suture zone (SSZ).The Sistan Suture Zone, in particular the less dense area of the Birjand Block, wraps around the central high-density body of the Lut Block (LBDB).For references and acronyms explanation, see Figure1.

Figure 17 .
Figure 17.Average crustal susceptibility with geological provinces and magmatic outcrops.It is difficult to identify correlations between magmatic outcrops and high crustal susceptibility.For references and acronyms explanation, see Figure 1.Faults database from Styron and Pagani (2020).

Figure 18 .
Figure 18.Shear modulus model at different depth slices.Geological provinces are also indicated (a-d).(e) Crustal average shear modulus model and the position of the superficial earthquakes, from 0 to 30 km hypocenter depth.For references and acronyms explanation, see Figure 1.Faults database from Styron and Pagani (2020).Earthquakes database from the U.S. Geological Survey (2021).

Table 1 Mantle
Chemical Composition Used for the ModelingFigure 9. Test of Brocher's relation