Crustal Structures From Receiver Functions and Gravity Modeling in Central Mongolia

3D forward gravity modeling combined with receiver function (RF) analysis characterizes the crustal structures of the southern part of the Mongolian collage. The seismic signals of the 48 stations of the MOBAL2003 and the IRIS‐PASSCAL experiments were analyzed to get the RFs. This analysis revealed a significant difference between the crustal structures of the Hangay dome and the tectonic zones in the south. In addition, seismic stations south of the Hangay dome display significant signals related to the occurrence of a low‐velocity zone at lower crustal level confirmed by the gravity anomalies. Finally, these seismic analysis inputs, the boundaries, the lithologies, and the density values from rock samples of the different tectonic zones constitute the starting points from the 3D forward gravity modeling. The resulting crustal density model indicates: (a) the likely absence of a Precambrian basement block beneath the Hangay dome, (b) an alternation of two low‐velocity/low‐density zones (LVLDZs) with high‐density zones in the Baydrag microcontinent interpreted as fragments of early Tonian plutons, (c) the occurrence of an LVLDZ at the lower crustal level beneath the Lake zone, the Mongol‐Altai Accretionary Wedge, and the Trans‐Altai Zone. Therefore, the combination of the seismic RF with gravity analysis and modeling reveals new crustal structures of the Mongolian collage and enhances the occurrence and the extent of an LVLDZ at lower crustal level. These LVLDZ may demonstrate the existence of the relamination of a hydrous material in southern Mongolian collage.

GUY ET AL.
10.1029/2023JB027614 3 of 28 gravity and seismic studies over our study area.Their lithospheric scale was mostly dedicated to the Mesozoic and Cenozoic volcanism and the investigation of the timing and the source of the doming in the Hangay.

Tectonic and Geological Settings
The representative tectonic units of the CAOB located south of the Siberian craton constitute the Mongolian collage (Xiao et al., 2018), whose southern part includes Central Mongolia.It can be divided into six tectonic zones (Figure 1): (a) the Mongol-Okhotsk domain in the north is composed of Silurian to Carboniferous oceanic sediments, which are surrounded by (b) the Mongolian Precambrian continental blocks including to the north the Tarvagatay, to the west the Zavkhan and to the south the Baydrag blocks; (c) the Bayankhongor Complex is intercalated between the southern part of the Mongol-Okhotsk domain and the Precambrian blocks, and it is composed of an ophiolitic mélanges obducted within Neoproterozoic and Cambrian accretionary wedges; (d) the Lake Zone (LZ) in the south of the Baydrag block, which is composed of Cambrian-Ordovician arc, also called Ikh-Mongol Arc, developed on Late Proterozoic accretionary complexes; (e) the Mongol-Altai Accretionary Wedge (MAAW) in the south being Cambrian-Ordovician in age; (f) and the Trans-Altai Zone (TAZ) in the southernmost part of the study area, which is an oceanic domain, composed of Devonian to Carboniferous volcano-sedimentary sequences.
The Bayankhongor Complex is composed of ophiolites formed at 569 ± 21 Ma (Kepezhinskas et al., 1991) and accreted at 540 Ma, intercalated in between a Silurian passive margin in the north (Buchan et al., 2001) and a Neoproterozoic to Ordovician accretionary prism in the south (Buchan et al., 2001;Jian et al., 2010).However, the tectonic context of the emplacement of the ophiolites is still under debate as Buchan et al. (2001) demonstrated that the ophiolites mark the suture between the Baydrag and the Hangay microcontinents, whereas Osozawa et al. (2008) said that the ophiolites were obducted in an accretionary complex.The south-dipping South Hangay sinistral fault zone is interpreted as a major crustal boundary between the Bayankhangor Complex and the Baydrag microcontinent (e.g., Buchan et al., 2001;Comeau et al., 2018;Walker et al., 2007).
The Baydrag microcontinent is part of the Mongolian Precambrian blocks.It contains Proterozoic-Grenvillian tonalitic granulites and granitic gneiss that are unconformably overlain by Neoproterozoic to early Cambrian metasediments, intruded by Cambrian-Ordovician metaigneous rocks, overlain by Devonian to Permian volcaniclastic rocks (Buriánek et al., 2017;Demoux, Kröner, Badarch, et al., 2009).It formed in Paleoproterozoic and recorded the formation of Rodinia, its likely break-up, the formation of Gondwana and the Palaeo-Pacific subduction initiation, with the succession of the early Tonian magmatic arc, followed by the late Tonian passive margin deposits, and the Cambrian-Ordovician continental magmatic arc developed on the Precambrian basement (Buriánek et al., 2017;Ling et al., 2021).At the western and northern periphery of our studied area stand, respectively, the Zavkhan and the Tarvagatay Precambrian microcontinents, which constitute the continuation of the Baydrag microcontinents.The Zavkhan microcontinent is also characterized by Proterozoic high-grade igneous rocks and Neoproterozoic to Carboniferous metasediments and sediments (Bold et al., 2016;Kozakov et al., 2014;Soejono et al., 2023) and as well as for the Tarvagatay (Badarch et al., 2002;Kozakov et al., 2011;Kröner et al., 2014).The Cambro-Ordovician Ikh-Mongol Arc developed at the southern continental margin of Baydrag and constitutes a part of the LZ (Janoušek et al., 2018).The relics of this arc are assumed to be spread over more than 1,800 km along the Precambrian continental blocks (Figure 1a).
The three major faults we briefly described above, namely, the South Hangay, the Bogd, and the Trans-Altai faults, have poor determinations in their depth continuities.However, many studies concluded that these faults are reactivated and developed along older Paleozoic weak zone structures, which probably imply lithospheric fault features (e.g., Calais et al., 2003;Cunningham, 2001;Cunningham et al., 2009;Walker et al., 2007).

