The influence of aftershocks on seismic hazard analysis: a case study from Xichang and the surrounding areas

In some instances, a strong aftershock can cause more damage than the mainshock. Ignoring the influence of aftershocks may lead to the underestimation of the seismic hazard of some areas. Taking Xichang and its surrounding areas as an example and based on the Seismic ground motion parameters zonation map of China (GB 18306-2015), this study used the Monte Carlo method to simulate synthetic mainshock sequences. Additionally, the Omi–Reasenberg–Jones (Omi–R–J) aftershock activity model is used to simulate the aftershock sequences that follow mainshocks above a certain magnitude threshold. Then, the mainshock and the aftershocks are combined to calculate the regional seismic hazard using ground motion prediction equations (GMPEs). Finally, the influence of aftershocks on seismic hazard analysis is examined and considered. The results show that in areas with moderate to strong seismic backgrounds, the influence of aftershocks on probabilistic seismic hazard analysis can exceed 50 %. These results suggest that the impact of aftershocks should be properly considered for future probabilistic seismic hazard analyses, especially in areas with moderate to strong seismic activity backgrounds and in areas prone to secondary disasters such as landslides and mudslides.

1 Introduction

Aftershocks are commonly removed from observed earthquake catalogs during probabilistic seismic hazard analyses, assuming that mainshocks follow a Poisson distribution. However, the strong aftershocks that follow an earthquake may cause more damage than the mainshock and should not be underestimated. As there is not enough time to repair damage between the mainshock and subsequent aftershocks, buildings can suffer cumulative damage from aftershocks, which can lead to additional casualties and property losses (Bi et al., 2022). For example, after the 1976 M_7.8 Tangshan earthquake in China, most houses in the aftershock zone collapsed, and the railway lines on the deck of a local bridge were damaged during the M_7.1 and M_6.9 aftershocks (Lv et al., 2007). The M_5.0 aftershock that followed the 2003 M_8.0 Hokkaidō earthquake in Japan caused a secondary fire disaster due to a spilled tank (Zhao et al., 2005). The M_6.3 aftershock that followed the 2010 M_7.1 Canterbury earthquake in New Zealand caused damage to buildings, 146 deaths, and over 300 people to go missing (Zhang et al., 2011). Lv et al. (2007) statistically analyzed the aftershocks that followed 21 M > 7.0 mainshocks in China and found that the proportion of peak ground accelerations caused by aftershocks that exceeded that of the mainshock was 76.2 %; that is, aftershocks may cause more severe damage than the mainshock. Therefore, ignoring the impact of aftershocks may lead to underestimation of the seismic risk in some areas. The cumulative damage-induced losses caused by strong aftershocks have attracted considerable attention in the field of disaster and catastrophe insurance modeling (Xiong, 2019).

Cornell (1968) proposed the classical probabilistic seismic hazard analysis (PSHA) method, and based on that work, Wiemer (2000) proposed the aftershock probabilistic seismic hazard analysis (APSHA). Gallovič and Brokešová (2008) combined the generalized form of the Omori law (Omori, 1894; Utsu, 1961; Utsu et al., 1995) that was given by Shcherbakov et al. (2004), refined the APSHA steps and parameterizations, and analyzed the seismic hazard probability of aftershocks that followed several earthquakes as case studies. Shen and Yang (2018) used the APSHA method established by Gallovič and Brokešová (2008) to analyze the aftershock seismic hazard probability of the 2017 M_7.0 Jiuzhaigou earthquake in China. In addition, many scholars have derived the influence of aftershocks on seismic hazard analysis using analytical solutions (Yeo and Cornell, 2009; Marzocchi and Taroni, 2014; Iervolino et al., 2014; Davoudi et al., 2020; Taroni and Akinci, 2021). Boyd (2012) and Xu and Wu (2017) used the epidemic-type aftershock sequence (ETAS; Ogata, 1988, 1998) model to generate catalogs with and without aftershocks. They used a spatially smooth seismicity model to calculate the impact of aftershock clusters for probabilistic seismic hazard analysis. Canales and Baan (2020) used the Poisson model to generate mainshock sequences and the ETAS model to generate aftershock sequences. Field et al. (2022) used the Third Uniform California Earthquake Rupture Forecast (UCERF3) ETAS model (UCERF3-ETAS) to evaluate the effects of declustering and Poisson assumptions on seismic hazard estimates. Wang et al. (2022) compared the ETAS-simulated hazard with approximations based on the declustered Poisson approach (DP), the nondeclustered Poisson approach (NDP), and the recently proposed sequence-based PSHA (Iervolino et al., 2014).

Based on the Seismic ground motion parameters zonation map of China (GB 18306-2015), this study uses the Monte Carlo method to simulate synthetic mainshock sequences. Then, according to the magnitude of the mainshock, the Omi–Reasenberg–Jones (Omi–R–J) aftershock activity model (Omi, 2013, 2016, 2019) is used to simulate the aftershock sequences that follow mainshocks for a certain magnitude threshold. Finally, the mainshocks and the aftershocks are combined to calculate the regional seismic hazard using ground motion prediction equations (GMPEs). Thus, the influence of aftershocks on seismic hazard analysis is analyzed.

Figure 1 Distribution of seismic events and the tectonic background in Xichang and its surrounding areas.

