Integrated Geophysical Evidence for the Middle-Lower Crust Melting of the Songpan-Aba Terrain, NE Tibetan Plateau

The Songpan−Aba region is located on the northeastern edge of the Tibetan Plateau. Tectonically, the area is surrounded by the West Qinling orogenic belt in the north, the Longmenshan orogenic belt in the southeast, and the East Kunlun and Sanjiang orogenic belts in the west and southwest, forming a triangle that provides an ideal location to study the crust-mantle structure and deep tectonics of the eastward extrusion of the Tibetan Plateau. In this study, the magnetic and electrical structures of the Songpan−Aba area were investigated by inversion using high-precision magnetic anomaly and magnetotelluric data to obtain the subsurface magnetization inversion intensity and resistivity of Songpan–Aba and adjacent areas. The results revealed a continuous magnetic layer up to 20 km below Songpan–Aba and its surrounding areas in the south, possibly originating from a magma root southwest of the Longmenshan massif. In the West Qinling, Songpan–Aba, and Longmenshan areas, pervasive low-resistance, weakly magnetic, or magnetic layers were identified below 20 km that might be formed from the molten mantle material extruded from the eastern edge of the Tibetan Plateau.


INTRODUCTION
Songpan-Aba is located at the northeastern edge of the Tibetan Plateau ( Figure 1) and situated at the interface of the east-west and north−south tectonic components of China. It is an important region for studying the eastward extrusion of crustal material from the Tibetan Plateau (Tapponnier et al., 1982).
Since the 1980s, several studies have proposed deformation models to explain the eastern Tibetan Plateau uplift. The first was the rigid block model, which concluded that the uplift of the eastern Tibetan Plateau could be primarily attributed to the deformation of brittle substances in the crust caused by thrust nappe (Tapponnier and Molnar, 1976, 1982, 2001Avouac and Tapponnier, 1993); the second was the continuous deformation model (England and Dan, 1985, 2010Houseman and England, 1986); and the third was the crustal flow model, which proposed that a low-viscosity asthenosphere in the lower crust of the eastern Tibetan Plateau slid eastward under plate movement and was blocked by the Sichuan Basin, resulting in the accumulation and lift of a lower crustal substance with low velocity and high conductivity (Bird, 1991;Royden et al., 1997;Royden et al., 1977;Clark and Royden, 2000).
Seismic tomography studies (Wang et al., 2005;Wang, 2007;Wu et al., 2009), receiver functions (Liu et al., 2015), and environmental noise imaging (Yang et al., 2013) of the Tibetan Plateau concluded that crustal and upper mantle structures in the eastern margin of the Tibetan Plateau were different from those in its surrounding regions, such as the Sichuan Basin and the South China Plate. The seismic wave velocity indicated a relatively steep gradient change at the junction of the Tibetan Plateau, its surrounding plates, and middle and lower crustal flows in the northeastern Tibetan Plateau (Pan et al., 2017).
Magnetotelluric (MT) is an important method to probe the deep structure of the lithosphere and can effectively reveal the distribution of subsurface fluids. Several MT studies conducted in the eastern periphery of the Tibetan Plateau have elucidated the electrical structure of the regional crust and mantle, showing expansive low-resistance layers widely distributed in the middle and lower crust of the northeastern margin (Sun et al., 2003;Xu et al., 2007;Zhao et al., 2008;Wang et al., 2009Wang et al., , 2013Zhang et al., 2012;Bai et al., 2013), indicating the prevalent presence of fluids and local melting. Their distribution characteristics were consistent with the "crustal flow" model, providing strong evidence for it from the perspective of conductive structures ).
The regional magnetic anomalies were primarily due to highly magnetized crustal rocks that provided information on crustal   (Wang et al., 2020). Previous studies in this area did not adopt the combined magnetic and MT approach to investigate subsurface tectonics. Therefore, this study aimed to apply high-precision magnetic survey data to conduct magnetization intensity inversion in Songpan-Aba and FIGURE 2 | Magnetic anomalies in study area (A) and Magnetotelluric survey lines (B) (pink triangle is Line 1, blue square is Line 2, red triangle is Line 3 and green square is Line 4).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 752693 3 adjacent areas to obtain three-dimensional (3D) regional magnetic properties. Meanwhile, the subsurface electrical structure was derived from four high-precision MT profiles collected over a long period to analyze the deep structure of the Songpan-Aba area from subsurface magnetic and electrical property perspectives, and to further examine the deep geological structures of the eastern margin of the Tibetan Plateau.

