The Tengchong volcano region, situated in the southeast of the India-Asia collision zone, hosts an active intraplate volcanic cluster (Fig.1)^[1], where the latest volcanic eruption occurred in 1906^{[2, 3]}. The young Holocene volcanoes north of Tengchong, such as the Ma’anshan, the Daying Mount and the Heikong Mount^[4], exhibit high K calc-alkali basalt and andesite on the surface, associated with the magma activity caused by mantle melting^[5]. However, the origin of the Tengchong volcano is still controversial. Some researchers ascribed its origin to the subduction of the Burma microplate down to 400 km depth^{[2, 3, 6]}, while others suggested that it is related to the subduction of the Indian plate down to 150 km^{[7, 8]} or 200 km depth^[9].

(Fig.1)

Fig.1 Map showing major faults, epicenters and volcanoes of the study area. Inset in upper left indicates location of the study area and tectonic setting

There are a series of NE-NNE axial fault basins and NS striking normal faults in the study area^[4], where frequent earthquakes attract concerns of the geoscience experts. M7.3 and M7.4 earthquakes successively occurred in the Longling area south of Tengchong on 29 May 1976. The 10 March 2011 earthquake (M5.8) killed 25, injured more than 250 people, and displaced about 300000. The widespread hot springs indicate the stronger geothermal activity and the possibility of volcanic activity in the future^[9]. Therefore, probing and research of the deep complex structure of this region are of great significance for predicting volcanoes and earthquakes.

According to geological and geophysical data, the Tengchong area is featured by high heat flow, low resistivity, low Q value and low velocity anomalies^{[9, 10, 11]}. The tomographic results of Lei et al. show the low velocity anomalies under the Tengchong volcanic area extend down to ~400 km depth^[2]. The velocity structure of the Sichuan-Yunnan region of Huang et al. suggest that the Pn-wave velocities are closely associated with the recent tectonic movement and geothermal activity^[12]. Jiang et al. analyze the focal mechanisms of small shocks in the Tengchong volcanic region, and find that earthquakes concentrate around the Ma’anshan and Atami area, and consider that earthquake activity is closely related to volcanic activity^[13]. Lei et al. suggest that the 21 July 21 (M6.1) and 16 October 16t (M6.2) 2003 Dayao earthquakes and the three 2011 Yingjiang earthquakes (M > 4.0) are related to the upwelling of hot and wet material caused by dehydration of the subducting Indian plate^{[2, 3]}. The above research results supply some explanation on the velocity structure and seismogenesis in the Tenchong volcanic area and surroundings. They, however, did not provide fine velocity structure of this region.

In order to obtain more detailed 3-D velocity structure in the study area, Institute of Geology and Geophysics, Chinese Academy of Sciences carried out an 1-year portable seismic array observations since December 2002 and acquired large amounts of earthquake data, funded by the National Natural Science Foundation of China. In this study, we used the finite difference tomography method to inverse the P-wave velocity structure beneath the Tengchong area, with the data recorded by mobile networks and fixed stations, and discussed the relationship between the deep structure and the genesis of the Tengchong volcanoes and the M7 Longling earthquakes.

The study region spans 24°N-26°N and 97.5°-100°E. The rectangular coordinate system is adopted in tomography calculation. Data is converted from latitude and longitude coordinates to rectangular coordinates about the origin (97.5°E, 24°N). After calculating, the data range is 0~253.5 km in east-west direction and 0~222.1 km in south-north direction, respectively. This area covers the main Tengchong volcanoes, the source region of the Longling earthquakes and main faults, such as the Nujiang, Tengchong, Dayingjiang and Longling faults. Mobile network was composed of 24 temporary stations (Fig.2), in which each wide-band seismography equipped with a RefTek data collector (72A or 130) and a three-component seismometer (Guralp CMG-3ESP). The continuous seismograms were recorded for more than one year since December 2002, with a sampling rate of 50 sps.

(Fig.2)

Fig.2 Distribution of stations and earthquakes used in this study

