Principles of Probabilistic Seismic Hazard Assessment (PSHA) and Site Effect Evaluation and Its Application for the Volcanic Environment in El Salvador

Abstract

This book chapter explains the fundamental concepts of the probabilistic seismic hazard and site effect evaluation. It is divided into four parts: firstly, the theoretical background of the probabilistic seismic hazard methods is explained to compute the earthquake loads used in structural analysis of buildings, namely, the rigid-zone method, the free-zone methods, and the characteristic models. We emphasize the physical meaning of the seismic coefficient prescribed in the seismic code regulations and its association with the return period of ground motion and spectral ordinates. The interconnection of the return period, the recurrence interval, and the lifetime concepts are explained to clarify misconceptions among these terms in connection with the probability of exceedance of motion. Secondly, the seismic hazard methods are applied employing volcanic chain seismicity data, and preliminary seismic hazard maps for rock site are presented for flat topography conditions along El Salvador. Thirdly, the site effects in terms of the amplification of ground motion are studied using soil profiles characterized by the interbedding of lava flows and volcanic ashes. Finally, we present a summary that highlights the most important concepts explained in this book chapter.

1. Introduction

The probabilistic seismic hazard assessment (PSHA) involves the production of contour maps that represent the levels of ground and structural shaking expected to be experienced over a particular return period (in years). Generally, the shaking is represented in terms of peak ground acceleration (PGA) that is generated by the effects of the rupture of a geological fault and the propagation and attenuation of seismic waves reaching the bedrock. The shaking can be also represented by the structural acceleration (SA), and it depends on the period of vibration and damping of the building, and it is related to the characteristics of the ground motion itself; the SA is also referred in the literature as a spectral ordinate, and the damping is set as 5% of the critical (e.g., reinforced concrete structures) in PSHA computations. The PGA is associated with a zero structural period or with an infinite rigid structure. Then the severity of the shaking computed in a PSHA is primarily calculated on rock site conditions (NEHRP B soil class with a Vs = 760 m/s) and assuming a flat topography. The amplification of ground motion due to the presence of sediments above the bedrock is taken into consideration as a site effect in a second phase of evaluation.

The conjunction of spectral accelerations is represented in a response spectrum that is indeed the seismic coefficient C_S that is always prescribed in the codes’ regulations as a fraction of the gravity g (9.8 m/s²), i.e., 1.19 g = 11.66 m/s² at 0.2 s period (Figure 1a, b). A common misunderstanding is to assume that the seismic coefficient is related only to the ground motion or PGA; such confusion comes from the fact that the shear base V of the building is related with the seismic coefficient employing the following equation:

$$V=C_{S}W$$E1

where W is the total weight of the building. What Eq. (1) really reflects is the reaction at the ground level in the foundation of the building due to the all inertial forces F that result at each mass story subjected to its correspondent spectral acceleration SA during an earthquake. When considering a one degree of freedom system or a one-story building, the inertia force F can be calculated by the second Newton law as

$$F=\mathit{ma}$$E2

where m is the mass and a is the spectral acceleration. For the one-story building case V = F, a = C_S, and m = W/g. Note that engineers use the weight instead of mass for structural calculations, so the seismic loads are directly calculated vanishing the gravity g constant when using Eq. (2): F = (W/g)*1.19 g. In the case of having a several degrees of freedom system or a multistory building, Eq. (1) also represents the reaction at the base of the building due to all inertia forces above; its application is still valid since spectral ordinates are always computed considering a single degree of freedom system that represents the seismic demand of one mode of vibration in the structure; in other words, the shear base in Eq. (1) retains the total inertia forces that are distributed later among all the stories along the building height when considering multi-degree of freedom systems. When considering the first overtones or harmonics in the calculations, the seismic demand of each mode of vibration can be combined employing classical methods of structural dynamics based on the same response spectrum (Figure 1b) derived by seismic hazard computations.

Figure 1.

(a) Geographical location of El Salvador. (b) Elastic design response spectrum for 5% of critical damping at San Vicente City setting several years of return period.

2. Principles of probabilistic seismic hazard assessment (PSHA)

In PSHA calculations, the return period RP is associated with the expected level of shaking C_S at one place and not with the recurrence interval RI of a single earthquake related with a specific magnitude in a seismogenic source. Despite the concepts of the return period RP and the recurrence interval RI are different—although sometimes confused in the literature—they complement each other to produce the seismic loads that depend of the selected lifetime L (in years) of the building. The seismic coefficientC_S in the modern code provisions is usually calculated setting 500- or 2500-year return period which is associated with the probability of exceedance: this is the probability that the building experiment at least one time—might be two or more—a specific level of shaking in its lifetime L. Employing a Poisson distribution, the probability of exceedance q for a specific motion level can be expressed as

$$q=1-e^{-\mathit{LE}\left(z\right)}$$E3

where E(z) is the mean annual rate (number of earthquakes per year) for a specific level of motion z. Setting q to 2 and 10% and L = 50 years, E(z) yields 0.0004 and 0.002 earthquakes per year after evaluating Eq. (3). The return period is the reciprocal of E(z); RP = 1/E(z) yields 2500 and 500 years in these cases. Note that the lifetime of 50 years is presumably set as the average time span for which a person can own a building property. Note that the multiplication of L by E(z) vanishes the unit of years allowing probability results q. The usefulness of the earthquake interval concept on hazard computations is presented in Section 2.1.

