1 Introduction

The Qujiang Fault is located in the southeastern part of the Sichuan-Yunnan Block and is one of the most seismically active regions in China (Kan et al. 1977; Wang et al. 2014). Situated on the southeastern edge of the Qinghai-Tibet Plateau, it is bounded by the stable South China Block to the east and the Indochina Block to the west. The collision of the Indian Plate with the Eurasian Plate has resulted in the compression of the Qinghai-Tibet Plateau, leading to the southeastward motion of the Sichuan-Yunnan Block (Deng et al. 2003; Wen et al. 2011). The total length of the Qujiang Fault is approximately 80 km, with a width of about 20 km, and it exhibits a general trend oriented N120°E (Zhang 1980; Liu et al. 1999). At the southeastern end of the fault, its orientation undergoes a directional change, forming a constrained bend before resuming its original strike (Fig. 1). The devastating Ms 7.7 Tonghai Earthquake on 5 January 1970 occurred on this fault, with the surface rupture extending over 48 km (Zhang and Liu 1978; Liu et al. 1999). This earthquake inflicted severe seismic disaster in the region, resulting in over 15,000 fatalities, nearly 30,000 injuries, with economic losses exceeding RMB 300 million yuan (at 1970 prices). Apart from the surface rupture along the fault trace, the Tonghai Earthquake triggered dozens of landslides, rockfalls, soil liquefaction, sand ejections, and other geological hazards (Liu et al. 1999). Numerous studies have indicated that the Qujiang Fault remains an active fault within the frequently active Sichuan-Yunnan Block, with a Quaternary right-lateral slip rate of approximately 3 mm per year (Wen et al. 2011; Wang et al. 2014). Historical records reveal that the Qujiang Fault experienced M 7 earthquakes in 1588 and 1913 (Liu et al. 1999). The most recent two M 7 earthquakes occurred more than 50 years apart, while it has been over half a century since the Tonghai Earthquake, which has once again drawn attention from the society (Yu et al. 2020). More importantly, the key concern for future seismic disaster assessment in the surrounding areas of active faults is not when the next destructive earthquake will occur, but rather the extent of damage and impact it will have on the society once it strikes. Therefore, it is imperative to conduct a seismic disaster loss and risk assessment of the Qujiang Fault. This assessment will provide scientific basis and data support for local government earthquake preparedness, disaster reduction, and emergency management planning. It will enable better preparation for future potential seismic disasters, reduce earthquake risks, and promote regional sustainable development.

Fig. 1
figure 1

Distribution of surface trace lines of the Qujiang Fault in western Yunnan Province, southwestern China. The beach ball represents the focal mechanism solution of the 1970 Ms 7.7 Tonghai Earthquake. The stars indicate the epicenter locations of six earthquake scenarios (red star is in the same position as the 1970 event). Background colors reflect the topography. The inset in the lower-left corner shows the location of the study area in China (blue rectangle)

Full size image

In traditional risk assessment, seismic hazard is calculated using the probabilistic seismic hazard assessment (PSHA) method (Cornell 1968; Li et al. 2020). This approach involves numerous assumptions to construct simplified mathematical models, aiming to approximate the true physical characteristics of earthquakes and ultimately calculate empirical ground shaking. These assumptions include: (1) Neglecting the three-dimensional (3D) complexity of seismic sources and replacing them with simple point, line, or planar sources (Xu and Gao 2007); (2) Ignoring the actual complex processes that generate seismic motions and using ground motion prediction equations (GMPEs) instead (Douglas and Edwards 2016; Gerstenberger et al. 2020). However, the data used to fit GMPEs are typically derived from moderate earthquakes at far distances, lacking sufficient data from near-fault events. This limitation makes it difficult for GMPEs to accurately account for the detailed effects of near-fault ground motions, which are often the primary areas of seismic damage. Moreover, GMPEs suffer from the “ergodic assumption” and fail to effectively consider “spatial correlation” between adjacent sites (Xin and Zhang 2021). The former leads to a lack of precision when predicting ground motions for specific regions based on multi-source data fitting (Infantino et al. 2020; Stupazzini et al. 2021), and the latter results in earthquake hazard and risk assessments that are typically underestimated (Erdik 2017; Chen and Baker 2019; Du and Ning 2021; Infantino et al. 2021).