The city of Xichang, one of the three major space launch facilities in China, is located in the Anning River valley in the southwest of Sichuan Province. The Anning River fault and the Zemu River fault run through the city. Historically, three M ≥ 7.0 earthquakes have occurred in the region: a M_7.0 event in 814, a M_7.5 earthquake in 1536, and a M_7.5 event in 1850 (Fig. 1). The Anning River fault is one of the main faults in the north–south Sichuan–Yunnan tectonic belt and is also an important fault in southwest China. According to regional geological data (Li, 1993; He and Ikeda, 2007), the Anning River fault zone is the boundary of different tectonic units with Paleozoic to Mesozoic ages. The west side of the fault contains a metamorphic complex and magmatic rock belt, and a Mesozoic–Cenozoic sedimentary basin lies on the east side. The Zemu River fault has been active throughout the Holocene (Li and Wang, 1977; Du, 2000) and is connected to the Anning River fault zone in the north and the Xiaojiang fault zone in the south. The fault has an overall strike of 330°, a fault plane dip angle of more than 60°, and a dip direction of northeast or southwest. Since the late Quaternary, the Anning River fault and Zemu River fault have been characterized by continuous strike-slip movements (Xu et al., 2003a, b). The Anning River fault and Zemu River fault are located at the boundary of the central Yunnan secondary block in the rhombus-shaped Sichuan–Yunnan block, which controls the focal positions of most nearby earthquakes with M≥7 (Lu et al., 2012).

Xichang is located in an area prone to strong earthquakes. Considering the impact of aftershocks in seismic hazard assessment, it is of critical importance to focus on fortifying areas subject to strong aftershocks, especially against landslides, debris flows, and other secondary geological disasters. However, these preparations require the development of accurate disaster prevention technologies.

2 Aftershock activity models and their parameters

2.1 Omi–R–J aftershock sequence model

After moderate or strong earthquakes, when direct information is available, the early activity characteristics of the aftershock sequences are used for sequence type determination (Jiang et al., 2007), strong aftershock prediction (Omi et al., 2013), and short-term aftershock probability prediction (Reasenberg and Jones, 1989; Gerstenberger et al., 2005). These characteristics have important scientific value and practical significance in earthquake relief, regional earthquake risk assessment, and understanding of the earthquake sequence itself. Reasenberg and Jones (1989) developed the R–J model to predict the occurrence rate of early aftershocks based on the Omori–Utsu formula (Omori, 1894; Utsu, 1961) and the Gutenberg–Richter (G–R) law (Gutenberg and Richter, 1944).

According to the R–J model, the aftershock intensity function, with a magnitude no less than M at time t in the earthquake sequence, can be expressed as

$$\begin{array}{}\text{(1)}& \mathit{\lambda }\left(t,M\right)={\displaystyle \frac{k}{(t+c{)}^p}}{\mathrm{10}}^{-bM},\end{array}$$

where t is the time after the mainshock. The parameter k controls the overall aftershock productivity. The parameter p represents the decay degree of the seismic sequence. Parameter c is used to adjust the incompleteness of the aftershock records within a very short time after the mainshock. This parameter is a positive and small constant and is negatively correlated with the focal depth (Shebalin and Narteau, 2017). Parameter b represents the stress accumulation level (Wiemer and Katsumata, 1999; Enescu et al., 2011). This model is simple in principle and suitable for estimating the parameters of moderate to strong earthquake sequences with simple decay laws. As a classical seismic sequence analysis method, it is widely used in aftershock prediction throughout the world and for earthquake hazard assessments by the Global Earthquake Model (GEM) project.

After the occurrence of moderate or strong earthquakes, many small aftershocks will be “submerged” in the early stage, resulting in a reduction in the completeness of the earthquake catalogs, making it difficult to assess many of the small earthquakes below the magnitude of completeness. Based on the R–J model, Omi et al. (2013) proposed the “Omi–R–J” model by considering the aftershocks below the magnitude of completeness during the early stage of the earthquake sequence in the model parameter fitting and in the aftershock occurrence rate prediction. Omi et al. (2013) used the expression of the detection rate function q(M) given by Ogata and Katsura (1993) (OK1993 model) to describe the detection rate of the incomplete part of the earthquake catalog. The actual recorded earthquake probability density function can be expressed as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}P\left(M\left|\mathit{\beta },\mathit{\mu },\mathit{\sigma }\right.\right)& ={\displaystyle \frac{e^{-\mathit{\beta }M}q\left(M\left|\mathit{\mu },\mathit{\sigma }\right.\right)}{{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{-\mathit{\beta }M}q\left(M\left|\mathit{\mu },\mathit{\sigma }\right.\right)dM}}\\ & ={\displaystyle \frac{e^{-\mathit{\beta }M}q\left(M\left|\mathit{\mu },\mathit{\sigma }\right.\right)}{e^{(-\mathit{\beta }\mathit{\mu }+{\mathit{\beta }}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}/\mathrm{2})}/\mathit{\beta }}}\\ & =\mathit{\beta }{e}^{-\mathit{\beta }(M-\mathit{\mu })+{\mathit{\beta }}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}/\mathrm{2}}q\left(M\left|\mathit{\mu },\mathit{\sigma }\right.\right),\end{array}\end{array}$$