Then, the main objective of seismic hazard computations is to obtain the mean annual rate E(z) for a specific level of motion z, and consequently the probability of exceedance is calculated in Eq. (3) setting a lifetime usually set as L = 50 years. For a specific period of vibration, a seismic hazard curve is calculated at several motion levels, and the design value of acceleration for a specific period of vibration is obtained at 2 or 10% of probability of exceedance via interpolation of the hazard curve (Figure 2). When considering another period of vibration in the hazard computations, it is possible to construct a response spectrum related to the engineering project at the return period of interest. Finally, the seismic hazard maps are constructed using a grid of points smoothing the acceleration values obtained at each point that constitutes the grid over the region under consideration.

Figure 2.

Seismic hazard curves for San Vicente City (Figure 1a) depicting the number of earthquake per year vs. the peak ground acceleration and spectral accelerations for 0.2 and 1.0 s structural period of vibration.

In the next sections, the main concepts for three PSHA classical methods, namely, the rigid-zone method, the free-zone methods, and the characteristic model are discussed.

2.1. The rigid-zone method

The rigid-zone method is the standard Cornell-McGuire approach [1, 2] based on a Poissonian probability distribution for which the occurrence of an event is independent of others. The first step in PSHA is to quantify the level of seismicity in the region of study answering three questions: (i) where do earthquakes happen? (ii) what magnitude do they have? and (iii) how do they frequently occur? Several ith seismogenic sources must be declared as polygons that enclose a group of earthquakes in particular geographical location that have similar geological characteristics, namely, the tectonic environment, the focal mechanism, the focal depth, and a maximum magnitude. As a result, that area constitutes a rigid zone (Figure 3). It is called the rigid zone because it is assumed that the level of seismicity is constant all over the area inside the polygon and cannot change with time and dimensions. Since the hazard analysis is based on a Poisson probability distribution, the compiled earthquake catalog must be declustered, that is, the foreshocks and the aftershocks must be removed from the list of events, and only the main shocks must remain in the catalog; this is usually accomplished employing decluster algorithms [3]. To quantify the seismicity in a specific seismogenic source, a Gutenberg-Richter (G-R) recurrence relationship is defined as

$$\mathit{log}{N}_i=A-\mathit{BM}$$E4

where N_i represents the cumulative number of earthquakes per year above magnitude M, and A and B set two constants derived from regression analysis; such constants are independent of the period of vibration of buildings. The number of earthquakes N_i within different magnitude ranges are usually computed employing the bounded G-R relationship based on [4]. The constant -B is the measure of the relative abundance of large to small shocks, and it includes the negative sign indicating the negative slope that must result from the regression analysis. Note that Eq. (4) must be derived after a completeness analysis of the earthquake catalog [5, 6]; what it means is that the cumulative number of earthquakes must be divided by the number of years for which a particular magnitude bin is completed. It is worth mentioning that the reciprocal value of N_i yields a recurrence interval RI and does not yield the return period for a certain magnitude M. For example, if A = 5.2 and B = 1 for a M 4.0 and M 7.0, N_i yield 16.0 and 0.016 earthquakes/year, and the recurrence interval are ¾ and 63 years, respectively. Then the return period is associated with the level of motion at a particular site, and the recurrence interval is related with the size of an earthquake [7].

Figure 3.

Epicenters (solid circles) and geometrical delimitation (polygon) of volcanic chain earthquakes in El Salvador and neighboring countries 1528–2009. The active quaternary volcanoes are depicted with white triangles. In this Figure the epicenters of main events, foreshocks, and aftershocks are included.

Although this information is valuable in terms of occurrence characteristics of the earthquakes sources, such recurrence interval is still not associated with a specific ground or structural motion needed to design a building. In the first instance, a ground motion prediction equation (GMPE) must be developed in order to relate the expected shaking at one place of interest depending on the magnitude and distance from the source at rock site conditions. Such GMPE must be compatible with the same tectonic environment of the seismogenic area ith to develop Eq. (4). A general form of a GMPE is

$$\text{log}\left(y\right)=a+\mathit{bM}+\mathit{cR}+\mathit{\text{dlogR}}\pm {\sigma }_{\text{log}\left(y\right)}$$E5

where y is the PGA or spectral acceleration SA in g, a and b are regression constants related to the source, and c and d are constants related with the anelastic and geometrical attenuation respectively (path term). The logarithm nature of this equation is presented in Section 3. M is the moment magnitude, R is the distance from the project site to any point of the seismogenic source i, and σ_log(y) is the standard deviation of the regression in terms of the logarithm value of y. Employing Eq. (5) and setting a = 1.68, b = 0.30, c = −0.01, and d = −1.0, the estimated PGA for an earthquake of M = 7.0 at a distance R = 10 km gives 479 cm/s² = 0.49 g. Despite that such value could be used for the design of a short period structure, it is not associated with any recurrence interval RI of this earthquake size and with any lifetime L and return period RP. Moreover, there is no information of other magnitude and distance pairs that can affect the project site; besides, the standard deviation σ is not being taken into account in the ground motion estimation at this stage.

