Excess Pore Water Pressure Ratio Comparison from Empirical and Numerical Methods to Determine Liquefaction Potential in Palu, Central Sulawesi, Indonesia

2.1. The Research Area

The research area is located downstream of the river in the eastern part of Palu City. This study is based on four soil investigation borehole points, BH 1, BH 2, BH 3, and BH 4, as shown in Fig. (3). The Palu 4 Bridge is being constructed in the research area, specifically on the eastern side abutment and embankment.

2.2. Geological and Geotechnical Condition

The research location consists of alluvial and coastal deposits. Its position downstream of the Palu River results in sediment transport quite far from the weathering source, leading to sand grains with rounded shapes. The rounded grain shape reduces the interlocking between particles. Additionally, these deposits formed during the Quaternary period, meaning the sediments' cementation process has not yet occurred [8].

Palu City is located along the Palu Koro Fault, a strike-slip fault that moves at a rate of 42 mm per year [9]. This movement is associated with the high seismic activity in the area. The research location is at the mouth of the Palu River, which is directly connected to Palu Bay, resulting in a shallow groundwater table. With sandy soil, high seismic activity, and a shallow groundwater level, there are early indications that the research location is prone to liquefaction.

Based on the sieve analysis results, the lithology of the four borehole points indicates that the lithology consists of clayey silt and sandy gravel up to a depth of 20 meters. The NSPT values at the four borehole points vary, ranging from 6 to 50, as shown in Fig. (4). The four borehole points have different groundwater levels: BH 1 has a groundwater depth of 1.1 m with a drilling elevation of 2.72 m above sea level, BH 2 has a depth of 1.05 m with a drilling elevation of 2.52 m above sea level, BH 3 has a depth of 0.6 m with a drilling elevation of 2.54 m above sea level, and BH 4 has a depth of 1.74 m with a drilling elevation of 3.08 m above sea level.

2.3. Data Acquisition

The data consists of soil investigation data from Standard Penetration Test (SPT) values from four boreholes: BH 1, BH 2, BH 3, and BH 4. This study also utilizes laboratory test data.

The ground motion data used are earthquake recordings from the Pacific Earthquake Engineering Research Center (PEER) via the website https://ngawest2. berkeley.edu/, with input parameters derived from deaggregation using the Indonesian Earthquake Hazard Deaggregation Map for Seismic-Resistant Infrastructure Planning and Evaluation [10]. Once determined, the ground motion matches the target design spectra created based on the Indonesian National Standard 2833:2016 [11].

2.4. Data Analysis

The field measurement data from the Standard Penetration Test (SPT) is corrected using several correction factors, such as overburden pressure correction, correction for the energy ratio produced by the hammer blow, borehole diameter correction, SPT rod length correction, fine grain content correction, and correction for the liner of the SPT sampler [12].

After applying the corrections to the SPT values, the next step is calculating the Cyclic Resistance Ratio (CRR) with a correction factor for effective overburden pressure correction (K_σ). The earthquake magnitude used is 7.5 Mw, referring to the Palu earthquake that occurred on September 28, 2018. The CRR value is calculated using Eq. (1) mentioned below [13].

(1)

where CRR_M,σ'v represents the CRR value that has been corrected with the magnitude scaling factor (MSF) and the effective overburden pressure correction (K_σ).

