Seismic Hazard Assessment using Monte Carlo Simulation¶

Prepared by: Cesar P. Malenab Jr.

Executive Summary¶

Statistics of significant earthquakes (Mw ≥ 5.0 and distance ≤ 400 km) are established for a building in Metro Manila based on an earthquake catalog in the past 122 years. Seismic hazard assessment utilizing conventional Probabilistic Seismic Hazard Analysis (PSHA) may be validated using the Monte Carlo Simulation (MCS) method and this notebook presents a framework for performing MCS in Python for seismic hazard assessments. The seismicity of the project site is determined through descriptive statistics of various earthquake parameters which will also serve as input variables for the simulation. This notebook will also highlight common methods, function calls, and visualization tools that are important in statistical studies and random sampling methods.

Introduction¶

Prediction of ground motion intensities of different hazard levels is essential, especially in highly urbanized areas like Metro Manila. In this regard, seismic hazard assessments are conducted to aid structural engineers in designing earthquake-resistant structures and ensure that structures are within a desired level of performance.

Probabilistic Seismic Hazard Analysis (PSHA) is a mathematical approach that accounts for the stochastic nature of earthquakes and gives the future ground motions that may occur at a given site. PSHA is composed of five steps (Baker, 2015): (1) identification of sources of potentially damaging earthquakes, (2) characterization of the distribution of earthquake magnitudes, (3) characterization of source-to-site distances of potential earthquakes, (4) prediction of the resulting ground motion intensity, and (5) combining uncertainties in earthquake size, location, and intensity using the total probability theorem.

From a different perspective, Musson (2012) pointed out that the results of PSHA may be validated using a quasi-observational approach. Monte Carlo simulation (MCS) can be used to generate a synthetic catalog that is compatible with the properties of the seismicity of a site and the known relationship between ground motion, magnitude, and distance. The simulation can be repeated with a large number of trials to produce pseudo-observational data and reduces the problem to a simple counting of outcomes.

A historical earthquake catalog from the years 1900 - 2022 was obtained from the United States Geological Survey (USGS) database for a building in Metro Manila and the statistics of major earthquakes (Mw ≥ 5.0 and distance ≤ 400 km) are presented in this notebook. The framework of Monte Carlo Simulation for seismic hazard assessment in Python is also highlighted to provide an alternative means to the conventional PSHA by which the results of MCS can serve as baseline values of earthquake intensities.

Conceptual Discussion¶

Monte Carlo simulation can be used to generate synthetic earthquake catalogs via random sampling of observed seismicity. The expected ground motion is estimated from the synthetic catalog and empirical attenuation models, from which seismic hazard is derived. Random variables, such as moment magnitude and source-to-site distance, are simulated by generating uniformly distributed numbers between 0 and 1 and obtaining the corresponding value of the random variable using its CDF plot. The figure below shows the block diagram of Monte Carlo Simulation for earthquake hazard assessment.

Kanno et al. (2006) developed attenuation models for peak ground acceleration (PGA), 5% damped response spectral acceleration (SA), and peak ground velocity (PGV). Regression models were developed from strong ground motion records of Japan from 1963 to 2003 and 12 earthquakes from other countries to improve the number of near-source data. The ground motion prediction equations (GMPE) were limited only to three parameters: moment magnitude (Mw), source-to-site distance (X), and focal depth (D).

Seismic hazard curves are developed for PGA and SA for periods of up to five seconds based on an idealized shear wave velocity (VS30) of 760 m/s. The uniform hazard spectrum for a return period of 43 (50% in 30 years), 475 (10% in 50 years), and 2475 (2% in 50 years) years are developed which can be eventually used for structural design, risk engineering, or comparative studies with those of conventional method.

Codes¶

Import Libraries¶

The following libraries are needed for this notebook:

• numpy, scipy, math - mainly used to handle statistical models, mathematical computations, and random number generators for the simulation

• pandas - data manipulation and storage of multi-dimensional arrays

• folium, seaborn, matplotlib - data visualization of geographical distribution, histograms, probability mass function, and cumulative distribution function of random variables