The mean annual rate E(z) for a specific motion level “z” can be obtained combining the seismicity evaluation and the GMPEs employing the following integrals considering all possible magnitudes and distances probabilities:

$$E\left(z\right)={\displaystyle \sum _{i=1}^{\mathit{Ns}}}{\nu }_i{\int }_{R=0}^{\alpha }{\int }_{M_{\mathit{min}}}^{M_{\mathit{max}}}{f}_i\left(M\right)f_i\left(R\right)P\left(Z>z|M,R\right)\mathit{\text{dRdM}}$$E6

where Ns = number of seismogenic sources; ν_i = mean rate of occurrence of earthquakes above the minimum magnitude for the “ith” source; f_i(M) = probability density distribution of magnitude within the “ith” source, which is obtained using the G-R relationship; f_i(R) is the probability density distribution of epicentral distance “R” between various locations within source “ith” and the site where hazard is estimated; and P(Z > z |M,R) = probability that a given earthquake of magnitude “M” and source-site distance “R” will exceed motion level “z”, which is obtained employing the selected GMPE and the standard deviation σ.

Events with a lower magnitude M_min (i.e., 5.0) might not cause damage to structures, so seismic hazard computation is performed establishing a lower bound limit in all cases. The mean rate of occurrence ν_i is evaluated using Eq. (4) at M_min and represents all earthquakes per year that could happen in the seismogenic source above the lower bound limit nearby the project location, including the maximum credible magnitude M_max that can produce each seismogenic source ith; M_max is considered as the upper-bound physical limit in hazard computations. The units of E(z) are in accordance with the units of ν_i and must be multiplied by the results of the integrals that contain the unitless quantities represented by the probabilities of occurrence of magnitudes, distances, and the level of shaking z produced by all possible magnitude-distance pairs that affect the site of interest. The multiplication of f_i(M), f_i(R), and P(Z > z |M,R) obeys to the fact that the occurrences of these probabilities are physically and statistically independent.

The probability distribution of magnitude M or the cumulative probability function CDF(M) can be expressed as

$${\mathit{CDF}}_i\left(M\right)=\frac{N_{\mathit{\text{Mmin}}}-N_M}{N_{\mathit{\text{Mmin}}}}$$E7

The CDF_i(M) is the quotient between the numbers of earthquakes of a magnitude M that exceeds the number of earthquakes of minimum magnitude M_min after evaluation in Eq. (4). However, the derivation of Eq. (7) with respect to the magnitude M is used in seismic hazard calculations, yielding the probability density distribution after Neperian logarithm transformation of the B value in the G-R relationship [8]:

$$f_i\left(M\right)=\beta e^{-\beta \left(M-\mathit{\text{Mmin}}\right)}$$E8

where β = (ln10)B and e is the Neperian number. Note that the probability density function is suitable for PSHA since it represents the chances of success of the occurrence of particular size M of an earthquake. In computational practice, Eq. (8) is discretized over the same bin magnitude classes (e.g., M1 5.5 to M2 6.0) used to develop the G-R recurrence relationship (Eq. 4); then the probability P(M) of a magnitude to fall in that bin can be written as

$$P\left(M\right)={\int }_{M1}^{M2}{f}_M\mathit{dM}$$E9

where the value of M is taken as the average of the bin magnitudes M1 and M2.

When performing PSHA calculations, the polygon that represents the rigid zone is usually discretized in squares of 1 km side which represent all possible earthquake sources. The probability f_i(R) that an earthquake source could happen at a specific distance R from the project site is the quotient between R and all the possible distances R_j:

$$P\left(R\right)=f_i\left(R\right)=\frac{R_j}{{\displaystyle \sum _{j=1}^n}{R}_j}$$E10

where n is the number of elements used to discretized the seismogenic source. Note that all the earthquake sources can have magnitudes between the lower- and upper-bound limits with related probabilities of occurrence explained above. So, f_i(M) and f_i(R) are independent of the level of motion being evaluated; f_i(R) is dependent on the location of the project relative to the seismogenic source ith into consideration. The probability P(Z > z |M,R) that a given earthquake of magnitude M and source-site distance R will exceed motion level z is evaluated using a standard normal variable G:

$$G=\frac{\text{log}\left(z\right)-\text{log}\left(y\right)}{{\sigma }_{\text{log}\left(y\right)}}$$E11

To estimate log (y) in Eqs. (5) and (11), the value of M to be used in Eq. (5) is taken as the average of M1 and M2 in Eq. (9). P(Z > z |M,R) is then calculated as

$$P\left(Z>z|M,R\right)=1-{\varphi }_{\left(G\right)}$$E12

where ϕ_(G) is the probability density function derived from a log-normal Gaussian distribution evaluated for the G value obtained in Eq. (11). For computational purposes, the exact integrals in Eq. (6) are substituted by discrete sums yielding:

$$E\left(z\right)=\sum _{i=1}^{\mathit{Ns}}{\nu }_i\sum _{R_{\mathit{min}}}^{R_{\mathit{max}}}\sum _{M_{\mathit{min}}}^{M_{\mathit{max}}}P\left(M\right)P\left(R\right)P\left(Z>z|M,R\right)$$E13