The factor of safety against liquefaction requires the analysis of the Cyclic Stress Ratio (CSR) using Eq. (2). The CSR is calculated using the Peak Ground Acceleration (PGA) value of 0.794g from the Indonesian Seismic Code, with a 7% probability of exceedance over 75 years. The total vertical stress influences the CSR value (σ_v), effective vertical stress (σ'_v), maximum ground acceleration due to the earthquake, gravitational acceleration (g), and the stress reduction factor (r_d). The equation is as follows:

(2)

The Safety Factor (SF) value is the ratio between the CRR and CSR values, as shown in Eq. (3) below.:

(3)

The Boulanger and Idriss (2014) method is commonly used for identifying liquefaction potential. This method can be applied using data that are typically available, such as N-SPT values, as well as laboratory test data commonly performed, such as specific gravity, plasticity index for silt and clay, and sieve analysis. Therefore, the Simplified Procedures method is the most suitable empirical method for this research.

The obtained safety factor results are then processed using the equation given by Yegian and Vittelli in 1981 [14] to determine the excess pore water pressure ratio (r_u), as shown in Eq. (4) below, where α and β values are 0.7 and 0.19, respectively.

(4)

The Yegian and Vittelli (1981) method presents the results of the simplified procedure calculations from Boulanger and Idriss (2014), which are expressed as the excess pore water pressure ratio (r_u) to facilitate comparison with the nonlinear and finite element methods. These comparison results will be part of the analysis in this study. Although this method was introduced many years ago, it is still widely used by several recent authors [15-17].

The target spectra are determined using the PGA value, the spectral value at the ground surface for a 1-second period (S_DS), and the spectral value at the ground surface for a short period of 0.2 seconds (S_D1). The value of S_DS is obtained from https://lini.binamarga.pu.go.id/ based on reference [9]. Several parameters are needed to establish the target spectra to define the curve boundaries. The first parameter required is the initial transition period from the phase of maximum acceleration (T ₀), followed by the transition period where spectral acceleration starts to decrease as the period increases (T_S). After the period T_S, the spectral acceleration (S_a) is determined. These parameters can be calculated using Eqs. (5-7), given as follows:

(5)

(6)

(7)

After determining the target spectrum, matching is performed so that the ground motion can be adjusted to field conditions. The matching results are shown in Fig. (5). The results of the ground motion matching that has been performed can be seen in Fig. (6).

2.5. Nonlinear Site Response Analysis

Numerical analysis is conducted as a comparison to the empirical method [14, 15]. Numerical analysis can be performed using 1-D Site Specific Response Analysis (SSRA) with the DEEPSOIL v7 application to estimate the excess pore water pressure. This analysis uses the Generalized Quadratic/Hyperbolic (GQ/H+PWP) approach.

The main parameters for numerical analysis using DEEPSOIL v7 with the Dorby and Matasovic Model Parameters for sand layers are normalized excess pore pressure (U_N), equivalent number of cycles (>N_eq), current reversal shear strain (>γ_c), threshold shear strain value (>γ_tvp), curve fitting parameters (P, s, and F), dimensionality factor (f), and degradation parameter (v).

In its calculations, this model is influenced by the Fines Content (FC) and shear wave Velocity (V_S) values [18], as shown in Eqs. (8 and 9) below:

(8)

(9)

The parameter v is calculated using Eq. (10) given below, where D_r represents the relative density in %.

(10)

For clay layers, the Matasovic and Vucetic (1995) model parameters are used, with key parameters including normalized excess pore pressure (U_N), equivalent number of cycles (N_eq), current reversal shear strain (γ_C), threshold shear strain value (γ_tvp), curve fitting parameters (r and s), and curve fitting coefficients (A, B, C, and D). The values for each parameter are based on Matasovic and Vucetic (1995) for marine clay material (OCR 1.0).

The r_u analysis using the nonlinear method is performed to account for the nonlinear behavior of the soil during an earthquake. The PGA value, which is one of the triggers of liquefaction, will be adjusted according to the soil response. This method is used to determine the amplification of earthquake vibrations, which in turn increases pore water pressure, leading to liquefaction.

2.6. Finite Element Method

Excess pore water pressure can be calculated numerically using the Plaxis 2D application with the PM4Sand model. Plaxis 2D utilizes the finite element method to analyze and estimate the increase in pore water pressure. PM4Sand simulates the behavior of sand under dynamic conditions, which generates the values of pore water pressure increase, liquefaction conditions, and displacement values following the dynamic event [19].

The main parameters in the PM4Sand model are atmospheric pressure (pA), contraction rate parameter (h_po), shear modulus (G ₀), and relative density (D_R) [20]. The atmospheric pressure value used is 101.3 kPa.

The relative density (D_R) parameter is obtained using the equation from Idriss and Boulanger (2008) [21], with parameters including the corrected SPT value at 1 atm overburden ((N₁)(60) and other values as shown in Eq. (11) below:

(11)

Shear modulus coefficient (G₀) controls the small strain shear (elastic), the expression covering a range of typical densities can be used to calculate the parameter directly G₀ as proposed by Boulanger and Ziotopoulou (2015) [22] in Eq. (12) given below:

(12)

The steps to obtain the contraction rate parameter require using the Cyclic Direct Simple Shear (CDSS) test. This is conductedto calibrate the CRR value under cyclic conditions. Liquefaction is considered to occur when the peak shear strain reaches 3% after undergoing cycles in the CDSS test.

The input parameters for the sand layer analysis using the PM4Sand model include unsaturated unit weight (γ_unsat), saturated unit weight (γ_sat), void ratio (e_init), Rayleigh damping with a value of 0.096 for Rayleigh α and 0.00079 for Rayleigh β, initial relative density (D_R), shear modulus coefficient (G ^₀), contraction rate (h_p0), atmospheric pressure (pA), maximum void ratio (e_max), minimum void ratio (e_min), bounding surface position according to ξ_R(n^b), bounding surface position according to ξξ_R(n^d), critical state friction angle (φ_CV), critical state line parameters for Q and R, and the earth pressure coefficient K ₀, with each value as shown in Table 1.

Table 2.

Parameter	Value	Unit
α	0.096	-
β	0.00079	-
	9,000	kN/m2
	9,000	kN/m2
	27,000	kN/m2
m	1	-
cref'	30	kN/m2
φ'	26	o
ψ	0	o
γ0.7	0.0007	-
	60,000	kN/m2
vur'	0.2	-
Pref	100	kN/m2
	0.5616	-
cinc'	0	kN/m2/m
yref	0	m
Rf	0.9	-
OCR	2	-

The parameters used for the clay layer utilize the HS small model. This model includes several input parameters such as secant stiffness in standard drained triaxial test tangent stiffness for primary oedometer loading power for the stress-level dependency of stiffness (m), cohesion friction angle (φ'), unloading/reloading stiffness dilatancy angle (ψ), shear strain at which G_s = 0.722 G₀ (γ_0.7) shear modulus at very small strains Poisson's ratio reference stress (P_ref), cohesion increment (c_inc'), reference coordinate (y_ref), failure ratio (R_f), and Over Consolidation Ratio (OCR), as demonstrated in Table 2.

The analysis of r_u using the finite element method with the help of Plaxis 2D is employed to provide a more detailed simulation of the pore water pressure ratio increase during an earthquake. The finite element method divides the entire soil layers into smaller elements, allowing for more precise and accurate analysis.

This study is conducted in stages, starting with an analysis using an empirical method that simplifies the calculation process. Next, a non-linear method is applied to analyze each soil layer individually, providing a more detailed soil response compared to the empirical method. Finally, the finite element method is used to divide the entire soil layers into smaller elements, yielding more detailed and accurate results.

3.1. Liquefaction Potential Analysis using the Empirical Method

Liquefaction potential analysis at four borehole points was conducted using the simplified procedure introduced by Boulanger and Idriss (2014) [13]. The earthquake parameter used in the study was a magnitude of 7.5 Mw, which occurred in Palu City on September 28, 2018, with a PGA value of 0.794g. Liquefaction is unlikely to occur at depths greater than 20 meters [23]. A study by Iwasaki in 1981 on the liquefaction potential index also limits liquefaction potential to a depth of 20 meters [24]. Based on previous research, empirical analysis was conducted up to a depth of 20 meters.

The results of soil data processing using the simplified procedure by Boulanger and Idriss (2014) [13] indicate that borehole BH 1 has liquefaction potential up to a depth of 17 m, borehole BH 2 up to a depth of 12 m, borehole BH 3 up to a depth of 15 m, and borehole BH 4 up to a depth of 13.5 m, as shown in Fig. (7).

Excess pore water pressure calculation using the Yegian and Vitelli (1981) equation employs parameters such as the liquefaction safety factor value. Thus, for layers with liquefaction potential, the limit with an r_u value > 0.8 is at the same depth as the SF value < 1, as shown in Fig. (8).

3.2. Liquefaction Potential Analysis Using Numerical Methods

In conducting liquefaction risk analysis using numerical methods, ground motion is required as the earthquake load to observe the soil's response to seismic loading. This study used ground motion from the Imperial Valley earthquake in 1979, with a magnitude of 6.5 Mw and a strike-slip fault mechanism. This ground motion is then matched to the target spectrum to approximate the ground motion conditions at the research location.

At point BH 1, data processing using DEEPSOIL v7 shows that layers with an r_u > 0.8 are located at depths of 6.5 m and 8.5 m, as shown in Fig. (9). The results also indicate an increase in Pore Water Pressure (PWP) ratio at 37.4 seconds for the 6.5 m layer and 37.2 seconds for the 8.5 m depth layer, as shown in Fig. (10). The liquefied layer is a gravelly sand layer with an average N-SPT value of 14 and 17.

At point BH 2, data processing with DEEPSOIL v7 shows an r_u > 0.8 at a depth of 7 m, as shown in Fig. (11). The results also indicate that r_u > 0.8 occurs at 21.6 seconds, as shown in Fig. (12). The liquefied layer is a silty sand to gravelly sand layer with an average N-SPT value of 14 and 17.

At point BH 3, data processing using DEEPSOIL v7 shows an r_u value > 0.8 at depths of 4.5 m and 7.5 m, as shown in Fig. (13). The results also indicate an increase in r_u > 0.8 at 20.8 seconds for the soil layer at a depth of 4.5 m and at 37.2 seconds for the soil layer at a depth of 7.5 m, as shown in Fig. (14). The liquefied layer is a silty sand to gravelly sand layer with average N-SPT values of 9 and 20.

At point BH 4, data processing with DEEPSOIL v7 shows an r_u > 0.8 at a depth of 10 m, as shown in Fig. (15). The results also indicate that r_u > 0.8 occurs at 35.7 seconds, as shown in Fig. (16). The liquefied layer is a gravelly sand layer with an average N-SPT value of 17.

The analysis results of the increase in pore water pressure ratio (r_u) using the nonlinear method show that liquefaction occurs in the silty sand to gravelly sand layers with an N-SPT range of 9 to 20. In the time history analysis, liquefaction begins after 21 seconds.

A numerical analysis was conducted using another tool, Plaxis 2D, with the PM4Sand model to compare the numerical method. The results show that borehole GA 04 has an r_u > 0.8 at 2-9 m depth, as shown in Fig. (17). The liquefied layer is a gravelly sand layer with an average N-SPT value of 9 and 17.

At borehole BH 2, liquefaction occurred down to a depth of 11 m. This is indicated by an r_u value greater than 0.8 at 3 m and 11 m depths, as demonstrated in Fig. (18). The liquefied layer consists of silty sand and gravelly sand layers with average N-SPT values of 6 and 17.

At borehole BH 3, liquefaction occurred down to a depth of 10 m. This is indicated by an r_u value greater than 0.8 at depths of 0-2 m, followed by depths of 4-5 m and 5-10 m, as shown in Fig. (19). The liquefied layer consists of silty sand and gravelly sand layers with average N-SPT values ranging from 6 to 20.

At borehole, BH 4, the Plaxis 2D calculations indicate a liquefied layer at a depth of 6.25 m. This is shown by an r_u value greater than 0.8 at a depth of 6.25 m, as illustrated in Fig. (20). The liquefied layer consists of gravelly sand with an average N-SPT value of 17.

The results obtained from the finite element method analysis indicate that the liquefied layer consists of silty sand to gravelly sand with N-SPT values ranging from 6 to 20. An interesting finding from the finite element method is that when a layer has a significant potential for liquefaction, with a thickness greater than the critical threshold, the r_u value exceeding 0.8 is not evenly distributed throughout the entire layer. Instead, the high r_u values are concentrated at the lower boundary of the layer.