Upper crustal structure of the Xinfengjiang reservoir from ambient noise double beamforming tomography and its implications for induced seismicity

SUMMARY

The Xinfengjiang Reservoir in Guangdong, China, triggered one of the largest reservoir-induced earthquakes in the world, with a magnitude of 6.1 in 1961 March. Frequent local seismic events have been recorded in the vicinity of the reservoir, posing a potential threat to the dam infrastructure and surrounding regions. In this study, we obtain a high-resolution S-wave velocity model of the upper crustal structure underlying the Xinfengjiang Reservoir, using double beamforming tomography method from newly deployed dense array across the reservoir. Our results reveal that several high-velocity structures are distributed beneath the Xinfengjiang Reservoir. These structures are robust, facilitating a greater accumulation of stress before fault slip. Earthquakes beneath the reservoir are primarily concentrated at the boundaries of the high-velocity bodies, indicating their controlling role on the location of seismic events. Low-velocity zones, acting as pathways for fluid migration, extend from the surface to the depths of seismic sources, thereby inducing earthquakes by elevating pore pressure within fault and fracture systems. Due to significant differences in the fluid diffusion coefficient, the delay times between the onset of earthquakes and the peak water levels vary considerably with the different crustal structures.

1 INTRODUCTION

Human activities, such as mining (McGarr 2002), wastewater disposal (Goebel et al. 2016; Improta et al. 2017; He et al. 2021), unconventional oil and gas production related to hydraulic fracturing (Schultz et al. 2020), geothermal extraction (Ellsworth et al. 2019) and surface water impoundment in reservoirs (Gupta 2002; He et al. 2018; Dong et al. 2022) have the potential to induce earthquakes by changing the stress or strain state of the subsurface. The growing demand for energy has driven the rapid expansion of water conservation and hydropower projects, leading to a rise in the number of reservoirs worldwide (Ding 1989). Typically, induced seismicity resulting from water impoundment in reservoirs tends to occur at shallow depths in the vicinity of the reservoirs, posing potential threats to both the reservoir dams and surrounding regions (e.g. Chen & Talwani 1998; Gupta 2022). For instance, the 1961 M 6.1 earthquake impacted the Xinfengjiang Reservoir in China, resulting in a 108-m-long crack in the reservoir dam embankment (Song et al. 2017). Similarly, the 1967 M 6.3 earthquake in the Koyna Reservoir resulted in the deaths of approximately 200 individuals and caused severe damage to nearby towns (Gupta 2022).

The Xinfengjiang Reservoir has been recognized as the site that triggered the world's largest reservoir-induced earthquake (M 6.1 in 1961 March). It is situated on a large granitic mass trending east–west and is bounded to the east by the Heyuan Basin (Tannock et al. 2019). The granitic intrusions beneath the Xinfengjiang Reservoir formed during the Jurassic to early Cretaceous crustal extension in South China, which undergone several episodes of deformation during the Cenozoic (Shu 2012; Li et al. 2017). Mountainous and hilly rocks of Devonian, Jurassic, and Cretaceous ages are widely distributed around the southern and northern boundaries of the reservoir (Ye et al. 2017). Well-developed NE-trending and NW-trending groups of fault systems have formed in the reservoir since the Meso-Cenozoic era (Fig. 1). Among these faults, the NE-trending faults system is the most prominent and well-developed, with associated surface expressions including the Heyuan fault (HYF), Renzishi fault (RZSF) and Daping–Yanqian fault (DYF) (Ding et al. 1983; Yang et al. 2013). These fault systems intersect, characterizing the complex geotectonic settings and profoundly influencing the occurrence of regional induced seismicity in the study area (He et al. 2018; Dong et al. 2022).

The Xinfengjiang area experienced relatively infrequent seismic activity until the impoundment of the reservoir began in 1959 (e.g. Shen et al. 1974; He et al. 2018; Gupta 2022). Following the commissioning of the reservoir, a significant increase in seismic activity was recorded within the reservoir (Ding et al. 1982; Guo et al. 2004). Seismic activity has persisted to the present, with over 50 earthquakes recorded that exceeded a magnitude of 4.0 (He et al. 2018). Seismicity in this area is primarily concentrated around the dam, but since 2012, earthquakes have become more prominent in the northwest region of the reservoir (Wang et al. 2018; Dong et al. 2022), with the latest being an ML 4.5 earthquake occurring near Wudong in 2023 March. Numerical simulations indicate that the pore pressure diffusion plays a vital role in triggering earthquakes, with pore pressure increases exceeding 0.01 MPa at a depth of ∼8 km in the Xinfengjiang Reservoir (Cheng et al. 2012). In addition, the five-year periodic changes in water levels (fluctuating between 90 and 115 m) have also proven to have a significant impact on seismicity in the Xinfengjiang Reservoir (He et al. 2018). Therefore, various factors, including the state of the regional stress field, fracture permeability, water level fluctuations, etc. control the occurrence of reservoir-induced seismicity (Zoback & Hickman 1982; Talwani 1997).