According to the local earthquake catalogue, we selected the seismic waves from the data recorded by the mobile network and picked up the first arrival times. The earthquakes satisfy the following criteria are selected from the above data sets: (1) The earthquake epicenters are located in the study area. (2) Magnitudes are greater than 2.0. (3) Seismic waveforms have high signal-to-noise ratio, so that seismic phases are distinguishable. (4) Focal depth is less than 25 km, and recorded by more than 5 stations. According to the criteria, a total of 1587 P-wave first arrivals from 106 local earthquakes were picked up. The picking accuracy of the P-wave first arrivals is estimated to be 0.02 s.

Considering the short time of recording by the mobile network and insufficient arrival time data, we also collected the events which occurred from 1990 to 2009 recorded by 7 fixed stations. Hypocenter parameters and P-wave arrival times were taken from the bulletin of seismological observation of Yunnan province. In total, we picked 4471 P-wave first arrivals from 1494 local earthquakes, each of which was recorded by more than 4 stations.

Next, we integrated the data by fixed stations and that by mobile stations. In order to ensure the uniform distribution of the earthquakes, we only selected the earthquakes recorded by the most stations in a sphere with radius of 5 km, according to the method of homogenization treatment proposed by Liang et al.^[14]. Finally, we got 2829 P-wave first arrivals from 582 local earthquakes, each of which was recorded by more than 4 stations. Fig.2 shows the distribution of stations and earthquakes used in this study.

We used the finite difference method to calculate travel times of seismic waves. This method was proposed and developed by Vidale in 1988^{[15, 16]}. Afterwards Hole improved it to suit for complex velocity structure^{[17, 18, 19]}. The inversion is achieved by using the LSQR algorithm^[20], with damping and first-order smoothing regularizations, to resolve the large and sparse system of equations. The LSQR algorithm is economical on regularized problems and appears to be more reliable if the residual is not small. It is suitable for numerical analysis of ill-posed problems in the iteration process. Additionally, the damping factor and the smooth factor are used to ensure the stability of solution, because uneven distribution of data could cause unreasonable local anomalies in inversion results. Compared with the ray tracing tomographic method, the images at the middle depth among the vertical grids can better reflect the velocity perturbation in the smoothed inversion results by finite difference tomography, instead of velocity perturbation values at grid nodes.

The model space is in the range of -30~280 km in the east-west direction, -30~250 km in the south-north direction and -4~40 km in depth, larger than the scope of the target area. The model was parameterized with a grid spacing of 1 km×1 km×1 km. With this model, we used the finite difference method to synthetic arrival times. In the inversion, we redivided the model space according to the situation of ray distribution. Because of the dense rays in the shallow position, the grids for inversion are smaller; whereas, if the rays are relatively sparse due to the fewer events in the deep crust, we set larger grids. Finally, the model was parameterized with a grid spacing of 20 km×20 km in horizontal directions and with grid nodes at depth of 2, 8, 14, 20, 26 and 33 km. There were 16×14×7 grid cells for inversion. Velocity perturbations at grid nodes were taken as unknown parameters that were solved in the inversion.

Our starting 1-D velocity model (Fig.3) was inferred from previous studies of deep seismic sounding^{[9, 21, 22]} and seismic tomography^{[10, 23, 24]} in the Tengchong area. In Fig.3, the solid line represents the starting 1-D velocity model. After adjusting on the basis of residual distribution of observational data, we obtained the adjusted 1-D velocity model shown as the dotted line in Fig.3. With the adjusted 1-D model, we used the finite difference tomographic method to invert the 3-D velocity perturbations by three times of iteration.

(Fig.3)

Fig.3 Starting velocity model

We carried out the method of the checkerboard resolution test (CRT) to assess the resolution, including the following three steps.

(1) In CRT, ±5% velocity perturbations are alternately assigned to the grid nodes along the three component directions. The model space is divided with a grid spacing of 20 km×20 km in the horizontal directions and with grid nodes at depth of 2, 8, 14, 20, 26, and 33 km. Alternating velocity perturbation boundaries are the same as the 16×14×7 divided grids for inversion, at whose nodes the velocity perturbations with the same amplitude with positive and negative values.

