Nonlinear analysis for liquefaction simulation: implication of constitutive soil model

This research is located on Lombok Island, especially in the Mandalika circuit project. The detailed location is presented in Figure 1.

FIGURE 1.

Research location in macroscale (a) and microscale (b)

Source: own work based on Ina-Geoportal maps (2025).

Regarding liquefaction susceptibility, the area is in the medium liquefaction risk zone based on the Indonesian liquefaction hazard map, which is layered in Figure 1. The site investigation was conducted using a standard penetration test (N-SPT) with five boreholes.

Sandy soil with a shallow groundwater table dominates the ground. Based on the site investigation, the five boreholes have different bedrock conditions, although all the properties are mixed with clay or silt. The groundwater level is generally 1.5 m, but in this analysis, the groundwater level is assumed to be 1 m deep for all boreholes. The bedrock condition is simulated at 20 m for BH1, BH2, and BH3, while at 10 m for BH4 and 17 m for BH5. To avoid misunderstanding, the unit weight of the soil is described as γ_S. Detailed information is presented in Table 1.

The soil properties that apply for the simulation are unit weight, shear wave velocity corresponding to the N-SPT value, and friction angle (ϕ). The relationship between shear wave velocity (V_s) and the SPT value (N₆₀) and initial effective vertical stress of the ground (σ′_v0) is described in Equation 1 (Brandenberg et al., 2010; Amalia et al., 2023): (1) lnVs=4.045+0.096lnN60+0.236lnσv0'. \[ \ln \left( {V_s } \right) = 4.045 + 0.096\ln \left( {N_{60} } \right) + 0.236\ln \left( {\sigma _{v0}^' } \right). \

This analysis computes the same layer consistency for the similar shear wave velocity values. Detailed N-values of SPT and V_s are also depicted in Figure 2.

FIGURE 2.

Results of standard penetration test (N-SPT) and shear wave velocity (V_s) test in designated boreholes: a – BH1, b – BH2, c – BH3, d – BH4, e – BH5

Source: own work.

In this study, the analysis of seismic behavior is based on site-specific response analysis through DEEPSOIL. This nonlinear calculation is applied to a one-dimensional model with input from the soil properties, using the seismic load from the Kobe earthquake ground motion. This is due to lack of real ground motion recorded at the location, although the area is highly impacted as a seismically prone site. The Kobe earthquake is embedded in DEEPSOIL and has been applied in several liquefaction cases in Indonesia, such as the 2018 Palu earthquake (Jalil et al., 2021; Hori et al., 2023) and the Yogyakarta airport project (Fitri & Pramana, 2025).

The dynamic properties of sand, such as soil properties, are estimated from Darendeli’s empirical model (Stokoe et al., 2004). This approach is calculated to simulate cyclic loading behavior, such as loading and unloading, and yields the stress–strain loop connection. Moreover, Phillips and Hashash (2009) employed the damping factor (D) and modulus reduction (G/G_max) for calibration and minimizing discrepancies between the targeted and resulting curves.

In general, numerous soil constitutive models for dynamic analysis have been proposed to account for the solid–liquid coupling phenomena in saturated sandy soils exposed to intense shaking, such as applying the effective stress method (Chiaradonna et al., 2018; Wang et al., 2014, 2021; Boccieri et al., 2024; Chen et al., 2024). Moreover, several constitutive models are available in DEEPSOIL for assessing the dynamic behavior of the ground. However, this research will apply two constitutive models for liquefaction analysis: General Quadratic/Hyperbolic (GQ/H) and Modified Kondner–Zelasko (MKZ). The detailed analysis, which focuses on shear strain and stress, is distinguished by the following approach.

The MKZ model is an upgrade from the hyperbolic approach model provided by Matasović and Vucetic (1993). This analysis applies β and s for fitting the shape of the curve where G₀ is the initial shear modulus, shear strength is τ, and shear strain is γ. The process ensures a stress–strain curve that complies with a user-specified shear strength at a specified strain while matching the analysis goal modulus-reduction (and related damping) curve over a selected strain range: (2) τ=G0γ1+βγγrS. \[ \tau = \frac{{G_0 \gamma }} {{1 + \beta \left( {\frac{\gamma } {{\gamma _r }}} \right)^S }}. \

While in DEEPSOIL, the model is augmented by γ_r as confining pressure according to Hashash’s and Park’s theory (2001). The value is calculated by dividing the vertical effective stress (σ′_v) by the reference stress (σ_ref): (3) γr=σv'σrefb. \[ \gamma _r = \left( {\frac{{\sigma _v^' }} {{\sigma _{{\text{ref}}} }}} \right)^b . \

