E-ISSN: 2528-388X
P-ISSN: 0213-762X
INERSIA
Vol.
No.
May 2026 Influence of Constitutive Models on Cyclic Pore Pressure and Liquefaction Assessment in Srandakan.
Bantul.
Indonesia Anisa Nur Amalina*.
Lalu Makrup and Elvis Saputra Department of Civil Engineering.
Universitas Islam Indonesia.
Yogyakarta, 555834.
Indonesia
ABSTRACT
Keywords:
Hardening Soil
Mohr-Coloumb
PLAXIS 2D
PSHA
PM4SAND
Liquefaction potential in earthquake-hazard areas is correlated with earthquake-induced cyclic This condition leads to an important selection in the constitutive soil model in numerical site response analysis.
This study aims to evaluate the effect of different constitutive soil models on cycling pore pressure by using a numerical model.
The case was taken in Srandakan.
Bantul.
Indonesia, which consists of sand and clay soil.
A representative ground motion was derived from a probabilistic seismic hazard analysis, with PGA OO 0.
31, in a one-dimensional dynamic analysis using PLAXIS 2D.
Three constitutive models.
MohrAeCoulomb (MC).
Hardening Soil with small-strain stiffness (HSsmal.
, and PM4Sand, were compared in terms of excess pore pressure development and pore pressure ratio (R.
The results showed that MC neglected a pore-pressure accumulation, while HSsmall captured only a limited nonlinear response.
Meanwhile.
PM4Sand predicted more significant cyclic pore pressure changes in sandy layers due to its state-dependent plasticity formulation.
Nevertheless.
Ru values remained well within acceptable limits at all depths, indicating that the applied seismic excitation did not trigger liquefaction.
The findings demonstrated that model selection substantially affected pore pressure response in seismic liquefaction analysis.
This is an open access article under the CCAeBYlicense.
Introduction Liquefaction is a phenomenon in which saturated, cohesionless soils .
ypically sandy soil.
lose strength and stiffness due to an increase in pore-water pressure, often triggered by seismic activity.
This increase in pore water pressure reduces the effective stress in the soil, causing it to behave like a liquid .
, .
, .
, .
, .
This process poses serious risks to structures, potentially leading to foundation failure and collapse.
Therefore, assessing liquefaction potential is essential, particularly in areas planned for development.
Earthquakes, through the propagation of seismic waves, are the primary triggers of Developing reliable models to evaluate liquefaction potential is thus critical for seismic risk Earthquake ground motions can be characterized using either deterministic or probabilistic approaches.
Probabilistic Seismic Hazard Analysis (PSHA) has become the standard tool for estimating site-specific *Corresponding author.
E-mail: anisa.
amalina@uii.
https://doi.
org/10.
21831/inersia.
Received 20 February 2026.
Revised 13 April 2026.
Accepted 13 April 2026 Available online 01 May 2026 ground motions.
In Indonesia and Southeast Asia.
PSHA has been extensively applied to produce seismic hazard and spectral maps and support microzonation studies and revisions of the national seismic building code .
, .
, .
, .
, .
, .
Ground motion modification techniques are essential in geotechnical analyses to evaluate the impact of seismic loading on soil and foundation systems, particularly for assessing liquefaction potential.
These techniques involve altering recorded ground motions to match a target response spectrum, which helps in reducing the number of ground motions required for analysis and ensures consistency in seismic evaluations .
, .
, .
Liquefaction potential has likewise been widely studied in Indonesia and abroad.
Previous works have evaluated liquefaction hazards using SPT data, state parameter analysis, and ground motions derived from probabilistic These studies consistently highlight the importance of integrating seismic hazard analysis with INERSIA.
Vol.
No.
May 2026 Anisa Nur Amalina, et.
geotechnical site response to better capture liquefaction risk .
, .
, .
2D in Figure 1.
The analysis focuses on differences in Ru evolution and the distribution of excess pore pressure to evaluate the sensitivity of liquefaction assessment to constitutive model selection.
By acknowledging model limitations, liquefaction-induced damage can be mitigated through remediation methods that depend upon a thorough evaluation of liquefaction.
Also, this study offers insights to improve the development of constitutive models for liquefaction analysis.
Site response analysis (SRA) is a critical tool in seismic design, particularly for understanding how local site conditions affect ground response during earthquakes .
This method is robust to manual comparisons of Cyclic Resistance Ratio (CRR) and Cyclic Stress Ratio (CSR), as it enables a more detailed, dynamic analysis of soil behavior, including the development of pore pressure under cyclic loading.
Methods Site response analysis using finite element platforms such as PLAXIS can help model the effect of cyclic loading on Moreover, it offers several constitutive models that offer insights into soil behaviour under various conditions .
Soil behavior in dynamic analysis can be represented by various constitutive models of varying The MohrAeCoulomb (MC) model is simple but does not explicitly capture cyclic pore pressure generation .
The Hardening Soil model with smallstrain stiffness (HSsmal.
improves the modeling of stiffness degradation under cyclic loading .
, whereas advanced critical-state-based models, such as PM4Sand, are specifically formulated to simulate pore pressure accumulation and cyclic mobility in sands .
Differences in model formulation may significantly influence predicted liquefaction indicators, particularly the pore pressure ratio (R.
This study evaluates the influence of constitutive soil models on liquefaction prediction through dynamic site response analysis.
Seismic input motions were selected based on a Probabilistic Seismic Hazard Analysis (PSHA) for Bantul.
Indonesia.
The PSHA results were used solely to define the input ground motion level and were not analyzed further in this study.
1 Probabilistic Seismic Hazard Analysis The Probabilistic Seismic Hazard Analysis (PSHA) explicitly accounted for uncertainties in earthquake occurrence, magnitude, location, and ground motion The magnitude distribution, describing the recurrence of earthquakes, is typically modelled with the GutenbergAeRichter relationship law .
To avoid unrealistic values, the magnitude range is truncated between a minimum magnitude .
CA) and an upper-bound magnitude .
A).
In most hazard analyses, mCA ranges between 3 and 5 .
This study compared MC.
HSsmall, and PM4Sand within a one-dimensional soil column subjected to PSHA-based input motion in a liquefaction-prone area using PLAXIS Figure 1.
Srandakan.
Bantul.
Yogyakarta.
Indonesia .
Anisa Nur Amalina, et.
INERSIA Vol.
No.
May 2026 In this study, the recurrence law for earthquake magnitude followed the truncated exponential GutenbergAeRichter model .
and its modification by Youngs and Coppersmith .
, known as the characteristic earthquake recurrence law.
This combined approach accounted for small-to-moderate events using the GutenbergAeRichter relation and large events using the characteristic model, where the probability distribution becomes uniform as it approaches the maximum magnitude.
spectrally matched to the Uniform Hazard Spectrum (UHS) obtained from the PSHA at the bedrock level.
Spectral matching was performed within the relevant engineering period range to ensure consistency with the target hazard spectrum while preserving realistic temporal The resulting modified acceleration time history had a peak ground acceleration (PGA) of approximately 0.
31 g, a total duration of about 30 s, and a time increment .
of 005 s.
The matched acceleration record was applied directly as a prescribed base acceleration input in PLAXIS 2D under compliant base boundary conditions, without additional scaling .
caling factor = 1.
The distance distribution described the probability of a rupture occurring at different distances from the site.
PSHA, experts generally use the rupture width equal to the rupture length .
This is commonly estimated using the relative rupture area on the fault plane, assuming equal likelihood of rupture along the fault length.
It should be noted that detailed PSHA derivations and spectral-matching procedures are not elaborated on in this paper, as the primary objective of this research is not hazard quantification itself but rather the comparative evaluation of different constitutive soil models in predicting site response and liquefaction behaviour under identical seismic loading conditions.
Ground-motion prediction equations (GMPE.
were used to evaluate this conditional probability across different tectonic environments affecting the study area.
Hazard deaggregation was performed to identify the dominant earthquake scenario contributing to the target hazard level.
The results indicate that the Opak Fault provides the largest contribution to the seismic hazard at the study The controlling scenario corresponds to a moment magnitude of Mw OO 6.
32 and a hypocentral distance of approximately 15.
5 km.
This magnitudeAe distance pair was subsequently used to guide the selection of a representative strong ground motion record.
2 Ground Motion Input The acceleration record was imported into PLAXIS 2D in tabulated format and applied directly as a prescribed base acceleration under compliant base boundary conditions.
No additional scaling was applied in the numerical model .
caling factor = 1.
The compliant base formulation allows wave transmission from the underlying bedrock while minimizing artificial reflection at the model The ground motion input is shown in Figure 2.
A representative strong ground motion record from the 1987 Whittier Narrows earthquake (Mw 5.
was selected based on similarity in magnitude and source-to-site distance to the governing seismic scenario.
The record was Figure 2.
Ground motion input INERSIA Vol.
No.
May 2026 Anisa Nur Amalina, et.
developed to simulate cyclic mobility and pore pressure accumulation in sandy soils.
Each constitutive model has their own input parameters.
3 Numerical Model Configuration This study employed a one-dimensional soil column approach to evaluate seismic site response and liquefaction potential.
The theoretical basis follows the classical one-dimensional shear wave propagation concept introduced by Bardet and Tobita .
, in which vertically propagating shear waves travel from bedrock to the ground surface within a horizontally layered soil The governing equation for one-dimensional wave propagation in a viscoelastic medium can be expressed as Equation .
C 2d Cd CA Ct Cz Ct All soil parameters were estimated based on available SPT data and empirical correlations.
For sandy layers, relative density and friction angle (I') were correlated with NSPT values using the empirical relationships proposed by .
Meanwhile, for the clay soil, the correlations were based on .
The basic input parameters for the MC model can be seen in Table 1.
These parameters consist of unsaturated and saturated unit weight .
nsat and sa.
, cohesion .
, dilatancy angle (O), modulus of elasticity (E), and PoissonAos ratio ().
Several parameters in MC were also used as input parameters for .
where A is the soil mass density, d is the horizontal displacement, z is the depth coordinate, t is time.
E is shear stress, and represents the damping coefficient.
Although the physical problem corresponds to onedimensional wave propagation, the numerical implementation was carried out in PLAXIS 2D under plane strain.
A narrow soil column was modeled with tied lateral degrees of freedom to suppress horizontal deformation incompatibility, thereby enforcing a condition equivalent to one-dimensional shear wave 4 Geometry and Soil Profile The analysis was conducted based on soil investigations.
The average N-SPT and soil classification are shown in Figure 3.
Based on SPT data, the soil classifications consisted of sand, silt, and clay.
N-SPT value ranged between 25 to 20 at 0 Ae 12 m and below 10 at the depths of 14 m and 18 m.
The dense sand was found below 20 m with N-SPT > 60.
The groundwater table was assumed at This soil profile was input to PLAXIS 2D, and the lateral boundaries were modelled using tied degrees of freedom to enforce one-dimensional wave propagation in Figure 4.
Figure 3.
Soil stratigraphy in the study area 5 Constitutive Models Three constitutive models available in PLAXIS 2D were adopted to evaluate their influence on liquefaction-related The first is the MohrAeCoulomb (MC) model, a linear-elastic, perfectly plastic model with constant The Hardening Soil model with small-strain stiffness (HSsmal.
is a nonlinear model accounting for stress-dependent stiffness and small-strain shear modulus degradation, and PM4Sand is a stress-ratio-controlled, critical-state-based Figure 4.
Soil stratigraphy in the study area Anisa Nur Amalina, et.
INERSIA.
Vol.
No.
May 2026 Table 1.
Input parameters for Mohr-Coulomb For HSsmall, the reference stiffness parameters (ECICA-ref.
Eoed-ref.
Eur-re.
were estimated using correlations with NSPT and overburden stress as suggested in .
The small-strain shear modulus GCA-ref was defined based on estimated Vs.
Silt .
Ae16.
Clay .
Ae20.
Sand .
Ae Model Undrained (A) Undrained (A) Undrained (A) Drained 79,000 90,000 43,000 180,000 Table 2.
Input parameters for Hardening Soil with small-strain Sand Silt Clay Sand Layer .
Ae6.
Ae16.
Ae20.
Ae28.
E50-ref For PM4Sand, this soil model is only used for the sandy soil, the cohesive layers were input as the HSsmall model.
The input parameters were calibrated using relative density and the typical ranges recommended by Boulanger and Ziotopoulou .
The contraction rate parameter and fabric parameters were selected within recommended bounds for medium-dense to dense sand conditions, consistent with the measured NSPT values.
The basic input parameters for the PM4Sand model can be seen in Table 3.
The input parameters include relative density (D.
, shear modulus coefficient (G.
, contraction rate parameter .
, minimum and maximum void ratios .
min and ema.
, critical state friction angle (Ic.
, dilatancy parameter (Q), and dilatancy surface parameter (R).
Eur-ref .
Eoed-ref .
G0ref .
79,000
90,000
43,000
180,000
237,000
270,000
129,000
540,000
79,000
90,000
43,000
180,000
350,000
300,000
120,000
600,000
Table 3.
Input parameters for PM4Sand Parameter Sand .
Ae6.
Icv (A) 6 Liquefaction Assessment Parameters Liquefaction potential was evaluated using timedependent excess pore water pressure and pore pressure ratio (R.
The pore pressure ratio was defined as .
A yuayc0 Sand .
Ae6.
N/m.
N/m.
c' .
I' (A) O (A) E .
For silt and clay layers, the HSsmall model was adopted under undrained conditions.
Undrained shear strength was estimated from empirical SPT correlations .
, and stiffness parameters were adjusted to reflect stressdependent behaviour.
The basic input parameters for the HSsmall model are shown in Table 2.
The input parameters include secant stiffness at 50% strength (E50re.
, unloading/reloading stiffness (Eur-re.
, eodometer stiffness (Eoed-re.
, small-strain shear modulus (G0re.
, reference shear strain .
, and stress-dependency exponent .
ycIyc .
= Layer Sand .
Ae28.
Additional response parameters evaluated include excess pore pressure evolution .
u vs dynamic tim.
Representative evaluation depths were selected at 3 m, 11 m, 18 m, and 25 m to enable consistent comparison of model behaviour across layers.
where iu.
is the excess pore pressure at time t, and EAv0 is the initial vertical effective stress obtained from the geostatic initial phase.
This equation follows the effective stress principle and liquefaction framework, which states that liquefaction occurs when excess pore pressure approaches the initial effective vertical stress (Ru OO .
When Ru remains below 1, the soil does not reach liquefaction conditions.
In cases where excess pore pressure becomes negative, the soil exhibits dilative behavior, reflected by a reduction in pore pressure during cyclic loading .
Results and Discussion 1 Time-Dependent Excess Pore Pressure Response The evolution of excess pore pressure .
mostly showed a negative response in both the shallow and deep sand layers, particularly at 3 m and 25 m.
This behaviour indicated a dilative behaviour.
The relatively small INERSIA.
Vol.
No.
May 2026 Anisa Nur Amalina, et.
magnitude of iu across all depths suggests that full liquefaction conditions were not reached under the applied ground motion intensity.
The results are shown in Figure 5 to Figure 8.
The excess pore pressure response .
showed predominantly negative values in the sand layers, which are considered dense and very dense sand.
This behaviour indicates dilative behaviour under cyclic loading.
Dense sands subjected to moderate seismic excitation can develop negative pore pressures due to dilation effects, and the shear modulus decreases slowly.
This behaviour enhances the undrained shear strength and reduces the likelihood of liquefaction.
HSsmall PM4Sand Dynamic time .
Figure 5.
iu vs Dynamic time for layer 1-Sand .
m-dept.
Across all layers, this behaviour is commonly captured by PM4Sand, which showed a gradual development of a negative excess pore pressure of nearly -5,5 kPa at the end of dynamic loading.
Even the excess pore pressure reached -18 kPa, indicating strong dilative behaviour due to the soil's stiffness.
On the other hand, the other models present limited behaviour of the development of excess pore pressure.
Hssmall PM4Sand iu .
For the cohesive layers (Figure 6 and Figure .
, the soil in the PM4Sand scenarios was still set to HSsmall.
For the clay layer, especially, the response to excess pore pressure differs from that of the observed sandy layers.
The MC model predicted a rapid drop of iu at the beginning of dynamic time, in contrast, the HSsmall model, for both scenarios, produces a minor change in excess pore The MC model assumes purely elastic behavior before failure, which can lead to an overestimation of pore pressure fluctuations in cohesive soils.
This is because the model does not account for pre-failure plastic deformations that can generate excess pore pressure under undrained conditions .
Dynamic time .
Figure 6.
iu vs Dynamic time for layer 2-Silt .
m-dept.
Furthermore, there is an inconsistency in the MC for modelling excess pore pressure in cohesive soil.
This condition relates to the concept that MC cannot predict pore pressure generation.
For the HSsmall, the behaviour difference occurred due to the main parameter, 0.
7 and The silt layer has a smaller 0.
7 value than the clay layer, thus its stiffness decreased faster and developed greater cyclic strain.
This led to a larger iu response than in the clay layer .
HSsmall PM4Sand Dynamic time .
Figure 7.
iu vs Dynamic time for layer 3-Clay .
m-dept.
Anisa Nur Amalina, et.
INERSIA.
Vol.
No.
May 2026 iu .
HSsmall PM4Sand HSsmall PM4Sand Dynamic time .
Dynamic time .
Figure 8.
iu vs Dynamic time for layer 4-Sand .
m-dept.
Figure 9.
Ru vs Dynamic time for layer 1-Sand .
m-dept.
2 Generation of Pore Pressure Ratio (R.
Time-dependent evolution of Ru was used to demonstrate the liquefaction potential along the profile model under the applied ground motion of 0.
31 g.
In general.
Ru values 2 indicate that full liquefaction is not triggered under dynamic motion.
In addition, iu is negative in several layers.
Ru appears positive in the PLAXIS output due to the compression-negative stress convention adopted in the software.
Because both iu and EAv0 are negative under compressive conditions, their ratio becomes positive.
The results are shown in Figures 9 to Figure 12.
Hssmall PM4Sand Dynamic time .
Figure 10.
Ru vs Dynamic time for Layer 2-Silt .
m-dept.
In the shallow sand layer .
, the PM4Sand model predicted an increase in Ru up to approximately 0.
indicating moderate pore pressure accumulation under cyclic loading.
A similar but lower trend is observed in the deeper sand layer .
, where Ru reached about 0.
In contrast, using the HSsmall model, the silt layer .
had limited Ru development, while the clay layer .
showed minor fluctuations with relatively small Ru values throughout the shaking duration.
HSsmall PM4Sand In the deep sand layer .
Ru reached only around 0.
under the PM4Sand model.
Even though the excess pore pressure in Figure 7 is relatively large in absolute terms, the initial effective stress at greater depth will also reduce the pore pressure ratio.
This result highlights that the liquefaction potential strongly depends on the confining stress level .
Dynamic time .
Figure 11.
Ru vs Dynamic time for layer 3-Clay .
m-dept.
INERSIA.
Vol.
No.
May 2026 Anisa Nur Amalina, et.
However, the maximum pore pressure ratio (R.
remained well below unity at all investigated depths, indicating that full liquefaction was not triggered under the applied seismic excitation.
Overall, the results demonstrate that constitutive model selection plays a crucial role in liquefaction-oriented site response analysis, as it directly affects the predicted magnitude and distribution of cyclic pore pressure generation.
HSsmall PM4Sand Dynamic time .
References