• warnings - to control warning messages that occasionally arise

The United States Geological Survey (USGS) is an agency that monitors natural hazards around the world including earthquakes. They have a vast data of earthquake records which seismologists, engineers, and scientists can use for statistical studies.

Data downloaded from: https://earthquake.usgs.gov/earthquakes/

Read csv file using pandas:

One of the initial steps when importing a dataset is to study its attributes to determine if there would be some form of data cleaning, filtering, or preprocessing before the actual analysis.

The data frame consists 817 earthquake records with 22 attributes or columns.

There are only certain attributes that are needed in the study including: time, latitude, longitude, depth, mag, magType, and place.

The information above shows that except for place, the attributes mentioned above does not contain null values and allows us to proceed without performing any data preprocessing.

Mapping earthquake records¶

The project site is located in Quezon City (Latitude: 14.628056ºN, Longitude: 121.068611ºE) and mapping the earthquake records will give the seismicity of the project location in terms of geographical distribution.

In the cell below, a new dataframe is created to isolate columns that will be used for mapping as well as the details that will be displayed.

A function named magnitude classifier will cluster the earthquakes into five categories. A new column is added to the map_df dataframe to classify the earthquakes into feature_groups that will enable the isolation of these specific groups in the folium map.

The .map() method will apply the function to a series in the dataframe.

Folium Map¶

Base Map

• to have a scalable circle marker in folium, scalable vector graphics (SVG) was defined using a snip of HTML code to define its properties. The svg rescales when zooming in or out of the folium map which provides better presentation than the basic circle marker.

Feature Groups

• to group similar terms, the .groupby() method of pandas will initially group the rows by the selected attribute by which unpacking and iteration can be performed.
• for each record, a scalable circle marker is created with varying sizes and colors depending on the magnitude group while assigning the magnitude as the tooltip and time, place, and depth for the popup. Each record is added to the feature group and the feature group is added to the map.

Layer Control

• the feature groups are aggregated in the layer control and the user can isolate magnitude groups of interest.

Statistics of Major Earthquakes¶

The following visualizations show the distributions of several earthquake parameters that will be later used for Monte Carlo simulation. The earthquake catalog is assumed to be declustered of foreshocks and aftershocks such that no dependent events exist in the records.

Earthquake Rate¶

Poisson distribution gives the discrete probability distribution of a given number of independent events with a constant mean rate over a specified time. Wang et al (2014) verified the Poisson model to be suitable for earthquake engineering practices. The poisson distribution is given by the formula:

where v is the mean annual rate of earthquakes and x is the number of earthquakes per year.

In the code below, mean_rate is determined by dividing the total number of earthquake records obtained by the number of years.

The probability mass function (PMF) and cumulative distribution function (CDF) of the earthquake rates are obtained using the poisson function of scipy library. Since the mean rate calculated is around 6 earthquakes per year, the probability for 0-15 earthquakes are obtained using the code below.

In the succeeding visualizations, seaborn library is used to graph the statistics either by barplot, lineplot, or regplot and will be modified using the matplotlib library.

Plot of Earthquake Rate¶

The sns.barplot() function is used to initially graph the PMF of earthquake rates and the succeeding codes are matplotlib modifiers.

The plt.subplots() function returns a tuple containing the figure and axes objects which can be unpacked. It provides a more concise assignment of figure and axes than the plt.figure() method call. Method calls such as ax.set_title(), ax.set_xlabel, and ax.set_ylim are used to modify the seaborn plot for better visualization. In line 14, ax2=ax.twinx() allows us to overlay a new plot by assigning ax2 as a twin axis of ax. In this way, a new plot may be used with the y-axis located on the other side of the plot.

The plot of earthquake rate shows that the most likely number of earthquakes per year in the project site is six earthquakes per year.

Magnitude¶

Magnitude scales in the dataset are varying and conversion of each value to a single scale is needed for consistency. From the results of .value_counts() method, the three magnitude scales used are surface-wave magnitude (ms), body-wave magnitude (mb), and moment magnitude (mw).