Crustal Architecture From Previous Seismic and Gravity Studies
The results of the MOBAL2003 seismic experiment mostly characterized the deep lithospheric structures of central and southern Mongolia.Its RF analysis demonstrated a constant crustal thickness of ∼44 km in the southern part of the Hangay dome (Mordvinova et al., 2007).Moreover, the asthenospheric upwelling beneath the Hangay dome explains its long-wavelength topography (Petit et al., 2008).A 3D forward modeling of the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite gravity gradient data also revealed a thick crust with a constant thickness south of the Hangay dome (Guy et al., 2017).The Hangay dome was studied by gravity-based, seismic-based, and topography-based models, which revealed strong velocity contrast within the crust and an anomalous low-density body at the Moho (i.e., in the uppermost mantle and/or lower crust) beneath the Hangay dome (Petit et al., 2002;Tiberi et al., 2008).These studies proposed that the Cenozoic uplift of the Hangay dome and its peripheral volcanism are most probably due to a mantle plume and is probably not driven by the India-Asia collision (Mordvinova et al., 2015;Petit et al., 2002;Tiberi et al., 2008).
The crustal structures and the major contacts between the different tectonic zones were determined along three N-S-oriented 2D gravity profiles (Guy et al., 2015;Zorin et al., 1993) and recently (Figure 2), the existence of a relaminant in southern Mongolia was demonstrated based on gravity and magnetic forward modeling constrained by whole-rock geochemical data (Guy et al., 2015(Guy et al., , 2021)).

Seismic Data Set
The MOBAL 2003 data (MOngolian-BAikal Lithosphere) were acquired in the frame of a collaborative seismic experiment, between Russian, Mongolian, and French institutions.A temporary seismic network of 19 broad-band seismic stations was deployed from the south of the Siberian craton to the Gobi-Altai range in southern Mongolia, building a 1,000-km-long profile (Figure 1a).The Central Mongolia Seismic Experiment (CMSE) consists of 112 broad-band seismic arrays that cover 900 × 600 km 2 in two separate deployments between 2012 and 2016 for a cumulative 21-month record.This experiment collaborates with Lehigh University, the University of Florida, and the Institute of Astronomy and Geophysics, Mongolian Academy of Sciences (Meltzer et al., 2019), and complements the MOBAL transect (Figure 3a).
The stations recorded teleseismic earthquakes and, by using traveltime tomography as well as P-to-S RF method, information on the lithospheric structures in Central Asia was obtained (Feng, 2021;Mordvinova et al., 2007;Petit et al., 2008;Tiberi et al., 2008;Zhao et al., 2021).The previous studies focused on the first-order signal and were able to characterize the main interfaces (Moho, LAB) beneath profiles or wider regions.Here, we mainly used the radial component of the RFs, while investigating the potential dip interfaces of the crust through the pattern of the transverse component.
We focused in priority on the relaminated arc, that is, the southern tip of MOBAL profile (TSET to ALTA, Figure 1) and the 44 easternmost stations of CMSE.Using those stations allows the analysis of the along-strike variations of structures in the lower crust (Figure 3a).

Gravity Data Set
The 2,847 ground gravity station data are provided by MonMap Company (convention between Czech Geological Survey and MonMap engineering company).The 10-25 km spaced gravity stations were acquired in the frame of Russian-Mongolian collaborative work in the 90s, in addition to the N-S oriented profile where the stations are 2 km spaced.The rest of the gravity data and specifically the densely measured profiles with 2-2.5 km spaced gravity stations were acquired by MonMap Company.Here, we present the map of the observed gravity data, leveled and corrected from the tide (Figure 3b).The gravity data need to be corrected from the topography effect to investigate the gravity anomalies at the crustal level.

Processing of the Receiver Function
RFs are time series that show the relative response of Earth's structure beneath the seismic station.They are computed from three-component seismograms, rotated into ZRT or LQT system coordinates (Kind et al., 2012;Langston, 1977).When a P wave encounters a velocity discontinuity at depth, part of its energy will convert into an S-wave (hereafter called Ps), mainly recorded on the radial R (ZRT) or Q (LQT) component, and which polarity reflects the sign of the velocity contrast: positive (negative) if the velocity increases (decreases) with depth (Langston, 1979).If the interface is horizontal and isotropic, the Ps phase will only be present on the radial component (R or Q).In case of anisotropy or dipping interface, however, the P-wave will also produce energy on the transverse (T) component (e.g., Jones & Phinney, 1998).Anisotropy and dipping interfaces result in timing and amplitude variations which are sinusoidal in the backazimuth for radial and transverse components.Both can also yield reverse polarity, where dips exceed the incidence angle, for example, examining the transverse component and its backazimuth variation is then of prime interest to identify potential dipping or anisotropic pattern in the crustal or subcrustal layers.
Moho depth for this area has been first estimated by Mordvinova et al. (2007) through different Vs depth profiles.We computed the RF for the seven southernmost stations of the French Mongolian Baikal Lithosphere seismic experiment (MOBAL2003) combined with 44 stations from the Central Mongolia Seismic Experiment in the study area (Figure 3a).We selected teleseismic events (Mb > 6.0) with epicentral distance ranging between 30° and 90° and filtered the original data between 0.08 and 0.8 Hz.The three components were rotated into the LQT system to ensure the maximum energy on the radial (Q) component, and we obtained the Q and the T RF using an iterative time-domain deconvolution with a Gaussian filter of 3.0 (Ligorría & Ammon, 1999).We select data with more than 80% of signal recovery, which eventually leads to minimum of 22 and maximum of 109 correct RFs per station (Table 1).