where β is equal to bln10, μ represents the corresponding magnitude when the detection rate is 50 %, σ is the corresponding magnitude dispersion, and μ+2σ or μ+3σ is usually used to approximate the minimum magnitude of completeness M_c. In the parameter estimation of Eq. (2), the “state-space” model developed by Omi et al. (2013) was used to estimate the time-varying factor μ(t). Specifically, μ(t) is set as the discrete distribution function corresponding to the aftershock time sequence $t_i\le t\le t_{i+\mathrm{1}}(i=$ 1, 2, …, n). The hyperparameter V is set to control the smoothness of the distribution, assuming a priori the distribution with a smooth constraint on μ(t). After the parameters β,σ, and V are optimized and the maximum a posteriori estimation is performed by the maximum expectation (EM) iterative algorithm, the parameter $\mathit{\mu }=({\mathit{\mu }}_{\mathrm{1}},{\mathit{\mu }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\mu }}_n{)}^T$ is obtained by Bayesian estimation. In the early period after a mainshock, the waveforms of small earthquakes are submerged by the waveforms of large earthquakes, making it difficult to identify small earthquakes and resulting in a lack of a catalog of small earthquakes. The EM algorithm is based on the estimation of hyperparameters in the Newton iterative algorithm. It can optimize the parameters in the case of missing small earthquakes in the early period, reducing the error in the Newton iterative algorithm and obtaining more objective parameters. The earthquake detection rate function considering incomplete earthquake records can be expressed as $\mathit{\nu }(t,M)=\mathit{\lambda }(t,M)q\left(M\left|\mathit{\mu }\left(t\right),\mathit{\sigma }\right.\right)$. The logarithmic likelihood function related to parameters $p,c,$ and k is

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\ln L(k,c,p)& ={\sum }_{M_i\ge M_{\mathrm{c}}}\ln \mathit{\nu }(t,M)\\ & -{\int }_{M_{\mathrm{c}}}^{\mathrm{\infty }}\mathrm{d}M{\int }_{\mathrm{0}}^T\mathrm{d}t\mathit{\nu }(t,M),\end{array}\end{array}$$

where t_i and M_i are the time and magnitude of the ith aftershock that occurred within the “learning period” [0, T] during model fitting. Then, the parameters p, c, and k and the standard deviation are determined by combining the Omori–Utsu formula and the maximum likelihood method.

2.2 ETAS time series model

The ETAS model introduces the idea of self-similarity and assumes that both background earthquakes and triggered earthquakes can stimulate their own aftershocks and that many direct aftershocks and indirect aftershocks (aftershocks of aftershocks) can be generated after a mainshock. Therefore, the ETAS model is constructed with branch point process characteristics (Ogata, 1988; Bi and Jiang, 2019). The conditional intensity function can be expressed as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{\lambda }\left(t\right)& ={\mathit{\mu }}_{\mathrm{ETAS}}+K_{\mathrm{ETAS}}\sum _{t_i<t}{\displaystyle \frac{e^{{\mathit{\alpha }}_{\mathrm{ETAS}}(M_i-M_{\mathrm{0}})}}{(t-t_i+c_{\mathrm{ETAS}}{)}^{p_{\mathrm{ETAS}}}}},\\ & M_i>M_{\mathrm{0}},\end{array}\end{array}$$

where t−t_i represents the elapsed time of seismic event i and K_ETAS is a normalized constant that determines the expected number of aftershocks directly triggered by the M_i event. The parameter α_ETAS represents the ability of a seismic event to stimulate secondary aftershocks (Ogata, 1989, 1992). Compared with isolated earthquakes and main aftershocks, the α_ETAS of a swarm-type earthquake sequence is smaller, generally α_ETAS < 1 (Ogata, 2001), and p_ETAS represents the decay degree of the seismic sequence. Parameter c_ETAS is used to adjust the incompleteness of aftershock records within a very short time following the mainshock. Parameter μ_ETAS indicates the occurrence rate of background earthquakes. In the calculation process, when the occurrence rate of background earthquakes in the area is low, μ_ETAS= 0 is used to better ensure the stability of parameter fitting.

The maximum likelihood estimation (MLE) method is used to estimate the parameters {$K_{\mathrm{ETAS}},c_{\mathrm{ETAS}},{\mathit{\alpha }}_{\mathrm{ETAS}},p_{\mathrm{ETAS}}$} in the ETAS model. The likelihood function L is expressed as

$$\begin{array}{}\text{(5)}& \log L=\sum _{i:S\le t_i<T}\lg \mathit{\lambda }\left(t_i\right)-\underset{S}{\overset{T}{\int }}\mathit{\lambda }\left(t\right)\mathrm{d}t.\end{array}$$

2.3 Aftershock sequence models in the Xichang area