To elucidate the mechanisms responsible for these induced earthquakes, a more detailed understanding of the structures beneath the Xinfengjiang Reservoir is essential. Ye et al. (2017) reported that the upper crust beneath the Xinfengjiang Reservoir exhibits an alternating distribution of high- and low-velocity anomalies at depths of 5–10 km. Moreover, seismic events were primarily recorded within the high-velocity bodies and the transition zone occurring between high- and low-velocity regions. Further insights were provided by He et al. (2018), who identified a low-velocity zone along the northwest–southeast direction at depths of 7–10 km beneath the HYF. Dong et al. (2022) proposed a model comprising channels for fluid infiltration in a 3-D S-wave velocity structure of the upper crust beneath the reservoir. According to their model, water infiltrates through the shallow faults (such as DYF and RZSF) and penetrates to greater depths.

Previous studies were limited by the relatively low resolution of the velocity models in Xichang and Heyuan Basin, primarily due to the sparse coverage of seismic recording stations in the area (e.g. Ye et al. 2017; He et al. 2018; Dong et al. 2022). Consequently, the relationship between the increased seismicity and subsurface structures in the Xichang and Wudong areas has not been clearly explained. At depth, seismic events are mainly distributed along the middle section of the Heyuan fault, located near the dam, while no seismic activity has been recorded at the tail end of the fault (Wang et al. 2018). The mechanisms underlying these variations in seismic activity remain ambiguous.

In this study, we conduct a high-resolution investigation of the structures beneath the Xinfengjiang Reservoir, utilizing ambient noise recorded from a newly installed temporary seismic network consisting of 603 stations. Our study takes advantage of the new ambient noise Rayleigh-wave phase velocity data to produce a high-resolution model of the fine structures beneath Xinfengjiang Reservoir. Our new models and results elucidate the relationship between seismicity and crustal structure, providing insights into the mechanisms underlying variations in seismic activity across the depth profile of the Heyuan fault.

2 DATA AND METHODS

To obtain the continuous seismic data used in this study, a total of 603 short-period three-component seismometers (Smartsolo IGU-16HR 3C) were deployed for about 50 d, from 2021 June to August. The dense seismic array features an average station interval of 2 km and commendably covers the Xinfengjiang Reservoir region (Fig. 1a). The elevations of the stations are relatively low (∼200 m), so the effect of topography is negligible in the subsequent analysis and inversion (Köhler et al. 2012; Ping et al. 2018; Wang & Sun 2018). The vertical components’ continuous time-series data are used to compute the ambient noise cross-correlation and subsequently extract Rayleigh-wave phase velocity.

2.1 Cross-correlations

We follow the procedure of Bensen et al. (2007) to calculate the cross-correlations between seismic station pairs. We first cut the continuous vertical component waveform data into smaller 1-hr segments, decimate the sampling rate to 25 Hz, remove mean and trend, and applied a band-pass filter to the data at periods of 0.2–10 s. Time-domain normalization and spectral whitening are applied to the pre-processed data, and the hourly cross-correlations for station pairs are then calculated. Thereafter, we linearly stack all cross-correlations for each station pair to obtain 181 503 interstation empirical Green's functions (Bensen et al. 2007). Fig. 2 shows examples of the vertical–vertical (ZZ) component of cross-correlations obtained as earlier described. Fundamental Rayleigh wave signals can be clearly seen on both negative and positive time coordinates of the cross-correlations. However, the amplitudes of positive and negative lags are significantly different, that is not time-symmetric for stations at the edge of the array (e.g. station 10 471 and 06 256 in Fig. 2), indicating the heterogeneous distribution of ambient noise sources in azimuth. Previous studies have suggested that the anisotropic distribution of noise energy has little effect on the extraction of the phase velocity (Yang & Ritzwoller 2008). Therefore, to fully utilize the causal and acausal parts of the cross-correlations, we compress the two-sided signals into one-sided signals by averaging the positive and negative sides of the cross-correlations. Finally, to mitigate the overlap between neighbouring sampling, narrow-band filtering (±10 per cent), utilizing a Butterworth filter, is applied for periods of 1.0–5.0 s. The final ‘symmetric’ signals with signal to noise ratios greater than 5 are retained for further analysis.

