# Spatiotemporal Characteristics of Reservoir-Induced Earthquakes UsingP-Wave Velocity Structures

Authors: Zhiren Feng; Bo Jin; Hongfu Lei; Di Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

In this study, through the seismic phase observation report obtained by a river reservoir seismic network from September 2012 to July 2015, the joint inversion method of source and velocity structure is employed, combined with the regional seismic geological environment data. Such method gives theP-wave velocity structure of different sections in the reservoir area and tries to find out the impact range of the reservoir water depth. And the impact of reservoir water infiltration on crustal medium is discussed. The obtained location results show that the earthquake activity is mainly concentrated in the first reservoir tail and the second reservoir head. As such, the spatial distribution characteristics of the recorded earthquakes are given, and the resolution results of detection board at different depths are speculated and discussed.


---

## Body

## 1. Introduction

The seismic tomography was proposed in 1984, which is a method to invert the three-dimensional velocity structure of the region by using seismic observation data. The earliest application of tomography is medical CT, and then it has been applied in the field of geophysics. Especially, inversion of three-dimensional underground structure by using seismic wave travel time tomography is widely used [1–4].On the one hand, many scholars use seismic tomography to study the velocity structure of strong earthquakes and then discuss the tectonic environment of strong earthquakes [5–9]. Sun and Liu [6] found that the major earthquakes in the Beijing-Tianjin-Tangshan region are mostly distributed on the high-speed blocks or on the side of the high-speed blocks at the intersection zone between the high-speed and low-speed blocks. Zhao et al. [7] considered that the major earthquakes in the crust of Japan from 1885 to 1999 are mostly located at the boundary between the low-velocity and the high-speed zones, which are explained by the results of the tomography. In addition, there are some low-velocity layers on the focal area of the major earthquakes. Zhou et al. [8] used the three-dimensional velocity structure evolution for the Magnitudes 6.2 and 6.1 earthquakes of 2003 in Dayao and demonstrated that there is an NNW (North, North, West)-direction high- and low-wave velocity boundary near the focal area before Dayao earthquake. The focal point is located on one side of the high wave velocity boundary. Pei et al. [9] obtained the lateral change and anisotropy of Pg wave velocity in Yushu area by two-dimensional tomography. The results show that the magnitude 7.1 Yushu earthquake in 2011 occurs in the high-speed abnormal area, and the rupture mainly expands to the low-speed abnormal direction, indicating that the velocity structure controls the occurrence and rupture expansion of Yushu earthquake. In addition, many scholars in China have also applied this method to study the distribution of strong earthquake sequence and the velocity structure characteristics of aftershock area [10–13].On the other hand, many researchers try to use seismic tomography to explore the three-dimensional velocity structure of reservoirs, water flooding, and other areas and discuss the state of water seepage [14–18]. Guo and Feng [14] inverted the three-dimensional velocity structure of Xinfengjiang reservoir, and the results show that the low-speed region corresponds to the surface fragmentation region, namely, the seismic zone region, and the high-speed region has fewer earthquakes, which are closely related to the water infiltration. Through studying the velocity structure of Shanxi Reservoir in Wenzhou, Zhong et al. [15] also found that most of the earthquakes in the reservoir area occurred in low-speed abnormal areas, which may be related to the water infiltration after the reservoir impoundment. Zhou [16] studied the three-dimensional velocity structure of Zipingpu reservoir, and the velocity profile shows that there are high-speed blocks for the depth from 0 to 8 km below the Beichuan-Yingxiu fault, which indicates that the largest depth of fluid along the fault plane penetration is about 8 km and speculates that the Zipingpu reservoir water osmosis may be related to Wenchuan earthquake occurred in 2008. Yang et al. [17] studied the three-dimensional structure tomography of the upper crust in Xinfengjiang reservoir area and found that the microearthquakes in the reservoir area are distributed in specific parts of the physical structure of the medium. The soft and hard interleave medium environment is a possible reason for the occurrence of the microearthquakes of the dip slip normal fault type. Zhang et al. [18] studied the P-wave velocity (Vp) structure in the Zigong-Rongchang flooding area. A transition zone of high- and low-wave velocity near NS strike is formed between the water injection area. The crustal medium P-wave velocity in the water injection area is obviously higher. This phenomenon may be caused by the higher saturation of the underground medium in the water injection area due to the presence of cavity saturation.The rock in the reservoir is porous medium, which can be filled with water or oil. The presence of fluid will affect the seismic parameter characteristics of the rock medium, and its velocity will change with the saturation. Many experiments prove that longitudinal and shear wave velocities are related to the degree of saturation. Gregory [19] believed that the saturation has more influence on the speed of low porosity rock than the speed of high porosity rock. Vp of the rock, which completely saturated with water, is bigger than that of the rock, which partially saturated with water obviously. Shear wave velocity (Vs) does not always decrease with the increase in the degree of saturation but related to the pressure, porosity, and chemical action between pore fluid and rock skeleton. Domenico [20] showed that Vp decreases with the increase in sample density because the saturation increases when the water saturation is low. However, when the saturation reaches a higher value, the porosity of rock increases significantly. This exceeds the velocity change caused by the increase in density, which leads to the increase in Vp significantly. For the S-wave, since there is almost no relation between the shear modulus and saturation, Vs decreases with the increase in saturation when the density increases. Through the rock influent experiment, when water enters the rock, Vp of the rock will decrease firstly and then remain at a low value, but when the saturation is high, it has a great influence on Vp, and Vs is basically unaffected [21, 22]. Kauster and Toksoz [23] believed that the elastic wave velocity has a great relationship with the fracture pore stiffness, and the increase in fracture pore stiffness will obviously increase the elastic wave velocity of rock.Based on the above experiment and previous results, two reservoirs and their adjacent area (27.8°∼28.8°N and 103.3°∼104.1°E) are taken as study area in this paper. Based on observations in the reservoir recorded stations from September 2012 to July 2015, 5230 earthquakes (M ≥ 1.0) are selected in the study area, and a total of 45667 P-wave radiation data source associated with the velocity structure inversion method is used to determine the Vp structure of research area for different sections. The influence area of the reservoir is analyzed to explore the influence of penetration effect of reservoir water on crustal medium.


