Liquefaction and Landslide#
Landslides and liquefaction are wellknown perils that accompany earthquakes. Basic models to describe their occurrence have been around for decades and are constantly improving. However, these models have rarely been incorporated into PSHA.
The tools presented here are implementations of some of the more common and appropriate secondary perils models. The intention is seamless incorporation of these models into PSH(R)A calculations done through the OpenQuake Engine, though the incorporation is a work in progress.
Tools for preparing the data for these models are also presented. This can be a nontrivial challenge, and consistent and correct data preparation is necessary for accurate secondary peril hazard and risk calculations.
Liquefaction models#
Several liquefaction models are implemented in the OpenQuake engine. One of them is the method developed for the HAZUS software by the US Federal Emergency Management Agency. This model involves categoriza tion of sites into liquefaction susceptibility classes based on geo technical characteristics, and a quanitative probability model for each susceptibility class. The remaining models are the academic geospatial models, i.e., statistical models that uses globally avai lable input variables as firstorder proxies to characterise satura tion and density properties of the soil. The shaking component is expressed either in terms of Peak Ground Acceleration (PGA) or Peak Ground Velocity (PGV).
HAZUS#
The HAZUS model classifies each site into a liquefaction susceptibility class (LSC) based on the geologic and geotechnical characteristics of the site, such as the sedimentological type and the deposition age of the unit. In addition to the LSC and the local ground acceleration at each site, the depth to groundwater at the site and the magnitude of the causative earthquake will affect the probability that a given site will experience liquefaction.
The equation that describes this probability is:
P(LPGA=a)
is the conditional probability that a site will fail based on the PGA and the LSC. P_{ml}
is the
fraction of the total mapped area that will experience liquefaction if P(LPGA=a)
reaches 1. These terms both have
LSCspecific coefficients; these are shown in Table 1.
K_m
is a magnitudecorrection factor that scales P(L)
for earthquake magnitudes other than M=7.5, potentially
to account for the duration of shaking (longer shaking increases liquefaction probability). K_w
is a groundwater
depth correction factor (shallower groundwater increases liquefaction probability).
LSC 
PGA min 
PGA slope 
PGA int 


very high 
0.09 
9.09 
0.82 
0.25 
high 
0.12 
7.67 
0.92 
0.2 
moderate 
0.15 
6.67 
1.0 
0.1 
low 
0.21 
5.57 
1.18 
0.05 
very low 
0.26 
4.16 
1.08 
0.02 
none 