Ground motion models are usually defined in terms of the moment magnitude which requires the conversion of other magnitude types. Tang et al. (2016) performed General Orthogonal Regression (GOR) method to develop empirical equations for the conversion of surface-wave and body-wave magnitude to moment magnitude based on the earthquake catalog of Western China from the China Earthquake Data Center (CEDC). The empirical equations are given as follows:

• $M_w = -0.55 + 1.16 * mb$
• $M_w = 1.61 + 0.69 * ms$

The function below is an implementation of the conversion formulas above.

The np.vectorize() method is used to transform functions to become numpy-aware and provides faster implementation compared to other methods (e.g. .apply()). Since the initial function expects a single value of magType and mag, an array can now be an input due to the vectorization of function.

Random samples are obtained using .sample() to see the new values of the magnitude in terms of Mw scale. The dataframe below shows that most of the magnitude increased when converted to Mw scale.

A frequency distribution of the dataset can be obtained by using the np.histogram() function which returns the count of instances for a specified interval. The intervals are assigned in the bin argument of the function and for this case, the bins are from magnitude 5.0 to 8.0 with 0.1 intervals.

The output of np.histogram() can be unpacked into freq and mag_bins in order to calculate the cumulative distribution function using the frequency.

The CDF of a random variable is calculated by dividing the cumulative sum over the total number of earthquake records.

Plot of Earthquake Rate¶

The magnitudes discretized at 0.1 are plotted using the sns.histplot() function. Initially, the tick marks of the graph are left-aligned and subtracting half of the interval (i.e. 0.1/2=0.05) in the bins argument will center the tick marks. The shrink argument reduces the size of the bars for better readbility of the graph.

The frequency histogram shows that at the event of an earthquake, there is a high probability that the earthquake would have a magnitude of 5.3.

Distance¶

The distance of an earthquake event from a specific site can be obtained using the Haversine formula. Haversine formula calculates the distance between two points on a sphere given the latitude and longitude of each.

The next lines of code are similar to the implementation of statistics of magnitude. Functions such as np.vectorize(), np.histogram(), and plt.subplots() are used to determine the PDF and CDF of earthquake distance.

Plot of Earthquake Distances¶

Focal Depth¶

The focal depth can be classified into two: shallow and deep events, which was adopted from the attenuation model in the study of Kanno et al (2006). The earthquake is a shallow event if the focal depth is less than 30km, otherwise it is a deep event.

The .map() method is used to classify each earthquake record in the dataframe by defining an anonymous function using the lambda function.

Plot of Focal Depth¶

The plot of focal depth shows that almost two-thirds of the earthquakes happened more than 30 km below the surface.

The series depth_df is redefined numerically in order to use the parameter for Monte Carlo Simulation.

Correlation of Magnitude and Distance¶

Before proceeding with the simulation, it is important to establish the correlation between random variables. If a correlation exists, there is a need to perform transformation techniques before proceeding with the simulation.

The code below concatenates the magnitude from the mag_df dataframe and distance from the dist_df dataframe using the pd.concat() method call. The argument axis is assigned with a value of 1 to indicate that the two series are added as columns and not as rows.

The correlation coefficient between the magnitude and distance can be calculated using the np.corrcoef() function. This function returns a matrix of correlations between variables.

The correlation coefficient is selected by indexing the array with the value rounded to four decimal places.

Plot of correlation between magnitude and distance¶

The figure above shows that there is no correlation between magnitude and distance as the value of R is close to 0.

Monte Carlo Simulation Function¶

After defining the statistics of earthquake parameters and obtaining their respective CDFs, the function that would perform the Monte Carlo Simulation is defined in the following lines of code.

The function would need three inputs namely, the period of the structure, the number of earthquake events, and the AVS30 which is a measure of a property of soil in the site. The output of the function will be a dataframe corresponding to the pseudo-earthquake records of the location.

Ground Motion Prediction Models coefficients¶

The regression coefficients of the prediction models are imported using the .read_csv() method of pandas with the period of the structure as the index.

Unpacking of coefficients from the dataframe¶

The index of each dataframe for the regression coefficients of the prediction models is the period of the structure, and the .loc[] method can be used to obtain the coefficients and unpack them to several variables.

