Earthquake Demand Energy Attenuation Model for Liquefaction Potential Assessment

Abstract

This study introduces an attenuation model based on the strain energy approach for estimating earthquake demand energy (EDE) to evaluate soil liquefaction potential. A new method is presented to estimate the EDE at a free-field site when only one record of the ground surface acceleration is available. This method was generated after analyzing the earthquake data of 18 downhole arrays in California. The developed method was later employed for calculating the EDE values of 328 earthquake records worldwide. Results showed that several parameters affected the EDE amount, including earthquake magnitude, faulting mechanism, site-to-source distance, shear wave velocity, and peak ground acceleration (PGA). These parameters were categorized by three main functions including source, distance, and site effect functions. An attenuation model was incorporated as a result of these three functions. Finally, the demand energy of two liquefaction array sites—Port Island (PI) site in Japan and Wildlife Liquefaction Array (WLA) site in California—were predicted by the proposed attenuation model and compared to the calculated capacity energies of these sites with satisfactory results.

Introduction

Since the 1970s, researchers (e.g., Trifunac 1976, Abrahamson and Shedlock 1997, Abrahamson et al. 2008, Gregor et al. 2014, Travasarou et al. 2003) have proposed many attenuation relationships for estimating ground motion parameters, such as peak ground acceleration (PGA), peak ground velocity (PGV), response spectra, and Arias intensity.

Earthquake energy is also attenuated from source to site, and can be implemented to evaluate the liquefaction potential of a soil layer based on the energy approach. Theoretically, the energy approach method is more sound compared to the other methods. Kokusho (2013) conducted a comparative study between energy-based and stress-based methods to show the importance of the energy approach. While the stress method also has been successful to model pore pressure build-up (Baziar et al. 2011a), the energy approach has been, for various reasons, the focus of many research projects in recent years, some of which are listed below.

Energy released from the earthquake source is quantified and estimated by seismological parameters (Richter 1936; Gutenberg and Richter 1954, 1956).

There is a strong relationship between energy loss during cyclic loading and generated pore water pressure. The amount of energy loss is independent of the loading and stress path amplitudes (Nemat-Nasser and Shokooh 1979, Towhata and Ishihara 1985, Baziar and Sharafi 2011).

The main advantage of the strain energy approach, compared to previous methods, is the manner in which it considers the effects of the entire shear stress and shear strain time histories simultaneously (Towhata 2008).

The earthquake energy reaching a soil element consists of both strain energy and kinetic energy (Kokusho and Suzuki 2011, 2012). Recently, the kinetic energy concept has also been employed for assessing soil liquefaction potential (Jafarian et al. 2014). In this alternative method, the amount of earthquake kinetic energy, estimated by soil particle velocity, is considered as a replacement for strain energy to determine the earthquake demand energy (EDE) value. However, the strain energy method has always been more popular among liquefaction researchers due to its increased compatibility with the unit of energy and direct relationship with pore pressure build-up.

In contrast to surface ground motion parameters (e.g., PGA and PGV), which can be recorded directly by instruments, the amount of EDE at any site needs to be calculated. Some researchers (Zeghal and Elgamal 1993, Davis and Berrill 1998) have proposed a series of linear and nonlinear methods that help determine the shear stress-strain behavior of the soil layer at the downhole array site during an earthquake, when at least two recorded accelerations are available. Since the downhole array sites are not common in different parts of the world, researchers have not been able to study the strain demand energy based on a reasonable amount of actual data. In this study, a new method is introduced to estimate the amount of EDE at a site using only surface recorded acceleration. Using this method, EDE values at wide-ranging sites are calculated, and the effect of several parameters, such as earthquake magnitude, faulting mechanism, site-to-source distance, average shear wave velocity, and PGA on the EDE are assessed. Additionally, to account for the near-field effect in EDE values, a threshold distance R_th, proportional to earthquake magnitude is introduced, and the collected databank is divided into near-field and far-field data. Eventually, an attenuation model for estimating EDE at a site is proposed based on the categorized databank.

Definition of Demand Energy

The main concept of the energy method is to determine the values of capacity and demand. The ratio of capacity and demand, like stress or strain methods, defines the safety factor against liquefaction potential. The capacity energy, studied experimentally by many researchers (Law et al. 1990, Figueroa et al. 1994, Green et al. 2000, Baziar and Jafarian 2007, Baziar and Sharafi 2011, Jafarian et al. 2012), is the minimum required energy for liquefying a soil element, and can be determined by the cumulating area of hysteresis loops up to liquefaction triggering. Contrary to capacity energy, determining the demand energy requires seismological and field study. The amount of demand energy can be calculated using the cumulating area of hysteresis loops during the shaking period, if the soil element has enough capacity to withstand the earthquake motion and does not liquefy.

Richter (1936) related the earthquake released energy at a source to earthquake magnitude. Later, Gutenberg and Richter (1954, 1956) developed the following relations in linear and nonlinear forms:

\log E = 1.8 + 1.5 M

(1)

\log E = 9.4 + 2.14 M - 0.054 M^{2}

(2)

where E is the total radiated energy (kJ) and M is the earthquake magnitude (Richter scale). Other researchers (Law et al. 1990, Davis and Berrill 1982, Berrill and Davis 1985, Trifunac 1995) have used general form of Equation 1 and have proposed several empirical relationships among demand energy, earthquake magnitude, site-to-source distance, and initial effective overburden stress. Furthermore, Jafarian et al. (2011) numerically investigated the suitability of the energy approach for evaluating the liquefaction potential. They introduced 32 earthquake records to the base of a centrifuge-simulated model, and produced an empirical equation to estimate the demand energy in a soil layer.

Overall, the empirical equations proposed by researchers for estimating the strain demand energy are based on simplified assumptions or numerical modeling. In other words, there is no comprehensive study on the evaluation of strain demand energy based on a wide range of actual strong motion data.