## 2. Seismic Activity

Based on records of observations from a river reservoir network in Sichuan (Figure1, dark blue triangles) in Southwest China, its storage capacity is about the half of billion km3, where the drainage area is about 406,100 km2, the maximum height of the dam is 270 m, and the installed capacity of power station is 10,200 MW. The spatial distribution of small earthquakes is given in the study area since September 2012. Earthquakes are mainly concentrated in two areas. Area 1 is the region of the first reservoir tail. Long axis of small earthquakes is expanded along NNW-direction (about 50 km), dense distribution area with about 25 km is in the reservoir water, and short axis is about 25 km. Area 2 is the region of the second reservoir head. Small earthquakes are distributed along the river in the NNE direction, with a long axis of about 45 km and a short axis of around 25 km (Figure 1).Figure 1
Spatial distribution of small earthquakes in the reservoir area.Figure2(a) shows the profile of the focal depth of small earthquakes along the long axis CC′ (Figure 1) for the first reservoir tail and adjacent area, and the distribution of the focal depth (Figure 2(b)) and magnitude (Figure 2(c)) of small earthquakes in this area is statistically analyzed. It can be seen from this figure that the focal depth of the earthquake in this region is mainly distributed within 10 km, accounting for 93% of the total number of earthquakes. The earthquakes for magnitude 2.0 and above are mainly distributed with 5–8 km. The area is mainly dominated by small earthquakes, with earthquakes of magnitude 2.0 or less, accounting for 85%.Figure 2
Earthquake statistics for the first reservoir tail and adjacent areas: (a) small earthquakes; (b) depth; (c) magnitude.
(a)(b)(c)Figure3(a) shows the focal depth profile along the long axis AA′ (Figure 1) of small earthquakes in the second reservoir head and adjacent area, and the distribution of focal depth (Figure 3(b)) and magnitude (Figure 3(c)) of small earthquakes in this area is also presented. It can also be seen from this figure that the focal depth of the earthquake in this region is also mainly distributed within 10 km, accounting for 93% of the total earthquake, among which earthquakes with 5 km account for 52.1%. The region is mainly dominated by small earthquake activities, and earthquakes with magnitude less than 2.0 account for 90%.Figure 3
Earthquake statistics for the second reservoir head and adjacent areas: (a) small earthquakes; (b) depth; (c) magnitude.
(a)(b)(c)Figure4 is the M-T and N-T graphs which indicate small earthquakes in two reservoirs, respectively. The first reservoir began filling on October 11, 2012, and the water level increased by 71 m to reach 352 m from October 17, 2012. Since then, the water level has been stable. The water was filled on June 27, 2013, and the water level increased by 15 m on July 6, 2013. The water level rose again in September 2013 when the high-water level was submerged to the part of the first reservoir tail [24]. To July 31, 2015, the reservoir seismic network has recorded a total of 964 earthquakes of magnitude 1.0 and above in this region, among which 823 times for M1.0–1.9, 128 times for M2.0–2.9, and 13 times for M3.0–3.9.Figure 4
Crustal mean velocity model in the study area.On November 16, 2012, the second reservoir dam started to hold the water, and the diversion tunnel overflowed the water. As a result, the water level in the reservoir area rose by 28 m. From May 4, 2013 to July 29, 2013, the water level raised from 440 m to 554.6 m and has remained at this level since then [25]. To July 31, 2015, the reservoir seismic network has recorded a total of 4072 earthquakes of magnitude 1.0 and above in this area, among which 3651 times for M1.0–1.9, 361 times for M2.0–2.9, 49 times for M3.0–3.9, 10 times for M4.0–4.9, and 1 time for M5.0–5.9. The largest earthquake was Yongshan M5.3 earthquake on August 17, 2014. The frequency of small earthquakes has increased significantly since May 2013.


