A new method to construct model structures from left Frobenius pairs in extriangulated categories

1.   Introduction

In the 1960s the great impact of earthquake-induced liquefaction on the sustainability of a modern society was first recognized when the Alaskan and Niigata earthquakes incurred severe damage to urban infrastructures and buildings constructed on poor soils of newly urbanized lowlands. Since then, liquefaction damage has occurred often in many countries around the world as one of the serious impacts on the sustainability and resilience of society. As populations and economies are growing rapidly in developing countries, urbanization in lowlands of poor soil conditions will make the liquefaction threats more serious because of the rising groundwater table due to surging sea levels deemed to be inevitable in global warming. Thus, the liquefaction mitigation of urban facilities is increasingly prioritized in engineering designs not only in developed but also in developing countries.

Moreover, from the dawn of the liquefaction research, a stress-based method (SBM), wherein undrained cyclic strength as liquefaction capacity compares with seismically induced shear stress as liquefaction demand, has been developed (Seed & Idriss 1971 ^[1]), and standardized in engineering practice for liquefaction potential evaluations in many design codes worldwide. In its simplified procedure wherein the earthquake effect is idealized by harmonic stress of a given amplitude and given cycles, the cyclic resistance ratio (CRR) is determined by in situ tests such as Standard Penetration Tests (SPTs) or Cone Penetration Tests (CPTs) using their depth-dependent penetration resistances plus pertinent soil parameters (fines content F_c and plasticity index I_p), while the earthquake-induced cyclic stress ratio (CSR) is calculated from peak ground horizontal acceleration (PGA) assuming depth-dependent stress decay curves or SH-wave propagation analyses using design acceleration motions. However, the liquefaction behavior is not governed solely by applied stress but also by earthquake-induced shear strain and associated soil dilatancy which builds up the pore pressure in undrained saturated soils. Hence, not only earthquake-induced stresses but also other pertinent parameters governing strains are needed to make relevant liquefaction evaluations.

Apart from that, a possibility of a unified liquefaction theory was first presented (Nemat-Nasser and Shokooh, 1979 ^[2]), which paved the way to develop the energy-based method for liquefaction evaluation (EBM). Instead of the SBM, it focuses on dissipated energy determined from the stress versus strain relationship causing dilatancy and pore-pressure buildup for liquefaction during cyclic loading. Since then, quite many experimental efforts have been made to confirm and demonstrate the superiorities of the EBM compared to the SBM.

As obvious merit of the EBM, Figure 1 exemplifies undrained torsional simple shear test results on medium loose clean sand (D_r ≈ 50%) cyclically loaded by two widely varied earthquake motions: (a) The 2011 Tohoku earthquake (M_J = 9.0: M_J = Japanese Meteorological Agency Magnitude similar to the surface-wave magnitude M_s) lasting more than 4 minutes, and (b) 1995 Kobe earthquake (M_J = 7.2) lasting only 20 seconds. This indicates that the cumulative dissipated energy normalized by initial effective confining stress ∑ΔW/σ_c^' in the horizontal axis can uniquely predict the pore-pressure buildup ratio r_u = Δu/σ_c^' despite the tremendous difference between the two motions in amplitudes, time durations, and waveforms. This implies that whatever earthquake motions are in terms of intensities, focal distances, durations, waveforms, predominant frequencies, irregularities, etc., the liquefaction behavior in terms of pressure buildup is determined exclusively by the cumulative energy dissipated in soils ∑ΔW during earthquake motions.

The SBM, in contrast, has to properly choose pertinent parameters critical to make a reliable prediction. Despite such advantages of the EBM, it has rarely been used in engineering practice partly because the SBM has long been authorized in design manuals as a standardized evaluation procedure in the historical background of the force-equilibrium using acceleration in earthquake engineering in general.

However, with the maturity of the EBM research and the increasing need for more reasonable and robust liquefaction evaluation tools, the energy-based liquefaction design method should be authorized so that it can be followed widely by practicing engineers to prepare for increasing liquefaction threats to resilient and sustainable societies.

In this review article, major research developments of the energy-based liquefaction evaluation method (EBM) are summarized by revisiting major previous research. Then, the EBM is discussed concerning the issues.

◆ How to determine capacity energy (dissipated energy) for liquefaction

◆ How to evaluate demand energy (earthquake wave energy) from design earthquake motions or relevant parameters

◆ How to compare the capacity energy with the demand energy to make a reasonable and simplified liquefaction evaluation.

Theoretical backgrounds and experimental data associated with these issues are also addressed to discuss in comparison with the SBM whenever necessary.

Finally, the evaluation procedures of the EBM in engineering practice are delineated step by step based on preceding discussions, which are followed by a couple of case studies to demonstrate its features and differences compared to the SBM.

2.   Previous research on energy-based methods

The energy-based liquefaction evaluation method (EBM) was first proposed by Davis and Berrill (1982) ^[3] and Berrill & Davis (1984) ^[4], followed by Law et al. (1990) ^[5]. In their model, the pore-pressure buildup for liquefaction was correlated with the capacity energy represented by dissipated energy expressed as a function of in situ penetration resistance. The demand energy was seismologically given as a function of earthquake magnitude and hypocenter distance by Gutenberg (1956) ^[6]. It was motivated by the pioneering research on energy capacity for liquefaction by Nemat-Nasser and Shokooh ^[2].

Since then, the uniqueness of dissipated energy or energy capacity in the liquefaction behavior has been experimentally demonstrated by quite many test data in various laboratory cyclic loading tests. Among them, Towhata and Ishihara (1985) ^[7] conducted hollow cylindrical torsional shear tests with various loading paths in which a unique relationship was found between shear work (dissipated energy) and excess pore-pressure buildup to be independent of the shear stress path. Yanagisawa and Sugano (1994) ^[8] conducted similar cyclic shear tests using irregular time histories and found again the uniqueness of cumulative energy on the pressure buildup. Laboratory soil tests were also conducted in a strain-controlled torsional shear test by Figueroa et al. (1994) ^[9], demonstrating the uniqueness of dissipated energy in pore-pressure buildup under different confining stresses. Baziar and Sharafi (2011) ^[10] conducted stress-controlled undrained cyclic torsional tests on silty sands under different confining stresses and found that the dissipated energy for pore-pressure buildup is independent of CSR but highly dependent on the effective confining pressure and fines content. Pan and Yang (2017) ^[11] and Azeiteiro et al. (2017) ^[12] both presented systematic cyclic triaxial test results with irregular time histories to find that the pore-pressure generation and the number of cycles for liquefaction are significantly influenced by the stress amplitudes and their sequence. Conversely, a unique relationship was found between the cumulative dissipated energy and the pore pressure wherein the effect of loading conditions can be ignored, indicating a superiority of the energy-based approach.

Empirical relationships correlating the dissipated energy with the pore-pressure buildup ratio were developed by Green et al. (2000) ^[13] and Jafarian et al. (2012) ^[14] to be used for liquefaction evaluations. Furthermore, the dissipated energy concept has been employed in evaluating liquefaction resistance in practical engineering problems such as plant-rooted soils (Karimzadeh et al. 2021) ^[15] and calcium-carbonate treated sands (Baziar and Alibolandi 2023) ^[16] by developing their energy versus CRR relationships.

Apart from those mainly concerned with pore-pressure buildup, Kazama et al. (2000) ^[17] carried out strain-controlled cyclic triaxial tests to focus on the dissipated energy evaluated in stress ~ strain hysteretic loops even after 100% pressure-buildup to take soil ductility into account in liquefaction design. Kokusho (2013) ^[18] summarized a series of their own cyclic triaxial test results by harmonic loading and found that the cumulative dissipated energy is uniquely correlated with not only pore-pressure buildup but also induced strain not only up to 100% pore-pressure buildup (initial liquefaction) but beyond. Hence, unlike the pressure buildup, the strain can serve as an index to evaluate the severity of liquefaction even after the onset of liquefaction. In torsional simple shear tests on loose sands without/with non-plastic fines, Kokusho and Kaneko (2018) ^[19] confirmed the uniqueness of ΣΔW in determining liquefaction-induced strain not only for harmonic motions but also a variety of earthquake motions.

In some of the research, using reconstituted sand specimens, ΣΔW in EBM and CRR in SBM were found to be uniquely correlated despite wide differences in relative density and fines content in triaxial tests ^[18] and in torsional tests ^[19]. Furthermore, Kokusho and Tanimoto (2021) ^[20] visited a set of cyclic triaxial test results on numerous intact soils sampled from various sites, wherein a unique correlation between the dissipated energy ΣΔW and CRR for the number of cycles for initial liquefaction N_L = 15 or 20 was developed irrespective of soil types.

Thus, a close correlation between the dissipated energy and the excess pore pressure or induced strain has been found and a great possibility of the EBM for liquefaction evaluation has been demonstrated in the previous research. Nevertheless, its application has been very limited so far in contrast to the conventional SBM in engineering practice. One of the reasons besides the historical background may be because concrete/detailed procedures for the EBM have not been discussed concerning how to evaluate the demand (earthquake wave energy) in particular and how to compare it with the capacity (dissipated energy), simply and reasonably.

In the earlier EBM papers ^[3,4,5], both capacity and demand energies were not explicitly quantified but intuitively represented respectively by SPT N-values and by well-known empirical formulas in seismology ^[6]. Namely, a liquefaction triggering curve was given as a boundary curve drawn on a chart of capacity versus demand energies segregating site-specific manifestation plots into occurrence and non-occurrence of liquefaction in case history records at various sites during past earthquakes. In a similar approach based on previous liquefaction case histories, Kayen and Mitchell (1997) ^[21] used Arias Intensity (1970) ^[22] as a demand for liquefaction potential evaluation, although the Arias Intensity was defined to be different from the energy in its physical meaning. Apart from the above, in the EBM proposed by Kazama et al. ^[17], the earthquake demand was not quantified to directly compare with the capacity energy, but a dynamic response analysis using design earthquake motions was carried out instead, to implicitly compare with the dissipated energy.

Moreover, Kokusho (2013) ^[18] developed another type of EBM wherein upward seismic wave energies are quantified as the demand energies and directly compared with dissipated energies at individual soil layers. Theoretical backgrounds of the comparison were discussed by Kokusho (2017) ^[23] on how to evaluate the energy of design motions and compare it with the dissipated energy based on laboratory soil tests. How to quantify the demand energy during earthquakes was first investigated based on an energy flow model of SH-wave using vertical array records during the 1995 Kobe earthquake by Kokusho & Motoyama (2002) ^[24] followed by Kokusho & Suzuki (2011, 2012) ^[25,26] using numerous vertical array data during strong earthquakes throughout Japan.

Kokusho and Mimori (2015) ^[27] conducted EBM studies on a hypothetical uniform soil deposit as well as actual liquefaction case histories where geotechnical data and recorded earthquake motions nearby were available. In the former, the significant impact of predominant frequency on the liquefaction potential was vividly shown. In the latter, the results demonstrated that, for several ground motions, the EBM tends to be mostly compatible with the SBM, if stress reduction coefficients r_n, similar to Magnitude Scaling Factors, MSF, in the USA (e.g. Idriss and Boulanger 2008) ^[28], are properly chosen in SBM. However, the gap between them tended to widen for ground motions with exceptionally large or small demand energy compared to the corresponding acceleration. In those cases, the EBM successfully replicated actual field performance whereas SBM could not properly appreciate the great impact of the demand energy.

Lau et al. (2019) ^[29] conducted comparative studies of the EBM developed by Kokusho ^[18] at several sites in Christchurch utilizing strong motion records during the 2010–2011 Canterbury Earthquake sequence and also at three sites in the Wellington and Marlborough for the 2013 Lake Grassmere and 2016 Kaikoura earthquakes in New Zealand. The method was compared with the SBM evaluation recommended by Idriss & Boulanger ^[28] concerning critical liquefaction depth and layer thickness, data scatter, and the number of false-negative (unconservative) predictions. As a result, the EBM was confirmed to make the least unconservative predictions compared to the CPT-based SBM, while yielding comparable results to the SBM as a whole.

Besides, the base-isolation mechanism of upward wave energy caused by soil liquefaction has also been investigated from the viewpoint of energy demand by Kokusho (2014) ^[30]. Furthermore, Kokusho (2020) ^[31] developed simplified evaluation steps to predict not only liquefaction potential but also associated induced strain and soil settlement, if liquefied, by assuming an equal allocation of the demand energy to potentially liquefiable layers. This reflects one of the excellent features of the EBM; total demand energy available is quantified to pursue the overall liquefaction behavior of a given site by allocating it among multiple potentially liquefiable layers. In contrast, no such interlayer interference of acceleration can be considered in the global liquefaction development of a site in the SBM.

Thus, after having revisited the previous research in this Chapter, the energy capacity is first dealt with in detail in Chapter 3 by utilizing various laboratory soil test data on reconstituted as well as intact natural soils to discuss the uniqueness of dissipated energy despite the variability of pertinent parameters in developing liquefaction. Then, in Chapter 4, wherein the energy demand, defined as cumulative upward energy by the SH-wave, is discussed concerning theories and earthquake records on how to compare with the capacity energy for liquefaction evaluation. Most of the data incorporated here are cited from previous papers in the references published by our research group.

3.   Energy capacity for liquefaction

In this Chapter, cyclic triaxial test results on reconstituted specimens are first addressed to summarize correlations of the dissipated energy with pore-pressure buildup and induced strain to examine the effects of relative density, fines content, and effective confining stress. Similar results on intact soils from in situ are also addressed to discuss the impact of natural soils on dissipated energy. In these tests, dissipated energy is compared with corresponding resistant stress CRR in the SBM for liquefaction onset to find a unique relationship between them. Then, torsional simple shear tests, which can more closely mimic in situ stress conditions during earthquakes, are addressed to recognize the similarity of the dissipated energy with that in the triaxial tests. Finally, the effect of wave irregularity of various earthquakes on the liquefaction behavior is discussed based on the torsional tests from the viewpoint of the uniqueness of dissipated energy.

3.1.   Triaxial tests on reconstituted specimens by harmonic loading

A series of stress-controlled undrained cyclic triaxial tests were carried out ^[18] using reconstituted specimens of Futtsu beach sand (along the Tokyo Bay), non-weathered sub-round particles with the mean grain size D₅₀ = 0.19 mm, and the uniformity coefficient C_u = 1.9. The size of the specimen was 10 cm in diameter and 20 cm in height. In some of them, low-plasticity fines, (the plasticity index I_p = 6) originated from decomposed granite was mixed with parametrically changing fines content F_c = 0~20%.

