1. Introduction
Hydropower has gradually become one of the important ways to solve the national energy crisis, and many countries have built large reservoirs. In the 1940s, Carder (1945) first found the correlation between impoundment and earthquakes in Lake Mead, USA. The water led to the reactivation of pre-existing faults and further induced earthquakes in the reservoir region. With the construction of reservoirs, researchers found that impoundment and periodical variations in water levels are strongly correlated with seismicity in the reservoir area. Seismic activity in the reservoir and surrounding areas became abnormal after impoundment, and in some cases, it even increased [
1,
2,
3,
4] (Carder, 1945; Simpson et al., 1988; Talwani, 1997; Cao et al., 2015). The earliest seismicity associated with reservoir impoundment occurred near the Kremasta Reservoir in north Greece [
5]. Reservoir earthquakes have occurred in more than 140 regions around the world, especially in the United States, China, and India [
4,
6].
In China, more than 90,000 reservoirs have been built, with a total capacity of nearly 900 billion cubic meters. Researchers have proven that about 20 earthquakes had been related to reservoir impoundment by 2012 [
7]. One example, the 1962
6.1 Xinfengjiang earthquake in Guangdong Province, is the earliest and most destructive reservoir-induced earthquake in Chinese history [
8]. Large reservoirs have been concentrated in southwest China since the 20th century, where strong earthquakes have occurred frequently. Therefore, the correlation between the construction of large-scale hydropower projects and seismicity changes in the surrounding area has drawn attention in academia for years, especially the relation between the occurrence of the 2008
Wenchuan earthquake and the impoundment of the nearby Zipingpu Reservoir [
9,
10,
11,
12]. The Jinsha River is located in the important south-north seismic zone in Southwest China (
Figure 1a,b). It is adjacent to the Sichuan-Yunnan Block in the west and the South China Block in the east, which has always been a region with high seismicity. In recent years, hydropower stations and reservoirs have been built along the Jinsha River to harness its significant elevation drop and abundant water resources. Four reservoirs, Wudongde, Baihetan, Xiluodu, and Xiangjiaba, were built along the river from the upper to the lower, and several notable earthquakes have occurred. For instance, Xiluodu and Xiangjiaba stations are located in the Mabian-Zhaotong seismic zone, where there have been
earthquakes [
13], and Wudongde and Baihetan Stations are located in the Dongchuan-Songming seismic zone, where there have been
earthquakes [
14]. Large earthquakes often result in significant losses to human society and the economy, posing a threat to people’s lives. Therefore, it is of great significance to study the possibility of earthquake occurrences and evaluate earthquake losses reasonably [
15,
16,
17]. Consequently, there is a research priority on whether the construction of hydropower stations will induce large earthquakes in active fracture zones.
Baihetan is located in the southern Daliangshan sub-block, near the junction of the Zemuhe and Lianfeng faults [
22]. Most of the nearby faults are high-angle strike-slip. There is also the thrust Longmenshan fault in the east-north direction (
Table 1). Zhang et al. (2021) [
23] suggested that the seismicity level around the reservoir was quite low before the impoundment of the Baihetan Reservoir but it significantly increased after the impoundment, and
earthquakes occurred at the same time (
Figure 1c). Therefore, it is very important to carry out further research on the relationship between reservoir impoundment and regional seismicity in the Baihetan Reservoir.
In general, seismicity results from the change of reservoir water level or impoundment, which is called a reservoir-induced earthquake. These earthquakes are different from normal tectonic earthquakes, and their occurrence times are strongly correlated with reservoir impoundment. Usually, microquakes occurred immediately or several months after the reservoir began to store water, while larger shocks occurred one or more years later. Different impoundment processes for reservoirs with different geological conditions affect the occurrence time and locations of larger shocks; Dou et al. (2022) [
24] studied and predicted the potential seismicity changes caused by the Baihetan Reservoir before its impoundment, using the poroelastic equation and an assumed water storage model. However, after the impoundment of the Baihetan Reservoir, there have been few studies on seismicity using its actual impoundment model. Furthermore, the long-term effect of reservoirs on seismicity has not been studied. There are still many problems to be discussed and studied, such as what will be caused by the continuous downward infiltration of reservoir fluid and whether long-term loading will affect seismicity.
This study proposed a new pore-viscoelasticity coupled mechanism, derived physical equations, and developed a finite element program to study the influence of the long-term effect of a reservoir on seismicity. We established a three-dimensional finite element model for the Baihetan area, calculated Coulomb failure stress on the major faults, analyzed seismicity changes using water level data as the boundary condition, and discussed the long-term influence of the Baihetan Reservoir on seismicity. This study may provide a theoretical basis and valuable reference for future reservoir-earthquake hazard mitigation.
2. Pore-Viscoelastic Method
Considering the long-term use of the reservoir after construction and the viscosity of the crustal medium, it is necessary to construct fully coupled pore-viscoelasticity equations to study the impact of reservoir impoundment on strata. According to the instantaneous elastic effect caused by water storage and the corresponding viscous relaxation effect for a long time, we derived a pore-viscoelasticity constitutive equation describing the evolution of underground media with time from fluid-saturated elastic porous media equations [
25]. Therefore, we can better analyze the early elastic response and long-term relaxation effect of the crust. In Rice and Cleary’s theory, the diffusion and the elastic constitutive equations are as follows:
where
denotes the stress tensor,
denotes the strain tensor, G is the shear modulus, and
is the diffusivity related to Poisson’s ratio
, undrained Poisson’s ratio
, Skempton coefficient
, and permeability
. Equation (2) reflects the diffusion of the pore pressure with time.
The term
in Equation (1) is the elastic stress caused by strain. We replaced it with a viscoelastic integral form constitutive equation based on the Boltzmann superposition principle (4), so we get the pore-viscoelastic constitutive equation with an integral form (5):
where
is a Kronecker delta function and K is the constant of bulk modulus, and
for maxwell viscoelastic medium, and
when
, in which
and
are the constants of shear modulus and viscosity, respectively.
Next, we transform Equation (5) to simplify it:
where
where
Then we make some changes to
, as
, and use the theorem
Thus, we can transform Equation (7) to
For quasi-static phenomena, we use
and
to express the governing equations, where the stress
must satisfy the equilibrium equation.
in weak form,
where
denotes the internal force such as gravity, and the
is the face force on the surface.
Using a forward difference scheme in the time domain
, and according to Equations (2) and (11)–(13), we can get the finite element weak solution form of the equations:
The superscript represents the result of the previous calculation step. Based on the finite element weak solution from Equations (14) and (15), we write the calculation program to simulate the problem of reservoir-induced earthquakes.
3. Model
We utilized the self-developed program to generate the three-dimensional pore-viscoelasticity finite element numerical model of the Baihetan Reservoir area (
Figure 2a), which mainly contains the Zemuhe Fault (ZMHF), Daliangshan Fault (DLSF), Lianfeng Fault (LFF), Ninghui Fault (NHF), Puduhe Fault (PDHF), north part of the Xiaojiang Fault (XJF1), Xiaojiang Fault (XJF2), and Zhaotong-Ludian Fault (ZLF). According to Hui et al. (2020) [
26], fault geometry has a certain influence on stress distribution in porous media. The fault geometry in the study area is incomplete. The fault geometry used in calculations is presented in
Table 1. We assigned a vertical strike-slip mechanism to the unspecified faults due to the high fault dip in this area. In this example, we utilized four nodes to generate grid files in parallel (
Figure 2), which can effectively improve the speed of model generation. Firstly, the study area was divided into multiple subregions (
Figure 2b takes four subregions as examples), and the vertical grid of each subregion was generated based on the thickness of each stratum and the predefined grid density. At the same time, the grid nodes were adjusted using real terrain data. Finally, the subregional grids were summarized to form an overall grid. The parallel mesh generation program can effectively enhance the calculation speed for complex terrains and boundary conditions, significantly improving the efficiency of mesh generation.
According to Liu et al. (2021) [
27], we divided the 30 km depth model into four layers, and the material parameters are shown in
Table 2. SRTM30 interpolates the surface, and the boundary conditions of the fixed normal and free tangential direction are set on the model’s surroundings and bottom. The variations of the strain and stress are calculated from the beginning of water injection for 200 years from 2021 to 2221, and the initial response to the impoundment is mainly calculated from April 26 to 23 September 2021, a total of 150 days after the initial water storage. Since the impoundment began, the water level of the reservoir has been steadily rising and has now stabilized at a high level. Thus, the part of the surface DEM (Digital Elevation Model) below the water level will be subjected to a pressure of size
and a face force of size.
, where
is the water density,
is gravity,
is the difference between DEM and water level, and
represents the unit area.
4. Results
The computed model predicts deformation, stress, and pore pressure for each time step. The results of surface deformation are shown in
Figure 3. The deformation is affected by both water gravity load and fluid pore pressure, and the continuous rise of the reservoir water level gradually increases the vertical compression, directly impacting the strata under the reservoir. In addition, the gradually increasing water load and the continuous downward diffusion of pore pressure cause the compressive stress under the reservoir to increase (
Figure 4). During the impoundment period of the reservoir, the water gradually squeezes towards both sides of the Jinsha River. This causes the principal stress in the reservoir center to be characterized by NW-SE compression and NE-SW tension, with near NS compression in the periphery. As the depth increases, the stress variation decreases. The stress begins to decrease due to the relaxation effect after the impoundment is completed. The stress changes share the same characteristics as the displacement changes.
Reservoir impoundment has different effects on the seismicity of different faults. Faults are affected by the pressure of the reservoir water gravity in the crust, and the stress imparted by the water gravity depends on the fault geometries. The closer to the reservoir center, the more significantly the compressive stress of the fault increases. According to the Coulomb failure stress change theory [
28,
29,
30], an increase in compressive stress will reduce Coulomb failure stress (CFS) and the seismic risk of the fault, while an increase in pore pressure will enhance the CFS and seismic hazard. Based on the fault geometries of
Table 1, we selected the reasonable values (
Table 3), and combined with the stress state (
Figure 4), we calculated the CFS changes (
) of all the faults at different times via Equation (16),
where
and
are the normal stress change and the shear stress change in the rupture direction, which are calculated based on fault geometries, respectively, and
is the pore pressure change. μ is the friction coefficient, and a small value represents the easily slipping fault. This study set used μ = 0.5 in the calculation because all fault seismicity in the study area is active.
The XJF1 and XJF2 are the two closest faults to the reservoir center, so pore pressure grows fast. PDHF and ZLF are far from the reservoir; thus, the pore pressure increases relatively slowly (
Figure 5a). The
on XJF1 and ZLF increases with time, indicating increasing seismic risks on the faults. However, the combined effects lead to the decrease of
on XJF2 and PDHF, which suggests a lower possibility of earthquake occurrence (
Figure 5b). The pore pressure and
on each fault become stable when the reservoir impoundment is completed, which indicates that the influence of reservoir impoundment on the underground pore pressure is stabilized. Regional stress also gradually reaches a new equilibrium state, which may lead to the stabilization of seismicity.
Since other faults were far away from the reservoir and the seismicity changes were not obvious after impoundment, the relationship between the
and seismicity on XJF1 and XJF2 was primarily studied and discussed. The
on the fault plane of XJF1 and XJF2 along the Jinsha River and the earthquake distribution within 10 km were drawn to explore the relation between the
and seismicity (
Figure 6). Compared with the background seismicity, the earthquake almost occurred near the area or below where the
was positive after the impoundment. One hundred fifty-two earthquakes occurred near the XJF1 after the impoundment, all of which occurred close to the area where the CFS increases, and there were 62 earthquake epicenters with
exceeding 10 kPa, compared to the data of 11 and 5, respectively, before impoundment. The XJF2 seismicity was more complex. A total of 112 earthquakes occurred near the XJF2 after the impoundment; 39 earthquakes were located near the negative area of
and 73 earthquakes near the positive area, where the
at the epicenter of 27 earthquakes exceeded 10 kPa, with a ratio of about 1.87. However, there were 20 earthquakes on the XJF2 before impoundment, and the ratio of the positive and negative areas of
was
. Therefore,
does have an impact on the seismicity. In particular, three
earthquakes occurred near the fault plane with larger
after the impoundment, which reflects the impoundment having a certain role in promoting the frequency and magnitude of earthquakes. The
was negative on the southern side of the XJF2, which corresponded to a decrease in seismicity in the area. According to the collected seismic records, there was less seismicity in this area, and a seismic gap was formed at the southern end, which further validated our calculation results.
5. Discussion
According to the seismic records, earthquakes were widely distributed before the reservoir impoundment, with many on each fault (
Figure 1). We have drawn the water level change and earthquakes’ M-T (
Figure 7a,b), and compared the focal depths of the background earthquakes before the impoundment with those of earthquakes during the first year of impoundment (the seismicity before impoundment is called “background earthquakes” in this paper) (
Figure 7d). The seismic b-value is one of the important indexes for seismological research and seismic risk assessment.
We get the Gutenberg–Richter (G-R) magnitude–frequency relationship in our study area and use Equation (17) to calculate the seismic b-value (Utsu, 1966) [
31] (
Figure 7c). Our seismic data show great integrity, and the earthquake frequency increases significantly after the reservoir begins to impound with b-value increases from 0.8179 to 0.8430. It may indicate that the accumulated stress is released due to the infiltration of pore fluid, and mostly small earthquakes with magnitudes smaller than
M3.0 (226 earthquakes smaller than
.0), and
earthquakes rarely occur (only 5 times) (
Figure 8). For the depth distribution, the earthquake numbers at depths shallower than 10 km have not increased, but the numbers at depths deeper than 10 km have increased significantly, from 84 to 197. Especially at depths between 10 km and 15 km, there were 61 earthquakes before, but the number of earthquakes increased to 139 after impoundment, even more than twice as much as before (
Figure 7d). This reflects that the previously unaffected deep faults are gradually opened by pore pressure as the reservoir water level gradually increases, resulting in shear stress on the fault plane exceeding the sum of its friction and cohesion, and earthquakes occur on the faults to release their energy. Because the fault stress is released before it is fully accumulated, there are usually more earthquakes than background earthquakes, with generally smaller magnitudes. There were 258 smaller and 5 bigger-than-
earthquakes before, and 390 smaller and 6 bigger-than-
M3.0 earthquakes after impoundment.
where
is the average magnitude,
is the minimum complete magnitude, and
is the minimum magnitude interval.
The
earthquake, which occurred between LFF and ZLF on 18 May 2022, brought a large number of aftershock activities (
Figure 8). It may have caused many earthquakes to occur in the NE direction away from the reservoir, and there seems to be no distinct change in seismic records between LFF and ZLF. Compared with background seismic records, the later earthquakes at all depths are concentrated in the area with the largest pore pressure increase below the reservoir center, and even a
earthquake occurs at a depth of 5 km (
Figure 8a). The seismicity with the focal depths of 5~10 km below the reservoir center increased, from 164 to 171, and the
earthquakes gradually migrated to the reservoir center after impoundment (
Figure 8b).
The number of earthquakes with focal depths between 10 and 15 km has significantly increased after impoundment, mostly distributed in the areas near the reservoir center, where the pore pressure increases the most (about 40 kPa). The seismic rates in the other areas far from the reservoir are constant or even decreased compared with the pre-impoundment seismicity level (
Figure 8c). This may be because the shallow pore pressure has increased and the accumulated stress has been released due to river seepage before the impoundment. In contrast, deep pore pressure also increases due to growth of the reservoir water level, and the fault has a stimulation effect on the earthquake.
The pore pressure below the reservoir center exceeds the typical static triggering threshold (10 kPa) after impoundment. These earthquakes can be considered induced earthquakes. Nevertheless, the impoundment does not induce earthquakes immediately in the areas far from the reservoir center. For example, the Qiaojia
earthquake on 30 March 2022, following another Qiaojia
earthquake on 18 May 2020 (
Figure 8b) might be considered a triggered earthquake resulting from the Baihetan Reservoir impoundment. However, previous studies have shown that a large earthquake can be triggered when the pore pressure increases by 0.1 bar [
32]. We find that the pore pressure at the epicenter of this
earthquake increases by ~1180 Pa according to the calculation results, and the stress caused by reservoir gravity is also on the order of kilopascals, which is much smaller than the typical static triggering threshold (10 kPa); thus, the reservoir may promote the occurrence of earthquakes in advance but the impact is limited. Therefore, the
Qiaojia earthquake can be considered a natural event.
According to a previous study [
33], the viscosity of the upper mantle is ~
, and the viscosity of the crust is less than
in the Sichuan-Yunnan region. Reservoirs are usually designed to serve as power supplies for hundreds of years. To study whether the long-term existence of the reservoir has an impact on regional seismicity, we have calculated the 200 years after the completion of the reservoir impoundment and calculated the fault stress changes resulting from the viscosity characteristics of the crust and mantle (
Figure 9). The results indicate that the long-term existence of the reservoir gravity load results in the relaxation effect in the surroundings, and the stress on the fault plane, will change constantly. However, the normal stress of the fault has changed slightly within 200 years. Even for the XJF1 below the reservoir center, the normal stress change caused by viscosity is no more than 20 Pa, which is far less than the stress changes imparted by tectonic loading. Consequently, the reservoir seismicity will be stabilized after the equilibrium of pore pressure and fluid field, and the long-term existence of the reservoir has little impact on seismicity.