Numerical Study on the Liquefaction Potential Analysis Using Constitutive Models for Sand and Silt in Mamuju, West Sulawesi

2.1. Geological and Geotechnical Condition

A two-story building located in Mamuju Regency, approximately 320 m away from the shoreline, was investigated. It is also part of the government initiative to support the rehabilitation and reconstruction of structures following the 2021 earthquake. Meanwhile, among the affected buildings was the grand mosque. Additionally, Fig. (2) shows the locations of soil investigation from the 2022 survey. For this study, soil investigation will be conducted using the Standard Penetration Test (SPT) at three borehole locations.

The geological conditions in Mamuju Regency comprised a combination of breakthrough, sedimentary, and volcanic rocks, as well as surface deposits. The sedimentary and volcanic rock formations are predominantly characterized by the following Adang (Tma) and Talaya volcanic rocks (Tmtv), as well as Mamuju formation (Tmm) and Tapalang members of Mamuju formation (Tmmt). The surface deposits in the coastal plains mainly constituted alluvial formation (Qa) and fan deposits (Qf). The study location near the coast was composed of Quaternary deposits in the form of alluvial formations (Qa) comprising boulders, crusts, pebbles, sand, silt, mud, and plant remains [21]. These quaternary deposits, subjected to weathering, were generally loose, soft, and less compact, thereby resulting in susceptibility to changes when subjected to external forces [11].

The soil layers were interpreted based on three drills and laboratory tests shown in Fig. (3). Moreover, it was reported several types of silt, clay, and sand layers characterized the study site. A silty sand layer (SM) was found at BH-02 and BH-03 with medium to dense density. The sandy silt layer (ML) at BH-01 and BH-03, ranging from 0 m to 8 m in depth, had a plasticity index (PI) value of less than 7 and was categorized as sand-like criteria [22]. However, the sandy clay (CL) layers at BH-01 and BH-02, with depths ranging from 8 m to 21 m and 0 m to 13 m, had PI values greater than 7, falling under clay-like criteria [22]. Hard soil was found starting at 19 m to 30 m depth. The study site had a potential for liquefaction due to the presence of a sand layer. Furthermore, silt and clay soils with CL, CL-ML, and ML classifications, exhibiting specific characteristics, also had the liquefaction potential [23].

Based on the result of soil investigation, the groundwater table at the study site is relatively shallow, ranging from 1.5 m to 2 m. This shallow depth is a significant factor in the potential occurrence of liquefaction, as a groundwater table less than 3 m deep had very high liquefaction susceptibility [24].

2.2. Developing Synthetic Ground Motion

Synthetic ground motion was generated because no ground motion recordings were available at the study location. This was realized using two methods, namely spectral matching and amplitude scaling [25]. In addition, spectral matching was performed because it enabled accurate matching of the target spectrum with the resulting time series. The method was also used to address discrepancies in frequency content between the initial time series and the target spectrum, which cannot be rectified through amplitude scaling [26].

The 2021 Mamuju earthquake was detected by 63 stations, and the closest to the epicenter was the Mamuju Sulawesi Transportation Agency Station (MMSN), situated 47.77 km away. Meanwhile, the Reis Watulimo Station (TDSR) was the farthest, located at 994.56 km. The study location was 4.5 km away from the nearest station, the MMSN, as shown in Fig. (4).

The ground motion data provided by BMKG was obtained only from the MMSN station, positioned in the highlands near the coast and characterized by a solid soil type with the site class SC. However, the motion was recorded in the following three components, horizontal East-West (HNE), horizontal North-South (HNN), and vertical (HNZ). The HNE component exhibited the highest maximum peak ground acceleration (PGA_max) value of 0.151 g, as shown in Fig. (5a, b) [27]. Ground motion parameters recorded at the MMSN station for the HNE component had a PGA_max, PGV_max, and PGD_max of 0.151 g, 15.38 cm/sec, and 4.47 cm, occurring at periods of 11.45 sec, 11.64 sec, and 11.42 sec, respectively.

Hanindya et al. [28] conducted a study on spectral matching using a target spectrum designed with DSHA. Meanwhile, Makrup et al. [29] performed spectral matching with a target spectrum based on the UBC 1997 code. This study adopted the Seismic Design Code of Indonesia (SNI 1726:2019). In addition, the target spectrum for the maximum earthquake considered risk-targeted (MCE_R) was assigned as the earthquake with a 2% probability of exceedance in 50 years and a return period of 2,475 years [30]. The target spectrum applicable to the study location, characterized by a soft soil (SE) site class, is shown in Fig. (6), corresponding to a PGA value of 0.619 g in bedrock [31].

