the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The influence of aftershocks on seismic hazard analysis: a case study from Xichang and the surrounding areas
Qing Wu
Guijuan Lai
Jian Wu
Jinmeng Bi
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 183062015), 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.
 Article
(12585 KB)  Fulltext XML
 BibTeX
 EndNote
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 damageinduced 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 epidemictype 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 (UCERF3ETAS) to evaluate the effects of declustering and Poisson assumptions on seismic hazard estimates. Wang et al. (2022) compared the ETASsimulated hazard with approximations based on the declustered Poisson approach (DP), the nondeclustered Poisson approach (NDP), and the recently proposed sequencebased PSHA (Iervolino et al., 2014).
Based on the Seismic ground motion parameters zonation map of China (GB 183062015), 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.
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 strikeslip 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 rhombusshaped 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.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 shortterm 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
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
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 “statespace” model developed by Omi et al. (2013) was used to estimate the timevarying 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}},\phantom{\rule{0.25em}{0ex}}\mathrm{\dots},\phantom{\rule{0.25em}{0ex}}{\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(M\left(\right)open="">\mathit{\mu}\left(t\right),\mathit{\sigma})$. The logarithmic likelihood function related to parameters $p,c,$ and k is
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 selfsimilarity 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
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 swarmtype 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
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 earlyseismicmonitoring 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.
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 183062015). 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 threelevel 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 moderatemagnitude earthquakes under different tectonic conditions; and the upper layer consists of the structural potential sources, which are used to reflect the smallscale 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.
The seismic zoning map of China (GB 183062015) 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.
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).

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 ν_{4} 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 ν_{4} as parameters, where L is the number of earthquakes in the seismic zone for the length of time T to be simulated.

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_{0} 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)}& \mathrm{log}N=abM,\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}^{abM}.\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.

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.

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 1year seismic sequence in this paper. When the time length is set to 10 years, the sequence is called a 10year 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 1year 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 50year 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 183062015), 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 183062015). Each synthetic sequence represents a 1year 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}}^{abM}$. 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
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: ${\mathrm{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 5a shows a schematic diagram of the spatial distribution of 1year 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 1year samplings.
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 183062015), the GMPEs of peak ground acceleration (PGA) suitable for the Xichang area are as follows (Xiao, 2011):
When M < 6.5,
and when M≥6.5,
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.
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 50year 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 50year exceedance probability of 10 % in Xichang and its surrounding areas calculated without and with aftershocks as well as the aftershock impact rate distribution map.
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
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_{7} 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 shows the PGA (gals) contour map of the 50year 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_{7} 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.
The ETAS model is a widely used statistical method for capturing shortterm 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:
in which A is the productivity; α is the coefficient of the exponential magnitudedependent 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.
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 physicsbased models to establish comparative bridges. However, this endeavor goes beyond the scope of the current study.
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 regionwide calibration approach is inadequate for constructing an ETAS model suitable for simulationbased PSHA. Generally, sequenceaveraged 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 183062015). 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 highmagnitude 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 wellcharacterized catalogs. These findings suggest the need to additionally investigate and improve the models through more sophisticated statistics and physicsbased models (Hardebeck et al., 2023).
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 decisionmaking and judgment, the implementation of catastrophe insurance, and other fields.
All raw data can be provided by the corresponding authors upon request.
All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by QW and GL. The first draft of the manuscript was written by QW, and all authors commented on draft versions of the manuscript. All the authors have read and approved the final paper.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank the compiling committee of the Seismic ground motion parameters zonation map of China for providing the seismic source zone model, seismicity model, and GMPE model data. Thanks are extended to Zongchao Li for his help in drawing the map showing the seismic event distribution and tectonic background in Xichang and its surrounding areas. Thanks are also extended to Jiawei Li for his assistance with the ETAS model. The authors thank Matteo Taroni and the two anonymous reviewers for their insightful comments, which helped improve the quality of this article.
This work was supported by the National Key Research and Development Program (grant no. 2022YFC3003502); an independent project initiated by the Institute of Geophysics, China Earthquake Administration (grant no. JY2022Z41); and the National Key Research and Development Program of the Xinjiang Uygur Autonomous Region (grant no. 2020B030064).
This paper was edited by Oded Katz and reviewed by Matteo Taroni and two anonymous referees.
Bi, J. M. and Jiang, C. S.: Distribution characteristics of earthquake sequence parameters in North China, Chinese J. Geophys., 62, 4300–4312, 2019 (in Chinese).
Bi, J. M., Jiang, C. S., Lai, G. J., and Song, C.: Effectiveness evaluation and constrains of early aftershock probability forecasting for strong earthquakes in continental China, Chinese J. Geophys., 65, 2532–2545, 2022 (in Chinese).
Boyd, O. S.: Including foreshocks and aftershocks in timeindependent probabilistic seismichazard analyses, B. Seismol. Soc. Am., 102, 909–917, 2012.
Canales, M. R. and Baan, M. V. D.: Are aftershock sequences pertinent to longterm seismic hazard assessments? Insights from the temporal ETAS model, J. Geophys. Res.Sol. Ea., 125, e2019JB019095, https://doi.org/10.1029/2019JB019095, 2020.
Cornell, C. A.: Engineering seismic risk analysis, B. Seismol. Soc. Am., 58, 1583–1606, 1968.
Davoudi, N., Tavakoli, H. R., Zare, M., and Jalilian, A.: Aftershock probabilistic seismic hazard analysis for Bushehr province in Iran using ETAS model, Nat. Hazards, 100, 1159–1170, 2020.
Du, P. S.: Recurrence interval on macroearthquake along the active fault of Zemuhe, Earthquake Research in Sichuan, 12, 102–118, 2000 (in Chinese).
Enescu, B., Enescu, D., and Ito, K.: Values of b and p: their variations and relation to physical processes for earthquakes in Japan and Romania, Rom. J. Phys., 56, 590–608, 2011.
Felzer, K. and Brodsky, E.: Decay of aftershock density with distance indicates triggering by dynamic stress, Nature, 441, 735–738, 2006.
Field, E. H., Milner, K. R., and Luco, N.: The seismic hazard implications of declustering and Poisson assumptions inferred from a fully timedependent model, B. Seismol. Soc. Am., 112, 527–537, 2022.
Gallovič, F. and Brokešová, J.: Probabilistic aftershock hazard assessment I: Numerical testing of methodological features, J. Seismol., 12, 53–64, 2008.
Gao, M. T. (Ed.).: Publicity and implementation textbook of the Seismic Ground Motion Parameters Zonation Map of China (GB183062015), V1.0, Standards Press of China, ISBN 9787506678889, 2015 (in Chinese).
Gerstenberger, M. C., Wiemer, S., Jones, L. M., and Reasenberg, P. A.: Realtime forecasts of tomorrow's earthquakes in California, Nature, 435, 328–331, https://doi.org/10.1038/nature03622, 2005.
Gutenberg, B. and Richter, C. F.: Frequency of earthquakes in California, B. Seismol. Soc. Am., 34, 185–188, https://doi.org/10.1038/156371a0, 1944.
Hardebeck, J. L., Llenos, A. L., Michael, A. J., Page, M. T., Schneider, M., and van der Elst, N. J.: Aftershock Forecasting, Annu. Rev. Earth Planet. Sc., 52, 2.1–2.24, https://doi.org/10.1146/annurevearth040522102129, 2023.
He, H. L. and Ikeda, Y.: Faulting on the Anninghe fault zone, southwest China in late quaternary and its movement model, Acta Seismologica Sinica, 29, 537–548, 2007 (in Chinese).
Iacoletti, S., Cremen, G., and Galasso, C.: Validation of the EpidemicType Aftershock Sequence (ETAS) Models for SimulationBased Seismic Hazard Assessments, Seismol. Res. Lett., 93, 1601–1618, 2022.
Iervolino, I., Giorgio, M., and Polidoro, B.: Sequencebased probabilistic seismic hazard analysis, short note, B. Seismol. Soc. Am., 104, 1006–1012, 2014.
Jiang, H. K., Zheng, J. C., Wu, Q., Qu, Y. J., and Li, Y. L.: Earlier statistical features of ETAS model parameters and their seismological meanings, Chinese J. Geophys., 50, 1778–1786, 2007 (in Chinese).
Japan Meteorological Agency (JMA): The IwateMiyagi Nairiku Earthquake in 2008, Reports coordnated Commuication on Earthquake Prediction, 81, 101–131, https://cais.gsi.go.jp/YOCHIREN/report/kaihou81/03_04.pdf (last access: 25 March 2024), 2009.
Li, P.: XianshuiheXiaojiang fault zone, Beijing: Seismological Press, ISBN 7502809872, 1993 (in Chinese).
Li, P. and Wang, L. M.: Discussion on the basic characteristics of seismogeology in Western SichuanYunnan, Beijing: Seismological Press, 1318030, 1977 (in Chinese).
Lu, P., Yuan, Y. F., Yuan, H. K., Liu, X. L., Yu, X. H., and Duan, Y. S.: Probabilistic estimate of strong earthquake risk in the AnningheZemuhe tectonic zone, Earthquake, 32, 62–72, 2012 (in Chinese).
Lv, X. J., Gao, M. T., Gao, Z. W., and Mi, S. T.: Comparison of the spatial distribution of ground motion between mainshocks and strong aftershocks, Acta Seismologica Sinica, 29, 295–301, 2007 (in Chinese).
Marzocchi, W. and Taroni, M.: Some thoughts on declustering in probabilistic seismic hazard analysis, B. Seismol. Soc. Am., 104, 1838–1845, 2014.
Mancini, S. and Marzocchi, W.: SimplETAS: A Benchmark Earthquake Forecasting Model Suitable for Operational Purposes and Seismic Hazard Analysis, Seismol. Res. Lett., 95, 38–49, https://doi.org/10.1785/0220230199, 2023.
Ogata, Y.: Statistical models for earthquake occurrences and residual analysis for point processes, J. Am. Stat. Assoc., 83, 9–27, 1988.
Ogata, Y.: Statistical model for standard seismicity and detection of anomalies by residual analysis, Tectonophysics, 169, 159–174, 1989.
Ogata, Y.: Detection of precursory relative quiescence before great earthquakes through a statistical model, J. Geophys. Res.Sol. Ea., 97, 19845–19871, 1992.
Ogata, Y. and Katsura, K.: Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues, Geophys. J. Int., 113, 727–738, 1993.
Ogata, Y.: Spacetime pointprocess models for earthquake occurrences, Ann. I. Stat. Math., 50, 379–402, 1998.
Ogata, Y.: Increased probability of large earthquakes near aftershock regions with relative quiescence, J. Geophys. Res., 106, 8729–8744, 2001.
Omi, T., Ogata, Y., Hirata, Y., and Aihara, K.: Forecasting large aftershocks within one day after the mainshock, Sci. Rep., 3, 2218, https://doi.org/10.1038/srep02218, 2013.
Omi, T., Ogata, Y., Shiomi, K., Shiomi, K., Enescu, B., Sawazaki, K., and Aihara, K.: Automatic aftershock forecasting: A test using realtime seismicity data in Japan, B. Seismol. Soc. Am., 106, 2450–2458, https://doi.org/10.1785/0120160100, 2016.
Omi, T., Ogata, Y., and Shiomi, K.: Implementation of a realtime system for automatic aftershock forecasting in Japan, Seismol. Res. Lett., 90, 242–250, 2019.
Omori, F.: On aftershocks of earthquakes, J. Coll. Sci. Imp. U. Tok., 7, 111–200, 1894.
Reasenberg, P. A. and Jones, L. M.: Earthquake hazard after a mainshock in California, Science, 243, 1173–1176, 1989.
Seif, S., Mignan, A., Zechar, J. D., Werner, M. J., and Wiemer, S.: Estimating ETAS: The effects of truncation, missing data, and model assumptions, J. Geophys. Res., 122, 449–469, https://doi.org/10.1002/2016JB012809, 2017.
Shebalin, P. and Narteau, C.: Depth dependent stress revealed by aftershocks, Nat. Commun., 8, 1317, https://doi.org/10.1038/s4146701701446y, 2017.
Shen, W. H. and Yang, F.: Probabilistic aftershock hazard assessment for Jiuzhaigou MS7.0 earthquake in 2017, Acta Seismologica Sinica, 40, 654–663, 2018 (in Chinese).
Shcherbakov, R., Turcotte, D. L., and Rundle, J. B.: Ageneralized Omori 's law for earthquake aftershock decay, Geophys. Res. Lett., 31, L11613, https://doi.org/10.1029/2004GL019808, 2004.
Šipčić, N., Kohrangi, M., Papadopoulos, A. N., Marzocchi, W., and Bazzurro, P.: The Effect of Seismic Sequences in Probabilistic Seismic Hazard Analysis, B. Seismol. Soc. Am., 112, 1694–1709, 2022.
Stallone, A. and Marzocchi, W.: Features of seismic sequences are similar in different crustal tectonic regions, B. Seismol. Soc. Am., 109, 1594–1604, https://doi.org/10.1785/0120180175, 2019.
Tahir, M., Grasso, J.R., and Amorèse, D.: The largest aftershock: How strong, how far away, how delayed?, Geophys. Res. Lett., 39, L04301, https://doi.org/10.1029/2011GL050604, 2012.
Taroni, M. and Akinci, A.: Good practices in PSHA: declustering, bvalue estimation, foreshocks and aftershocks inclusion; a case study in Italy, Geophys. J. Int., 224, 1174–1187, 2021.
Utsu, T.: A statistical study on the occurrence of aftershocks, Geophys. Mag., 30, 521–605, 1961.
Utsu, T., Ogata, Y., and Matsu'ura, R. S.: The centenary of the Omori formula for a decay law of aftershock activity, J. Phys. Earth., 43, 1–33, 1995.
Wang, S., Werner, M. J., and Yu, R.: How well does Poissonian probabilistic seismic hazard assessment (PSHA) approximate the simulated hazard of epidemictype earthquake sequences?, B. Seismol. Soc. Am., 112, 508–526, 2022.
Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, 1994.
Wiemer, S.: Introducing probabilistic aftershock hazard mapping, Geophys. Res. Lett., 27, 3405–3408, 2000.
Wiemer, S. and Katsumata, K.: Spatial variability of seismicity parameters in aftershock zones, J. Geophys. Res.Sol. Ea., 104, 13135–13151, 1999.
Wu, Q. and Gao, M. T.: A preliminary study on the correlativity of seismic hazard between Beijing area and Xiong'an New Area, Seismology and Geology, 40, 935–943, https://doi.org/10.3969/j.issn.02534967.2018.04.015, 2018 (in Chinese).
Wu, Q., Wu, J., and Gao, M. T.: Correlation analysis of earthquake impacts on a nuclear power plant cluster in Fujian province, China, Environ. Res., 187, 109689, https://doi.org/10.1016/j.envres.2020.109689, 2020.
Xiao, L.: Study on the attenuation relationship of horizontal ground motion parameters near the source of rock site, PhD, Beijing: Institute of Geophysics, China Earthquake Administration, 120, 2011 (in Chinese).
Xiong, Z. H.: The study and application of catastrophe model for earthquake insurance, PhD, Beijing: Institute of Geophysics, China Earthquake Administration, 2019 (in Chinese).
Xu, W. J. and Wu, J.: Effect of temporalspatial clustering of aftershocks on the analysis of probabilistic seismic hazard, Chinese J. Geophys., 60, 3110–3118, 2017 (in Chinese).
Xu, X. W., Wen, X. Z., Zheng, R. Z., Ma, W. T., Song, F. M., and Yu, G. H.: The latest tectonic variation pattern and its dynamic source of active blocks in SichuanYunnan region, Sci. China Ser. D, 33, 151–162, 2003a (in Chinese).
Xu, X. W., Cheng, G. L., Yu, G. H., Song, F. M., Xiang, H. F., Zhang, L. F., Ron, H., Wang, Y. L., and Wen, X. Z.: Tectonic and paleomagnetic evidence for the clockwise rotation of the SichuanYunnan rhombic block, Seismol. Geol., 25, 61–70, 2003b (in Chinese).
Yeo, G. L. and Cornell, C. A.: A probabilistic framework for quantification of aftershock groundmotion hazard in California: methodology and parametric study, Earthq. Eng. Struct. Dyn., 38, 45–60, 2009.
Zhao, J. B.: Damage assessment of reinforced concrete frame structures under main and aftershocks, MS thesis, Beijing: Institute of Geophysics, China Earthquake Administration, 2005.
Zhang, C. J., Hou, Y. Y., Hu, B., Xu, H. H., and Wang, D. B.: Analysis on the seismic activities and hazards of M_{7.1} earthquake, 2010 and M_{6.3} earthquake, 2011 in New Zealand, Recent Developments in World Seismology, 4, 44–51, 2011 (in Chinese).
Zhou, B. G., Chen, G. X., Gao, Z. W., Zhou, Q., and Li, J. Y.: The technical highlights in identifying the potential seismic sources for the update of national seismic zoning map of China, Technol. Earthq. Disaster Prev., 8, 113–124, 2013 (in Chinese).
Zhuang, J. C., Ogata, Y., and Wang, T.: Data completeness of the Kumamoto earthquake sequence in the JMA catalog and its influence on the estimation of the ETAS parameters, Earth Planets Space, 69, 1–12, https://doi.org/10.1186/s4062301706146, 2017.
 Abstract
 Introduction
 Aftershock activity models and their parameters
 Probabilistic seismic hazard analysis method considering aftershocks
 Influence of aftershocks on probabilistic seismic hazard analysis
 Comparison with the ETAS model
 Discussion
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Aftershock activity models and their parameters
 Probabilistic seismic hazard analysis method considering aftershocks
 Influence of aftershocks on probabilistic seismic hazard analysis
 Comparison with the ETAS model
 Discussion
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References