3-D Inversion of Total-Field Magnetic Anomalies
The non-uniqueness and instability of the solutions have typically focused on previous magnetic anomaly inversion studies, as they are important factors restricting magnetic inversion development. Tikhonov and Arsenin (1977) proposed the regularization theory to solve such problems. The forward sum function, A, decreases rapidly with increasing depth, yielding an inversion result near the ground (Li and Oldenburg, 1998;Zhdanov, 2002). The depth weighting of this regularization method depends on the selection of the regularization factor, α. The adaptive regularization factor selection method proposed by Hu et al. (2019) was applied in this study. We followed Li and Oldenburg (1996), using the reciprocal of the vertical surface component of the total magnetic anomaly data as the dataweighting factor. If the model sensitivity matrix is directly incorporated into data-fitting (Zhdanov, 2002;Hu et al., 2019;Zhao et al., 2020), this problem can be solved. After model weighting, the total sensitivities of the weighted model parameter of all grid cells are equal, and the contribution of observed data is essentially identical. The objective function is defined as follows: where ϕ(m w ) is the data-fitting function and S (m w ) is the modelstabilizing function. The model weighted matrix (W m ) is: where m w apr is the weighted prior reference model, A w AW −1 m is the weighted forward operator, and m w W m m is the weighted model parameter.
If there is no remanence, and assuming that the inclination (I 0 ) and declination (D 0 ) of each data grid are known, then the oblique magnetization forward matrix A can be directly calculated. In this study, the inclination and declination were directly incorporated into the total magnetic intensity inversion process, ΔT. Zhdanov (2002) reported the iteration process of the conjugate gradient.