The probabilities in Eq. (13) are calculated after selecting a segmentation of magnitude and source-site distance pairs. Usually such segmentation is done for every 0.25–0.5 magnitude units of magnitude and 1–10 km distance. R_min and R_max constitute the minimum and maximum source-site distances.

2.2. The free-zone methods

Woo [9] proposed an alternative approach to compute the PSHA in order to remove the uncertainties inherent in the definition of the rigid seismogenic sources above. The arbitrariness in delineating zone boundaries may be artificially expanded to dilute activity around a site or on the other hand artificially contracted to concentrate activity around a site. Besides, the assumption that the activity rate is constant in a seismogenic source is rarely tested statistically [9]. In this method, the seismic hazard is calculated directly based on the characteristics of the earthquake catalog rather than defining rigid zones considering a spatial uniformity seismicity; then earthquake sources are not defined according to a geological and tectonic criteria. The Woo method is called the kernel estimation method since it implements spatial smoothing directly to earthquake epicenters contained in the catalog depicting the spatial nonuniformity of seismicity in the region of study. Then the earthquake catalog must contain the main events, the foreshocks, and aftershocks when performing hazard calculations.

The first step in the method consists in defining a grid of nodes that constitute the point earthquake sources around the project site. Each point of the grid would have an activity rate that is magnitude-distant dependent based on fractal characterization. A multivariate probability density function K is expressed in terms of magnitude and distance yielding:

$$K\left(M,x\right)=\frac{n-1}{\pi H^2}{\left[1+{\left(\frac{R}{H}\right)}^2\right]}^n$$E14

where n is the exponent of the power law (usually between 1.5 and 2.0), R is the source-site distance, and H is the bandwidth for normalizing distances that is a function of magnitude M representing the average distance between epicenters of the same magnitude bin:

$$H=c{e}^{\mathit{dM}}$$E15

where c and d are constants to be derived by a regression analysis. By summing over all cataloged events with historical epicenters x_i, the cumulative activity rate density is computed at a general point x in the grid for each magnitude bin class from the minimum magnitude to the highest value, as follows:

$$\lambda \left(M,x\right)={\displaystyle \sum _{i=1}^{n_o}}K\left(M,x-x_i\right)/T\left(x_i\right)$$E16

T(x_i) is the effective observation time period which is equivalent to the completeness period, and n_o is the number of earthquakes in the catalog; the contribution of each event is inversely weighted by its completeness period. Once λ(M,x) is calculated for each node in the mesh, a GMPE is employed, and the seismic hazard is calculated by summing over each node source for an specific motion z. Note that the probability P(Z > z |M,R) will be evaluated using Eqs. (11) and (12) and multiplied by λ(M,x) to calculate E(z) and the correspondent probability of exceedance q in Eq. (3) setting a specific lifetime L. The observational uncertainties in event magnitude and epicentral location can be accounted for in this method employing Gaussian error distributions [10]. The basic principle underlying the Woo method is that the epicenter of each past shallow event is smoothed geographically to generate a spatial probability distribution for event recurrence. As a conclusion, the method intends to represent both the fractal clustering of moderate magnitude epicenters and the haphazard migration of major events [9].

Frankel [11] also proposed a free-zone method used for the US National Hazard Maps developed by the US Geological Survey (USGS). In this method a region is partitioned into cells by a grid of 0.1° (11 km on a side). Within the ith grid cell, a count is made of the number of events with a magnitude above a reference value M_ref . This cumulative count of events is converted to incremental values from M_ref + ΔM or for a specific magnitude bin. The grid of n_i values is smoothed spatially in a simple Gaussian way with correlation distance c_d. The smoothed value of n_i is obtained as follows:

$$n_i=\frac{\sum _j{n}_j{e}^{-{\Delta }^2\mathit{ij}/{c_d}^2}}{\sum _j{e}^{-{\Delta }^2\mathit{ij}/{c_d}^2}}$$E17

where Δ_ij is the distance between the ith and jth cells. For a specific site, the values of n_i are binned by their distance from the site, and a value of N_r is calculated as the total value of n_i for cells within a certain distance increment of the site.

The mean annual rate E(z) for a specific motion level “z” at the project site is determined as follows:

$$E\left(z\right)=\sum _r\sum _m{10}^{\left[\text{log}\left(\frac{N_r}{T}\right)+B\left(M_m-M_{\mathit{ref}}\right)\right]}P\left(Z>z|M,R\right)$$E18

where r and m are the counters for the distance and the magnitude bins, respectively. It is worth mentioning that the term ${10}^{\left[\text{log}\left(\frac{N_k}{T}\right)+B\left(M_i-M_{\mathit{ref}}\right)\right]}$ is the value of 10^A in Eq. (4) for the completeness period T; $P\left(Z>z|M,R\right)$ is calculated as in the previous methods. Since a Poisson probability distribution is used in this method, the earthquake catalog must be declustered removing the foreshocks and aftershocks retaining only the main events for PSHA calculations.

2.3. The characteristic earthquake model

The characteristic earthquake model is employed when detailed fault information is collected by a geologist in the field and where there is a clear evidence of the slip via geodetic surveys (GPS) and/or paleoseismological investigations [12, 13, 14]. Instead of declaring an area as a rigid seismogenic source or employing fractal geometry concepts in free-zone methods, a fault is defined employing the following input data:

