Study of 3-D velocity structure characteristics in Dadu river Houziyan reservoir area at different impoundment stages

1 Introduction

Since the initiation of ‘The Development of the Western Region in China’ and ‘Electricity Transmission from West to East’ more than 2 decades ago, numerous reservoirs have been constructed in Western China. As China advances swiftly towards its goals of ‘Peak Carbon Dioxide Emissions’ and ‘Carbon Neutrality,’ the hydropower industry, a vital component of energy transformation, is expected to flourish (Fang et al., 2010; Li et al., 2022; Liu et al., 2022; Ma and Nie, 2022). However, the geological conditions in Western China are highly intricate, characterized by the presence of complex geological structures at high altitudes. Consequently, artificial projects such as reservoir dam construction and reservoir impoundment have elevated the risk of seismic activity (Su et al., 2020; Dou et al., 2022; Guo et al., 2022; Liu et al., 2023). Reservoir-induced seismicity (RIS) is a distinctive form of seismic activity that occurs within a defined timeframe in or around reservoir areas, primarily triggered by the processes of reservoir impoundment and drainage (Li et al., 2004; Rinaldi et al., 2020). To date, approximately 130 confirmed cases of RIS have been reported worldwide (Van Eijs et al., 2006; Basbous et al., 2022; Zheng et al., 2022). Despite most RIS events having low magnitudes, they can be highly destructive, causing damage to nearby dams and structures and even giving rise to secondary disasters like landslides and collapses, posing a threat to downstream safety (Che et al., 2009; Zhang et al., 2019). Consequently, RIS plays a crucial role not only in hydraulic and hydroelectric engineering research but also in the realms of seismology, regional structural stability, and environmental engineering geology research (Gaucher et al., 2015; Chang et al., 2022; Chen, 2022).

The impounded water level is a pivotal factor in triggering RIS events. Scholars widely concur that the likelihood of an earthquake induced by a high dam and a large storage reservoir is higher. Numerous studies have demonstrated a significant increase in the frequency of small and medium seismic activities around reservoirs following impoundment. Changes in reservoir water levels often serve as triggers for earthquakes, with sharp fluctuations in water levels capable of instigating major seismic events. Additionally, both the frequency and magnitude of reservoir-induced earthquakes tend to peak as the reservoir approaches its maximum water level (Liu and Zhang, 2016; Lizurek et al., 2021; Li and Luo, 2022). Research by Simpson has shown that an increase in nearby seismic activity is likely if the water level of a reservoir surpasses its original maximum level. Moreover, the occurrence of earthquakes is closely linked to the rate of change in the reservoir water level: seismic activity decreases when the rate is stable but surges when the rate changes (Simpson and Negmatullaev, 1981). Generally, seismic activity increases with rising water levels and decreases with falling water levels. However, in some cases, water levels exhibit a negative correlation with seismic activity, with induced earthquakes occurring after drainage. Prior to significant RIS events, small seismic enhancements, akin to ‘foreshocks,’ are often observed. Increased seismic activity typically follows a peak in water levels or a sudden change in the rate of water injection into the reservoir. Based on the timing of impoundment and earthquake occurrence, reservoir-induced earthquakes can be categorized as ‘fast response type earthquakes’ or ‘delayed response type earthquakes’. The elastic stress resulting from the water load leads to pore compression in the bottom rocks of the reservoir, causing an increase in pore pressure and a reduction in rock strength, which in turn triggers fast-response type earthquakes. In the case of delayed response types, their activity and lag time are associated with the diffusion of pore water pressure (Simpson et al., 1988). Research has indicated that the mechanism of reservoir-induced seismicity comprises reservoir load and pore pressure diffusion categories, with reservoir load mechanisms typically leading to fast responses and pore pressure diffusion mechanisms tending to induce delayed responses (Kinscher et al., 2023). A variety of factors, including the duration of the impoundment process and variations in geological configurations, can influence the lag time of each reservoir-induced earthquake (Kim et al., 2022).

Research into RIS events has shown that changes in underground media are reflected by tracking the wave velocity ratio in a reservoir area. This phenomenon can be attributed to three key factors: in addition to the aforementioned effects of fluid loading and pore pressure diffusion, another key factor is that the filling of the original ‘dry’ fractured faults with water reduces the frictional resistance of fault activity. Therefore, even small stress disturbances can easily cause rock fracturing. When seismic waves traverse the source area of a future significant earthquake, the transmission speed may fluctuate by a variable percentage, typically ranging from a few to more than ten percent. Changes in the stress state within the source area lead to alterations in the elastic modulus, causing variations in wave velocity (Wang et al., 2021; Bashir et al., 2022; Jin, 2022). Consequently, the study of the three-dimensional velocity structure, wave velocity ratio, and attenuation structure in a reservoir area can provide valuable insights into the subsurface fluid conditions (Xiong et al., 2013; Ersoy et al., 2019).