2.2 Double beamforming tomography

In this study, we apply the double beamforming (DBF) method (Krüger et al. 1993; Wang et al. 2019a, b; Wu et al. 2022) to directly measure the Rayleigh wave phase velocity of the array. In recent years, the DBF has been widely used to enhance weak body and surface wave signals from cross-correlations (Boué et al. 2013; Nakata et al. 2016; Castellanos et al. 2020), reconstruct ray paths (Boué et al. 2014), and measure azimuthal anisotropy (Wu et al. 2022). The array-based DBF method provides a direct approach to measure phase velocity compared to the traditional ambient noise tomography method. This eliminates errors introduced by assumptions and intermediate steps that are inherent in conventional phase velocity inversion. Moreover, the uncertainty of local phase velocity can be estimated effectively and conveniently through repeated measurements among different source–receiver pairs (Wang et al. 2019a, b; Wu et al. 2022, 2023).

Given a source beam centred at $X_s^c$ and a receiver beam centred at $X_r^c$⁠, the stations within D/2 from the centre of source/receiver beam are selected as the source/receiver stations, where ‘D’ is defined as the beam width (Fig. 3a). Cross-correlations between source stations and receiver stations are used for subsequent DBF analysis. To avoid inadequate stacking of cross-correlations for each beam pair, ensure the stability of the computation, and optimally improve the transverse resolution of the results (Fig. S1, Supporting Information), we set D at 6 km (Wang et al. 2019a; Wu et al. 2023). A far-field criterion (Yao et al. 2006; Luo et al. 2015) is adopted to remove beam pairs with distances less than 1.5 wavelengths (Wang et al. 2019a; Wu et al. 2022, 2023). To isolate the fundamental Rayleigh wave signals, the cross-correlations are windowed between the minimum time (⁠$d/V_{\max }-T_0$⁠) and the maximum time (⁠$d/V_{\min }+T_0$⁠). Here, d represents the interstation distance, and $V_{\min }$ and $V_{\max }$ are determined by the apparent propagation velocity of the cross-correlations for each period $T_0$⁠. The truncated waveforms are then normalized by their maximum amplitudes (Figs 3b and c) and stacked in the frequency domain (Boué et al. 2013, 2014; Wang et al. 2019b; Wu et al. 2022). The shifting and stacking processes conducted in the frequency domain can minimize errors caused by the sampling rate in time domain, and the stacked waveforms are then converted into time domain.

Assuming great-circle propagation, we employ a 2-D grid search method to compute the optimal local slowness $u_{\mathrm{s}}$ and $u_{\mathrm{r}}$ in each source beam and receiver beam, so that the energy of the envelope of the stacked waveform is maximum. In addition, a coarse-to-fine grid-search strategy is adopted to save computation time meanwhile ensure the computation accuracy (Fig. 4). The slowness for each beam centre is obtained through the iteration of different source–receiver beam pairs, and the final local slowness is determined by averaging all slowness within two standard deviations (Wu et al. 2023).

2.3 Shear-wave velocity inversion

In order to obtain the shear-wave velocity structure of the study area, an iterative least-square 1-D inversion algorithm (Herrmann 2013) is employed. The 1-D initial model for the entire region is derived from the average estimation from He et al. (2018), and the final Vs velocity model is computed and obtained after 30 iterations. As the Vs is updated in each iteration, Vp is calculated using a fixed ratio of Vp/Vs = 1.7, and the density is determined based on an empirical relationship between density and Vp (Brocher 2005). Fig. S2 (Supporting Information) shows an example of the 1-D shear-wave velocity inversion and fitted dispersion curve at a depth range of 0–11.0 km. Based on our observations, the data fitting is generally satisfactory. The resolvable depth range of the model is determined by analysing the sensitivity kernels for fundamental Rayleigh wave phase velocity at different periods (Fig. S3, Supporting Information). The 3-D shear-wave velocity model of the whole region can then be obtained by interpolating the 1-D shear-wave velocity results inverted from all the beam centres.

3 RESULTS

3.1 Checkerboard tests