1. The slip rate (mm/year)

2. A B value (usually set as 1.0)

3. Precise geographical coordinates of a series of points that define the fault trace

4. The dip, the sense of slip, and depth of the fault

5. The minimum magnitude M_min (e.g., 6.0) to be considered in the analysis

Two alternate characteristic models are applied actually in a PSHA:

1. Considering only the maximum magnitude M_max and its associated activity rate or number of earthquakes per year

2. Classical Gutenberg-Richter (G-R) relationships

For both models the first step in a PSHA is to calculate the number of earthquake per year E(z) for the M_max in model (i) and for the M_min in the model (ii).

The idea of using two alternate models is that it is unknown whether future earthquakes will be large fault-filling ruptures of smaller ruptures which only occupy a relatively small portion of the fault [15]. Generally a weight of 0.5 is assigned to each model if the mean recurrence time for M_max in the first model is less than the return period of motion.

2.4. Probabilistic seismic hazard application for the volcanic chain earthquakes in El Salvador

In this section, the free-zone method, the rigid-zone methods, and the characteristic models are applied for the volcanic chain earthquakes in El Salvador. For the Cornell-McGuire method, an own program was developed and used the KERFRACT program for the Woo method. The programs developed by the USGS [15] were employed for the Frankel and characteristic models. The seismic activity of this source is concentrated in upper part of the crust with focal depth less than 25 km and within a continuous belt of 20–30 km width along the active volcanoes in Central America (Figure 3). The events can appear in clusters. Faults with right and left lateral focal mechanisms are located parallel and perpendicular to the volcanic axis, respectively. The last destructive earthquakes inside the volcanic chain in El Salvador dated on October 10, 1986, M5.7, and February 13, 2011, M 6.5 [16].

Seismic hazard maps for PGA were developed employing a 0.1 × 0.1${}^{°}$ grid assuming a flat topography setting 500-year return period of ground motion in terms of PGA. The homogenized earthquake catalog compiled by Salazar et al. [17] was used for the moment magnitudes above 5.0 (M_min) and the GMPE developed by Sadigh et al. [18] for shallow crustal earthquakes and a truncation value of 3σ in the hazard computations. We updated the work of Bommer et al. [19] which also presents PSHA maps for El Salvador in terms of normalizing PGA to the maximum computed value; here the absolute values of PGA are presented based on new seismological information [17] applying the methods explained above.

A M_min of 5.0 and M_max of 6.9 were set when using all methods with a constant seismogenic depth of 5 km. The G-R relationship for the rigid-zone method is presented in Figure 4 after the completeness analysis. A tri-linear trend that reflects interesting features of the seismicity along the volcanoes is observed. The first part (B1) suggests a large number of small earthquakes related to big shocks, the middle part yields a very low B2 value reflecting the effects of clusters or the abundance or moderate size events, and the third part (B3) reflects the small number of destructive events related to small shocks. The tri-linear trend of the G-R relation yields the probability density function when applying Eqs. (7) and (8). The sum of all magnitude bin probabilities of the tri-linear trend yields 76, 9, and 15% for the first (B1), second (B2), and third (B3) segments, respectively, which sum must be 100%. Although, it is interesting to note that the probability of occurrence of an earthquake of a specific magnitude M 5.6 and M 6.3 is the same of about 3%, presumably, the contribution to the seismic hazard from both events in terms of magnitude probability is the same despite the size differences of these events (Figure 5). Generally, the greater the magnitude, the lower the probability. The probability of occurrence for magnitudes between 5.7 and 6.1 yield similar (2%) due to the low B2 value in the G-R relationship.

Figure 4.

Gutenberg-Richter (G-R) relationships for the volcanic chain earthquakes. Solid line and solid circles: tri-linear trend employed in the rigid-zone Cornell-McGuire method. Broken line: linear regression employed in the free-zone Frankel method with B value of 0.97; see Eq. (4).

Figure 5.

Magnitude probability density function for the volcanic chain earthquakes based on the tri-linear G-R relationships depicted in Figure 4; see Eq. (8). The percentages are the sum of all magnitude bin probabilities for each segment of the tri-linear trend with a total of 100%.

Carr et al. [20] suggest that the volcanoes in Central America cluster into centers, whose spacing is random but averages about 27 km. Table 1 and Figure 6 present the average distance of earthquakes within magnitude bins of 0.25; the all average distance for all size bins yields the same 27 km as the volcano centers suggested by Carr et al. [20]. What it means is that earthquake epicenters and volcano centers have the same spatial distribution among them. Regarding with the Woo free-zone method, the kernel only depends on the epicentral distance rather than magnitude, so the kernel function must be unique for the volcanic chain events (Figure 7) and independent of the size of the earthquakes; then c = 27 and d must be zero yielding H = 27 km in Eq. (15).

Magnitude range	Number of earthquakes	Distance (km)
5.00 ≤ M < 5.25	68	15.84
5.25 ≤ M < 5.50	13	35.19
5.50 ≤ M < 5.75	20	22.33
5.75 ≤ M < 6.00	14	24.79
6.00 ≤ M ≤ 6.25	22	19.44
6.25 ≤ M ≤ 6.50	11	22.64
6.50 ≤ M ≤ 6.75	5	47.98
6.75 ≤ M ≤ 7.00	2	27.75
Total	155	Average H = 27