Seismic tomography is not only a crucial method for detecting transverse inhomogeneous structures in media at various scales within the Earth but also a major technology for constraining the distribution of geophysical parameters, including elastic, non-elastic, heterogeneity, and anisotropy parameters, which influence seismic wave propagation (Rawlinson et al., 2010; Zhao et al., 2016; Tan et al., 2022). Since the mid-1970s, seismic tomography has advanced rapidly and became an indispensable tool in modern seismology (Zhang and Thurber, 2006). In recent years, with increased attention to reservoir earthquakes and the establishment of reservoir seismic networks, seismic tomography has been widely utilized to invert the velocity structure of reservoir areas. Many scholars have gained a more accurate understanding of the fine velocity structures of the upper crust of reservoirs by using seismic tomography to analyze the velocity structure characteristics of seismic waves in a reservoir area, study the impact of reservoir water on crustal structure, and investigate the relationship between seismic activity and the distribution of velocity anomalies. In this paper, we use the Fast Marching Tomography Package (FMTOMO) developed by Rawlinson to invert the three-dimensional velocity structure of a reservoir area (Rawlinson and Sambridge, 2003; Rawlinson and Sambridge, 2004a; Rawlinson and Sambridge, 2004b). FMTOMO addresses the forward problem of travel time prediction through a multi-stage fast marching method based on the Fast Marching Method (FMM) and adjusts model parameters using a subspace inversion scheme to align with data observations. Employing FMTOMO for velocity structure inversion offers several advantages, including the ability to invert and track the travel time of various seismic phases, effectively simulate complex velocity and interface structures, perform velocity structure tomography calculations with combined seismic data from nearby, regional, and remote earthquakes, and simultaneously or separately invert various types of velocity structures, layer interfaces, or source locations.

The Dadu River Houziyan Reservoir Seismic Network in Kangding County, Ganzi Prefecture, Sichuan Province, commenced operations in October 2015. Up to November 2021, we have collected a substantial amount of seismic data from monitoring 6330 seismic events using this network. In this paper, we have used FMTOMO tomography to invert the underground velocity structure of this reservoir at different impoundment stages, enabling us to analyze the temporal and spatial characteristics of the velocity structure in the reservoir area at different impoundment stages, changes in the underground velocity structure due to reservoir impoundment, and the overall impact of impoundment. Furthermore, the data we have collected can provide valuable support for future research into the mechanisms of seismic activity in the Houziyan reservoir area, the influence of reservoir water infiltration on seismic activities and medium properties in the reservoir area, as well as earthquake prevention and disaster reduction efforts, among other aspects.

2 Data and methods

2.1 Structural setting and seismic activities in the reservoir area

The Dadu River Houziyan Hydropower Station was located in an alpine valley region in Kangding County, Sichuan Province, at the southeast tip of the Tibetan Plateau, the southernmost tip of the strong uplift area in northwest Sichuan, with alpines and rivers as the main landform (Tang and Yong, 2012). The Wenchuan-Miaoxian fault zone, trending in a northeast direction, and the Xianshuihe fault and Yuke fault zones, trending in a northwest direction, are developed around the reservoir area. Among them, the seismic activity of the Xianshuihe fault zone is the most intense. The Yuke fault is closer to the dam site. The regional geological structure and seismic geological setting of this reservoir area are very complicated. Seismic hazard analysis shows that this reservoir area is greatly affected by the seismic activities in the Xianshuihe fault zone. The area near the dam site is the main part of the Northwest Sichuan Plateau, which has been lifted to become a component of the Tibetan Plateau since the Quaternary. The faults near the dam site are relatively small in scale, and they are all early to mid-Pleistocene faults, with relatively weak activity and good preservation of the plateau surface. Also, the seismicity near the dam site area is relatively weak, with the maximum magnitude at M_S5.7 (which occurred on 30 September 1972, about 19 km away from the site) and a few earthquakes at M_S2.0.

There are 8 field observation stations in the Houziyan Reservoir Seismic Network, consisting of short-period seismometers (FSS-3M), seismic data collectors (EDAS-24GN), and other monitoring equipment. The completeness magnitude in the range of the reservoir area is M_L0.5, which conforms to the data precision required for this research. The network was put into operation in October 2015. As of November 2021, 6330 earthquakes were monitored within the monitoring range (101.5°–102.5° east longitude, 30.2°–31.2° north latitude), which includes 3 earthquakes with magnitude M_L4.0 or higher. See Table 1 for detailed earthquake magnitude statistics. The seismic location method in Sichuan-Yunnan region (Loc3dSB), developed by science and technology special project of earthquake industry ‘Compiling earthquake travel time table in Sichuan-Yunnan Region’, was used to locate the 6330 seismic events in the Houziyan Network, and The Chinese Sichuan-Yunnan three-dimensional velocity model is used for positioning velocity model. There is a total of 27799 Pg travel time data and 23961 Sg travel time data from 6330 earthquakes.