In contrast, in recent years, deterministic seismic hazard analysis methods have rapidly developed, overcoming the issues of ergodic assumption and spatial correlation. These methods focus solely on investigating worst-case earthquake scenarios for specific seismic sources (Graves et al. 2011; Smerzini and Pitilakis 2018; Zhu 2018; Bradley 2019; Zhu and Yuan 2020; Riaño et al. 2021; Wang, Wang, et al. 2023). They take into account various factors, such as the fault geometric structure, rupture propagation directionality, heterogeneity of underground media, and site effects, using physical equations like elastic wave equation and friction law, and numerical simulation techniques to model the complete physical processes of source rupture, seismic wave propagation, and surface damage. As a result, these methods provide seismic ground motion distributions that better align with physical laws (Xin and Zhang 2021). However, most existing studies using deterministic approaches are limited to simulating ground motion distributions, with scarce cases that quantitatively assess population fatalities, economic losses, and earthquake risk under specific earthquake scenarios. Therefore, this study aimed to expand the current scope of research by integrating deterministic seismic hazard analysis methods with earthquake loss estimation methods to achieve seismic disaster loss and risk assessments under specified earthquake scenarios.

Specifically, in this study, we began by setting various physical parameters for earthquake simulations of the Qujiang Fault based on the research by Yu et al. (2020) on the recurrence of the 1970 Ms 7.7 Tonghai Earthquake. These parameters include the 3D complex fault geometry, regional stress field, friction law, media properties, and topography. Subsequently, we designed different nucleation locations and conducted a series of dynamic rupture simulations and strong ground motion simulations to obtain earthquake scenarios for the Qujiang Fault under different conditions. Next, based on the intensity distribution provided by different earthquake scenarios, we combined earthquake fatality estimation model and economic loss estimation model to quantitatively estimate the human and economic losses under the set earthquake scenarios. This assessment allowed us to determine the extent of damage and impact caused by various earthquake scenarios on the society. Furthermore, by analyzing the simulated earthquake motion parameters and human and economic loss data, we conducted a seismic risk analysis to identify areas at risk of seismic damage and areas likely to experience significant social impact. Through this seismic risk assessment, we aimed to gain a more scientific and in-depth understanding for potential seismic disaster characteristics posed by the Qujiang Fault in the future. This will provide essential insights for earthquake preparedness, disaster reduction, and emergency management decisions by local authorities.

2 Model Setup and Method

This section describes the parameter settings of the models and data sources, numerical simulation and disaster loss assessment methods, algorithm flow and the key equations, and validation of the models and methods.

2.1 Dynamic Rupture and Ground Motion Modeling Setup

To comprehensively understand the entire physical process of an earthquake, from fault rupture initiation to seismic wave propagation through the subsurface, and finally to the distribution of seismic ground motion, dynamic rupture simulation and strong ground motion simulation are essential. Various model parameters and physical properties need to be determined, including the complex fault geometry, regional stress field magnitude and direction, friction coefficients, 3D heterogeneous media, terrain, and earthquake nucleation locations. Based on field geological investigations (Zhang and Liu 1978; Zhu 1984, 1985; Wang et al. 2014), Yu et al. (2020) determined the 3D geometry of the Qujiang Fault while retaining the detailed variations of the fault strike as much as possible. The overall strike of the fault is approximately N120°E, with a length of about 77 km along its strike. The fault width is approximately 20 km, inferred from the distribution of aftershocks of the Tonghai Earthquake (Zhang 1980; Liu et al. 1999). Based on field geological investigations and numerical simulation experiments (Zhou et al. 1995; Wang et al. 2014), Yu et al. (2020) determined that the dip of the Qujiang Fault is relatively complex. The northwestern segment of the fault dips at an angle of 80° towards the southwest, the middle segment is vertical, and the southeastern segment dips at an angle of 75° towards the northeast. Moreover, the actual fault surface exhibits certain irregular roughness, and studies suggest that the roughness of the fault surface contributes significantly to the generation of high-frequency seismic waves during dynamic rupture (Andrews and Barall 2011; Shi and Day 2013). Therefore, the fault model considers the roughness of the fault surface, with a Hurst exponent (H) of 0.4 and a roughness amplitude of 200 m, following the work of Andrews and Barall (2011), Dunham et al. (2011), and Shi and Day (2013). Figure 2a displays the complex 3D fault model constructed in this study.