Following the procedure of Wu et al. (2023), we conduct checkerboard tests to estimate the lateral resolution of the phase velocity model obtained from DBF. Each cell is set to a size of 6 × 6 km. The test model has a background velocity of 3.0 km s^-1 and a peak perturbation of ±6 per cent. Unlike the traditional tomography approach based on ray-tracing strategy, the DBF method directly extracts the local phase velocity by analysing the waveforms of each source–receiver pair. Therefore, we use Ricker wavelets to simulate the synthetic cross-correlation functions for different frequencies. Subsequently, we apply time-shift to the synthetic cross-correlations of all station pairs with respect to their traveltimes. Traveltimes are calculated using the fast-marching method (Sethian 1996). Gaussian random noise, with a mean of 0 s and a standard deviation of 0.01 s, is added to the synthetic data to simulate the noise level present in the observed data (Rawlinson et al. 2014). To better estimate the lateral resolution of each frequency, the cross-correlations analysed in each frequency are identical to those in the actual calculations. As shown in Fig. 5, there is a high level of velocity model recovery as revealed across each period. In addition, the tests show high resolution across the majority of the study area, including Xichang, the Heyuan Basin, and the centre of the reservoir.

3.2 Phase velocities and uncertainties

For each period, the slowness and corresponding uncertainties of each beam are used to construct the phase velocity maps (Fig. 6). Generally, phase velocity in our study area increases with period. At short periods (<3 s), a low-velocity anomaly is predominantly observed across the Heyuan Basin, and its morphology aligns well with the geological structure (Figs 6a–d). At longer periods (>3 s), high-velocity anomalies are observed beneath most parts of the Xinfengjiang Reservoir (Figs 6e and f). The uncertainty of the estimated phase velocity is generally low, indicating that the DBF analysis is stable (Fig. S4, Supporting Information). However, as the period increases, the minimum distance between beam pairs also increases. In this case, the number of cross-correlations involved in the calculation decreases, resulting in a slight increase in the corresponding uncertainty of phase velocity, particularly for periods greater than 3.6 s (Fig. S4e and f, Supporting Information).

3.3 Shear-wave velocity structure

Fig. 7 presents the horizontal slices of results from shear-wave velocity models at depths of 1, 2, 3, 4, 6 and 8 km. At a depth of ∼1 km, which is represented by Fig. 7(a), strong low-velocity shear-wave anomalies are observed to be laterally distributed beneath Heyuan Basin and Xinfengjiang Reservoir, and high-velocity anomalies are observed to correlate with the surface mountainous zones surrounding the reservoir. At depths of 2–4 km, a consistent low-velocity anomaly observed within this range is typical of sediments which characterize the Heyuan Basin (Figs 7b–d). Similarly, low-velocity zones are observed beneath Xinfengjiang Reservoir and decrease further as they approach ∼4 km (Figs 7c and d). At deeper depths of about 6–8 km, high-velocity bodies are distributed in the northwest, extending from the Heyuan Basin to Xichang (Fig. 7f). Furthermore, a total of 1528 earthquakes occurring between 2012 and 2015 were relocated using the HypoDD method (Waldhauser & Ellsworth 2000) in a previous study (He et al. 2018). We subsequently overlaid these earthquakes on to the horizontal velocity-depth slices (Fig. 7). The results indicate that earthquakes located in Xichang and the dam are primarily concentrated at the boundaries between the low- and high-velocity zones at depths ranging from 6 to 8 km (Figs 7e and f), whereas a small number of earthquakes are distributed within the shallow low-velocity zone beneath the centre of the reservoir (Figs 7c–e).

Fig. 8 presents four cross-sections: AA’, BB’, CC’ and DD’ of the absolute and relative shear-wave velocity models, which traverse the study area as shown in Fig. 1. The profiles AA’, BB’ and CC’ are oriented NE and transect the earthquake clusters in Xichang, Wudong, and the dam area, while profile DD’ is oriented NW and is projected to traverse through the three earthquake clusters (Fig. 1). At shallow depths (<4 km), all depth profiles reveal a wide distribution of low-velocity anomalies. Particularly notable at this depth is a thicker low-velocity layer beneath the Heyuan Basin, with a maximum thickness of up to 3 km. The low-velocity anomalies beneath Xichang and the main reservoir section have a thickness of approximately 1 km (Figs 8g and h). High-velocity zones correspond to structures beneath the hills, where velocities are ∼2 per cent greater than those in the surrounding areas (Figs 8b, d and f). The relative velocity model beneath the Heyuan Basin shows significant low anomalies (−5 to −8 per cent) that extends from the surface to a depth of ∼8 km. At deeper depths (>4 km), three distinct zones of high-velocity anomalies are observed at depths range from 6 to 11 km, beneath the reservoir. As illustrated in the relative velocity models (Fig. 8h), the high-velocity anomalies beneath Xichang and Heyuan faults are the most significant, with velocities 8 per cent higher than those of the surrounding rocks. These anomalies are overlain by significantly low-velocity anomalies (Fig. 8h). The seismogenic zone, shown to be located beneath Xichang and the dam, is located mainly at the edge of the high-velocity bodies, with depths greater than 6 km. Two M 4.8 earthquakes recorded near Xichang on 2012 February 16 (ID1) and 2013 February 22 (ID2), and the M 4.3 earthquake on 2023 February 11 (ID3) near the dam, support the identification of this seismogenic zone. On the other hand, earthquakes recorded directly beneath the reservoir are predominantly situated within the low-velocity zones that overlie the high-velocity zone.