TABLE 1. The statistical table of earthquake magnitudes in the houziyan network (M_L).

We have obtained the distribution diagram of epicenters (Figure 1) and the diagram of M-T of reservoir area earthquakes and water level (Figure 2) based on the monitoring data from the Houziyan Reservoir Seismic Network and data on water levels in the reservoir area. The diagrams reveal the following observations: 1) Seismic activity is relatively weak before impoundment. The number of earthquakes significantly increases during the initial impoundment stage, with most of them being weak earthquakes with a magnitude of M_L2.0, and only a few are felt earthquakes; 2) During the initial impoundment stage, there is a clear correlation between induced seismic activities increases with the increase of reservoir water levels. However, as the impoundment time lengthens, seismic activities in the reservoir area appear to be influenced by the disturbance of the regional stress field; 3) The epicenters of the earthquakes are generally widely distributed around the reservoir area, particularly near the Xianshuihe fault in the southwest corner of the monitoring area, to the south of the monitoring range, and other areas. Based on the regional geological structure, our analysis indicates the presence of several relatively dense seismic regions (zones) in this area: There is an inclined seismic concentration region (zone) along the Xianshuihe fault, where 15 earthquakes with magnitudes of M_L3.0∼3.9 were located, accounting for 60% of earthquakes in the entire monitoring area in that magnitude range. Three earthquakes above M_L4.0 were also recorded, representing 100% of earthquakes at M_L4.0 and higher. There is an earthquake region (zone) in the upstream area of the dam, which is also a seismic concentration region with magnitudes below M_L1.0. Currently, two events have been monitored at M_L3.0∼3.9. In the northeast of the dam, near the monitoring boundary area, there is an earthquake region with a few earthquakes, mainly within M_L3.0. Currently, one event has been monitored with M_L3.0∼3.9. Just below the dam, near the monitoring boundary area, there is an earthquake region with earthquakes mainly within M_L3.0. Currently, two events have been monitored with M_L3.0∼3.9. There is an earthquake region in the NW-SE direction between the Xianshuihe fault and Yuke fault and a horizontal earthquake region in the downstream area of the dam. Both regions have dispersed earthquakes, mainly within M_L3.0. Currently, no event has been monitored with M_L3.0 and higher.

FIGURE 1. The earthquake epicenter distribution in the Houziyan network. The dashed lines represent the latitude and longitude grid.

FIGURE 2. The M-T diagram in the reservoir area and water level.

2.2 Data processing

Using the seismic travel time data recorded by the Houziyan Network from October 2015 to November 2021, the three-dimensional velocity structures of V_P, V_S, and V_P/V_S in this region were inverted. In order to keep a big enough data volume for the research and guarantee the precision of imaging inversion, according to the monitoring capabilities, station distribution, and waveform recording of the Houziyan Seismic Network, some obviously erroneous travel time data were first removed. Then, earthquakes with M_L0.5 or higher recorded by at least four stations simultaneously were selected for inversion. After the data were screened and processed, the inversion ultimately utilized 14,093 Pg travel time data from 2,627 earthquakes and 10,681 Sg travel time data from 2,045 earthquakes. The time distance curves of the Pg and Sg phases are shown in Figure 3. According to the figure, the errors of Pg and Sg travel time dates are respectively within 0.5s and 1.0s.

FIGURE 3. Pg and Sg travel time data for inversion.

2.3 Impoundment stage division

To analyze changes of the underground velocity structure at different impoundment stages, the following five stages of the Houziyan reservoir area impoundment period are divided according to the impoundment time, water level cycle, and seismicity (Figure 2; Table 2): before impoundment (10/2015–11/2016), first impoundment stage (12/2016–12/2017), second impoundment stage (01/2018–12/2018), third impoundment stage (01/2019–12/2019), and fourth impoundment stage (01/2020–11/2021). See Figure 2 for details of impoundment stage division and changes in water level. Table 3 shows the number of filtered seismic events and Pg and Sg travel time data used for inversion at different impoundment stages.

TABLE 2. The statistical table of impoundment stage division and changes of water levels for the dadu river houziyan reservoir area.

TABLE 3. The statistical table of seismic data volume used for inversion at different impoundment stages.

2.4 FMTOMO tomography