Determination of the Moho Depth
We proceed with the data set with the stacking method of Zhu and Kanamori (2000) to obtain the Moho depth and a Vp/Vs ratio estimate beneath each of the 48 sites.This inversion allows a combined estimate of crustal depth (H) and Vp/Vs ratio using a grid search on the radial component.It takes advantage of the Ps-converted wave and its multiples to reduce the inherent trade-off between H and Vp/Vs ratio.We used an a priori crustal P-wave velocity of 6.4 km.s −1 similar to the value of Mordvinova et al. (2007) and close to the one in He et al. (2016).The converted phase weighting is set to 0.6, 0.3, and 0.1 for Ps, PpPs, and PpSs + PsPs, respectively (Zhu & Kanamori, 2000).We estimate the deviation of our results relative to the data quality and the a priori P-wave crustal velocity with a bootstrap algorithm technique (Tiberi et al., 2005).We first calculate standard error from the inversion of 200 random subsets of data to estimate the variability of the results within the data set for a given crustal velocity (Vp = 6.4 km.s −1 , two last columns in Table 1).This helps to assess the data set coherency and to  (Mordvinova et al., 2007) and IRIS-PASSCAL (Meltzer et al., 2019) experiments superimposed on the topography map from Earth2014 (Hirt & Rexer, 2015).(b) Location of the gravity stations and the related raw gravity anomaly map leveled and corrected from the tide effect.sort out azimuthal variation or noisy data as sources of disparity.A second source of discrepancy within the inversion is tested through 200 inversions with crustal P-wave velocities randomly picked between 6.0 and 6.8 km.s −1 .
In Table 1 (uncertainty from Vp column), we report the dependency of the result on the ad hoc velocity model.These two bootstrap analyses show that uncertainty in our results mainly comes from crustal velocity influence (Table 1).In our case, a 5% decrease (respectively increase) of velocity will reduce (respectively increase) the crustal thickness by about 6%-8%.The Vp/Vs ratio will only change by ∼1%.The data set accounts for about 2% in average on the H uncertainty, with a maximum of 9% for HD68.
In most cases, the crustal thickness is retrieved within 2 km of uncertainty, except for HD25 and HD68 station, which presents a higher standard deviation (>3 km) due to the signal azimuthal variability.For these two stations, Moho multiples are hardly detectable, probably due to interference with either lithospheric interface or intracrustal multiples.
For most of the stations, the Moho signal is clearly associated with a positive and sharp velocity gradient (at ∼6-7 s), and its first and second multiples are identified between ∼19 and 21 s and ∼24-26 s, respectively (Figure 4; Figure S1 in Supporting Information S1).For some stations, in particular HD09, Ps converted wave at the Moho exhibits a weaker amplitude that may be flattened by interference with strong intracrustal conversions.From the stacking method, the Moho depth varies from 42.5 to 52.8 km (Figure 5) for an average value of 47.2 km.Vp/Vs ratio ranges from 1.70 to 1.82 with a mean value of 1.75.Those values are coherent with He et al. (2016) for their westernmost area (near 104° longitude).
Except for OVGO and TUSG, the relative Moho variations along the profile are similar to Mordvinova et al. (2007) study.Total crustal thickness variation deduced from the stacking method along the profile (6 km) is less than the 9 km identified by Mordvinova et al. (2007).As no error estimate is available from their study, we estimated their crustal thickness directly from their Vs velocity profiles and we cannot rigorously compare.The differences can also be due to the inversion techniques themselves, which differ.In our case, we consider a common and constant crustal P-wave velocity of 6.4 km.s −1 , whereas they consider a complex multi-layer crustal structure.Thus, we do not account for lateral geological heterogeneities, which can explain the discrepancy for the northern stations of the profile.The highest difference with Mordvinova et al. (2007) study is found for TUSG and OVGO stations, where an LVZ is proposed at ∼30-35 km depth.However, those two stations are reprocessed and complemented in our study with HD11 and HD15 from CMSE, which show remarkable coherency in the seismic signal.The uncertainty drops to ± 3 km for these two sites (Table 1).
A depth-migrated P-wave RF profile is presented (Figure 6).The P RFs are back-projected along the ray paths using IASP91 velocity model, and stacked within 3 × 3 km 2 bins (Kind et al., 2012).We identify continuous P-to-S conversion at the Moho, and several other crustal interfaces from South to North through the different tectonic zones.The migrated profile and the normalization of its Ps amplitude reveal two main trends of the seismic signal with a boundary located in the middle of the Baydrag microcontinent.The MOO and the northern part of Baydrag crust (HD08-HD07) have a rather homogenous seismic signal, whereas the southern Baydrag, the LZ, and the MAAW display a contrasting and alternating seismic signal (HD09 and HD69; Figure 6).
The Moho clearly deepens beneath the Hangay dome (north of HD06 and South Hangay fault, >50 km), following the topographic variation and coherent with previous seismic analyses (Feng, 2021).The crust thins southwards to reach ∼45 km with a Moho signal blurred and unfocused in the Baydrag zone beneath HD09 station.This apparent thickening of the crust-mantle boundary is certainly coming from a backazimuth variation effect, which is clearer for the non-smoothed version of the migrated profile (Figure 6), and which also appears in the  data (Figure 4b).In addition, HD09 displays contrasted and high amplitude intracrustal reflections, which are also retrieved beneath HD08 and ALTA in the MAAW.Those interfaces are associated with strong negative signals that we relate to negative velocity contrasts within the crust (low-velocity zones).

Intracrustal Interfaces Determination
Aside from a clear Moho signal, some of the analyzed stations exhibit a distinct and coherent negative polarity signal at ∼2-3 s (Figure 4).This is particularly the case for HD09, where this intracrustal signal is stronger in amplitude than the Ps Moho conversion.These negative peaks are often associated with transverse azimuthal variations (Figure S2 in Supporting Information S1), and may reflect a decreasing velocity with depth associated with a dip and/or anisotropic material.This behavior is particularly clear for the southern part of the migrated profile (MAAW, Lake, and Baydrag zones, Figure 6), with an enhanced amplitude below HD09 (Baydrag).
We investigate the characteristics of these interfaces and potential LVZ through a series of inversions following the method of Frederiksen et al. (2003).The inversion is stochastic and offers the possibility to retrieve both dipping layers (thickness, Vp, Vs, dip, and strike) and their anisotropy behavior (%Vp, %Vs, trend, and plunge).However, the increase of model space parameters will result in a less constrained solution.We decided to look only for layer thickness, Vp, Vs, dip, and strike within the crust.Even if methods allow to investigate the anisotropic behavior of crustal part (Feng, 2021;Feng & Diaz, 2023;Shen et al., 2013), the anisotropy is too complex in our multilayer model to be constrained.Besides, we are looking for intracrustal information we can combine with gravity signal, and anisotropy will be of no help in our case.
We apply the inversion on HD09 station, which presents the most coherent radial component and which overlays potential relaminated material (Guy et al., 2015).Among all our tests, the best results are obtained for Even if the transverse component displays a coherent pattern (Figure S2 in Supporting Information S1), the number of layers prevents from having a simple well-identified solution due to the large degree of freedom and to the complex interferences between the conversions.From our preferred models, all dips are less than 15° and more likely South or Southwest.