## 3. Inversion Analysis of Crustal Velocity Structure in Reservoir Area

### 3.1. Inversion Method

Seismic tomography is a geophysical method which uses inversion analysis to deduce the important information such as velocity structure of underground medium and other physical property parameters by analyzing the kinematics (such as travel time and ray path) and dynamics (such as waveform and amplitude) characteristics of observed seismic wave phases.In the parametric representation of 3D velocity structure, grid method is a common method. Its advantage is that the grids can be divided into different spaces according to the density of the distribution of the epicenter and the station, to ensure sufficient ray cross coverage in each grid and reduce the impact of the prior assumptions. In this method, the velocity structure is expressed as a continuous function, the unknown parameter is the velocity value on the 3D mesh node, and the velocity at any point in the model is calculated by interpolation.In the joint inversion of source location and velocity structure, the travel time residualδt is caused by the perturbation of source parameters and velocity. According to the relevant literature [26–28], the problem can be expressed by the following linear equation:(1)δt=Δt+∂t∂xΔx+∂t∂yΔy+∂t∂zΔz+∑n=1N∂t∂vΔvn.Δt, Δx, Δy, and Δz and Δvn represent the source of original times, longitude, latitude, depth of disturbance, and speed disturbance, respectively; N is the total number of velocity parameters. For l earthquakes and K stations, equation (1) can be written in the compact form as follows:(2)δt=Aδv+Bδx,where δt is the m-dimensional travel time residual vector; δv is the n-dimensional node velocity perturbation vector; δx is the 4l-dimensional source parameter perturbation vector; A is the m × n-dimensional travel time partial derivative matrix with respect to velocity; and B is the m × 4l-dimensional travel time partial derivative matrix with respect to source parameters.According to the basic equation (2) of joint inversion, velocity parameters and source parameters are coupled to each other. For inverting two parameters of different dimensions in the same equation at the same time, the numerical instability of the algorithm not only will increase but also a large amount of computer memory and machine time are needed in practice. As such, parameter separation must be carried out [29–31]. In this study, the orthogonal projection operator proposed by Liu et al. [31] is adopted to decompose equation (2) into the following two equations for solving velocity parameters and source parameters, respectively:(3)I−PBAδv=I−PBδt,Bδx=PBδt−Aδv,where PB, which is related to source parameters, is the orthogonal projection operator on the image space R (B) from cluster Rm to B. The coupling analysis of velocity parameters and source parameters shows that the determination of velocity disturbance is not related to the source location disturbance but only to its initial value. In the plane direction, the study area is divided into a uniform grid of 0.1°×0.1°. The velocity distribution in the model is expressed as a continuous function, and the velocity at any point in the grid is calculated by interpolation [32, 33]. The one-dimensional reference velocity model is shown in Table 1.Table 1
Crustal mean velocity model ofP wave in the target area.
Depth (km)P-wave velocity (km/s)0.05.683.05.706.05.8810.05.9515.06.0020.06.3530.06.4340.06.8055.07.80Based on the crustal upper mantle velocity model of the Sichuan-Yunnan region obtained by Liu and Wang [34] and combined with crustal velocity research results of this region in recent years [35, 36], Figure 4 is finally adopted as the crustal velocity model of this study, and detailed parameters are shown in Table 1.


### 3.2. Resolution Analysis of Solutions

In this study, the resolution of the solution is estimated using the test board method [37, 38]. The test board method is a kind of the method of seismic velocity and attenuation imaging. The basic principle of this method is that the same disturbance is applied from plus or minus direction on each node based on given speed model parameters. And then according to the actual distribution of ray, the travel time data are obtained by forward modeling calculation theory. The theoretical travel time data with a certain random error are considered as observation data, and inversion is carried out. Meanwhile, requirements in the process of the inversion method and the actual imaging method are consistent. Finally, the results from inversion and test board are compared, and their similar is regarded as the estimates of the reliability. The disturbance value is ±3% of the normal value. As 90% of the focal depth of small earthquakes is above 10 km, the resolution of solutions at different depths, such as 1 km, 3 km, 5 km, 7 km, and 9 km, is given in Figure 5. Among them, for some regions with a depth of 3 km or above (the first reservoir tail and the second reservoir head), the resolution of solution is satisfactory. Only partial nodes with a depth of 5 km still have a certain resolution, while the resolution of solution with a depth of 7 km or below is poor, which may be affected by the small sample size of the earthquake below 7 km.Figure 5
Resolution detection results ofP-wave test plate at different depths: (a) 1 km; (b) 3 km; (c) 5 km; (d) 7 km; (e) 9 km.
(a)(b)(c)(d)(e)The river reservoir seismic network is fully founded. It is a digital earthquake observation station network and has been monitoring seismic activities very well. The resolution of the solution is estimated by the test board method, which requires the depth within 9 km. Furthermore, we mainly focus on the effect of shallow earthquakes onP-wave structure. The data are enough to analyze the P-wave structures from near-field shallow earthquakes.