The Fast Marching Tomography Package (FMTOMO) is a tomography software package developed by Rawlinson and Sambridge in 2004 (Rawlinson and Sambridge, 2004a; Rawlinson and Sambridge, 2004b). FMTOMO utilizes the Fast Marching Method (FMM) for predicting travel times in the forward step and employs a subspace inversion scheme to adjust model parameters to fit observed data. In solving the forward problem, FMTOMO relies on the multi-stage Fast Marching Method based on the FMM. The FMM is a grid-based eikonal solver that implicitly tracks the evolution of wavefronts in 3-D layered media. To solve the inverse problem, it requires the adjustment of model parameters in order to satisfy the data, subject to regularization constrains. A common problem with tomographic inversion is that not all model parameters will be well constrained by the data alone. To help combat this problem, several regularization terms are added to the objective function. The damping parameter and smoothing parameter control these regularization terms and govern the trade-off between how well the solution satisfies the data, the proximity of the solution model to the starting model, and the smoothness of the solution model (Rawlinson et al., 2008).

FMTOMO can simultaneously invert one or more types of travel time datasets for variations in wave speed, interface depth, and source location. The initial model parameterization in the FMTOMO program consists of three types of model grids: a velocity grid, an interface grid, and a propagation grid. Interface are represented by a regular grid of points in latitude and longitude, the depths of which can vary, with cubic B-splines used to define a continuous surface. Propagation grid represents a discrete sampling of the velocity fields for implicitly tracking seismic wavefronts to solve the forward problem. To ensure the resolution and precision of inversion results, various model parameterizations are tested before inversion using checkerboard detection plates. These tests help to determine the model parameterization to be used for the final inversion. For the Houziyan network, the following model parameterizations are chosen based on tests involving travel time data, station distribution, and source locations: a velocity grid of 10 km × 10 km × 2.5 km; a propagation grid of 2 km × 2 km × 0.5 km; and an interface grid of 10 km × 10 km (horizontal node spacing 10 km).

The FMTOMO program introduces damping and smoothing coefficients during the inversion process to control the convergence rate of the model solution and to constrain the smoothness of the results. A trade-off analysis of damping coefficients and smoothing coefficients is performed using a compromise curve (see Figure 4) to select suitable values for these coefficients in the inversion process. The analysis indicates that the model is relatively smooth, and the mean square error of relative travel time is small when the damping coefficient is set to 0.5 and the smoothing coefficient to 1.0. Therefore, these coefficients are chosen in the inversion.

FIGURE 4. The compromise curve.

Reasonable selection of the initial velocity model helps to obtain reliable inversion results. The 1-D P -velocity initial velocity models used as the input model of the 3-D P-wave tomography are mainly the crustal velocity models from experts such as Zhaozhu (Zhao and Zhang, 1987) and are finely adjusted according to trial results of travel time data and research outcomes of other scholars (Zhao et al., 1997; Wang et al., 2002). The initial V_P/V_S was obtained by fitting the travel-time relation between P waves and S waves. The slope of the fitting line (Figure 5) shows that the average V_P/V_S of the study area is 1.69. The 1-D S-wave velocity model used as the input model of the 3-D S-wave tomography comes from the P-wave velocity model divided by 1.69. The 1-D initial velocity model was illustrated in Figure 6.

FIGURE 5. The Fitting line of the Ratio of V_P/V_S Calculated from the Travel Times of P and S Waves.

FIGURE 6. The 1-D initial velocity model for inversion.

3 Results

Through FMTOMO tomography inversion, we got the 3-D velocity structure, including values of V_P, V_S, and V_P/V_S ratios at different impoundment stages of the Houziyan reservoir area. To better reflect the characteristics of the velocity structure in the reservoir area, the horizontal section result maps including different depths of 1 km, 2 km, 3 km, 4 km, 5 km, 6 km, 8 km, and 10 km and the vertical profile result maps including different directions of 30.50°N, 30.55°N, 30.60°N, 30.65°N, 30.70°N, 30.75°N, 101.95°E, 102.00°E, 102.05°E, 102.10°E, 102.15°E, 102.20°E, WN-ES have been derived, as shown in Figure 7. At each stage, we chose the horizontal section maps with depths of 3 km and 6 km as well as the profile result maps in the directions of 30.65°N and WN-ES to show different wave velocities. All results were referenced in the actual analysis and interpretation.

FIGURE 7. The cross-sections profile projections in the horizontal planes.

3.1 Reliability tests

In order to ensure data volume and the reliability of velocity structure tomography results, we examined the distribution of rays in horizontal sections and vertical sections at different water storage stages. The horizontal ray distribution of the P-wave are shown in Figure 8. Meanwhile, we conducted a checkerboard test to assess the degree to which the initial checkerboard pattern could be restored using the data. This allowed us to evaluate the resolution of the data, the impact of model parameterization, and the suitability of damping and smoothing parameters. Firstly, perturbations were added to the initial velocity model and then the perturbed model was used to calculate the theoretic travel times. Subsequently, the synthetic travel times were taken as input data to invert the velocity model when the hypocentral parameters were kept constant. The inversion results after the checkerboard test are shown in Figure 9 and Figure 10. The resolution of the V_S tomography is similar to that of V_P tomography because of the similar ray coverage and the same model parameterizations and the same damping and smoothing parameter. Thus, the checkerboard test results and ray distribution of V_S are not presented here.