0.0 
0.0 
0.0 
Table 1: Liquefaction values for different liquefaction susceptibility categories (LSC). PGA min is the minimum ground
acceleration required to initiate liquefaction. PGA slope is the slope of the liquefaction probability curve as a
function of PGA, and PGA int is the yintercept of that curve. P_{ml}
is the Map Area Proportion, which gives the area
of liquefaction within each map unit conditional on liquefaction occurring in the map unit.
Geospatial models#
Zhu et al. (2015)#
The model by Zhu et al. (2015) is a logistic regression model requiring specification of the Vs30 [m/s]
, the
Compound Topographic Index (CTI), a proxy for soil wetness or groundwater depth gwd [m]
, the PGA_{M,SM} [g]
experienced at a site, and the magnitude of the causative earthquake.
The model is quite simple. An explanatory variable X
is calculated as:
and the final probability is the logistic function:
The term PGA_{M,SM}
is the PGA corrected by magnitude scaling factor (MSF) that serves as proxy for earthquake
duration. The MSF
is calculated as per Idriss et al. (1991):
Both the CTI
and the Vs30
may be derived from digital elevation data. The Vs30
may be estimated from the
topographic slope through the equations of Wald and Allen (2007), which uses a very low resolution DEM compared to
modern offerings. As topographic slope tends to increase with increased DEM resolution, the estimated Vs30
does too;
therefore a lowresolution DEM (i.e., a 1 km resolution) must be used to calculate Vs30
, rather than the 30 m DEM
that is the current standard. This results in a more accurate Vs30
for a given slope measurement, but it also means
that in an urban setting, subkmscale variations in slope are not accounted for.
The CTI (Moore et al., 1991) is a proxy for soil wetness that relates the topographic slope of a point to the upstream drainage area of that point, through the relation:
where d_a
is the upstream drainage area per unit width through the flow direction (i.e. relating to the DEM
resolution). It was developed for hillslopes, and is not meaningful in certain very flat areas such as the valley
floors of major lowgradient rivers, where the upstream drainage areas are very large. Unfortunately, this is exactly
where liquefaction is most expected away from coastal settings.
Model’s prediction can be transformed into binary class (liquefaction occurrence or nonoccurrence) via probability threshold value. The authors proposed a threshold of 0.2.
Bozzoni et al. (2021)#
The parametric model developed by Bozzoni et al. (2021), keeps the same input variables (i.e., PGA_{M,SM}
, CTI
,
Vs30
) and functional form as in Zhu et al. (2015). Regression parameters are calibrated based on the liquefaction
case histories observed during seismic events in Europe. The implemented model is associated with the ADASYN sampling
algorithm. The explanatory variable \(X\) is computed as:
and the probability of liquefaction in calculated using equation (3).
The adopted probability threshold of 0.57 converts the probability of liquefaction into binary outcome.
Zhu et al. (2017)#
Two parametric models are proposed by Zhu and others (2017), a coastal model (Model 1), and a more general model
(Model 2). A coastal event is defined as one where the liquefaction occurrences are, on average, within 20 km of the
coast; or, for earthquakes with insignificant or no liquefaction, epicentral distances less than 50 km.The implemented
geospatial models are for global use. An extended set of input parameters is used to describe soil properties (its
density and wetness). The ground shaking is characterised by PGV [cm/s]
. Soil density is described by Vs30 [m/s]
.
Soil wetness in Model 1 is chatacterised by a set of features: mean annual precipitation precip [mm]
, distance to
the coast d_{c} [km]
, and distance to the river d_{r} [km]
. Distance to the coast also indicates the geologic
age  younger deposits are found near the coast. Soil wetness in Model 2 is characterised by closest distance to the
water body d_{w} [km]
which is determined as \min(d_{c}, d_{r})
, and the water table depth wtd [m]
. Mean
annual precipitation is from a global layer developed by Hijmans et al. (2005). Distance to the nearest river is
calculated based on the HydroSHEDS database (Lehner et al. 2008). Water table depth is retreived from a global dataset
by Fan et al (2013).Distance to the nearest coastline data was computed from https://oceancolor.gsfc.nasa.gov.
The explanatory varibale \(X\) is calculated as:
Model 1:
Model 2:
and the probability of liquefaction is calculated using equation (3). Zero probability is heuristically assigned if PGV < 3 cm/s
.
The proposed probability threshold to convert to class outcome is 0.4.
Another model’s outcome is liquefaction spatial extent (LSE). After an earthquake LSE is the spatial area covered by surface manifestations of liquefaction reported as a percentage of liquefied material within that pixel. Logistic regression with the same form was fit for the two models, with only difference in squaring the denominator to improve the fit. The regression coefficients are given in Table 2.:
Parameters 
Model 1 
Model 2 

a 
42.08 
49.15 
b 
62.59 
42.40 
c 
11.43 
9.165 
Table 2: Parameters for relating probabilities to areal liquefaction percent.
Rashidian et al. (2020)#
The model proposed by Rashidian et al. (2020) keeps the same functional form as the general model (Model 2) proposed by
Zhu et al. (2017); however, introdu cing two constraints to address the overestimation of liquefaction extent. The mean
annual precipitation has been capped to 1700 mm. No liquefaction is heuri stically assign when pga < 0.1 g
as an
additional measure to decrease the overestimation of liquefaction. Additional novelty introduced in this model is the
magnitude scaling factor MSF
to multiply the PGV to mitigate the potential overprediction in earthquake with low
magnitude.
The explanatory variable X
is evaluated using the equation (7) that corresponds to the general model of Zhu et al.
(2017). The spatial extent is evaluated identically using the equation (8).
The proposed probability threshold to convert to class outcome is 0.4.
Akhlagi et al. (2021)#
Expanding the liquefaction inventory to include 51 earthquake, Akhlagi et al. (2021) proposed two candidate models to
predict probability of liquefaction. Shaking is expressed in terms of PGV [cm/s]. Soil saturation is characterised using
the set of proxies: distance to the nearest coastline d_{c} [m]
, distance to the closest river d_{r} [m]
,
elevation from the closest water body Z_{wb} [m]
. Soil density is characterised either by Vs30 [cm/s]
or
topographic roughness index TRI
which is defined as the mean difference between a central pixel and its eight
surrounding cells. The explanatory variables of two candidate models are:
Model 1:
Model 2:
and the probability of liquefaction is calculated using equation (3). Zero probability is heuristically assigned if
PGV < 3 cm
or Vs30 > 620 m/s
.
The proposed probability threshold to convert to class outcome is 0.4.
Allstadth et al. (2022)#
The model proposed by Allstadth et al. (2022) uses the model proposed by Rashidian et al. (2020) as a base with slight
changes to limit unrealistic extrapolations. The authors proposed capping the mean annual precipitation at 2500 mm, and
PGV at 150 cm/s. The magnitude scaling factor MSF
, explanatory variables X
, probability of liquefaction P(L)
,
and liquefaction spatial extent LSE
are calculated using the set of equations previously shown. The proposed
probability threshold to convert to class outcome is 0.4.
Todorovic et al. (2022)#
A nonparametric model was proposed to predict liquefaction occurrence and the associated probabilities. The general
model was trained on the dataset including inventories from over 40 events. A set of candidate variables were considered
and the ones that correlate the best with liquefaction occurrence are identified as: strain proxy, a ratio between
pgv [cm/s]
and Vs30 [m/s]
; distance to the closest water body d_{w} [km]
, water table depth gwd [m]
,
average precipitation precip [mm]
.
Permanent ground displacements due to liquefaction#
Evaluation of the liquefaction induced permanent ground deformation is conducted using the methodology developed for the HAZUS software by the US Federal Emergency Management Agency. Lateral spreading and vertical settlements can have detrimental effects on the built environement.
Lateral spreading (Hazus)#
The expected permanent displacement due to lateral spreading given the susceptibility category can be determined as:
Where: E[PGD(PGA/PL_{SC})=a]
is the expected ground displacement given the susceptibility category under a
specified level of normalised shaking, and is calculated as:
(PGA/PGA(t))
PGA(t)
is theminimum shaking level to induce liquefaction (see Table 1) K_{\Delta}
is the
displacement correction factor given thhat modify the displacement term for magnitudes other than M7.5
:
Vertical settlements (Hazus)#
Ground settlements are assumed to be related to the area’s susceptibility category. The ground settlement amplitudes are given in Table 3 for the portion of a soil deposit estimated to experience liquefaction at a given ground motion level. The expected settlements at the site is the product of the probability of liquefaction (equation 1) and the characteristic settlement amplitude corresponding to the liquefaction susceptibility category (LSC).
LSC 
Settlements (inches) 