There are several basic approaches to spectral matching, including frequency domain methods, frequency domain methods with Random Vibration Theory (RVT), and time domain methods [32]. Furthermore, the time domain method has been extensively studied and improved by various experts [26, 33-37]. This included adjusting the time series through the addition of wavelets. The method reportedly had good convergence properties and could maintain the non-stationary character of the initial time series [26]. Subsequently, the SeismoMatch [38] program was used to carry out spectral matching. This tool had the capability to modify the earthquake accelerogram to correspond with the response spectrum of a specific target by using the wavelet algorithm developed by [26, 37].

Based on SNI 1726:2019, the spectral matching requirement proposed that each component of the ground motion must be matched within the range of 0.8 T_lower to 1.2 T_upper. In this case, T_lower and T_upper correspond to 0.031 sec and 0.805 sec, depicting the period for spectral matching ranging from 0.025 sec to 0.966 sec. The mismatch tolerance for the two algorithms were established at 0.3.

The spectral matching results of the 2021 Mamuju earthquake accelerogram, considering the three components with the target spectrum, are shown in Figs. (7 and 8). According to Table 1, while the Al Atik and Abrahamson [26] algorithm had a smaller misfit value, it required a higher number of iterations compared to the Hancock algorithm. The mean misfit for the three records using the Al Atik and Abrahamson [26] algorithm is 2.36%, while the Hancock [37] algorithm had a mean misfit value of 3.68%. Based on the misfit tolerance value, the ground motion served as the outcome of spectral matching using the Al Atik and Abrahamson [26] algorithm on the HNE component. The time history in Fig. (9) shows the spectral matching results used to analyze the liquefaction potential at the study location, considering PGA, PGV, and PGD values of 0.478 g, 34.71 cm/sec, and 7.25 cm, respectively.

2.3. Constitutive Model

Constitutive soil models are crucial elements of analytical frameworks in the geotechnical field. These were designed to reflect the most representative conditions of the original settings, particularly for liquefied soils. The models aided in the formulation of a mathematical framework to comprehend and forecast soil response to earthquake or dynamic loads. The selection of a constitutive soil model depended on factors such as the soil type, availability of existing data, and ease of implementation. Meanwhile, two constitutive models were selected, namely PM4Sand [39] and PM4Silt [40]. Considering the soil behavior criteria proposed by [22], PM4Sand was applied to simulate soil layers consisting of sandy silt and silty sand soil types with sand-like behavior. However, PM4Silt was used to represent soil layers with sandy clays and sandy silt soil types with clay-like behavior [41-43].

2.3.1. PM4Sand Constitutive Model

PM4Sand is a plasticity model specifically designed to analyze the behavior of soils under earthquake loading, with a focus on liquefaction. It was developed based on the fundamental principles of stress-ratio controlled and critical-state-based boundary surface plasticity models initially introduced and later refined by [44] and [45], respectively.

The PLAXIS 2D program was used for modeling. Furthermore, calibration was performed using four primary and nine secondary parameters with recommended default values that can be adjusted under specific circumstances [46]. The recommended default values for the secondary parameters were proposed in a previous study [39]. Moreover, the four primary parameters included relative density (D_r), shear modulus coefficient (G₀), contraction rate parameter (h_p0), and atmospheric pressure (p_A).

The relative density (D_r) influences the dilatancy and stress-strain response characteristics, and is expressed as a ratio rather than a percentage. D_r values tend to be generally be estimated based on CPT or SPT penetration resistance. According to a previous study [14], the general form for D_r correlation based on SPT was represented by Eq. (1), where (N₁)₆₀ is the 1 atm overburden corrected SPT value, and the C_d value is determined based on the recommendation of 46.

(1)

- The shear modulus coefficient (G₀) controls the elastic (or small-strain) shear modulus. According to [39], the value of G₀ can be estimated based on the (N₁)₆₀ as in Eq. (2).

(2)

- The contraction rate parameter (h_p0) was the last to be calibrated once all other have been adjusted. This parameter was used to modify the contractiveness, enabling calibration of the model to a specific value of the Cyclic Resistance Ratio (CRR). In addition, PLAXIS 2D soil test feature enabled calibration using the CDSS test. Liquefaction was presumed to occur when the peak shear strain reaches 3% after undergoing 15 uniform loading cycles during the CDSS test [39].

- The atmospheric pressure (p_A) used was 101.3 kPa.

2.3.2. PM4Silt Constitutive Model

The previously discussed PMSand constitutive model is effective for describing sand-like behavior, and not clay-like behavior of soils. This led to the need for a model that can accurately represent the clay-like behavior of soils under earthquake loads. In another study, the PM4Silt constitutive model was developed based on the PM4Sand framework, with modifications [40] to enhance its ability to replicate the undrained monotonic and cyclic responses of silt and clay.