Fig. 2
figure 2

3D view of the fault geometry and stress distributions. a The nonplanar Qujiang Fault geometry. The six pentagrams represent the epicenters of the six earthquake scenarios (from left to right, they are recorded as Scenarios 1–6, respectively), with the red pentagram indicating the epicenter of the 1970 Ms 7.7 Tonghai Earthquake (Scenario 3). b The normal stress (Tn), and c shear stress (Ts) distributions on the fault plane. d The absolute value of the ratio of shear stress to normal stress \({\mu }_{0}\)

Full size image

The initial stress and friction parameters are consistent with Yu et al. (2020). Figures 2b and c show the distribution of normal and shear stresses obtained through the stress field resolution method. The magnitudes of the three axial stresses and the orientation of maximum principal compressive stress are listed in Table 1. In this study, we defined the absolute value of the ratio of shear stress to normal stress on the fault plane as \({\mu }_{0}\), which reveals the propensity for fault rupture (Zhang, Zhang, Xin, et al. 2020). A larger \({\mu }_{0}\) indicates a higher likelihood of fault rupture (Fig. 2d). We adopted the slip-weakening friction law (Ida 1972), with the parameters listed in Table 1. We employed an over-strength nucleation approach, setting the center of the nucleation region on the fault plane at a depth of 13 km (Fig. 2a) based on historical earthquake hypocenter depths. The radius of nucleation region is set to 2 km, and the initial shear stress within the nucleation region is slightly larger than the fault strength (0.1%) to trigger the rupture. Once the rupture starts within the nucleation region, it propagates to locations outside the nucleation region, and the stress conditions and frictional properties will determine whether the rupture will spontaneously propagate along the fault surface. See Yu et al. (2020) for more details.

Table 1 Parameter values set in the earthquake simulations
Full size table

In our simulations, we adopted the S-wave velocity (Vs) structure reference model for the Chinese region, which was obtained using the surface wave dispersion method by Shen et al. (2016). We also used the empirical relationship proposed by Brocher (2005) to derive the P-wave velocity and density. Additionally, we incorporated velocity structures from the sedimentary basins near the Qujiang Fault, including the Tonghai Basin and Quxi Basin (Yu et al. 2020), in which the minimum Vs velocity reaches 600 m/s. The spatial distribution of the Vs model (computational spatial extent) is shown in Fig. 3. The study area encompasses vast mountainous regions, and the undulating terrain significantly influences the propagation of seismic waves (Zhu et al. 2013; Khan et al. 2020). Therefore, to account for real topographic conditions, we used the Shuttle Radar Topography Mission (SRTM) terrain dataset (Reuter et al. 2007). This study did not consider the issues of visco-elastic and Q factor. Ultimately, the fault model was discretized into a grid of 1551×400 points, and the 3D spatial discretization for calculation was set at 200×1648×450 grid points, with a grid spacing of 50 m. The time step for the simulation was set to 0.003 s, with a total simulation time of 30 s. For the strong ground motion simulation, seismic wave propagation calculations were carried out at a larger spatial scale. Due to computational resource constraints, the grid spacing was set to 300 m, and the calculation spatial grid size was 1800×1400×300. The time step was 0.01 s, and the total simulation time was 150 s. All parameters used in the dynamic rupture simulation and strong ground motion simulation are listed in Table 1. Eight Nvidia Tesla A100 graphic processing units (GPUs) were used to conduct dynamic rupture and ground motion simulations.

Fig. 3
figure 3

Spatial distribution of the S-wave velocity (Vs). The area range of the model displayed is the computational space for the numerical simulation of strong ground motion in this study

Full size image

2.2 Numerical Simulation Method and Disaster Loss Assessment Method

The earthquake simulation method employed in this study is the curved-grid finite-difference method (CG-FDM), originally proposed by Zhang and Chen (2006) to address seismic wave propagation on irregular surfaces. Subsequently, Zhang et al. (2014) incorporated it into dynamic rupture simulations on irregular fault surfaces. A series of earthquake simulations have confirmed the efficiency and stability of this method in simulating non-planar fault dynamic ruptures and seismic wave propagation in complex media (Wang, Li, et al. 2022; Wang, Li et al. 2023; Xu et al. 2023; Yu et al. 2023). We employed the improved GPU-based CG-FDM dynamic rupture program developed by Zhang, Zhang, Li, et al. (2020) to simulate the dynamic rupture process. We used the GPU-accelerated seismic disaster rapid response platform CGFDM3D-EQR developed by Wang, Zhang, et al. (2022) to simulate the seismic wave propagation process.