(2) For this model we calculated the synthetic travel times which were taken as the observation travel times.

(3) Then, we utilized the above tomographic method to reconstruct the checkerboard model using the 1-D initial model without any velocity anomalies. The results at the small grid nodes are not smoothed, so they are equal. Fig.4 shows the results of the CRTs at depth 1, 5, 11, 17, 23 and 30km (any depth between grid boundaries). For depth 33~40 km, the CRT results are not listed because of the poor resolution.

(Fig.4)

Fig.4 Plots of recovered checkerboard velocity models

Figure4 shows the results of the CRTs with a horizontal grid spacing of 20 km. The resolution is considered to be good for areas at depth of 0~33 km where the checkerboard model has good ray coverage, according to the distribution of earthquakes and stations. The resolution is high in the central of the study area for this test, especially in the Tengchong volcanic area. In this study, we added the seismic data recorded by more than 20 mobile stations for one year to improve the uniformity and rationality of ray distribution to enhance resolution of the results.

Figure5 shows the P-wave velocity perturbations of the obtained tomographic images in six horizontal sections at depth 0~33 km, relative to a 1-D velocity model as shown in Fig.3.

(Fig.5)

Fig.5 P-wave velocity perturbations on horizontal sections

At 1 km depth, obvious low-V anomalies occur around Tengchong-Yingjiang, Tuantian-Luxi, the east of Baoshan and the southeast of Shidian area. A high-V anomaly exists east of Tengchong to Shidian. Another small high-V anomaly is visible to the north of Gudong.

At 5 km depth, prominent high-V anomalies exist on the area of Tengchong-Yingjiang, the south of Luxi and Daojie-Baoshan, while the low-V anomaly distributes mainly between Gudong and Tuantian.

At 11 km depth, low-V anomalies are visible in the area of Xincheng-Mengyang-Luxi, Tuantian-Shidian, the north of Yongde and the north of Baoshan. There exist high-V anomalies to the north of the Dayingjiang fault, to the west of the Tengchong fault and in the Yongping-Shidian area. The two Longling earthquakes lie on the boundary between the high and low velocity anomalies.

At 17~30 km depth, the scopes of the high-V anomalies to the north of the Dayingjian fault and to the west of the Tengchong fault become smaller, even the anomalies turn into the low-velocity anomalies with increasing depth. A low-V anomaly to the southeast of Shidian gets smaller with increasing depth, meanwhile turns into the high-V anomaly on a large scale at 30 km depth. Another low-V anomaly exists under the Tengchong volcano and Tuantian-Shidian, which extends to the deep subsurface.

Figure6 shows three vertical cross sections. Section A passes through the hypocenters of the Longling earthquakes. Section B and C pass through the Tengchong Volcanic area along the east-west and north-south trends.

(Fig.6)

Fig.6 Vertical cross sections of P-wave velocity perturbation tomography

In Section A, a high-V anomaly is visible under the Tengchong Block and the Baoshan Block. A low-V anomaly forms a band down to depth from west to east between the Longling fault and the Nujiang fault. The values on both sides of these two fractures are close to or higher than that in the initial model. The Longling earthquakes occurred near the velocity boundary beneath the faults. There is a low-V anomaly beneath the Shidian basin at 7~8 km depth, which extends to the east to connect with the low-V anomaly beneath the hypocenters of the Longling earthquakes. A similar low-V anomaly exists under the Yingjiang basin down to about 10 km depth, then turns into a high-V anomaly in the middle-lower crust.