FIGURE 8. The Horizontal Ray Distribution of P-wave in different water storage stages. (A–E) represent before impoundment, impoundment stage 1, impoundment stage 2, impoundment stage 3, and impoundment stage 4 respectively.

FIGURE 9. The Horizontal Results of Checkerboard test of P-wave in different Water Storage Stages. (A–J) represent before impoundment in 3 km, before impoundment in 6 km, impoundment stage 1 in 3 km, impoundment stage 1 in 6 km, impoundment stage 2 in 3 km, impoundment stage 2 in 6 km, impoundment stage 3 in 3 km, impoundment stage 3 in 6 km, impoundment stage 4 in 3 km, and impoundment stage 4 in 6 km, respectively.

FIGURE 10. The Vertical Results of Checkerboard test of P-wave in different Water Storage Stages. (A–J) represent before impoundment along longitude, before impoundment along latitude, impoundment stage 1 along longitude, impoundment stage 1 along latitude, impoundment stage 2 along longitude, impoundment stage 2 along latitude, impoundment stage 3 along longitude, impoundment stage 3 along latitude, impoundment stage 4 along longitude, and impoundment stage 4 along latitude, respectively.

The results of these tests indicate that, before impoundment, the ray distribution was relatively sparse, resulting in low detection plate resolution. However, after impoundment, the ray distribution density significantly increased, leading to improved resolution. The selected data exhibited high precision and provided dependable results. Overall, the imaging results within the depth range of 1–12 km show good resolution, effectively recovering the velocity model (Figure 10). We notice that although the test result can recover the pattern of the velocity structure at 12 km depth, the resolution decreases beyond 12 km depth. The results of Checkerboard test are consistent with the distribution of earthquake depths. In addition, we can notice that there is a loss of resolution close to the boundaries of the investigated area. Therefore, this article only discusses regions with good resolution.

3.2 Imaging results of V_P

According to V_P in the section result maps (Figure 11) and profile result maps (Figure 12), there is local lateral heterogeneity before impoundment. The manifestations include that a high wave velocity zone of 5700–6100 m/s appears from the Houziyan dam to the upstream reservoir area and in a 4 km depth area from the surface to underground; overall, the P-wave velocity increases as the depth becomes deeper, with a value of about 5300 m/s at the shallow area of 1 km to about 6100 m/s at the deep area of 10 km.

FIGURE 11. The horizontal section maps of V_P at different impoundment stages. (A–J) represent before impoundment in 3 km, before impoundment in 6 km, impoundment stage 1 in 3 km, impoundment stage 1 in 6 km, impoundment stage 2 in 3 km, impoundment stage 2 in 6 km, impoundment stage 3 in 3 km, impoundment stage 3 in 6 km, impoundment stage 4 in 3 km, and impoundment stage 4 in 6 km, respectively.

FIGURE 12. The vertical profile maps of V_P at different impoundment stages. (A–J) represent before impoundment along longitude, before impoundment along distance, impoundment stage 1 along longitude, impoundment stage 1 along distance, impoundment stage 2 along longitude, impoundment stage 2 along distance, impoundment stage 3 along longitude, impoundment stage 3 along distance, impoundment stage 4 along longitude, and impoundment stage 4 along distance, respectively.

There are local lateral changes in V_P from stages 1–4 after impoundment. At stage 1 after impoundment, a high wave velocity zone of 5800–6400 m/s appears from the Houziyan dam to the upstream reservoir area and in a 5 km depth area from the surface to underground; at stage 2, a low wave velocity zone of 5100–5600 m/s appears from the dam to the upstream reservoir area and in a 5 km depth area from the surface to underground. Compared with the cases before and at first stage of impoundment, the original high value zone becomes a low value zone with a variation of 600–800 m/s, which indicates that impoundment has resulted in changes of V_P of the underground media and there is dynamic adjustment of underground wave velocity. At stage 3, a high wave velocity zone of 5800–6400 m/s appears from the Houziyan dam to the upstream reservoir area and in a 6 km depth area from the surface to underground; at stage 4, a high wave velocity zone of 6200–6500 m/s appears from the Houziyan dam to the upstream reservoir area and in a 6 km depth area from the surface to underground.

Throughout all impoundment stages except stage 2, V_P changes from 5300 m/s at the shallow area of 1 km–6500 m/s at the deep area of 10 km; compared with the case before impoundment, the position of the high wave velocity zone does not change greatly, with the velocity increased by 200–500 m/s. Compared with the cases before and at first stage of impoundment, the original high value zone of V_P at stage 2 becomes to a low one with a variation of 300–600 m/s, which indicates that impoundment has resulted in changes of V_P of the underground media.