Table 1.

Average distances for magnitude bins for volcanic chain earthquakes.

Figure 6.

Bandwidth for normalizing epicentral distances in El Salvador volcanic chain and neighboring countries where H = 27 km in Eqs. (14) and (15). Carr et al. [20] suggest the volcanoes in Central America clusters into centers, whose spacing is random but averages about 27 km. It is noted that the all average distance for all earthquake epicenters also yields the same 27 km as the volcano centers (see Table 1).

Figure 7.

Kernel function for upper-crustal earthquakes in the volcanic chain. Note that the kernel probability function K is unique for H = 27 m in Eqs. (14) and (15).

We model the characteristic earthquake for the San Vicente Fault employing the following parameters:

1. The slip rate: 4.1 mm/year [21] based on paleoseismological investigations

2. A B value equal to 1.0

3. The geographical coordinates that define the fault trace are taken from [21]

4. The dip = 90°, the sense of slip as a pure right lateral and depth of the fault = 15 km

5. The minimum magnitude M_min as 6.0 to be considered in the analysis.

We present in Figures 8–11 the seismic hazard maps employing the four methods explained above. As a first observation, it is noted that the rigid-zone method yields a uniform level of ground motion along the volcanoes and that the free-zone methods can capture the nonuniformity of hazard spatial distribution especially at the western part of the country. However, at San Salvador City (see location in Figure 1a), the PGA yields about 0.45–0.5 g for all methods. At San Vicente City (see location in Figure 1a) when applying the rigid zone, the Woo free zone, and the characteristic model, the PGA yields 0.45 g; however, a lower value of about 0.35–0.40 g yields for the Frankel free-zone method. This difference is attributed to the fact that the Woo method employs the unclustered catalog and the Frankel method employs the declustered one that it only contains the main events; it is noticed also that for the Frankel method, we employed a unique B value of 0.97 for a G-R relationship (Figure 4); however, such hazard estimates are not significantly different from each other. When applying the characteristic model employing the paleoseismological slip rate at San Vicente City of 4.1 mm/year, a near value of 0.45 g is obtained as well. The coincidence of the PGA values revealed the quality of our earthquake catalog in the volcanic chain. The actual state of the art suggests to treat the epistemic uncertainties in the PSHA through a logic-tree framework considering all the methods explained in this chapter and compute the covariance or error associated with the result of each branch in the logic tree [22, 23, 24, 25, 26]. The PSHA could be performed for spectral ordinates as well; Figure 1b depicts an elastic response spectrum for San Vicente City employing the Woo free-zone method for 13 periods of vibration between 0.0 and 4.0 s setting 5% of critical damping.

Figure 8.

Seismic hazard map for the peak ground acceleration (PGA) as a fraction of gravity “g” employing the rigid-zone Cornell-McGuire method setting 500 years of return period.

Figure 9.

Seismic hazard map for the peak ground acceleration (PGA) as a fraction of gravity “g” employing the free-zone Woo method setting 500 years of return period.

Figure 10.

Seismic hazard map for the peak ground acceleration (PGA) as a fraction of gravity “g” employing the free-zone Frankel USGS method setting 500 years of return period.

Figure 11.

Seismic hazard map for the peak ground acceleration (PGA) as a fraction of gravity “g” employing the characteristic model at San Vicente Fault setting 500 years of return period and 4.1 mm/year of slip rate.

Figure 12 shows a disaggregation chart depicting the contribution to the mean annual rate of exceedance (or number of earthquakes per year) for a specific ground motion value at a site due to scenarios of a given magnitude and source-to-site distance. Then it is possible through the disaggregation analysis to determine the magnitude-distance pair with the largest contribution to the hazard [27]. This analysis is performed setting a PGA of 0.46 g observed at San Vicente City (Figure 15) when applying the rigid-zone method. It is interesting to note that small earthquakes of M 5.0 have the largest contribution and that an earthquake of M 6.2 has the same contribution of an earthquake of M 5.4 to the seismic hazard in the volcanic chain of El Salvador. We attribute such especial features due to the tri-linear trend observed in the seismicity evaluation represented by the G-R recurrence relationships and its correspondent magnitude probability density function presented in Figure 5.

Figure 12.

Example of magnitude-distance pair disaggregation charts for the peak ground acceleration (PGA) of 0.46 g at 500-year return period employing the rigid-zone Cornell-McGuire method at San Vicente City (see location in Figure 1a). Note that an earthquake of M 6.2 has the same contribution to seismic hazard of an earthquake of M 5.4 in the volcanic chain of El Salvador.

3. Site effects related with the interbedding of sediments and lava flows

Different studies carried out after the occurrence of several destructive seismic events have confirmed that the presence of sedimentary deposits exerts a particularly strong influence on the distribution of damages. These studies have also corroborated that the spatial distribution of damages correlates on the magnitude of the amplifications of the arriving seismic waves caused by the upper sedimentary layers. Consequently, we can regard the distribution of the ground-shaking intensity as a phenomenon closely related to the filtering effects of the soil profile above the bedrock level [28].

