By considering the characteristics of the time-frequency curve of ground acceleration at liquefaction sites, this research introduces the concept of a step function, defines an error function, and establishes criteria for identifying site liquefaction time. A novel method for site liquefaction time identification is proposed. One-dimensional and two-dimensional non-linear site response analysis models are established, and real liquefaction cases are studied to compare and analyze the two methods. The results demonstrate that the error function effectively identifies the moments when the frequency of the time-frequency curve of ground acceleration undergoes rapid changes, enabling the identification of site liquefaction time.
- Site liquefaction time depends on the site response.
- Time-frequency distribution of ground surface motion could be obtained based on HHT.
- Liquefaction time could be identified based on the phase shift or rapid change of frequency contents.
The deformation of horizontally liquefied soil can be divided into three stages: small-magnitude seismic deformation before liquefaction, moderate-large deformation induced by seismic shaking after liquefaction, and permanent deformation after the earthquake . Therefore, identifying the liquefaction time is important for evaluating the dynamic response of the site .
Several studies have been conducted to obtain the liquefaction time of the liquefied site, Yang and Kavazajian  used FLAC 2D to construct detailed profiles of five layers and obtained the liquefaction time of the site using the PM4sand constitutive model. Özener et al.  used power spectra, combined with Stockwell spectrograms to predict the liquefaction time. Zhang et al.  developed a liquefaction time prediction model by training a CNN using time-frequency spectral. These studies [3-5] have made significant contribution to the prediction of liquefaction time, while the method used in  requires highly detailed geotechnical investigation profile, the method proposed in  is relatively complicated, and the accuracy of the method in  depends on the number of the seismic motions used in the training of convolutional neural network.
In this research, based on the distribution of seismic records in the time and frequency domains, a new method for predicting liquefaction time is proposed. The real liquefaction time of the site is taken as the benchmark for comparison, and the accuracy of the proposed method as well as non-linear site response analysis method in determining liquefaction time are analyzed.
2. Site response analysis
The present study focuses on the analysis of the Treasure Island site during the 1989 Loma Prieta earthquake ( 6.9). This earthquake was characterized by various liquefaction phenomena, such as sand boiling and ground deformation . Nonlinear analysis models of the site were established using DEEPSOIL and PLAXIS 2D. The deconvolution results of the bedrock seismic waves were used as inputs for the seismic waves. The liquefaction time of the site was determined based on the changes in effective stress.
2.1. Site liquefaction case
The soil profile and material parameters of the site are shown in Fig. 1 and Table 1 . During the earthquake, bedrock motion was obtained near Treasure Island in Yerba Buena Island and deconvolution was performed on it. Based on the deconvolution analysis results, the horizontal ground motion records with rich energy content were selected using Arias intensity as inputs for DEEPSOIL and PLAXIS 2D analyses, as shown in Fig. 2.
Fig. 1The soil profile of Treasure Island 
Fig. 2Input earthquake records
Table 1Basic size and style requirements 
Young Bay mud
Old Bay mud
dry density (kg/m3)
Angle of internal friction (°)
Shear wave velocity (m/s)
DEEPSOIL is a one-dimensional site response analysis program developed by the University of Illinois, USA. It offers three types of site response analysis methods: linear, equivalent linear and non-linear . According to reference , the sand fill material is identified as a liquefiable soil layer. To investigate the liquefaction time of the site, a soil column shown in Fig. 3 is constructed based on the blue column in Fig. 1. The stress-strain models for all layers are set as the MKZ model. The excess pore pressure model for the sand fill material below the groundwater level adopts the Vucetic/Dobry model, while the Holocene and Pleistocene Bay Mud layers use the Matasovic excess pore pressure model. The bedrock is modeled as an elastic half-space with of 1500 m/s, a unit weight of 25 kN/m3, and a damping ratio of 5 %. The dynamic shear modulus and damping ratio are determined based on experimental results from reference . The monitoring point for excess pore pressure is located at point A, at an elevation of –11.6 m, at the bottom of the sand fill material layer.
PLAXIS 2D is selected as the software for two-dimensional non-linear site response analysis. Based on the soil parameters of Treasure Island shown in Fig. 1 and Table 1, and referring to the research findings from reference , the following parameters are adopted. The Mohr-Coulomb model is used to generate initial stresses in the soil layers based on the available data. The UBC3D-PLM sand liquefaction model is applied to the liquefiable soil layer, while the remaining non-liquefiable soil layers are modeled using small-strain soil hardening models. The UBC3D-PLM model includes several key parameters, peak friction angle (), friction angle at constant volume (), effective cohesion () and the normalized SPT blow count [()60]. Among them, the value of ()60 is determined through back-analysis based on Eq. (1) . For Sandy fill, 31.29°, 30.45°, 1kPa and ()60 = 8.45:
Fig. 3The parameters of DEEPSOIL soil column model
To ensure that the shear strength at the lateral boundaries is not completely lost, as this could lead to excessive deformations, a Mohr-Coulomb constitutive model is applied with a 1-meter-wide drained boundary at the lateral boundaries. Additionally, to mitigate the influence of high-frequency components in the seismic motion, the Rayleigh damping method proposed in reference  is implemented. Considering the propagation theory of seismic waves in an elastic half-space, a horizontal displacement component of 0.5 meters is assigned.
Based on the above considerations, a numerical model is established with 6651 elements and 54510 nodes. The grid of model is depicted in Fig. 4.
Fig. 4Finite element analysis model
The input for DEEPSOIL and PLAXIS 2D is the bedrock seismic wave shown in Fig. 2, which is used for non-linear site response analysis. The excess pore water pressure time history at location A is obtained and depicted in Fig. 5. The black solid line represents the calculation results from DEEPSOIL, while the red dashed line represents the results from PLAXIS 2D. The liquefaction time predicted by DEEPSOIL is 11.5 s, while the liquefaction time predicted by PLAXIS 2D is 11.9 s.
3. Prediction method of site liquefaction time
Earthquake signals are non-stationary time-domain signals, and the frequency components of the ground motion will change due to liquefaction. If the phase shift of the motions could be captured, then the corresponding liquefaction time could be obtained. Thus, this paper proposed a method for predicting the site liquefaction time as follows.
1. Select a seismic record with abundant energy in the horizontal direction based on the Arias intensity. Apply the Hilbert-Huang Transform (HHT) to process the ground acceleration time history. From the Hilbert spectrum, extract the excellent frequency curve by identifying the frequency corresponding to the maximum energy at each time instant. The specific physical interpretation of the excellent frequency is discussed in detail in reference .
2. Use the 90 % significant duration  to capture the excellent frequency time history within the corresponding time interval. Define the average frequency on the left and right beam sections at any given time as the step function, represented by and .
Fig. 5The excess pore water pressure ratio time history
3. Define the difference between the excellent frequency and the step function at any given time as the error function , as shown in Eq. (2). This error function effectively identifies the moments of rapid frequency change in the acceleration-time-frequency curve of the seismic record:
where, represents the average frequency at time , represents the excellent frequency at time , and and represent the start and end times, respectively, of the time interval within which the 90 % significant duration occurs.
4. Normalize the error function, and the time corresponding to 1 is considered as the estimated liquefaction time for the site.
3.1. Case analysis of Treasure Island
The proposed method was applied to process the ground seismic records monitored at the Treasure Island station during the 1989 Loma Prieta earthquake. The 90 % significant duration of the seismic record is shown in Fig. 6(a), from the PEER-NGA-West2 database . Additionally, the corresponding normalized error function time history is shown in Fig. 6(b).
Fig. 6a) Treasure Island acceleration time history , b) error function curve
The minimum value of the normalized error function corresponds to the estimated liquefaction time. From Fig. 6(b), the normalized error function occurs at approximately 13.6 s, indicating that liquefaction likely occurred at parts of the Treasure Island site at around 13.6 s.
3.2. Comparison of three liquefaction prediction time methods
As indicated by the published reference , it is established that the liquefaction triggering time for the Treasure Island site during the 1989 Loma Prieta earthquake event was 14.3 s. It was found that the liquefaction time for Treasure Island using the one-dimensional and two-dimensional non-linear site response analyses was 11.5 s and 11.9 s, respectively. In this paper, the proposed liquefaction time prediction method based on the prominent frequency time series resulted in a predicted liquefaction time of 13.6 s, which is very close to the true liquefaction time (with a difference of only 0.7 s). This indicates that the proposed method in this research demonstrates a high level of accuracy in predicting liquefaction time.
4. Validation of time-frequency site liquefaction prediction method
Based on the seven liquefaction sites with known liquefaction times from the 1995 Hyogoken-Nambu earthquake in Japan, the 1987 Superstition Hills earthquake in the United States, and the 2011 Tohoku-Oki earthquake in Japan, this paper further validates the reliability of the method for determining liquefaction times. To assess the correlation and reliability between the predicted liquefaction times and the actual liquefaction trigger times, the coefficient of determination (), as shown in Eq. (3), can be introduced. Based on a total of 7 cases in Table 2 and Treasure Island site, the calculated for the method is 0.99. This suggests that the proposed method exhibits a high level of reliability in predicting liquefaction time for different seismic events and sites:
where, represents the true liquefaction time, represents the predicted liquefaction time, and represents the mean of the predicted liquefaction times.
Table 2Comparison of predicted and real values of site liquefaction time
6.9, 1995 Hygoken-Nambu, Japan
6.6, 1987 Superstition Hills, USA
9.1, 2011 Tohoku-Oki, Japan
The present study proposes a new method for predicting liquefaction initiation time based on the time-frequency characteristics of ground acceleration. The defined step function and error function in this research effectively identifies the moments of rapid frequency changes in the ground seismic records, enabling the prediction of liquefaction initiation time. A 0.99 during the validation for 8 documented liquefaction sites associated with 4 seismic events indicates this method is reliable.
It should be noted that the proposed method only utilizes horizontal ground motion records. The inclusion of vertical ground motion records in the prediction of liquefaction initiation time may potentially improve the accuracy of the results.
J. M. Zhang., “Effect of large horizontal post-liquefaction deformation of level ground on pile foundation,” (in Chinese), Journal of Building Structures, No. 5, pp. 75–78, 2001, https://doi.org/10.14006/j.jzjgxb.2001.05.016
S. L. Kramer, S. S. Sideras, and M. W. Greenfield, “The timing of liquefaction and its utility in liquefaction hazard evaluation,” Soil Dynamics and Earthquake Engineering, Vol. 91, pp. 133–146, Dec. 2016, https://doi.org/10.1016/j.soildyn.2016.07.025
Y. Yang and E. Kavazanjian, “Numerical evaluation of liquefaction-induced lateral spreading with an advanced plasticity model for liquefiable sand,” Soil Dynamics and Earthquake Engineering, Vol. 149, p. 106871, Oct. 2021, https://doi.org/10.1016/j.soildyn.2021.106871
P. T. Özener, M. W. Greenfield, S. S. Sideras, and S. L. Kramer, “Identification of time of liquefaction triggering,” Soil Dynamics and Earthquake Engineering, Vol. 128, p. 105895, Jan. 2020, https://doi.org/10.1016/j.soildyn.2019.105895
W. Zhang, F. Ghahari, P. Arduino, and E. Taciroglu, “A deep learning approach for rapid detection of soil liquefaction using time-frequency images,” Soil Dynamics and Earthquake Engineering, Vol. 166, p. 107788, Mar. 2023, https://doi.org/10.1016/j.soildyn.2023.107788
R. W. Boulanger, L. H. Mejia, and I. M. Idriss, “Liquefaction at Moss Landing during Loma Prieta Earthquake,” Journal of Geotechnical and Geoenvironmental Engineering, Vol. 123, No. 5, pp. 453–467, May 1997, https://doi.org/10.1061/(asce)1090-0241(1997)123:5(453)
Y. Yang, J. Ma, S. Yuan, and W. Chen, “Numerical simulation of lateral spreading using an advanced plasticity model,” in International Conference on Geotechnical and Earthquake Engineering 2018, pp. 1–9, Jan. 2019, https://doi.org/10.1061/9780784482049.001
Y. M. Hashash et al., “Deepsoil 7.0 User Manual,” Urbana, IL, Board of Trustees of University of Illinois at Urbana-Champaign, 2020.
N. Puri, A. Jain, P. Mohanty, and S. Bhattacharya, “Earthquake response analysis of sites in state of Haryana using DEEPSOIL software,” Procedia Computer Science, Vol. 125, pp. 357–366, 2018, https://doi.org/10.1016/j.procs.2017.12.047
H. Chen, R. Sun, X. Yuan, and J. Zhang, “Variability of nonlinear dynamic shear modulus and damping ratio of soils,” in 14th World Conference on Earthquake Engineering, 2008.
H. Lu, Z. Lin, X. Zhan, Y. Yang, and D. Wu, “Numerical Simulation of seismic liquefaction for Treasure Island Site,” in Lecture Notes in Civil Engineering, pp. 407–415, 2023, https://doi.org/10.1007/978-981-99-1748-8_36
M. Hatanaka and A. Uchida, “Empirical correlation between penetration resistance and internal friction angle of sandy soils,” Soils and Foundations, Vol. 36, No. 4, pp. 1–9, Dec. 1996, https://doi.org/10.3208/sandf.36.4_1
B. van der Kwaak, “Modelling of dynamic pile behavior during an earthquake using PLAXIS 2D: Embedded beam (row),” Ms. thesis, Delft University of Technology, 2015.
Y. Hu, Y. Zhang, J. Liang, and R. R. Zhang, “Recording-based identification of site liquefaction,” Earthquake Engineering and Engineering Vibration, Vol. 4, No. 2, pp. 181–189, Dec. 2005, https://doi.org/10.1007/s11803-005-0001-3
M. D. Trifunac and A. G. Brady, “A study on the duration of strong earthquake ground motion,” Bulletin of the Seismological Society of America, Vol. 65, No. 3, pp. 581–626, Jun. 1975, https://doi.org/10.1785/bssa0650030581
T. D. Ancheta et al., “NGA-West2 database,” Earthquake Spectra, Vol. 30, No. 3, pp. 989–1005, Aug. 2014, https://doi.org/10.1193/070913eqs197m
About this article
This research was supported by Guangxi Natural Science Foundation under Grant (No. 2022GXNSFBA035569) and Innovation Project of Guangxi Graduate Education under Grant (No. YCSW2023297).
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.
The authors declare that they have no conflict of interest.