Due to the limitation of computing resources and the lack of high-resolution media models, the maximum frequency calculated based on numerical methods is usually low. The calculation formula of the maximum frequency is \(f_{\max } = {{V_{\min } } \mathord{\left/ {\vphantom {{V_{\min } } {\left( {ppw \cdot dh} \right)}}} \right. \kern-0pt} {\left( {ppw \cdot dh} \right)}}\). In this study, the Vmin = 600 m/s. grid size dh = 300 m, points per wavelength (ppw) = 6. The maximum computational frequency of this study is 0.33 Hz. Although engineering applications require broadband ground motions, Chen et al. (2001) indicated that due to the uncertainty of building-based inventory methods in earthquake loss assessment, the results estimated of empirical methods based on macroscopic parameters usually fall between the upper and lower limits of the inventory method. The loss estimation based on individual buildings and broadband seismic motion time history curves is indeed more accurate, but it requires a detailed database and a large number of engineering experiments, which are beyond the scope of this study. Therefore, the empirical methods based on macroscopic parameters were selected in this study.

For earthquake disaster loss estimation, we used the earthquake fatality estimation model and earthquake economic loss estimation model developed by Li et al. (2021, 2023) for China’s Mainland. These models enable the estimation of human and economic losses in the study area using reasonable seismic intensity distribution data, population exposure data, and fixed assets exposure data. The basic equations for earthquake disaster loss estimation are as follows:

$$E_{i} { = }\sum\limits_{Int}^{{}} {r(Int) \times \left( {\frac{{HDI_{regression - year} }}{{HDI_{event - year} }}} \right)^{h} \times {\text{per capita Exposure}}_{(Int,event - year)} \times {\text{PE}}_{{\left( {Int,event - year} \right)}} }$$
(1)
$$P\left( {a < {\text{loss}} \le b} \right) = \Phi \left[ {\frac{{\ln b - \ln (E_{i} )}}{\zeta }} \right] - \Phi \left[ {\frac{{\ln a - \ln (E_{i} )}}{\zeta }} \right]$$
(2)
$$\zeta { = }\sqrt {\frac{1}{N - 2}\sum\limits_{i = 1}^{N} {\left[ {\ln \left( {O_{i} } \right) - \mu_{{\ln O_{i} |\ln E_{i} }} } \right]^{2} } }$$
(3)

where, \(E_{i}\) represents human or economic losses estimated for an earthquake. \(r(Int)\) represents a human or economic loss ratio function based on seismic intensity (Int). The Int was derived from the relationship between peak ground velocity (PGV) and intensity in the Chinese Seismic Intensity Scale (SAC 2020). The formula is \(Int{ = }3.00\log_{10} ({\text{PGV}}) + 9.77\). The HDI, also known as the Human Development Index, is the ultimate indicator proposed by the United Nations Development Programme in 1990 for the comprehensive evaluation of the level of social development. In this study, we used the HDI term to adjust the loss ratio function at the time scale (building fortification, emergency preparedness level, people’s seismic awareness, and other social impacts). \({\text{per capita Exposure}}_{(Int,event - year)}\) represents the per capita exposure quantity at each seismic intensity in an earthquake, which can be ignored for human losses, and for economic losses, it is the per capita fixed asset value. \({\text{PE}}_{{\left( {Int,event - year} \right)}}\) represents the number of people exposed to each seismic intensity in an earthquake. Φ represents the normal cumulative distribution function. \(\ln \left( {E_{i} } \right)\) denotes the expectation of the distribution, which is calculated by the loss estimation model with Eq. 1. ζ denotes the standard deviation of the distribution, which represents the uncertainty of the loss estimation model. \(\mu_{{\ln O_{i} |\ln E_{i} }}\) denotes the predicted economic loss based on the linear function regressed from the recorded (\(O_{i}\)) and estimated (\(E_{i}\)) loss pairs of all earthquake events following the functional form of \(\ln \left( {O_{i} } \right) = a + b\ln \left( {E_{i} } \right)\). N is the total number of earthquakes. In this study, the Qujiang Fault belongs to the Qinghai-Tibet Plateau region model in Li et al. (2021, 2023). In this region, the standard deviation ζ and h corresponding to the population and economic loss assessment models are 1.8474, 1.0, and 1.4418, 2.0, respectively. The human and economic loss ratio functions are as follows:

$$\log \left( {r_{human} } \right) = - 10.365 + 0.892 \cdot Int$$
(4)
$$r_{{economic}} = \Phi \left[ {\frac{1}{{0.133}} \cdot \ln \left( {\frac{{Int}}{{7.630}}} \right)} \right]$$
(5)

where \(r_{human}\) and \(r_{economic}\) represent the human loss ratio and economic loss ratio, respectively. The specific process of constructing model parameters and HDI, normality test, and uncertainty analysis of parameters can be referred to Li et al. (2021, 2023).

The approach has been successfully applied to actual earthquake cases in China in recent years (Li et al. 2022; Wang, Li, et al. 2023). The population distribution data used in this study was obtained from the LandScan population spatial distribution dataset, developed by the Oak Ridge National Laboratory (ORNL) in 1997 to meet the urgent need for improved population assessment for various impact assessments. The dataset is updated every year and provides global population spatial distribution data in ESRI grid format, with a spatial resolution of 1 km (30″ × 30″) (Sims et al. 2022). For the economic loss estimation, we referenced fixed asset data at the prefecture level from the research outcomes of Wu et al. (2014).

It is worth noting that the estimated losses in this study are only statistically significant expected values and cannot be used as accurate prediction results. Therefore, we used a normal cumulative distribution function based on model expectation and standard deviation to calculate the probability of different loss intervals, thereby reflecting the reliability of the expected loss and relative loss scale. The purpose of this study was not to accurately predict potential losses for potential earthquake scenarios, but rather to characterize the scale of disaster losses and probability, which may have more practical significance for future regional earthquake defense.

2.3 Validation of the Models and Methods

The 3D numerical simulation methods (CG-FDM) for dynamic rupture and strong ground motion simulation have been verified by many benchmark models in the SCEC/USGS Spontaneous Rupture Code Verification Project (Harris et al. 2018). For the earthquake scenario simulations of Qujiang Fault, the physical parameters (see Table 1) used in this study are consistent with Yu et al. (2020), which successfully realizes the recurrence of the 1970 Ms 7.7 Tonghai Earthquake and is in good agreement with the observations. To further verify our models and methods, we estimated the human and economic losses based on the strong ground motion simulation results of the 1970 Ms 7.7 Tonghai Earthquake (Scenario 3 in this study). The estimation results (Fig. 4) show that the maximum probability interval of human and economic losses estimated by this method has a good match with the actual reported value (deaths: 15,621; economic losses: RMB 300–3800 million yuan) (Zhang et al. 2013). Therefore, our methods are applicable in the study area.

Fig. 4
figure 4

This study estimated the 1970 Ms 7.7 Tonghai Earthquake disaster losses. a Estimated earthquake death expectation and the probability of occurrence of different loss degrees; b Expected economic losses assessed and the probability of occurrence of different levels of loss (at 1970 RMB prices)

Full size image

3 Earthquake Scenario Simulation Results

This section describes the earthquake scenario simulation results, including dynamic rupture and strong ground motion simulations of the six earthquake scenarios, and analyzes the similarities and differences.

3.1 Dynamic Rupture Simulation Results

We designed six nucleation zone along the Qujiang Fault to explore different earthquake scenarios that may occur in the future. Among them, the third nucleation zone corresponds to the hypocenter of the 1970 Ms 7.7 Tonghai Earthquake, with the remaining nucleation zones extending from it towards both sides (see Figs. 1 and 2a). The six earthquake scenarios triggered by the six nucleation zones (marked with pentagrams) from left to right in Fig. 2a are defined as Scenarios 1 to 6. After conducting dynamic rupture simulations for the six earthquake scenarios, the slip distributions and rupture time contours (1 s interval) on the fault surface are depicted (Fig. 5). It can be observed that earthquakes nucleated on the northwest side of the fault can rupture the entire fault surface and can all have moment magnitudes exceeding Mw 7. Particularly, Scenarios 2–4 release the highest energy, with moment magnitudes reaching Mw 7.32 and a maximum slip exceeding 4 m. However, on the southeast side of the fault, the geometric morphology changes significantly for the nucleation regions of Scenario 5 and on the northwestern side of Scenario 6. The fault exhibits curvature, resulting in a region with lower \({\mu }_{0}\) between the nucleation regions of Scenarios 5 and 6 (see Fig. 2d). Consequently, earthquakes initiated from the southeast side of the Qujiang Fault are influenced by the complex geometric structure, restricting the rupture to the southeast end of the fault, and preventing it from propagating across the entire fault. As a result, the moment magnitudes in this region are noticeably smaller, with Mw 5.92 and Mw 6.34 for Scenarios 5 and 6, respectively.