Magnetic Data
High-precision magnetic anomaly data were acquired in the study area using five G-856 and four G-858 magnetometers. We set up nine magnetic observation stations for diurnal variations in Jiangyou, Shuijingbao, Sedi, Lianghekou, Hongyuan, Jiuzhi, Maqu, Diebu, and Dangchang to adequately and effectively control the measurement area and improve the accuracy of each measurement point. The locations were free from human interference and active magneticity, with the convenience of use, clear markings, and accessible storage. Joint measurements were conducted at these stations from 5-8 am and 5-10 pm (when the magnetic field was relatively stable) to achieve magnetic field normalization. The data FIGURE 4 | Three-dimensional joint inversion of oblique magnetization. M 0.1 A/m iso-surface image (Z axis is the number of grids, and one grid represents 2 km). recorded during these two periods were used to obtain the basic magnetic field value of the observation stations. The radius for the joint measurements was controlled to within 50 km. The data collected included the total magnetic anomaly △T, with a total root-mean-square error of the magnetic measurements <5nT. The field survey was set up with 13 survey lines spaced 30 km apart. The measurement points were arranged on each line at 2km intervals, resulting in 1,985 points (see Table 1).
Data from other areas were collated and supplemented using Satellite Magnetic Anomaly Data (https://www.ngdc.noaa.gov/ geomag/emag2.html) (Figure 2A) from the Earth Magnetic Anomaly Grid (EMAG2) database published by the Commission for the Geological Map of the World (Maus et al. , 2009). This database of crustal magnetic anomalies was accumulated from airborne, offshore, and satellite magnetic surveys over 50 years at a 2' × 2' global grid resolution. As shown by the magnetic anomaly map of the Songpan-Aba area, the object under investigation in this study was characterized by stable weak magnetic anomalies, varying from −30 to −50 nT in a steady, smooth magnetic field. The magnetic field in the eastern part of the survey area varied drastically with local anomalies. The measurement area was divided into the following zones according to the strength of the magnetic anomalies: 1) North China plate: This region exhibited varying positive and negative magnetic anomalies. Numerous rock outcrops primarily formed in the Indosinian period, and a some formed during the Hercynian and Yanshanian periods, exist in this region. The distribution of exposed rocks was well correlated with the magnetic anomalies. The magnetic anomaly was in the north-northeast direction and positive, with the highest value >200 nT. The rocks formed during the Indosinian period exhibited various magnetic properties, ranging from weak to strongly magnetic, with their magneticity related to the magnetic minerals content of the rocks. The ferromagnetic mineral content is low in acidic rocks, which exhibit weak to moderate magneticity, whereas intermediate mafic rocks are more magnetic and show high magnetic anomalies. 2) West Qinling area: This area revealed broad, smooth negative magnetic anomalies stacked with local magnetic anomalies. The magnetic anomaly in this area was characterized as low, with some local high magnetic spots. The background value of the low negative magnetic anomaly was approximately −40 nT and the stacked local anomalies were higher overall than the background anomaly value. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 752693 in this area was north−east trending, with dense anomaly contours and rapid gradient changes between the negative anomalies. 4) The Songpan-Aba crust was a stable zone of low, negative magnetic anomalies of approximately −40 nT, with ∼10 nT fluctuations. The anomalies decreased to approximately −50 nT eastwards at Nanping.

Three-Dimensional Magnetic Structure
First, a theoretical model was designed to verify the practicality and stability of the proposed algorithm, then a model consisting of three blocks of different depths and sizes was designed, as shown in ( Figure 3A). The model space was 40 × 40 × 21 grid cells, that were each 1 × 1 × 1 km in size, and the magnetization of the model was 1 A/m. Figure 3B shows the synthetic data with the I 0 52°and D 0 −1°, where the data space was 40 × 40 grid cells, and each cell was 1 × 1 km. A profile selected in this study is indicated by the pink line in Figure 3B. Regularized inversion was then performed, and the results are shown in Figure 3C. The profile was selected for a detailed comparison of the model and inversion results, as indicated in Figure 3D. Further, the number of grids in the east-west, north-south, and vertical directions were 97, 81, and 40, respectively, for the 3D inversion of the measurement data obtained from the 5 × 5 km horizontal-dimension grids in the study area. The vertical depth of the forward modeling and inversion was on average, divided into 40 grids covering 40 km (the vertical interval for the grid is 1 km). One hundred iterations were performed in the inversion, in which the data-fitting error was reduced from an initial value of six root-mean-square (RMS) to approximately two RMS. Figure 4 shows the results of the regularized inversion based on the weighted model. We derived the trend and distribution of the magnetic anomaly in the study area from the iso-surface images.

MT-Regularized Inversion
A method for solving the non-uniqueness of geophysical inverse problems is described in the Tikhonov and Arsenin (1977) regularization theory that uses both data misfit and a modelstabilizing functional to construct the parametric objective functional. The objective functional of the regularized inversion can be written as where p (m,d) is the parametric objective functional, m is a vector consisting of model parameters, A(m) is the forward operator, d is a vector containing the observation data, ||·|| is the L2-norm, φ(m) is the model-stabilizing functional, and α is the regularization parameter. A stabilizing functional is used to select the appropriate class of models from the model space of all possible models that are a correction set for solving the inverse problem (Zhdanov, 2002). The minimum support (MS) stabilizing functional was proposed to obtain a model that minimizes the volume with anomalous model parameters to improve the resolution of the sharp electrical interface (Last and Kubik, 1983). It focuses on where/how the model parameters are different with respect to the given a priori model. Here, the MS stabilizing functional can be expressed as Xiang et al. (2017) and Zhang et al. (2020) analyzed the advantages and disadvantages of several commonly used families of stabilizing functionals, such as the smooth, MS or minimum gradient support (MSG). We then applied a new stabilizing functional, the MSG that combines smooth and sharp-boundary constraints with a unified form for the objective functional in the 2-D regularized inversion. The MSG stabilizing functional is defined as: The spatial gradient is calculated after the MS stabilizer is obtained using Eq. 3, thus, we refer to this new stabilizing functional as the MSG, which uses an MS value instead of m-m apr, leading to a stable constraint that focuses on sharp boundaries for inversion. This functional acts on the class of models such as the domain with an anomalous parameter distribution that occupies the minimum volume, while avoiding any instability caused by simultaneously focusing Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 752693 8 parameter tuning. The MSG stabilizing functional is minimally affected by the a priori model. Zhdanov (2002) proved that the MS and MGS satisfy the Tikhonov criterion for a regularization stabilizer. Since the MSG is a spatial gradient of the MS, it also satisfies the regularization criterion.

Data Sources
Four survey lines (Lines 1, 2, 3, and 4) were identified for MT data acquisition in the Songpan-Aba area. Lines one and two were parallel and trending in a SE−NW direction with an ∼45°a zimuth, whereas Lines three and four were trending in an SW−NE direction with the same azimuth. The survey lines spanned the entire Songpan-Aba area. Field testing was conducted using Model V-2000 devices in a five-component tensor-impedance format, in 2002. In total, 809 points were collected from the 4 MT profiles in the region. The detailed line layout and workload requirements are shown in Figure 2B and Table 2, respectively. The survey was fully implemented using fixed-station and remote-reference methods throughout the survey area. The effective recording bands were 320-0.001 Hz, the electrode distance was 100 m, and the ground resistance was <1 KΩ. The horizontal magnetic rods were buried below 50 cm and the vertical magnetic rods were inserted approximately three-quarters into the ground. Differential GPS positioning was applied to determine the locations of the survey points, maintaining the precision of the plane coordinates within ±10 m. The measurement network was maintained in a regular form, with point distance and lateral deviations within 400 m. Maintenance was conducted to ensure that at least 85% of the physical points were working during the survey.

Electrical Conductivity Structure
First, a theoretical model was designed ( Figure 5A) with the resistivity of wall rock at 100 Ω m and the resistivity of a highconductivity strip at 5 Ω m. The above-mentioned twodimensional MT-regularized inversion method, with inversion parameters based on the constraints of the MSG model, was used to perform the inversion for this theoretical model. The horizontal grids were set up on a 10-km scale and the vertical grids were divided into 86 layers, which were equally spaced in a log scale up to 300 km. The model response incorporated 5% Gaussian random noise, with 40 frequencies from 320-0.00055 Hz. The initial model was an infinite half-space of 100 Ω m. The inversion results are shown in Figure 5B, which illustrates that the resistivity and distribution of the high-conductivity strip can be accurately delineated in the model, with the fitting error of the inversion data reduced to 1.06. Further, this two-dimensional MT-regularized inversion method was applied to the inversion of the four profiles in the Songpan area. The inversion parameters were as follows: the horizontal grid-scale was the spacing of the measurement points, and the vertical grids were divided into 86 layers, which were equally spaced in a log scale up to 300 km. The maximum number of iterations in the inversion was 50, and the error level was 5% (i.e., errors with an initial value < 5% were set to 5%, while >5% remained unchanged). The initial model was an infinite half-space of 100 Ω m. The four profiles converged to 1.5, 1.8, 1.3, and 1.2, respectively, using the chi-square test error as the convergence criterion. As shown in Figures 6−9 B, the resistivity structure of the four profiles revealed a vertical distribution of high-resistance-low-resistance-high-resistance. A high-resistance layer ∼20 km thick existed in the shallow ground, composed of sedimentary materials and substances from the upper crust. A high-conductivity layer was located 20-60 km below the study area, while heterogeneous highresistance and high-conductivity anomalies appeared below 60 km. A connection was found between the highconductivity layer (20-60 km) and the high-conductivity anomalies below 60 km; therefore, it was concluded that the high-conductivity layer between 20 and 60 km might be formed by a substance from the deep mantle which was brought to the crust through upwelling channels. A comprehensive elaboration on this conclusion, integrating the magnetic anomaly inversion results will be presented in the following section.

ANALYSIS AND INTERPRETATION
First, the planar distribution of magnetization intensity was extracted from the 3D magnetic results at different depths (10, 15, 20, and 30 km) (Figure 10) for the comparison of resistivity  −9) to analyze the magnetic and electrical resistivity structure of the study area and its geological features. The 10-km magnetization map showed that highly magnetic matter existed in the South China plate, the eastern part of the West Qinling Mountains, Sichuan, and the southern part of the Longmenshan massif, where a wide range of primarily Indosinian period rock outcrops were shown on the surface, with a small quantity from the Hercynian and Yanshanian age. The distribution of exposed rock adequately corresponded to the magnetic anomaly. The distribution of weakly magnetic substances was identified in the southwestern parts of the Songpan-Aba massif, Bayankara, and the Longmenshan massif, whereas the major parts of the West Qinling and Longmenshan massif did not contain substantial magnetic structures.
Next, from the 15-30 km magnetization map, it was observed that the magnetic bodies became significantly fewer and smaller, with only one magnetic body detected in the southern part of the Longmenshan massif at 20 km deep. The magnetic layers in the study area were distributed above 20 km. The lower boundary of a magnetic body is determined by the depth of the Curie isotherm. Temperature increases with depth, and when the temperature exceeds the Curie point of the enclosed ferromagnetic minerals (primarily magnetite in the crust), they become paramagnetic and no longer generate magnetic anomalies. The Curie temperature of magnetite in the crust is ∼580°C; therefore, the 580°C isotherm is approximately located at the bottom interface of the deep magnetic anomaly. The upper and lower boundary conditions based on the upper mantle temperature are estimated by surface temperature and the S-wave velocity from ground stations. Sun et al. (2013) estimated that the typical temperature range was Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 752693 11 800-1,000°C for the Moho below the Tibetan Plateau, while the thermal profile of 90°E, 30°N in southern Tibet indicated that the Curie isotherm depth was ∼20 km, consistent with our magnetic inversion results. Figure 4 illustrates a magnetic connection linking the upper and lower compartments in the southern part of the Longmenshan massif. It was tentatively deduced that the magnetic source in the Songpan-Aba area originated from a magma root in the southwestern part of the Longmenshan massif (red circles in Figures 10C,D).
Beneath the northeastern part of the Tibetan Plateau (West Qinling, Songpan−Aba, and Longmenshan Mountains), we found an extensive low-resistance layer below 20 km that showed no magneticity, possibly due to partial melting of the crust, which might have been sourced from the substance extruded from the eastern margin of the Tibetan Plateau that was pushed eastwards under the south−north compression from the Indian and the North China plate. Subsequently, blocked by the Sichuan Basin, the lower crustal low-resistance substance congregated and lifted near the Longmenshan area, contributing to the crustal uplift, consistent with the crustal flow model (Royden et al., 1977;Clark and Royden, 2000).

CONCLUSION
In this study, the magnetic and electrical structures of the Songpan-Aba area were obtained from 3D oblique magnetization and MT regularization inversion using high-precision total magnetic anomalies and four long-term MT survey lines. Furthermore, we explored the origin of the magnetic layer in the southern part of the study area and the formation of zones in the northern area with or without weak magneticity. 1) A continuous magnetic layer exists above depths of 20 km under Songpan-Aba and adjacent areas on the south; it is located above the low-resistance layer in the crust and can be traced to the west side of Songpan, deepening in the northwestern part of the Longmenshan area, and possibly originating from a magma root in the southwestern part of the Longmenshan massif. 2) There are extensive low-resistance and weakly magnetic layers below 20 km in the West Qinling, Songpan-Aba, and Longmenshan areas that may be partially molten crustal material sourced from substance extruded from the eastern edge of the Tibetan Plateau.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CZ performed the experiment, the data analyses and wrote the article. LZ contributed significantly to analysis and article preparation. PY performed the analysis with constructive discussions. XX performed the analysis with constructive discussions and plotted the figures.