### 3.3. Results ofP-Wave Velocity in Study Areas

Based on the seismic observation from September 2012 to July 2015 in the two reservoir areas, the joint inversion method of source and velocity structure is used to obtain the small earthquake distribution of crustalP-wave velocity at different depths in the reservoir area (Figures 6 and 7). In this region alone, the P-wave velocity presents different images at different depths. The characteristics of P-wave velocity will be analyzed in detail for each layer in the following. At 1 km depth, a large area of low-speed zone (polygonal zone) appears along the river basin and its vicinity, especially at the first reservoir tail and the second reservoir head. The low-speed boundary of two sections with dense small earthquakes is about 13-14 km away from the river basin.Figure 6
Inversion results based on crustalP-wave velocity (a) at 1 km depth (depth of small earthquake: 0.0–1.9 km); (b) at 3 km depth (depth of small earthquake: 2.0–3.9 km); (c) at 5 km depth (depth of small earthquake: 4.0–5.9 km).
(a)(b)(c)Figure 7
Inversion results based on crustalP-wave velocity (a) at 7 km depth (depth of small earthquake: 6.0–7.9 km); (b) at 9 km depth (depth of small earthquake ≥8.0 km).
(a)(b)Shi et al. [39] and Dong and Teng [40] carried out seismic sounding and high-resolution seismic refraction observation experiments along the Yanyuan-Xichang-Mahu line and presented the high-resolution seismic refraction Pg wave travel time data by using the finite difference seismic travel time tomography algorithm. The fine structure of P-wave velocity is obtained in the crust along the boundary of the active block in western Sichuan region. The Daliangshan area, where the two reservoirs are located, is an uneven high-speed area, which indicates that the wave velocity in this area is relatively high before the impoundment.At 3 km depth, the low-speed zone is more obvious along the river basin and its vicinity, and the range is obviously expanded. The long axis of the low-speed zone extends from Wagang District in Leibo County to Longqiao in Pingshan County, showing a NE strike, about 100 km, and the boundary of the low-speed zone is slightly different from the short wheelbase of Jinsha river (Figure6(b)). In the northern section, the distance from Pingbian to Suijiang in Pingshan County is about 43 km long. In the middle section, the distance from Guihua in Leibo County to Tuanjie in Yongshan County is around 33 km long. In the southern section, the distance from Lami in Lebo County to Malan in Yongshan County is approximately 33 km long. On the west side of the water area, the low-wave velocity areas are mainly distributed in Xinzhen, Pingbian, Xining town, Guihua, Lami, and Wagang district of Pingshan County, while on the east side, the low-wave velocity areas are mainly distributed in Suijiang Pass, Leibo, Gumi, Yongshan, Tuanjie, and Malan. The spatial distribution of small earthquakes within the depth of plus or minus 1 km is similar to that of 0–2 km and is also distributed in Maoshui District of Pingshan County to Huanglang District of Leibo County and Jingxin Town to Shangtian Dam district at the junction of Yongshan County and Leibo County.At 5 km depth, the low-velocity zone gradually decreases along the Jinsha river basin and its vicinity, which is obviously weaker than the result of 0–3 km. The low-velocity zone in the water area at the first reservoir tail is mainly distributed in the areas of Anhe and Maoshui in Pingshan County, Shuanghekou in Leibo County, Suijiang County, and Guixi in Yongshan County. The low-wave velocity area in the first section of the second reservoir is mainly distributed in Wenshui, Shangtian Dam area, and Fotan area of Yongshan and Leibo Counties. The small earthquakes within the depth of plus or minus 1 km mainly occur in the low-wave velocity area (Figure6(c)). At 7 km depth, a wide range of the low-speed phenomenon gradually disappeared along the river basin and its vicinity. The large-range low-speed phenomenon in the reservoir area at 9 km basically disappears (Figure 7).To sum up, the obvious low-valueP-wave phenomenon is presented along the river basin and its adjacent areas, which is obviously different from that on the remote areas. Especially, the phenomenon of low-speed is the most obvious around 3 km, and the range of low-speed zone is the largest. With the increase in depth, the range of low-speed zone gradually decreases, and the influence of water at about 7 km depth and below gradually disappears. Based on the experimental results of water flooding in rocks made by Shi et al. [21], the low-speed bodies along the river basin and its surrounding areas may be caused by water infiltration, and the influence range and water depth may be estimated based on the low-speed P-wave region.To further analyze theP-wave velocity structure and small earthquake distribution characteristics at the regions of first reservoir tail and second reservoir head, four vertical velocity profiles are made along the long axis and the short axis of the two regions, respectively. The P-wave velocity structure and small earthquake distribution are obtained in the first area of the second reservoir. Herein, AA′ (Figure 1) is the long axis section along the Jinsha river basin in the first area of the second reservoir and BB′ (Figure 1) is the short axis section perpendicular to the Jinsha river basin. There is an obvious low-speed zone below the water area of the second reservoir. The southwest boundary of the affected area is Siguxi in Leibo County, and the affected area may extend to Jiaba reservoir in the northeast. The affected depth is mainly within 5 km, and the low-speed zone near the dam (Jincheng Town) is slightly deeper than other areas. The dam and nearby small earthquakes are mainly distributed near the interface of low and high-speed zones. However, the epicentral depth, about 17 km away from the dam, is below 5 km in low-speed zone, which indicates that there are both reservoir earthquakes and regional seismic activities in this region. So far, there is no mature theoretical approach to tell the difference between the reservoir earthquakes and regional seismic activities. Nevertheless, an analysis of focal mechanism in Xiangjiaba reservoir region is referred [24].P-wave velocity structure and small earthquake distribution are obtained in the first reservoir tail, in which CC′ (Figure 1) is the long axis profile of small earthquake in this area and DD′ (Figure 1) is the short axis profile of small earthquake. Away from Suba zone in Mabian Yi nationality autonomous county to Wanjiagou zone, P-wave velocity in this zone is higher than that in the water area, and the focal depth in small earthquakes is greater than 5 km, whose influencing zone can reach the border of Pingshan county and Mabian Yi nationality autonomous county. The focal depth of small earthquakes in the low-speed region is mainly concentrated within 5 km, and there are also small earthquakes in other regions above 5 km, which may indicate that there are both reservoir earthquakes and regional seismic activities in this region.According to the measurement of China Earthquake Networks Center, M5.3 and M5.0 earthquakes occurred in Yongshan County, Yunnan Province, on April 5 and August 17, 2014, respectively. Focal mechanism solutions of M6 ≥4.0 earthquakes are obtained based on CAP (Cut and Paste) inversion, as shown in Figure8. The CAP method is to inverse the focal mechanism of the earthquake and the magnitude of moment by waveform fitting. Then, the optimal solution of focal parameters is obtained by fitting the theoretical seismograms with the observed waveforms through grid search [41, 42]. The results show that the focal mechanism of these two earthquakes is different, and the dislocation type of Yongshan M5.3 earthquake on April 5 is the thrust type. However, the focal mechanism of the Yongshan M5.0 earthquake on August 17 was of strike slip type. Four aftershocks of the earthquake with a magnitude of 4 are in good agreement with main earthquake. Whether the difference in focal dislocation type between two earthquakes is related to reservoir water storage is worth for further study.Figure 8
Focal mechanism distribution of M5.3 and 5.0 earthquakes and M4 aftershocks in Yongshan.According to the error change of hypocenter mechanism inversion with depth, the focal depth of Yongshan M5.3 earthquake on April 5, 2014, is relatively shallow, around 2 km. Combined with theP-wave velocity along the long axis section of second reservoir given in Figure 9, this earthquake occurs in the low-wave velocity region. On August 17, 2014, the focal depth of Yongshan M5.0 earthquake is 6 km, a relatively deep depth. The earthquake occurs at the intersection of high- and low-speed zones.Figure 9
Change of focal mechanism with depth for two earthquakes: (a) Yongshan M5.3 earthquake, April 5, 2014; (b) Yongshan M5.0 earthquake, August 7, 2014.
(a)(b)