It is convenient to introduce the most general and widely accepted formulation of the ground motion phenomenon, departing from assuming that in the frequency domain, the S-wave can be expressed in terms of Fourier amplitude spectra A_ie(f) recorded at a site i due to an earthquake e as

$$A_{\mathit{ie}}\left(\mathrm{ }f\mathrm{ }\right)=S_e\left(\mathrm{ }f\mathrm{ }\right)P_{\mathit{ie}}\left(\mathrm{ }f\mathrm{ }\right)G_i\left(\mathrm{ }f\mathrm{ }\right)$$E19

where f denotes the frequency domain and Se(f) represents the source, P_ie(f) the path, and G_i(f) the site effect term (Figure 13). When applying logarithm to Eq. (19) yields

$${\mathit{\text{logA}}}_{\mathit{ie}}\left(\mathrm{ }\mathit{f}\mathrm{ }\right)=\mathit{log}{\mathit{S}}_{\mathit{e}}\left(\mathrm{ }\mathit{f}\mathrm{ }\right)+\mathit{log}{\mathit{P}}_{\mathit{ie}}\left(\mathrm{ }\mathit{f}\mathrm{ }\right)+\mathit{log}{\mathit{G}}_{\mathit{i}}\left(\mathrm{ }\mathit{f}\mathrm{ }\right)$$E20

where the amplitude term A_ie(f) is related to PGA or and spectral ordinate SA, the source term can be converted to a + bM and the path term to cR + dlogR in Eq. (5). The convolution of the source and path effects yields the input bedrock motion B_n; in the first instance the site effect term can be defined as

$${\mathit{G}}_{\mathit{i}}\left(\mathrm{ }\mathit{f}\mathrm{ }\right)=\frac{{\mathit{B}}_{\mathit{m}}+{\mathit{A}}_{\mathit{m}}}{{\mathit{B}}_{\mathit{n}}}=\frac{\mathbf{2}{\mathit{B}}_{\mathit{m}}}{{\mathit{B}}_{\mathit{n}}}$$E21

where B_n and A_n are the input and refracted bedrock motion, respectively, and B_m and A_m are the incident and refracted wave at the surface, respectively (Figure 13). Note that the free-surface effects yield an amplification of 2.0 at the outcrop engineering bedrock since B_n = A_n. Such formulation is useful when there are strong motion recordings (accelerograms) for both sediment and outcrop sites. In such cases some authors [29, 30, 31, 32, 33] have proposed several empirical techniques to evaluate the magnitude of site effects on ground motion without relying on specific soil profile information.

Figure 13.

Illustration of theoretical transfer functions at soil sites and outcrop site (reference site) taking into account the free-surface amplification. B_n and A_n: input and refracted bedrock motion, respectively. B_m and A_m: incident and refracted wave at the surface, respectively.

The estimation of the magnitude of site effects applying theoretical methods requires knowledge of the physical properties of the layers composing the soil profile, namely, shear wave velocity, density, and thickness of the layers of the upper 30 meters. In most practical cases involving thick sedimentary deposits, knowledge of the thickness and physical properties of the upper 30 meters of the soil profile appears to provide a rather poor representation of the problem.

The amplification function employing 1-D SH wave propagation can be obtained when solving the following partial equation of wave motion:

$$\frac{{\partial }^{2}u}{\partial t^2}=v_s^2\frac{{\partial }^{2}u}{\partial v_v^2}$$E22

where u is the horizontal displacement, v_s is the shear wave velocity, t is time, and v is the vertical direction [8].

A theoretical 1D SH amplification function is derived solving Eqs. (21) and (22) for a soil profile volcanic environment at the capital of El Salvador, San Salvador. The work of Salazar and Seo [34] is revisited analyzing a volcanic ash and rigid lava flows’ interbedding patterns (Figures 14 and 15). A value of v_s = 225 m/s is used for the volcanic ashes in the upper part (Tierra Blanca) and a value of v_s = 500 m/s for the soils under laying the stiff lava flow (consolidated tuffs). In the case of the stiff lava flow, we used a v_s = 2100 m/s. The interbedding of soft and rigid materials is the consequence of different volcanic eruptions through historic times. The results of the 1D amplification functions are presented in Figure 16.

Figure 14.

Interbedding of lava flow and volcanic ashes in San Salvador City (see location in Figure 1a).

Figure 15.

Profiles of lava flow and volcanic ashes interbedding in San Salvador City; see location in Figure 1a (after Salazar & Seo (2002)). The value of thickness of the volcanic ash (white part) and the lava flow (black part) is written to the right of the profiles in meters.

Figure 16.

1D SH wave transfer functions for interbedding of lava flow and volcanic ashes in San Salvador City; see location in Figure 1a. (a) Site 1 and (b) Site 2 in Figure 15. Note that the second mode is predominant at 0.4 s rather than the first one for the Site 1 and that similar amplification factors of about 10–12 are obtained for the first and second modes for Site 2.