Gravity Analysis
The gravity anomalies are used to investigate the density contrasts and decipher the crustal structures, particularly in lower crustal parts.Before correlating the gravity signal with the different tectonic units and performing the density modeling, the measured data has to be corrected from the effect of topography to get the complete Bouguer gravity anomalies.Then, the geometry and density distribution of the crust can be investigated by building a 3D forward gravity model constrained by the geological and RF data.

Gravity Processing
The complete Bouguer gravity anomaly grid is obtained from the free-air anomalies by computing the terrain correction of each gravity station (Data Set S1).The terrain correction was performed using a mean crustal density of 2,670 kg/m 3 by the GEEC Matlab-based code (Cattin et al., 2015;Saraswati et al., 2019), which takes into account the sphericity of the Earth and the global effect of the Earth's gravity field.Three different digital elevation models (DEMs) were used to test the sensitivity of the results according to the resolutions of the DEM and the gravity data: (a) Earth 2014 global topography model (Hirt & Rexer, 2015) of 1 min arc resolution (∼2 km); (b) Shuttle Radar Topography Mission 30 (SRTM) of 30 s arc resolution (∼1 km); and (c) SRTM1 (Farr et al., 2007) of 1 s arc resolution (∼30 m).
The data sets and the complete Bouguer gravity results were extracted along the N-S profile, where most of the seismic stations are located, to examine the similarities and the discrepancies between the topography models, the results of the complete Bouguer anomalies, and the different seismic Moho depth (Figure 7).The three topography models are equivalent from the southern part of the Mongol-Okhotsk domain to the MAAW but display significant discrepancies from 0 to 165 m in the rest of the Mongol-Okhotsk domain, which corresponds to the high-frequency topography relief of the Hangay Dome.Thus, the Earth2014 model has a lower amplitude in this area than the SRTM30 and the SRTM1.The raw gravity data are logically anti-correlated with the topography.The differences in the frequency of the topography between the topography models are in return impacting the terrain correction results, as the complete Bouguer anomalies obtained with the Earth2014 model have lower amplitude than those obtained with SRTM30 and SRTM1.The comparison of the resulting complete Bouguer anomalies also shows that the use of SRTM1 does not necessarily bring strong differences in the terrain correction taking into account the resolution of the gravity data.Thus, for future computations, using SRTM30 may be enough and for sure a gain of time.In the last section, we display and discuss the differences between the Moho depths obtained by GOCE satellite gravity gradient modeling, the previous seismic RF analyses and our seismic results (Figure 7).The discrepancies observed between the seismic Moho and the gravity Moho depth from GOCE gravity gradient data are quite significant.In addition, the different seismic studies show rather poor agreement, particularly for the seismic stations OVGO and TUSG, which are the seismic stations located in the  Hangay Dome.These differences are probably due to the methods and assumptions used to get the Moho depth, as explained in the seismic section.
The resulting complete Bouguer gravity anomaly map with a grid spacing of 10 km ranges from −310 to −185 mGal (Figure 8).The Trans-Altai fault displays a mild gravity gradient, whereas the South Hangay fault correlates with a strong gravity gradient.The Bogd fault also displays gravity contrasts, which are not really linear but in the shape of "boudins."Finally, the gravity gradient located at the northeastern part of the map corresponds to a fault, which presents significant seismic events from 0 to 15 km depth (Meltzer et al., 2019) and identified as an Early Mesozoic rift system (Yarmolyuk & Kuzmin, 2012).

Constraints and Implementation of the Model
The gravity modeling of the central Mongolian crust is constrained by the depths of the Moho obtained by the RF analyses of this study, complemented by the topography of the Moho obtained by the 3D forward modeling of the GOCE gravity gradients (Guy et al., 2017).In addition, the modeling is constrained by geological data such as the surface boundaries of the different tectonic zones and their characteristic lithological composition (e.g., metamorphic, magmatic, or sedimentary rocks) (Figure 1), the dip of the contacts determined by previous studies (Badarch et al., 2002;Buchan et al., 2001;Cunningham et al., 1996;Sukhbaatar et al., 2022;Walker et al., 2007;Zorin et al., 1993).Thus, these benchmarks bring constraints on the shape and the depth of the main crustal density units.The RF inversion brought only poor constraints on the dip of the contact, and the signal's frequency content does not permit high detail level in the structure.We thus only took first-order information from these inversions: average depth location of low-velocity zones and interfaces.
The model also benefits from density values measured on 343 rock samples collected in the different tectonic zones (Figure 9; Table S1).The density values range from 2,430 to 3,312 kg/m 3 and the mean rock density of 2,739 kg/m 3 is obtained with a standard deviation of 145 kg/m 3 .However, the Gaussian of the frequency diagram gives a pick at 2,680 kg/m 3 .The highest density average distribution correlates with the location of the Precambrian metaigneous garnet-bearing amphibolites in northern Baydrag and the lowest density average distribution corresponds with the Devonian-Carboniferous turbiditic basin in the Mongol-Okhotsk domain.The intermediate density values from 2,680 to 2,750 kg/m 3 are found in the MAAW and the TAZ containing volcano-sedimentary sequences related to their oceanic environments.However, the constraints on the density values remain poor for a given lithology, surely metamorphosed in the deeper part of the crust and thus, the seismic velocities and the geological model are integrated to assess the most suitable density.Table 3 summarizes the density values that were assigned in our 3D forward gravity model.
The 3D forward gravity modeling was performed with the IGMAS+ software package (Anikiev et al., 2020;Gotze & Lahmeyer, 1988;Schmidt et al., 2011).The 3D block model is built by 32 north-south oriented 2D parallel sections with spacing of 15 km, which contains the density units linked by polyhedrons with triangular surfaces connecting the vertices between the successive sections.The model covers an area of ∼480 km × 400 km and reaches 70 km in depth.First, as the topography of the Moho is constrained by the seismic RF results and the previous GOCE gravity gradient model (Guy et al., 2017), a simple homogeneous crust-mantle model with standard density values was generated in order to observe the contribution of the variations of the Moho within the long-wavelengths of the gravity signal (Figure 10a).Thus, we can see that the Moho density contrast does not explain the entire gravity signal.Second, the crust was divided into six tectonic zones, each split into the conventional upper, middle, and lower crusts (Christensen & Mooney, 1995), with densities increasing with depth (Rudnick & Fountain, 1995).During this process, the dips of the main contacts between the distinct tectonic zones were maintained and the average density values were attributed according to the characteristic surface lithologies observed for each tectonic zone.At this step, the gravity signal is still not correctly modeled (Figure 10b).Finally, the model gained complexity by including the geological models (Buchan et al., 2001;Guy et al., 2015Guy et al., , 2021;;Zorin, 1999;Zorin et al., 1993) and the interfaces obtained by RF analyses.As the target is to model the structures of the crust, we did not take into account the possible variation of densities in the upper mantle and we assumed a laterally homogeneous upper mantle of 3,200 kg/m 3 .

