MCPSHA: A New Tool for Probabilistic Seismic Hazard Analysis Based on Monte Carlo Simulation

Abstract

The utilization of the Monte Carlo method in conjunction with probabilistic seismic hazard analysis (PSHA) constitutes a compelling avenue for exploration. This approach presents itself as an efficient and adaptable alternative to conventional PSHA, particularly when confronted with intricate factors such as parameter uncertainties and diverse earthquake source models. Leveraging the Monte Carlo method and drawing from the widely adopted Cornell-type seismicity model in engineering seismology and disaster mitigation, as well as a seismicity model capturing temporal, spatial, and magnitude inhomogeneity, we have derived a formula for the probability of earthquake intensity occurrence and the mean rate of intensity occurrence over a specified time period. This effort has culminated in the development of a MATLAB-based program named MCPSHA. To assess the model’s efficacy, we selected Baoji City, Shaanxi Province, China, as our research site. Our investigation delves into the disparity between occurrence probability and extreme probability (a surrogate commonly employed for occurrence probability) in the Baoji region over the next 50 years. The findings reveal that the Western region of Baoji exhibits a heightened hazard level, as depicted in the maps, which illustrate a 10% probability of exceedance within a 50-year timeframe. The probability of earthquake occurrence under various intensities (VI, VII, and VIII) over 50 years follows a declining trend from west to east. Furthermore, the likelihood of seismic intensity exceeding VI, VII, and VIII indicates the lowest exceeding probability in the northeast and the highest in the northwest. Notably, for intensities VI-VII, the difference between occurrence probability and extreme probability approaches twice, gradually diminishing with increasing intensity. This study underscores the MCPSHA model’s efficacy in providing robust technical support for mitigating earthquake risk and enhancing the precision of earthquake insurance premium rate calculations.

1. Introduction

China, as one of the world’s most seismically active countries [1], faces significant earthquake risks, with over half of its cities situated in zones prone to earthquakes of intensity VII or higher [2,3]. Recognizing the pivotal role of seismic risk analysis in disaster prevention and reduction, the fundamental aim is to calculate and analyze ground motion characteristics over time, forming an integral part of engineering seismic fortification [4]. Given the escalating global trend in earthquake disasters, conducting medium- and long-term seismic risk assessments becomes crucial. Simultaneously, these assessments provide a scientific foundation for disaster prevention and reduction planning [5,6,7], as well as informing decisions related to disaster reconstruction and site selection [8,9].

The probabilistic seismic hazard analysis (PSHA) result is generally expressed by ground-motion with a varying probability of exceedance at a site [10]. Initially introduced by Conell in 1968 [4], this methodology has been extensively applied for seismic risk assessment in numerous regions [11,12,13]. However, the Cornell–McGuire method, a prominent approach in PSHA, is often perceived as non-transparent by end-users due to its lack of clarity regarding the relative contributions of different sources and magnitudes to seismic hazards. Many scholars have made notable contributions and enhancements to the program. For instance, Slejko et al. [14] applied a logic tree approach to calculate robust seismic hazard estimates for the northeastern regions of Italy. Their methodology involves various seismotectonic knowledge levels, seismicity rate computation methods, maximum magnitude estimation approaches, and PGA attenuation relations. The resulting hazard maps, processed through a GIS system, highlight elevated hazard areas in central Friuli and around Vittorio Veneto. Additionally, D’Amico and Albarello [15] introduced a method specifically designed for estimating seismic hazards through a probabilistic analysis of macroseismic data. This method, available freely upon request, facilitates the computation of hazard curves, intensity, or PGA values with fixed exceedance probabilities. It leverages locally available seismic effects information, integrating both actual and “virtual” intensity values from neighboring sites while considering uncertainties and reliability factors. The method also incorporates adaptable, region-dependent coefficients suitable for international applications.