All the samples with various target relative densities D_r were consolidated to effective stress of ${\sigma '_c}$ = 98 kPa, with back-pressure of 196 kPa, and cyclically loaded with frequency 0.1 Hz under an undrained condition with constant axial stress amplitudes ${\sigma _d}$. Following the SBM practice, Figure 2 shows cyclic stress ratios CSR = σ_d/2σ'_c in isotropically consolidated triaxial tests versus the number of load cycles N_c plotted on the log-log chart to attain double-amplitude (DA) axial strain ε_DA = 5% and pore-pressure buildup ratio Δu/σ'_c = 1.0 for all the tests; (a) clean sands (F_c = 0) with D_r ≈ 30, 50, 70%, (b) F_c = 0~20% with D_r ≈ 50% and (c) F_c = 0~20% with D_r≈70%. Then, the plots for ε_DA = 5% in the chart are approximated by lines regressed by the following empirical formula with positive constants a and b, which are listed in the chart;

Despite data scatters, the CSR-values for ε_DA = 5% and r_u = 1.0 tend to increase systematically with increasing D_r and decreasing F_c. From the regression lines for ε_DA = 5%, cyclic resistance ratio CRR = σ_d/2σ'_c in isotropically consolidated triaxial tests is determined as CRR = CSR at the number of cycles N_c = N_L = 15 as summarized in Figure 2 for the initial liquefaction to be used in the normal SBM practice.

3.1.1.   Dissipated energy versus pore-pressure and shear strain in triaxial shear tests

Figure 3(a) exemplifies a typical cyclic stress-strain relationship for D_r = 51% and F_c = 0%. Dissipated energy per unit volume in the test specimen is calculated from a hysteretic area ABCD of a thick dashed curve for a k-th stress cycle and summed up to have a cumulative value from the start to that cycle as;

The details of this energy calculation from the stress versus strain hysteresis curve will be addressed again in torsional shear tests in Eq. (7).

In the top of Figure 3 (b), the cumulative dissipated energy per unit volume ∑ΔW calculated by Eq. (2) for the stress ~ strain curve in Figure 3 (a) is plotted versus the stress cycle N_c, together with axial stress σ_d, axial strain ε and excess pore-pressure Δu. The energy staying minimal in the earlier loading is followed by a drastic increase after the pore-pressure approaches σ'_c = 98 kPa at the onset of liquefaction. Note that the energy also keeps increasing steeply with increasing strain amplitude, though the pore pressure stops rising.

Based on a series of similar tests, the maximum pore pressure in each cycle is plotted against corresponding cumulative dissipated energy together with double-amplitude axial strain ε_DA for each loading cycle in Figure 4 for sands of F_c = 0. Here, the pore pressure is normalized by the initial effective confining stress which is named the pore-pressure buildup ratio r_u = Δu/σ_c^', and the energy is also normalized as ∑ΔW/σ_c^' in the horizontal axis. The plots are shown with different symbols for different nominal relative densities D_r = 30, 50, and 70%. The pore-pressure buildup ratio r_u correlates well with the dissipated energy and becomes r_u = 1.0 at around ∑ΔW/σ_c^' = 0.01 ~ 0.04 with small D_r-dependent differences.

In good contrast, a dominant effect of D_r is clear for the strain amplitude ε_DA plotted in the right vertical axis versus the normalized energy ∑ΔW/σ_c^', despite some data scatters caused by gaps between the targeted and actual D_r values. The induced strain ${\varepsilon _{DA}}$ is almost in proportion to the normalized energy ∑ΔW/σ_c^' for each relative density D_r up to a strain around ε_DA = 10%. The strain development can thus be correlated consistently with dissipated energy not only up to the initial liquefaction (ε_DA = 5%) but also beyond and serves as an indicator for the severity of liquefaction.

In Figure 5, the normalized dissipated energies ∑ΔW/σ_c^' to attain specific values of strain amplitudes, ε_DA = 2, 5, 10% as well as the pressure buildup ratio r_u = Δu/σ_c^' = 1.0 are plotted in the vertical axis versus numbers of cycles N_c corresponding to achieve the specific values in the horizontal axis on the log-log chart. The plots in (a) are for clean sands (F_c = 0) with parametrically varying D_r, and those in (b) are for sands of parametrically changing fines content (F_c = 5~20%) of nominal D_r = 50% and 70%. There are groups of 2 to 4 data points with the same symbols corresponding to the same strain ${\varepsilon _{DA}}$ or r_u = 1.0 achieved at different cycles N_c, which were read off from multiple tests carried out with different CSR for specimens of the same F_c and nominal D_r. The lines connecting the same symbols are essentially flat without consistent N_c-dependent up/down trends of ∑ΔW/σ_c^', despite some different trends in higher strains ${\varepsilon _{DA}}$ for dense sands.

The essentially flat lines seem to indicate that the dissipated energy ∑ΔW/σ_c^' almost uniquely determines the induced strain amplitude ${\varepsilon _{DA}}$ or pore-pressure buildup ratio r_u, no matter how many cycles N_c and how large the CSR-values are needed to attain those specific ${\varepsilon _{DA}}$ or r_u. This further indicates that a CSR ~ N_c line corresponding to a particular strain ${\varepsilon _{DA}}$ or r_u drawn in Figure 2, which normally serves as a basis for liquefaction potential evaluation in the SBM and is interpreted as the lines of equal damage in the fatigue theory (e.g., Annaki and Lee 1977 ^[32], Green and Terri 2005 ^[33]), which also represents the line of equal dissipated energy corresponding to specific conditions ^[18]. This observation paves the way for the EBM using the CSR ~ N_c data in the SBM.

3.1.2.   Dissipated energy versus CRR in isotropically consolidated triaxial tests

From the full-logarithmic CSR~N_c chart in Figure 2, individual cyclic resistant ratios CRR for N_c = 15 (CRR₁₅) can be determined from individual plots of the same symbols using the constant b in Eq. (1) for ${\varepsilon _{DA}}$ = 2, 5, 10%, and Δu/σ_c^' ≈ 1.0. This can be done by drawing parallel lines of the same gradient b passing through individual plots to determine the corresponding CRR_{Nc = 15}. Thus, the CRR₁₅-values determined are directly plotted versus corresponding dissipated energies ∑ΔW/σ_c^' to develop the CRR₁₅ versus ∑ΔW/σ_c^' chart as shown in Figure 6. Despite some data scatters, the CRR₁₅–values for ${\varepsilon _{DA}}$ = 5% (open circles) seem to be uniquely correlated with ∑ΔW/σ_c^' despite widely varying relative densities and fines contents. The correlation may be approximated by a parabolic function Eq. (3) for $CR{R_{20}} \geqslant$0.1; a practically meaningful condition in normal liquefaction problems with the determination coefficient R² = 0.86 ^[18].

Note that the ∑ΔW/σ_c^' -values in the vertical axis correspond to the dissipated energies to attain the axial strain ${\varepsilon _{DA}}$ = 5% regardless of CSR and N_c, while CRR₁₅–values in the horizontal axis represent the stress amplitudes at N_c = 15. For other strains ${\varepsilon _{DA}}$ = 2% and 10%, similar curves with the vertical coordinate 0.4 times and twice that of 5%, respectively, are drawn in the chart. The curves show a fair fitting with the corresponding plots because the strain ${\varepsilon _{DA}}$ is almost in proportion to the normalized energy ∑ΔW/σ_c^' up to around ${\varepsilon _{DA}}$ = 10% as already indicated in Figure 4.

3.2.   Triaxial tests on intact specimens from in situ by harmonic loading

In most of the experimental studies conducted so far on the dissipated energy for liquefaction, reconstituted clean sands or those mixed with non-plastic fines have been used. Experimental data on intact soils recovered in situ wherein the energy for liquefaction was focused have rarely been available. Here, a series of cyclic undrained triaxial tests conducted on intact soils with various soil properties sampled from different sites are addressed to examine the uniqueness of the cumulative dissipated energy ΣΔW in determining liquefaction behavior in terms of induced soil strain.

3.2.1.   Intact specimens addressed

A systematic test program was implemented by the Public Works Research Institute (PWRI) of the Ministry of Land, Infrastructure and Transportation (MLIT) of the Japanese Government after the 2011 Tohoku earthquake (M_w = 9.0) in Japan to review the liquefaction damage in the eastern part of Japan (Sasaki et al. 2016 ^[34]). Eventually, 190 specimens tested from 49 intact soils sampled from alluvium and hydraulic/land-fills at different sites were used. Among them, softer soils were sampled by fixed-piston thin-wall samplers, while stiffer soils were sampled by rotary triple-tube samplers. The sampling holes were drilled near pilot borings wherein SPT blow-counts and P/S-wave velocities were measured.

Stress-controlled undrained cyclic triaxial tests of intact specimens with a diameter of 5 cm and a height of 10 cm were carried out, where soil specimens were fully saturated and isotropically consolidated by in situ effective mean stresses and cyclically loaded in the undrained condition. Among the soil samples recovered by tube sampling at various depths of different sites comprising 51 specimens altogether, 14 of them named here as S1 ~ S14 were selected as of higher reliability. The samples were very variable even in the same sample tube in density, F_c, and PI (plasticity index) ^[34]. In each sample, 3 ~ 4 specimens were tested with different CSRs leading to different N_L. Those from the same soil sample were assumed uniform, though actually, a certain degree of heterogeneity reflecting soil stratifications was inevitable.

In Figure 7, all the triaxial test results in terms of CSR versus the number of cycles N_c to attain ε_DA = 5% (corresponding to the initial liquefaction) are shown together with pertinent parameters of the 14 sample S1~S14 in the legend. Note that the soils are sandy, silty and clayey, and very variable with fines content F_c = 0~100%, plasticity index PI = 0~50, and in situ SPT N₁ = 1~11. The N₁-values were calculated from SPT N-values conducted in pilot borings nearby using the formula ${N_1} = {{170N} / {\left({{{\sigma '}_v} + 70} \right)}}$ where ${\sigma '_v}$ = effective overburden in kPa. (Meyerhof 1957 ^[35]). The CSR ~ N_c plots are approximated on a full-logarithmic diagram with the straight line by Eq. (1) with the constants a and b also listed in the legend. The liquefaction strengths are widely distributed and the gradients of the lines, b, tend to decrease with increasing CRR₁₅ except for a few exceptions (S2, S9, S11).

3.2.2.   Hysteresis curves and dissipated energy of intact soils

On the left in Figure 8 (a)–(d), four hysteresis curves of cyclic axial stress σ_d versus axial strain ε are exemplified (S3, S4, S8, S10 out of S1 ~ S14), individually. At first glance, one may notice that, in most, the cyclic strain tends to develop larger in the extension side (left) than in the compression. This trend is generally observed in triaxial cyclic loading tests. It may presumably be attributed to the asymmetry in the three-dimensional stress condition inherent to cyclic triaxial tests between compression loading (one-axis vertical compression and two-axes horizontal extension) and tension loading (vice versa), as well as to anisotropy in the soil fabric.

The cumulative dissipated energies per unit area ΣΔW were calculated from the axial stress ${\sigma _d}$ versus strain $\varepsilon$ curves using Eq. (2) and nondimensionalized by initial effective confining stresses ΣΔW/σ'_c in the same way as the reconstituted specimens. In the ${\varepsilon _{DA}}$ versus ΣΔW/σ'_c correlations individually shown on the right of Figure 8 (a)–(d), the energy dissipated to attain a given strain amplitude tends to be considerably larger as the CRR of the sample becomes greater from S2 to S10. In S8 of CRR = 0.303, for instance, about 10 times difference in energy can occur among the 3 specimens to attain the same strain ${\varepsilon _{DA}}$ = 5%. Also pointed out is that the ${\varepsilon _{DA}}$-value is not uniquely correlated with the energy ΣΔW/σ'_c for larger ${\varepsilon _{DA}}$ in particular, while the uniqueness tends to hold better for soils with lower CRR such as S4 of CRR₁₅ = 0.204. Thus, the uniqueness of energy in determining liquefaction-induced strain, recognized in Figure 5 for reconstituted specimens, does not hold in these triaxial tests on intact samples for larger strains except for the low CRR soils ^[20].

In Figure 9 (a) to (d), the number of cycles N_L to attain reference DA strains ${\varepsilon _{DA}}$ = 2, 5, and 10% are plotted in the horizontal log axis versus corresponding dissipated energies ΣΔW/σ'_c in the vertical axis. The plots represented by different symbols for different samples and reference strains are connected with dashed lines for easy identification of grouping. In (a) the non-cohesive soils S1~S6 of lower CRR (0.153$\leqslant CR{R_{15}} \leqslant$0.252), the dissipated energies ΣΔW/σ'_c are relatively small and almost unchanged along increasing N_L. In contrast, the non-cohesive soils S7~S12 (0.285$\leqslant CR{R_{15}} \leqslant$0.357) in (b), ΣΔW/σ'_c tends to increase with increasing N_L for larger reference strains in particular despite some data scatters. The same trend as (b) can be observed more clearly in (c) the cohesive soils of S13~S14 (0.346$\leqslant CR{R_{15}} \leqslant$0.495).

Apart from these intact samples, Figure 9 (d) shows cyclic triaxial test results on reconstituted Futtsu sand of D_r ≈ 30, 50, 70%, and F_c = 0% or 20% with CRR₁₅ = 0.086~0.206, already addressed in Figure 5. Their dissipated energies ΣΔW/σ'_c are smaller than the intact samples in Figure 9 (a), and almost unchanged with increasing N_L, presumably reflecting the absence of the aging effect in the reconstituted samples.

In Figure 10, CSR values are plotted versus corresponding dissipated energies ΣΔW/σ'_c in full logarithmic diagrams to attain the reference DA strains ${\varepsilon _{DA}}$ = 2, 5, and 10%. In the non-cohesive soils of lower CRR, S1~S6 in (a), the energies are mostly constant and independent of CSR despite large data scatters. In contrast, the non-cohesive soils of higher CRR, S7~S12 in (b), exhibit decreasing trends of ΣΔW/σ'_c with increasing CSR (except abnormal S8). In the cohesive soils, S13~S14 in (c), the decreasing trend of ΣΔW/σ'_c with increasing CSR is more visible.

From the test results in Figures 9 and 10, it may be inferred that intact non-cohesive soils of low CRR in (a) are similar to reconstituted sands, wherein the energy ΣΔW/σ'_c almost uniquely determines induced strain with no regard to the number of cycles N_L or CSR. However, intact non-cohesive/cohesive soils of higher CRR which presumably reflect stronger aging appear to lose the uniqueness of energy in liquefaction behavior. Thus, the N_L/CSR-dependency of dissipated energy for liquefaction seems more conspicuous for intact soils of higher CRR in cyclic triaxial tests ^[20].

In order to quantify the extension strain bias thus observed, peak-to-peak mid-strain ${\varepsilon _{mid}} =$ ${{\left({{\varepsilon _{\max }} + {\varepsilon _{\min }}} \right)} / 2}$ is calculated from the maximum and minimum strains in each loading cycle, ${\varepsilon _{\max }}$ and ${\varepsilon _{\min }}$, respectively, wherein the extension strain is to be positive. The ${\varepsilon _{mid}}$-values for DA strains ${\varepsilon _{DA}}$ = 2, 5, 10% in the triaxial tests are plotted versus corresponding CSR in Figure 11 with different symbols. Note that the plots are overwhelmingly on the positive side of ${\varepsilon _{mid}}$, indicating that almost all cyclically induced strains are biased on the extension side. The ${\varepsilon _{mid}}$-values are widely varied for the non-cohesive lower CRR soils in (a) presumably reflecting variability of in situ natural soils, but they show no meaningful correlations with CSR. For the reconstituted Futtsu sand of D_r ≈ 50 & 70% and F_c = 0% plotted on the same chart, and the ${\varepsilon _{mid}}$-values are also positive and essentially flat versus increasing CSR. In the non-cohesive higher CRR soils in (b), however, ${\varepsilon _{mid}}$ tends to increase with increasing CSR for higher strains ${\varepsilon _{DA}}$, and the same trends are more clearly recognized in the cohesive soils of higher CRR in (c). Thus, the strain bias on the extension side seems intrinsic in triaxial tests in general, which tends to be more pronounced with increasing CSR for intact in situ soils of higher CRR with a larger DA strain in particular, while this trend is not so evident in reconstituted sands.

In Figure 12 (a)–(c), the dissipated energy ΣΔW/σ'_c is directly plotted versus the peak-to-peak mid-strain ${\varepsilon _{mid}}$ for the three groups of samples in terms of their averages for reference DA strains ${\varepsilon _{DA}}$ = 2, 5, 10% with the large plots connected with solid lines, while the smaller plots with dashed lines are individual data with different CSRs corresponding to individual ${\varepsilon _{DA}}$-values. The energy ΣΔW/σ'_c tends to increase with increasing ${\varepsilon _{mid}}$ if the average values are concerned, though ΣΔW/σ'_c for the same ${\varepsilon _{DA}}$ is also dependent on CSR. That dependency is small for the soils of low CRR in (a) though getting stronger with increasing CRR in (b) and (c) as already addressed. Thus, the strain bias in extension is one of the major causes of the apparent non-uniqueness of energy specific in the triaxial stress system ^[20].

Unlike cyclic triaxial tests, however, soil layers in situ are cyclically sheared left and right during earthquakes by SH waves on a horizontal plane, wherein no such three-dimensional stress asymmetry in compression/extension exists. The best way to reproduce the in situ stress system in the laboratory is (hollow cylindrical torsional) simple shear tests.

The above experimental finding regarding triaxial tests may have a significant impact not only on the evaluation of liquefaction energy capacity but also on numerical liquefaction analyses by effective stresses, wherein soil constitutive models are calibrated to CSR~N_c correlations often obtained by laboratory element tests. Considering the significant effect of strain bias in the stress versus strain hysteretic relationship, it is recommended to use simple shear tests rather than triaxial tests for better reliability of the analyses.

3.2.3.   Dissipated energy versus CRR correlation for intact soils

Considering that CRR used in SBM is normally defined as CSR at N_L = 15, the dissipated energy corresponding to that particular CSR or N_L may be focused and correlated with CRR. For that goal, ΣΔW/σ'_c versus CSR plots of individual samples for ${\varepsilon _{DA}}$ = 5% (corresponding to the initial liquefaction) in triaxial tests depicted in Figure 10 are approximated here by a straight line on the log-log diagram as conceptually illustrated in Figure 13 (a) and formulated by a power function of Eq. (4) ^[20].

Here, CRR₁₅ is the CSR-value in isotropically consolidated cyclic triaxial tests at N_L = 15 as normally accepted in the SBM liquefaction evaluation, wherein a constant α is equal to ΣΔW/σ'_c when CSR = CRR₁₅, and −β is a negative gradient of the line. If the energy were uniquely determined irrespective of N_L or CSR, the dashed flat line with β = 0 in Figure 13 (a) would be obtained. Thus, it seems practically meaningful to represent the dissipated energies corresponding to various CSR-values in cyclic triaxial tests by the energy for CRR = CSR at N_L = 15 by using Eq. (4). It could be N_L = 20 or other values in place of N_L = 15, whichever possible in engineering practice, as long as the corresponding dissipated energy ΣΔW/σ'_c is postulated to uniquely determine the strain amplitude for liquefaction irrespective of CSR or N_L.

Figure 13 (b) depicts the β-values in Eq. (4) plotted versus CRR₁₅ for all the 49 intact samples or the 14 selected samples with open or closed plots, respectively, where the symbols are grouped by stepwise F_c-values. For CRR around 0.25 or smaller, β is nearly zero, indicating in Eq. (4) that the energy for liquefaction can be uniquely determined independent of N_L or CSR as can be seen in Figures 9(a) and 10(a). With increasing CRR, β tends to increase from nearly zero to around 5, indicating that the non-uniqueness of energy or dependency on N_L or CSR in cyclic triaxial tests is pronounced as demonstrated in Figures 9 and 10.

The α-values, identical to the dissipated energy (∑ΔW/σ'_c)_{εDA = 5%} for CSR = CRR in Eq. (4), are calculated from Figure 10 and plotted versus CRR₁₅ in Figure 13 (c). The open and closed dots represent the global and selected data, respectively, and different symbols represent different stepwise F_c-values, again.

For the selected data of CRR < 0.4, a good correlation can be obtained between CRR₁₅ and the dissipated energy ΣΔW/σ'_c to attain ε_DA = 5% (the initial liquefaction) as;

If CRR₂₀ for N_L = 20 is chosen instead of 15 in SBM such as in the road-bridge design code in Japan, it becomes;

The determination coefficient for both equations is R² = 0.92 ^[20].

Remarkably, the liquefaction strength CRR used in SBM is almost uniquely correlated with the energy ∑ΔW/σ'_c despite the considerable variations in densities, fines contents, plasticity indices, and ages of in situ soils. Star symbols superposed in Figure 13 (c) are the cyclic triaxial test results of the reconstituted Futtsu sand already addressed in Figure 6, which also showed a unique relationship of Eq. (3) between the energy and CRR₁₅ for the initial liquefaction (${\varepsilon _{DA}}$ = 5%). Note that the close symbols of the selected data from the intact samples are essentially compatible with the star symbols, though the absolute energies are slightly larger presumably due to the aging effect of the intact soils.

Thus, the widely varied ΣΔW/σ'_c versus CSR plots in Figure 10 obtained at different CSR-values in triaxial tests on intact soils from in situ have yielded a well-correlated ΣΔW/σ'_c versus CRR₁₅ plots as shown in Figure 13 (c). Hence, Eqs. (5) or (6) may be conveniently used to evaluate the energy capacity ΣΔW/σ'_c for EBM from CRR to be determined in SBM from in situ penetration tests.

3.3.   Torsional shear tests by harmonic loading

As shown above, the triaxial tests despite their good access from practical engineers have intrinsic difficulties in replicating simple shear conditions for evaluating liquefaction-induced dissipated energy. Figure 14 exemplifies stress-strain hysteresis curves of reconstituted clean sands by the torsional simple shear tests of (a) Futtsu sand of D_r = 45%, and (b) Toyoura sand of D_r = 80%, to compare with (c) the triaxial test result on reconstituted Futtsu sand of D_r = 51%. In the simple shear tests, shear strain develops almost symmetrically, whereas, in the triaxial test, the strain tends to be biased on the extension side with increasing numbers of cycles even in medium-density reconstituted sands.

In the following, typical torsional undrained simple shear test results on reconstituted samples are addressed to examine the uniqueness of energy ^[19]. Harmonic and irregular cyclic stresses are given to the specimens without and with non-plastic fines to correlate excess pore pressure and induced strain with the dissipated energy. Effects of cyclic stress ratio CSR, number of cycles N_c, and wave irregularity on the energy-dependent strain evaluation are particularly focused.

The specimen size was r_i = 30 mm and r_o = 50 mm in inner and outer radius, respectively, and H = 100 mm in height. The Futtsu beach sand was used with non-plastic fines parametrically changing F_c = 0 to 30%. The test specimens were prepared by moist tamping to target relative densities, D_r = 30 and 50%. The specimens were saturated and isotropically consolidated by effective confining stress ${\sigma '_c}$ = 98 kPa in most cases, or 48 or 196 kPa in some other cases and cyclically sheared by torsional stress ${\tau _d}$. The cyclic stress ${\tau _d}$ was either the harmonic motion of frequency 0.1 Hz or irregular seismic motion with its time scale elongated 10 times.

Following the SBM practice, Figure 15 shows cyclic stress ratios $CSR = {{{\tau _d}} \mathord{\left/ {\vphantom {{{\tau _d}} {{\sigma _c}^\prime }}} \right. } {{\sigma _c}^\prime }}$ for γ_DA = 7.5% (the initial liquefaction corresponding to axial strain ε_DA = 5% in triaxial tests) versus the number of cycles N_c on the log-log chart for all the torsional shear tests on reconstituted Futtsu sand specimens; (a) D_r ≈ 50%, F_c = 0–30%, σ'_c = 49, 98,196 kPa and (b) D_r ≈ 30%, F_c = 0%, σ'_c = 49, 98,196 kPa. Then, the plots are approximated by straight lines drawn after the empirical formula Eq. (1) with constants a and b regressed from the plots. Note that the plots are aligned almost linearly, and increasing F_c tends to reduce CSR considerably under the same D_r, while the effective confining stress seems to have only a minor effect.

3.3.1.   Dissipated energy by harmonic loading

Figure 16 (a) and (b) exemplify the stress τ_d versus strain γ relationship or effective stress path (τ_d versus p' = σ_c'−Δu), respectively, for clean sand of D_r = 45% and F_c = 0% sheared by $\operatorname{CSR} = \tau_d / \sigma_c^{\prime}$ = 0.236 in the effective confining stress ${\sigma _c}^\prime$ = 98 kPa. Dissipated energy in the cycle ABCDEFGHA′ in (a) can be calculated as;

where ∑_sa stands for the summation of sliced areas such as PQQ′P′ for a strain increment $\Delta \gamma$, and $\oint {}$ for the integral in terms of shear strain $\gamma$ in a single loading cycle. Obviously, the increment of dissipated energy in the sliced Area (PQQ'P') can be automatically calculated in Eq. (7) because τΔγ alternates signs plus/minus between Area (PQRS) and Area (P'Q'RS) for increasing and decreasing strains, respectively.

The cumulative dissipated energy to a k-th cycle ∑ΔW is obtained as the summation of ΔW_k or expressed eventually as a seamless integral from the 1st to the arbitrary cycle as;

Eq. (8) is valid not only for harmonic waves but also for irregular motions of variable amplitudes.

Figure 17 (a) exemplifies typical time histories of (ⅰ) cyclic shear stress ${\tau _d}$, (ⅱ) excess pore-pressure normalized by effective confining stress Δu/σ_c^', (ⅲ) shear strain $\gamma$, (ⅳ) normalized cumulative dissipated energy ∑ΔW/σ_c^', (ⅴ) axial stress ${\sigma _a}$, and (ⅵ) axial strain ${\varepsilon _a}$. The axial stress and strain in the hollow cylindrical specimen tend to change slightly after the initial liquefaction as the specimen softens, though the associated energy is negligibly small. ∑ΔW/σ_c^' increases its gradient with N_c when the initial liquefaction for ${\gamma _{DA}}$ = 7.5% occurs at N_L = N_c = 12.7 in this case.

Figure 17 (b) shows typical variations of excess pore-pressure ratio ${{\Delta u} \mathord{\left/ {\vphantom {{\Delta u} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}$ and shear strain $\gamma$ versus normalized cumulative dissipated energy ∑ΔW/σ_c^' exhibited in the same test as (a). Larger close symbols in the chart represent pore-pressure buildup ratios r_u = (${{\Delta u} \mathord{\left/ {\vphantom {{\Delta u} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}$)_peak at the end of the individual stress cycles, while larger open symbols are double-amplitude shear strains ${\gamma _{DA}}$, peak to peak, in individual cycles. The r_u-value correlates with the dissipated energy in an asymptotic curve with its maximum r_u ≈ 1.0 for ∑ΔW/σ_c^' ≈ 0.03 at N_c ≈ 13, corresponding to double-amplitude shear strains ${\gamma _{DA}}$ = 7.5 % for the initial liquefaction.

The ${\gamma _{DA}}$ versus ∑ΔW/σ_c^' curve changes its gradient at two inflection points, though it is not so strongly nonlinear. In this example, the first inflection occurs at ∑ΔW/σ_c^' ≈ 0.01 around N_c = 9th cycle when the pore pressure starts fluctuating up/down due to cyclic mobility. In the τ_d ~ p^' diagram of Figure 16 (b), the effective stress path around there is in a thick solid curve which is positioned about half in the dilative zone crossing the phase transformation line (PTL) (Ishihara et al. 1975 ^[36]). Moreover, in the same diagram, we observe that the second inflection appears at around N_c = 15th cycle where the effective stress path moves up and down along the failure line (FL) shown with another thick curve through the origin. Thus, a ${\gamma _{DA}}$ versus ∑ΔW/σ_c^' curve in general may be divided into the next three sections. Section Ⅰ: From the start to the 1st inflection point, where the effective stress path is contractive in cyclic loading with low strain versus energy gradient. Section Ⅱ: From the 1st to 2nd inflection point, where the effective stress path is highly affected by PTL; semi-dilative in loading versus contractive in unloading, with higher strain versus energy gradient. Section Ⅲ: The effective stress path moves along the failure line, exhibiting again a low strain versus an energy gradient.

Similarly, the correlations of pore-pressure ratio r_u and double-amplitude strain ${\gamma _{DA}}$ are plotted versus dissipated energy ∑ΔW/σ_c^' in Figure 18 (a) and (b), respectively, for clean sands of not only D_r ≈ 50% (closed symbols) but also D_r ≈ 30% (open symbols) both tested under σ_c^' = 98 kPa. The r_u-value is almost uniquely correlated with ∑ΔW/σ_c^' despite largely diverted CSR or N_L-values. The pressure buildup for D_r≈30% tends to occur slightly faster and attains r_u = 1.0 (the initial liquefaction) at lower dissipated energy than that for D_r ≈ 50% around ∑ΔW/σ_c^' = 0.02.

In the ${\gamma _{DA}}$ versus ∑ΔW/σ_c^' correlation of Figure 18 (b), the double-amplitude strain ${\gamma _{DA}}$ is distinctively larger for D_r ≈ 30% than D_r ≈ 50%, indicating a significant effect of D_r on the induced strain. Among the same D_r group, the dissipated energy looks almost uniquely correlated with ${\gamma _{DA}}$ up to the initial liquefaction (${\gamma _{DA}}$ = 7.5%) and beyond. In the correlations, the energy ∑ΔW/σ_c^' at the 1st inflection point around 0.01 or less, tends to be larger with increasing N_L. Hence, the γ_DA versus ∑ΔW/σ_c^' correlation may not be precisely unique presumably due to the inflection point mentioned above. Nevertheless, as shown later, the dissipated energy ∑ΔW/σ_c^' achieving a given γ_DA -value for harmonic waves is found to be almost uniquely determined independent of N_L or CSR, presumably because the N_L-dependent effects are not strong enough to overcome the data scatters.

Figure 19 shows similar correlations of r_u or γ_DA plotted versus ∑ΔW/σ_c^' for D_r ≈ 50% sands containing fines tested under σ'_c = 98 kPa. Plots of sands for fines content F_c = 10~30% are compared with solid plots for F_c = 0 already addressed in Figure 18. In Figure 19 (a), ∑ΔW/σ_c^' is fairly well correlated with r_u in individual F_c -values, although the pressure buildup tends to occur earlier for the same energy as F_c increases from 0 to 30% and attain r_u≈1.0 at smaller energy of ∑ΔW/σ_c^' < 0.02. In (b), the γ_DA -values for soils of F_c = 10~30% tend to grow larger for the same energy ∑ΔW/σ_c^' than those for the clean sand F_c = 0%, indicating a significant effect of F_c particularly due to the increase from F_c = 0 to 10%. Thus, the effect of fines content on the ∑ΔW/σ_c^' versus γ_DA correlation is considerable even for the same D_r.

3.3.2.   Dissipated energy versus CRR in torsional shear tests

From the CSR ~ N_c chart for the SBM obtained by the torsional test on isotropically consolidated specimens in Figure 15, the cyclic resistant ratio CRR₁₅ for ${\gamma _{DA}}$ = 7.5% for N_c = 15 can be determined by Eq. (1) for each test specimen using the value b available in the legend. Thus, the determined CRR₁₅-values are directly correlated with corresponding dissipated energies ∑ΔW/σ_c^' to construct CRR₁₅ versus ∑ΔW/σ_c^' correlations for ${\gamma _{DA}}$ = 7.5% in Figure 20. The similar correlations analogously obtained for ${\gamma _{DA}}$ = 3% and 15%, are also plotted in the chart. Despite that the plots included here widely differ in relative density D_r ≈ 30 ~ 50%, fines content F_c = 0 ~ 30% and confining stress = 49 ~ 198 kPa, the CRR₁₅-values for the strain level ${\gamma _{DA}}$ = 7.5% plotted by the circles seem to be correlated fairly well with the dissipated energy ∑ΔW/σ_c^' and can be approximated by a parabolic function;

for CRR ≥ 0.1 or larger with the determination coefficient R² = 0.79.

A similar ΣΔW/σ'_c ~ CRR approximation formula Eq.(3) previously obtained in the triaxial tests under ${\sigma '_c}$ = 98 kPa using the same reconstituted Futtsu sand of D_r = 30, 50, 70% for axial strain ${\varepsilon _{DA}}$ = 5% (corresponding to ${\gamma _{DA}}$ = 7.5% in shear strain) is overlaid with a dashed curve in Figure 20. Note that the dissipated energy in Eq.(3) is defined differently in the chart as 3/4 of that in Eq.(8), because the shear stress and strain in the simple shear stress condition are expressed as $\tau = {\sigma \mathord{\left/ {\vphantom {\sigma 2}} \right. } 2}$ and $\gamma = {{3\varepsilon } \mathord{\left/ {\vphantom {{3\varepsilon } 2}} \right. } 2}$ where $\sigma$ and $\varepsilon$ are the axial stress and strain in triaxial shear, respectively. Despite the significant differences in the test conditions, the solid curve in Eq. (9) looks similar to the dashed curve in Eq. (3).

3.3.3.   Effect of confining stress on dissipated energy

Figure 21 shows the influence of effective confining stress changing stepwise, σ'_c = 49, 98, and 198 kPa, on the correlations of ΣΔW/σ'_c versus ratio r_u in (a) or induced shear strain ${\gamma _{DA}}$ in (b) for the sand of D_r ≈ 50% and F_c = 0%. The pressure ratio r_u for the same ∑ΔW/σ'_c -value tends to be lower with increasing σ'_c in (a) while a similar trend can be seen for the γ_DA versus ∑ΔW/σ_c^' correlation in (b) only for smaller energy values despite larger data scatters.

In Figure 22, the energies normalized by confining stress, ∑ΔW/σ_c^', to attain (a) stepwise r_u- values or (b) stepwise γ_DA-values are plotted versus the effective confining stresses for σ'_c = 49 to 196 kPa to examine how σ'_c impacts the ∑ΔW/σ_c^'-value. The inclined dashed lines and horizontal solid lines drawn there assume that ∑ΔW/σ_c^' is either dependent on (σ'_c)^0.5 or independent of the confining stress (σ'_c)⁰, respectively.

As for the effect of σ'_c on the dissipated energy, let us assume that a stress τ versus strain γ hysteresis loop of cyclically loaded sands is idealized by the hyperbolic model f() often employed in liquefaction analyses (e.g. Kokusho 2017) ^[37] as;

where G₀ = initial shear modulus, γ_r = reference strain for G/G₀ = 0.5. Then, the normalized cumulative dissipated energy is calculated as;

Because G₀ and γ_r are dependent on the effective confining stress σ'_c as;

where G₀^* = shear modulus under a unit pressure p₀, ${\tau _f} = {\sigma '_c}\tan \phi$ = ultimate shear strength, $\phi$ = internal friction angle, and p₀ = unit pressure. Then;

It shows that the normalized cumulative dissipated energy ∑ΔW/σ_c^' tends to increase proportionally to (σ'_c /p₀)^0.5, the square root of the confining stress.

From experimental observations shown in Figure 22 (a), the normalized energy ∑ΔW/σ_c^' looks to be proportional to (σ'_c)^0.5 as in Eq. (14) in terms of pressure buildup r_u. However, in terms of induced strain γ_DA in (b), it seems more appropriate to assume that ∑ΔW/σ_c^' is σ'_c-independent. In the latter case, the energy ∑ΔW/σ_c^' can be normalized dimensionally and quantitatively. Then, the normalized dissipated energy in determining liquefaction-induced strain becomes almost independent of soil depth, to be beneficial for in situ liquefaction evaluations based on energy.

3.4.   Torsional shear tests by irregular loading

3.4.1.   Irregularity of input stress waves

A variety of earthquake motions A to F, covering intra-plate as well as plate-boundary earthquakes recorded in recent decades in Japan widely varied in terms of magnitudes (M_J = 6.8 to 9.0), hypocenter distance (R = 24~380 km), duration, and number of cycles, were used for the liquefaction tests as tabulated in Table 1. The undrained cyclic torsional shear tests were conducted twice using the same waves with different amplitudes on Futtsu sands of D_r ≈ 50%, F_c = 0% under initial isotropic effective confining stress σ_c^' = 98 kPa ^[19].

In Figure 23 (a), six acceleration time histories are illustrated with dash lines superposed by corresponding stress time histories with solid lines; A-1 to F-1, reproduced in the undrained cyclic torsional simple shear tests. The duration and number of cycles are considerably differentiated among the waves depending on the magnitudes M_J and the hypocenter distances R.

Sequences of double amplitudes of acceleration/stress obtained by zero-cross calculations of the time histories are illustrated in Figure 23 (b) with open/close columns, respectively, from larger to smaller DA amplitudes normalized by the first maximum value (starting from unity) versus the sequential number of cycles (a half of zero-crossing numbers). This reveals that the stress time histories are good enough to represent wave irregularities of widely varied earthquakes, though they are not in exact proportion to the acceleration amplitudes due to insufficient performance of the test device in reproducing smaller amplitudes. The total number of cycles N_c-total of the stress waves with amplitudes exceeding 5% of the maximum listed in Table 1 tends to increase from Wave A to F with corresponding M_J and R-values, which seems to divide the motions into two groups A ~ C and D ~ F of smaller and larger N_c, respectively.

Also indicated in Table 1 are relative damage levels D_L for liquefaction calculated on Wave A to F in the fatigue theory sometimes addressed as a theoretical background of the stress-based liquefaction evaluation (e.g. Annaki and Lee ^[32]).

Postulate here that substituting the maximum cyclic stress ratio CSR = CSR_max = 1.0 and the number of cycles N_c = N_L = 1 into Eq. (1) gives a relative damage level D_L = 1.0, then a single cycle of an arbitrary stress ratio CSR_j smaller than CSR_max = 1.0 gives a relative damage level given as D_Lj = 1/N_Lj = (CSR_j)^−1/b, and the sum of the relative damage levels for individual amplitudes can be calculated as;

The D_L-values summarized in Table 1 are larger in wave group D ~ F than in group A~C, though Wave D shows unexpectedly large D_L despite its smaller earthquake magnitude and epicenter distance. However, despite the large differences among the D_L-values, the test results indicate that the liquefaction potentials by the six motions are almost solely dependent on the dissipated energy ∑ΔW/σ_c^' as shown later in Figure 28.

3.4.2.   Effect of wave irregularity on dissipated energy for liquefaction

Figure 24 shows 12 stress-strain curves from the start to double-amplitude strain γ_DA = 20~30% obtained by the torsional tests using the 6 waves A to F tested twice each with different amplitudes. From each chart, the cumulative dissipated energy ΣΔW to a given cycle addressed in the following figures was calculated in the same way as in harmonic waves using Eq. (5). Each hysteresis curve here is defined in terms of stress, starting upward from τ_d = 0 and ending with the next τ_d = 0.

The solid dots connected with solid lines overlaid on the diagrams represent the means of maximum and minimum stresses and corresponding strains in individual cycles of the irregular motions. Unlike harmonic loading where stress and strain are both symmetrical over the origin, the stress-strain loops are not only irregular but their centers tend to deviate from the origin cycle by cycle.

In Figure 25 (a), the plots of the pressure ratio r_u versus energy ∑ΔW/σ_c^' are superposed for all the irregular waves to compare with those of harmonic waves shown in Figure 18. The correlations for the former are more dispersed than the latter presumably because of poor reproductions of the targeted D_r-values, though their mean values are essentially coincidental. The ${\gamma _{DA}}$ versus ∑ΔW/σ_c^' plots in Figure 25 (b) look similar between the irregular and harmonic motions, and the two inflection points mentioned before in the harmonic motions can be identified in the irregular motions, too. The correlations may be recognized to be almost unique up to the energy around ∑ΔW/σ_c^' ≈ 0.02~0.04, corresponding to the initial liquefaction but increasingly dispersed beyond that.

3.4.3.   Uniqueness of dissipated energy for liquefaction by torsional shear test

In Figure 26 (a) ~ (c), dissipated energies ∑ΔW/σ_c^' for given strains γ_DA = 3, 7.5, 15% in the harmonic wave tests read off from Figures 18, 19 and 21 are plotted versus the corresponding number of cycles N_c in the horizontal axis. Among them, the chart (a) shows plots for D_r ≈ 50%, F_c = 0, σ_c^' = 49, 98,196 kPa, while the charts (b) for D_r ≈ 30%, F_c = 0, and the charts (c) for D_r ≈ 50%, F_c = 10~30%, respectively, both under σ_c^' = 98 kPa. In Figure 26 (a), ∑ΔW/σ_c^' may be considered N_c-independent for a wide range of N_c ≈ 2 to 100 and DA strains γ_DA = 3, 7.5, 15%. If the difference of confining stress σ_c^' = 49~196 kPa on ∑ΔW/σ_c^' is ignored as suggested in Figure 22 (b), the average dissipated energies for D_r ≈ 50% are ∑ΔW/σ_c^' = 0.017, 0.028, and 0.047 as indicated with the thick dashed lines corresponding to γ_DA = 3, 7.5, 15%, respectively. For other test conditions in Figure 26 (b) and (c), no systematic N_c-dependent variations of energy outrivaling the data dispersions may be observed, indicating that the dissipated energy is almost uniquely correlated with double-amplitude shear strain γ_DA for sands regardless of the number of cycles N_c or the cyclic stress ratio CSR.

In Figure 27 (a), another torsional shear test using clean Toyoura sand conducted by Sako (2019) ^[38] (the data was processed later by the present author ^[20]) is depicted in terms of the strain amplitude γ_DA versus the energy ∑ΔW/σ_c^' relationships. Because the test was carried out with greater care in precisely reproducing targeted relative densities, the energies ∑ΔW/σ_c^' are very uniquely correlated with the double amplitude strain γ_DA, and those for γ_DA = 3, 7.5% read off are surprisingly flat versus the number of loading cycles N_c in Figure 27 (b). Though the D_r-value reported as D_r = 50% ^[38] seems to be too small considering ∑ΔW/σ_c^' or CRR_{NL = 15} shown in Figure 27 (c), the almost perfect N_L-independency of ∑ΔW/σ_c^' shown here serves as good experimental evidence demonstrating the uniqueness of energy for induced shear strain in the simple shear stress condition if sand density is kept precisely the same.

In Figure 28, the cumulative dissipated energies ∑ΔW/σ_c^' by the irregular waves A~F on the clean Futtsu sand with D_r ≈ 50% under σ_c^' = 98 kPa, which are read off from Figure 25 (b) for γ_DA = 3, 7.5, 15%, are plotted with solid symbols versus the number of cycles N_c in a log scale corresponding to those particular strains. To count N_c from the start of the irregular waves, major cycles with DA stress amplitudes exceeding 5% of the maximum were selected as already explained in Figure 23 (b).

On the same diagram, the ∑ΔW/σ_c^' versus N_c plots by harmonic wave for D_r ≈ 50% and σ_c^' = 98 kPa shown in Figure 26 (a) are superposed, by assuming similarity between the numbers N_c in harmonic and irregular motions, respectively. The irregular and harmonic test results are almost compatible in terms of the absolute value of ∑ΔW/σ_c^' for γ_DA = 3, 7.5, and 15% despite data scatters. It also shows that the energy ∑ΔW/σ_c^' is essentially flat against N_c for γ_DA = 3 and 7.5%, indicating that the dissipated energy uniquely determines the initial liquefaction regardless of the motions, harmonic or irregular.

Thus, Figure 28 indicates that the energies for γ_DA = 3% and 7.5% are essentially independent of N_c, implying that the performance up to the initial liquefaction is uniquely determined by energy even in the irregular waves A ~ F widely varied in terms of earthquake magnitude and epicenter distance, CSR, number of cycles, frequency contents, duration, and the damage level D_L already explained.

Figure 28 also indicates that the energies ∑ΔW/σ_c^' at larger strains such as γ_DA = 15% appear to decrease with increasing N_c though not so drastically. It is shown that, in irregular waves unlike harmonic waves, the center points of stress ~ strain hysteresis loops tend to move around the origin as already pointed out with solid lines & dots in Figure 24. This indicates that the deviated shear stress may increase the liquefaction resistance due to the cyclic mobility of dilative sands (e.g. Kokusho 2020 ^[39]). This will increase the dissipated energy for irregular motions of smaller N_c or larger CSR particularly, wherein the off-origin deviations of the mean stress tend to be greater.

To evaluate the stress-deviation effect quantitatively, the additional energy due to the stress bias may be roughly assessed from the movements of the solid dots in Figure 24, schematically illustrated by the vector k to l in Figure 28 (b). Then, the additional dissipated energy may be roughly evaluated by shaded areas E_LC, trapezoidal or triangular depending on whether the vector kl is on one side or crossing the τ = 0 line. Thus, the energy ∑E_LC/σ_c^', using E_LC calculated individually and summed up for all the major loops of an irregular wave and normalized by σ_c^', may evaluate the dissipated energy by the stress bias.

In Figure 28 (a), the dissipated energy after subtracted by the biased energy, (∑ΔW/σ_c^'−∑E_LC/σ_c^'), is also overlaid with open symbols. Note that the energy ∑E_LC/σ_c^' seems to have an insignificant impact regarding the calculated energy for γ_DA = 3 and 7.5. Also note for γ_DA = 15% that, the N_c-dependency of dissipated energy becomes milder by this modification, resulting in the energy trend almost flat against N_c and compatible with the harmonic motion except for the two motions C and F.

Thus, the test results indicate that the initial liquefaction can be uniquely predicted using the dissipated energy irrespective of input waves. However, the energy-based prediction for larger strains in high precision seems to be difficult and less meaningful from an engineering point of view, considering the variabilities of soil properties under ultra-low effective confining stress during liquefaction and irregular design motions. Hence, uncertainties in predicting strains far beyond the initial liquefaction should be accepted in design to a certain extent.

Such an acceptance seems inevitable not only in the energy-based evaluation but also in other sophisticated numerical tools when larger post-liquefaction strains need to be evaluated in Performance-Based Design. This also indicates that a deterministic approach for post-liquefaction residual deformation based on near-zero effective stress properties of liquefying sands using a single earthquake motion is irrelevant. Instead, a probabilistic approach taking account of the variability of soil properties and irregular waves has greater relevance in evaluating post-liquefaction large strains. Regarding this, the energy approach can serve as a simple and effective tool to quantify not only the liquefaction potential but also post-liquefaction strain development quite reasonably.

4.   Energy demand for liquefaction to compare with capacity

Earthquake energy was first dealt with quantitatively by seismologists (Gutenberg and Richter 1942 ^[40], 1956 ^[6,41]) to evaluate the total energy released from a seismic source based on observed earthquake records assuming the spherical energy radiation for body waves. From the viewpoint of engineering design where acceleration and stress are always focused on, very few researchers have tried to interpret earthquake motions as energy demand. Among them, Sarma (1971) ^[42] calculated site-specific demand energies from velocity records and compared them with spherically radiated energy from the earthquake source. In liquefaction evaluations, energy-based methods have been proposed ^[3,4,5], where the energy capacity for liquefaction triggering was compared with earthquake energy demand ^[6].

While the wave energy defined seismologically seems to correspond to demand energy at hard bedrock neglecting the effect of softer surface layers, how the demand energy transmits as earthquake waves from hard bedrock at a depth to ground surface had not been discussed for a long time. Kokusho and Motoyama (2002) ^[24] performed a basic study on the energy flow of seismic waves in surface layers based on the one-dimensional multi-reflection theory of SH waves utilizing vertical array records during the 1995 Kobe earthquake, which was followed by a theoretical study on the same topic by Kokusho et al., (2007) ^[43]. Similar studies using numerous vertical array data were further carried out to understand general trends of the energy demand in surface layers ^[25,26].

In this chapter, the energy demand in surface soils is calculated assuming 1D multi-reflection of SH-wave using vertical array records in actual soil profiles. After touching upon basic theories on wave energy, empirical formulas already developed using a vertical array strong motion database are explained. Then, how to determine the demand energy from a given earthquake and compare it with the capacity energy in the liquefaction evaluation is discussed considering the effect of free ground surface.

4.1.   Energy flow of a 1D propagating SH wave

The horizontal displacement u in the SH wave propagating to the positive direction of the z-axis, as illustrated in Figure 29, can be expressed as;

Here, t = time, z = upward coordinate, V_s = S-wave velocity, A = wave amplitude, and f () is an arbitrary waveform function. Then, shear strain γ is readily obtained.

where $\dot u = {{\partial u} \mathord{\left/ {\vphantom {{\partial u} {\partial t}}} \right. } {\partial t}}$ is particle velocity (e.g., Kokusho 2017 ^[37]).

As for the wave energy carried by the upward SH wave passing through a horizontal plane A-A' of a unit area, kinetic energy increment ΔE_k can be written as Eq. (18) for a soil element of a unit horizontal area times a small thickness dz = V_sΔt (a travel distance in a time increment Δt) having particle velocity $\dot u$.

Strain energy ΔE_γ is expressed similarly by shear stress τ = Gγ and shear strain γ, and using Eq. (17) as:

Hence, $\Delta {E_k} = \Delta {E_\gamma }$ holds in the same soil element, and the wave energy passing through the unit area in the time increment Δt is their sum expressed as:

Cumulative energy in a time interval t = t₁ ~ t₂ can be expressed as the sum of the kinetic and strain energies, E_k and E_e, of an equal amount (Timoshenko and Goodier 1951 ^[44], Bath 1956 ^[45], and Sarma 1971 ^[42]).

Note that the unit of E is Energy divided by Area, and kJ/m² will be used hereafter. The time derivative of the energy called energy flux or energy flow rate is written as:

Thus, the seismic wave energy depends on the particle velocity $\dot u$ and the S-wave impedance ρV_s of the soil where the ground motion is recorded. In this context, it seems unreasonable from the viewpoint of energy to define design motions in terms of acceleration or velocity without specifying the associated impedance value ρV_s. Hence, when a design motion with a given amplitude is discussed, it is essential from the viewpoint of energy to identify the impedance value or soil condition where the motion is recorded.

4.2.   Energy flow calculated by vertical array records

Subsurface energy flows were calculated utilizing vertical array records during nine strong earthquakes acquired in recent decades in Japan by assuming the vertical propagation of SH waves to know depth-dependent energy demands and apply them in energy-based liquefaction evaluation.

4.2.1.   How to calculate energy flow

First, the level ground is idealized by a model of horizontal soil layers where the SH wave propagates vertically as shown in Figure 30. Then, it is necessary to evaluate the energy flow by separating a recorded subsurface record into upward and downward waves. Let E_u,_m and E_d,_m denote the upward and downward energies at the upper boundary of the m-th layer, and corresponding energies at the lower boundary of the (m−1)-th layer as E′_u,_m−1, E′_d,_m−1, respectively. Considering the energy balance at the boundary between the m-th and (m−1)-th layer, the next equation holds, where E_T means the total energy passing through the boundary during a given earthquake.

From this, the following equation is readily derived.

Here, E_w is the energy dissipated in soil layers above the layer boundary during an earthquake, because all the energy computed here is assumed to transmit vertically.

Based on the multiple reflection theory, the upward and downward SH waves and corresponding wave energies at arbitrary levels can be evaluated from a single record at a given level using the condition that the ground surface is free from shear stress (e.g., Schnabel et al., 1972 ^[46]). If downhole records are available, however, they will considerably improve the energy flow evaluation, which may not fully comply with the simple 1D SH-wave propagation theory. Suppose seismic records at two subsurface levels, B and C, are available in addition to surface A, as illustrated in Figure 30. The energy flow between B and C can be calculated from earthquake records at the two levels ^[24]. Between the ground surface, Point A, and Point B, conversely, two sets of energy flow can be calculated using the earthquake record either at A or B. The two sets of energy may then be averaged with the weight of relative proximity to A and B to have plausible values (Kokusho and Suzuki 2011 ^[25]).

Nine earthquakes (Kobe EQ. & EQ1 to EQ8) are addressed here with the moment magnitude M_W = 6.6 to 7.9 or JMA magnitude (magnitude in Japan Meteorological Agency scale, similar to Richter scale) M_J = 6.7 to 8.0. In total, 30 vertical array sites were selected during the nine earthquakes with focal distances ranging from 9 to 227 km. The deepest depth of the vertical arrays of three-dimensional accelerometers spanned from 83 to 260 m, though most of them were nearly 100 m. The S-wave velocities at the base were widely diverged as V_s = 380~2800 m/s due to geological differences, while the surface velocities were mostly V_s = 90~430 m/s. Four vertical array sites for the 1995 Kobe EQ. and one site for EQ7 consist of accelerometers at three or more different levels, including the ground surface, while all others belonging to the KiK-net consist of only two levels, a surface and base ^[25].

The scalar sum of the wave energies calculated from the two orthogonal horizontal acceleration records was used for the energy flow evaluations. Equivalent linear soil properties, S-wave velocity V_s and damping ratio D both assumed as non-viscous or frequency independent (Ishihara 1996 ^[47]), were optimized for the main shock records and incorporated in the evaluations.

4.2.2.   General trends of energy flow in vertical arrays

Figure 31 depicts the variations of upward energy E_u along the depth z calculated for the nine earthquakes at 30 vertical array sites ^[25]. Considering the large differences in the absolute energies among the records, the horizontal axis is taken in log scale. The energy E_u shows an obvious decreasing trend in most sites with decreasing depth regardless of the energy magnitude. In some sites, it tends to decrease to less than one-tenth from the base to the surface, and the trend is more pronounced near the surface in contrast to the depth of 50~100 m or below.

There used to be a traditional view in engineering seismology that the wave energy ( = square of velocity amplitude × wave impedance) is kept constant as seismic waves propagate underground (e.g., Joyner and Fumal 1984 ^[48]). Hence, the velocity amplitude was generally considered inversely proportional to the square root of the impedance ρV_s. Figure 31 indicates that this view is not true, at least in shallow depths at most sites despite a few exceptions (KNK and KKNPS-SH).

Closer scrutiny of Figure 31 has revealed (^[25], Kokusho 2022 ^[49]) that the major cause of the drastic decrease of upward energy with decreasing depth is layer boundaries of clear impedance ratio near the ground surface, where a large amount of upward wave energy is reflected downward to return to the earth before arriving at surface soils. Also revealed was that internal energy dissipation is another cause in those sites where soft soils are prevalent near the surface.

4.3.   Empirical formulas for upward energy

4.3.1.   Upward energy ratio versus impedance ratio

In order to evaluate the decreasing rate of the upward energy as it approaches the ground surface, an empirical formula has been developed, wherein ratios of upward energies between layers are correlated to corresponding impedance ratios using the dataset of vertical array records addressed above ^[26]. Among the depth-dependent upward energy variations, 23 sites are used where the mutual difference is within 25% between the two upward energies at the deepest level calculated from measured motions at the ground surface and the deepest level.

Impedance ratios α and upward energy ratios β defined between two neighboring layers, m and m+1, in a given soil profile in Figure 30 are calculated at individual sites from surface to base as;

The soil density was assumed depending on the S-wave velocity as: ρ = 1.6~2.0 t/m³ for V_s≤300 m/s, ρ = 2.0~2.2 t/m³ for 300 m/s ≤ V_s ≤ 700 m/s, ρ = 2.3~2.4 t/m³ for 700 m/s ≤ V_s≤ 1000 m/s, and ρ = 2.5~2.7 t/m³ for 1000 m/s ≤ V_s < 3000 m/s. The energy ratios β are plotted versus the corresponding impedance ratios α in Figure 32 (a) for all layers above the deepest levels with different symbols in the 23 vertical array sites. For most data points, α ≤ 1.0 holds because the impedance ratio is normally less than unity (ρV_s is getting larger in deeper layers). It is quite reasonable to assume that β = 0 for α = 0 at a ground surface, and β = 1 for α = 1 within a uniform base layer. Hence, a next simple power function may be practically used for 0 ≤ α ≤ 1.0 to approximate the plots as the thick solid curve in Figure 32 (a), and the power n = 0.70 is obtained from the least mean square method with the determination coefficient R² = 0.81.

Furthermore, the impedance ratio α and the upward energy ratio β may be redefined differently from Eq. (25), between an arbitrary m-th layer and the deepest base layer in a vertical array as follows.

In Figure 32 (b), data points for all the layers at the 23 vertical array sites are plotted on the α versus β diagram following the definition in Eq. (27), wherein the symbols are connected with thin dashed lines for individual sites and differentiated by four steps of V_s-values at the base layers. Though the plots are more dispersed than in Figure 32 (a), the curve by Eq. (26), using α and β redefined in Eq. (27) seems to approximate the plots to a certain degree. As shown in the chart (b), some of the V_s-values of base layers are very high, 2400–3000 m/s, almost equivalent to seismologically defined stiff bedrocks. This indicates that Eq. (26) may be used for practical purposes to evaluate the upward energy in a shallow soil layer directly from the incident energy at a deep seismological bedrock using the impedance ratio defined between the two in Eq. (27).

4.3.2.   Upward energy at seismological bedrock

Upward energies at the deepest levels (base layers) denoted here as (E_u)_base calculated from the downhole records and the associated impedance values are plotted with close symbols in the full-logarithmic scale versus hypocenter distances R (m) in Figure 33 at all 30 vertical array sites for nine earthquakes ^[25,26]. Although the plots are widely scattered, the decreasing trends in (E_u)_base with increasing R for individual earthquakes are not difficult to identify. Among the nine earthquakes, the plots of EQ3 with M_J = 8.0 (2003 Tokachi-Oki earthquake) are relatively higher on the upper-right side of the diagram reasonably as a larger magnitude subduction earthquake, while other crustal earthquakes with M_J around 7.0 are lower on the lower-left side.

Based on the finding that Eq. (26) may be conveniently used to roughly evaluate the energy ratio between arbitrary two layers, the same equation is further used here to estimate the incident or upward energies at the seismological bedrocks (E_u)_sbr from those at the base layers (E_u)_base of individual vertical arrays. The impedance for the seismological bedrock is calculated thereby using V_s = 3000 m/s and ρ = 2.7 t/m³. It may well be justified here that major wave energy is carried by vertically propagating SH waves even in seismological bedrocks as stiff as V_s = 3000m/s, as substantiated by some of the vertical array data in Figure 32 (b), despite the potential involvement of SV waves.

In Figure 33, the upward energies at the calculated seismological bedrock (E_u)_sbr are also plotted with open symbols for the nine earthquakes. They are positioned higher than the corresponding close plots because of higher ρV_s-values than those at the bottom of the vertical arrays, reflecting the upward energy decreasing as it goes up from the seismological to engineering bedrock. Straight lines drawn in Figure 33 represent the following formula of incident energy E_IP (kJ/m²) versus hypocenter distance R (m) by assuming the spherical energy radiation from the center of energy release, assumed here as the hypocenter (e.g. ^[6]).

Total released wave energy E_Total in kJ is calculated by the next empirical formula ^[6] from earthquake magnitudes M,

Concerning magnitude M, it should be reminded that in the original paper ^[6] the surface wave magnitude M_S was employed as M in Eq. (29). According to Utsu (1982) ^[50], M_S may be similar to M_W (Moment Magnitude) and M_J (Japanese Meteorological Agency Magnitude) as well with a small difference. However, the actual magnitude values of the nine earthquakes (Kobe EQ, EQ1~EQ8), which were announced by JMA in Japan and USGS in the USA, are all distinctively larger for M_J than the corresponding M_W. We also found that the incident wave energies estimated by the empirical formula Eqs. (28) and (29) using M = M_J tend to be more compatible than those by M = M_W with the upward energies (E_u)_sbr calculated at the seismological bedrock using Eqs. (26) and (27) ^[49]. Thus, the JMA magnitude may be suitable to estimate the incident energy E_IP at the seismological bedrock for engineering purposes. Certainly, more detailed studies will be needed to upgrade the empirical equations by incorporating fault and path mechanisms of individual earthquakes.

4.4.   Energy flow in wave propagation

One of the highlights of the energy-based liquefaction evaluation developed here is to grasp the total amount of demand energy and directly compare it with the capacity to determine the extent of liquefaction at a site. The upward energy going up through a soil profile will reflect at the ground surface and come down again as the downward energy except that dissipated in the upper layers ^[43]. Thus, the upward energy represents all the demand energy to be compared with the capacity energy in the EBM, because the downward energy is originally a part of the upward energy. In this regard, how to compare the demand energy of design earthquake waves with the capacity or dissipated energy for liquefaction considering wave attenuation and surface reflection of the SH-wave is discussed below.

Let us consider a harmonic SH-wave shown in Figure 34 (a) propagating in the z-axis (upward direction) with time t in a viscoelastic medium as;

Here, B = wave amplitude, $\omega$ = 2πf = angular frequency, and V_s^* = complex S-wave velocity considering non-viscous damping (Ishihara 1996 ^[47]) which can be written as;

where $G + iG'$ = complex shear modulus with real and imaginary parts and ρ = soil density. δ is the phase delay angle of the stress-strain hysteresis obtained in cyclic loading of soils, using G and G' and also correlated with damping ratio D (e.g. ^[37,47]) as;

Eq. (30) is also written in the following form;

where modified S-wave velocity is defined as

which can be approximated as ${V'_s} \approx V_s^*$ for δ < < 1.0. The wave attenuation coefficient β is defined ^[37,47] as;

This can be approximated for δ < < 1.0 as

If Eq. (33) is substituted into Eq. (21), and integrated for one period ${{T = 2\pi } \mathord{\left/ {\vphantom {{T = 2\pi } \omega }} \right. } \omega }$ of the harmonic wave, the energy in one wavelength $\lambda = {{2\pi {V_s}} \mathord{\left/ {\vphantom {{2\pi {V_s}} \omega }} \right. } \omega }$ can be obtained as;

Let the amplitude of particle velocity be expressed as

and the shear strain amplitude is assumed $\delta < < 1.0$ as

Then, the upward wave energy in Eq. (37) can be expressed as

Hence, the energy density per unit volume in the upward wave is

Here, W is equal to the maximum strain energy defined in a cyclic loading test of an ideally viscoelastic material of shear modulus $G = \rho {V_s}^2$ with the amplitude ${\gamma _a}$ shown in Figure 34 (b).

The wave energy E' arriving at z+z₀ shown in Figure 35 (a) is formulated using the energy E at z in Eq. (40) as;

Hence, the energy dissipated in travelling from z to z+z₀ in one wavelength z₀ = λ is;

Thus, the ratio of dissipated energy to the total energy in one wavelength z₀ = λ is written as;

It can be approximated by using the Taylor expansion for $\beta = {{\omega D} \mathord{\left/ {\vphantom {{\omega D} {{V_s}}}} \right. } {{V_s}}} < < 1.0$ or D < < 1.0 as;

This relationship (Kokusho 2017) ^[37] resembles the next equation which correlates the dissipated energy ΔW with the elastic strain energy W using the damping ratio D in cyclic loading tests as indicated in Figure 34 (b).

Hence, the same mechanism is supposed to govern the wave propagation as in the laboratory cyclic loading test, wherein E and ΔE correspond to W and ΔW, respectively. If the dissipated energy ΔW is negligibly small, strain energy W given to the system in the first half cycle is all recovered and used in the subsequent loading. It is thus clear that the wave energy dissipates with the ratio defined in Eq. (45) for small D during wave propagation.

4.5.   Energy demand near free surface

The upward SH-wave is reflected at a ground surface and becomes a downward SH-wave as illustrated in Figure 35 (a). In this condition, the displacement of the harmonic wave with angular frequency $\omega$ may be formulated as;

where the z-axis is taken downward to be positive, and A and B are amplitudes of downward and upward waves, respectively. The internal damping is assumed zero in this case for simplicity. Due to zero shear stress τ = 0 at the free ground surface,

or A = B, indicating that the wave reflection occurs with the surface displacement $u = 2B{e^{i\omega t}}$. Accordingly, the depth-dependent variations of displacement u, particle velocity $\dot u = du/dt$ and shear strain γ = du/dz during steady-state vibration by the harmonic wave are written respectively as;

The displacement in Eq. (49) illustrated with a dashed curve in Figure 35 (a) indicates that nodes and antinodes of vibration, ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$ apart in between, appear regularly at every ${\lambda \mathord{\left/ {\vphantom {\lambda 2}} \right. } 2}$, starting from the top antinode at the surface. The kinetic and strain energies per unit volume E_k/λ and E_γ/λ averaged over the period $T = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } \omega }} \right. } \omega }$ can be obtained as Eqs. (52) and (53) ^[31].

The sum of the two energy densities becomes quite reasonably a constant twice larger than the upward energy density as;

As shown in Figure 35 (a), the strain energy normalized as (E_γ/λ)/2W is 0% at the antinodes and 100% at the nodes, while the kinetic energy (E_k/λ)/2W is vice versa. Thus, the two kinds of energy are distributed with a fixed rate of 0 to 100% depending on the depth and are not convertible to each other because of the constraints of the nodes and antinodes appearing at every quarter wavelength λ/4 for the steady-state harmonic vibration. This is a significant difference from the one-way wave propagation without any boundary in Figure 34 (a), where both the kinetic and strain energies are eligible to compensate for the dissipated energy.

Figure 35 (a) is, however, an extreme case by a steady harmonic motion of a single frequency and may not represent a realistic seismic response by transient irregular waves, wherein nodes and antinodes are difficult to appear steadily. Nevertheless, even in irregular seismic motions, the free ground surface consistently serves as an antinode for waves of any frequencies because of the stress-free condition. Such a surface boundary effect may well be assumed to fade away with increasing depth from the surface, difficult to transmit beyond a certain depth.

In Figure 35 (b), where the shallow depth of ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$ from the surface is zoomed in, the normalized strain energy (E_γ/λ)/2W is depicted using Eq. (53), indicating that the ratio of strain energy E_γ imparted from the upward energy E_u to compensate for the liquefaction-induced dissipated energy is increasing from zero at the ground surface to the maximum value 1.0 at z = ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$, and the ratio is fixed depending on z during harmonic vibration as;

For earthquake motions, however, the effect of the surface boundary may be limited within a shallow depth because of the irregularity of motions and the nonlinear soil properties during strong earthquakes, too. Hence, it is assumed here that, beyond the depth of a quarter wave-length z = ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$ for λ = V_sT corresponding to predominant periods T of earthquake motions, the SH wave propagates up and down as in Figure 34 (a) without any surface boundary effect where E_k and E_γ can freely compensate for the dissipated energy. In other words, if E_u^* stands for the wave energy that can compensate for the dissipated energy, and ${E_u}$ is the total wave energy available at the same depth z, then:

In between the surface and z = ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$, the boundary effect gradually weakens with increasing depth from E_γ = 0 or E_u* = 0 at z = 0 to Eq. (56) at z = ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$. Hence, it may be possible to take the average of Eq. (55) in between considering the effects of irregular earthquake motions as well as soil nonlinearity and simplify the energy ratio as the next equation for practical purposes.

Considering that the wavelength $\lambda = {V_s}T$ may be 80 ~ 160 m in most liquefiable site conditions with average shallow S-wave velocity V_s > 160 m/s for z < 20 m and predominant periods of earthquake motions T > 0.5 s, the depth of ${\lambda \mathord{\left/ {\vphantom {\lambda 4}} \right. } 4}$ means around 20 m or deeper, indicating that for normal liquefaction evaluation depths z < 20 m, the upward energy may be halved as in Eq. (57) to compare with the liquefaction energy capacity ^[31].

To examine how well the simplification mentioned above holds, we conducted an additional numerical study on irregular earthquake motions. A uniform horizontal soil layer of 20 m thick with ρ = 1.8 t/m³, V_s = 160 m/s, and the damping ratio D = 0, 20 & 40% is vibrated by surface input motions as illustrated in Figure 36 (a). Four acceleration records EQ1 to EQ4 listed in Figure 36 (b) with different durations and irregularities obtained at soft soil sites in Japan have been used. Their predominant periods read off in velocity response spectra in Figure 36 (c) are varied from T = 0.26 s to 2.5 s. Induced strains due to upward and downward waves, ${\gamma _u}\left({t, z} \right)$ and ${\gamma _d}\left({t, z} \right)$, have been calculated at different depths z in the 1D wave propagation analysis, respectively, from which cumulative strain energies E_γ (z) are determined by the next equation derived from Eqs. (17) and (19).

The energy E_γ thus calculated is plotted versus depth z with open symbols in Figure 37 to compare with the demand energy of close symbols of constant E_u along the depth z for the case of D = 0%. If D > 0%, the demand energy is the average of upward and downward energies, (E_u+E_d)/2 along z considering wave attenuation. Here, a quarter wavelength λ/4 calculated from the predominant period T as λ/4 = V_sT/4 for V_s = 160 m/s is indicated with a horizontal dash line for each earthquake. The quarter wavelength λ/4 tends to decrease along with equivalent linear V_s lowering with increasing soil nonlinearity. Furthermore, the damping ratio will also increase as D = 0%, 20%, and 40 %, for example, which shows the effect of damping on the strain energy E_γ, appearing insignificant in shallower depths as illustrated.

According to Figure 37, the depth-dependent increase of E_γ seems variable among the four earthquakes: First, due to the difference in the characteristic wavelength; second, due to the wave irregularity. Nonetheless, the strain energy E_γ to compensate for the dissipated energy has been confirmed to start from zero at the ground surface, implying that liquefaction is difficult to occur near the surface as long as the SH-wave field is assumed. With increasing depth, E_γ tends to be nearly half of (E_u+E_d)/2 for the same D-value at the same depth but never goes far beyond.

Thus, as assumed for the harmonic motion in Figure 35 (b), the strain energy E_γ averaged over the depth of a quarter wavelength λ/4 cannot be larger than 1/2 of the demand energy for strong earthquake motions causing soil nonlinearity. Hence, the simplification by Eq. (57) may be justified in determining the energy ${E_u}^*$ compensating for liquefaction-induced dissipated energy within the depth of a quarter-wavelength or any depth for normal liquefaction evaluations.

4.6.   How to compare upward energy with dissipated energy

As mentioned, the key of the present EBM is that the upward wave energy is quantified and directly compared with the dissipated energy for liquefaction. The comparison should be made so that liquefaction occurs when the cumulative wave energy in the field attains the same amount of dissipated energy to liquefy the soil in laboratory cyclic loading tests.

In the laboratory test of a viscoelastic material in Figure 38 (a), shear stress $\tau = {\tau _a}\sin \omega t$ induces strain $\gamma = {\gamma _a}\sin \left({\omega t - \delta } \right)$ with a phase-delay angle $\delta$. The dissipated energy ΔW is formulated from the area of the τ~γ hysteresis loop ACDA'C'D'A as

Thus, the angle δ is a key parameter determining the dissipated energy during cyclic loading. In the same loading, twice the elastic strain energy W = Area (OAB) = Area (OA'B') is given to the material. Of the energy W = Area (OAB) in the first 1/2 cycle, the energy corresponding to Area (ABC) can be recovered and recycled in the second 1/2 cycle filling Area (A'B'C'). The Area (ABC) is calculated as

Hence, the energy supply in one cycle loading denoted here as 2W₋ considering the energy recycling by the Area (ABC) is obtained as ^[23,37];

Then, the ratio of the dissipated energy ΔW to the supplied energy 2W₋ in one cycle is written as;

In Figure 38 (b), the variation of the energy ratio ΔW/2W₋ in Eq. (62) is shown versus the damping ratio D and compared with ${{\Delta E} \mathord{\left/ {\vphantom {{\Delta E} E}} \right. } E} = 1 - {e^{ - 4\pi D}}$ in Eq. (45). The two curves on the chart are very similar to each other, both have the same initial tangent and tend to approach to the asymptote ΔW/2W₋ = ΔE/E = 1.0 with increasing D ^[23]. This indicates that the energy dissipation mechanism during wave propagation is very similar for not only infinitely small but also large damping and reproduced by laboratory cyclic loading tests despite a small gap of a maximum of 10%. The gap may be attributed to the difference in loading mechanism; simultaneous cyclic loading on a whole soil specimen versus time-delayed loading on in situ soils during wave propagation ^[23].

In conventional data reduction of cyclic loading tests, the elastic strain energy W = Area (OAB) is normally employed to divide the dissipated energy ΔW, as ΔW/W. Figure 36 (b) shows that the wave energy ratio ${{\Delta E} \mathord{\left/ {\vphantom {{\Delta E} E}} \right. } E} = 1 - {e^{ - 4\pi D}}$ is represented very well by ${{\Delta W} \mathord{\left/ {\vphantom {{\Delta W} W}} \right. } W} = 4\pi D$ when D is less than a few percent. This is because of the energy 2W supplied in one cycle, W is recycled without being dissipated on account of small damping D. With increasing D, ${{\Delta W} \mathord{\left/ {\vphantom {{\Delta W} W}} \right. } W} = 2\pi D$ may become better to approximate ΔE/E (e.g. Kazama et al. ^[17]) though there is a distinctive difference from ΔE/E. In contrast, Eq. (62) seems to successfully simulate the energy dissipation ratio during wave propagation for wide D-values as shown in Figure 38 (b).

The above formulation in one-cycle loading may be extended to a similar relationship for cumulative energies in liquefaction tests by multi-cycle loading if the damping ratio D can be approximated as a constant averaged over the whole developing process of liquefaction. We will demonstrate that the damping ratio stays D = 0.1 to 0.2 with the average D = 0.15 during liquefaction tests of saturated sands in Figure 39. Hence, the following equation may be valid in place of Eq. (62).

Instead of the above equation where the soil is assumed ideally viscoelastic, the value ∑ΔW/∑2W₋ may be calculated on actual soil materials using laboratory cyclic loading test data already addressed. Figure 39 (a) exemplifies a typical stress-strain relationship in the undrained cyclic triaxial test on saturated Futtsu sand. Namely, using the dissipated energy per cycle, ΔW = Area (ABCDA), and the energy corresponding to Area (ABB'CDD'A) as well, the strain energy 2W₋ supplied per cycle considering the recycling effect of elastic energy can be formulated as;

In Figures 39 (b) and (c), the cumulative energies Σ2W_— thus calculated are plotted versus the cumulative dissipated energy ΣΔW (^[31], JGS Committee Report 2018 ^[51]) from the 1^st to the last cycle in triaxial and torsional liquefaction tests, respectively. Note that the gradient of the plots, (Σ2W₋/σ'_c)/(ΣΔW/σ'_c), evaluates very similarly among the test results of different conditions; 1.03 for the reconstituted Futtsu sand, 1.07~1.34 for intact soils both tested in triaxial tests, and 1.10~1.23 for the Futtsu sand tested in torsional tests. Also noted is the evident proportionality of Σ2W₋ to ΣΔW in all the plots, indicating that the D-value in Eq. (63) is almost constant throughout the loading cycles as already assumed. This implies that the energy dissipation mechanism in soils will not significantly vary with the development of liquefaction.

The finding that the gradients (Σ2W₋/σ'_c)/(ΣΔW/σ'_c) read off from all the plots are near unity can lead to a simple result as ΔW = 2W_−, so that the soils tend to dissipate almost all the energy supplied cycle by cycle throughout the liquefaction process during cyclic loading tests. Hence, it may be justified to make use of such simplification in the EBM of liquefaction evaluation as;

Thus, ${{\sum {\Delta W} } \mathord{\left/ {\vphantom {{\sum {\Delta W} } {\sum {2{W_ - }} }}} \right. } {\sum {2{W_ - }} }} \approx 1$ in Eq. (65) and ${{{E_u}^*} \mathord{\left/ {\vphantom {{{E_u}^*} {{E_u}}}} \right. } {{E_u}}} = 1/2$ in Eq. (57) are two essential formulas on energy ratios in the EBM of liquefaction evaluation developed here. Because the dissipated energy for liquefaction in the laboratory, ∑ΔW, should be compared with the wave energy density ${{{E_u}^*} \mathord{\left/ {\vphantom {{{E_u}^*} \lambda }} \right. } \lambda }$ available in situ, ∑2W₋ has to be compared with E_u/(2λ) in evaluating in situ liquefaction behavior. Consequently, E_u is to be compared with 2 × ∑2W₋ as already stated in the previous paper ^[18,23], which is denoted as ∑W₋^* ≡ 2 × ∑2W₋ and also used in this paper hereafter.

5.   Evaluation steps and examples of EBM

After visiting basic research findings on the energy-based method (EBM) for liquefaction evaluation in terms of capacity, demand, and their mutual comparison, it seems necessary to show practical steps to implement the EBM in engineering practice. Thus, the evaluation steps are as simplified and practical as possible, reliable, and presented and followed by some examples. Some of the equations and methods used here have been updated from previous publications by the author and his coworkers reflecting more recent research findings.

The evaluation steps of the current EBM ^[31] quantifying the upward wave energy as the energy demand are illustrated in Figure 40, composed of two stages. Stage-Ⅰ comprises Step-1 to Step-6 to make a preliminary evaluation of the liquefaction potentials of individual layers. Stage-Ⅱ comprising Step-7 to Step-10 is to make a higher-level liquefaction prediction including induced shear strains and settlements based on the results of Stage-Ⅰ.

5.1.   Stage-Ⅰ: Preliminary Evaluation Steps for Initial Liquefaction

Step-1:

At a given site, a soil profile is divided into layers of an identical thickness H = 1 to 2 m throughout the soil model following in situ penetration test data available. The CRR- values for an equivalent number of cycles N_L for the initial liquefaction in SBM, typically CRR₁₅ for N_L = 15, are determined in individual layers based on the penetration resistances or other soil test data.

Step-2:

Normalized cumulative dissipated energy (abbreviated as dissipated energy hereafter) $\sum {{{\Delta W} \mathord{\left/ {\vphantom {{\Delta W} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}}$ for soil to attain the initial liquefaction (corresponding to maximum double amplitude axial strain ${\varepsilon _{DA\max }}$ = 5% in triaxial tests or maximum double amplitude shear strain ${\gamma _{DA\max }}$ = 7.5% in simple shear tests) is determined in each layer from CRR-value; such as CRR₁₅ for N_L = 15;

Here, the dimension of ΣΔW is energy per unit volume, and non-dimensionalized by effective confining stress ${\sigma '_c}$. Eq. (66) is identical to Eq. (5) developed from a series of triaxial liquefaction tests on intact soils of various properties. Thus, the correlation is available to uniquely determine the dissipated energy ΣΔW/σ'_c corresponding to the initial liquefaction (${\varepsilon _{DA\max }}$ = 5% or ${\gamma _{DA\max }}$ = 7.5%) regardless of earthquake motions from CRR in SBM for widely varying in situ soils irrespective of soil types ^[20].

Step-3:

The strain energy ∑2W₋ considering the energy-recycling effect is identical to the dissipated energy ∑ΔW all through the liquefaction development as already given in Eq. (65) based on the cyclic loading test results shown in Figure 39. Because the upward energy compensating for the dissipated energy is approximated as half of the upward wave energy E_u in surface soils down to a quarter wavelength λ/4 (normally larger than 20m) as in Eq. (57), ∑2W₋/σ'_c is doubled in advance and denoted as $\sum {{{W_ - ^*} \mathord{\left/ {\vphantom {{W_ - ^*} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}}$ to directly compare with the demand energy E_u later on in STEP-6.

Step-4:

Energy capacity (energy per area in dimension) in a layer of thickness H is calculated as ΣW_—^*H from Eq. (67) using ΣΔW/σ'_c corresponding to the initial liquefaction, where σ'_c = (1+2K₀) σ'_v/3 is effective confining stress, K₀ = earth-pressure coefficient at rest, and σ'_v = effective overburden.

Step-5:

The upward wave energy E_u (energy per area in dimension) is calculated from an equivalent linear one-dimensional SH-wave propagation analysis in a soil profile at a given site using $\dot u$ = particle velocity of upward wave, ρV_s = impedance of SH wave wherein the degradation of V_s due to soil nonlinearity is reflected, and energy at the end of motion E_uf is determined in each layer by integrating Eq. (68) or Eq. (21) in terms of time t.

For those sites where a design earthquake motion is unavailable, the energy E_uf may be roughly estimated by incorporating empirical formulas, Eqs. (26) to (29) as delineated. From earthquake magnitude M and focal distance R (unit: m) as well as the impedance ratio α = (ρV_s)_m/(ρV_s)_sbr between m-th layer and a seismological bedrock, the upward energy at a layer m, (E_u)_m (unit: kJ/m²), is determined using;

where the upward energy at the seismological bedrock (E_u)_sbr (unit: kJ/m²) is computed as;

If the upward energy is estimated using these empirical equations, note that the ultimate upward energy E_uf equivalent to Eq. (68) at m-th layer should be halved as E_uf = (E_u)_m/2 to compare with the dissipated energy ∑ΔW, because (E_u)_sbr in Eq. (69) is evaluated as the sum of two horizontal directions while ∑ΔW is normally determined from soil tests by one-directional shearing.

Step-6:

The energy ratio ΣW_—^*H/E_uf is calculated in each layer using the values ΣW_—^*H and E_uf above. A layer with a smaller energy ratio has higher and earlier liquefaction potential than others among a given soil profile ^[18]. The energy ratios of individual layers are numbered sequentially starting from the lowest ratio (j = 1) toward the higher ones and summed up in terms of j, which is named AER (Accumulated Energy Ratio).

Liquefaction occurs in that sequence exclusively in those layers of AER$\leqslant$1.0, because the upward energy can liquefy individual layers in the above sequence until it is completely dissipated by liquefaction. Thus, unlike the conventional Stress-Based Method, the EBM has the potential to decide the extent of liquefaction to be within the availability of the total energy demand. To make the most of this potential, the steps in Stage-Ⅱ are taken for a higher-level liquefaction evaluation as below.

5.2.   Stage-Ⅱ: Steps for Strain and Settlement Evaluation

As mentioned, a remarkable feature of the present EBM is that the total energy demand for liquefaction can be prescribed at a site if a design earthquake is given. This will enable one to evaluate not only if liquefaction occurs or not but also if it occurs how large liquefaction-induced strain and settlement will develop.

This is possible in a simple case, where only a single layer has the potential to liquefy, to evaluate the induced strain by demand energy exclusively given there. In the case of multiple liquefiable layers, it may also become possible to calculate the strains if a rule is adopted on how the upward energy is shared among the layers. However, such a universal rule is never easy to establish because the energy distribution seems to be dependent on soil profiles, input seismic motions, and nonlinear soil properties. Sophisticated nonlinear effective stress analyses based on rigorous constitutive laws may provide the solution, though that is out of the scope of simplified and efficient liquefaction evaluations investigated here.

Instead of postulating the same induced strain of γ_DA = 7.5% corresponding to the initial liquefaction in all layers in Stage-Ⅰ, a simple rule may be introduced in Stage-Ⅱ that the wave energy is equally shared among those layers where liquefaction has been judged to occur (AER < 1.0) in Stage-Ⅰ in Eq. (71) ^[31]. It may well be expected that this rule despite its simplicity will lead the evaluation not exactly but closer to actual behavior than Stage-Ⅰ. Let us assume for example the number of liquefied layers in Stage-Ⅰ as m, and count other sequential numbers for Stage-Ⅱ from the top to the bottom of the soil layers, k = 1~m, as indicated in Figure 40.

Maximum induced DA shear strain in one of the liquefiable layers k may be calculated assuming the proportionality to dissipated energy as;

Here, (ΣW₋^*H)_k = the energy capacity (determined by dissipated energy for the initial liquefaction) corresponding to γ_DAmax = 7.5%, and (E_uf)_k/m = the energy demand (upward energy) allocated to the k-th layer among the m liquefiable layers. This equation stems from experimental observations such as in Figure 25 (b) or Figure 4 in cyclic torsional or triaxial tests using harmonic and earthquake motions ^[18,31] despite the data scatters between double amplitude shear strain γ_DAmax and dissipated energy ΣΔW/σ_c'. Thus, the uniqueness of energy in terms of induced strain seems valid regardless of different loading histories as a great merit of the energy-based evaluation.

The maximum strain γ_DAmax in Eq. (72) is accepted to be closely correlated with liquefaction-induced volumetric strain ε_v by many investigators (Silver and Seed (1971) ^[52], Tatsuoka et al. ^[53], Nagase and Ishihara ^[54], Tsukamoto et al. ^[55]). Quite a lot laboratory test data are available for the volumetric strain by undrained cyclic loading. Among them, Ishihara and Yoshimine ^[56] simplified their data as a bilinear relationship between volume change ε_v and single or double amplitude maximum shear strain, γ_max or γ_DAmax = 2γ_max, schematically shown in Figure 41 (a), where ε_v is proportional with γ_DAmax up to the upper limit ε_vmax corresponding to γ_DAmax = 20% (γ_max = 10%), and the volumetric strain in the vertical axis is normalized by the corresponding upper limit as ε_v/ε_vmax.

To evolve the ${\varepsilon _v}$ versus ${\gamma _{DA\max }}$ correlations to be more robust and applicable to in situ soils of various physical properties, the database ^[34] of in situ intact soils disseminated by PWRI is utilized. This consists of numerous cyclic triaxial liquefaction tests of intact soils sampled from various natural and manmade ground in liquefied areas during the 2011 Tohoku earthquake in Japan together with pertinent in situ test data mainly by Standard Penetration Test including N₁ (corrected blow counts for effective overburden 98 kPa), F_c (fines content), and G_c (gravel content). A multiple regression analysis for the volumetric strain ε_vmax corresponding to the maximum shear strain γ_DAmax = 20% by utilizing the database has yielded the following formula, though the determination coefficient R² = 0.458 is not high ^[31].

In Figure 41 (b), the values ε_vmax on the vertical axis are plotted versus N₁ and F_c horizontally in the three-dimensional diagram. The plots are projected on the two-dimensional plane of ε_vmax versus N₁, and compared with the star plots connected with solid lines for clean sand of F_c = 0, which are from a design chart ^[56] based on reconstituted clean sands for γ_max = 10% (γ_DAmax = 20%). The two plots essentially share a similar trend of ε_vmax decreasing with increasing N₁-value, wherein the widely varying fines content in the natural soils (F_c = 1~97%) obviously has a great influence.

Based on the above test data and interpretations, the following steps to evaluate liquefaction-induced strains and settlements can be developed as Stage-Ⅱ of the energy-based liquefaction evaluation procedure.

Step-7:

The upward energy at a base layer is divided by m (the number of liquefied layers in Stage-Ⅰ), and the energy at a k-th layer (E_uf)_k /m comes up independently to layer k to be dissipated for liquefaction exclusively there.

Step-8:

The maximum induced shear strains assumed in Stage-Ⅰ as γ_DAmax = 7.5% and identical in all those layers that were judged liquefiable are now recalculated in Eq. (72) to yield individually different (γ_DAmax)_k-values caused by the allocated upward energy (E_uf)_k/m using the strain energy (Σ2W_—^*H)_k corresponding to the initial liquefaction. Consequently, the number of layers beyond the initial liquefaction (γ_DAmax = 7.5%) tends to be smaller because the equally-allocated upward energy tends to induce higher strain than γ_DAmax = 7.5% in weaker layers, and lower strain in stronger layers. Thus, the higher-level liquefaction evaluation in Stage-Ⅱ can be implemented differently from Stage-Ⅰ, which may be more plausible in actual liquefaction performance.

Step-9:

Using the induced shear strains in individual liquefiable layers determined in the previous step, the corresponding volumetric strains ε_v are calculated in proportion to γ_DAmax if γ_DAmax ≤ 20% as;

or a constant if γ_DAmax > 20%.

where (ε_vmax)_k is determined in Eq. (73) in individual layers.

Step-10:

One-dimensional ground surface settlement S can be obtained as the sum of settlements in individual layers k = 1~m involved in liquefaction as;

Here, not only those beyond the initial liquefaction (γ_DAmax > 7.5%) but also those before are involved in calculating the settlement as long as they are judged liquefiable in Stage-Ⅰ.

5.3.   Examples & case histories

Four soil models; one hypothetical uniform sand deposit and three case history sites liquefied and are evaluated using the EBM steps delineated above. The results are compared with those of the conventional stress-based method (SBM). The ground surface settlements in the case history sites are compared with the observation if available during the earthquakes.

5.3.1.   Uniform sand model

A hypothetical uniform sand model 10 m thick with corrected SPT blow counts N₁ = 8 and F_c = 0%, underlain by a stiff base shown in Figure 42 (a) was studied ^[18]. The cyclic resistance ratio corresponding to N₁ = 8 and F_c = 0% for the equivalent number of cycles N_eq = 20 is given by a Japanese design code (Japan Road Bridge Association ^[57]) as CRR₂₀ = (γ_DAmax)_{triax.Neq = 20} = 0.191 determining the energy for the initial liquefaction ΣΔW/σ'_c = 0.0372 in Eq. (6). Due to the uniqueness of energy, CRR₁₅ can be derived from Eqs. (5) and (6) as;

CRR₁₅ = 0.204 thus calculated is used hereafter in the SBM because N_eq = 15 is more widely accepted internationally than N_eq = 20 corresponding to M = 7.5 earthquakes.

The soil conditions of the model are listed on the left side of Table 2. The model is sliced into 5 layers of H = 2 m thick, L1 to L5, with the top (L1) unsaturated and the rest saturated. The SPT N at each layer was calculated from N₁ = 8 using N₁ = 1.7N/(σ_v^'/p₀ + 0.7) ^[35] where σ_v^' = effective overburden, p₀ = unit pressure, and the S-wave velocity V_s (m/s) was determined from the N-value using an empirical formula V_s = 80 N^{1/3 [57]}.

The dynamic ground response was calculated by an equivalent linear 1D response program ^[45] using the hyperbolic model (Hardin & Drnevich 1972 ^[58]) modified for better fitting to laboratory data ^[18]. The acceleration motion in Figure 42 (b) recorded at K-NET Urayasu (EW-direction) during the 2011 Tohoku earthquake in Japan (offshore subduction event of M_W = 9.0) was given to the model surface either in the real-time scale (RT: Duration 236 s) or in a compressed half time scale (RT/2: Duration 118 s) with the same acceleration amplitude. Note that the cumulative upward energy E_u and its ultimate value at the end of earthquake E_uf calculated by Eq. (68) differs tremendously (8 to 1 ratio) between RT and RT/2 despite the same absolute accelerations as depicted in the bottom of (b).

In Table 2, Stage-Ⅰ evaluations of EBM are shown together with Stage-Ⅱ for the RT and RT/2-motions, where layers judged to liquefy are shaded in the columns (note that summations Σ for cumulative energy by cyclic loading in EBM are all abbreviated in the following tables). From the dissipated energy for liquefaction onset ΣΔW/σ'_c = 0.0372, the strain energy per unit volume is given as ΣW_—^*/σ'_c = 0.0744 from Eq. (67). The liquefaction energy capacities WH for H = 2 m thick layers to liquefy are calculated using the average confining stresses σ'_c = σ'_v (1+2K₀)/3, assuming the lateral earth-pressure coefficient at rest K₀ = 0.5. Energy ratio ΣW_—^*H/E_uf is calculated in each layer using the values ΣW_—^*H and E_uf above, and summed up as AER = Σ_i (ΣW_—^*H/E_uf)_i in terms of i from the lowest ratio (i = 1) toward the higher ones in Eq. (71). In the RT-motion, all four saturated layers show AER-values below unity and hence liquefy in the EBM evaluation. In contrast, for the RT/2-motion, Σ_i (ΣW_—^*H/E_uf)_i < 1.0 only for L2, indicating the upward energy is insufficient to liquefy more than that. Thus, a clear difference exists between the liquefaction potentials of the two input motions, reflecting the tremendous energy reduction in the RT/2-motion as mentioned before.

The corresponding SBM evaluation results on the same model are also listed in Table 2 in terms of a safety factor F_s = CRR/CSR. So far, the terms "CSR" and "CRR" have been defined to represent stress ratios between cyclic shear stress τ_d = σ_d/2 and isotropic effective confining stress σ'_c considering laboratory soil tests on isotropically consolidated specimens as CSR = (σ_d/2σ'_c)_triax. and CRR = (σ_d/2σ'_c)_triax. for N_L, respectively. In calculating F_s in SBM, CRR has to be redefined as a ratio of cyclic shear stress to effective vertical stress CRR = σ_d/σ'_v and further written as CRR = 0.9 × (1 + 2K₀)/3 × (σ_d/2σ'_c)_triax., using the stress ratio obtained in isotropically consolidated cyclic triaxial tests (σ_d/2σ'_c)_triax., and K₀ = earth-pressure coefficient at rest, as well as a reduction constant 0.9 for horizontal two-directional shaking (Tokimatsu & Yoshimi 1983 ^[59]). The CSR is also redefined hereafter to represent the cyclic stress ratio of a given earthquake in the field CSR = (τ_d/σ'_v)_field = r_n × (τ_max/σ'_v) where τ_max = the maximum shear stress of an earthquake motion and σ'_v = effective overburden in the field. The index r_n = τ_d/τ_max is a stress reduction coefficient empirically proposed as r_n = 0.1 × (M−1.0) ^[59] to convert maximum shear stress τ_max to equivalent shear stress τ_d of a harmonic motion of the equivalent number of cycles N_L = 15 using earthquake magnitudes M. The index r_n is correlated with Magnitude Scaling Factor used in North American practice; MSF, as MSF = 0.65/[0.1(M-1)] = 0.65/r_n, ^[28]. Hence, r_n = 0.8 and 0.65 correspond to M = 9.0 and 7.5 for the RT and RT/2-motion, respectively.

The results in the SBM and the EBM Stage-Ⅰ (Step-1 to 6) l for the RT-motion do not look mutually contradictory in that all the saturated layers liquefy. For the RT/2-motion, however, only L2 liquefies in EBM because of the drastic reduction of upward energy from the RT-motion. In contrast, the SBM results are not markedly different between RT and RT/2 with only small increases in F_s. As an interesting trend generally observed in the comparison, liquefaction tends to occur more easily in shallower layers in EBM, while vice versa in SBM in the uniform soil model.

Hereafter, the EBM can evaluate the maximum strain γ_DAmax in Steps 7 and 8 of Stage Ⅱ as listed in bold letters in the columns of strain evaluation in Table 2. In the RT-motion, where the four layers liquefied in Stage-Ⅰ, the upward wave energy is equally divided by m = 4, from which γ_DAmax can be calculated in each layer by Eq. (72) as 17.0% ~ 8.4% and beyond 7.5% of liquefaction onset. Note that the strain tends to be higher in weaker layers where liquefaction occurs earlier in the sequence. For the RT/2-motion where L2 is a single layer to liquefy, the strain can easily be determined as 9.3% because all the wave energy is exclusively consumed there.

Based on the strain evaluated above, the 1D surface settlement is computed by following Step-9 and 10 in the settlement evaluation columns with bold letters. The upper limit strain ε_vmax is determined as 3.40% from Eq. (73) for N₁ = 8.0, F_c = G_c = 0, and the volumetric strains ε_v are calculated from γ_DAmax in individual layers by Eqs. (74) or (75). Then, they are multiplied by the layer thickness H = 2 m and summed up to determine the surface settlement S in Eq. (76). Eventually, the settlements S = 16.2 cm and S = 3.2 cm are obtained for the RT and RT/2-motion, respectively, which are widely varied despite the same acceleration. Thus, the EBM proposed here highlights a significant impact of wave energy or frequency content of earthquake motions other than acceleration intensity on the liquefaction-induced settlement without resorting to sophisticated effective stress numerical analyses.

5.3.2.   Liquefaction case of loose sand fill by small PGA earthquake

A gently inclined farmland artificially filled with loose volcanic sandy soil liquefied and fluidized, leaving surface depression behind as shown in the photographs of Figure 43 (a) in Kitami city in Hokkaido, Japan, during the 2003 Tokachi-oki earthquake (M = 8.0) ^[27]. The site was 230 km from the hypocenter of the offshore subduction earthquake, and the maximum acceleration near the site was only 0.056 g.

The strange liquefaction with underground fluidization seems to have occurred along an old shallow valley filled with sand relatively thick, 4~7 m, and the water table was GL. −1~−2 m ^[23,27]. The upslope portion of the farmland, 200 m long and 50 m wide, was depressed by 3.5 m maximum, but the ground surface remained intact with marginal fissures and no sand boils there. On the downslope side, the boiled sand erupted collectively at several ejection holes and flowed 1 km away along a ditch. Figures 44 (a) and (b) show in situ soil investigation points P1 to P12, and depth-dependent SPT N-values converted from SWS (Swedish Weight Sounding) carried out there after the earthquake. The sandy fill contained 33% non-plastic fines with the mean grain size D₅₀ = 0.02 mm.

In Table 3, representative soil models at the two points, P1 and P7, are tabulated. The liquefiable sand fill was sliced into 5–7 layers 1 m thick each with soil density ${\rho _{sat}}$ = 1.8 t/m³. The converted SPT N-values in individual layers were corrected to N₁ corresponding to effective overburden stress ${\sigma '_v}$ = 98 kPa using ${N_1}{\text{ = 1}}{\text{.7}}N/{\text{(}}\sigma _v^{'}/{p_0}{\text{ + 0}}{\text{.7)}}$, and the S-wave velocity V_s (m/s) was determined from a corresponding N-value using V_s = 80 N^1/3.

The weak acceleration motion with a long duration of more than a hundred seconds in Figure 43 (b) recorded at K-NET Kitami (EW-direction), 10 km far from the site, was given at the surface of the soil model. The velocity response spectrum in (c) indicates a very long predominant period of about 2.5 s presumably due to the large magnitude of M_J = 8.0 earthquake and long hypocenter distance R = 230 km. Figure 43 (d) exemplifies upward energies calculated at Point 7 in the one-dimensional equivalent linear analysis, which merely exceeds a few kJ/m² in the fill.

In the SBM, values of N_L = 15 were first determined from CRR₂₀ in the design code ^[57] using N₁ and F_c listed in Table 3. It was converted to the corresponding dissipated energy for the initial liquefaction $\sum {{{\Delta W} \mathord{\left/ {\vphantom {{\Delta W} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}}$ using Eq. (5) for EBM, and the normalized strain energy $\sum {{{W_ - ^*} \mathord{\left/ {\vphantom {{W_ - ^*} {{{\sigma '}_c}}}} \right. } {{{\sigma '}_c}}}}$ was calculated from it by Eq. (67) in each soil layer. The liquefaction energy capacity $W_ - ^*H$ for H = 1 m thick soil layer to liquefy was calculated using corresponding average confining stress ${\sigma '_c} =$${\sigma '_v}{{\left({1 + 2{K_0}} \right)} \mathord{\left/ {\vphantom {{\left({1 + 2{K_0}} \right)} 3}} \right. } 3}$, assuming the earth-pressure coefficient at rest K₀ is 0.5.

Table 3 compares the liquefaction potentials of EBM with those of SBM at P1 (a) and P7 (b). In the EBM, liquefaction is to occur from the groundwater level GL-1 to -2 m down to GL -3 to -4 m at P1 and P2, respectively, as shaded in the table columns. Moreover, the liquefaction potential is higher (AER is smaller) for shallower soils as seen in the uniform sand layer. In contrast, F_L-values in SBM (for two stress reduction coefficients r_n = 0.65 for M = 7.5 as a default value and 0.70 for the 8.0 earthquake, respectively) are considerably larger than 1.0 at all depths both for r_n = 0.65 and 0.70 predicting no possibility of liquefaction at all.

Figure 44 (c) depicts liquefaction potentials AER by EBM and F_s by SBM (r_n = 0.70 for the M_J = 8.0 earthquake) directly compared at multiple points near the periphery of the subsided zone. Obviously, liquefaction is very unlikely to occur at all the points according to SBM, while it is highly probable in EBM at shallower depths beneath the water table. The fundamental reason for this gap seems to rest on very small PGA (only 5% of g) versus non-negligible energy demand due to long predominant period and long duration as indicated in Figure 43 (b) to (d).

Thus, this liquefaction case study suggests better applicability of EBM to a large-magnitude far-field earthquake, in which the acceleration was too small to be accountable for the liquefaction onset by adjusting the coefficient in SBM. In contrast, the seismic wave energy was sufficient to liquefy sands in EBM because of the long duration and predominant period of the ground motion.

5.3.3.   Takasu elementary school during the 2011 Tohoku earthquake

During the 2011 Tohoku earthquake (M = 9.0), reclaimed areas along Tokyo Bay, Japan, more than 350 km from the offshore causative fault, underwent widespread liquefaction. Takasu Elementary School in Urayasu City, Chiba prefecture, was focused as one of the typical study sites ^[27], where extensive liquefaction occurred and induced large subsidence accompanying huge sand boils of more or less 10 centimeters thick as photographed in Figure 45 (a). The soil consisted of landfill (B1, B2), hydraulic fill (F), alluvial sand (As1, As2), and clay (Ac, Nac) underlain by stiff gravelly base Dg at GL-43.7 m as illustrated in Figure 45 (b) together with their layer-by-layer soil properties. Fines contained in the soils in large percentages are assumed here to be non-plastic, because the ejecta coming out from the liquefied layers contained lots of fines that were non-plastic; however, according to some soil investigation reports before the earthquake, the plasticity of in situ fines were very spatially variable (Kokusho et al. 2014 ^[60], Kokusho 2015 ^[61]).

An equivalent linear one-dimensional SH-wave propagation analysis was carried out using V_s-logging data and a downhole acceleration record in the EW direction at 14 m below the ground surface provided by the Earthquake Research Institute, the University of Tokyo ^[27]. As depicted in Figure 45 (c), the recorded peak acceleration was 0.096 g downhole, and the upward wave energy calculated based on the SH-wave propagation analysis was E_uf = 146 kJ/m² in the gravelly base layer at the end of the major motion.

Table 4 shows the EBM liquefaction evaluation compared with SBM in the top 16 m soils from the ground surface, comprising 16 layers of H = 1 m each (with the top layer unsaturated), together with pertinent parameters such as N₁-values and fines content F_c, from which (σ_d/2σ'_c)_triax. for N_L = 15 was determined using ^[57], and the dissipated energy ΣΔW/σ_c' was further derived from it by Eq. (5).

The layers with AER smaller than 1.0 in EBM Stage-Ⅰ and F_s smaller than 1.0 in SBM are judged liquefiable and shaded in the table. The SBM employing r_n = 0.80 for the M = 9.0 earthquake estimates all the saturated layers to liquefy, while the EBM Stage-Ⅰ predicts not all but quite many layers belonging to hydraulic fill (F), B2, and As1 to liquefy.

Following the EBM results obtained in Stage Ⅰ, induced strains are evaluated in the EBM Stage-Ⅱ. As the total number of layers to liquefy is m = 11 in Step-6, the upward energies allocated to individual layers are 1/11 of the total energy, from that the maximum induced shear strain γ_DAmax can be determined in Eq. (72) as listed in the right of Table 3 with bold letters. The induced strains γ_DAmax, assumed as 7.5% in all the liquefied layers in Stage-Ⅰ, are recalculated quite differently, as a maximum of 56.6% to a minimum of 4.4%. Consequently, the number of layers beyond the initial liquefaction (γ_DAmax > 7.5%) reduces to 6 (fewer than 11 in Stage-Ⅰ) as a result of strain concentration in fewer layers, which is seemingly closer to actual liquefaction behavior than Stage-Ⅰ.

Maximum volumetric strains ε_v are calculated by Eq. (74) as tabulated in Table 3 following Step-9 using ε_vmax = 3.61 ~ 4.91% determined by Eq. (73) according to N₁ and F_c (G_c can be assumed to be 0) in those layers, which were judged to liquefy in Stage-Ⅰ. In one layer with γ_DAmax = 56.6% exceeding the upper strain limit γ_DAmax = 20%, ε_v is calculated by Eq. (75). Then, the ground surface subsidence can be obtained in Eq. (76) as the sum of settlements of individual layers of thickness H = 1.0 m. Here, not only those beyond the initial liquefaction (γ_DAmax > 7.5%) but also those before are involved in calculating the subsidence if they have been judged liquefiable in Stage-Ⅰ.

The calculated subsidence of 24.8 cm at the bottom of Table 3 seems to be smaller than the soil settlement, which was around 40 cm observed by the present author during his reconnaissance at a pile-supported building of the Takasu elementary school as photographed in Figure 45 (a). However, the difference between the two values may probably be narrowed because the ejecta in large volumes (more or less 10 cm thick on average) will help increase the ground subsidence, which is defined solely in Eq. (76) from the volume contraction by undrained cyclic loading.

5.3.4.   Maihama Ni-Chome during the Tohoku earthquake

In another liquefied site during the same Tohoku earthquake, Maihama Ni-Chome, also in a reclaimed area of Urayasu City (Azuno & Kokusho 2020) ^[62], the soil consisted of landfill (B), hydraulically reclaimed clay and sand (Rc, Rs), alluvial sand (As, with high silt content in the lower part), and alluvial clay (Ac1, Ac2), underlain by stiff gravelly base layer at GL.-51 m, as illustrated in Figure 46 (a). The acceleration motion in (b) was obtained by deconvoluting a record a few km distant at a KIK-NET Urayasu station (EW) (NIED 2021 ^[63]), and given at the base layer of this site. The equivalent linear SH-wave propagation analysis was conducted to calculate dynamic shear stress τ_d, CSR and upward energies E_u as depicted in the bottom of (b) at different levels of the soil profile.

Table 5 shows the soil model developed in the research ^[62] for evaluating liquefaction potential in the shallower soils with the water table GL.-3.1 m, wherein 11 layers of H ≈ 1 m each are chosen excluding clayey soils. From N₁-values and fines contents F_c, (σ_d/2σ'_c)_triax for N_L = 15 was determined using ^[57], and the dissipated energy ΣΔW/σ_c' was further derived from it by Eq. (5).

The evaluation results by the EBM Stage-Ⅰ in Table 5 indicate that 5 shaded layers (2 layers of Rs and 3 layers of upper As) are to liquefy ^[62]. In the SBM, two layers shaded each at Rs and As are to liquefy for r_n = 0.8 corresponding to the M = 9.0 earthquake according to ^[59] to be slightly similar to the EBM in terms of liquefied depth.

Near this site, Taira et al. (2012) ^[64] carried out post-earthquake drilling/sampling of soils and X-ray CT-scan imaging of core samples for the depth GL.−2~−13 m, and observed that, in the hydraulically reclaimed soils of GL. −6.15~−8.85 m, the original sedimentary structures were completely obliterated because of liquefaction. The observation seems compatible with the present evaluation of EBM or SBM, despite the slight difference in soil profiles due to the horizontal distance of 450 m between the two points.

At this site, a 1D effective stress liquefaction analysis was also conducted to compare the results utilizing the triaxial liquefaction data on in situ intact specimens ^[62]. Despite sophisticated constitutive laws of soil properties incorporated in a commercially available widespread computer code, a pressure buildup ratio of r_u = 0.90 was recommended there as a criterion to determine the onset of liquefaction. While the analytical result using the r_u = 0.90 criterion was mostly compatible with the corresponding EBM & SBM results, r_u = 0.90 was found to correspond to the axial strain only ε_DA = 0.5% and the normalized dissipated energy only ΣΔW/σ_c' = 0.005 ^[62] in the program, which is seemingly far smaller than normally anticipated. Hence, from the viewpoint of the uniqueness of energy, it seems recommendable to pay more attention to the compatibility with the energy concept to examine the reliability of sophisticated but quite variable/tricky numerical nonlinear analyses often employed in practical design for liquefaction.

Such research efforts to calibrate nonlinear numerical results in the energy concept have begun (e.g., ^[65,66]) for sophisticated numerical analyses of various structures. In other words, the EBM is expected to serve as a reasonable/common scale to measure the reliability of sophisticated nonlinear numerical evaluations without no closed-form rigorous solutions.

Following the EBM Stage-Ⅰ results, induced strains are evaluated in Step-7 and 8. As the energies allocated to individual layers are 1/5 of the total because of m = 5, the maximum induced shear strains γ_DAmax can be individually determined by Eq. (73) as listed in one of the right columns of Table 5 as 39.7%~4.8%. Among them, only those with the sequence number 1~4 exceed the strain γ_DAmax = 7.5% for the initial liquefaction to make the number of liquefied layers smaller from 5 to 4 in the Stage-Ⅱ evaluation.

As for the surface settlement, the volumetric strains ${\varepsilon _v}$ are calculated by Eq. (74) or (75) in those layers that were judged to liquefy in Stage-Ⅰ following Step-9 using (ε_v)_{γDAmax = 20%} = 2.92~3.96% determined by Eq. (73). In the three layers with their γ_DAmax-values exceeding 20%, ε_v is calculated by Eq. (74). The calculated subsidence of 11.8 cm by Step-10 at the bottom of Table 5 is slightly smaller than 15~17 cm reported ^[34] as the ground settlements relative to tip-supporting piles nearby. The gap between the two values may be narrowed in this site again by ejecta thicker than a few centimeters on average observed in the same report ^[34].

6.   Summary and conclusions

In this review article, the energy-based liquefaction evaluation method (EBM) has been addressed in two different aspects in comparison with the conventional stress-based method (SBM): Energy capacity and Energy demand for liquefaction. How to compare the capacity and demand energies to make simplified and practical liquefaction evaluations has been also discussed. Furthermore, steps for simple EBM liquefaction predictions have been presented and demonstrated by example-case histories. Major outcomes from the reviews are as follows.

Energy Capacity for Liquefaction:

(1) The capacity of liquefaction can be uniquely determined by cumulative dissipated energy ΣΔW regardless of earthquake motions much better than applied stress histories employed in the SBM as demonstrated in various lab tests conducted on reconstituted and intact soil samples.

(2) The energy ΣΔW is almost uniquely correlated with pore-pressure buildup ratio r_u or induced strain amplitude γ_DA with no regard to cyclic stress amplitudes or the number of cycles, and waveforms. With increasing ΣΔW, the pore-pressure ratio r_u = Δu/σ'_c tends to climb up to r_u = 1.0 corresponding to initial liquefaction, and stay there for larger ΣΔW. The strain γ_DA keeps monotonically increasing almost in proportion to the energy with a slight nonlinearity even beyond the initial liquefaction. Hence, the energy can serve as a good index to evaluate the liquefaction intensity such as induced strain even beyond the initial liquefaction.

(3) The uniqueness of dissipated energy for liquefaction can be confirmed in cyclic simple shear tests which can best mimic in situ stress conditions during earthquakes. For cyclic triaxial tests in contrast, care is needed in dealing with the measured dissipated energy because the uniqueness of energy appears to fail due to the stress anisotropy inherent to the cyclic loading triaxial test system.

(4) The energy capacity in EBM is often expressed in the cumulative dissipated energy normalized by initial effective confining stress σ'_c as ΣΔW/σ'_c. This is reasonable because ΔW (energy per unit volume; the same dimension as stress) can be nondimensionalized by σ'_c (stress), though the correlation of ΣΔW/σ'_c versus pore-pressure ratio r_u or induced strain ε may not be fully independent of σ'_c in soil tests.

(5) The CSR versus number of cycles N_c curve employed in the SBM can be interpreted as an equal-energy line in the EBM corresponding to a particular pore-pressure buildup ratio or induced strain. Hence, the liquefaction performance can be uniquely determined in the EBM solely by the dissipated energy or the capacity energy without referring to an equivalent number of cycles for liquefaction N_L or the stress reduction coefficient r_n, unlike the SBM.

(6) Thus, an empirical correlation has been established to evaluate the capacity energy for liquefaction from the CRR of the SBM for particular N_L values. Using triaxial tests on intact samples recovered from various sites, Eqs. (5) & (6) have been obtained connecting the CRR for N_L = 15 or 20 to the corresponding energy ΣΔW/σ'_c uniquely with no regard to the difference in soil properties.

(7) Concerning the influence of earthquake wave irregularities investigated by cyclic torsional liquefaction tests using six recorded motions with widely spanned damage level D_L of the fatigue theory, the energy ΣΔW/σ'_c is effective in uniquely determining the induced strain at least up to the initial liquefaction of γ_DA = 7.5%. This indicates that the energy concept is superior to the fatigue theory, a theoretical basis of the SBM in dealing with wave irregularity, as the liquefaction susceptibility can be governed solely by the dissipated energy despite the widely varied D_L-values.

(8) For strain ranges beyond initial liquefaction (γ_DA > 7.5%), however, the energy for e.g. γ_DA = 15% or larger tends to decrease slightly with an increasing number of cycles or decreasing stress amplitudes in irregular motions. This is probably due to stress bias temporarily working at irregular stress peaks where larger biased stress tends to make the soil more resistant and require larger dissipated energy in dilative soils considered here.

Energy Demand for Liquefaction to Compare with Capacity:

(9) Compared to the energy capacity, the energy demand for liquefaction has scarcely been discussed and more often dealt with implicitly through numerical analyses directly using design earthquake motions. Besides, the seismologically given empirical formula by Gutenberg ^[6] was utilized in early-time liquefaction evaluations by several researchers, though the energy demand was not directly compared with the capacity but used only as an index for judging liquefaction susceptibility based on previous liquefaction case histories.

(10) One of the highlights of the energy-based liquefaction evaluation developed here is the capability to quantify the total amount of demand energy from upward design earthquake waves and to directly compare them with the energy capacity for liquefaction, not only to predict liquefaction potential, but also to evaluate the extent of liquefaction at a glance. This includes induced strains in individual layers and ground settlements without resorting to complicated nonlinear effective stress analyses.

(11) In the EBM developed here, the demand energy of a given earthquake E_uf is determined in Eq. (68) as cumulative upward wave energy using a site-specific design motion. Numerous vertical array strong earthquake motion data during strong earthquakes of recent decades in Japan demonstrated that the E_uf-value tends to monotonically decrease with decreasing soil depth depending on the corresponding impedance ratio.

(12) In those sites where specific design motions are not available, empirical formulas, Eq. (69), derived from the vertical array database may be used to roughly evaluate Eu at particular ground depths by employing the well-known formula seismologically given in Eq. (70).

(13) The energy dissipation rate during wave propagations in situ ΔE/E_u is very similar to that during cyclic loading tests ΔW/2W_－ and is mostly reproducible in laboratory cyclic loading tests if the energy-recycling effect during cyclic loading is taken into account when evaluating 2W_－. Since the energy dissipated for liquefaction, in situ ΔE and laboratory ΔW, respectively, should be the same; ΔE = ΔW × H, for a soil layer of thickness H, the cumulative wave energy E_u coming up through unit horizontal area is to be compared with cumulative strain energy Σ2W_－×H in the liquefaction evaluation.

(14) Because of 100% reflection of the SH-wave at a free ground surface of zero shear stress, the strain energy cannot be present near the surface to compensate for the dissipated energy for liquefaction. In this regard, wave propagation studies using harmonic and earthquake waves have led to a simple approximation that only half of the upward wave energy is available to compensate for the dissipated energy for liquefaction in soils shallower than a quarter-wavelength λ/4 from the surface. Considering that λ/4 is normally larger than 20 m in practical liquefaction problems, the upward energy Eu should be compared with the strain energy twice larger; ΣW_－^* ≡ 2×Σ2W_－.

Evaluation Steps and Case Study Examples:

(15) Ten steps of the energy-based liquefaction evaluation divided into Stage-Ⅰ and Stage-Ⅱ are proposed. In Stage-Ⅰ (Step 1~6), the capacity and demand energies are quantified from site/soil and earthquake conditions and compared layer by layer individually to know the spatial extent of initial liquefaction corresponding to γ_DAmax = 7.5% (pore-pressure buildup ratio r_u = 100%) and their sequence of occurrence in terms of AER (Accumulated Energy Ratio) in Eq. (71) where liquefaction occurs only for AER$\leqslant$1.0. Thus, unlike the conventional SBM, the EBM has the potential to decide the extent of liquefaction within the availability of the total energy demand.

(16) In Stage-Ⅱ (Step 7~10), the same upward energy E_u is equally allocated exclusively to those layers judged as liquefiable in Stage-Ⅰ to evaluate induced γ_DAmax and post-liquefaction ground settlements S. The equal energy allocation may not be theoretically rigorous but tends to approach nearer to actual liquefaction performance than Stage-Ⅰ by concentrating larger induced strains in a smaller number of liquefied layers.

(17) The impact of frequency content in input motions on liquefaction potential has been demonstrated by a simple uniform model where, despite the same input motion of the same acceleration amplitude given to a uniform sand model, the liquefaction potential, induced strain, and ground settlement are evaluated significantly lower when the time axis is compressed by 1/2 in the EBM, while the corresponding difference was found marginal in the SBM. Also found was a general trend that liquefaction tends to occur earlier in shallower layers on uniform sand layers in the EBM, while it is reversed in the SBM.

(18) A liquefaction case history in a loose sandy fill with its recorded acceleration of only 0.056 g during a far-field earthquake has been successfully replicated using a recorded motion nearby in the EBM, presumably because the wave energy was large enough due to the large magnitude of M_J = 8.0, while the SBM predicted no possibility of liquefaction, quite different from the actual behavior.

(19) Two severe liquefaction cases in hydraulically filled residential lands during a far-field M_J = 9.0 earthquake have been studied using acceleration records nearby. Both the EBM (Stage-Ⅰ) and SBM have similarly predicted heavy liquefaction involving most layers in shallow depth, though, in the EBM Stage-Ⅱ, the number of liquefiable layers tends to decrease from Stage-Ⅰ, while their induced strains tend to concentrate to a smaller number of layers with higher liquefaction susceptibility. The post-liquefaction settlements calculated in the EBM Stage-Ⅱ have been compared with the performance observed at the two sites and found a fair agreement if the thickness of ejected sands is considered.

Consequently, the energy-based liquefaction evaluation has been comprehensively reviewed from various aspects to recognize that it has been mature enough for practical use in engineering design. It can consider a variety of earthquake motions exclusively in terms of demand energy. The capacity energy, uniquely connected with the resistant stress CRR in the SBM, is compared with the demand energy to readily evaluate liquefaction potential, induced strain, and settlement without resorting to sophisticated numerical analyses.

Moreover, considering the uniqueness of energy, it can serve as a common scale to measure the reliability of nonlinear numerical liquefaction evaluation tools which attract designers because of the modeling capability of complicated boundary value problems but have no rigorous solutions to depend upon.

There are limitations in the EBM that need to be overcome in future studies. The most crucial is the accumulation of in situ liquefaction case histories in terms of surface manifestations including lateral deformations and settlements to demonstrate the applicability of the EBM in comparison with the SBM for various seismic and geotechnical conditions. Induced strains and settlements calculated in the EBM may be compared first with in situ observation or model shaking table tests, and also with sophisticated stress-based numerical tools to know how well this simplified evaluation works in engineering designs. Another challenge is how the effect of initial shear stress should be considered so that the EBM can be reasonably applied to liquefaction designs concerning slope failures and soil structure interaction problems. With these efforts, the EBM is hoped to serve as a major player in liquefaction-related engineering designs and also a common scale to calibrate numerical results of practical problems in terms of energy.

Acknowledgments

Professor Motoki Kazama of Tohoku University, Japan, and all members of the research committee on Energy-Based Liquefaction Evaluation are gratefully appreciated for their valuable contributions to the committee activated from 2015 to 2018 in JGS (the Japanese Geotechnical Society). Various outcomes in the committee have been utilized in this review paper to demonstrate the usefulness of the EBM compared with the SBM. Among them, Shunsuke Tanimoto, a Senior researcher at the Public Works Research Institute in Japan, dedicated to publicizing a set of cyclic triaxial test results conducted on intact soils sampled in numerous sites liquefied during the 2011 Tohoku earthquake (M_W = 9.0) to be available for data analyses in terms of dissipated energy. Keigo Azuno, Senior Engineer of Chuo Kaihatsu Corporation, Japan, contributed a lot to applying the EBM to case histories and comparing it with the SBM and the effective stress numerical analyses. NIED (National Research Institute for Earth Science and Disaster Prevention) Japan, who disseminated numerous K-NET and KiK-net data is gratefully acknowledged. Last but not least, the great devotion of graduate & undergraduate students of the Civil & Environmental Engineering Department at Chuo University who conducted cyclic triaxial tests as well as cyclic torsional shear tests over many years, without which this paper could never have been realized, is also very much appreciated.

Use of AI tools declaration

The author declares that Artificial Intelligence (AI) tools have not been used in the creation of this article.

Conflicts of interest

The author declares no conflict of interest concerning this article.

Data availability statement

All data presented in this study are available for research purposes upon request directly to the present author.