3D Forward Gravity Modeling Results
The residual map of the final 3D model (Figure 11a) shows that the highest amplitudes of the residuals are at the western edge of the model.This border effect can be explained by the changes of lithologies in the western periphery of the model from Mongol-Okhotsk oceanic domain to the Zavkhan Precambrian microcontinent, whereas the geological units continue to the eastern periphery of the model (Figure 11b).The average misfits between the measured and the calculated anomalies have a correlation coefficient of 09.975 and a standard deviation of 2.35 mGal.
The N-S trending profile extracted from the 3D model is 530 km long (Figure 11c).Its location has been defined as it perpendicularly crosses the main tectonic zones and it runs along the area where most of the seismic stations from the MOBAL2003 and IRIS-PASSCAL experiments are located.The density structure profile starts at the southern margin of the Tarvagatay microcontinent, crosses the Mongol-Okhotsk oceanic domain, the Bayankhangor Complex, the Baydrag microcontinent, the LZ, the MAAW, and ends in the northern part of the TAZ.The complete Bouguer gravity anomalies range from ∼−215 mGal in the south to −290 mGal in the center of the Mongol-Okhotsk domain.The long wavelength gravity signal along the profile can be divided into two main trends: a high amplitude and low-frequency gravity over the MOO and the northern Baydrag, and a low amplitude and medium frequency gravity anomalies over the southern Baydrag, the LZ, the MAAW, and the TAZ (Figure 11c).The Moho depth, constrained by the RF analysis and the GOCE gravity gradient modeling (Guy et al., 2017), varies from ∼−45 to 50 km depth along the profile and can be considered as rather flat in its southern part as it varies from −46.6 to −48 km.The upper crust of the southern Paleozoic tectonic zones, namely the LZ, the MAAW, and the TAZ, is modeled, taking into account the compositions of the different lithologies, as laterally variable nearly horizontal density bodies, which reproduce the large-scale folding and thrusting of the accretionary system cut off by steep contacts between each tectonic zone (Guy, Schulmann, Clauer, et al., 2014;Lehmann et al., 2010;Sukhbaatar et al., 2022).These features correspond to the short-wavelength anomalies observed in the southern part of the profile.The nearly constant long-wavelength anomaly is interpreted as the Table 3 Densities of Rocks Used in the Model lack of density variation in the lower crust, since the Moho depth does not significantly oscillate.Thus, the lower crust underneath the LZ, the MAAW, and the TAZ has a homogeneous gravity signature and is characterized by a felsic density value of 2,730 kg/m 3 (Guy et al., 2015(Guy et al., , 2021)).The Bogd fault located in the LZ can be identified by a small gradient in the gravity signal and is modeled accordingly.The gravity analysis and modeling point toward a slight crustal discontinuity between the Baydrag microcontinent and the LZ as exemplified by the short wavelength gravity low at the contact, which can however correspond to the presence of Ordovician metasediments in the Proterozoic basement of the LZ or to the presence of the Ikh-Mongol arc relics east and west of the profile (Figure 11c).Geological constraints concerning the Baydrag microcontinent are rather poor so far.We can just assert that there is a significant difference between the southern, seismically heterogeneous, and the northern part, seismically homogeneous (Figure 6).Thus, the Baydrag microcontinent is modeled as a multi-layered basement crust in its southern part, where two low-density layers alternate with higher-density layers according to the resulting seismic interfaces, and a rather homogeneous to two-layer crust with densities of 2,700 and 2,850 kg/m 3 in its northern part, with an undetermined steep boundary in the middle of the microcontinent.The Bayankhongor Complex is 20 km wide in the N-S direction and the gravity data have 10 km spaced measurements in this area.Thus, it has a poor resolution and could not have been modeled indeed.Nevertheless, according to geological studies (Buchan et al., 2001;Osozawa et al., 2008) and to maintain a tectonic coherency, the Bayankhangor Complex is modeled as three vertical density units reproducing the ophiolitic mélange intercalated in volcano-sedimentary to sedimentary clastic sequences.Finally, the Mongol-Okhotsk domain is modeled as a two-layer sedimentary crust, which represents the Silurian sedimentary sequences overlain by the Devonian-Carboniferous turbiditic basins together folded into a synclinorium pattern, later intruded by Permian-Triassic granitoids.The Cenozoic sedimentary basins in central Mongolia are supposed to be approximately 0-2 km thick and short-wavelength anomalies do not spatially correlate with them at the location and scale of our 3D modeling area.Therefore, they were not included in the model.

Discussion
The following section examines the seismic and gravity patterns of the crust in Central Mongolia by first discussing the similarities and differences of our crustal model with the electromagnetic profiles acquired in the same area and the significance of these geophysical anomalies.Then, we analyze the distribution of the seismic interfaces from the RFs and the geometry of the density structures from the forward modeling, related to surface geology of the different tectonic zones.These correlations enhance new crustal structures, such as middle and lower crust low-velocity/low-density zones (LVLDZs) in the Precambrian block and the Paleozoic domains, and challenge some boundaries between the tectonic zones.Thus, we demonstrate the complementarity of seismic and gravity methods to decipher the crustal structures in complex, reworked, and long-lasting tectonic settings.