Fig. 5
figure 5

Slip distributions and rupture time contours (gray line, 1 s interval) on the fault surface for the six earthquake scenarios. a to f correspond to earthquake scenarios 1 to 6. The epicenter locations correspond to the pentagrams from left to right in Fig. 2a

Full size image

3.2 Strong Ground Motion Simulation Results

Based on the dynamic rupture models obtained from the dynamic rupture simulations (Fig. 5), we conducted strong ground motion simulations, taking into account the topographical variations and 3D complex media. This allows us to calculate the entire process of seismic wave propagation and obtain the distribution of ground motion in the study area. The distributions of peak ground velocity (PGV) for earthquake scenarios 1–6 are shown in Fig. 6. From these simulation results, an intriguing phenomenon emerged: earthquakes with the same magnitude and similar slip distribution but different epicenters (Fig. 5a–d) lead to significantly varied ground motion patterns (Figs. 6a–d). This observation highlights the significant impact of the seismic source location and the rupture directivity on the distribution of ground motion. Among the earthquake scenarios, Scenario 3 exhibits the highest PGV at 1.90 m/s, followed by Scenario 1 at 1.74 m/s and Scenario 4 at 1.58 m/s. According to the relationship between ground motion and seismic intensity in the Chinese Seismic Intensity Scale (SAC 2020), these three scenarios reach a maximum seismic intensity XI. Scenario 2 reaches a maximum PGV of 1.35 m/s, corresponding to seismic intensity X. As the epicenter moves from the northwest to the southeast, the concentrated release of seismic energy shifts from the southeast to the northwest. However, for Scenarios 5 and 6 at the southeast end of the fault, the rupture scales are relatively small, leading to limited energy propagation to the free surface (Fig. 6e and f). Consequently, these scenarios result in maximum PGVs of only 0.13 m/s and 0.19 m/s, corresponding to seismic intensities VII and VIII, respectively.

Fig. 6
figure 6

Peak ground velocity distribution of the six earthquake scenarios. The white solid line represents the surface trace of the Qujiang Fault, while the white pentagrams indicate the epicenter locations. The background map shows the topography. The color bar is saturated at 1.41 m/s

Full size image

4 Seismic Disaster Loss Estimation and Risk Analysis

In this section, we present the disaster loss estimation (including fatalities and economic losses) and seismic risk analysis for all potential earthquake scenarios on the Qujiang Fault. We used the Chinese Seismic Intensity Scale (SAC 2020) to convert the seismic ground motion into seismic intensity distributions for each earthquake scenario, as shown in Fig. 7. We estimated fatalities and economic losses for each earthquake scenario in every prefecture-level city within the study area (located in the upper part of each subplot in Fig. 7). We also estimated overall losses and probability distribution of loss ranges for each earthquake scenario (located in the lower part of each subplot in Fig. 7, with the left side indicating fatality estimation and the right side indicating economic loss estimation).

Fig. 7
figure 7

Intensity distribution and human and economic loss assessment (at current RMB prices) distributions for prefecture-level cities under the six earthquake scenarios (upper part of each subplot), as well as the overall human and economic loss values and probabilistic distributions for each earthquake scenario (lower part of each subplot, with the left subplot indicating human fatalities and the right subplot indicating economic losses). Subplots (a) to (f) correspond to earthquake scenarios 1 to 6

Full size image

