# MCPSHA: A New Tool for Probabilistic Seismic Hazard Analysis Based on Monte Carlo Simulation
## 2 Method
### Strategy for Generating a Stochastic Earthquake Catalogue
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)=\sum_{l^{\prime}=1}^{l}f_{hl^{\prime},mj},\left(l=1,2, \ldots,N_{hl}\right) \tag{9}\]  
The value of \(F\left(l,m_{j}\right)\) is monotonous between [0, 1]. \[F\left(l,m_{j}\right)=r \tag{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. Figure 2: Cumulative distribution function (**a**): Magnitude Cumulative Distribution Function. (**b**) the cumulative probability curves of potential sources under different magnitudes (mj). 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(l,m_{j})\). According to Equations (9) and (10), given the random value r, the value of \(F(l,m_{j})\), 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} =(x_{max}-x_{min})r_{x}+x_{min} \tag{11}\] \[y_{i} =(y_{max}-y_{min})r_{y}+y_{min} \tag{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. 1. The seismic azimuth angle (\(\theta\)) 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. 2. 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 (\(\theta\)). 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{aligned} \frac{x_{c}^{2}}{R_{a}^{2}}+\frac{y_{c}^{2}}{R_ {b}^{2}}=1\\ I(m,R_{a})=I(m,R_{b})\\ I(m,R_{a})=a1+a2\times m+a3\times m^{2}+a4\times log[R_{a}+a5\times \exp(a6\times m)]+a7\times R_{a}+\varepsilon_{a}\\ I(m,R_{b})=b1+b2\times m+b3\times m^{2}+b4\times log[R_{b}+b5\times \exp(b6\times m)]+b7\times R_{b}+\varepsilon_{b}\end{aligned} \right.