Downhole Array Database

In this section, an attempt is made to develop a new method that estimates the shear stress and shear strain time histories caused by an earthquake motion, when only one accelerometer on the ground surface is available. Earthquake data from 18 downhole arrays in Northern and Southern California were utilized. Table 1 presents these earthquakes with relevant information.

Table 1

List of downhole array events used in this study (data source: CESMD)

No.	Event	Year	M _w	Station name	Dist.^∗ (km)	V_S30 (m/s)	PGA (g)	Depth of Accelerometers (m)
1	South Napa	2014	6	Corckett	19.6	345	0.960	Surface, 20.4
2	South Napa	2014	6	Vallejo	11.6	509	0.194	Surface, 17.9
3	South Napa	2014	6	Rohnert-Park	37.5	340	0.055	Surface, 10.7
4	South Napa	2014	6	Antioch	53.4	150	0.044	Surface, 6, 15, 24
5	South Napa (A)^∗∗	2014	4	Corckett	15.7	345	0.053	Surface, 20.4
6	Encino	2014	4.4	Los Angeles La Cienega	14.7	241	0.104	Surface, 18.3
7	Berkeley	2011	4	Alameda	8.6	208	0.054	Surface, 5.2, 12.2
8	Ferndal Area	2010	6.5	Eureka	54	194	0.192	Surface, 19, 33
9	Calexico Area	2010	7.2	El Centro	59.1	182	0.226	Surface, 30
10	Inglewood	2009	4.7	Los Angeles La Cienega	11.4	241	0.107	Surface, 18.3
11	Inglewood	2009	4.7	Los Angeles Vincent	21.7	149	0.057	Surface, 15, 30
12	Chino Hills	2008	5.4	Corona	21.6	349	0.156	Surface, 7.9, 21.6
13	Chino Hills	2008	5.4	San Bernardino	44.8	254	0.048	Surface, 10.7
14	Piedmont	2007	4.2	Alameda	8.1	208	0.052	Surface, 5.2, 12.2
15	Parkfield Area	2004	6.2	Turkey Flat 1	7.1	900	0.191	Surface, 24
16	Parkfield Area	2004	6.2	Turkey Flat 2	7.7	417	0.286	Surface, 11, 23
17	SanCelemente	2004	5.4	San Diego-Coronado	73	210	0.057	Surface, 11.6, 23
18	Whollywood	2001	4.4	Los Angeles La Cienega	4.3	241	0.481	Surface, 18.3
max		2014	7.2		73	900	0.960	33
min		2001	4		4.3	149	0.044	Surface

∗

Epicentral Distance.

∗∗

Aftershock.

Using the trigonometric interpolation method (Davis and Berrill 1998), the shear stress-strain hysteresis loops were extracted in two horizontal directions (E-W and N-S components) for different depths of downhole array sites between the surface and downhole accelerometer levels. Figure 1 depicts the E-W component of the Corona site hysteresis loops at a depth of 3.95 m, obtained by the trigonometric interpolation method. The cumulative area bounded by the hysteresis loops was later calculated using the Trapezoidal method (see Green et al. 2000). Figure 2 illustrates the EDE values, which are calculated from various depth hysteresis loops of downhole array sites. The EDE value increases with increasing depth for all downhole array sites, as seen in Figure 2. This behavior is justifiable as more energy is needed to overcome the friction between soil particles under higher confining pressure; EDE is therefore, expected to increase with increasing depth.

Figure 1.

Shear stress-strain hysteresis loops at a depth of 3.95 m for the Corona site during the 2008 Chino Hills earthquake.

Figure 2.

Earthquake demand energy (EDE) versus depth for the 18 downhole arrays.

Estimation of Shear Stress Time History

Seed and Idriss (1971) proposed a simplified procedure to assess the induced cyclic shear stress at depth (see Figure 3). They assumed a rigid behavior for the soil column above an element of soil at depth h, and obtained the maximum shear stress on a horizontal plane by the following equation:

τ_{\max (r i g i d)} = ρ h a_{\max}

(3)

where ρ is the mass density and a_max is the surface peak acceleration. Since a soil column is not a rigid body, they introduced a nonlinear shear mass participation factor, r_d, which is the ratio of maximum shear stress of deformable soil, to the rigid body. As shown in the following equation, Cetin and Seed (2004) used the results of 2,153 seismic site response analyses to determine r_d based on site conditions and ground motion excitation characteristics for depths of less than 20 m.

