The use of the SPT-based seismic soil liquefaction triggering evaluation methodology in engineering hazard assessments

Graphical abstract


Method details
A summary of proposed SPT-based probabilistic and deterministic liquefaction triggering methodology: background A new set of probabilistic and deterministic seismic soil liquefaction triggering relationships is presented in Cetin et al. [1], on the back analyses of standard penetration test liquefaction triggering case histories, which are fully documented in Cetin et al. [6]. The use in forward (design) assessments of these new relationships requires i) the correct understanding of the protocols behind case history processing, and ii) the consistent use of these protocols in design assessments. This manuscript is intended to discuss these protocols, and to guide engineers through the correct and consistent use of them, along with a discussion on "tricks" alongside these protocol. For this purpose, an illustrative soil profile shaken by a scenario earthquake is used to outline the use of the proposed methodology. These new relationships are given in Eqs. (1) and (2).
In Eq. (1), P L is the probability of liquefaction in decimals (i.e. P L = 30% is input as 0.30), CSR s 0 v ;a¼0;Mw is not "adjusted" for vertical effective stress or magnitude/duration effects (corrections are executed within Eq. (1) itself), FC is percent fines content (by dry weight) expressed as an integer (e.g.: 12% fines is input as FC = 12) with the limit of 5 FC 35, P a is atmospheric pressure (1 atm. = 101.3 kPa = 2116.2 psf) in the same units as the in-situ vertical effective stress (s 0 v ), and F is the standard cumulative normal distribution. The cyclic resistance ratio (CRR), for a given probability of liquefaction can be expressed as given in Eq. (2), where F À1 P L ð Þ is the inverse of the standard cumulative normal distribution (i.e. mean = 0, and standard deviation = 1). For spreadsheet construction purposes, the command in Microsoft Excel for this specific function is "NORMINV (P L ,0,1)". In Fig. 1, factor of safety (FS) values corresponding to probabilities of liquefaction 5, 20, 50, 80, 95% are also presented.
The new probabilistic boundary curves, as shown in Fig. 1, are estimated by considering the uncertainty due to model error only. In this figure, the dots and circles represent liquefied and nonliquefied case histories compiled from available literature. As stated earlier, a complete documentation of these case histories including their source references is presented in Cetin et al. [6]; hence will not be repeated herein. These new triggering relationships, which are discussed in detail in Cetin et al. [1], will be referred to as CEA2018 (Cetin et al.), hereafter. The values of the model coefficients (i.e. u i ) as defined in Eqs. (1) and (2), are listed in Table 1. Approximate factors of safety, FS values can be estimated by Eq. (3) based on the assumption that P L of 50% corresponds to a best-estimate factor of safety value of 1.0.
It should be noted that Eqs. (1) and (2) In Eq. (4), s input represents the consolidated uncertainty of the liquefaction engineering input parameters, and as discussed in CEA2018, a scaling factor of u 7 needs to be systematically applied. This factor (u 7 ) is one of the regressed parameters of the overall triggering relationship. If the higher-order terms are eliminated, then s input can be estimated as given Eq. (5).
The estimation of the uncertainties in input parameters will be discussed later in this manuscript.

Recommended use of the proposed liquefaction triggering assessment methodology
The proposed new probabilistic correlations can be used in two ways. They can be used directly, all at once, as summarized in Eqs. (1) and (2), oralternatively, they can be used "in parts" as most of the previous similar conventional methods. A flow chart is presented in Fig. 2, which summarizes the required assessment protocols for the consistent use of recommended methodology.
For illustrating the use of the proposed methodology and details of assessment steps described in Fig. 2, a generic site is to be assessed. As presented in Fig. 3, this illustrative soil site is composed of three layers: There exists a 3 m thick high-plasticity clay (CH) layer at the surface, which is underlain by a 5 m thick, potentially liquefiable, medium-dense silty sand (SM) layer. At and below 8 m depth, there exists another highly plastic clay (CH) layer. The earthquake scenario is deterministically selected as a M w = 6.8 event, which is assumed to produce a maximum peak ground acceleration (a max ) of 0.28 g. The coefficient of variation of a max is assumed as 0.15. The ground water table depth (h w ) is assumed to be 3 m with standard deviation of 1 m. The other necessary input parameters are also presented along with the soil profile given in Fig. 3. Site and soil specific moist and saturated unit weights are estimated based on laboratory test results as given in the same figure.
Step 1: Determination of Mean Input Parameters and Their Uncertainties The input parameters required for seismic soil liquefaction triggering assessments will be discussed under two separate sections: the ones related to a) soil profile and soil characteristics, and b) seismic scenario.