Random generation of magnitude, source, and depth¶

A uniformly distributed number between 0 and 1 is generated using the np.random.uniform() method as shown below.

Initially, the CDF of each random variable was determined and plotted in the above visualizations. In the code below, the CDF will be used as the bins for np.digitize() method wherein np.digitize() returns the index of the random variable corresponding to that CDF.

Finally, the value of the random variable is obtained by passing the index into their respective bins and using np.around() to round off the value.

For the standard normal error, the norm.ppf() method of scipy is used to return the inverse CDF of a given value.

Monte Carlo Simulation¶

In the simulations below, the Peak Ground Acceleration (i.e. period= 0) is simulated and the effect of the number of trials is studied by varying the specified number of years.

Trial 1: 10,000 years¶

To determine the number of earthquakes that would occur within a 10,000 year period, the number of earthquake per year is randomized using the poisson distribution model. The np.random.poisson is used to perform the randomization with the mean_rate and no_of_years as arguments.

The array below shows that on the first year, there would be eight earthquakes, eleven earthquakes on the second year, eight earthquakes on the third year until the four earthquake events on the 100,000th year.

The array above is then added to determine the total number of earthquake events for a specific time period.

Since the mean rate of earthquake is around 6 earthquakes per year as defined in the above statistics, the 66915 earthquake events for a 10,000-year period is a reasonable result.

The Monte Carlo Simulation is initiated by entering the period, no_of_eq, and AVS30 in the function.

Trial 2: 100,000 years¶

Trial 3: 1,000,000 years¶

The spectral acceleration (SA) from 0 to 4.0-g (g-gravitational acceleration) is discretized by 0.02 as shown in the code below.

The exceedance rate or the probability by which a structure will experience an acceleration is determined by counting the number of earthquakes that exceeds a specific spectral acceleration.

In the code below, a lambda function is defined inside the .apply() method wherein the function will count the number of instances where a specific (SA) is exceeded in the simulations above for all values of SA. If the value of the sum is 0, a null value is assigned instead of 0 to avoid an undefined value since the plot below will make use of a logarithmic scale in the y-axis.

For better presentation, the precision of columns for the final dataframe is assigned with the .style.format() method.

Seismic Hazard Plot of Peak Ground Acceleration¶

Mulitple lineplots are added to the figure below by calling the sns.lineplot() for each data column of mean exceedance rate. The y-axis of a seismic hazard plot is usually in logarithmic scale and this is achieved with ax.set_yscale('log') method.

The plot of hazard curves with varying number of years show that at lower probabilities of exceedance (more rare events), the tail end of the graph is not as defined as those with higher probabilities. Thus, a large number of samples is required to obtain more data points at the regions of lower probability. In the above simulation, the sampling of 10,000 years of earthquake still gives a bad approximation for probabilities lower than $10^{-3}$ and the total number of years should be extended. A smoother curve is obtained when the number of trials was extended to 1,000,000 years. Thus, the same number of years will be considered in the succeeding plots for other structural periods.

Seismic Hazard Plot for Different Structural Periods¶

The Monte Carlo simulation was iterated for varying structural periods to obtain the annual hazard curves of the site. A number of structural periods were considered between 0 and 5 seconds which will be used to construct the Uniform Hazard Spectra (UHS).

Structures are designed with varying levels of earthquake hazard that are defined in terms of return period. Return period is the reciprocal of the annual frequency of exceedance and the code below adds a new column in the hazard dataframe that calculates the return period of each record.

Plot of Seismic Hazard Curves¶

Uniform Hazard Spectra¶

A Uniform Hazard Spectra (UHS) corresponds to spectral accelerations that have the same annual frequency of exccedance. Three levels of hazard are of particular interest in PSHA namely, Service-level Earthquake (43-year return period), Design-level Earthquake (475-year return period), and Maximum Considered Earthquake (2475-year return period).

In order to interpolate the spectral accelerations with equal return periods in the hazard dataframe, the pivot_table command is used to assign the structural periods as columns with the return periods as values.