As an extended method from MKZ, the GQ/H method is a further approach that combines ultimate shear strength and modulus reduction factor with θ_τ as a parameter for fitting the curve, and maximum shear stress (τ_max): (4) τ=τmax1θτ1+γγr−1+γγr2−4θτγγr, \[ \tau = \tau _{\max } \left[ {\frac{1} {{\theta _\tau }}\left\{ {1 + \left( {\frac{\gamma } {{\gamma _r }}} \right) - \sqrt {\left\{ {1 + \frac{\gamma } {{\gamma _r }}} \right\}^2 - 4\theta _\tau \frac{\gamma } {{\gamma _r }}} } \right\}} \right], \ where θ_τ is calculated by the equation: (5) θτ=θ1+θ2xθ4x(γγr)θ5θ3θSθ4x(γγr)θ5. \[ \theta _\tau = \theta _1 + \theta _2 x\frac{{\theta _4 x\left( {\frac{\gamma } {{\gamma _r }}} \right)^{\theta _5 } }} {{\theta _3^{\theta _S } \theta _4 x\left( {\frac{\gamma } {{\gamma _r }}} \right)^{\theta _5 } }}. \

DEEPSOIL utilizes the MRDF (Phillips & Hashash, 2009) technique that improves the fit to damping behavior and modulus reduction. The correction factor is calculated as Equation 5. Parameters θ₁, θ₄, and θ_τ are employed to alter the large strain values according to the designated large strain shear strength and to maintain the modulus reduction curves derived from reference studies as much as possible. After generating the layered domain, the GQ/H curve fitting procedure generates the shear strength correction for the soil shear strength–shear strain graph and delivers the parameters θ₁ through θ₅. All further detailed descriptions are provided in the DEEPSOIL user manual (DEEPSOIL, 2024).

In order to improve the simulation of cyclic strain calculation, the GQ/H model further specifies a thorough hysteresis formulation for including unloading and reloading approaches. The reaction to hysteretic stress is provided by Equation 6. The τ_rev and γ_rev are the point of reversal at stress and strain point, and τ acts as a reduction factor, and Gγ_max is the tangent modulus that associated with the particular γ_max: (6) τ=Fγmax−τmaxθτ1+γ−γrev2γr−−1+γ−γrev2γr2−4θrγ−γrev2γr−Gγmax(γ−γrev)+Gγmax(γ−γrev)+τrev. \[ \begin{gathered} \tau = F\left( {\gamma _{\max } } \right)\left[ { - \frac{{\tau _{\max } }} {{\theta _\tau }}\left\{ {1 + \frac{{\gamma - \gamma _{{\text{rev}}} }} {{2\gamma _r }}} \right\} - } \right. \hfill \\ \left. { - \sqrt {1 + \left( {\frac{{\gamma - \gamma _{{\text{rev}}} }} {{2\gamma _r }}} \right)^2 {-4{{\theta _r }}} \left( {\frac{{\gamma - \gamma _{{\text{rev}}} }} {{2\gamma _r }}} \right) - G\gamma _{\max } (\gamma - \gamma _{{\text{rev}}} )} } \right] + G\gamma _{\max } (\gamma - \gamma _{{\text{rev}}} ) + \tau _{{\text{rev}}} . \hfill \\ \end{gathered} \

The pore water pressure generation is calculated by Equation 7 and considers the Dobry–Matasovic (Dobry et al., 1982) model approaches when it is called normalized excess pore pressure (r_u) as u_n, a parameter for fitting the curve (p, s, and F), f is the dimensionality factor, N_c is the number of cycles, the shear strain reversal value is γ_c, and γ_tvp is the threshold for shear strain: (7) un=pfNcFγC−γtvpS1+fNCFγC−γtvpS. \[ u_n = \frac{{pfN_c F\left( {\gamma _C - \gamma _{tvp} } \right)^S }} {{1 + fN_C F\left( {\gamma _C - \gamma _{tvp} } \right)^S }}. \

The normalized excess pore water pressure ratio is also defined as the u_n parameter in the following equation: (8) ru=u′ σv', \[ r_u = \frac{{u'}} {{\sigma _v^' }}, \ where r_u is the pore water pressure ratio, u′ is excess pore water pressure, and σ′_v is the effective vertical stress.

The output for r_u from each layer is presented in Figure 3 for all layers. All the results are included in two constitutive models, the MKZ and GQ/H models. The Darendeli model (Stokoe et al., n.d.) was employed to evaluate the modulus reduction and damping contours in this investigation for each layer, with the coefficient at rest (K₀) applied from Jacky’s formula (Diaz-Segura, 2016). Different behaviors are depicted in all grounds since the model is on layered ground with varying soil properties. The r_u value occurs after 1 m depth because the modeling of the groundwater table is below 1 m depth.

FIGURE 3.

Pore water pressure generation from DEEPSOIL analysis: a – BH1, b – BH2, c – BH3, d – BH4, e – BH5

Source: own work.

In general, the MKZ and GQ/H results describe almost similar values in the several maximum ranges of r_u (BH1, BH2, BH3, BH5), although they have differences around the r_u under 0.5 (BH2, BH3, BH5). In addition, all the boreholes depict the deeper layer of the ground, and the values of the two constitutive models tend to be similar but have slight differences. The maximum value of r_u is in BH1 at 18 m, BH2 at 16 m, BH3 at 12 m, BH4 at 4 m, and BH5 at 10 m.

This study only focuses on modeling liquefaction simulation based on pore water pressure (r_u) generation without any supporting data from the physical laboratory tests. Based on the centrifuge analysis (Groholski et al., 2016), the MKZ typically generates the r_u effectively at a moderate level of liquefaction but lacks prediction accuracy under full liquefaction conditions. Furthermore, Raza and Ahmad implied that r_u ≥ 0.9 indicates complete liquefaction condition while marginal liquefaction is 0.8 ≤ r_u ≤ 0.9 (Raza & Ahmad, 2024). Based on the results, the smaller the value of r_u from the MKZ approach, the larger the gap between the two constitutive models. This means that in the no liquefaction to marginal liquefaction conditions, the MKZ predicts less safe results than GQ/H. This threshold also combines with the strain value to predict the real phase of the liquefied soil layer more accurately.

GQ/H ignores the multiple, occasionally subjective processes frequently required to calibrate MKZ for strength under high stresses in cyclic loading. This condition leads to a more straightforward and transparent model calibration and more dependable output for crucial liquefaction situations. In addition, some studies also declare that similar degradation indicators for modulus and strength are included in the GQ/H model, which can more accurately represent large-strain shear strength than MKZ (Mei, 2018; Di Buccio & Pagliaroli, 2020).

Due to the large gap between the two constitutive models in the r_u output, especially at shallow depths, the hysteresis loops (relation between stress and strain) are focused on the −5 m depth in all boreholes. This depth is chosen because the r_u value tends to be under 0.8 at this particular depth except in BH4, which means it is not in the full liquefaction condition. All the results are shown in Figure 4.

FIGURE 4.

Shear stress and shear strain relation output in depth of −5 m: a – BH1, b – BH2, c – BH3, d – BH4, e – BH5

Source: own work.

All the boreholes present a similar behavior in that strain results from MKZ tend to be smaller than GQ/H, although the values were reversed under conditions of full liquefaction (BH4). This output is appropriate for the study that proved MKZ underestimates the rate of shear strength degradation at high pore water pressures and the loss of cyclic stiffness, but it may be overestimated at large strains. Meanwhile, the GQ/H approach improves estimates of how soil deteriorates before and after the first phase of liquefaction by providing more realistic hysteresis that appropriately represents the significant decrease in cyclic resistance and energy dissipation (Shamsher et al., 2025).

The liquefaction analysis is not only about predicting the behavior of the excess pore water pressure ratio (r_u), but also the output for shear and strain of the modeling. The combination of the failure process in the soil due to strain stress, especially in shear force, is considered. The maximum r_u value in each borehole with the MKZ model is described in Figure 5.

FIGURE 5.

Relationship among excess pore water pressure ratio (r_u) and strain in the DEEPSOIL model: a – BH1, b – BH2, c – BH3, d – BH4, e – BH5

Source: own work.

The strain output shows that the peak of r_u does not always produce the maximum strain ratio. However, the pore water pressure of the MKZ model is lower than of GQ/H because the degradation model considers a reduction in soil material stiffness at each time step. For a more accurate simulation, Chiaradonna et al. (2015) suggested that the recalculation of stiffness and soil strength during the post-cyclic consolidation phase should be considered, as should the dissipation and redistribution of excess pore pressure inside a soil deposit.

The study covers the numerical analysis through modeling soil properties in Indonesia to obtain the liquefaction behavior of the ground. Although the cyclic load does not fully represent the actual seismic motion, the behavior of the softening soil is expected to be the basic threshold for disaster management in the Lombok region. The two constitutive models offer benefits and drawbacks to the output that a geotechnical expert should consider in engineering practices. In addition, laboratory physical testing, such as centrifuge and 1 g shaking table tests, is necessary for future study because these tests can replicate the soil model and specific site profiles, allowing for advanced validation of the MKZ and GQ/H constitutive models.