In earthquake loss estimation, the extreme probability of intensities is commonly used to replace the probability of occurrence at a certain intensity, which is then combined with the vulnerability matrix to obtain seismic risk results [16,17]. Although the assessment results are expressed as a percentage, the PSHA results give the probability of exceeding a certain ground motion intensity, and the seismic vulnerability analysis is the expected degree of losses within a defined area under earthquakes with various seismic intensities, often expressed in matrix form [18]. In recent years, only a few studies about the occurrence probability of intensity have been carried out [19]. Fu et al. [20] use the Bernoulli repeated test to simulate and calculate the probability of ground motion occurrence. Attempts have been made to improve the seismic intensity integral interval of the Poisson model, and the complex integral formula is simplified using the Taylor formula by Gao [19]. Based on the actual ground motion, Peresan et al. [19] define a reasonable and reliable design earthquake, and the deterministic method is used. However, only by increasing the constraints and simplifying the seismic activity model or mathematical expression can the algorithm be implemented. The intensity exceedance probability is calculated using the extreme value distribution, and the final expression is thoroughly considered by taking into account various factors such as the potential source and its magnitude upper limit, as well as the temporal and spatial distribution unevenness of seismic activity. Using hypothetical conditions to reduce the mathematical expressions of the model may result in large deviations in the calculation results, especially in high intensity, and addressing the problem of the uncertainty of various factors caused by model simplification is also challenging. As a result, understanding how to build the resulting occurrence probability of the intensity model is critical for the study of earthquake loss assessment.

The data required for a seismic hazard assessment using traditional probabilistic seismic hazard assessment methods (PSHA) can also be used for probabilistic hazard analysis using Monte Carlo simulation methods [21,22,23]. Monte Carlo simulation methods stand out due to their exceptional flexibility and transparency, proficient handling of parameter uncertainties, and implementing complex seismic activity models [24]. Bolotin [25] systematically explores random factors in seismic risk assessment for structures using the Monte Carlo method, examining the randomization of crucial parameters in numerical simulations and employing asymptotic distributions of extremes to enhance the reliability of estimates, particularly in the context of rare seismic events with multiple response restrictions. Emmi and Horton [26] employ a CIS-based Monte Carlo simulation to assess the impact of random perturbations in earthquake ground shaking intensity zone boundaries and consider sensitivity to intensity estimation errors, revealing that variations within a 1560 m corridor have minimal effect on property loss estimates, while a five percent error in ground shaking intensities leads to an 11 to 15 percent error in property loss estimates, demonstrating the use of GIS to simulate error propagation and identify strategic error reductions in risk assessments. Numerous established software packages exist [24,27,28,29]. These approaches rely on Monte Carlo simulations for hazard calculation and have found extensive application in various seismic risk assessments [22,30,31]. Nevertheless, the package software may not be suitable for the Chinese area. CPSHA (China Probability Sesmic Hazard Assessment) is a modified PSHA method that has advantages in characterizing the non-uniform spatial distribution of intraplate earthquakes in China [32,33]. The effect of the Monte Carlo method with CPSHA is unknown on the resulting probabilistic results in China.

To fill this gap, we present an exploration of the framework for Chinese probabilistic seismic hazard analysis utilizing the Monte Carlo approach (MCPSHA) implemented in MATLAB. This approach allows for the determination of earthquake intensity, occurrence, and exceedance probabilities. The three integral components of MCPSHA, namely the seismic source model, earthquake recurrence model, and ground motion attenuation model, are comprehensively addressed. The Baoji area is selected as the research site. To assess the compatibility of our proposed method with conventional approaches, seismic hazard maps (with a 10 percent probability of exceedance in 50 years) are computed on a 0.05-degree grid using both our program and a classic seismic hazard assessment program (Cornell–McGuire approach). Furthermore, the occurrence and extreme probability of VI, VII, and VIII seismic intensities over 50 years are calculated and analyzed. The development of the framework and the results derived from MCPSHA serve as a benchmark for comprehensive seismic hazard and seismic risk analysis.

2. Method

2.1. Monte Carlo Simulation

Monte Carlo methods, considered numerical techniques, can be broadly defined as statistical simulation methods, where statistical simulation refers to any physical model employing sequences of random numbers for simulation purposes. An integral aspect of a Monte Carlo simulation involves modeling the physical process using one or more probability density functions (pdf). By employing random sampling, a substantial number of samples that satisfy the pdf are generated. Consequently, by compiling the statistics of the occurrence time “n” in “N” simulation runs, the probability of event occurrence can be determined using Equation (1) [34].