Gao (2015) divided the Chinese mainland and its adjacent areas into 29 seismic belts, of which 25 seismic belts are located inside mainland China. Since 1970, the Xianshuihe East–Yunnan seismic belt, where the Xichang area is located (see Fig. 2), has experienced six M ≥ 7.0 earthquakes; the 1970 M_7.8 Tonghai earthquake in Yunnan, 1973 M_7.6 Luhuo earthquake in Sichuan, 1973 M_7.2 Nima earthquake in Tibet, 1974 M_7.1 Daguan earthquake in Yunnan, 1997 M_7.4 Nima earthquake in Tibet, and 2010 M_7.1 Yushu earthquake in Qinghai. As early-seismic-monitoring ability in Tibet is limited and the number of recorded aftershocks is low, the two M_7.0+ earthquakes in Tibet cannot be used to fit aftershock parameters. We estimated the aftershock sequence parameters of the other four M_7.0+ earthquakes using the ETAS model and Omi–R–J model. The results are shown in Tables 1 and 2. To obtain more samples of aftershock sequence parameters, we use the Omi–R–J model to calculate the aftershock sequence parameters of 40 M_4.5–7.0 earthquakes. The results are also shown in Table 2.

Figure 2 The spatial distribution of the potential sources for the Xianshuihe East–Yunnan seismic belt where the Xichang area is located and six M≥7.0 earthquakes in the belt since 1970. Publisher's remark: please note that the above figure contains disputed territories.

Table 1 Basic information on four mainshocks with M≥7.0 and aftershock sequence parameters, as calculated by the ETAS model, from the Xianshuihe East–Yunnan seismic belt, where the Xichang area is located. The date format is year/month/day.

Table 2 Basic information on 44 mainshocks with M≥4.5 and aftershock sequence parameters, as calculated by the Omi–R–J model, from the Xianshuihe East–Yunnan seismic belt, where the Xichang area is located. The date is given in the format year/month/day.

3 Probabilistic seismic hazard analysis method considering aftershocks

3.1 Probabilistic seismic hazard analysis considering only the Poisson distribution

Wu et al. (2020) used the Monte Carlo method to simulate synthetic earthquake catalogs for probabilistic seismic hazard analysis based on the Seismic ground motion parameters zonation map of China (GB 18306-2015). The seismic source zone model used by this map is based on seismological and geological data for China. To reflect the heterogeneity in potential seismicity and describe the structural complexity more faithfully, the model adopts a three-level delineation of seismic belts, uses background and structural sources, and considers the tectonic differences between eastern and western China (Zhou et al., 2013). The spatial relationship of the three source levels is as follows (see Fig. 2): the base layer is the seismic belt (seismic statistical area), which is used to reflect the overall statistical characteristics of seismicity; the middle layer is the background potential sources, which are used to reflect the differences in seismic characteristics of small- and moderate-magnitude earthquakes under different tectonic conditions; and the upper layer consists of the structural potential sources, which are used to reflect the small-scale spatial seismic heterogeneity caused by the differences in local seismic structural conditions. This is a peculiar property of the seismicity model used for probabilistic seismic hazard assessment in China (CPSHA). Figure 2 shows the spatial distribution of the potential sources for the Xianshuihe East–Yunnan seismic belt where the Xichang area is located.

Figure 3 The calculation sites and the potential seismic sources in and around Xichang. The meaning of the color of the potential seismic source is the same as that in Fig. 2.

The seismic zoning map of China (GB 18306-2015) has established the corresponding probability model and spatial distribution model of earthquake occurrence and gives the basic parameters of each seismic zone. Figure 3 shows the potential seismic sources in and around Xichang.

Table 3 List of seismicity parameters of potential seismic sources in and around Xichang.

According to the basic assumptions and seismicity parameters of the zoning map (Table 3), the following steps are used to synthesize the sets of earthquake sequences (Wu and Gao, 2018; Wu et al., 2020).

1. Based on the assumption that the occurrence of earthquakes in seismic zones satisfies the Poisson distribution, the time length T of the simulated earthquake sequence and the average annual occurrence rate ν₄ of earthquakes with magnitude 4 and above in the seismic zone should be determined first. Then, a Poisson distribution random number L is generated with T and ν₄ as parameters, where L is the number of earthquakes in the seismic zone for the length of time T to be simulated.

2. Based on the assumption that the magnitude distribution of seismic zones satisfies the truncated Gutenberg–Richter relationship (magnitude–frequency relationship), with the minimum magnitude level M₀ and the maximum magnitude level M_UZ, the magnitudes of earthquakes to be simulated are determined.

   The magnitude–frequency relationship is represented as

   $$\begin{array}{}\text{(6)}& \log N=a-bM,\end{array}$$

   where a and b are coefficients, N is the number of earthquakes whose magnitude are equal to or greater than M, and the initial magnitude of the zoning map is 4. The cumulative number of earthquake events is

   $$\begin{array}{}\text{(7)}& N\left(M\right)=e^{a-bM}.\end{array}$$

   If we take ΔM=0.1, then

   $$\begin{array}{}\text{(8)}& N\left(M\right)>N(M+\mathrm{\Delta }M).\end{array}$$

   Based on M= 4.1, 4.2, 4.3, …, M_UZ, a random number u that satisfies a uniform distribution between 0 and 1 is generated. Then, the following is determined:

   $$\begin{array}{}\text{(9)}& u\in {\displaystyle \frac{N(M+\mathrm{\Delta }M)}{N\left(\mathrm{4}\right)}}\sim {\displaystyle \frac{N\left(M\right)}{N\left(\mathrm{4}\right)}}.\end{array}$$

   If the above formula is true, the magnitude M of an earthquake event is determined.