Soil profile
Estimation of mean and standard deviation of N 1,60 The critical layer has multiple and consistently scattered SPT blow counts (N 60 ). The field N values are corrected for effective normal stress (C N ), hammer energy (C E ), rod length (C R ), sampler (C S ), borehole diameter (C B ), and procedural effects to fully standardized N 1;60 values, as given in Eq. (6).
The corrections for C N , C R , C S , C B and C E correspond closely to those recommended by NCEER Working Group (NCEER [7], also given in Youd et al. [8]), and they are summarized in Table 2 for the sake of completeness.
On the basis of Taylor's expansion, the first-order approximations of the mean and variance of N 1,60 are given in Eqs. (7) and (8) with the assumption that correction factors are all exact (i.e.: uncertainties in correction factors are ignored).
For design purposes it is recommended that the field N-values within a critical layer from one or more boreholes at a soil site to be corrected as recommended in Eq. (6) to estimate N 1;60 values. These are then plotted vs. the topography-corrected depth below the ground surface. In many cases, a given soil layer will be found to contain an identifiable critical sub-stratum based on a group of localized low N 1;60 -values. Occasional high values, assumed to be gravel and not apparently representative of the general characteristics of the critical stratum, are removed if they sit outside the main cluster of points. Similarly, though less often, very low N 1;60 values that are much lower than the apparent main cluster of points representing the stratum are typically associated with locally high fines content and eliminated if the soils have significant fines and if it appears that the fines are plastic in nature. The remaining, corrected N 1;60 values are then used to evaluate both the mean of N 1;60 within the critical stratum, and the variance in N 1;60 . When applicable, SPT blowcounts above the groundwater table may be also used to characterize the relative density state of the critical layer.  Table 3 presents the corrected SPT blowcounts along with the shear wave estimations, details of which will be given in Step1.A.4. After applying the corrections for the SPT blowcounts, the critical silty sand layer is characterized by the mean and the standard deviation of scattered SPT blow counts, as 11.4 blows/ft and 2.3 blows/ft, respectively. For the estimation of rod length, a 1.2 m of stick up is assumed.