$$\mathrm{P}=\frac{n}{N}$$

(1)

where N is the number of experiments and n is the occurrence frequency of the event in N times of experiments.

Theoretically, within a given physical or mathematical model, the relationship between a core parameter(s) and its probability distribution function (F(s)) is established, where F(s) is a monotonically increasing function ranging from 0 to 1, as expressed in Equation (2).

F(s) = r  for  0 ≤ r ≤ 1

(2)

where F(s) is continuous and monotonic, it assumes that r follows a uniform distribution in [0, 1].

$$s={\mathrm{F}}^{-1}\left(r\right)$$

(3)

If a random value r is generated according to a uniform distribution, we can obtain the variable s according to Equation (3). This method is a direct sampling technique (DST), and we generate a random number ($r$) using the Mersenne Twister algorithm [35,36].

2.2. Strategy for Generating a Stochastic Earthquake Catalogue

The seismic zoning map of China establishes the corresponding earthquake occurrence probability model, a tri-class seismic source model (Figure 1), and provides the basic parameters based on geophysical data obtained from each region. China’s earthquake activity model is based on a three-class seismic source model, which includes the seismic belt, the background seismicity zone, and the tectonic feature zone [37,38].

Based on the basic assumptions and seismic activity parameters of the Chinese ground motion parameter zoning map, the following steps can be used to synthesize the seismic event set:

(1)

Youngs and Coppersmith [39] propose the exponential magnitude probability distribution approach to deal with the possible maximum ($m_{uz}$) and minimum ($m_0$) magnitudes, as expressed in Equation (4). Based on the basic assumptions of Chinese probabilistic seismic hazard analysis, in the seismic belt, the magnitude distribution of earthquakes satisfies the doubly-truncated Gutenberg-Richter recurrence relation Equation (4) expressed below.