r_{d} = \frac{[1 + \frac{- 23.013 - 2.949 a_{\max} + 0.999 M_{w} + 0.0525 V_{S 12}}{16.258 + 0.201 e^{0.341 (- h + 0.0785 V_{S 12} + 7.586}}]}{[1 + \frac{- 23.013 - 2.949 a_{\max} + 0.999 M_{w} + 0.0525 V_{S 12}}{16.258 + 0.201 e^{0.341 (0.0785 V_{S 12} + 7.586}}]} \pm σ_{ε r_{d}}

(4)

σ_{ε r_{d}} = {\begin{matrix} 0.0 198 h^{0.850} h < 12 m \\ 0.0 198 12^{0.850} h \geq 12 m \end{matrix}

(5)

where a_max is the maximum surface ground acceleration (g), M_w is the earthquake magnitude, V_S12 is the average shear wave velocity of the top 12 m of the site (m/s), h is the considered depth (m), and σ_εrd is the standard deviation value.

Figure 3.

A schematic illustration for determining maximum shear stress, τ_max, at depth and stress reduction coefficient, r_d (after Seed and Idriss 1971).

In the present study, the authors employed Equation 4 and replaced the maximum surface ground acceleration, a_max, with the acceleration time history, a_(t). This modified coefficient was referred to as the nonlinear time series factor, r_d*(t), and was used to estimate the shear stress time history at a depth of h, based on Equation 3:

τ_{(h, t)} = ρ h a_{s (t)} r_{d (t)}^{*}

(6)

where τ_(h,t) is the shear stress time history at depth h, and a_s(t) is the surface ground acceleration time history. Figure 4a illustrates the calculation using the measured vertical array records, and predicted the shear stress time histories of the Eureka site, at a depth of 9.5 m and during the 2010 Ferndale Area earthquake, by the simplified method. Moreover, the cross-correlation function, a criterion for measuring the similarity between two signals (Bendat and Piersol 1980), was used to assess the accuracy of this estimation. As shown in Figure 4b, the cross-correlation function reaches a major peak of 0.98 at a time delay of 0.00 s. It indicates a 98% similarity between the two shear stress time histories in Figure 4a, without any time delay. In addition, the obtained cross-correlation number for all downhole array sites ranged between 96% and 99%, suggesting that the estimation of shear stress time histories at shallow depths (2.5 m to 16 m depth) is accurate. It should be noted that the high accuracy of these estimation is due not only to the accuracy of Equation 6, but also to the non-significant changes in the acceleration time histories at shallow depths for these sites. In other words, the proposed nonlinear interpolation method determines the shear stress time history, using the average surface and downhole acceleration time histories. It is therefore, expected that the results of Equation 6, which were obtained using only the surface acceleration time history, to be almost similar to the results of the interpolation method that used two acceleration records.

Figure 4.

Comparison between calculated, using measured vertical array records, and predicted shear stress time histories at a depth of 9.5 m for the Eureka site during the 2010 Ferndale Area earthquake: (a) Shear stress, (b) Cross-correlation between calculated and predicted time series.

Estimation of Shear Strain Time History

In this study, it was observed that the time histories of displacement differences between the surface and shallow depths of the sites in all downhole array data were greatly related to the combination of the smoothed acceleration time history and the velocity time history at ground surface. Figure 5 illustrates the differences between the displacement time histories at the ground surface and at a depth of 21.6 m for the Corona downhole array site during the 2008 Chino Hills earthquake. In addition, this figure shows the surface acceleration, smoothed using the LOESS method, and the surface velocity time histories of the site during the earthquake. As seen in this figure, the waveforms of these time series are largely similar, which was almost true for all the events in Table 1. The following relation was thus suggested to estimate differences in the displacement of time histories:

Δ d_{(t)} = K_{1} a_{s m o o t h (t)} + K_{2} v_{(t)}

(7)

where Δd_(t) is the displacement difference time history between the surface and downhole accelerometer levels, a_smooth(t) is the smoothed acceleration time history at the ground surface, v_(t) is the surface velocity time history, and K₁ and K₂ are the empirical factors which should be calculated. In the present study, a relationship for estimating the maximum displacement differences was generated using a nonlinear regression analysis of all the downhole array data, as given below:

Δ d_{(\max)} = 0.332 a_{s (\max)} . r_{d} {(\frac{h}{V_{S 30}})}^{2} + 0.00 {484 h (d}_{s (max)})^{(\frac{1 . 748 – 6 .0 13}{\ln V_{S 30}})}

(8)

where a_s(max) is the maximum surface acceleration (g), V_S30 is the average shear wave velocity of the top 30 meters of soil layer (m/s), d_s(max) is the maximum displacement on the ground surface (m), h is depth (m), and r_d is obtained using Equation 4. In the second term of Equation 8, the interaction between the parameters can compensate the effect of soil nonlinear behavior, which was ignored as a simplifying assumption. The following relations yield the amount of factors K₁ and K₂, so that the maximum points of the real time histories coincide with the estimated time histories:

K_{1} = {\begin{cases} \begin{matrix} 0.003 + 0.776 \frac{Δ d_{(\max)}}{a_{s (\max)}} & a_{s (\max)} < 0.1 g \end{matrix} \\ \begin{matrix} 0.905 \frac{Δ d_{(\max)}}{a_{s (\max)}} & a_{s (\max)} \geq 0.1 g \end{matrix} \end{cases}

(9)

K_{2} = 23.11 \frac{Δ d_{(\max)}}{h^{1.25}}

(10)

where the dimensions of the parameters are the same as Equation 8. By substituting Δd_(t) (calculated using Equation 7) in the trigonometric interpolation method (Davis and Berrill 1998) or linear interpolation method (Zeghal and Elgamal 1993), the shear strain time history can be made available.

Figure 5.

Comparison between time histories of displacement difference, surface acceleration, and surface velocity of the Corona site during the 2008 Chino Hills earthquake: (a) displacement difference between the surface and the depth of 21.6 m, (b) smoothed surface acceleration, (c) surface velocity.

As evident in Figures 6a and 6b, Equations 7–10 are well suitable for estimating the shear strain time history. Figure 7 illustrates the predicted EDE using the proposed simplified method versus the calculated values, for 18 vertical arrays. In addition, Figure 8 demonstrates the degree of coincidence between the predicted and calculated results of 72 points (four different depths of 4, 7, 10, and 14 m for 18 vertical arrays).

Figure 6.

Comparison between calculated, using measured vertical array records, and predicted shear strain time histories at a depth of 9.5 m for the Eureka site during the 2010 Ferndale Area earthquake: (a) Shear strain, (b) cross-correlation between calculated and predicted time series.

Figure 7.

Calculated EDE, using measured vertical array records (presented in Figure 1), and predicted EDE, by simplified method, for 18 vertical arrays.

Figure 8.

Predicted EDE, using simplified method (with only surface acceleration), against calculated EDE for different depths of 18 vertical arrays (with accelerations in depth and surface).

Ground Motion Database

In this study, 328 records of ground motion events from 27 shallow crustal earthquakes were utilized. These data were extracted from the Pacific Earthquake Engineering Research Center (PEER), and recorded in free-field stations, mainly in California, Japan, New Zealand, Taiwan, Iran, Italy, and Turkey. In addition, the database was collected so that the records for each faulting mechanism (reverse, normal, and strike slip) cover uniformly a wide range of parameters such as magnitude, distance, shear wave velocity, and PGA. Furthermore, these ground motion data should be recorded before soil liquefaction occurs during the same motion. Table 2 presents relevant information pertaining the earthquakes. The three types of data, categorized according to fault mechanism, were normal events (52 records), reverse events (154 records), and strike-slip events (122 records). According to the 1997 UBC (Dobry et al. 2000), the sites corresponding to these data have been classified into class C, D, and E. Figure 9 shows the distribution of earthquake data with respect to moment magnitude and PGA against the closest distance to the rupture plane.

Figure 9.

Distribution of earthquake data with respect to magnitude (left) and PGA (right).

Table 2

List of earthquakes used in development of attenuation model (data source: PEER)

No.	Event	Year	Faulting mechanism	M _w	Distance range (km)	Number of data
1	San Fernando	1971	Reverse	6.61	19.3–92.59	16
2	Tabas (Iran)	1978	Reverse	7.35	2.05–120.81	5
3	Mammoth Lakes	1980	Strike Slip	4.73	5.17	1
4	Irpinia (Italy)	1980	Normal	6.2	8.83–122.35	10
5	Borah Peak	1983	Normal	5.1	18.51–22.35	2
6	Morgan Hill	1984	Strike Slip	6.19	9.81–45.47	5
7	Lazio-Abruzzo (Italy)	1984	Normal	5.8	31.23–51.29	5
8	SMART1 (Taiwan)	1985	Normal	5.8	41.71–41.82	8
9	Kalamata	1986	Normal	5.3	5.6–28.86	6
10	Whittier Narrows	1987	Reverse	5.99	13.28–64.27	32
11	Manjil (Iran)	1990	Strike Slip	7.37	12.55–93.62	5
12	Landers	1992	Strike Slip	7.28	2.19–148.47	32
13	Little Skull Mtn	1992	Normal	5.65	16.06–99.58	4
14	Northridge	1994	Reverse	6.69	12.44–68.62	23
15	Kobe	1995	Strike Slip	6.9	1.47–95.72	19
16	Dinar (Turkey)	1995	Normal	6.4	3.36–86.31	4
17	Northwest (China)	1997	Normal	5.93	27.86–46.24	4
18	Chi-Chi	1999	Reverse	7.62	3.12–89.48	27
19	Chi-Chi 03	1999	Reverse	6.2	8.48–105.48	28
20	Chi-Chi 04	1999	Strike Slip	6.2	6.2–124.74	21
21	Tottori	2000	Strike Slip	6.61	5.86–49.82	23
22	Anza	2001	Normal	4.92	20.04–120.14	5
23	Niigata	2004	Reverse	6.63	6.63–22.74	2
24	Iwate	2008	Reverse	6.9	22.75–25.16	3
25	L'Aquila (Italy)	2009	Normal	6.3	73.34–148.97	4
26	Darfield	2010	Strike Slip	7	1.22–95.18	16
27	Christchurch	2011	Reverse	6.2	3.26–91.32	18

Development of Attenuation Model

The proposed method was employed to estimate the amount of EDE at different depths for the 328 ground motion records. Because the 10 m-interval (4–14 m) is usually the most important depth for a liquefaction study, the depths selected for this study were 4, 7, 10, and 14 m.

Effective Parameters

In this study, moment magnitude (M_w) was used to define the earthquake magnitude. Figure 10 illustrates the calculated EDE of several sites at a depth of 10 m for the earthquakes with different magnitudes. These sites have been located at approximately the same distance from the rupture plane, and all the earthquakes were generated from the reverse fault type. As shown in this figure, the EDE for all the soil site classes increases as the moment magnitude of earthquakes increases.

Figure 10.

Variation in earthquake demand energy (EDE) with earthquake magnitude, M_w, for different site classes.

In the present study, the closest distance to the rupture plane (R_rup) was selected as the site-to-source distance. Figure 11 depicts the calculated EDE at a depth of 10 m for the sites with the different site-to-source distances. In this figure, all the sites were subjected to the earthquakes with the same magnitude and the same fault mechanism. This figure indicates that the EDE at any site decreases with increasing distance between the source and the site for all the soil site classes. In addition, Figure 11 shows that V_S30 also has an effect on the amount of EDE. According to this figure, dense soils, with higher shear wave velocities, are able to transmit more seismic energy, and less input energy is dissipated through these soils compared to looser soils.

Figure 11.

Earthquake demand energy (EDE) versus closest distance, R_rup, for different site classes.

The effects of faulting mechanism on the amount of EDE for a constant earthquake magnitude within the same site class are presented in Figure 12. From this figure, it can be understood that reverse mechanism produces and propagates more seismic energy than normal and strike-slip mechanisms. McGarr (1984) used theoretical and empirical studies to conclude that since strike-slip and reverse-faulting are in the compressional stress tectonic regime, they are associated with higher ground motion, rather than normal-faulting earthquakes located in the extensional stress tectonic regime. His findings are consistent with Figure 12 results.

Figure 12.

Earthquake demand energy (EDE) versus closest distance, R_rup, for different faulting mechanisms.

In addition to the parameters discussed above, the PGA was selected as the surface amplitude of seismic waves, which represent subsurface conditions in the attenuation model. Figure 13 illustrates the relationship between EDE and PGA of 23 sites in California and Japan. This figure shows that despite having similar characteristics (M_w, R_rup, V_S30, and fault type), the PGA and the corresponding EDE of the 23 sites have different values. Thus, since the PGA is affected by the site condition, it can also be a significant parameter for estimating EDE at a site.

Figure 13.

Earthquake demand energy (EDE) versus PGA for 23 sites, related to the 1971 San Fernando and 1994 Northridge earthquakes in California, and the 2004 Niigata, 2007 Chuetsu-oki, and 2008 Iwate earthquakes in Japan.

Proposed Attenuation Model

A simple attenuation model was proposed to estimate earthquake demand energy using regression analysis of 1,312 data points (four different depths for 328 earthquake records). This model includes three main functions and the model error, ε_log(EDE), as follows:

\log (E D E) = F_{1} (M_{w}, F) + F_{2} (R_{r u p}, R_{t h}) + F_{3} (V_{S 30}, P G A, z) + ε_{\log (E D E)}

(11)

where EDE is the earthquake demand energy at the site (J/m³). Also, F₁, F₂, and F₃ are the source, the distance, and the site effect functions, respectively (see below).

Source function:

F_{1} (M_{w}, F) = F * c_{1} M_{w}^{c_{2}}

(12)

where M_w is the earthquake magnitude and F is the faulting mechanism factor, obtained at 1.00 for normal faulting, 1.12 for strike-slip faulting, and 1.25 for reverse faulting using regression analysis.

Distance function:

F_{2} (R_{r u p}, R_{t h}) = {\begin{cases} \begin{matrix} c_{3} Ln {(R_{rup})}^{c_{4}} & R_{r u p} < R_{t h} \end{matrix} \\ \begin{matrix} c_{3} Ln {(R_{th})}^{c_{4}} + c_{5} {(\frac{R_{r u p} - R_{t h}}{R_{t h}})}^{c_{6}} & R_{r u p} \geq R_{t h} \end{matrix} \end{cases}

(13)

where R_rup is the closest distance to the rupture plane (km) and R_th is the threshold distance (km). To determine the threshold distance, R_th was assumed proportional to the square of the magnitude, and the calculated EDE at a depth of 10 m was plotted against K (ratio of R_rup and

M_{w}^{2}

) in the arithmetic space (Figure 14). According to this figure, K = 0.3 was assumed as an approximate boundary between two regions, which seem different considering the rate of EDE. Therefore, the following empirical equation was proposed for deriving R_th value:

R_{t h} = 0.3 M_{w}^{2}

(14)

Site effect functions:

F_{3} (V_{S 30}, P G A, z) = c_{7} L n (V_{S 30}) + c_{8} P G A^{c_{9}} + c_{10} L n (z)

(15)

where z is the depth under consideration at the site (in meters), V_S30 is the average shear wave velocity of the top 30 m of the site (m/s), and PGA is the peak ground surface acceleration (g).

Figure 14.

Earthquake demand energy (EDE) versus an empirical parameter, K (ratio of closest distance (km) and magnitude-squared).

The coefficients c₁ to c₁₀ are defined for three functions as shown in Table 3. The regression procedure for determining these coefficients included three steps. In the first step, the data points were categorized by fault mechanism, and the logarithm of EDE was estimated by incorporating source and distance functions for any mechanism. Then, the coefficient c₁, which was obtained using normal mechanism data, was selected for other mechanisms, and the coefficient c₂ and the faulting mechanism factor were determined for each of the three faulting mechanisms. In this step, the effect of near-field ground motions was ignored. In the second step, c₁ and c₂ were fixed and the coefficients of the near and far distance functions were estimated. In the third step, all coefficients, except the coefficient of the site effect function were fixed, and the amounts of c₇ to c₁₀ were estimated.

Table 3

Regression coefficients and faulting mechanism factors for attenuation model

Model coefficients										Faulting mechanism factor (F)
c ₁	c ₂	c ₃	c ₄	c ₅	c ₆	c ₇	c ₈	c ₉	c ₁₀	Normal	Strike-Slip	Reverse
0.171	1.791	-0.389	0.125	-0.019	2.119	-1.021	2.652	0.101	0.582	1.00	1.12	1.25

Since the input energy values, used in the regression model, are associated with some degree of uncertainty, the uncertainty propagation in the model error must be considered. The standard deviation of the proposed simplified method was estimated using results of the predicted energy values of 18 vertical arrays (see Figure 8). Because the standard deviation, 0.243, is completely uncorrelated with the independent variables of the model, it cannot influence the model parameters, and only affects the standard deviation of the regression model, which is calculated to be 0.635. Eventually, the total model error is obtained from the square root of the sum of variances:

ε_{L o g}_{(EDE)} = \sqrt{{0.243}^{2} + {0.635}^{2}} = 0.680

(16)

where this error is conditioned on the knowledge of independent variables. Figure 15 illustrates the residual values imposed by the estimation of the total data points, which include the function of magnitude, closest distance, V_S30, and PGA. As shown in this figure, no significant relation is observed between residual and estimator parameters.

Figure 15.

Distribution of estimation residuals with respect to magnitude, closest distance, V_S30, and PGA.

Model Predictions

To assess the accuracy of the proposed attenuation model, two liquefied downhole array sites were selected (see Table 4), and their computed capacity energies were compared with the predicted demand energies (when the recorded accelerations from these two sites were not included in the databank, and implemented for the model).

Table 4

Information of two liquefied downhole array sites used for validation of the EDE model

Station name	Event name	Year	Magnitude (M_w)	Mechanism	R_rup (km)	V_S30 (m/s)	PGA (g)
Port Island site	Hyogoken-Nambu	1995	6.9	Strike-Slip	3.31	198	0.341
Wildlife Liquefaction	Elmore Ranch	1987	6.22	Strike-Slip	17.59	179	0.13
Array site	Superstition Hills	1987	6.54	Strike-Slip	23.85	179	0.21

Port Island Site, Japan

The Port Island (PI) is a reclaimed island on the South West coast of Kobe, Japan. In August 1991, several instruments including accelerometers at the ground surface and at the depths of 16, 32, and 83 m were installed to record future earthquakes (Iwasaki 1995).

On January 17, 1995, the PI site was severely shaken by the Hyogoken-Nambu earthquake, M_w = 6.9, and shallow layers were liquefied during this earthquake. The surface evidences of liquefaction such as ground settlement, lateral spreading, and ground fissures were observed after the earthquake (O'Rourke 1995). The recorded acceleration of the earthquake site at four different depths are illustrated in Figure 16. This figure also illustrates the computed cumulative dissipated energies at the depths of 4, 7, 10, and 14 m, using a nonlinear interpolation method. Zeghal and Elgamal (1996) characterized these recorded accelerations in three stages: (1) an initial phase of low amplitude vibration, (2) a strong shaking phase, and (3) low amplitude vibration at the surface. Using finite element analyses, they showed that an abrupt rise in excess pore pressure occurred mainly during stage 2.

Figure 16.

Horizontal (N-S) accelerations of Hyogoken-Nambu earthquake recorded at different depths of the PI site (Zeghal and Elgamal 1996), and computed cumulative dissipated energy at different depths during the earthquake.

In this study, it was assumed that the liquefaction triggering point falls into stage 3, and that the amount of cumulative dissipated energy at 7 s and 40 s were considered for the lower and the upper bounds of capacity energy for this soil layer, respectively. Figure 17 illustrates the capacity energy bounds, calculated using the nonlinear interpolation method, along the top 14 m of the site. As seen in Figure 17, the profile of the predicted demand energy falls on the right side of the upper bound, indicating liquefaction occurrence at this site. It should be mentioned that the closest distance of this site (3.31 km), which was considered a near-field site, was smaller than the threshold distance of the Hyogoken-Nambu earthquake (14.28 km).

Figure 17.

Calculated capacity energy region and predicted demand energy for Hyogoken-Nambu earthquake along the top 14 m of the PI site.

Wildlife Liquefaction Array Site, California

The Wildlife Liquefaction Array (WLA) site was established in 1982 on the plains of the Alamo River in the Imperial Valley desert of Southern California (Holzer and Youd 2007). The United States Geological Survey (USGS) instrumented this site with two triaxial accelerometers and six piezometers (Youd and Wieczorek 1984).

On 24 November 1987, two significant earthquakes occurred along the southern San Jacinto fault zone and the related structural elements in Southern California (Hanks and Allen 1989). The first event was the Elmore Ranch earthquake, M = 6.2, and was followed by the Superstition Hills earthquake, M = 6.6, which occurred approximately 12 hours later. Holzer et al. (1989) reported the occurrence of liquefaction at the Wildlife site during the Superstition Hills earthquake in the sandy silt and silty sand layers (2.5–6.8 m). In addition, other researchers (e.g., Dobry et al. 1989, Hushmand et al. 1992, Matasovic and Vucetic 1994) conducted various studies on the Wildlife site after the 1987 earthquakes, and reached a consensus that while the Superstition Hills earthquake did cause liquefaction at the Wildlife site, the Elmore Ranch earthquake did not cause liquefaction. Figure 18 depicts the recorded downhole array accelerations, excess pore pressure ratio, and the cumulative dissipated energy of the WLA site during the Superstition Hills earthquake. Zeghal and Elgamal (1996) characterized these recorded downhole array accelerations in four stages (stage 1: 0–13.7 s; stage 2: 13.7–20.6 s; stage 3: 20.6–40 s; and stage 4: 40–97 s). Stage 4 was selected to calculate the possible amounts of capacity energy at the WLA site. Figure 19 illustrates the lower and upper bounds of the capacity energy area and the energy demand amount of the Elmore Ranch and Superstition Hills events. These demand energies were predicted using the EDE attenuation model along the top 6.6 m of the WLA site. Since the threshold distances of the Elmore Ranch and Superstition Hills earthquakes (11.61 and 12.83 km, respectively) were smaller than the site-to-source distances (17.59 and 23.85 km, respectively), the motions were not considered near-field.

Figure 18.

Horizontal (N-S) accelerations and pore pressure ratio at 2.9 m (Zeghal and Elgamal 1996), besides cumulative dissipated energy at different depths of the WLA site during the Superstition Hills earthquake.

Figure 19.

Calculated capacity energy region and predicted demand energies for the Elmore Ranch and Superstition Hills earthquakes along the top 6.6 m of the WLA site.

Using a wide-range of experimental database, Baziar et al. (2011b), proposed an equation to calculate the amount of capacity energy by soil parameters, including initial mean effective confining pressure, relative density, fine content, uniformity coefficient, and mean grain size. They calculated the capacity energy of the liquefied layer of the WLA site using their equation (see Figure 19). According to the findings, namely the comparison of the two predicted demand energies with the calculated capacity energy region and the empirical study of Baziar et al. (2011b), it can be concluded that this site did indeed liquefy during the Superstition Hills earthquake, but did not liquefy during the Elmore Ranch earthquake.

Conclusions

In this research, an attenuation model was developed to estimate earthquake demand energy (EDE) at a site. This research was conducted in two parts: (1) providing a new method for determining shear stress and shear strain time histories at shallow depth of a site when only one surface accelerometer is available; (2) using the obtained method to calculate the demand energy for the wide-ranging ground motion events, and creating a suitable database for developing a rigorous relationship among EDE, earthquake parameters, and soil characteristics.

Initially, recorded data from 18 downhole arrays in California were analyzed, and several semi-empirical equations for estimating shear stress and shear strain time histories due to earthquake loading at the shallow depth of a site were generated. The demand energies of 328 worldwide earthquake records at various depths of the free-field sites (4, 7, 10, and 14 m) were then calculated using the new method. Subsequently, an attenuation model was proposed using a nonlinear regression analysis. The obtained EDE values were based on the source mechanism, the distance, and the site effect functions. Finally, the proposed attenuation model was used to predict the demand energies at two well-known liquefied sites.

The developed method to calculate EDE using only one recorded acceleration at the surface sites, referred to in the first part of this study, can largely remove the limitation of the required two recorded surface and downhole accelerations. Also, analyzing downhole array records showed that the demand energy increases with increasing depth at 18 downhole arrays. The proposed attenuation model takes into consideration all the following, namely the role of source mechanism, distance and site effects. In the source function, the faulting mechanism factor (F) was considered as a constant factor, and the amount of F was determined using regression analysis. In the distance function, a simple energy-based criterion was used to consider the effects of near-field area on the EDE value. In the site effects function, besides the V_S30, the PGA was used as representative of subsurface characteristics. The presented attenuation model is suitable for predicting the EDE at the top 15 m of free-field sites and earthquakes with M_w > 4.5 and R_rup < 150 km.

Footnotes

Acknowledgments

The downhole array data were provided by the Center for Engineering Strong Motion Data (CESMD). Furthermore, the worldwide ground motion data was made available by the Pacific Earthquake Engineering Research Center (PEER). In addition, the Wildlife Liquefaction Array (WLA) data were obtained through the Network for Engineering Earthquake Simulation (NEES) at the University of California Santa Barbara. The authors also thank Professor I. Towhata, University of Tokyo, for providing data of the Port Island (PI) site in Japan. Useful comments provided by Dr. H. Shahnazari and Dr. H. Sharafi are also gratefully acknowledged.

References

Abrahamson

N. A.

, and Shedlock

K. M.

, 1997. Overview, Seismol. Res. Lett 68, 9–23.

Abrahamson

N. A.

, Atkinson

, Boore

, Bozorgnia

, Campbell

, Chiou

, Idriss

I. M.

, Silva

, and Youngs

, 2008. Comparisons of the NGA ground-motion relations, Earthquake Spectra 24, 45–66.

Baziar

M. H.

, and Jafarian

, 2007. Assessment of liquefaction triggering using strain energy concept and ANN model, capacity energy, Soil Dynamics and Earthquake Engineering 27, 1056–1072.

Baziar

M. H.

, and Sharafi

, 2011. Assessment of silty sand liquefaction potential using hollow torsional test: an energy approach, Soil Dynamics and Earthquake Engineering 31, 857–865.

Baziar

M. H.

, Shahnazari

, and Sharafi

, 2011a. A laboratory study on the pore pressure generation model for Firouzkooh silty sands using hollow torsional test, International Journal of Civil Engineering 9, 126–134.

Baziar

M. H.

, Jafarian

, Shahnazari

, Movahed

, and Tutunchian

M. A.

, 2011b. Prediction of strain energy-based liquefaction resistance of sand–silt mixtures: an evolutionary approach, Computers and Geosciences 37, 1883–1893.

Bendat

J. S.

, and Piersol

A. G.

, 1980. Engineering Applications of Correlation and Spectral Analysis, Wiley, New York.

Berrill

J. B.

, and Davis

R. O.

, 1985. Energy dissipation and seismic liquefaction of sands: revised model, Soils Found. 25, 106–118.

Cetin

K. O.

, and Seed

R. B.

, 2004. Nonlinear shear mass participation factor r_d for cyclic shear stress ratio evaluation, Soil Dynamics and Earthquake Engineering 24, 103–113.

10.

Davis

R. O.

, and Berrill

J. B.

, 1982. Energy dissipation and seismic liquefaction in sands, Earthquake Engineering and Structural Dynamics 10, 59–68.

11.

Davis

R. O.

, and Berrill

J. B.

, 1998. Rational approximation of shear stress and shear strain based on downhole acceleration measurements, International Journal of Numerical Analysis Methods and Geomechanics 22, 603–619.

12.

Dobry

, Elgamal

A-W.

, Baziar

M. H.

, and Vucetic

, 1989. Pore pressure and acceleration response of Wildlife site during the 1987 earthquake, in The Second US–Japan Workshop on Soil Liquefaction ( O'Rourke

T. D.

and Hamada

, editors), 145–160.

13.

Dobry

, Borcherdt

R. D.

, Crouse

C. B.

, Idriss

I. M.

, Joyner

W. B.

, Martin

G. R.

, Power

M. S.

, Rinne

E. E.

, and Seed

R. B.

, 2000. New site coefficients and site classification system used in recent building seismic code provisions, Earthquake Spectra 16, 41–67.

14.

Figueroa

J. L.

, Saada

A. S.

, Liang

, and Dahisaria

M. N.

, 1994. Evaluation of soil liquefaction by energy principles, J Geotech Eng, ASCE, 120, 1554–1569.

15.

Green

R. A.

, Mitchell

J. K.

, and Polito

C. P.

, 2000. Energy-based excess pore pressure generation model for cohesionless soils, Proceedings of the John Booker Memorial Symposium, Sydney, NSW, Australia.

16.

Gregor

, Abrahamson

N. A.

, Atkinson

G. M.

, Boore

D. M.

, Bozorgnia

, Campbell

K. W.

, Chiou

, Idriss

I. M.

, Kamai

, Seyhan

, Silva

, Stewart

J. P.

, and Youngs

, 2014. Comparison of NGA-West2 GMPEs, Earthquake Spectra 30, 1179–1197.

17.

Gutenberg

, and Richter

C. F.

, 1954. Seismicity of the Earth, Princeton University Press.

18.

Gutenberg

, and Richter

C. F.

, 1956. Magnitude and energy of earthquakes, Annali de Geofisica 9, 1–15.

19.

Hanks

T. C.

, and Allen

C. R.

, 1989. The Elmore Ranch and Superstition Hills earthquakes of 24 November 1987: Introduction to the special issue., Bull. Seismol. Soc. Am. 79, 231–239.

20.

Holzer

T. L.

, Youd

T. L.

, and Hanks

T. C.

, 1989. Dynamics of liquefaction during the Superstition Hills Earthquake (M = 6.5) of November 24, 1987, Science 244, 56–60.

21.

Holzer

T. L.

, and Youd

T. L.

, 2007. Liquefaction, ground oscillation, and soil deformation at the Wildlife array, California, Bulletin of the Seismological Society of America 97, 961–976.

22.

Hushmand

, Scott

R. F.

, and Crouse

C. B.

, 1992. In-place calibration of USGS transducers at Wildlife liquefaction site, California, Proc. 10^th World Conf. on Earthquake Engrg., A.A. Balkema Publishers, Rotterdam, Netherlands, 1263–1266.

23.

Iwasaki

, 1995. Geological and geotechnical characteristics of Kobe area strong ground motion records by 1995 Kobe earthquake, Tsuchi-to-Kiso, Japanese Society of Soil Mechanics and Foundation Engineering 43, 15–20.

24.

Jafarian

, Sadeghi Abdollahi

, Vakili

, Baziar

M. H.

, and Noorzad

, 2011. On the efficiency and predictability of strain energy for the evaluation of liquefaction potential: a numerical study, Computers and Geotechnics 38, 800–808.

25.

Jafarian

, Towhata

, Baziar

M. H.

, Noorzad

, and Bahmanpour

, 2012. Strain energy based evaluation of liquefaction and residual pore water pressure in sands using cyclic torsional shear experiments, Soil Dyn. Earthquake Eng. 35, 13–28.

26.

Jafarian

, Vakili

, Sadeghi Abdollahi

, and Baziar

M. H.

, 2014. Simplified Soil Liquefaction Assessment Based on Cumulative Kinetic Energy Density: Attenuation Law and Probabilistic Analysis, International Journal of Geomechanics 14, 267–281.

27.

Kokusho

, and Suzuki

, 2011. Energy flow in shallow depth based on vertical array records during recent strong earthquakes, Soil Dynamics and Earthquake Engineering 31, 1540–1550.

28.

Kokusho

, and Suzuki

, 2012. Energy flow in shallow depth based on vertical array records during recent strong earthquakes (supplement), Soil Dynamics and Earthquake Engineering 42, 128–131.

29.

Kokusho

, 2013. Liquefaction potential evaluations: energy-based method versus stress-based method, Can Geotech J 50, 1–12.

30.

Law

K. T.

, Cao

Y. L.

, and He

G. N.

, 1990. An energy approach for assessing seismic liquefaction potential, Can Geotech J 27, 320–329.

31.

Matasovic

, and Vucetic

, 1994. Analysis of seismic records obtained on November 24, 1987 at the Wildlife Liquefaction Array, (same as Matasovic J, M. Sc. Thesis, Civil Engineering Department, University of California, Los Angeles, California 1993).

32.

McGarr

, 1984. Scaling of ground motion parameters, state of stress, and focal depth, Journal of Geophysical Research 89, 6969–6979.

33.

Nemat-Nasser

, and Shokooh

, 1979. A unified approach to densification and liquefaction of coheshionless sand in cyclic shearing, Can Geotech J 16, 659–678.

34.

O'Rourke

T. D.

, 1995. Geotechnical Effects, Preliminary report from the Hyogoken-Nambu earthquake of January 17, 1995, National Center for Earthquake Engineering Research Bulletin, SUNY, Buffalo.

35.

Richter

C. F.

, 1936. An instrumental earthquake magnitude scale, Bulletin, the Seismological Society of America 25, 1–32.

36.

Seed

H. B.

, and Idriss

I. M.

, 1971. Simplified procedure for evaluating soil liquefaction potential, J Soil Mech. Found. Div., ASCE, 97(SM8), 1249–1274.

37.

Towhata

, and Ishihara

, 1985. Shear work and pore water pressure in undrained shear, Soils and Foundations 25, 73–84.

38.

Towhata

, 2008. Geotechnical Earthquake Engineering, Springer Series in Geomechanics and Geoengineering, Springer-Verlag, Berlin, 684 pp.

39.

Travasarou

, Bray

J. D.

, and Abrahamson

N. A.

, 2003. Empirical attenuation relationship for Arias intensity, Earthquake Engng Struct. Dyn. 32, 1133–1155.

40.

Trifunac

M. D.

, 1976. Preliminary analysis of the peaks of strong earthquake ground motion dependence of peaks on earthquake magnitude, epicentral distance, and recording site conditions, Bull. Seismol. Soc. Am. 66, 189–219.

41.

Trifunac

M. D.

, 1995. Empirical criteria for liquefaction in sands via standard penetration tests and seismic wave energy, J Soil Dynamic and Earthquake Engineering 14, 419–426.

42.

Youd

T. L.

, and Wieczorek

G. F.

, 1984. Liquefaction during 1981 and Previous Earthquakes Near Westmorland California, USGS Open-file Rep. 84-680, U.S. Geological Survey, Washington, D.C.

43.

Zeghal

, and Elgamal

A-W.

, 1993. Lotung Site: Downhole Seismic Data Analysis, Report, Department of Civil Engineering Rensselaer Polytechnic Institute (to appear as Electric Power Research Institute Report TR-106093).

44.

Zeghal

, and Elgamal

A-W.

, 1996. Analysis of site liquefaction using downhole array seismic records, 11th World Conference on Earthquake Engineering, Paper No. 371, Acapulco, Mexico.