Comparison of Seismic and Gravity Results With Electromagnetic Profile
To evaluate our crustal model based on seismic and gravity data sets, we compare it with recent N-S-oriented electromagnetic profiles (Comeau et al., 2018(Comeau et al., , 2020)), whose data were acquired nearly at the same location where most of the seismic stations from the MOBAL2003 and IRIS-PASSCAL experiments are distributed (Figure 12).In general, the electrical resistivity models of the crust reveal a high resistive upper crust (0-20 km) with several conductive features, that are interpreted as fault or suture zones, and a low resistive lower crust, interpreted as a presence of fluids and a weakened lower crust.
Specifically, the Mongol-Okhotsk upper crust displays a variable resistivity (10-2,000 Ωm) in its first 0.5 km (C1 (2018) , Figure 12b) probably due to the occurrence of porous sediments.Below, the upper crust presents rather highly resistive zones (∼10,000 Ωm) interspersed with vertical channels (C4 ( 2018) ) of lower resistivity (400-1,500 Ωm), which correlate with location of volcanism and present-day hydrothermal activity.The MOO lower crust has a low resistivity zone between 30 and 50 km (20-80 Ωm), indicating the presence of fluids and a weakened lower crust (C2 ( 2018) ).These two electrical layers of the MOO, each rather homogeneous, concur with the homogeneous seismic signal of the MOO crust (Figure 6) and the two density layers of the 3D forward model (Figure 12a).The electrical resistivity profile reveals a major crustal boundary in the shape of a low resistivity zone (20-40 Ωm) concerning the South Hangay fault system and the ophiolitic belt (F1 (2018) ) interpreted as the marker of hydrothermal alteration along fossil fluid pathways due to the location of a suture zone (Comeau et al., 2021).Ophiolites have to provide distinct gravity and seismic signals due to their compositions (Christensen, 1978).The Bayankhangor ophiolitic complex is a NW-SE belt of cca.300 km long and cca.20 km wide presumably obducted within two accretionary systems, so that a strong seismic and gravity signal should be expected.However, the Bayankhongor complex does not exhibit strong seismic interfaces (HD6, HD7, and BUMB; Figure 6).Moreover, only the southern part of the complex correlates with a gravity high centered on the ophiolitic mélange; otherwise, no distinct gravity high can be distinguished along its entire length (Figure 8).This lack of seismic interface and gravity highs may be related to the spacing of gravity and seismic stations or because the ophiolites may be only superficial and not deeply rooted (Figure 12a).Similarly, the South Hangay fault is not associated with strong seismic and gravity signal anomalies, which may indicate a shallow depth of the contact and a low-density contrast between the juxtaposed lithologies, respectively.Comeau et al. (2018Comeau et al. ( , 2020) ) interpreted the shallow conductive features (<5 km) in the southern Baydrag microcontinent as a significant mineralization and past hydrothermal fluid alteration leading to the Tsagaan Tsahir Uul gold deposit.Their electrical model also shows a low resistivity zone starting from 20 km in the northern and 10 km in the southern parts of the Baydrag microcontinent with a vertical heterogeneous resistivity up to 20 km depth in between (G1 (2020) ).Overall, the northern part of the Baydrag microcontinent appears more resistant than the southern part.This trend is similar to our seismic and gravity results showing that the Baydrag microcontinent is geophysically divided into two parts: the northern part does not display strong seismic interfaces, whereas the southern part presents an alternation of low-velocity with high-velocity zones respectively corresponding to low-density with high-density layers.The Bogd fault in the LZ correlates with a low resistivity anomaly (30-100 Ωm) dipping to the south (F2 ( 2018) ).However, the apparent dip angle of the fault given by the electrical model is ∼30° (Comeau et al., 2018), whereas the gravity modeling results in a subvertical dipping fault (Figure 12a).Moreover, no seismic contrast enables the detection of the Bogd fault, although it was identified as a deeply rooted (up to 20 km depth) active fault (Kurushin et al., 1998;Rizza et al., 2011).This poor resolution of the fault may result from the seismic coverage and the orientation of the fault, which might be too vertical to generate an RF signal.Finally, the distinct zones of low resistivity material (50-100 Ωm) located at 30-45 km depth beneath the LZ and the MAAW (C1 ( 2020) ) and the low resistivity zone (∼100 Ωm) at 15-30 km depth in the TAZ (C2 ( 2020) ) are in agreement with the identification of the LVLDZ in our seismic and gravity models (Brownish density body of 2,730 kg/m 3 , Figure 12a), although the RF results give a shallower low-velocity zone starting from 20 km up to the Moho in the southern part of the profile.Therefore, the electrical resistivity model and our model based on seismic and gravity data converge regarding the boundaries of the different tectonic zones and their upper and lower crustal structures.