In Sections B and C, the tomographic images show a small-scale low-V anomaly beneath the volcanic area down to 5 km depth, where a high-V anomaly is visible at 5~15 km depth. Below 15 km depth, an obvious low-V anomaly extends down to deep. In Section B, our model shows a prominent low-V anomaly at depth of 5~15 km to the east of the volcanoes. The Baoshan basin is characterized by a low-V layer at 0~8 km depths and a high-V anomaly at 8~20 km depth, to the west of which a low-V band extends to the deep crust. In Section C, a low-V anomaly exists and extends to the north in the middle-lower crust beneath the volcanoes. A high-V zone is present at 5~15 km depth under Gudong. The velocity perturbation between the Yingjiang fault and the Longling fault is lower at 0~14 km depths and higher below the depth of 14 km than that under the volcanoes.

The Tengchong area has experienced high-temperature metamorphism, magma intrusions and volcanic activities since Cenozoic time. In the recent years, earthquakes have occurred frequently. The deep structure and the magmatic activity beneath the Tengchong area have attracted broad attention of the geoscientists. Meanwhile, the structural genesis of the Longling earthquakes is also a focused issue in geosciences. The seismic velocity structure of the upper-middle crust can be better constrained when the data by mobile stations and fix stations are used together. Our tomographic model provides new information on the magmatic activity in the Tengchong volcanic area and the structural genesis of the Longling earthquakes.

It is commonly considered that low and high velocities beneath the volcanoes represent the unsolidified magma capsule or the partial melting magma and cooled and solidified magmatic intrusions or the less volatile ultramafic remnants, respectively^{[25, 26]}. The examples include the Hawaii Volcano, New Zealand’s Mount Taranaki volcano, Ecuador’s Tungurahua Volcano, the Yellowstone Volcano, the French Mont Dore, and so on^{[27, 28, 29, 30, 31, 32]}. Our results suggest that a low-V layer with small thickness between the surface and a depth of 5km may be caused by Tertiary and Quaternary basin sedimentary and volcanic deposits, consistent with previous tomographic images^{[9, 10, 21, 22, 24]}. The high-V anomaly at depth of 5~15 km is observed at approximate locations in other tomographic results, which is visible between 9 and 15 km depth, in the lower crust and at depth of 5~10 km in the results by Qin et al., Wang et al. and Yang et al., respectively^{[21, 23, 24]}. Yang et al. inferred this high-V anomaly as the cooled and solidified magmatic intrusions or the less volatile and high-density remnants^[24].

There is a lot of evidence that suggest the existence of magma chambers under the Tengchong volcano. Kan and Zhao have proved that a magma chamber underlies the Tengchong volcano on the basis of recent volcanic movement, different physical properties and focal mechanisms associated with magma impact^[1]. The tomographic results by geophysicists show the existence of magma chambers or the paths of magma intrusions in the deep crust and mantle of this area^{[2, 9, 10, 21, 22, 24]}. There are abundant surface manifestations in the Tengchong area, including the high-K calc-alkaline basalt, to indicate that the low-V band may be closely associated with the path through which the mantle magma upwells^[33]. From our tomographic images, we can see a high-V anomaly between 5 and 15 km depth under the Tengchong volcano related to the high-density remnants, and a prominent low-V anomaly extending to deep in the middle-lower crust in Fig.6(b-c), which is inferred to be the molten or semi-molten magma. In Fig.5f, an obvious low-V anomaly is located between the Tengchong fault and the southern segment of the Nujiang fault, and joined to another one beneath the Tengchong volcanic area. Our study provides a new piece of evidence that both of the lava flows represented by the low-V anomalies come from the same magma source in the mantle. Lei et al. suggested that the Tengchong volcano is a rift-related volcano caused by the subduction and dehydration of the subducting Indian plate and corner flow in the mantle wedge^[2], which better supports the shape of the distribution of the low-V body along the faults.

The M7.3 and M7.4 earthquakes occurred successively in Longling, Yunnan province on 29 May 1976. According to the National Earthquake Network, these two earthquake’s epicenters were at coordinates (24.5°N, 99°E) and (24.6°N, 98.7°E), nearby the Nujiang fault and the Longling fault which are strike-slip faults with steep inclination^[34]. Different studies provided different information about the depth of the Longling earthquakes, such as 12 km by Chen et al.^[35], 20 km according to China Earthquake Network Center and Yunnan seismic station network, and the 15 km depth determined through focal mechanism solutions by Harvard University.