The distribution characteristics of P-wave velocity from different impoundment stages demonstrate the impact of reservoir water on the velocity structure of the upper crust in the reservoir area. Before impoundment, the velocity of the reservoir area to the northeast (upstream of the Houziyan Dam, with the Yuke fault as the boundary) increased slightly, reflecting the influence of regional geological structures on the velocity structure. This characteristic still exists in the first stage of impoundment. In the first stage of impoundment, a large number of small induced earthquakes near the reservoir shore reflect the influence of reservoir water infiltration and its dynamic effects. Changes in gravity before and after impoundment also indicate a significant impact of reservoir water infiltration and load changes on the deformation field of the reservoir area near the shore. The water level of the reservoir has been fluctuating quasi-periodically between 6 and 40 m since 2018. From a physical perspective, this large-capacity, large water-level difference fluctuation will cause changes in the physical structure of the underground medium, such as pore pressure and permeability, making it more conducive to reservoir water infiltration. The low-velocity anomaly distribution in the second stage of impoundment is related to reservoir water erosion and infiltration, a phenomenon consistent with the results of P-wave velocity structure imaging studies conducted by Li Qiang et al. (2009) in the Three Gorges reservoir area (Li et al., 2009). Long-term large fluctuations in reservoir water level will change pore pressure, leading to deeper reservoir water infiltration. Therefore, in the third and fourth stages of impoundment, the area of high-velocity anomalies in the 1–6 km depth layer of the reservoir area and its surroundings increases, and the degree of anomaly change intensifies, reflecting an expanding influence of reservoir water on the velocity structure of this depth layer. This is consistent with the results obtained by Wu Haibo et al. (2018) in their study of the upper crustal velocity structure imaging in the Three Gorges reservoir area (Wu et al., 2018).

3.3 Imaging results of V_S

According to V_S in the section result maps (Figure 13) and profile maps (Figure 14), there is local lateral heterogeneity before impoundment of the Houziyan reservoir area. The manifestations include that a high wave velocity zone of 3300–3700 m/s appears from the Houziyan dam to the upstream reservoir area and in a 4 km depth area from the surface to underground; overall, the velocity increases as the depth becomes deeper, with a value of about 3300 m/s at the shallow area of 1 km to about 3800 m/s at the deep area of 10 km.

FIGURE 13. The horizontal section maps of V_S at different impoundment stages. (A–J) represent before impoundment in 3 km, before impoundment in 6 km, impoundment stage 1 in 3 km, impoundment stage 1 in 6 km, impoundment stage 2 in 3 km, impoundment stage 2 in 6 km, impoundment stage 3 in 3 km, impoundment stage 3 in 6 km, impoundment stage 4 in 3 km, and impoundment stage 4 in 6 km, respectively.

FIGURE 14. The vertical profile maps of V_S at different impoundment stages. (A–J) represent before impoundment along longitude, before impoundment along distance, impoundment stage 1 along longitude, impoundment stage 1 along distance, impoundment stage 2 along longitude, impoundment stage 2 along distance, impoundment stage 3 along longitude, impoundment stage 3 along distance, impoundment stage 4 along longitude, and impoundment stage 4 along distance, respectively.

In general, the characteristics of the velocity structure inversion results for P-waves and S-waves are very similar. Similarly, there are local lateral changes in V_S from stages 1–4 after impoundment. At stage 1, a high wave velocity zone of 3400–3900 m/s appears from the Houziyan dam to the upstream reservoir area and in a 5 km depth area from the surface to underground; at stage 2, a low wave velocity zone of 3000–3300 m/s appears from the Houziyan dam to the upstream reservoir area and in a 3 km depth area from the surface to underground; at stage 3, a high wave velocity zone of 3400–3900 m/s appears from the Houziyan dam to the upstream reservoir area and in a 5 km depth area from the surface to underground; at stage 4, a high wave velocity zone of 3400–3900 m/s appears from the Houziyan dam to the upstream reservoir area and in a 6 km depth area from the surface to underground.

Throughout all impoundment stages except stage 2, V_S changes from about 3300 m/s at the shallow area of 1 km to about 3800 m/s at the deep area of 10 km; compared with the case before impoundment, V_S has increased by 100–300 m/s. Compared with the cases before and at the first stage of impoundment, the original high-value zone of V_S at stage 2 becomes a low one with a variation of 400–600 m/s, which indicates that impoundment has resulted in changes in the V_S of the underground media.