The six earthquake scenarios exhibit significant differences in the damage caused to the society. Scenarios 5 and 6 result in limited impacts, causing minor losses to Yuxi City and Honghe Hani and Yi Autonomous Prefecture (hereafter referred to as Honghe). Conversely, Scenarios 1–4 may lead to severe social and economic destruction in the study area. Among these scenarios, Scenario 1 results in the most severe damage. The maximum probable fatality range is 10,000–100,000 people (with 41.2% occurrence probability), and the mathematical expectation of fatalities estimated is 11,694 people (Fig. 7a). Among the affected cities, Honghe is the most heavily affected, with potential fatalities exceeding 6,515, accounting for more than 55% of the total deaths under the earthquake scenario. Yuxi and Wenshan Zhuang and Miao Autonomous Prefecture (hereafter referred to as Wenshan) follow, potentially causing around 2800 and 2200 deaths, respectively. These three cities account for over 97% of the total human losses under this earthquake scenario, which is closely related to the fault rupture directivity. In terms of economic losses for Scenario 1, the maximum probability of economic loss range in the study area is RMB 100 to 1000 billion yuanFootnote 1 (with 46.3% occurrence probability), with an estimated value of 912 billion yuan (Fig. 7a). Among the affected regions, Honghe suffers the most severe economic losses, with a potential loss of 377.55 billion yuan, followed by Wenshan with 286.04 billion yuan and Yuxi with 197.30 billion yuan. It is worth noting that while Wenshan’ human fatality is significantly lower than Honghe and Yuxi, its economic losses are comparable to Honghe and significantly higher than Yuxi. This discrepancy may be attributed to the fact that much of Wenshan falls under the coverage of moderate to strong seismic intensity, indicating that the correlation between economic losses and fatalities may not be entirely positive.

Under Scenarios 1 and 2, despite the epicenters being located within Yuxi, the directional rupture and propagation of the earthquake result in Honghe, to the east of the epicenter, being the most severely affected area. This distribution is similar to the pattern observed in the 2008 Ms 8.0 Wenchuan Earthquake, where the epicenter and the most severely affected area were not the same (Xu et al. 2009). This reveals the critical defense areas that are often overlooked during the disaster prevention and mitigation process for active faults. In this case, the focus should be on Dali Bai Autonomous Prefecture (hereafter referred to as Dali) and Chuxiong Yi Autonomous Prefecture (hereafter referred to as Chuxiong) to the northwest of the fault, and Honghe and Wenshan to the southeast, with Honghe being of particular importance.

Scenarios 1 to 4 all result in earthquakes with Mw 7.3. Although the total economic losses they incur are comparable, there exists a significant difference in the number of fatalities. The reason for the above phenomenon is the exceptional sensitivity of fatality in Honghe to the distribution of ground motion. There is a significant variation in fatalities across different earthquake scenarios, ranging from 166 to 6515 people. Under Scenarios 1 and 2, the fatalities in Honghe account for over half of the total fatalities under each scenario, respectively. Therefore, if future disaster prevention measures in Honghe could address the seismic hazards posed by the aforementioned earthquake scenarios, it could substantially reduce the seismic risk brought about by the Qujiang Fault.

To better assess the earthquake risk caused by the Qujiang Fault in the region and provide a scientific reference for seismic fortification and emergency management planning based on future extreme earthquake scenarios, we calculated two earthquake risk scenarios based on seismic intensity. The first scenario (Fig. 8a) considers no seismic fortification cost and synthesizes the maximum ground shaking intensity from all earthquake scenarios to obtain the distribution of the maximum ground shaking intensity at each grid point on the surface in the study area. This represents the potential distribution of the maximum seismic intensity that the Qujiang Fault may cause in the study area in the future. Intensity levels VII to XI are concentrated along the NW-SE direction in the study area, with intensity levels X–XI focused near the surface trace of the Qujiang Fault, and the maximum intensity level XI occurs in Yuxi and Honghe. Under the second scenario (Fig. 8b), seismic fortification costs are considered, and we select the maximum seismic intensity value that occurs with a probability exceeding 50% at each grid point based on the results of all earthquake scenarios. This provides the final intensity risk distribution in the study area. The overall intensity range under the second scenario is slightly narrower compared to the first scenario, especially with the maximum intensity level XI only occurring within Yuxi. Additionally, Dali and Wenshan are not exposed to the seismic risk of intensity IX under this second scenario.

Fig. 8
figure 8