4 DISCUSSION

Our model indicates that the upper-crust beneath the reservoir is characterized by relatively high velocity, while the Heyuan Basin exhibits low velocity, consistent with previous studies (Ye et al. 2017; He et al. 2018; Dong et al. 2022). Due to improved station distribution and a denser array (Fig. 1a), our models achieve higher resolution and reveal more detailed structures underlying the Xinfengjiang Reservoir (e.g. Fig. S5, Supporting Information). At shallow depths, the distribution and thickness of sediments in the study area are well imaged (Figs 7a–c). Although our model does not include a specific sediment-bedrock velocity jump, a transition velocity of 2.8–2.9 km s⁻¹ is reasonable (Xiong et al. 2024). Therefore, the maximum depth of the Heyuan Basin reaches 3 km (Figs 8e and g), which is deeper than the values obtained from ambient noise tomography using group velocity (Dong et al. 2022). In contrast to previous studies that proposed a single bedrock unit (Dong et al. 2022), our models reveal a strong lateral heterogeneity of the granitic bedrock, with three distinct high-velocity bodies located beneath the dam, Xichang, and the centre of the reservoir (Fig. 8g, h). Notably, earthquakes are predominantly concentrated at the edges of the high-velocity bodies rather than within them (Figs 8g and h).

Structural features indicate lateral heterogeneity within the internal structure of the granitic bodies (Fig. 8), with high-velocity anomalies representing intact granitic bedrock, while relatively low-velocity bodies indicate highly fractured granitic structures (e.g. Moos & Zoback 1983; Lees & Nicholson 1993; Clarke & Burbank 2011). The intact granitic bodies are sufficiently robust to facilitate stress accumulation (Tenthorey et al. 2003). The Xinfengjiang Reservoir is situated in a northeast–southwest compressional state (Ding et al. 1983), and such tectonic properties predispose the high-velocity bodies to significant stress, thereby becoming potential nucleation points for earthquakes. Earthquakes beneath the reservoir are mainly concentrated at the boundaries of the two high-velocity bodies, indicating their influence on the location of earthquakes (Fig. 8).

The velocity model shows there are strip-shaped low-velocity bodies overlying the high-velocity structures, particularly the inclined low-velocity strip beneath the HYF, which extends from the surface of the Heyuan Basin to a depth of ∼8 km beneath the RZSF (Figs 8g and h). These low-velocity strips are 7 per cent lower relative to the surrounding rocks. This sharp velocity contrast has been attributed to the presence of fluids (Chen et al. 2021), suggesting that the low-velocity anomalies may act as channels for downward fluid migration to deeper depths. We propose that the large volume of water impoundment in the Xinfengjiang Reservoir creates substantial downward pressure, facilitating fluid infiltration along fractures associated with the Heyuan fault, extending to a depth range of 6–8 km. Fluid diffusion within these structures leads to an increase in pore pressure at the boundaries of the high-velocity body, rendering the fracture susceptible to sliding and thus inducing earthquakes. Similar mechanisms underneath Xichang may also induce significant earthquakes.