The velocity perturbation images in Fig.6 show the heavily distorted velocity contours and complex deep structure in the source region. On the low-V and high-V boundary, there exist the exceptional structural conditions for earthquake generation and occurrence. A large number of magmatic intrusive rocks in different periods have outcropped in the Cambrian basement between the Nujiang fault and the Longling fault, underlain by a low average velocity in the crust in the vertical sections. This evidence demonstrates that multiperiodic hot-tectonic events have occurred in history^{[24, 34]}. In the western edge of the Baoshan block, the basement strata to the east of the Nujiang fault formed in Cambrian time. In the eastern edge of the Tengchong block, to the west of the Longling fault exposes the Gaoligong high-grade metamorphic granite at the surface. The geological features in these two regions correspond to the high-V structures in Fig.6, with large strain strength. Thus it can be seen that the variation of velocity structure in the crust reflects the heterogeneity of strain strength of rocks, which causes the uneven distribution of cumulative stress. Strong earthquakes can occur in this tectonic setting. The distribution of velocity anomalies in our model is in agreement with that in the tomographic image by Yang et al.^[24]. Comparing with the results by Bai et al.^[22], there is some difference in the depth of hypocenters. However, there is a common understanding to the information that the important tectonic condition of earthquake preparation and occurrence is the inhomogeneous distribution of cumulative stress and energy around the boundary between low- and high-V anomalies.

In addition, the intersection of the Longling fault and the Tengchong fault near the source region is underlain by a low-V anomaly tending to extend downward to the lower crust and upper mantle. Because controlled by multi-deep fractures, this area is characterized by the development of fractures and the invasion of liquids. It is probable that the magmatic activity beneath the volcanoes and fluid intrusions in the deep crust of the Longling earthquake zone are from the same magma source.

A high-resolution tomographic model of the crust under Tengchong has been obtained using earthquake data recorded by a temporary network and other permanent stations, and the finite difference tomographic method. The results provide new insights into the magma activity in the Tengchong volcanic area and the origin of the Longling earthquakes with M7.3 and M7.4.

The average velocity under the Tengchong volcanic region is relatively low. With increasing depth, the velocity varies strongly. The low-V zone in the shallow subsurface is mainly related to the basin deposits, volcanic deposits, fluid fractures and hot springs. The high-V anomaly between the depth of 5 and 15 km is inferred to associate with the cooled and solidified magmatic intrusions, the less volatile ultramafic remnants, or the remnants of the ultramafic accumulative rocks in the old volcanic pipe. The low-V body at depth represents the molten or semi-molten magma, where the present magmatic activity is ongoing. Besides, the magma from the same source in the mantle probably intruded upward along the intersection of the Tengchong fault and the Yingjiang fault and the intersection of the Longling fault and the Nujiang fault, which may cause the volcanic eruption and the Longling earthquakes.

The Longling earthquakes occurred near the intersections of the faults, recognized as a boundary between the low- and high-V anomalies in the tomographic images. The magmatic activity at depth of this area is controlled clearly by the Longling fault and the Nujiang fault which are deep fractures. The magmatic activity results in the decrease of the rheological strength of rocks in the crust. But the strain strength is higher in the crust on both sides of the faults, where stress accumulates. We conclude that the energy in the lower crust and upper mantle could not be released after the termination of volcanic activity, so the above-mentioned tectonic position became the locus of mass stresses. Therefore, the strong earthquakes in the Longling area can be attributed to the lateral variation of the crustal medium strength, the heterogeneity of stress accumulation and the deep magmatic activity.

We thank Liu Jianhua and Yang Xiaotao at the Institute of Geology and Geophysics, Chinese Academy of Sciences for their valuable suggestions. This work was supported by the National Natural Science Foundation of China (40974027), the National High Technology Research and Development Program (2009AA093401), and Postdoctor Foundation of Shandong Province (201103087).