very high 
12 
high 
6 
moderate 
2 
low 
1 
very low 
0 
none 
0 
Table 3: Ground settlements amplitudes for liquefaction susceptibility categories.
Landslide models#
Landslides are considered as one of the most damaging secondary perils associated with earthquakes. Earthquakeinduced
landslides occurs when the static and inertia forces within the sliding mass reduces the factor of safety below 1.
Factors contributing to a slope failure are rather complex. The permanentdisplacement analysis developed by Newmark
(1965) is used to model the dynamic performance of slopes (Jibson 2020, 2007). It considers a slope as a rigid block
resting on an inclined plane at an angle \alpha
(derived from Digital Elevation Model, DEM). When the input motion which
is expressed in terms of acceleration exceeds the critical acceleration a_{c}
, the block starts to move. The crtical
acceleration accounts for the shear strength and geometrical characteristics of the sliding surface, and is calculated
as:
The lower bound of a_{c}
is set to 0.05 to avoid unrealistically large displacements. The static factor of safety is
calculated as:
where: c' [Pa]
is the effective cohession with typical values ranging from 20 kPa for soils up to 20 MPa for
unfaulted rocks. \alpha [\degrees]
is the slope angle. \phi' [\degrees]
is the effective friction angle with
typical values ranging from 30 to 40 degrees. \gamma [kg/m^3]
is the dry_density of the soil or rock. It ranges
from ~1500 kg/m^3
for soils to ~ 2500  3200 kg/m^3
. t [m]
is the slopenormal thickness of a failure slab
with the default value of 2.5 meters. m
is the proportion of slab thickness that is saturated with default value of
0.1. \gamma_{w} [kg/m^3]
is the unit weight of water which equals to 1000 kg/m^3
.
Note that the units of the input parameters reported in this document corresponds to the format required by the Engine to produce correct results. The first and second term of the the equation corresponds to the cohesive and frictional components of the strength, while the third component accounts for the strength reduction due to pore pressure.
A variety of regression equations can be used to estimate the Newmark displacements, and within the engine, Newmark displacement as a function of critical acceleration ratio and moment magnitude is implemented. The displacement is in units of meters.
The computed displacements do not necessarily correspond directly to measurable slope movements in the field, but the modeled displacements provide an index to correlate with field performance. Jibson (2000) compared the predicted displacements with observations from 1994 Northridge earthquake and fit the data with Weilbull curve. The following equation can be used to estimate the probability of slope failure as a function of Newmark displacement.
The rockslope failures are the other common effect observed in earthquakes. The methodology proposed by Grant et al., (2016) captures the brittle behavior associated with rockslope failures and discontinuities common in rock masses. The static factor of safety is computed as:
where: c [Pa]
is the cohession with typical values ranging from 20 kPa for soils up to 20 MPa for unfaulted rocks.
c_{r}
is the cohesion provided by the root systems of vegetated hillslopes. Here, we adopted the default value of 0
root cohesion. \alpha [\degrees]
is the slope angle. \gamma [kg/m^3]
is the dry_density of the soil or rock. It
ranges from ~1500 kg/m^3
for soils to ~ 2500  3200 kg/m^3
. h [m]
is the vertical height of the failure mass
and it corresponds to 1/4 of the local relief H calculated based on the moving window analysis. \phi' [\degrees]
is
the effective friction angle with typical values ranging from 30 to 40 degrees. \beta
is the slope’s critical angle
calculated as:
The critical acceleration is computed similarly to equation (14). For rock slope failures, the \(\alpha\) term is replaced with \(\beta\).
Finally, the coseismic displacements are estimated using Jibson’s (2007) sliding block displacement regression equation:
Reference#
[1] HAZUSMH MR5 Earthquake Model Technical Manual (https://www.hsdl.org/?view&did=12760)
[2] Youd, T. L., & Idriss, I. M. (2001). Liquefaction Resistance of Soils: Summary Report from the 1996 NCEER and 1998 NCEER/NSF Workshops on Evaluation of Liquefaction Resistance of Soils. Journal of Geotechnical and Geoenvironmental Engineering, 127(4), 297–313. https://doi.org/10.1061/(asce)10900241(2001)127:4(297)
[3] I. D. Moore, R. B. Grayson & A. R. Ladson (1991). Digital terrain modelling: A review of hydrological, geomorphological, and biological applications. Journal of Hydrological Processes, 5(1), 330. https://doi.org/10.1002/hyp.3360050103
[4] Wald, D.J., Allen, T.I., (2007). Topographic Slope as a Proxy for Seismic Site Conditions and Amplification. Bull. Seism. Soc. Am. 97 (5), 1379–1395.
[5] Zhu et al., 2015, ‘A Geospatial Liquefaction Model for Rapid Response and Loss Estimation’, Earthquake Spectra, 31(3), 18131837.
[6] Bozzoni, F., Bonì, R., Conca, D., Lai, C. G., Zuccolo, E., & Meisina, C. (2021). Megazonation of earthquakeinduced soil liquefaction hazard in continental Europe. Bulletin of Earthquake Engineering, 19(10), 4059–4082. https://doi.org/10.1007/s10518020010086
[7] Zhu, J., Baise, L. G., & Thompson, E. M. (2017). An updated geospatial liquefaction model for global application. Bulletin of the Seismological Society of America, 107(3), 1365–1385. https://doi.org/10.1785/0120160198
[8] Rashidian, V., & Baise, L. G. (2020). Regional efficacy of a global geospatial liquefaction model. Engineering Geology, 272, 105644. https://doi.org/10.1016/j.enggeo.2020.105644
[9] Allstadt, K. E., Thompson, E. M., Jibson, R. W., Wald, D. J., Hearne, M., Hunter, E. J., Fee, J., Schovanec, H., Slosky, D., & Haynie, K. L. (2022). The US Geological Survey ground failure product: Nearrealtime estimates of earthquaketriggered landslides and liquefaction. Earthquake Spectra, 38(1), 5–36. https://doi.org/10.1177/87552930211032685
[10] Baise, L. G., Akhlaghi, A., Chansky, A., Meyer, M., & Moeveni, B. (2021). USGS Award #G20AP00029. Updating the Geospatial Liquefaction Database and Model. Tufts University. Medford, Massachusetts, United States.
[11] Todorovic, L., Silva, V. (2022). A liquefaction occurrence model for regional analysis. Soil Dynamics and Earthquake Engineering, 161, 1–12. https://doi.org/10.1016/j.soildyn.2022.107430
[12] Newmark, N.M., 1965. Effects of earthquakes on dams and embankments. Geotechnique 15, 139–159.
[13] Jibson, R.W., Harp, E.L., & Michael, J.A. (2000). A method for producing digital probabilistic seismic landslide hazard maps. Engineering Geology, 58(34), 271289. https://doi.org/10.1016/S00137952(00)000399
[14] Jibson, R.W. (2007). Regression models for estimating coseismic landslide displacement. Engineering Geology, 91(24), 209218. https://doi.org/10.1016/j.enggeo.2007.01.013
[15] Grant, A., Wartman, J., & Grace, A.J. (2016). Multimodal method for coseismic landslide hazard assessment. Engineering Geology, 212, 146160. https://doi.org/10.1016/j.enggeo.2016.08.005