Crustal Structures of the Mongol-Okhotsk Oceanic Domain
The Mongol-Okhotsk oceanic domain is the Paleozoic tectonic zone, also corresponding to the Cenozoic Hangay Dome, which is an intracontinental high-elevation low-relief landform in Central Mongolia.The Hangay Dome is surrounded by Mesozoic and Cenozoic intraplate volcanism attributed to the asthenospheric upwelling (Ancuta et al., 2018;Hunt et al., 2012;Sheldrick et al., 2018).Recent geochemical and geophysical studies were undertaken to decipher the mechanisms of such a dome formation and intraplate volcanism in intracontinental tectonic settings (e.g., Ancuta et al., 2018;Comeau et al., 2021;Tiberi et al., 2008).The thickness of the southern Mongolian crust is rather constant with the seismic and gravity Moho located at ∼45 km depth (Figures 6 and 11), which is an average continental crustal thickness considered above the estimated standard of around 39.5 km (Christensen & Mooney, 1995).This steady thickness of the crust in this part of the Mongolian collage was also demonstrated by previous seismic and gravity studies (Guy et al., 2017;Mordvinova et al., 2007;Petit et al., 2002).Then, the thickness of the crust increases in the southern part of the MOO beneath the Hangay Dome, whereas the asthenospheric upwelling reveals the thinning of the lithosphere leading to the Cenozoic magmatism (Guy et al., 2017;Petit et al., 2002Petit et al., , 2008)).However, the flexure of the Moho is not enough to explain the high amplitude gravity low over the MOO (Figure 10a).
Interpretations of geochemical and gravity data led some authors assuming that the Devonian-Carboniferous sedimentary basin of the MOO developed on a Precambrian basement, called Hangay microcontinent (Cunningham, 2001;Osozawa et al., 2008;Purevjav & Roser, 2012;Zorin et al., 1993).Taking into account the homogeneity of the seismic signal and the strong gravity low, these geophysical characteristics in the MOO can be explained by the combination of the Devonian-Carboniferous sedimentary layers on the Silurian metasediments and their magmatic reworking in Permian-Triassic with the intrusion of mantle and crustal source granitoids (Yarmolyuk et al., 2016).Although at the scale of the 3D model, no real close correlation between the location of the Permian-Triassic granitoids in the MOO and the short wavelength gravity low is generally observed, we integrated these granitoids to the model to avoid any possible bias.However, the involvement of a hypothetical Precambrian basement will not help to model the gravity low and this sedimentary layer structure is in agreement with the absence of a distinct seismic interface within the MOO, which could have been interpreted as the interface between the Paleozoic sedimentary basin and a Precambrian continental fragment.Therefore, our model may reveal the absence of a Precambrian basement in the MOO (Figure 13

Significance of the Low-Velocity and Low-Density Zones in the Baydrag Microcontinent
As previously mentioned, contrasting seismic signals are observed between the northern and the southern parts of the Precambrian Baydrag microcontinent (Figure 6).Its northern seismic structure does not show strong interfaces, which implies either a relative homogeneity of the lithologies or highly verticalized structures due to a high convergence rate.Hence, this was modeled as a homogeneous density body with a thin denser lower crust.Its southern seismic structure is characterized by sub-horizontal interfaces showing the alternance of low-velocity and high-velocity zones, which were modeled as a tectonic layering with the alternance of low-density and high-density bodies accordingly (Figure 11).
Geological investigations leading to a geodynamic reconstruction of the Baydrag microcontinent are rather scarce, particularly in its southern part due to the occurrence of the Cenozoic sedimentary basin, which conceals the Proterozoic and Paleozoic tectonic units.Thus, the tectonic interpretation of the Baydrag southern margin mostly comes from an assumed continuation of the tectonic trends from the Zavkhan southern margin (Soejono et al., 2023), as it belongs to the Precambrian microcontinent belts located in the center of the Mongolian collage.Nevertheless, studies converge with the probability of contrasting geodynamics affecting the northern and southern Baydrag microcontinent.The northern part of the Baydrag microcontinent is affected by an early Tonian (890-855 Ma) extension of the continental crust revealed by HT-LP metamorphism and interpreted as the consequence of the subduction roll-back of the Mirovoi Oceanic plate beneath the Rodinia (Štípská et al., 2023).This event is followed by the late Tonian thickening of a back-arc or arc domain due to the subduction advancing system.The southern part comprises a succession of a Tonian magmatic arc related to the subduction of the Mirovoi Ocean beneath Rodinia, followed by a Cambro-Ordovician continental magmatic arc related to the subduction of the Paleo-Pacific Ocean beneath Gondwana (Buriánek et al., 2017;Soejono et al., 2023).This Cambro-Ordovician magmatic event is not observed in the northern part.In addition, the Tonian magmatism related to an arc was not described in the northern part of Baydrag microcontinent and if it is, it is rather subordinate and will probably not be distinguished at crustal scale.
Therefore, based on the geological studies, the difference of crustal structures between the northern and southern parts of Baydrag revealed by seismic and gravity data analyses may be explained by the distinct tectonic events affecting the Baydrag microcontinent.In the southern part, the LVLDZs may constitute the remnants of the batholiths resulting from the early Tonian or the Cambro-Ordovician continental magmatic arc (Figure 13) or may result from late Tonian anatexis processes due to the extension during the break-up of Rodinia that would have affected southern Baydrag.However, this late Tonian extension related to the break-up of Rodinia was identified in all the other Precambrian microcontinents in the Mongolian collage except in Baydrag (Soejono et al., 2023).Thus, we rather assume the hypothesis of remnants of batholiths, which can partly trigger the alternation of high velocity-high density layers with low-velocity/low-density layers at the crustal scale, as it was demonstrated in the Archean Karelia and Kola provinces, where the acoustic transparent regions and low seismic reflection zone were interpreted as fragments of batholiths (Kostyuchenko et al., 2006;Mints et al., 2020).

Relaminated Material Under the Mongol Altai Accretionary Wedge and the Trans-Altai Zone
It is accepted that the present continental crust is formed by magmatic processes at volcanic arcs above subduction zones (e.g., Hacker et al., 2015;Hawkesworth & Kemp, 2006).Volcanic arcs have a geochemical composition very close to the continental crust, at least in its upper part, since the deep crust of the arcs (from ∼20 km below the surface) has a very different geochemical composition from the continental crust at equivalent depth.
The denser part of the arc can undergo delamination, which is a loss of part of the denser lower lithosphere that sinks into the mantle and is recycled (e.g., Bird, 1979;Magni et al., 2013).However, this process is not sufficient to transform a lower arc crust into a continental crust.Delamination does occur, but for it to be the only mechanism responsible for crustal differentiation, a complex process of crustal thickening and repeated metamorphic events would be required.Thus, a new process of crustal reworking was proposed (Hacker et al., 2011).During subduction followed by continental collision, a part of the crustal lower-plate material becomes buoyant and can be redistributed within the crustal upper-plate inducing its compositional and structural transformations.The process of such transformation of the continental crust is called relamination (Hacker et al., 2011;Kelemen & Behn, 2016) and the resulting material of redistributed meta-sedimentary and meta-igneous rocks is referred to as the relaminant.Recent petrological data (Hacker et al., 2011(Hacker et al., , 2015)), geochemical data (Guy et al., 2015; Kelemen The seismic interfaces observed in the LZ and the MAAW crusts show an alternation of increase and decrease of velocities and can be interpreted as a succession of high-velocity and low-velocity zones together with a low Vp/Vs ratio (Figures 4 and 6).The geological cross-sections proposed in these tectonic zones provide evidence of the succession of folding and thrusting of the different lithologies stemming from the magmatic and tectonic events, and sedimentary deposits affecting the LZ Cambrian passive margin and the Mongol-Altai Ordovician Accretionary Wedge (Buriánek et al., 2017;Sukhbaatar et al., 2022;Štípská et al., 2010).These folds and thrusts were roughly reproduced in the gravity model in the shapes of folded sedimentary layers in the MAAW and thrusts of metasediments on the Proterozoic basement in the LZ, which enables it to fit the observed with the calculated gravity signal (Figure 11).The geologically constrained upper parts of the LZ, MAAW, and TAZ allow us concluding that the weak amplitude of the gravity signal (∼25 mGal) over these zones indicates a lack of deep crustal discontinuities (Guy, Schulmann, Munschy, et al., 2014, 2015), which is modeled as a homogeneous felsic lower crust (Guy et al., 2015(Guy et al., , 2021)).In the lower crustal part of the LZ and MAAW, the LVLDZ is interpreted as the occurrence of an allochthonous hydrous or felsic material (Figure 13), which explains the low Vp/Vs ratio (1.71 below ALTA).Our 3D geophysical model also provides important constraints for the possible extent of the relaminant at the scale of southern Mongolia, and consequently, may indicate a huge amount of material being relaminated.

Conclusions
The combination of RF with gravity data analyses and modeling, constrained by geological data, not only refines the topography of the Moho in central Mongolia, but also provides substantial insights into the crustal structures of central Mongolia and its tectonic evolution.
• The significant difference of the crustal seismic and gravity anomalies between the Mongol Okhotsk Oceanic domain and the southern tectonic zones, namely the Baydrag microcontinent, the LZ, the MAAW, and the TAZ, involves different crustal structures and testifies to the contrasting geodynamics.• The seismic homogeneity of the Mongol Okhotsk Oceanic crust and the 3D gravity model demonstrate the absence of Precambrian microcontinent nuclei beneath the Hangay dome.• The contrasting seismic signatures between the northern and southern parts of the Baydrag microcontinent confirm the distinct geodynamic evolution between its northern and southern margins.Moreover, the LVLDZs within the Baydrag Precambrian block middle crust can be interpreted as the remnants of the Grenvillean or Cambro-Ordovician batholiths.
• The LVLDZ at the lower crustal level beneath the LZ, the MAAW, and the TAZ reveal the occurrence of a felsic material, and, thus the existence of an allochthonous crust of continental affinities beneath an oceanic crust.This is interpreted as the result of a relamination, which is probably an underestimated tectonic process of continental crust formation.

Data Availability Statement
Seismic data (Meltzer et al., 2019) used in this study are open-access and obtained from https://www.iris.edu/hq/.The density data used in the 3D model are available as supplementary material.Satellite gravity data (Guy et al., 2017) are available from web services at https://earth.esa.int/eogateway/missions/goce.Land gravity data supporting this research are available at the MonMap engineering company and were obtained upon a formal convention between Czech Geological Survey and MonMap engineering company.The Complete Bouguer anomaly grid used in the modeling is available in supplementary material.The grid was calculated using the GEEC software (Saraswati et al., 2019; https://github.com/RodolpheCattin/GEEC). The 3D forward gravity modeling was performed using the IGMAS+ software (Anikiev et al., 2020; https://igmas.git-pages.gfz-potsdam.de/igmaspages/).The other data can be obtained in the tables and the references.

Figure 1 .
Figure 1.Tectonic and geological maps (modified from Parfenov et al. (2003) and Wang et al. (2019)).(a) Global tectonic map of the Mongolian collage with the location of the 3D forward modeling area highlighted by the red box.The black curves are the state borders.The dark red area underlines the assumed extent of the Ikh-Mongol Arc remnants.(b) Geological map of Central Mongolia (modified from Parfenov et al. (2003) and Guy, Schulmann, Munschy, et al. (2014)) with the superimposed tectonic zones.

Figure 2 .
Figure 2. Spatial coverage of the previous seismic and gravity studies superimposed on the gray-shaded topography.

Figure 3 .
Figure 3. Central Mongolia seismic and gravity data sets.(a) Location of the seismic stations from MOBAL 2003(Mordvinova et al., 2007) and IRIS-PASSCAL(Meltzer et al., 2019) experiments superimposed on the topography map from Earth2014(Hirt & Rexer, 2015).(b) Location of the gravity stations and the related raw gravity anomaly map leveled and corrected from the tide effect.

Figure 4 .
Figure 4. Radial receiver functions for four stations associated with pivotal tectonic zones (a-Mongol-Altaï, b-Baydrag, and c-Hangaï).The origin time (t = 0 s) refers to the direct P-wave arrival time.The receiver functions are sorted by backazimuth, and the P to S converted phase at the Moho (Ps) is surrounded by the black rectangle.The alternance of negative and positive pics between t = 0 and the Ps arrival time indicates crustal interlayering with potential low velocity layers.

Figure 6 .
Figure 6.Depth migrated profile (middle panel) and CCP profile (bottom panel) along the transect N12° centered in HD10 station.Ps amplitude has been normalized.Top panel represents the topography along the profile (thick black line), and the gray area corresponds to the envelope for topography 30 km apart from the transect.Stations (black triangles) are projected along the transect.

Figure 7 .
Figure 7.Comparison of the different data sets, their resulting processing, and previous studies.The N-S profile extracted in all data set.along with Mordvinovaet al. (2007).Profiles comparing the topography of three different topography models, the raw gravity data, the resulting complete Bouguer anomalies, and the seismic Moho depth from our study and from previous studies.

Figure 8 .
Figure 8. Complete Bouguer anomaly map with the lithological contacts and the tectonic zones, the main faults, and the location of seismic stations superimposed.

Figure 9 .
Figure 9. Density sample locations and average distributions according to the distinct tectonic zones used for the modeling.Values in kg/m 3 .

Figure 10 .
Figure 10.Modeling constraints and implementation of the model.The profile presented here is extracted from the 3D block along the N-S profile, where most of the seismic stations are located.The resulting residual maps of the corresponding 3D forward gravity modeling are shown for each implementation.Densities are in kg/m 3 .(a) Crust-mantle gravity model with standard density values.The resulting Moho depth from seismic receiver function and in between the Moho depth from GOCE satellite gravity gradient 3D modeling are depicted.(b) Gravity model with conventional three-layer crust and density values increasing with depth.The tectonic boundaries are displayed.

Figure 11 .
Figure 11.Results of the 3D gravity model.(a) Illustration of the final 3D block model.(b) Residual anomaly map of the model.(c) N-S profile extracted along the seismic stations.The characteristics of the density bodies are given in the legend.Densities are in kg/m 3 .

Figure 13 .
Figure 13.3D block diagram illustrating the results and interpretation.

Table 1
Results for the 48 Stations Stacking Process of the Receiver Functions H (km) Vp/Vs H (km) Vp/Vs

Table 2
Preferred Models From the Inversion of HD09 Receiver Functions