$$F\left(m\right)=\left\{\begin{array}{cc}\frac{1-{\mathrm{e}}^{-\beta \left(m-m_0\right)}}{1-{\mathrm{e}}^{-\beta \left(m_{uz}-m_0\right)}}\hfill & , m_0\le m\le m_{uz}\\ 0& , other conditions\end{array}\right.$$

(4)

The value of $F\left(m\right)$ is monotonous in [0, 1]. Let $F\left(m\right)=r$ (r follows uniform distribution in [0, 1]). We can obtain Equation (5) through the inversion transformation of Equation (4).

$$m=m_0-\ln \left\{1-r\times \left[1-e^{-\beta \left({\mathrm{m}}_{uz}-{\mathrm{m}}_0\right)}\right]\right\}/\beta$$

(5)

Therefore, given a random number r from [0, 1], the magnitude of the earthquake $m$ can be determined (Figure 2) according to Equation (5).

(2)

The seismic activity in the seismic province follows the Poisson distribution [37]. In [0, t], the probability of a total of $n$ earthquakes will be observed, as shown in Equation (6).

$$P\left(n=k\right)=\frac{{\left(v_0\times t\right)}^k}{k!}·e^{-v_0\times t} ,\left(k=1,2,3\dots \dots \right)$$

(6)

where k is the number of earthquakes; $v_0$ describes the mean activity rate of earthquake occurrence with magnitude equal or greater than $m_0$.

We can acquire the cumulative distribution function of Equation (6) as shown in Equation (7):

$$P\left(1\le n\le k\right)={{\displaystyle \sum }}_{n=1}^k\frac{{\left(v_0\times t\right)}^k}{k!}.e^{-v_{0}t}$$

(7)

The value of $P\left(1\le n\le k\right)$ is continuity and monotonicity, and thus we can get the value of k through the DST.

If $t_i\left(i=1,2,\dots ,k\right)$ is set to the original time of an earthquake, the time interval of adjacent earthquakes (${\tau }_i$=$t_i-t_{i-1}$) follows the exponential distribution [40], as shown in Equation (8):

$$F{\tau }_i=\left\{\begin{array}{c}0 , {\tau }_i<0\\ 1-e^{-v_0\times {\tau }_i}, {\tau }_i\ge 0\end{array}\right.$$

(8)

According to the DST, the time interval ${\tau }_i$ of the earthquake event $i$ can be obtained. Let $t_i=0$, then, the sequence of earthquake events can be obtained.

(3)

Determination of the location of the epicenter. The non-uniform spatial distribution of earthquake activity within the seismic belt is manifested by the activity of earthquakes at the potential source (which consists of the background seismicity zone and the tectonic features zone). In practice, the magnitude (m) is divided into $N_m$ intervals, and $m_j$ ($m_j-\frac{1}{2}△m\le m_j\le m_j+\frac{1}{2}△m$) represent the mid-value of the magnitude interval $j$. The spatial distribution function ($f_{hl,mj}$) implies the occurrence probability of an earthquake generated by an individual potential source l with a magnitude interval mj.

According to the value of the $f_{hl,mj}$, we can obtain the cumulative distribution function, as shown in Equation (9)

$$F\left(l,m_j\right)={{\displaystyle \sum }}_{l^{\prime }=1}^l{f}_{h{l}^{\prime },mj},\left(l=1,2,\dots ,N_{hl}\right)$$

(9)

The value of $F\left(l,m_j\right)$ is monotonous between [0, 1].

$$F\left(l,m_j\right)=r$$

(10)

r follows a uniform distribution in [0, 1]. With the same method as the DST, given the random value r, the value of $F\left(l,m_j\right)$, $l$ can be determined according to Equations (9) and (10). For example, Figure 2b shows the spatial cumulative distribution function of the Liupanshan-Qilianshan seismic belt. The Liupanshan-Qilianshan seismic belt refers to a seismic zone in China that spans the Liupanshan Mountains and Qilianshan Mountains. These mountain ranges are located in the north-central part of the country. The seismic belt is characterized by tectonic activity, including the interaction of the Eurasian Plate and the Indian Plate. This geological region is known for its seismicity and has experienced earthquakes over the years. The graph’s horizontal axis provides a unique number for every background seismicity zone and tectonic feature zone, and the vertical axis shows the value of $F\left(l,m_j\right)$. According to Equations (9) and (10), given the random value r, the value of $F\left(l,m_j\right)$, which represents the number of $l$ can be determined.

Within the background seismicity zone and tectonic feature zone, the earthquake activity is supposed to be uniform. In order to generate the epicenter, we build up the minimal bounding rectangle of every background seismicity zone and tectonic feature zone and subsequently generate a two-dimensional random number ($r_x$, $r_y$) that follows a uniform distribution between 0 and 1. According to Equations (11) and (12), we can obtain the point $(x_i,y_i)$ within the minimal bounding rectangle. If the point is located in the background seismicity zone or tectonic feature zone, it is considered the epicenter (Figure 3).

$$x_i=\left(x_{max}-x_{min}\right)r_x+x_{min}$$

(11)

$$y_i=\left(y_{max}-y_{min}\right)r_y+y_{min}$$

(12)

where $x_{max}$ and $x_{min}$ are respectively the maximum and minimum longitudes of the minimal bounding rectangle, and $y_{max}$ and $y_{min}$ are the maximum and minimum latitudes of the minimal bounding rectangle, respectively.

(4)

The seismic azimuth angle (θ) is established by sampling the probability function associated with the long-axis azimuth of potential source region attenuation. By taking into account the comprehensive influence of seismic belts on the site, we iteratively execute these procedures until the specified number of earthquakes is met within the seismic belt, completing a sampling cycle. Consequently, seismic event sets can be generated through extensive sampling simulations.

(5)

The conversion of latitude and longitude coordinates to Cartesian coordinates is performed for each earthquake within the event set. Using the epicenter as the coordinate origin, a Cartesian right-handed rectangular coordinate system is established, with the x-axis pointing to the true east and the y-axis pointing to the true north. This coordinate system undergoes a counterclockwise rotation by a specific angle (θ). The coordinates of the site in the new Cartesian coordinate system, denoted as ($x_c,y_c$), can be determined through simultaneous equations, as illustrated in Equation (13).

$$\left\{\begin{array}{c}\frac{x_c^2}{R_a^2}+\frac{y_c^2}{R_b^2}=1\\ I\left(m,R_a\right)=I\left(m,R_b\right)\\ \begin{array}{c}I\left(m,R_a\right)=a1+a2\times m+a3\times m^2+a4\times log\left[R_a+a5\times \exp \left(a6\times m\right)\right]+a7\times R_a+{\epsilon }_a\\ I\left(m,R_b\right)=b1+b2\times m+b3\times m^2+b4\times log\left[R_b+b5\times \exp \left(b6\times m\right)\right]+b7\times R_b+{\epsilon }_b\end{array}\end{array}\right.$$

(13)

where a1–a7, b1–b7 are the regression constants, and $I\left(m,R_a\right)$ and $I\left(m,R_b\right)$ are the intensity along the major axis and minor axis, respectively. $R_a$ and $R_b$ are the semi-major axis and semi-minor axis of an ellipse, respectively.

Utilizing the iterative method [41] to solve the transcendental equation enables the determination of intensity values. Furthermore, the normalized residual (ε) is typically assumed to follow a normal distribution with a mean of zero and a standard deviation (σ). By incorporating ε into the intensity, the final intensity can be obtained. Considering the potential influence of all seismic zones on the site, the process is continued until the specified number of earthquakes is reached within the seismic zone, completing a sampling cycle. Subsequently, through large-scale sampling simulations, seismic event sets can be generated.

2.3. Estimation of Occurrence Probability of Earthquake Intensity

Based on the Monte Carlo method, for a site $s$, our goal is to predict the probability of the seismic intensity exceeding specified seismic intensity in T years. $Let N_r$ represents the number of earthquake catalog groups simulated in the forthcoming T years. The number of simulated earthquake events in the g-th $\left[g=\overline{1,N_r}\right]$ group is denoted as $k\left(g\right)$ and the ground motion generated by the i-th $\left[i=\overline{1,k\left(g\right)}\right]$ earthquake event on the site is $A_{gi}$.

The total number affected by earthquake magnitude $M\ge m_0$ during T years is:

$$K_{sum}={{\displaystyle \sum }}_{g=1}^{g=N_r}k\left(g\right)$$

(14)

The number of ground motions ($A$) which is greater than $a_0$ and less than $a_{0+1}$ is shown in Equation (15).

$$N\left(a_{0+1}\ge A\ge a_0\right)={{\displaystyle \sum }}_{g=1}^{g=N_r}count\left[a_{0+1}\ge A_{gi}\ge a_0|i=\overline{1,k\left(g\right)}\right]$$

(15)

Then, in T years and within the influence range of the site s, an earthquake with m $\ge m_0$ occurs, and the probability that the ground motion value is greater than $a_0$ and less than $a_{0+1}$ is shown in Equation (16).

$$p_1\left(a_{0+1}\ge A\ge a_0\right)=\frac{N\left(a_{0+1}\ge A\ge a_0\right)}{K_{sum}}$$

(16)

In T years, the mean rate of earthquake intensity occurrence at site s is shown in Equation (17).

$$V_s\left(a_{0+1}\ge A\ge a_0\right)=\frac{N\left(a_{0+1}\ge A\ge a_0\right)}{Nr}$$

(17)

In T years, the average number of earthquakes ($Nmean$) with a magnitude of $M\ge m_0$ is shown in Equation (18).

$$Nmean=\frac{K_{sum}}{Nr}$$

(18)

Therefore, the occurrence probability of earthquake intensity ($p\left(a_{0+1}\ge A_{con}\ge a_0\right)$) in T years is shown in Equation (19):

$$\mathrm{p}\left(a_{0+1}\ge A_{con}\ge a_0\right)=1-{\left[1-p_1\left(a_{0+1}\ge A\ge a_0\right)\right]}^{Nmean}$$

(19)

2.4. Estimation of Exceedance Probability of Earthquake Intensity

Based on the Monte Carlo method, in the g-th random independent repeat test, the maximum value of the ground motion generated by $k\left(g\right)$ events on the site s is:

$$A_{max}\left(g\right)=max\left\{A_{gi}|i=\overline{1,k\left(g\right)}\right\}$$

(20)

The number of ground motions that is greater than the set value $a_0$ and less than $A_{max}\left(g\right)$ is written as:

$$N\left(A_{max}\ge a_0\right)=count\left[A_{max}\left(g\right)\ge a_0|g=\overline{1,N_r}\right]$$

(21)

$$P\left(A_{max}\ge a_0\right)=\frac{N\left(A_{max}\ge a_0\right)}{Nr}$$

(22)

Correspondingly, the extreme probability of the intensity of the specified seismic intensity ($a_0$) in T years is:

$$P\left(A_{max}=a_0\right)=P\left(A_{max}\ge a_0\right)-P\left(A_{max}\ge a_{0-1}\right)$$

(23)

MCPSHA represents a relatively straightforward approach for seismic hazard assessment, depicted schematically in Figure 4. The program comprises three modules: (1) simulating earthquake catalog generation based on seismicity parameters; (2) calculating ground motion at the designated site; and (3) determining ground motion probability at the site.

3. Geological Setting and Seismic Activity Parameters

Baoji City, situated in Shaanxi Province, is encompassed by the Loess Plateau, Qinling Mountain, and Weihe Basin (Figure 5). The Weihe River traverses the region from west to east. The Baoji area is influenced by three tectonic units, namely the QinLing fold belt, the Fen-Wei fault-depression zone, and the southern margin of the Ordos block [42]. Extensive seismic events have been recorded in this locale, with data sourced from the Catalogue of Shaanxi Earthquakes, 1189BC-2001AD (2015) (http://news.ceic.ac.cn/index.html?time=1524648589, accessed on 13 May 2019). Figure 5 shows the historical earthquake distribution around the Baoji area, spanning from 780 BC to 2015 AD. Notably, since 1973, the region has experienced numerous moderate and small earthquakes characterized by periodicity, migration, and recurrence.

A total of 7 seismic belts may affect the Baoji area. In each background seismicity zone, tectonic feature zones are recognized and classified (Figure 6). The seismic activity parameters and the upper limit of magnitude for each seismic belt are measured based on the results of the “seismic microzonation of 4 cities in the Guan-Zhong” project (Figure 6). The background seismicity zone and tectonic feature zone parameters are obtained based on the results of the 5th National Seismic Zoning Map of China [43]. Seismic motion attenuation relations are based on the results of the “seismic microzonation of four cities in the Guan-Zhong” project (Figure 7) [44,45].

4. Results and Analyses

4.1. Earthquake Catalog Simulation

To assess whether the generated seismic catalog adheres to the predefined seismic activity parameters and satisfies the assumed Poisson distribution hypothesis, we individually tabulate the annual occurrences of earthquakes with a magnitude of 4 or higher for each catalog group. A total of 50 samples (T = 50 years) were obtained. In the case of seismic events occurring $i$ times $\left(i=0,1,2,3,\dots \right)$, the probability density $p\left(i\right)$ of earthquake occurrences per year is calculated.

$$p\left(i\right)=\frac{n_i}{50}$$

(24)

Figure 8 presents the Poisson distribution fitting of annual earthquake occurrences with a magnitude of 4 or higher, obtained through nonlinear fitting methods, based on four sets of simulated seismic catalogs from the Fenwei seismic zone (T = 50 years).

We compute the ${\mathit{v}}_{\mathit{h}}$ and b values for each subgroup of earthquake catalogs following the aforementioned procedure and then determine the mean values for each seismic belt. Figure 9 illustrates the parameters for the seven seismic belts obtained through simulation tests, while

Table 1. Simulation catalogs (1,000,000 groups) of $v_h$ and b value test parameter.

Seismic Province	b	${\mathit{v}}_{\mathit{h}}$
	Standard Deviation	Standard Deviation
Middle-reach of Yangtze River	0.13	0.43
North China plain	0.11	0.39
Ordos	0.12	0.26
Fenwei	0.11	0.40
Yinchuan-Hetao	0.14	0.37
Longmenshan	0.15	0.30
Liupan-qilian Mountain	0.16	0.13

provides details on the error parameters. The findings indicate that a sufficient number of repetitions guarantees the accuracy of the parameters.

4.2. Estimation of the Probability of an Earthquake Intensity

Using the method described in this paper, 1,000,000 groups of simulated earthquake catalogs over the next 50 years are generated. Figure 10 shows the distribution of earthquake epicenters (5.0 ≤ M ≤ 8.5) in 1 group of simulated earthquake catalogues in the next 50 years.

Using a grid size of 0.05° × 0.05°, we compute the seismic intensity for the Baoji area at a specific hazard level. Macro-seismic intensity uncertainty is randomly drawn from a distribution with a mean of zero and a standard deviation of 0.574, truncated at two standard deviations. Figure 11 illustrates the mean intensity incidence rates for VI, VII, and VIII sites in the region over 50 years based on Equation (18). The results reveal the lowest incidence rate in the northeast and the highest in the northwest. One contributing factor may be the presence of more tectonic feature zones with a m_uz greater than 7.5 in the northwest, as compared to the northeast and southeast of Baoji city.

Figure 12 displays the variation in intensity with the probability of occurrence and exceedance over 50 years in the region, calculated using our program. Notably, the difference between the two probabilities is more pronounced at intensity levels VI~VII, whereas at intensity level VIII, the difference is smaller.

Figure 13 illustrates a comparison of the mean curves for the two probabilities (VI, VII, and VIII) in the region. The results reveal that as the intensity value increases, the disparity between the two gradually diminishes. Baoji is categorized as a moderate earthquake area, where occurrences of intensity above VIII are infrequent, and the occurrence probability of the intensity aligns closely with the extreme probability of the intensity value. However, in seismic loss estimation, specific regions necessitate individualized analysis owing to their diverse tectonic environments to determine their fundamental consistency. Despite minor distinctions, these differences may still introduce deviations in the evaluation results. Therefore, in the context of general seismic damage loss estimation, the exceedance probability of intensity cannot substitute for the occurrence probability of intensity.

5. Discussion

To verify the accuracy of the results generated by our program, 300 points are randomly chosen in the Baoji area, and the hazard (return period T = 476 years) for each site is computed using both the traditional CPSHA method and the MCPSHA program (shown in Figure 14). Although the results from both methods exhibit similarities, some differences persist (Figure 14). The correlation between the two intensity results is high, at 0.98 (Figure 15a). To further assess the disparity between the CPSHA and MCPSHA results, we calculate the difference by subtracting the two outcomes (Figure 15b). The average difference is 0.16, with a minimum of 0.08, a maximum of 0.2, and a standard deviation of 0.03. The observed difference may stem from the excessive discretization employed by the CPSHA program, specifically in variables such as earthquake magnitude and epicenter distance, intended for simplified integration. In MCPSHA, the utilization of the Monte Carlo method based on Bayer’s large number theorem ensures that, with an ample number of random simulations, the frequency of events precisely reflects the probability of their occurrence. Additionally, a comprehensive series of repeated experiments assures the occurrence of all possible conditions. Consequently, as the number of simulations increases, even events with low probability are taken into account.

The most appropriate and scientifically sound approach for validating the results of the PSHA method is to rely on historical data specific to the region. The seismic activity history spanning approximately 3000 years in the studied area serves as a valuable validation tool. We have presented the expected occurrence probability for a mean return period of 476 years, a period that repeats more than six times within the historical information of the region. However, due to a lack of data, we were unable to locate several or a few sites with recorded macroseismic intensities. In future research, we intend to compare the modeling of observed intensity occurrences with those estimated by MCPSHA. Such a comparison is sufficient to either substantiate or negate the robustness of the proposed methodology.

The current seismic hazard analysis method is generally utilized for engineering seismic fortification; thus, only the extreme intensity of the site is taken into consideration. This method is not applicable to earthquake damage estimation [46,47]. The seismic loss of a site is not only caused by an earthquake, which corresponds to its extreme intensity. Although the original extreme value (only considering one maximum case) correctly handles the randomness of future earthquake occurrences, if the duration is long and the site is located in a region with frequently occurring strong earthquakes, the loss caused by multiple earthquakes with the same intensity and that associated with one earthquake may obviously differ [48]. Therefore, when calculating the expected earthquake-induced loss, both the occurrence probability and the occurrence frequency of the seismic intensity should be considered. For example, the site may suffer multiple VI-intensity ground motions within a given number of years, and the loss produced in this situation is twice that generated by an earthquake with an intensity of VI. The estimation method based on extreme intensity ignores the above effects. Therefore, the seismic hazard analysis method of the engineering site cannot be directly used.

Both the probability of occurrence and exceedance are important indicators of seismic hazards for a given site, but the meaning of the two indicators is distinctly different. Furthermore, the exceedance probability itself contains the meaning of extreme intensity. For any site, it is assumed that the frequency of strong earthquakes will be extremely high in the next 50 years, and the possibility of earthquakes with a magnitude of around 6 is very high (even more earthquakes with a magnitude less than 6). Under this assumption, it is almost inevitable that the site will be affected by intensities VI and VII within a given number of years. The occurrence probability of intensity VI is high, but the probability of intensity VI as an extreme value is relatively small (Figure 12). The method of subtraction adjacent exceedance probability within a given period of time can be considered the exceedance probability of the extreme intensity at the site within a specified number of years. Numerically, it is lower than the occurrence probability of the intensity (Figure 13). From a statistical point of view, if the extreme intensity is a fixed value, the probability of the site being affected is zero. There are two possibilities in this situation: The first is that the intensity values that affect the site are all less than the fixed value within a given period of time; the second is that the probability of an extreme intensity value in a given period of time is 100% equal to 1 (i.e., it is inevitable). In the second case, although the subtraction adjacent exceedance probability is zero, the probability that the intensity will affect the site may be large. Therefore, the exceedance probability of intensity cannot be used to replace the occurrence probability of intensity.

The Monte Carlo method poses computational challenges, particularly for intricate seismic models, leading to high computational costs due to the necessity for an extensive number of random simulations. Moreover, the accuracy of the method hinges on both the quantity and quality of random sampling; insufficient or inappropriate sampling could introduce biases to the results. In the seismic zone earthquake simulation of this study, the Poisson model employed for earthquake occurrence times often deviates from real seismic activity, typically following a clustered model. Addressing the issue of mutual triggering between earthquakes (such as Poisson + clustering), future research can refine the sampling function and incorporate a seismicity model (e.g., ETAS) to yield more precise results, especially in regions with limited historical earthquake records and complex geological structures. Meanwhile, Bolotin [25] systematically investigates random factors in seismic risk assessment for structures, discussing the randomization of key parameters in numerical simulations and applying asymptotic distributions of extremes for more reliable estimates in rare seismic events with multiple response restrictions. In this study, the uniform setting of key parameters, such as seismic intensity attenuation, to a single value using the Monte Carlo method in tandem with spatiotemporal heterogeneous seismic hazard analysis in China may impact the accuracy of the seismic hazard analysis results.

6. Conclusions

In this paper, we derived the formula for the occurrence probability of earthquake intensity and the mean rate of earthquake intensity occurrence within a specified time period, and we have developed a program named MCPSHA. The three fundamental components of MCPSHA are carefully considered, including the seismic source model, magnitude recurrence model, and ground motion attenuation model. This approach offers greater flexibility and transparency, making it easily adaptable to modifications in the seismicity model. For instance, the introduction of complex seismic activity models requires minimal reprogramming—only a new probability density function (pdf) function needs to be provided. Additionally, the modeling of uncertainty in input parameters is highly flexible, accomplished solely through the use of a standard deviation rather than relying on the discrete branches of a logic tree.

Over the next 50 years, we will investigate the disparity between occurrence probability and extreme probability (i.e., the difference in exceeding probability, often used as a substitute for occurrence probability) in the Baoji area. The findings reveal that, for intensities within VI~VII, the difference is nearly twofold, gradually diminishing with increasing intensity. This study underscores that exceeding probability and occurrence probability for specific intensities cannot be interchangeable. Accurate seismic intensity occurrence probability stands as a crucial foundation for seismic risk analysis and loss estimation.