3. The epicenter location is determined as follows. First, the potential source area H where the earthquake is located should be determined. According to the magnitude M determined in the previous step, the magnitude range d to which the earthquake belongs is determined. Because the probability P_d(h) of each magnitude range in each potential source area is known, a random number u is generated that satisfies a uniform distribution between 0 and 1. The following is then determined:

   $$\begin{array}{}\text{(10)}& u\in \sum _{h=\mathrm{1}}^{H-\mathrm{1}}{P}_d\left(h\right)\sim \sum _{h=\mathrm{1}}^H{P}_d\left(h\right).\end{array}$$

   If so, the potential source area H where the earthquake event is located is determined. Based on the assumption that the epicenter is evenly distributed in the potential source area, a point is randomly selected in the potential source area H as the epicenter location of an earthquake.

4. According to the azimuth of the potential source area, the azimuth of the earthquake is determined.

At this point in the calculation, the basic elements of an earthquake have been determined. Steps (2)–(4) are repeated until the required number L of earthquakes in the seismic zone is obtained, accounting for all possible seismic zones that may affect the site and thus determining a seismic sequence and completing one sampling.

If the time length T is set to 1 year, the seismic sequence obtained by one sampling is called a 1-year seismic sequence in this paper. When the time length is set to 10 years, the sequence is called a 10-year earthquake sequence.

For each earthquake in seismic sequences, the peak ground acceleration (PGA) of each site is calculated by the optimal ellipse search algorithm through ground motion prediction equations (GMPEs).

For 5 000 000 simulations of a 1-year earthquake sequence, if a site is affected by ground motions exceeding specific values, the sequence is assigned a value of 1. The sum of earthquake sequences identified as 1 is counted and is divided by the total number of earthquake sequence simulations (i.e., 5 000 000), resulting in the annual exceedance probability of specific ground motions. Through the annual exceedance probability, the 50-year exceedance probabilities of 10 % and 2 % can be calculated.

This is how the probabilities of seismic hazard are obtained.

3.2 Probabilistic seismic hazard analysis considering aftershocks

Since there are only a few strong earthquakes with M ≥ 7 in the Xianshuihe East–Yunnan seismic belt, the Omi–R–J model is selected as the aftershock parameter model. According to the spatial division of the Seismic ground motion parameters zonation map of China (GB 18306-2015), the median values of the p, c, K, and b values (see Table 2) used in the Omi–R–J model (Omi et al., 2013, 2016, 2019) for the aftershock sequence samples from the Xianshuihe East–Yunnan seismic belt are 0.8747, 0.0187, 0.0133, and 0.8361, respectively. The aftershock sequences are generated according to these median values and the following steps.

The mainshock sequences are simulated by the Monte Carlo method based on the Seismic ground motion parameters zonation map of China (GB 18306-2015). Each synthetic sequence represents a 1-year possible distribution of earthquakes in the region that is consistent with the seismicity model (Wu et al., 2020). Considering the destructiveness of the earthquake, when the magnitude threshold for the mainshock is met (i.e., M ≥ 6.0 in this study considering a sufficiently large potential impact on the site, and the value can be adjusted as needed), the aftershock sequence sampling is started.

The minimum magnitude of the aftershock sequence is set to 4.0, and the maximum magnitude is equal to the magnitude of the mainshock. In fact, the magnitude of aftershocks can be greater than that of the mainshock. In this study, we focus on the “aftershocks”, so we adopted the assumption of Iervolino et al. (2014). That method assumes foreshocks do not contribute exceedances, aftershocks do not trigger their own aftershocks, and aftershocks are smaller than the mainshocks. The aftershock sequence satisfies the magnitude–frequency relationship $N\left(M\right)={\mathrm{10}}^{a-bM}$. The aftershock occurrence time t is within 30 d of the mainshock and follows the Omori–Utsu formula $N\left(t\right)=\frac{K}{(t+c{)}^p}$. The time interval between a strong aftershock and the mainshock varies from a few seconds to several years, but most strong aftershocks occur a few days or even a day after the main shock (JMA, 2009; Tahir et al., 2012). A length of 30 d is taken as the duration for a simplified calculation and can be changed as needed. According to the median values of p, c, k, and b and the upper limit of magnitude of the potential sources, the magnitude and time series of aftershocks with M ≥ 4 are simulated.

According to the empirical relationship between the magnitude of the mainshock and the rupture scale (Wells and Coppersmith, 1994), the rupture length and width are calculated by

$$\begin{array}{}\text{(11)}& {\displaystyle }L={\mathrm{10}}^{(-\mathrm{3.22}+\mathrm{0.69}{M}_{\mathrm{w}})},\text{(12)}& {\displaystyle }W={\mathrm{10}}^{(-\mathrm{1.01}+\mathrm{0.32}{M}_{\mathrm{w}})}.\end{array}$$

The rupture strike is taken as the direction of the potential source area where the mainshock is located, and the model of Felzer and Brodsky (2006) is adopted; that is, the aftershock density decays exponentially with increasing distance r from the fault, $\mathit{\rho }\left(r\right)=c{r}^{-n}$, where n is 1.37 and c is a constant. Thus, the locations of the aftershock epicenters can be determined.