Estimation of mean and standard deviation of FC
The fines (silt and clay particles) of soils are widely expressed by fines content, which is defined as the portion of soil particles by mass finer than the No. 200 sieve (0.074 mm). Following the similar procedure outlined for the estimation of mean and standard deviation of N 1;60 values, the mean and standard deviation of FC values are estimated as 14 AE 5.3, respectively. Note that the FC values are limited to be within the range of 5-35 %. If any outlier data exists, it needs to be excluded or critical layer needs to be divided into more homogenous soil sub-layers. In the literature there exist different opinions regarding the effects of fines on penetration (Cubrinovski and Ishihara [10], Shahien [11]) and liquefaction resistances (Cetin et al. [1], Youd et al. [8], Tokimatsu and Yoshimi [12], Seed et al. [13], Table 2 Recommended corrections for SPT equipment, energy and procedures.
where the effective stress, s 0 v , and reference stress,Pa, are in the same units.
where d = rod length (or "stick-up") from the top of the SPT sampler to the striking point at the top of the rod. C S For samplers with an indented space for interior liners, but with liners omitted during sampling, where ER (energy efficiency ratio) is the fraction or percentage of the theoretical SPT impact hammer energy actually transmitted to the sampler, expressed as % The best approach is to directly measure the impact energy transmitted with each blow with instrumented rod. When available, direct energy measurements were employed. The next best approach is to use a hammer and automatic (mechanical) trip hammer release system that has been demonstrated to deliver repeatable energy, and which has been calibrated based on direct (-instrumented) energy measurements.
Otherwise, ER must be estimated. For good field procedures, equipment and monitoring, the following approximate guidelines for SPT performed with rope and cathead are suggested: For lesser quality fieldwork (e.g.: irregular hammer drop distance, excessive sliding friction of hammer on rods, wet or worn rope on cathead, etc.) further judgmental adjustments are needed.
Notes: (1) Based on rope and cathead system, two turns of rope around cathead, "normal" release (not the Japanese "throw" (Seed et al. [9]), and rope not wet or excessively worn.
(2) Rope and cathead with special Japanese "throw" release. (See also Note 4. (3) For the ranges shown, values roughly central to the mid-third of the range are more common than outlying values, but ER and C E can be even more highly variable than the ranges shown if equipment and/or monitoring and procedures are not good. (4) Common Japanese SPT practice (Seed et al. [9]) requires additional corrections for Borehole diameter and for frequency of SPT hammer blows. For "typical" Japanese practice with rope and cathead, donut hammer, and the Japanese "throw" release, the overall product of CB x CE is typically in the range of.1.0-1.3.

Estimation of mean and standard deviation of total and effective overburden stresses
Total and effective vertical stress estimations require the estimation of soil unit weights along with the depth of ground water table. If soil specific unit weight data is missing, the recommended values, which were also used in CEA2018 for back analyses of field performance case histories, as presented in Table 4, can be used. For this illustrative site, the mean values of moist (g moist Þ and saturated (g sat Þ unit weights are estimated based on laboratory test results as given in Fig. 3.; whereas, their standard deviations are assumed as 0.5 kN/m 3 . In addition to the uncertainty of the mean estimates of soil unit weights, the inexact estimation of the depth to water table affects the accuracy of vertical effective stress estimations. In the literature, clear definitions for uncertainty estimations of phreatic surface depths are not available. Based on expert opinions the following simple procedure is proposed: The mean values for total and effective vertical overburden stresses are estimated for the middepth of the critical layer as given in Eqs. (9) and (10), respectively. The uncertainty of these input parameters are estimated as given in Eqs. (11)-(13) on the basis of first order approximation as Cov In these equations, g 1 and g 2 represent the moist and saturated unit weights, respectively. On the other hand, m h and s h are the mean and standard deviation of the critical depth, respectively. For the illustrative site discussed herein, m h can be estimated as 5.
Estimation of mean shear wave velocity, V s;12m for r d calculations Cetin and Seed [22] r d relationship requires the estimation of a representative shear wave velocity for the upper 12 m of the soil site. If in-situ Vs measurements are available they can be directly used, otherwise Eqs. (21) and (22) from the Design Specification for Highway Bridges, Japan Road Association [23] can be used as one of many N vs. V s correlations. Note that the V s estimations presented in Table 3 were also determined via these equations.
Assuming that the average shear wave velocity values are 250, 170 and 320 m/s for the highlyplastic clay, silty sand, and claystone layers, the mean shear wave velocity for the upper 12 m is estimated as 220 m/s as shown in Eq. (24). Earthquake scenario Estimation of mean and standard deviation of ln(PGA) Estimating peak ground acceleration (a max ), at soil sites requires the understanding of both the seismicity (magnitude, source mechanism, travel path, directivity effects, etc.) and the response (both geological and geotechnical) characteristics of the site. In the literature, for the assessment of liquefaction triggering case histories, peak ground acceleration has been evaluated in the order of decreasing accuracy by using: 1 Strong ground motion recordings obtained directly at the site of interest (e.g.: Wildlife site, U.S.A., and Port Island site, Japan), 2 Site response analysis tools with a good, "representative" input motion developed from an eventspecific nearby ground motion record, 3 Site and earthquake specific attenuation relationships, derived from available strong ground motion data recorded on similar nearby soil sites, or on rock sites where the amplification or deamplification of soil sites are incorporated separately, and with reasonable azimuthal accounting for directivity and travel path effects, etc., 4 Generalized attenuation relationships (e.g.: Abrahamson and Silva [24], Idriss [25], etc.), with modifications to account for the effects of local site conditions, 5 Generalized attenuation relationships without local records to provide event specific calibration, 6 Intensity scales. (e.g.: modified Mercalli scale).
Within the scope of Cetin et al. [1] studies, case histories where a max cannot be estimated by one of the first three methods were eliminated from further consideration. This also defines an upper boundary on the uncertainty of the mean estimates of a max as to be less than that predicted by generalized attenuation relationships. In all cases, a max estimations adopted for these studies are based on the geometric mean of the two orthogonal components of available recordings. Typical attenuation relationships (e.g.: Abrahamson and Silva [24], Idriss [25]) can estimate peak ground acceleration at soil sites for a wide range of earthquake magnitudes and distances with an error term, which is dependent on a number of additional factors (e.g.: event magnitude, event type and mechanism, rupture distance, site stiffness, etc.) A typical coefficient of variation term varies in the range of 30% to 40%. Any relevant information other than magnitude and distance should improve the accuracy of the estimations, which in turn should decrease the coefficient of variation (c.o.v.) to a value less than $ 0.35.
Similarly, comparisons of the actual recorded a max values with site response analysis predictions based on "good" site characterization and seismic data revealed that the discrepancy in the matches is more typically in the range of c.o.v.10-20 %. The error represented by c.o.v. of a max reduces to < 10% for the case history sites where actual strong ground motion recordings are available at the site of interest.
These outlined recommendations define guidelines for the estimation of mean and standard deviation for ln(a max ). For illustration purposes, ln(a max )AE s amax values are chosen as ln(0.28) AE 0.15.

Estimation of mean of moment magnitude
Due to use of different earthquake magnitude scales, a conversion factor may be needed to express the magnitude in moment magnitude scale. However, this conversion may not be exact. On the other hand, the reported magnitude might be in terms of moment magnitude but due to uncertainties in the estimations of the fault rupture dimensions, the rigidity of the ruptured material or for some other reason, the documented moment magnitude itself may not be exact. Due to the relative minor importance of these input parameters in the overall model, moment magnitude of the earthquake is incorporated as a deterministic value with mean of 6.8. Step

2: Estimation of Mean and Standard Deviations for Capacity and Demand Terms
Demand term: cyclic stress ratio, CSR

Seismic site response analyses
In-situ equivalent uniform CSR s 0 v ;a;Mw can be evaluated either based on direct seismic site response analyses, or direct seismic site response and soil-structure-interaction analyses, as given by Eq. (25a,25b).
Using the "Simplified" approach Estimation of mean and standard deviation of nonlinear shear mass participation (Stress reduction) factor, r d . The stress reduction or nonlinear shear mass participation factor, r d , is estimated as given by Cetin and Seed [22] in the form of Eqs. (26) and (27).
The standard deviation of the model error term (s er d ) is defined as given in Eq. (28): For d< 12 m ($40 ft): For d ! 12 m ($40 ft): In Eqs. (26)-(28a,28b), "d" is in meters and corresponds to the depth of interest (critical depth of 5.5 m for this case), a max is in gravitational acceleration (in g's), V Ã s;12m is the time-averaged shear wave velocity over the top 12 m in m/sec calculated in the same manner as V s;30 , and "e" is the exponential symbol. A full explanation of the development of the probabilistic r d relationship is presented in Cetin and Seed [22]. By inputting the corresponding values of the illustrative case, r d and error term's standard deviation are calculated as 0.97 and 0.084, respectively as presented in Eqs. (29) and (30). Estimation of mean and standard deviation of lnðCSR s 0 v ;a;Mw Þ by simplified procedure. In-situ equivalent uniform CSR s 0 v ;a;Mw can be evaluated based on the "simplified" approach by employing Eq. (31) along with the r d relationships given by Cetin and Seed [22] as given in Eq. (29).
The mean CSR s 0 v ;a;Mw for the critical layer of illustrative case can be estimated as 0.23 as presented in Eq. (32).
Similarly, the uncertainty in CSR, where only the total stress and effective stress terms are assumed to be correlated, can be estimated as given in Eq. (33).
Coefficient of variation (d) of each term can be calculated by the ratio of standard deviation to mean (s=m). Thus, Eq. (33) can be re-written as given in Eq. (34); whereas, the correlation coefficient between total and effective stress terms is defined as given in Eq. (35).
Capacity term: cyclic resistance ratio, CRR s 0 v ;a;Mw or CRR s 0 v ¼100 kPa;a¼0;Mw¼7:5 Determination of CRR s 0 v ;a;Mw by closed form Eq. (2) CRR term corresponding to any vertical effective stress ratio and moment magnitude can be estimated by using the equation given in Eq. (2). As recommended by CEA2018, this equation will be solved for probability of liquefaction of 50%,for which F À1 P L ¼ 50%  Determination of P L from N 1;60;CS vs. CRR curves given in Fig. 1 Use of chart solution requires series of corrections which will be introduced in following sections. Estimation of K s . For the estimation of stress-scaling factor, K s , either the closed form solution given in Eq. (44) or the chart solution presented in Fig. 6 can be used. This relationship was developed based on the regression of the liquefaction performance field case history database. The histogram of the vertical effective stresses of these case histories is also presented in Fig. 6. Hence the use of the proposed K s relationship should be limited over the effective vertical stress range of 0.25 atm.
s' v 1.8 atm. Extrapolation to higher vertical effective stresses beyond 1.8 atm. for forward engineering analyses is controversial, and the readers are referred to the discussion presented in Cetin et al. [26].
Estimation of K a . As the illustrative soil site is free field level site, there exists no shear stresses acting on the horizontal plane (i.e.: a = 0). Hence no correction for K a effects is needed (i.e.: K a = 1.0).
However, it should be noted that at sloping soil sites, and at soil sites where a super-structure is overlying potentially liquefiable soils, K a correction needs to be applied as recommended in Harder and Boulanger [27], Boulanger [28], Cetin and Bilge [29].
Estimation of CSR s 0 v ¼100 kPa;a¼0;Mw¼7:5 . CSR s 0 v ;a;Mw is then adjusted with K s , K a and K Mw correction factors to convert the field and event specific CSR value to the reference CSR value, valid for s 0 v = 1 atm, a = 0 and M w = 7.5, as given in Eq. (48). Note that for level sites, K a = 1.0.
lim: CSR s 0 v ¼1atm;a¼0;Mw¼7:5 0.6 CSR s 0 v ¼1atm;a¼0;Mw¼7:5 ¼ 0:23 Á 1 1:08  Estimation of probability of liquefaction, P L . The resulting, fully adjusted and normalized values of N 1;60;CS and CSR s 0 v ¼1atm;a¼0;Mw¼7:5 can then be used as shown in Fig. 8 to estimate the probability of liquefaction triggering as 96%. It should be noted that these curves were developed by considering the uncertainty due to model error only. If the uncertainty of input parameters needs to be incorporated into the assessment, then Fig. 1 cannot be used. For N 1;60;CS = 13 blows / ft and CSR s 0 v ¼1atm;a¼0;Mw¼7:5 ¼ 0:17 , P L is estimated from Fig. 8 as 96%.
As discussed earlier, Eqs. (1) and (2) are only applicable for cases where the input parameters are exact. However, usually there exists significant uncertainty associated with these parameters in typical design applications. The proposed framework of CEA2018 allows assessment of the parameter uncertainty in forward engineering problems. For the cases where input parameters are uncertain, then in Eq. (1), instead of model error (s e ), the overall standard deviation term (s tot ) needs to be used along with Eqs. (4) and (5). By substituting the corresponding mean input parameters, the overall (consolidated) uncertainty of the illustrative case (s input ) can be calculated as presented in Eq. (50).
When the uncertainties of the input parameters are considered, then the probability of liquefaction triggering is estimated as 90% as presented in Eq. (51c).
Determination of factor of safety against liquefaction triggering, FS For the "deterministic" evaluation of liquefaction triggering, Eq. (52), which defines the ratio of capacity to demand terms, can be used to estimate the factor safety against liquefaction triggering. While using Eq. (52), it is vital to adjust the CSR and CRR terms to same stress and magnitude (duration) states. Thus, FS can be calculated by either using the values corresponding to the site and event specific state or adjusted to the reference state (i.e. s 0 v ¼ 1atm; a ¼ 0; M w ¼ 7:5) as given by Eqs.