Density structure of the crust in the Emeishan large igneous province revealed by the Lijiang‐ Guiyang gravity profile

The Emeishan large igneous province (hereafter named by its acronym ELIP) is the first accepted large igneous region in China. The current study tries to reconstruct the density structure of the crust in this region. For this purpose, we conducted the gravity survey along an 800‐km‐long profile, which stretched laterally along the latitude 27°N from Lijiang (Yunnan province) to Guiyang (Guizhou province). The fieldwork included 338 gravity measurements distributed from the inner zone to the outer zone of the mantle plume head. After a series of gravity reductions, we calculated the Bouguer gravity anomaly and then constructed the density model for ELIP by iterative forward modeling from an initial density model tightly constrained by wide‐angle seismic reflection data. The topography of the Moho, here physically interpreted as a density discontinuity of ~0.4 g·cm‒3, gradually rises from the inner zone (~50 km deep) to the outer zone (~40 km), describes a thicker crust in the inner zone than in any other segment of the profile and largely reproduces the shape of the Bouguer gravity anomaly curve. Both the Bouguer gravity and the density structure show significant differences with respect to the inner zone and the other two zones of ELIP according to the commonly accepted partition of the Emeishan area. A thicker and denser middle‐lower crust seems to be the main feature of the western section of the profile, which is likely related to its mafic magmatic composition due to magmatic underplating of the Permian mantle plume.


Introduction
The Emeishan large igneous province (ELIP) is generally recognized as the first accepted large igneous region in China (Xu YG et al., 2001). ELIP has drawn strongly the attention of the scientific community for several reasons, among them by its possible synchrony with the eruption of the Siberian Traps and its possible relationship with the mass extinctions during the Late Permian (Chung and Jahn, 1995;Chung et al., 1998;Wignall, 2001;Xu YG et al., 2001;Ali et al., 2002;Courtillot and Renne, 2003;Zhou MF et al., 2002;Lo et al., 2002). Chung and Jahn (1995) early invoked the plume model to explain the outpouring of the Emeishan flood basalts. In more recent years, a rigorous evaluation of the role played by a mantle plume in the generation of this large igneous province was addressed Xu YG et al., 2004;Zhang ZC et al., 2006). Xu YG et al. (2007) summarized the achievements reached in the research of ELIP and listed up to eight "diagnostic" arguments in favor of the involvement of a mantle plume in ELIP, namely: (1) uplifting prior to volcanism; (2) chemical characteristics of the magma; (3) high-temperature magma; (4) thermal zonation of the plume; (5) short duration of flood volcanism; (6) dimension and volume of magmas; (7) physical volcanology; (8) radiating dike swarm and age progression along hotspot tracks. These arguments come mainly from geochemistry, stratigraphy and tectonics. The volcanism generated by the ancient mantle plumes would lead to thermal anomalies in the crust and the mantle lithosphere, and also to the intrusion and eruption of magmas, so that the material composition of the crust and mantle would change correspondingly, which are the imprints of the events provoked by ancient mantle plumes. Nevertheless, the snapshot obtained through geophysics provide the most direct images of the deep structure of the Earth. The main features can be recognized by comparing variations in seismic velocity, Vp/Vs ratio, density, anisotropy, rheology, and magnetism with the properties in the surrounding areas. Recently, the mantle origin in ELIP has been studied from the analysis of residual gravity anomalies, so it has suggested the melting of ultramafic rocks at high temperatures intruded by a mantle plume and an abundant transport of basalts (Deng YF et al., 2014b). However, the detailed density structure across Emeishan regions is still unknown. If a mantle plume emerged beneath Emeishan in the past, at least some physical traces should remain in the region and might be observed from a snapshot.
To meet this target, we used the gravity method to obtain the density structure of the crust. The data acquired by means of a gravity profile conducted along the Guiyang-Lijiang transect in ELIP were used to adjust a density model through forward modeling, that we used later as a basis to explain the geophysical features and the geodynamics of the region.