The number of aftershocks. We have accounted for the number of M_4.0+ aftershocks for M_5.0+ mainshocks in the Chinese mainland and its surrounding area, and we have found that when the mainshock is greater than 6.0, the number of M_4.0+ aftershocks within a month (30 d) increases with the magnitude of the mainshock, yielding the following statistical relationship: ${\log }_{\mathrm{10}}\left(N\right)=\mathrm{0.84}M-\mathrm{4.57}$ (shown by the red line in Fig. 4). This relationship fluctuates within the range of ±0.8 (shown by the two dashed red lines) and obeys the normal distribution under linear coordinates. The number of aftershocks corresponding to a certain magnitude is generated according to this relationship.

Figure 4 The M_5.0+ mainshocks and the number of their M_4.0+ aftershocks for the Chinese mainland and its surrounding areas.

Figure 5a shows a schematic diagram of the spatial distribution of 1-year mainshocks with M≥6 sampled 100 times. The yellow star in the figure is the location of the mainshock. Figure 5b is a schematic diagram of the spatial distribution of the corresponding aftershocks. The small blue dot in the figure is the aftershock corresponding to the mainshock. The distribution direction of the aftershocks refers to the strike of the potential source area where the mainshock is located. In this study, considering the destructiveness of the earthquake, when the magnitude of the mainshock is ≥ 6.0, random sampling of the aftershock sequence begins, and the sampling time is set to within 30 d of the mainshock. The model program user interface can be used to adjust and refine the aftershock model to account for random aftershock sequences in the future that may have different requirements. To ensure the stability of the results, we conducted 5 million 1-year samplings.

Figure 5 Distribution of 1-year mainshocks with M≥6 and their corresponding aftershocks in 100 samplings: (a) mainshocks and (b) aftershocks. Publisher's remark: please note that the above figure contains disputed territories.

After the aftershocks are obtained, the mainshocks and aftershocks are combined, the ground motion value of the site is calculated using the ground motion prediction equations (GMPEs), and the exceedance probability for a specific case is counted; thus, a probabilistic seismic hazard analysis that considers aftershocks can be carried out. Figure 6 shows the calculation process for this analysis. According to the Seismic ground motion parameters zonation map of China (GB 18306-2015), the GMPEs of peak ground acceleration (PGA) suitable for the Xichang area are as follows (Xiao, 2011):

Figure 6 Flowchart of seismic hazard analysis when considering aftershocks.

When M < 6.5,