## 3.1. Inversion Method

Seismic tomography is a geophysical method which uses inversion analysis to deduce the important information such as velocity structure of underground medium and other physical property parameters by analyzing the kinematics (such as travel time and ray path) and dynamics (such as waveform and amplitude) characteristics of observed seismic wave phases.In the parametric representation of 3D velocity structure, grid method is a common method. Its advantage is that the grids can be divided into different spaces according to the density of the distribution of the epicenter and the station, to ensure sufficient ray cross coverage in each grid and reduce the impact of the prior assumptions. In this method, the velocity structure is expressed as a continuous function, the unknown parameter is the velocity value on the 3D mesh node, and the velocity at any point in the model is calculated by interpolation.In the joint inversion of source location and velocity structure, the travel time residualδt is caused by the perturbation of source parameters and velocity. According to the relevant literature [26–28], the problem can be expressed by the following linear equation:(1)δt=Δt+∂t∂xΔx+∂t∂yΔy+∂t∂zΔz+∑n=1N∂t∂vΔvn.Δt, Δx, Δy, and Δz and Δvn represent the source of original times, longitude, latitude, depth of disturbance, and speed disturbance, respectively; N is the total number of velocity parameters. For l earthquakes and K stations, equation (1) can be written in the compact form as follows:(2)δt=Aδv+Bδx,where δt is the m-dimensional travel time residual vector; δv is the n-dimensional node velocity perturbation vector; δx is the 4l-dimensional source parameter perturbation vector; A is the m × n-dimensional travel time partial derivative matrix with respect to velocity; and B is the m × 4l-dimensional travel time partial derivative matrix with respect to source parameters.According to the basic equation (2) of joint inversion, velocity parameters and source parameters are coupled to each other. For inverting two parameters of different dimensions in the same equation at the same time, the numerical instability of the algorithm not only will increase but also a large amount of computer memory and machine time are needed in practice. As such, parameter separation must be carried out [29–31]. In this study, the orthogonal projection operator proposed by Liu et al. [31] is adopted to decompose equation (2) into the following two equations for solving velocity parameters and source parameters, respectively:(3)I−PBAδv=I−PBδt,Bδx=PBδt−Aδv,where PB, which is related to source parameters, is the orthogonal projection operator on the image space R (B) from cluster Rm to B. The coupling analysis of velocity parameters and source parameters shows that the determination of velocity disturbance is not related to the source location disturbance but only to its initial value. In the plane direction, the study area is divided into a uniform grid of 0.1°×0.1°. The velocity distribution in the model is expressed as a continuous function, and the velocity at any point in the grid is calculated by interpolation [32, 33]. The one-dimensional reference velocity model is shown in Table 1.Table 1
Crustal mean velocity model ofP wave in the target area.
Depth (km)P-wave velocity (km/s)0.05.683.05.706.05.8810.05.9515.06.0020.06.3530.06.4340.06.8055.07.80Based on the crustal upper mantle velocity model of the Sichuan-Yunnan region obtained by Liu and Wang [34] and combined with crustal velocity research results of this region in recent years [35, 36], Figure 4 is finally adopted as the crustal velocity model of this study, and detailed parameters are shown in Table 1.