Characterization of the Emeishan Basalts
The Emeishan basalts are erosion remnants of mafic composition rocks generated by successive volcanic episodes that occurred during the Late Permian in the western margin of the Yangtze Craton, southwest China. They cover roughly a 250000 km 2 rhombic area that extends east of the Sanjiang tectonic belt, bounded by the Longmenshan thrust fault and the Ailaoshan-Red River (ASRR) slip fault to the northwest and southwest, respectively ( Figure 1). Nevertheless, some basalts and mafic complexes existing in the Simao basin, North Vietnam (west of the Ailaoshan fault), Qiangtang terrain, Lijiang-Yanyuan belt and in the Songpan active fold belt (northwest of the Longmenshan fault) are possibly an extension of ELIP (Chung et al., 1998;Xiao L et al., 2003Xiao L et al., , 2004Hanski et al., 2004). Some Emeishan-type basalt traces in southwest Yunnan and northern Vietnam might be related to the continental extrusion of Indochina along the Ailaoshan-Red River fault zone (Tapponnier et al., 1990;Chung et al., 1997).
The thickness of the Emeishan basalts ranges from about 5000 m in the west to a few hundred meters in the east. From sedimentological and paleogeographic data , ELIP has been divided into three zones, namely: the inner zone, according to the extent of the erosion of the Maokou Formation composed of Mid-Late Permian carbonates, and the intermediate and outer zones (Figure 1), which reveal a circular dome-shaped geometry (Deng YF et al., 2014b). The Emeishan volcanic successions overlie unconformably to the Permian carbonates, which appear covered by Permian sedimentary rocks to the east and Late-Triassic sedimentary rocks in the central part of ELIP.
Based on the analyses of the samples from Binchuan and Ertan, the Emeishan flood basalts have been divided into two major magma-types: high-Ti (Ti/Y>500) basalts and low-Ti (Ti/Y<500) basalts (Xu YG et al., 2001). Generally, the basalts belonging to the high-Ti (HT) group have higher Ti/Y (>500) and TiO 2 (>3.7 wt%) than those of the low-Ti (LT) group, which have low Ti/Y (<500) and TiO 2 (<2.5 wt%). The large volume of data obtained from basalts samples in Binchuan and Jinping has confirmed this classification (Xiao L et al., , 2004. Moreover, the available data on characteristic trace elements permit a finer subdivision of the LT group into the LT1 and LT2 types. In general, the LT1-type basalts have higher Mg# (67-61) than those of the LT2-type (54-48). For the basalts of the HT group basalts Mg# varies between 58 and 44. The distinction between LT1-and LT2-type lavas is clearly justified when plotting the Th/Nb ratio versus the Ti/Y ratio, which highlights the abundant content of highly incompatible elements in the LT1-type lavas. The LT1 lavas have a ratio Th/Nb >0.17 that is higher than in the case of LT2 lavas and HT basalts (Th/Nb<0.17). Furthermore, the measured and age-corrected 87 Sr/ 86 Sr and 143 Nd/ 144 Nd ratios found in the Emeishan basalts define an array that lies near the "mantle correlation line", although which is somewhat displaced on the right (Xu YG et al., 2001). In general, the HT basalts have relatively higher ε Nd (t) and lower 87 Sr/ 86 Sr(t) values than the LT basalts.
While the age of the Emeishan basalts is constrained by the above and below strata, its absolute age is probably 260.4±0.4 Ma (Gradstein et al., 2004). Courtillot and Renne (2003) reviewed the age of the main flood basalt events on Earth, and advanced the hypothesis that most of these events have the duration of less than 1 Ma and are linked to major bioclimatic events. Recent Ar-Ar radiometric dating on Permian basalts exposed in western Guangxi (Fan WM et al., 2004) and volcanic rocks near Qiaojia yielded 254-256 Ma and 258 Ma, respectively (Xu YG et al., 2007). These estimates are quite consistent with the stratigraphic constraints and strongly suggest that the emplacement of the Emeishan basalts took place prior to the Permian-to-Triassic transition.
Petrogenetic assessment suggests that the sediments in the Xuanwei Formation come from the erosion of the Emeishan volcanic rocks in the inner zone of Emeishan . This indicates that the Emeishan volcanism is a boundary event that has important implications for the duration of this large igneous province. Thus, the age of 259 Ma for the Middle-to-Late Permian transition (Zhong YT et al., 2014) can be taken as the time of the main phase of the Emeishan flood volcanism. Moreover, the short duration of the Emeishan volcanism is another feature that also

Model
In the explored area, the Bouguer gravity anomaly is mainly negative and is spatially distributed as shown in Figure 2 on a map extracted from the original Bouguer gravity map for South China and surrounding regions, well known as the Earth Gravitational Model (EGM2008), which was computed by Pavlis et al. (2008Pavlis et al. ( , 2012. The natural tectonic boundaries of the zones that make up ELIP are roughly parallel to the contours of the regional gravity anomaly, which decreases gradually (in absolute value) from the inner zone to the outer zone. A larger gravity anomaly of negative sign can be observed in the middle of the inner zone. Specifically, the extent of the anomaly is larger (up to -200 mGal) within the inner zone and then decreases to -300 mGal in the middle zone. High gravity anomaly often means that the Moho is shallow, or mantle material upwelling. In this area, the Moho is not shallow, so we performed a gravity experiment to obtain the density structure of the crust.

Data Acquisition
In this frame, we carried out a fieldwork along an 800-km-long profile that was deployed across the inner zone, intermediate zone and outer zone of ELIP from July to August 2012. The spacing gap between each two consecutive gravity measure points was about 2.2 km on average. The gravity measurements, 338 in total (Figure 1), were made with an accuracy of 10 μGal using a Burris Gravity Meter (No. B65). The gravity readings were recorded with respect to two base points belonging to the National Gravity Network of China: one near Guiyang (eastern end of the profile), and another near Lijiang (western end of the profile). Both points were used to later calibrate the drift error. Each observation point was identified by its GPS coordinates with a level of accuracy given by the instrument used, that was Trimble GeoXM.

Data Processing
First, we corrected all gravity values for tide and drift errors. The tide correction was carried out from the time data and coordinates annotated during the data acquisition process; the drift correction was performed assuming that the meter had a linear drift rate. We used the gravity values provided by the two base points to calculate the total drift value in the course of the experiment, and thereafter we calculated the drift correction assigned to each measurement point by linear interpolation.
Once the time effect was eliminated, we converted the relative gravity values to absolute values. For this process we compared the value obtained at each point of the profile with the gravity value at the base point in Guiyang, and added the difference between both values to the existing gravity value at the base point to determine the absolute gravity value at each point of the profile.
Finally, we made the classical reductions of gravity measures (latitude correction, free-air correction, Bouguer correction, and topographic or terrain correction) to obtain the set of Bouguer gravity anomalies. We took into account the same geographical coordinates of the observation points that those taken for tide correction.
In particular, we considered the typical density of 2670 kg/m 3 for gravity reduction. Figure 3 shows the Bouguer gravity anomaly correlated with the topography along the reference profile.

Initial Density Model
In order to reduce the multiple solutions inherent to any inversion method, we stated the initial density model from previously obtained wide-angle seismic reflection data . According to the position of the deep seismic sounding conducted by Xu T et al. (2015) to investigate the seismic velocity structure of the crust in ELIP, we selected a great part of the Bouguer gravity anomaly curve along a distance of 650 km for subsequent observation-calculation adjustments. The data that we used later for forward modeling are shown in Figure 3b.
Recently, Deng YF et al. (2014a) have proposed different laws for the tectonic units that from South China, all formulated as linear relationships between P-wave velocity and density for the tectonic units forming South China. However, because our study area differs significantly from this large region and particularly from the Yangtze Block, we used the relationship given by Christensen and Mooney (1995) to convert P-wave velocity into density of the crust and the relationship given by Zhu JS et al. (2005)    To start the inversion process, we initially took a density model consisting of a 1 km ×1 km-sized grid model, and associated density values using the known P-wave velocity data and the equations written above. In this way, the reference model is based on the available seismic reflection data . Here we used a three-layer reference model with interfaces defined by the seismic data at hand. Nevertheless, the final density model obtained by inversion is presented with a geometry constituted by irregular blocks of uniform density.

Results
Starting from the initial density structure shown in Figure 4b, we manually alter the density model until the calculated gravity fit to the observed gravity. Both the calculated gravity curve in the first stage and the observed gravity curve are shown in Figure 4a for    comparison, along with the residuals (differences) between them. It is obvious that these two curves do not match. The final density model for the crust and uppermost mantle obtained by forward modeling is shown in Figure 5b. Gravity is not solely a function of the density structure. An infinite number of models could produce the same observations. The fit to data does nothing to make the model more reliable than another model that fits the data. Therefore, the reliability of the model does not come from the coincidence between measured and calculated gravity (Figure 5a), although this is a necessary condition that partially support the inversion result, but from the fact that it honors the structures observed seismically.
The lateral density variation is practically non-existent in ELIP, except in the lower crust in the inner zone where is ~0.2 g·cm -3 with respect to the adjacent zones. The Moho, whose topography (Figure 5b) largely reproduces the shape of the Bouguer gravity an-omaly curve (Figure 5a), is interpreted as a density discontinuity of 0.4 g·cm -3 . In ELIP the Moho turns out to be deeper in the inner zone where it reaches the average depth of 50 km and depths of 45 km and 40 km in the intermediate and outer zones, respectively. The Moho gradually rises from the inner zone to the outer zone where it reaches a depth of ~38 km, and describes a thicker crust in the inner zone than that in any other transect of the profile.
The Bouguer gravity curve (Figure 5a) follows a similar trend to that of the Moho topography (Figure 3a) in the intermediate and outer zones; but the Bouguer gravity anomaly is high (in relative value) in the inner zone while the Moho is deeper (Figure 5b). According to the theory of isostasy, a deeper Moho correlates with a higher topographic elevation, but the topography is not high in the inner zone, so theoretically there must be a denser and thicker middle-lower crust to reach the isostatic equilibrium. He B et al. (2003) divided ELIP into three zones according to the extent of the erosion of the Maokou limestone (Figure 1). Our results based on gravity demonstrate that zonation of ELIP exists not only at surface but also at crustal depths. If we contrast the Bouguer gravity with the terrain elevation, we can divide our reference gravity profile into six different zones ( Figure 6). In the west end of the profile the topography is high (above 2000 m) and the Bouguer gravity anomaly is comparatively low (below -300 mGal). Part of the negative correlation between topographic elevation and Bouguer anomaly could be due to a reduction density that is too high. Tectonically, this segment of the profile is halfway between a part of the Sanjiang belt and the Yangtze Craton, and both the terrain elevation and the Bouguer anomaly take values clearly different from those in other segments of the profile (Figure 6b, points labeled as 1).

Zonation of the Bouguer Gravity
The segments on both sides of the inner zone of ELIP are two other cluster-type zones (Figure 6b, points labeled as 2 and 4) where the elevation is below 1650 m and the gravity anomaly ranges between -290 and -250 mGal. In particular, the zone near the Xiaojiang Fault (Figure 6b, points labeled as 4) that has the smallest elevation of the profile (below 1400 m) shows Bouguer gravity anomaly around -250 mGal. This low anomaly has its origin in a density loss of the upper crust (Figure 5b), which could be the result of faulting movement. The other small cluster-type zone, elevation of 1500 m and Bouguer anomaly of around -280 mGal (Figure 6b, points labeled as 2) could also be most likely the result  Figure 5. Upper plot (a): fit between the observed (red line) and calculated (blue line) gravity curves (the measured gravity was smoothed using a five-points average method twice); residuals (differences) between both curves are also plotted. Lower plot (b): final density structure (in g/cm 3 ) by blocks. Names and acronyms included on top of this plot allow identify places with respect to the inner, middle and outer zones of ELIP.
of faulting.
Most of the profile reflects roughly a correlation between the Bouguer anomaly and the terrain elevation, having a mid-high anomaly of more than -250 mGal (in relative value) in correspondence with a gradual decrease of the topographic elevation (Figure 6b, points labeled as 3, 5 and 6). It is possible to distinguish two zones according to the observed slope: one corresponds to the inner zone ( Figure 6, points labeled as 3) and another corresponds to the intermediate and outer zones of ELIP ( Figure 6, points labeled as 5 and 6). For the same elevation, the anomaly in the inner zone (cloud formed by points marked as 3) is lower (in relative value) than that in the other two zones (clouds formed by points marked as 5 and 6). This difference could be the result of the change in Moho topography.

Middle-lower Crust in the Inner Zone
Our density model for ELIP shows the existence of a comparatively denser and thicker middle-lower crust in the inner zone (Figure 5b), which is the cause of the higher gravity (in relative sense) observed in comparison with the gravity on both sides of this zone (Figure 5a). This is consistent with the circular-shaped high gravity found in the inner zone (Deng YF et al., 2014b), as shown in the map of the regional Bouguer gravity anomaly based on EGM2008 model (Figure 2), which can be the result of the active magma channel during the formation of ELIP. It is also consistent with the results on magmatic underplating and crustal growth in ELIP recently revealed by receiver-function-based passive source seismology ; and also with the results provided by gravity stripping and residual gravity (Deng YF et al., 2016). The thickening of the middle-lower crust detected by a virtual deep seismic sounding (Liu Z et al., 2017) also reflects the accumulation of mafic material in the central part of the plume, and the uplift of the lithosphere above the mantle plume head that triggered relatively high amounts of melting by decompression (Xu YG et al., 2007). The high density may be the result of the solidification of the basaltic magma in the middle-lower crust, which remained there after the eruption. The greater volume and higher density of the middle-lower crust resulted in a huge mass integrated by mantle rocks that caused the curved bending of the Moho in the inner zone.

Conclusions
In order to obtain the gravity response and density structure of the crust in ELIP, a gravity experiment was performed along the Guiyang-Lijiang profile which included 338 gravity measurements. The density model was determined by adjusting the measured Bouguer gravity anomaly starting from an initial density model constrained by wide-angle seismic reflection data. The density model provides new data to better understand the deep structure and features of ELIP. Both the Bouguer anomaly and the density structure show significant differences among the inner, intermediate and outer zones that make up ELIP. The lateral density variation hardly exists in ELIP, except in the lower crust in the inner zone where is ~0.2 g·cm -3 higher than the adjacent zones. The Moho, which is interpreted as a density discontinuity of ~0.4 g·cm -3 , turns out to be deeper in the inner zone where it reaches the average depth of 50 km. There is evidence of a thicker and denser middle-lower crust in the inner zone. Our results demonstrate that it is feasible to zonate the region at crustal depths and even distinguish up to six different zones west-to-east in ELIP; in particular, two narrow fault zones (Chenghai fault and Xiaojiang fault) on both sides of the inner zone contribute to infer a gradual uplift of this last zone driven by thickening of the middle-lower crust up to construct a horst. The thickening reflects the accumulation of mafic materials in the central part of the plume as well as vertical crustal growth. The high density would be the result of the solidification of the basaltic magma in the middle-lower crust in the inner zone. The high density, along with high velocity from seismic data in the inner zone means the material is mafic/ultramafic. The increase of mass and the lifting of the middle-lower crust would be the mechanism that triggers the current buckling up of the Moho in the inner zone of ELIP.