The distribution characteristics of S-wave velocity from different impoundment stages show similar variations to those of P-wave velocity, but with some differences in detail. During the first impoundment stage, S-wave velocity shows a significant increase, with an expanded range of horizontal and depth variations. In the second impoundment stage, S-wave velocity also shows a low anomaly phenomenon, but the depth of the velocity anomaly is shallow, within 3 km. In the third and fourth impoundment stages, the changes in S-wave velocity are consistent with those of P-wave velocity, but the high-velocity anomaly of S-wave (2–5 km) is more pronounced than the deep anomaly of P-wave (1–4 km). The S-wave phenomenon indicates that the influence of reservoir water infiltration on S-wave velocity structure is minimal, but the gravitational effect of the reservoir water is significant, suggesting that water has a relatively small impact on the shear modulus (He et al., 2018).

3.4 Imaging results of V_P/V_S

According to V_P/V_S ratios in the section result maps (Figure 15) and profile result maps (Figure 16), there is local lateral heterogeneity before impoundment. A positive anomaly zone between 1.69 and 1.72 appears in the dam and its downstream as well as a 4 km depth area from the surface to underground; a negative anomaly zone between 1.67 and 1.69 appears in the east and west sides of the central reservoir area and a 4 km depth area from the surface to underground. The values of the remaining areas are between 1.68 and 1.70.

FIGURE 15. The horizontal section maps of V_P/V_S ratios at different impoundment stages. (A–J) represent before impoundment in 3 km, before impoundment in 6 km, impoundment stage 1 in 3 km, impoundment stage 1 in 6 km, impoundment stage 2 in 3 km, impoundment stage 2 in 6 km, impoundment stage 3 in 3 km, impoundment stage 3 in 6 km,impoundment stage 4 in 3 km, and impoundment stage 4 in 6 km, respectively.

FIGURE 16. The vertical profile maps of V_P/V_S ratios at different impoundment stages. (A–J) represent before impoundment along longitude, before impoundment along distance, impoundment stage 1 along longitude, impoundment stage 1 along distance, impoundment stage 2 along longitude, impoundment stage 2 along distance, impoundment stage 3 along longitude, impoundment stage 3 along distance, impoundment stage 4 along longitude, and impoundment stage 4 along distance, respectively.

At stage 1, a positive anomaly zone between 1.69 and 1.72 appears in the dam and its downstream as well as a 3 km depth area from the surface to underground; a negative anomaly zone between 1.66 and 1.69 appears in the central reservoir area and a 3 km depth area from the surface to underground. The values of the remaining areas are between 1.68 and 1.70. At the second impoundment stage, there are local lateral changes in V_P/V_S ratios. The manifestations include that a negative anomaly zone between 1.63 and 1.67 appears in the upstream of the Houziyan dam (northeast of Yuke fault as a dividing line) and a 6 km depth area from the surface to underground; a positive anomaly zone between 1.71 and 1.74 appears in the west of the Houziyan dam (southwest of Yuke fault as a dividing line) and a 2–8 km depth area from the surface to underground. At stage 3, there are local lateral changes in V_P/V_S ratios. The manifestations include that a negative anomaly zone between 1.66 and 1.68 appears in the tail area of the Houziyan reservoir and a 4 km depth area from the surface to underground; a negative anomaly zone between 1.66 and 1.68 appears in the east and west of the Houziyan dam and a 3–6 km depth area from the surface to underground; a positive anomaly zone between 1.70 and 1.73 appears from the central part of the reservoir to the downstream of the dam and a 4 km depth area from the surface to underground; the values of the remaining areas are between 1.68 and 1.70. At stage 4, there are local lateral changes in V_P/V_S ratios. The manifestations include that a negative anomaly zone between 1.66 and 1.68 appears in the tail area of the Houziyan reservoir and a 4 km depth area from the surface to underground; a positive anomaly zone between 1.70 and 1.73 appears from the central part of the reservoir to the downstream of the dam and in a 4 km depth area from the surface to underground; and the values of the remaining areas are between 1.68 and 1.70.

Throughout all impoundment stages except stage 2, V_P/V_S ratios change from about 1.70 at the shallow area of 1 km to about 1.75 at the deep area of 10 km; compared with the case before impoundment, this value does not change significantly, but the range of abnormal horizontal position changes is large, which means the variations of V_P and V_S are different after impoundment. At stage 2, in a 4 km-deep area from the surface to underground, the values of V_P/V_S ratios at each side of the Yuke fault are one high and one low, which clearly shows the position and depth of the fault.

The V_P/V_S velocity structure imaging results show different anomaly features compared to P-wave and S-wave. In the upstream reservoir area, the V_P/V_S velocity anomaly changes are consistent with P-wave and S-wave, but the anomalies are mainly concentrated at shallow depths, indicating that reservoir impoundment has different effects on V_P, V_S, and V_P/V_S. On both sides of the Yuke fault, the V_P/V_S velocity anomalies show different variations at different impoundment stages. In the southwest of the Yuke fault, the V_P/V_S velocity anomalies show consistent positive anomalies, while in the northeast of the Yuke fault (upstream of the dam), V_P/V_S shows both positive and negative anomalies, indicating that the reservoir impoundment does not penetrate the fault, reflecting the inhibitory effect of the fault on the underground velocity structure changes. At the same time, the different impoundment stages exhibit different V_P, V_S, and V_P/V_S velocity structure anomaly characteristics, reflecting the dynamic adjustment of the reservoir impoundment’s influence on the underground velocity structure with changes in the reservoir water level.