For the site 1, the highest amplification is observed at the second mode of vibration, and it does not take place at the first mode as can be usually expected, so in this case, the predominant period is the second mode of soil vibration at 0.4 s rather that the first one. For the site 2, similar amplification factors of about 10–12 are observed at the first mode and the subsequent overtone. Note that for long period components above 2.0 s the amplification factors yield a value of 2.0 confirming that the earthquake motion at these periods yields no amplification. These amplification features might have serious implications in the distribution of damages since several modes can be highly excited during an earthquake affecting medium- and low-rise buildings in San Salvador City. Moreover, if the building foundations are being constructed on the first lava flow encountered as a recommendation from a standard penetration test (SPT), the site would be characterized as a rock yielding a lower seismic coefficient despite that volcanic ashes might underlay the lava several tens of meters.

4. Summary

To finalize this book chapter, the most important concepts are summarized below when performing a PSHA:

1. The necessary input data to start a PSHA is the conformation of a homogeneous earthquake catalog in terms of moment magnitude and a strong motion data base. The earthquake catalog serves to compute the classical G-R relationships that depend on the completeness analysis of the catalog itself for a specific seismogenic source when applying the rigid-zone method. New methods are oriented to capture the uncertainty of the seismicity rate in the G-R relationships [35]. The free-zone methods test directly the earthquake catalog itself containing in the case of the Woo method, the foreshocks, the main events, and the aftershocks, and its application also depends of the analysis of completeness of the catalog. On the other hand, the strong motion data base serves to conform ground motion prediction equations (GMPEs) for both, the peak ground acceleration and spectral ordinates, or to test the applicability of another GMPE derived for similar tectonic environments.

2. In all seismic hazard methods, the first objective is to compute the number of earthquakes per year that a single site can expect during the lifetime of the structure. The combination of a lifetime and the probability of exceedance yield the return period (in years) that is employed to establish a design seismic load of the building. Then, the return period, the recurrence interval, and the lifetime of the structure are different concepts that connect each other through probability concepts [7]; however, it is a must to employ GMPEs to estimate the seismic coefficient prescribed in the building code regulations to take into account all the magnitude-distance pairs that can affect the project site. Actual research frontiers for state of the art in PSHA suggest to test the hazard curves against record strong motions and Mercalli intensity observations [36, 37].

3. The rigid-zone methods have the advantage to assessing a seismicity level when there is a geological evidence that shocks can occur within specific area; however, one must have caution to spread such seismicity in arbitrary manner. Instead, the free-zone methods are based solely on the quality of the earthquake catalogs, and they do not depend on a priori geological information neither of a particular seismogenic source delimitation.

4. The characteristic models are usually applied when there is clear information about geological faults, namely, the slip rate (mm/years), the geographical trace coordinates, the dip, and the lock depth of the fault. However, in the Woo free-zone method, it is possible to establish a background seismicity setting in a certain area, a recurrence interval for the maximum credible magnitude that could be directionally dependent applying anisotropy features, which is indeed an equivalent approach of the characteristic model.

5. A robust seismic hazard map must be the result of a rigorous computation employing a logic-tree formulation based on all methods explained in this chapter and the combination of several GMPEs [23, 38, 39, 40], possible maximum magnitudes and slip rates. Besides, it is advisable to check the covariance or the associated error or dispersion of the final results [17, 24, 25]. When employing a logic tree, the source-to site distance definitions and the ways to use the two horizontal components of motion in the GMPEs must be compatible among them [41, 42, 43, 44, 45]. The effect of truncation of ground motion distributions at a specific number above the median (i.e., 3σ) might have a substantial influence in the results when related to a very low probability of exceedance [46, 47].

6. The PSHA is generally performed at rock site conditions (NEHRP B soil class with a Vs = 760 m/s) and a flat topography. However, the final response spectrum must take into account the site effects through the amplification of ground motion due to the presence of sediments above the bedrock; however, the topographic effects might be also influence the final motion for design at highland sites. The 1D SH wave amplification function (or transfer function) can be obtained via theoretical methods based on shear wave velocity test at the sites or employing nondestructive geophysical methods as the microtremor arrays and refraction methods.

7. When studying the amplification of ground motion, the compilation of a strong motion data base [48] and the mobile microtremor measurement surveys [29] can be combined to validate the transfer function obtained via inversion analysis of earthquake data. As a result, the elastic response spectra calculated in the PSHA at rock site condition can be amplified using suitable amplification factors for the final design of a building [49].

8. Before applying any method for PSHA or site effect computation, the inherent geological characteristics of the earthquakes and the soil profiles in the region of study must be checked carefully. We have demonstrated in this chapter that the volcanic chain earthquakes in El Salvador and the intercalation of rigid lava flows and volcanic ashes yield special patterns of seismicity [50, 51] and vibration of soil deposits, respectively, that must be included in the future seismic code regulations and design foundation guidelines for the country.

Acknowledgments

I firstly thank IntechOpen editorial board for inviting me to write this book chapter and to be a part of this important global educational project. I openly thank the book editor Valentina Svalova for the detailed review that helps to improve the manuscript. This article was sponsored by the Catholic University of El Salvador UNICAES in terms of research grants for international publishing. I thank Gordon Woo from RMS London who allows me to use his Fortran code KERFRACT. I also thank the assistance of Stephen Harmsen at the USGS who guides me to use the USGS Fortran codes available at GitHub: https://github.com/usgs/nshmp-haz-fortran. Vanessa Rivera at UNICAES helped to correct the manuscript.