$$\begin{array}{}\text{(13a)}& \left\{\begin{array}{l}{\log }_{\mathrm{10}}{\stackrel{\mathrm{‾}}{G}}_l(M,R)=\mathrm{2.331}+\mathrm{0.646}M\\ -\mathrm{2.431}\log \left(R+\mathrm{2.647}\exp \left(\mathrm{0.366}M\right)\right)\\ {\log }_{\mathrm{10}}{\stackrel{\mathrm{‾}}{G}}_s(M,R)=\mathrm{1.017}+\mathrm{0.614}M\\ -\mathrm{1.866}\log \left(R+\mathrm{0.612}\exp \left(\mathrm{0.457}M\right)\right)\end{array}\right.,\end{array}$$

and when M≥6.5,

$$\begin{array}{}\text{(13b)}& \left\{\begin{array}{l}{\log }_{\mathrm{10}}{\stackrel{\mathrm{‾}}{G}}_l(M,R)=\mathrm{3.846}+\mathrm{0.413}M\\ -\mathrm{2.431}\log \left(R+\mathrm{2.647}\exp \left(\mathrm{0.366}M\right)\right)\\ {\log }_{\mathrm{10}}{\stackrel{\mathrm{‾}}{G}}_s(M,R)=\mathrm{2.499}+\mathrm{0.388}M\\ -\mathrm{1.866}\log \left(R+\mathrm{0.612}\exp \left(\mathrm{0.457}M\right)\right)\end{array}\right.,\end{array}$$

where G(M,R) is the peak ground acceleration (PGA), M is the magnitude, R is the epicentral distance, and the other coefficients are obtained by regression.

4 Influence of aftershocks on probabilistic seismic hazard analysis

To calculate the impact of aftershocks on seismic hazard analysis, Xichang and its surrounding areas were divided into a 0.1° × 0.1° grid (see Fig. 3), and the PGA values of the 50-year exceedance probability of 10 % and 2 % were calculated for each grid point. The results of the calculation with and without aftershocks were compared.

Figure 7 shows the PGA (gals) contour map of the 50-year exceedance probability of 10 % in Xichang and its surrounding areas calculated without and with aftershocks as well as the aftershock impact rate distribution map.

Figure 7 Comparison of the aftershock impacts on PGA (gals) with a 10 % exceedance probability over 50 years in Xichang and its surrounding areas. (a) PGA (gals) contour map of the 50-year exceedance probability of 10 % considering only the mainshocks. (b) PGA (gals) contour map of the 50-year exceedance probability of 10 % considering the mainshocks and aftershocks simultaneously. (c) Distribution map of the aftershock impact rate; aftershock impact rate = (calculation results of model with aftershocks − calculation results of model without aftershocks) $/$ calculation results of model without aftershocks.

To calculate the aftershock impact rate, we take the difference between the calculation results of the aftershock model and the calculation results of the model without aftershocks and divide that value by the calculation results of the model without aftershocks. That is

$$\text{impact rate}={\displaystyle \frac{\begin{array}{c}(\text{results with aftershocks}-\\ \mathrm{\text{results without aftershocks}})\end{array}}{\text{results without aftershocks}}}.$$

The maximum impact rate of aftershocks in Xichang and its surrounding areas is 55 %, the minimum is 0 %, and the average is 10 %. Aftershocks have the largest impact in the Xichang urban area, where there was a M₇ earthquake in 814, a M_7.5 earthquake in 1536, and a M_7.5 earthquake in 1850. The upper limit of magnitude of the potential source area is 8.0.

Figure 8 Comparison of the aftershock impacts on PGA (gals) with a 2 % exceedance probability over 50 years in Xichang and its surrounding areas. (a) PGA (gals) contour map of the 50-year exceedance probability of 2 % considering only the mainshocks. (b) PGA (gals) contour map of the 50-year exceedance probability of 2 % considering the mainshocks and aftershocks simultaneously. (c) Distribution map of the aftershock impact rate, aftershock impact rate = (calculation results of model with aftershocks − calculation results of model without aftershocks) $/$ calculation results of model without aftershocks.

Figure 8 shows the PGA (gals) contour map of the 50-year exceedance probability of 2 % in Xichang and its surrounding areas calculated with and without aftershocks. Additionally, this figure also shows the aftershock impact rate distribution map. The maximum impact rate of aftershocks in Xichang and its surrounding areas is 72 %, the minimum is 0 %, and the average is 10 %. The greatest impact of aftershocks is also in the Xichang urban area, where there was a M₇ earthquake in 814, a M_7.5 earthquake in 1536, and a M_7.5 earthquake in 1850. The upper limit of magnitude of the potential source area is 8.0. In this calculation, only mainshocks with M ≥ 6.0 generate aftershocks. Therefore, the calculated results are consistent with the aftershock model. The seismic hazards for sites with different seismic backgrounds are affected by aftershocks to different degrees.

5 Comparison with the ETAS model

The ETAS model is a widely used statistical method for capturing short-term spatiotemporal earthquake clustering. However, its application is occasionally impeded by the challenge of estimating a substantial number of unknown parameters. Recent advancements in ETAS formulations introduce spatial and temporal variability into certain parameters, further complicating their estimation process. Mancini and Marzocchi (2023) introduced a simple ETAS method called SimplETAS. The basic idea behind SimplETAS is that the earthquake clustering process in crustal regions is time- and space- independent, a premise substantiated by empirical analyses conducted by Stallone and Marzocchi (2019).

The functions adopted in SimplETAS are defined as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}\mathit{\lambda }& (t,x,y,m)=\left[\mathit{\nu }\mathit{\mu }\right(x,y)\\ & +\sum _{i:t_i<t}\mathit{\kappa }\left(m_i\right)g(t-t_i)f(x-x_i,y-y_i;m_i)\left]s\right(m),\end{array}\text{(15)}& {\displaystyle }\mathit{\kappa }\left(m_i\right)=A{e}^{\mathit{\alpha }(M_i-M_{min})},\text{(16)}& {\displaystyle }g(t-t_i)={\displaystyle \frac{c^{p-\mathrm{1}}(p-\mathrm{1})}{(t-t_i+c{)}^p}},\text{(17)}& {\displaystyle }\begin{array}{rl}f(x-x_i,y-y_i\left|m_i\right.)& ={\displaystyle \frac{q-\mathrm{1}}{\mathit{\pi }D{e}^{\mathit{\gamma }(M_i-M_{min})}}}\\ & [\mathrm{1}+{\displaystyle \frac{(x-x_i{)}^{\mathrm{2}}+(y-y_i{)}^{\mathrm{2}}}{D{e}^{\mathit{\gamma }(M_i-M_{min})}}}{]}^{-q},\end{array}\end{array}$$

in which A is the productivity; α is the coefficient of the exponential magnitude-dependent productivity law; c and p are the time constant and the exponent of the modified Omori law, respectively; q and D define the spatial distribution of triggered events; γ accounts for the correlation between the aftershock area and the magnitude of the triggering event; and t_i, x_i, and y_i are the temporal and spatial distances of the ith past earthquake from the present time t and from the considered location (x,y), respectively.

In the model, the seven parameters in the conventional ETAS formulation governing earthquake clustering, namely $\mathit{\{}\mathit{\alpha },p,c,D,\mathit{\gamma },q,\mathit{\beta }\mathit{\}}$, are predetermined. Only the total background seismicity rate (ν) and the seismic productivity (A), which exhibit significant variations depending on the region, remain to be estimated. The SimplETAS model can work as well as the ETAS model. Therefore, in this study, we follow the SimplETAS model and use its codes to simulate 10 000 sets of earthquake catalogs in Xichang and the surrounding areas for comparison analysis. Referring to Mancini and Marzocchi (2023), Table 4 shows the corresponding parameters. The background seismicity spatial probability density function (PDF) is shown in Fig. 9. When estimating the PDF, the M_4.0+ earthquake catalog from 1 January 1970 to 23 August 2023 is used. The primary catalog spans 1 January 1975 to 23 August 2023, while the auxiliary catalog covers the period from 1400 to 23 August 2023. Figure 10 shows the distribution map of the simulated 10 000 earthquake catalogs by SimplETAS.

Table 4 The SimplETAS parameters used for simulations, with M_min= 3.95, which is similar to Mancini and Marzocchi (2023).

Figure 9 The estimated background seismicity spatial PDF μ (x,y) for the Sichuan–Yunnan region. The black rectangle represents the area where Xichang is located.

From Fig. 10, it can be seen that there is a significant difference between the location of earthquake clusters simulated by SimplETAS and the distribution of the mainshocks and aftershocks in Sect. 3.2.

The potential source models we employed to simulate earthquake catalogs in Sect. 3 comprehensively consider various data, including paleoearthquakes, historical earthquakes, seismogenic structures, and stress–strain fields. These data help constrain the locations of earthquakes, especially those of high magnitude. However, it is important to note that the ETAS model is an empirical statistical model, relying on earthquake catalogs as its fundamental data. This distinction makes it challenging to draw direct comparisons between the two models. To address this limitation, it is essential for future research to explore the incorporation of more physics-based models to establish comparative bridges. However, this endeavor goes beyond the scope of the current study.

Figure 10 The distribution map of the simulated 10 000 earthquake catalogs by SimplETAS. The black rectangle represents the area where Xichang is located.

6 Discussion

Several studies have shown that the estimates of the ETAS parameters are highly susceptible to the assumptions made, such as the magnitude cutoff, time dependency of the background rate, anisotropic aftershock triggering, and aftershock incompleteness (Seif et al., 2017; Zhuang et al., 2017). Bearing this in mind, it becomes apparent that comparing parameter values across different studies using diverse catalogs (with variations in quality, magnitude of completeness, and spatial and temporal windows) is not a straightforward task. Moreover, the inherent statistical correlation among the parameters further complicates the comparison process.

The outcomes of the study conducted by Iacoletti et al. (2022) suggest that the traditional region-wide calibration approach is inadequate for constructing an ETAS model suitable for simulation-based PSHA. Generally, sequence-averaged ETAS models prove to be more acceptable, exhibiting both a higher number of aftershocks and consistent spatial and magnitude–frequency distributions. Nevertheless, numerous regions (such as Xichang) face a challenge due to an insufficient occurrence of active sequences within the required (short and recent) period according to the method's criteria.

Šipčić et al. (2022) conducted a comparison of three alternative models (Poisson, Omori, and ETAS) under two different initial conditions: an “unconditional” case, with initial conditions characterized by average seismicity, and a “conditional case”, incorporating initial conditions of an ongoing active earthquake sequence. As expected, the traditional Poissonian approach for earthquake occurrence modeling tends to provide lower hazard estimates. The inclusion of aftershocks in the Omori model and consideration of all events in the ETAS model significantly enhance hazard estimates, providing more realistic values by not solely accounting for the effect of the largest events, as seen in the case of the Poissonian approach.

In our study, we have examined the classical Poissonian model that considers only mainshocks and the model that combines the Poissonian model for mainshocks and the Omi–R–J model for aftershocks, which is considered an approach for clustered seismicity modeling that is less complicated than ETAS, and the Omi–R–J model is sensitive to the identification of mainshocks.

The significant feature of our study is the simulation of the mainshocks based on the potential source model and the seismicity model of the Seismic ground motion parameters zonation map of China (GB 18306-2015). These models comprehensively consider various data, such as paleoearthquakes, historical earthquakes, seismogenic structures, and stress–strain fields, and provide probability functions for the spatial distribution of earthquakes with different magnitude ranges (Gao, 2015), thereby limiting the location of mainshocks (especially high-magnitude earthquakes). After the determination of the mainshocks, the aftershocks are distributed around the mainshocks. However, the ETAS model is an empirical statistical model, and the fundamental data are only earthquake catalogs. Therefore, the accuracy of the ETAS model depends on having well-characterized catalogs. These findings suggest the need to additionally investigate and improve the models through more sophisticated statistics and physics-based models (Hardebeck et al., 2023).

7 Conclusions

In this study, a probabilistic seismic hazard analysis based on the Monte Carlo method was combined with the Omi–R–J model to systematically study how aftershocks impact seismic hazard analyses in the city of Xichang and the surrounding areas. The results show that in areas with moderate to strong seismic backgrounds, the influence of aftershocks on probabilistic seismic hazard analysis can exceed 50 %. Aftershocks are typically ignored in traditional probabilistic seismic hazard analyses, which underestimate the seismic hazard to some extent and may cause potential risks. Our results suggest that the impact of aftershocks should be properly considered during future probabilistic seismic hazard analyses, especially in areas with moderate to strong seismic activity backgrounds and in areas prone to secondary disasters such as landslides and mudslides.

The model settings adopted for the calculation processes presented in this study can be modified according to the actual situation and specific requirements. The Monte Carlo method is highly adaptable and can take into account different parameters in different models. In future work, we can attempt to adjust the initial magnitude of the mainshock and the aftershock. Additionally, we can adjust the duration of the aftershock and use different mainshock models and aftershock models to study how aftershocks impact probabilistic seismic hazard analysis.

This work provides a scientific basis for governmental departments to minimize disaster losses and formulate corresponding earthquake prevention and disaster mitigation measures. Furthermore, this work plays very important roles in engineering decision-making and judgment, the implementation of catastrophe insurance, and other fields.