4 Discussion and conclusion

Through tomography, we have gained insights into the subterranean velocity structure and its variations in the Houziyan reservoir area before and after impoundment, as well as during different impoundment stages. The following conclusions can be drawn through a comprehensive analysis and comparison of the results:

1) Anisotropy in the underground velocity structure is evident in the Houziyan reservoir area before impoundment. Specifically, from the dam to the upstream reservoir area and within a depth range of 4 km from the surface to underground, there is a distinct high wave velocity zone of V_P、V_S and V_P/V_S, reflecting the anisotropic nature of the underground velocity structure in the pre-impoundment reservoir area.

2) The area of high wave velocity increases in stage 1, stage 3 and stage 4. The area of low wave velocity increases in stage 2, especially in depth, indicating that impoundment has altered the underground velocity structure. These changes vary significantly across different impoundment stages. Notably, high wave velocity zones of both V_P and V_S appear from the dam to the upstream reservoir area and extend to a depth of 5 km from the surface to underground in stage 1. Similar high wave velocity zones appear in stages 3 and 4, with depths extending to 6 km from the surface to underground. In contrast, low wave velocity zones of V_P and Vs are observed in stage 2, extending from the dam to the upstream reservoir area and to a depth of 5 km from the surface to underground. The influence depth of these changes is greater than that observed before impoundment.

3) The rule of the underground velocity structure changes in the Houziyan reservoir differs from the patterns observed in other reservoirs after impoundment. For instance, at stage 2, low wave velocity zones for V_P and Vs are observed from the dam to the upstream reservoir area, reaching a depth of 5 km from the surface to underground. This pattern contradicts the results observed before and during the first stage of impoundment. At stage 2, the values of V_P/V_S ratios are bounded by the Yuke fault, decreasing in the northeast area of the fault, with a change depth of approximately 4 km from the surface to underground. Conversely, V_P/V_S ratios increase in the southwest area of the fault, with a change depth ranging from 2 to 6 km. This indicates that the changes in P-wave and S-wave velocities for underground media differ with increased impoundment time. The depth of this fault can be inferred from the value changes on both sides.

4) In general, the underground velocity structure exhibits an overall increasing trend after impoundment, but there are instances of velocity decrease, reflecting the continuous dynamic adjustment of the underground velocity structure following impoundment.

The study of underground velocity structure before and after reservoir impoundment is helpful for us to gain a deeper understanding of the process and extent of reservoir water infiltration. It should be combined with the hydrological conditions and geological structures of the reservoir for detailed analysis. Regarding the anomalies in the second impoundment stage compared to other stages, it is speculated that the relatively less permeation of reservoir water may be due to the rapid rise and fall of the reservoir water level in the early stage of impoundment, and the underground velocity structure is mainly influenced by the gravitational changes of the reservoir water in dynamic adjustment. Currently, due to the lack of detailed hydrological and geological data, more seismic geological data will be used to further improve the research results in the later stage of the study. Overall, there exists a correlation between the underground velocity and reservoir water storage in the Houziyan reservoir area. The changes in the underground velocity structure vary across different impoundment stages, underscoring the influence of reservoir impoundment on the subterranean velocity structure.

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 authors.

Author contributions

ZY: Investigation, Writing–original draft. JC: Conceptualization, Writing–original draft. YW: Data curation, Writing–review and editing. JQ: Methodology, Writing–review and editing. HH: Formal Analysis, Writing–review and editing. XD: Validation, Writing–review and editing. MY: Resources, Writing–review and editing. JW: Resources, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was sponsored by the Key Scientific Research Projects of the Power China Guiyang Survey Design and Research Institute Co., Ltd. (YJ2022-09; DJ-ZDXM-2020-55); Chongqing Technology Innovation and Application Development Project of Chongqing Science and Technology Bureau (CSTB2022TIAD-KPX0197); Guizhou Provincial Science and Technology Projects [QKHJC-ZK(2022)-General 186]. Projects of National Natural Science Foundation of China (No. 41602275 and 41977298); Provincial Key Disciplines of Guizhou Province-Geological Resources and Geological Engineering [ZDXK(2018)001].

Conflict of interest

Authors ZY, YW, and XD were employed by Power China Guiyang Survey Design and Research Institute Co., Ltd.

Author HH was employed by Guodian Dadu River Basin Hydropower Development Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References