The PM4Silt model was also implemented into PLAXIS 2D as a user-defined model, presented in the form of a Dynamic Link Library (DLL) [47]. It required three primary parameters to calibrate and 17 secondary parameters. The PM4Silt model also recommended that default values modified in special circumstances, with the secondary parameter given [40]. Meanwhile, the three primary parameters include undrained strength at a critical state or undrained strength ratio (S_u) or (S_{u, ratio}), shear modulus coefficient (G₀), and contraction rate parameter (h_p0).

- The undrained strength in a critical state (S_u) was determined using various methods such as consolidated undrained triaxial and DSS tests, including empirical correlations. Due to data limitations in this study, the S_u value for soil types with PI ≤ 20 was estimated using empirical correlations based on corrected SPT values (N₆₀), water content, PI, and liquid limit proposed by [48] as stated in Eq. (3).

(3)

The undrained strength ratio (S_{u, ratio}) was calculated using Eq. (4), where σ'_vc represents the effective vertical stress.

(4)

- The shear modulus coefficient (G₀) was estimated using Eq. (5), where Eq. (6) was applied to determine the value of G_max. In this case, ρ represents saturated density, p depicts effective stress, and n_G has a value of 0.75 [40].

(5)

(6)

- Similar to PM4Sand, in PM4Silt, the contraction rate (h_p0) was the final parameter to be calibrated after all others had been established. The calibration of the h_p0 value for the PM4Silt model in this study was equated with the PM4Sand model using the CDSS test.

2.4. Numerical Modeling

Numerical modeling was conducted in this study using PLAXIS 2D software. The soil profile was modeled with a one-dimensional soil column with various layers in each borehole. The excess pore water pressure was determined by observing the mechanism of the soil column under cyclic loading [46, 47, 49].

Sand-like and clay-like soils were modeled using PM4Sand and PM4Silt, respectively. The Hardening Soil model was applied to soils assumed to have no possibility of liquefaction. The unavailability of bedrock depth information led to the assumption that the layer was located at the bottom of each borehole and therefore modeled with linear elastic material. The shear wave velocity value obtained for the engineering bedrock is 760 m/s [50]. The soil layer used for simulation, as well as the elevation of the groundwater table in each borehole, were adjusted based on the soil investigation results in Fig. (3).

The Rayleigh damping ratio values are α = 0.096, β = 0.00079 [46, 47], respectively. The values of G₀, D_r, S_{u, ratio} and h_p0 in each soil layer are shown in Table 2.

The h_p0 value in Table 2 was calibrated using the CDSS test in the PLAXIS 2D program. This test, which was conducted under undrained conditions, was essential for analyzing the liquefaction phenomenon. Meanwhile, the initial vertical stress value was set at 100 kPa. For the initial calibration, the number of cycles (N_M) were determined using the Eq. (7) where MSF represented the magnitude scaling factor, N_M=7.5 is 15, and coefficient b equivalent to 0.34 [46].

(7)

Correlation of (N₁)_60CS values in each soil layer to determine the cyclic stress ratio value proposed by [14] was performed to obtain the CRR value. Additionally, this was required to induce liquefaction in an equivalent number of uniform loading cycles, particularly when cyclic laboratory data are unavailable [46].

The results of the PM4Sand model calibration were validated using the curve in Fig. (10), showing the CSR value required to induce liquefaction or 3% shear strain against a uniform number of loading cycles at different D_r values. The curve conformed with the calibration performed by [39], depicting that larger D_r values require more cycles to achieve a specific CSR. It also proved that higher D_r values corresponded to greater soil resistance to liquefaction under equal loading. Additionally, the results of the calibration with the PM4Silt model were validated using the curve in Fig. (11), where the residual strength ratios of 0.45 closely matched the calibration curve performed by [40].

The numerical simulation was carried out in two stages. In the first, the K ₀ procedure was used to initialize static stresses by fabricating the model in layers and then raising the water level to the desired level. Meanwhile, at this stage, the earthquake load was not activated.

In the second stage, the earthquake dynamic load in the form of ground motion, as shown in Fig. (9), was applied along the bottom of the numerical model. Furthermore, to prevent the amplitude of the motion from doubling when it reached the free surface, only half of the input motion was required by assigning a value of 0.5 m to the x component of the specified displacement while the y direction was fixed [46, 47, 51]. The deformation boundary was selected freely in the x direction, while in the y direction, the default from PLAXIS 2D was used. Tied degrees of freedom were applied on the left and right vertical boundaries, including compliant base conditions at the bottom, to perform a 1D wave propagation analysis. The period of the dynamic analysis was based on the duration of the ground motion. In order to enable excess pore water pressure analysis, undrained A soil drainage was used. The cavitation cut-off was set at 100 kPa during the undrained dynamic analysis to reduce unrealistic excess pore pressure and stress concentrations.