At large dam sites, such as the Koyna Reservoir in India and the Aswan Reservoir in Egypt, dominant seismicity tends to occur after the water level has reached its peak, despite rapid increases in pore pressure as the water level rises to the peak following initial impoundment (Simpson et al. 1988; Gahalaut & Hassoup 2012; Gupta 2022). Similarly, the Xinfengjiang Ms 6.1 mainshock occurred about 6 months after the water level peaked at about 113.8 m (Gupta 2022). The observed time delay between the seismicity and peak water level is attributed to the time required for the pore pressure to diffuse from the reservoir to the hypocentres (Parotidis et al. 2003; Talwani et al. 2007). Moreover, this time delay varies considerably across different locations. For example, two M > 4 earthquakes were recorded at the dam site (ID3: M 4.3) on 2023 February 11, and in Wudong (ID4: M 4.5) on 2023 March 8, following the reservoir's peak water level of 113.3 m reached earlier in 2022 July. From the equation presented by Talwani et al. (2007):

where c is the hydraulic diffusivity and r is the distance. It has been demonstrated that, for two earthquakes with similar focal depths, the time delay $\tau$ is inversely related to the hydraulic diffusivity at different locations. Fig. 8 gives a pictorial illustration of the presence of a thick sedimentary layer near the dam. The low-velocity zones with high-fluid diffusion coefficients are also observed to extend directly towards the location of ID3. In contrast, a 2 km thick high-velocity layer overlays ID4 (Fig. 8h), indicating a more structurally intact rock with lower hydraulic diffusivity. This is suggested to account for the nearly one-month delay between the occurrences of ID4 and ID3. Therefore, it is important to consider the seismic time delays resulting from structural disparities in different areas of the reservoir.

Additionally, low-velocity anomalies are extensively distributed along the HYF, indicating fluid dispersion throughout the fault, potentially extending to greater depths (Figs 8e and h). Further evidence includes the abundant presence of hot springs and hot water wells along the HYF (Zhang et al. 2015; Qiu et al. 2018). High-velocity structures that serve as potential earthquake nucleation points are also located at the end of the HYF. However, earthquakes are concentrated along the middle section of the HYF in the study area, which is associated with a north-trending strike, while no earthquakes occur at the northern and southern edges of the fault (Fig. 1). Therefore, we propose that the fault's strike may influence the initiation of earthquakes in this area. According to Coulomb's criterion, fault rupture occurs when the Coulomb stress on the fault exceeds a threshold, such as the shear strength of the fault. In the study area, the maximum compressive horizontal stress is oriented in the northwest–southeast direction (Ding et al. 1983; He et al. 2018), approximately perpendicular to the NE-trending fault systems (Fig. 1). The angular difference between the maximum principal stress and the strike of the fault is relatively larger at the northern and southern edges of the HYF compared to its middle section (Fig. 1a). As a result, the resulting Coulomb stress is generally low at the northern and southern edges of the HYF. This may account for the variability in the distribution of earthquakes along the Heyuan fault. To adequately characterize the seismotectonic properties beneath the Xinfengjiang Reservoir, greater emphasis should be placed on detecting additional N-oriented faults.

5 CONCLUSIONS

In this study, we obtain a high-resolution S-wave velocity model of the upper crust underlying the Xinfengjiang Reservoir and surrounding areas, from newly acquired dense array ambient noise data using double beamforming method. Our observations reveal that several high-velocity structures are located beneath the Xinfengjiang Reservoir, with recorded earthquakes mainly occurring at the edges of these high-velocity regions, indicating their controlling role in seismic processes. Low-velocity zones, which serve as channels for fluid infiltration, are observed beneath Xichang and Heyuan faults and suggested to elevate the pore pressure at source depths, thereby inducing earthquakes. In addition, the structural differences at various locations beneath the Xinfengjiang Reservoir may have been contributed to the delays in the timing of the earthquakes. Given the seismicity along different sections of the Heyuan fault, we propose that more attention should be focused on N-oriented faults to adequately conduct the seismotectonic characterization of the Xinfengjiang Reservoir.

ACKNOWLEDGMENTS

This work was supported by the National Key Research and Development Program of China (Grant 2023YFB2905300), Shenzhen Sustainable Development Science and Technology Project (No. KCXFZ20201221173608023), NSFC (grants 42122027 and 41974052), SEA (grant SSKP202203) and the Institute of Earthquake Forecasting of China Earthquake administration (grant 2021IEF0601).

CONFLICT OF INTEREST

The authors acknowledge that there are no conflicts of interest recorded.

DATA AVAILABILITY

All plots are produced using the Generic Mapping Tools version 6.2.0 (Wessel et al. 2019) from https://mirrors.ustc.edu.cn/gmt/. The data used and produced in this study including the linearly stacked cross-correlations and inverted shear-wave velocity models are available at 10.5281/zenodo.10926703. All websites were last accessed in 2024 April. The supplemental material to this article includes 3 figures.

REFERENCES