## 3.2. Resolution Analysis of Solutions

In this study, the resolution of the solution is estimated using the test board method [37, 38]. The test board method is a kind of the method of seismic velocity and attenuation imaging. The basic principle of this method is that the same disturbance is applied from plus or minus direction on each node based on given speed model parameters. And then according to the actual distribution of ray, the travel time data are obtained by forward modeling calculation theory. The theoretical travel time data with a certain random error are considered as observation data, and inversion is carried out. Meanwhile, requirements in the process of the inversion method and the actual imaging method are consistent. Finally, the results from inversion and test board are compared, and their similar is regarded as the estimates of the reliability. The disturbance value is ±3% of the normal value. As 90% of the focal depth of small earthquakes is above 10 km, the resolution of solutions at different depths, such as 1 km, 3 km, 5 km, 7 km, and 9 km, is given in Figure 5. Among them, for some regions with a depth of 3 km or above (the first reservoir tail and the second reservoir head), the resolution of solution is satisfactory. Only partial nodes with a depth of 5 km still have a certain resolution, while the resolution of solution with a depth of 7 km or below is poor, which may be affected by the small sample size of the earthquake below 7 km.Figure 5
Resolution detection results ofP-wave test plate at different depths: (a) 1 km; (b) 3 km; (c) 5 km; (d) 7 km; (e) 9 km.
(a)(b)(c)(d)(e)The river reservoir seismic network is fully founded. It is a digital earthquake observation station network and has been monitoring seismic activities very well. The resolution of the solution is estimated by the test board method, which requires the depth within 9 km. Furthermore, we mainly focus on the effect of shallow earthquakes onP-wave structure. The data are enough to analyze the P-wave structures from near-field shallow earthquakes.


## 3.3. Results ofP-Wave Velocity in Study Areas

