To evaluate the liquefaction damage of soil, this study uses the Monte Carlo simulation to generate the combination of earthquake magnitude (M) and peak ground acceleration (PGA), which are the seismic parameters for the assessment. The simulation methods and the seismic source database are consistent with those used in probabilistic seismic hazard analysis, so the resulting parameters are also hazard-consistent. Moreover, the hazard-consistency is maintained in the fragility curves of the soil liquefaction analysis.
3.1. Procedures for Seismic Hazard Analysis
Probabilistic seismic hazards can be divided into four steps [
12], as shown in
Figure 1.
The first step is to identify the earthquake sources. As
Figure 2 illustrates, this study models the sources in the deep seismic zone (hypocentral depth > 35 km) as point sources [
13]. On the other hand, in the shallow seismic zone (hypocentral depth < 35 km), the sources are modelled as area sources or line sources (such as active faults) [
14].
The next step is to describe how often and how strongly each source produces earthquakes, or, in other words, the relationship between the frequency and the magnitude of seismic events. A maximum possible earthquake is selected for each source.
The third step is to assess the impact of the earthquake. To calculate how the peak ground acceleration (PGA) decreases with the magnitude (
M) and the distance from the site to the source (R), this study uses the Campell's formula [
15]:
where
is the horizontal
PGA value for hard site conditions;
b1~
b5 are regression coefficients with values of
b1=0.00369,
b2=1.75377,
b3=2.05644,
b4=0.12220 and
b5=0.78315, respectively [
16]. σ
lnE is the standard deviation of the natural logarithm of the model error.
The next step is to estimate the site's seismic hazard, which means calculating how likely it is that different levels of ground shaking will occur at the site in a given year. The seismic hazard analysis is based on the assumption that earthquakes follow a stationary Poisson process, which means that they occur randomly and independently over time.
If go through the above four steps, the annual probability of peak ground acceleration, PGA
R, exceeding a value
a can be expressed as [
14]:
where
n is the number of potential seismic sources in the region of the site;
is the event that the earthquake with
M ≥
M0 occurs in source
i ;
is the annual average occurrence rate of earthquakes with
M ≥
M0 in source
i; and
M0 is the smallest magnitude of concern to engineering.
3.2. Simulation of Hazard-Consistent Seismic Parameters
To produce the pair of seismic parameters that match the hazard level, the earthquake magnitude (M) and the peak ground acceleration (PGA), we first simulate the magnitude, M, of a seismic event and the distance from the site to the source, R. Then, we use Equation (14) to calculate the horizontal peak ground acceleration, PGAR, at the hard site for the given scenario of M and R. This way, we obtain the pair of (M, PGAR) as well.
Table 4 lists the random variables that the Monte Carlo simulation uses. The seismic sources can be either line sources or area sources, depending on whether they are in the shallow zone or not. For well-defined faults (Type I sources [
14]), the model uses line sources with known parameters such as position, direction and length of the fault line, for example, Type I active fault in Taiwan. The random variables for this type of source are magnitude (
M), epicentral distance (
Re), hypocentral depth (
h), and the model error term (
E). The future seismic events may start a rupture at any location along the fault line. The rupture length of the future events depends on magnitude
M.
For area sources, we can use a line source to represent them, but we have to account for the uncertainty of the line source's location and orientation (Type III source [
14]). The fault line's length depends on the magnitude
M. The random variables in this case are the magnitude (
M), the depth of the hypocenter (
h), the direction of the fault rupture (
θ), and the model error term (
E). We can divide the source areas into grid points and assume that each grid point can be the origin of a fault rupture. Therefore, we only use the depth of the hypocenter (
h) as a random variable to capture the uncertainty of the seismic event's position.
If the seismic sources are located in a deep zone, they can be approximated as point sources. In this case, the position error may exceed the rupture length, and the hypocentral distance is large enough to make the rupture length negligible. Thus, the only random variables to consider are the magnitude (M), the hypocentral depth (h), and the model error term (E).
The aforementioned random variables follow different distributions. The probability density function of earthquake magnitude,
M, is [
14]
where b is the seismic region coefficient; and are the smallest magnitude of concern to engineering (5.0 is used in this study) and the upper bound magnitude, which depends on the region, respectively.
We assume that the model error term, E, follows a log-normal distribution. The fault rupture azimuth, θ, has a uniform distribution over 2π. For the epicentral distance of Re and hypocentral depth of h, we also assume that they are uniformly distributed over the range given by the seismic source database.
3.3. Liquefaction Hazard and Fragility Curves
A diagram that illustrates the steps for assessing the risk of liquefaction is presented in
Figure 2. The process begins by randomly creating a set of seismic parameters that are compatible with the hazard level, namely the earthquake magnitude, M, and the corresponding peak ground acceleration at stiff sites, . These parameters are then used to evaluate the liquefaction potential.
In the following, after considering the soft site effects and the uncertainties of associated parameters of geological data, the factor of safety against liquefaction (FS) and the liquefaction probability (PLiq) for every depth considered are first computed. Then, the liquefaction potential index (IL), liquefaction-induced settlement (St), and the liquefaction probability index (PW) are calculated as follows to evaluate the liquefaction damage of the site subject to earthquake.
The system will only process the samples of PGAR ≥ 0.05g to assess the liquefaction damage level of the site. The soil is considered to be non-liquefiable if it has a factor of safety above 1.0 and a zero liquefaction probability for the sample of PGAR < 0.05g. The simulation runs until the number of samples of PGAR ≥ 0.05g, denoted by Ng, exceeds 100,000.
If we denote the damage index of liquefaction as
X, and assume that the occurrences of earthquakes follow the Poisson process, the annual probability of the damage index exceeding a specified state
x can be expressed as
The symbols above are the same as those in Equation (16).
Based on the flowchart of the liquefaction hazard analysis shown in
Figure 3 and the Equation (18), the damage state of the liquefaction for any specified return period T can be obtained as a reference for the engineering design. On the other hand, to meet the requirements of seismic scenario simulation and the early warning system for liquefaction damages, the fragility curves of the liquefaction of the site are another important part to be determined. As to how to build the fragility curves, the partial procedures are the same as those in
Figure 3. When the pair of seismic parameters (
M, PGA
R) is generated, the soft site effects and the uncertainties of the associated parameters of the geological data are considered, and then the damage indices (
IL,
St, and
PW) are calculated, the following equation is used to estimate the liquefaction fragility corresponding to a specified peak ground acceleration value of
a:
In which IMi and Xi are the ith Monte Carlo simulation values of PGAR and the corresponding damage index, respectively. Δa is the interval of PGAR considered. The value of 0.02 g is adopted in this study.