Seismic risk scenarios based on intensity in the study area. a The intensity at each grid point in the study area is represented by the maximum seismic intensity value across all seismic scenarios. b The intensity at each grid point in the study area is represented by the maximum seismic intensity value among scenarios with occurrence probabilities exceeding 50%

Full size image

Overall, the seismic risk in the intensity X–XI zones is very high. Taking unreinforced brick-masonry structures as an example, according to the descriptions in the Chinese Seismic Intensity Scale (SAC 2020), more than 60–90% of buildings in this zone would be destroyed (that is, beyond repair). These extreme high-risk zones are mainly concentrated near the surface trace of the Qujiang Fault, extending about 30 km southeastwards along the fault towards the eastern end in Honghe and about 60 km southwestwards in Yuxi, causing a large number of human and economic losses in Yuxi and Honghe. The seismic risk level in the intensity IX zones is high. Pedestrians may fall, and 40–70% of buildings may suffer moderate to severe damage (that is, requiring general to difficult repairs). Some areas in Chuxiong, Yuxi, and Honghe fall into this seismic risk level. In addition, the areas of Wenshan and Dali in the second seismic risk scenario are no longer classified under this high-risk level as compared to the first risk scenario. The areas in the intensity VIII range are at a moderate seismic risk level, where people may experience swaying, and 40%–70% of buildings may sustain minor to moderate damage (that is, requiring simple to general repairs). Areas within the intensity VII range are at a low seismic risk level, with 40–70% of buildings remaining basically intact or sustaining only minor damage. Areas below intensity level VI are at a very low seismic risk level, where more than 60–90% of buildings would remain basically intact.

5 Conclusion

In this study, we employed physical methods to successfully simulate the entire seismic process of the Qujiang Fault, from dynamic rupture to the distribution of ground shaking intensity. This led to the production of a series of earthquake scenario simulation results based on the Qujiang Fault. Furthermore, we assessed earthquake-related human and economic losses within the study area. As a result, we obtained preliminary earthquake risk assessment results that are specifically applicable to the Qujiang Fault region, revealing high-risk earthquake zones.

If a future earthquake nucleated in the northwest segment of the Qujiang Fault (covering approximately two-thirds of the total fault length), it is highly likely that earthquake scenarios with Mw 7 or higher will occur. The maximum seismic intensity of these earthquake scenarios could reach X–XI, resulting in significant earthquake disaster losses, with economic losses estimated at around RMB 1 trillion yuan. The human losses will vary significantly depending on the epicenter, but all are expected to be in the range of over 5000 fatalities. Conversely, if an earthquake nucleated in the southeast segment of the Qujiang Fault (covering approximately one-third of the total fault length), the risk is lower. Due to the physical properties such as fault geometry and stress distribution, the rupture may not propagate throughout the entire fault, resulting in an earthquake magnitude possibly not exceeding Mw 6.5. The maximum seismic intensity would be limited to VII–VIII, with fatalities in the range of about 10 people and economic losses estimated between RMB 10 and 20 billion yuan. The simulation and loss assessment results for various earthquake scenarios show that earthquakes of the same magnitude can result in significantly different distributions of ground shaking intensity and disaster losses, emphasizing the pivotal influence of seismic source location and fault rupture directivity on surface seismic disasters.

In the study area, regions with intermediate to high seismic risks are mainly concentrated in the NW-SE direction, following the Qujiang Fault. The most vulnerable cities include Dali, Chuxiong, Yuxi, Honghe, and Wenshan. Among them, certain areas in Yuxi and Honghe are at extremely high seismic risk, making them the critical focus for future seismic mitigation efforts related to the Qujiang Fault. If a strong earthquake occurs, regardless of the earthquake scenario, the scale of disaster losses in Yuxi remains relatively consistent, with around 2000 fatalities and economic losses of around  RMB 200 billion yuan. However, the disaster losses in Honghe would exhibit significant variations due to the rupture directivity. In extreme scenarios, the seismic risk in Honghe is significantly higher than in Yuxi, resulting in greater human and economic losses. Similar to the disaster pattern observed in the 2008 Ms 8.0 Wenchuan Earthquake, even though the Qujiang Fault is located within Yuxi, it is highly likely that Honghe would suffer the most severe societal damages. Implementing effective seismic retrofitting and emergency management planning for Yuxi and Honghe could significantly mitigate the seismic risk posed by the Qujiang Fault to the society in the future.