Based on the seismic observation from September 2012 to July 2015 in the two reservoir areas, the joint inversion method of source and velocity structure is used to obtain the small earthquake distribution of crustalP-wave velocity at different depths in the reservoir area (Figures 6 and 7). In this region alone, the P-wave velocity presents different images at different depths. The characteristics of P-wave velocity will be analyzed in detail for each layer in the following. At 1 km depth, a large area of low-speed zone (polygonal zone) appears along the river basin and its vicinity, especially at the first reservoir tail and the second reservoir head. The low-speed boundary of two sections with dense small earthquakes is about 13-14 km away from the river basin.Figure 6
Inversion results based on crustalP-wave velocity (a) at 1 km depth (depth of small earthquake: 0.0–1.9 km); (b) at 3 km depth (depth of small earthquake: 2.0–3.9 km); (c) at 5 km depth (depth of small earthquake: 4.0–5.9 km).
(a)(b)(c)Figure 7
Inversion results based on crustalP-wave velocity (a) at 7 km depth (depth of small earthquake: 6.0–7.9 km); (b) at 9 km depth (depth of small earthquake ≥8.0 km).
(a)(b)Shi et al. [39] and Dong and Teng [40] carried out seismic sounding and high-resolution seismic refraction observation experiments along the Yanyuan-Xichang-Mahu line and presented the high-resolution seismic refraction Pg wave travel time data by using the finite difference seismic travel time tomography algorithm. The fine structure of P-wave velocity is obtained in the crust along the boundary of the active block in western Sichuan region. The Daliangshan area, where the two reservoirs are located, is an uneven high-speed area, which indicates that the wave velocity in this area is relatively high before the impoundment.At 3 km depth, the low-speed zone is more obvious along the river basin and its vicinity, and the range is obviously expanded. The long axis of the low-speed zone extends from Wagang District in Leibo County to Longqiao in Pingshan County, showing a NE strike, about 100 km, and the boundary of the low-speed zone is slightly different from the short wheelbase of Jinsha river (Figure6(b)). In the northern section, the distance from Pingbian to Suijiang in Pingshan County is about 43 km long. In the middle section, the distance from Guihua in Leibo County to Tuanjie in Yongshan County is around 33 km long. In the southern section, the distance from Lami in Lebo County to Malan in Yongshan County is approximately 33 km long. On the west side of the water area, the low-wave velocity areas are mainly distributed in Xinzhen, Pingbian, Xining town, Guihua, Lami, and Wagang district of Pingshan County, while on the east side, the low-wave velocity areas are mainly distributed in Suijiang Pass, Leibo, Gumi, Yongshan, Tuanjie, and Malan. The spatial distribution of small earthquakes within the depth of plus or minus 1 km is similar to that of 0–2 km and is also distributed in Maoshui District of Pingshan County to Huanglang District of Leibo County and Jingxin Town to Shangtian Dam district at the junction of Yongshan County and Leibo County.At 5 km depth, the low-velocity zone gradually decreases along the Jinsha river basin and its vicinity, which is obviously weaker than the result of 0–3 km. The low-velocity zone in the water area at the first reservoir tail is mainly distributed in the areas of Anhe and Maoshui in Pingshan County, Shuanghekou in Leibo County, Suijiang County, and Guixi in Yongshan County. The low-wave velocity area in the first section of the second reservoir is mainly distributed in Wenshui, Shangtian Dam area, and Fotan area of Yongshan and Leibo Counties. The small earthquakes within the depth of plus or minus 1 km mainly occur in the low-wave velocity area (Figure6(c)). At 7 km depth, a wide range of the low-speed phenomenon gradually disappeared along the river basin and its vicinity. The large-range low-speed phenomenon in the reservoir area at 9 km basically disappears (Figure 7).To sum up, the obvious low-valueP-wave phenomenon is presented along the river basin and its adjacent areas, which is obviously different from that on the remote areas. Especially, the phenomenon of low-speed is the most obvious around 3 km, and the range of low-speed zone is the largest. With the increase in depth, the range of low-speed zone gradually decreases, and the influence of water at about 7 km depth and below gradually disappears. Based on the experimental results of water flooding in rocks made by Shi et al. [21], the low-speed bodies along the river basin and its surrounding areas may be caused by water infiltration, and the influence range and water depth may be estimated based on the low-speed P-wave region.To further analyze theP-wave velocity structure and small earthquake distribution characteristics at the regions of first reservoir tail and second reservoir head, four vertical velocity profiles are made along the long axis and the short axis of the two regions, respectively. The P-wave velocity structure and small earthquake distribution are obtained in the first area of the second reservoir. Herein, AA′ (Figure 1) is the long axis section along the Jinsha river basin in the first area of the second reservoir and BB′ (Figure 1) is the short axis section perpendicular to the Jinsha river basin. There is an obvious low-speed zone below the water area of the second reservoir. The southwest boundary of the affected area is Siguxi in Leibo County, and the affected area may extend to Jiaba reservoir in the northeast. The affected depth is mainly within 5 km, and the low-speed zone near the dam (Jincheng Town) is slightly deeper than other areas. The dam and nearby small earthquakes are mainly distributed near the interface of low and high-speed zones. However, the epicentral depth, about 17 km away from the dam, is below 5 km in low-speed zone, which indicates that there are both reservoir earthquakes and regional seismic activities in this region. So far, there is no mature theoretical approach to tell the difference between the reservoir earthquakes and regional seismic activities. Nevertheless, an analysis of focal mechanism in Xiangjiaba reservoir region is referred [24].P-wave velocity structure and small earthquake distribution are obtained in the first reservoir tail, in which CC′ (Figure 1) is the long axis profile of small earthquake in this area and DD′ (Figure 1) is the short axis profile of small earthquake. Away from Suba zone in Mabian Yi nationality autonomous county to Wanjiagou zone, P-wave velocity in this zone is higher than that in the water area, and the focal depth in small earthquakes is greater than 5 km, whose influencing zone can reach the border of Pingshan county and Mabian Yi nationality autonomous county. The focal depth of small earthquakes in the low-speed region is mainly concentrated within 5 km, and there are also small earthquakes in other regions above 5 km, which may indicate that there are both reservoir earthquakes and regional seismic activities in this region.According to the measurement of China Earthquake Networks Center, M5.3 and M5.0 earthquakes occurred in Yongshan County, Yunnan Province, on April 5 and August 17, 2014, respectively. Focal mechanism solutions of M6 ≥4.0 earthquakes are obtained based on CAP (Cut and Paste) inversion, as shown in Figure8. The CAP method is to inverse the focal mechanism of the earthquake and the magnitude of moment by waveform fitting. Then, the optimal solution of focal parameters is obtained by fitting the theoretical seismograms with the observed waveforms through grid search [41, 42]. The results show that the focal mechanism of these two earthquakes is different, and the dislocation type of Yongshan M5.3 earthquake on April 5 is the thrust type. However, the focal mechanism of the Yongshan M5.0 earthquake on August 17 was of strike slip type. Four aftershocks of the earthquake with a magnitude of 4 are in good agreement with main earthquake. Whether the difference in focal dislocation type between two earthquakes is related to reservoir water storage is worth for further study.Figure 8
Focal mechanism distribution of M5.3 and 5.0 earthquakes and M4 aftershocks in Yongshan.According to the error change of hypocenter mechanism inversion with depth, the focal depth of Yongshan M5.3 earthquake on April 5, 2014, is relatively shallow, around 2 km. Combined with theP-wave velocity along the long axis section of second reservoir given in Figure 9, this earthquake occurs in the low-wave velocity region. On August 17, 2014, the focal depth of Yongshan M5.0 earthquake is 6 km, a relatively deep depth. The earthquake occurs at the intersection of high- and low-speed zones.Figure 9
Change of focal mechanism with depth for two earthquakes: (a) Yongshan M5.3 earthquake, April 5, 2014; (b) Yongshan M5.0 earthquake, August 7, 2014.
(a)(b)


## 4. Conclusions

Based on the seismic observation report recorded by the reservoir seismic network from September 2012 to July 2015, this study uses the joint inversion method of source and velocity structure to obtainP-wave velocity structure for different sections in the reservoir area. The impact range of reservoir water storage is analyzed, and the impact of reservoir water infiltration on crustal medium is discussed.(1)
The earthquake location results show that the earthquake activities are mainly concentrated in two areas, namely, the first reservoir tail and the second reservoir head(2)
The focal depths of small earthquakes are mainly distributed within 10 km, among which the focal depths of small earthquakes in the first reservoir tail and the second reservoir head are mainly within 5 km(3)
Small earthquakes above 5 km are mainly distributed in the area far from the reservoir water area, which is consistent with both reservoir earthquakes and regional seismic activities in this areaThe resolution results of detection board at different depths show that the solution resolution in the study area within 5 km is satisfactory. Affected by the impaction of reservoir water storage, the low-speed phenomenon appears along the Jinsha river basin and its surrounding areas. Especially at 3 km, the low-speed phenomenon is the most obvious and the low-speed zone has the largest range. With the increase in depth, the low-speed zone gradually reduces, and the influence of water at 7 km depth and below gradually disappears. According to the distribution range of low-speed bodies, the influence range and depth of reservoir water storage on medium velocity in the surrounding area can be speculated. At present, the maximum depth of seepage in the first reservoir may be 5 km.In fact, observed data in this study are from September 2012 to July 2015. So far, a lot of new recorded data are obtained and analyzed in future. Nevertheless, the observed data from 2012 to 2015 contain the whole process from initial to full storage typically.



---

*Source: 6250086-2021-06-30.xml*