smt_model

The SMT modeling team developed probabilistic susceptibility, triggering, and manifestation models that are implemented using the functions below. The SMT model is described by Ulmer et al. (2024)

function	description
smt_model(depth, qt, fs, amax, m, pa=101.325, **kwargs)	Computes probabilities of susceptibility, triggering, and manifestation
smt_model_fragility(depth, qt, fs, amax, m, pa=101.325, **kwargs)	Computes a set of fragility functions representing probability of profile manifestation as a function of amax and m
get_pfs(Ic)	Computes probability factor for liquefaction susceptibility
get_pfts(csrm_hat, crr_hat)	Computes probability factor for liquefaction triggering conditional on susceptibility
get_pfmt(ztop, Ic)	Computes probability factor for manifestation of a layer conditional on triggering of the layer
get_pmp(pfmt, pfts, pfs, Ksat, t)	Computes probability of manifestation at surface of profile

The functions developed by the SMT use the following helper functions.

function	description
get_FC_from_Ic(Ic, epsilon=0.0)	Computes fines content from soil behavior type index following the method by Hudson et al. (2024)
get_Ic_Qtn_Fr(qt, fs, sigmav, sigmavp, pa=101.325, maxiter=30)	computes soil behavior type index, Ic, overburden-corrected dimensionless cone tip resistance, Qtn, and dimensionless friction ratio, Fr
get_qc1N_qc1Ncs(qt, fs, sigmav, sigmavp, FC, pa=101.325, maxiter=30)	computes overburden-corrected dimensionless cone tip resistance qc1N that uses qc1Ncs to set exponent for overburden correction, and overburden- and fines-corrected cone tip resistance qc1Ncs
box_cox(x, lam)	Computes Box-Cox transform. $\hat{x} = \frac{x^\lambda - 1.0}{\lambda}$
inv_box_cox(x_hat, lam)	Computes inverse Box-Cox transform $x = \left(\hat{x} \cdot \lambda + 1.0\right)^{\frac{1.0}{\lambda}}$
get_crr_hat(qc1Ncs)	Computes Box-Cox transformed cyclic resistance ratio
get_csrm(amax, m, sigmav, sigmavp, depth, qc1Ncs, pa=101.325)	Computes cyclic stress ratio corrected for magnitude and overburden

Note

For Box-Cox transforms, $\lambda_{CSR} = -0.361$ for CSR and CRR, and $\lambda_{DR} = 1.226$ for relative density.

smt_model(depth, qt, fs, amax, m, pa=101.325, **kwargs)

This function computes the probabilities of susceptibility, triggering, and manifestation and returns a tuple of relevant quantities. It is a wrapper function that parses user inputs, then calls other functions contained in the SMT software package. Equations are from Ulmer et al. (2024).

Arguments

argument	format	description
depth	Numpy array [dtype = float]	Depth values.
qt	Numpy array [dtype = float]	Cone tip resistance values in same units as pa.
fs	Numpy array [dtype = float]	Cone sleeve friction values in same units as pa.
amax	float	Maximum horizontal surface acceleration in g.
m	float	Earthquake magnitude.
pa	float	Atmospheric pressure. Optional. Default = 101.325 kPa.

Keyword Arguments (kwargs)

kwarg	format	description
gamma	float	Soil total unit weight
dGWT	float	Depth to groundwater table.
gammaw	float	Unit weight of water. Optional. Default = 9.81 kN/m3
sigmav	Numpy array [dtype=float]	Vertical total stress in same units as pa.
sigmavp	Numpy array [dtype=float]	Vertical effective stress in same units as pa.
Ksat	Numpy array [dtype=float]	Saturation factor to multiply probability of triggering.

Returns

The function returns a tuple containing the following variables.

variable	format	description
Ic	Numpy array [dtype=float]	Soil behavior type index
qt_inv	Numpy array [dtype=float]	Inverse-filtered cone tip resistance.
fs_inv	Numpy array [dtype=float]	Inverse-filtered sleeve friction.
qc1Ncs_inv	Numpy array [dtype=float]	Inverse-filtered overburden- and fines-corrected cone tip resistance.
Ic_inv	Numpy array [dtype=float]	Inverse-filtered soil behavior type index.
FC	Numpy array [dtype=float]	Fines content computed from Ic_inv.
qc1Ncs	Numpy array [dtype=float]	Overburden- and fines-corrected cone tip resistance.
ztop	Numpy array [dtype=float]	Depth to top of layers.
zbot	Numpy array [dtype=float]	Depth to bottom of layers.
qc1Ncs_lay	Numpy array [dtype=float]	Overburden- and fines-corrected cone tip resistance for layers.
Ic_lay	Numpy array [dtype=float]	Soil behavior type index for layers.
Ksat_lay	Numpy array [dtype=float]	Ksat value for layers.
pfs	Numpy array [dtype=float]	Probability factor for susceptibility of layers.
pfts	Numpy array [dtype=float]	Probability factor for triggering conditioned on susceptibility.
pft	Numpy array [dtype=float]	Probability factor for triggering of layers.
pfmt	Numpy array [dtype=float]	Probability factor for manifestation conditioned on triggering.
pfm	Numpy array [dtype=float]	Probability factor for manifestation of layers.
pmp	Float	Probability of profile manifestation.

Note

1. Either dGWT and gamma must be specified, or sigmav and sigmavp must be specified.
2. If sigmav, sigmavp, and gamma are specified, gamma is ignored.
3. If sigmav and sigmavp are specified and dGWT is not specified, dGWT is computed as the deepest point where sigmav = sigmavp.
4. After the layering algorithm is completed, a layer interface is imposed at dGWT if the following conditions are met: (i) dGWT does not already coincide with a layer interface, (ii) dGWT is below the ground surface, and (iii) dGWT is above the bottom of the profile.
5. If sigmav, sigmavp, and dGWT are specified, no checks are performed to make sure they are consistent.
6. If Ksat is not specified, it is be computed as 0.0 above dGWT and 1.0 below dGWT.
7. If Ksat and dGWT are both specified, no checks are performed to make sure they are consistent.
8. All input Numpy arrays must have the same length.
9. The following output Numpy arrays have the same length as the input Numpy arrays: qt_inv, fs_inv, Ic_inv, FC, qc1Ncs.
10. The following output Numpy arrays have a length equal to the number of layers identified by the layering algorithm: ztop, zbot, qc1Ncs_lay, Ic_lay, Ksat_lay, pfs, pfts, pft, pfmt, pfm.

smt_model_fragility(depth, qt, fs, amax, m, pa=101.325, **kwargs)

This function computes liquefaction fragility defining P[M_P] as a function of a_max and M. It is a wrapper function that parses user inputs, then calls other functions contained in the SMT software package, and returns a 2D array of P[M_P] values. Equations are from Ulmer et al. (2024).

Arguments

argument	format	description
depth	Numpy array [dtype = float]	Depth values.
qt	Numpy array [dtype = float]	Cone tip resistance values in same units as pa.
fs	Numpy array [dtype = float]	Cone sleeve friction values in same units as pa.
amax	Numpy array [dtype = float]	Maximum horizontal surface acceleration in g.
m	Numpy array [dtype = float]	Earthquake magnitude.
pa	float	Atmospheric pressure. Optional. Default = 101.325 kPa.

Keyword Arguments (kwargs)

kwarg	format	description
gamma	float	Soil total unit weight
dGWT	float	Depth to groundwater table.
gammaw	float	Unit weight of water. Optional. Default = 9.81 kN/m3
sigmav	Numpy array [dtype=float]	Vertical total stress in same units as pa.
sigmavp	Numpy array [dtype=float]	Vertical effective stress in same units as pa.
Ksat	Numpy array [dtype=float]	Saturation factor to multiply probability of triggering.

Returns

variable	format	description
pmp	Numpy ndarray [dtype=float]	Probability of profile manifestation corresponding to arrays of amax and m.

Note

Dimensions of pmp are (len(m), len(amax)).

get_pfs(Ic)

This function computes the probability of liquefaction susceptibility conditioned on soil behavior type index. Equations are from Ulmer et al. (2024).

Parameters

argument	format	description
Ic	float or Numpy array [dtype = float]	Soil behavior type index

Returns

variable	format	description
pfs	float or Numpy array [dtype = float]	Probability factor for susceptibility.

Note

$PF_S = \frac{1.0}{1.0 + exp\left[-14.8\left(\frac{I_c}{2.635} - 1.0 \right) \right]}$

get_pfts(csrm_hat, crr_hat)

This function returns the probability factor of liquefaction triggering conditional on susceptibility. Equations are from Ulmer et al. (2024).

Parameters

argument	format	description
csrm_hat	float or Numpy array [dtype=float]	Box-Cox transformed overburden- and magnitude-corrected cyclic stress ratio
crr_hat	float or Numpy array [dtype=float]	Box-Cox transformed cyclic resistance ratio

Returns

variable	format	description
pfts	float or Numpy array [dtype=float]	Probability factor for liquefaction triggering conditional on susceptibility

Note

$PF_{T|S} = \frac{1.0}{1.0 + exp\left[-2.020\left(\widehat{CSR_M} - \widehat{CRR} \right) \right]}$

get_pfmt(ztop, Ic)

This function returns the probability factor for manifestation conditional on liquefaction triggering. Equations are from Ulmer et al. (2024).

argument	format	description
ztop	float or Numpy array [dtype=float]	Depth to the top of the layer of interest in meters
Ic	float or Numpy array [dtype=float]	Median soil behavior type index of the layer of interest

Note

$PF_{M|T} = \frac{1.0}{1.0 + exp\left[-\left(7.613 - 0.338 z_{top} - 3.042 I_c \right) \right]}$

Returns

variable	format	description
pfmt	float or Numpy array [dtype=float]	Probability factor for manifestation of a layer conditional on liquefaction triggering of the layer

get_pmp(pfmt, pfts, pfs, Ksat, t)

This function returns the probability of manifestation of a profile. Equations are from Ulmer et al. (2024).

Parameters

argument	format	description
pfmt	Numpy array [dtype=float]	Probability factor of manifestation conditional on triggering
pfts	Numpy array [dtype=float]	Probability factor of triggering conditional on susceptibility
pfs	Numpy array [dtype=float]	Probability factor of susceptibility
Ksat	Numpy array [dtype=float]	Saturation factor. Ksat = 0.0 above groundwater table and 1.0 below groundwater table
t	Numpy array [dtype=float]	Layer thickness in meters

Returns

variable	format	description
pmp	float	Probability of manifestation of the profile in percent

Note

$P[M_P] = 1.0 - \prod_{i=1}^{N_{layer}} \left(1.0 - PF_{M|T} PF_{T|S} PF_S K_{sat} \right) ^ {t / t_c}$
$t_c = 2.0 m$

get_FC_from_Ic(Ic, epsilon=0.0)

This function computes fines content from soil behavior type index following the method by Hudson et al. (2024).

Parameters

argument	format	description
Ic	float, or Numpy array [dtype = float]	Soil behavior type index
<epsilon>	float, or Numpy array [dtype = float]	Number of standard deviations. Optional. Default = 0.0

Returns

variable	format	description
FC	float, or Numpy array [dtype = float]	Fines content in percent

Note

FC will be a scalar valued float if Ic and epsilon are scalar valued floats, and will be a Numpy array if either Ic and / or epsilon are Numpy arrays. If Ic and epsilon are both Numpy arrays, they must have the same length or an error will be returned.

get_Ic_Qtn_Fr(qt, fs, sigmav, sigmavp, pa=101.325, maxiter=30)

This function computes Ic, Qtn, Fr, qc1N, and qc1Ncs. Iterations are required to solve for these parameters because relationships among Ic, Qtn, Fr, qc1N, and qc1Ncs are implicit. Equations are from Robertson (2016).

Parameters

argument	format	description
qt	Numpy array [dtype = float]	Cone tip resistance
fs	Numpy array [dtype = float]	Cone sleeve friction
sigmav	Numpy array [dtype = float]	Vertical total stress
sigmavp	Numpy array [dtype = float]	Vertical effective stress
<pa>	float	Atmospheric pressure in same units as qt, fs, sigmav, and sigmavp. Optional. Default = 101.325 kPa
<maxiter>	int	Maximum number of iterations

Returns

variable	format	description
Ic	Numpy array [dtype = float]	Soil behavior type index.
Qtn	Numpy array [dtype = float]	Dimensionless cone tip resistance.
Fr	Numpy array [dtype = float]	Dimensionless cone sleeve friction.

Note

Values of qt, fs, sigmav, and sigmavp that are less than or equal to zero will cause errors due to the logarithms that must be computed. Such values are set to 0.001*pa before computing Ic, Qtn, and Fr. A warning is not issued when this occurs.

get_qc1N_qc1Ncs(qt, fs, sigmav, sigmavp, FC, pa=101.325, maxiter=30)

This function computes qc1N, and qc1Ncs. Iterations are required to solve for these parameters. The equations are implicit because the stress normalization exponent depends on the fines correction. Equations used herein are from Boulanger and Idriss (2016).

Parameters

argument	format	description
qt	Numpy array [dtype = float]	Cone tip resistance
fs	Numpy array [dtype = float]	Cone sleeve friction
sigmav	Numpy array [dtype = float]	Vertical total stress
sigmavp	Numpy array [dtype = float]	Vertical effective stress
FC	Numpy array [dtype = float]	Fines content in percent
<pa>	float	Atmospheric pressure in same units as qt, fs, sigmav, and sigmavp. Optional. Default = 101.325 kPa
<maxiter>	int	Maximum number of iterations

Returns

variable	format	description
qc1N	Numpy array [dtype = float]	Overburden corrected cone tip resistance, where exponent for overburden correction is a function of qc1Ncs.
qc1Ncs	Numpy array [dtype = float]	Overburden- and fines-corrected cone tip resistance

Note

Values of qt, fs, sigmav, and sigmavp that are less than or equal to zero will cause errors due to the logarithms that must be computed. Such values are set to 0.001*pa before computing qc1N, and qc1Ncs. A warning is not issued when this occurs.

box_cox(x, lam)

This function computes the Box-Cox transformation of a variable following Box and Cox (1964).

Parameters

argument	format	description
x	float or Numpy array [dtype = float]	Variable to be transformed
labmda	float	Exponent term used in transform

Returns

variable	format	description
x_hat	float or Numpy array [dtype = float]	Transformed variable

Note

x_hat will have the same type as x.

inv_box_cox(x_hat, lam)

This function computes the inverse Box-Cox transformation of a variable following Box and Cox (1964).

Parameters

argument	format	description
x_hat	float or Numpy array [dtype = float]	Variable to be inverse transformed
labmda	float	Exponent term used in transform

Returns

variable	format	description
x	float or Numpy array [dtype = float]	Inverse transformed variable

Note

x will have the same type as x_hat.

get_crr_hat(qc1Ncs)

This function computes the Box-Cox transformed value of cyclic resistance ratio based on overburden- and fines-corrected cone tip resistance. qc1Ncs is converted to relative density using equations from Idriss and Boulanger (2008), and cyclic resistance is computed using equations from Ulmer et al. (2024).

Parameters

argument	format	description
qc1Ncs	float or Numpy array [dtype = float]	Overburden- and fines-corrected cone tip resistance

Returns

variable	format	description
crr_hat	float or Numpy array [dtype = float]	Box-Cox transformed cyclic resistance ratio

Note

crr_hat will have the same type as qc1Ncs.

get_csrm(amax, m, sigmav, sigmavp, depth, qc1Ncs, pa=101.325)

This function returns magnitude- and overburden corrected cyclic stress ratio. Depth-reduction factor, $r_d$ is from Lasley et al. (2016), magnitude scaling factor, $MSF$, is from Lasley et al. (2017), and overburden correction factor, $K_\sigma$, is from Ulmer et al. (2024).

Parameters

argument	format	description
amax	float	Peak horizontal surface acceleration in g
m	float	Earthquake magnitude
sigmav	Numpy array [dtype = float]	Vertical total stress
sigmavp	Numpy array [dtype = float]	Vertical effective stress
depth	Numpy array [dtype = float]	Depth values
qc1Ncs	Numpy array [dtype = float]	Overburden- and fines-corrected cone tip resistance
<pa>	float	Atmospheric pressure in same units as sigmav and sigmavp. Optional. Default = 101.325 kPa

Note

*sigmav, sigmavp, depth and qc1Ncs must be Numpy arrays of the same length. An error will result if they are not Numpy arrays or if they are Numpy arrays of different length.

Returns

variable	format	description
csrm	Numpy array [dtype = float]	Peak horizontal surface acceleration in g

Example: smt_model()

To download the files below, right click the link and use "Save Link As".

Data file example_cpt.csv

Jupyter notebook smt_model_example.ipynb

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import pandas as pd
import ngl_tools.smt as smt

def interleave(v1, v2):
    vout = np.empty(2 * len(v1), dtype=float)
    vout[0::2] = v1
    vout[1::2] = v2
    return vout 
    
# Ready example CPT file and extract data into Numpy arrays
df = pd.read_csv('example_cpt.csv')
depth = df['depth'].values
qt = df['qt'].values
fs = df['fs'].values

# Define earthquake loading conditions
amax = 0.3
m = 7.0
gamma = 19
dGWT = 2.0

Ic, qt_inv, fs_inv, qc1Ncs_inv, Ic_inv, FC, qc1Ncs, \
ztop, zbot, qc1Ncs_lay, Ic_lay, Ksat_lay, \
pfs, pfts, pft, pfmt, pfm, pmp = smt.smt_model(depth, qt, fs, amax, m, dGWT=dGWT, gamma=gamma)

fig = plt.figure(figsize=(10,10))

gs = GridSpec(2, 12, figure=fig, wspace=0.8)
ax = []
ax.append(fig.add_subplot(gs[0, 0:3]))
ax.append(fig.add_subplot(gs[0, 3:6], sharey=ax[0]))
ax.append(fig.add_subplot(gs[0, 6:9], sharey=ax[0]))
ax.append(fig.add_subplot(gs[0, 9:12], sharey=ax[0]))

ax.append(fig.add_subplot(gs[1, 0:2], sharey=ax[0]))
ax.append(fig.add_subplot(gs[1, 2:4], sharey=ax[0]))
ax.append(fig.add_subplot(gs[1, 4:6], sharey=ax[0]))
ax.append(fig.add_subplot(gs[1, 6:8], sharey=ax[0]))
ax.append(fig.add_subplot(gs[1, 8:10], sharey=ax[0]))
ax.append(fig.add_subplot(gs[1, 10:12], sharey=ax[0]))

for i, a in enumerate(ax):
    if((i!=0) and (i!=4)):
        plt.setp(ax[i].get_yticklabels(), visible=False)

ax[0].plot(qt/1000, depth)
ax[0].plot(qt_inv/1000, depth)
ax[1].plot(fs/1000, depth)
ax[1].plot(fs_inv/1000, depth)
ax[2].plot(Ic, depth, label='measured')
ax[2].plot(Ic_inv, depth, label='inverse filtered')
ax[3].plot(qc1Ncs, depth)
ax[3].plot(qc1Ncs_inv, depth)
ax[4].plot(interleave(pfs, pfs), interleave(ztop,zbot))
ax[4].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(pfs, pfs), facecolor='C0', alpha=0.2)
ax[5].plot(interleave(pfts, pfts), interleave(ztop,zbot))
ax[5].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(pfts, pfts), facecolor='C0', alpha=0.2)
ax[6].plot(interleave(pft, pft), interleave(ztop,zbot))
ax[6].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(pft, pft), facecolor='C0', alpha=0.2)
ax[7].plot(interleave(pfmt, pfmt), interleave(ztop,zbot))
ax[7].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(pfmt, pfmt), facecolor='C0', alpha=0.2)
ax[8].plot(interleave(Ksat_lay, Ksat_lay), interleave(ztop,zbot))
ax[8].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(Ksat_lay, Ksat_lay), facecolor='C0', alpha=0.2)
ax[9].plot(interleave(pfm, pfm), interleave(ztop,zbot))
ax[9].fill_betweenx(interleave(ztop, zbot), np.zeros(2*len(ztop)), interleave(pfm, pfm), facecolor='C0', alpha=0.2)
ax[2].axvspan(1,1.31,alpha=0.3,color='darkorange')
ax[2].axvspan(1.31,2.05,alpha=0.3,color='peru')
ax[2].axvspan(2.05,2.6,alpha=0.3,color='palegreen')
ax[2].axvspan(2.6,2.95,alpha=0.3,color='mediumaquamarine')
ax[2].axvspan(2.95,3.6,alpha=0.3,color='navy')
ax[2].axvspan(3.6,4.0,alpha=0.3,color='saddlebrown')
ax[0].set_ylabel('depth [m]')
ax[4].set_ylabel('depth [m]')
ax[0].set_xlabel(r'$q_c$ [MPa]')
ax[1].set_xlabel(r'$f_s$ [MPa]')
ax[2].set_xlabel(r'$I_c$ [-]')
ax[3].set_xlabel(r'$q_{c1Ncs}$ [-]')
ax[4].set_xlabel(r'$PF_S$')
ax[5].set_xlabel(r'$PF_{T|S}$')
ax[6].set_xlabel(r'$PF_T$')
ax[7].set_xlabel(r'$PF_{M|T}$')
ax[8].set_xlabel(r'$K_{sat}$')
ax[9].set_xlabel(r'$PF_M$')
for a in ax:
    a.grid(True, alpha=0.5)
ax[2].plot(interleave(Ic_lay,Ic_lay), interleave(ztop,zbot), color='black', zorder=2, linewidth=0.5, label='layered')
ax[3].plot(interleave(qc1Ncs_lay,qc1Ncs_lay), interleave(ztop,zbot), color='black', zorder=2, linewidth=0.5)
ax[0].invert_yaxis()
ax[2].set_xlim(1.0,4.0)
ax[5].set_xlim(0.0,1.0)
ax[6].set_xlim(0.0,1.0)
ax[7].set_xlim(0.0,1.0)
ax[8].set_xlim(0.0,1.0)
ax[9].set_xlim(0.0,1.0)
ax[2].legend(loc='upper center', bbox_to_anchor=(0.9,1.15), ncol=3, fancybox=True, shadow=True)
#fig.suptitle(r'$P[M_P]$ = ' + str(np.round(pmp*100, 1)) + ' %')
ax[9].annotate(r'$P[M_P]$=' + str(np.round(pmp*100, 1)) + '%', xy=(0.03,0.03), xycoords='axes fraction')

Example

To download the files below, right click the link and use "Save Link As".

Data file example_cpt.csv

Jupyter notebook smt_model_fragility_example.ipynb

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import pandas as pd
import ngl_tools.smt as smt
    
# Ready example CPT file and extract data into Numpy arrays
df = pd.read_csv('example_cpt.csv')
depth = df['depth'].values
qt = df['qt'].values
fs = df['fs'].values

# Define earthquake loading conditions
amax = np.logspace(-2, 1, 100)
m = np.asarray([5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0])
gamma = 19
dGWT = 2.0

pmp = smt.smt_model_fragility(depth, qt, fs, amax, m, dGWT=dGWT, gamma=gamma)
fig, ax = plt.subplots()
for i in range(len(m)):
    ax.plot(amax, pmp[i], label=m[i])
ax.set_xscale('log')
ax.grid(True, alpha=0.5, which='both')
ax.legend(title='M')
ax.set_xlabel(r'$a_{max}$ [g]')
ax.set_ylabel(r'$P[M_P]$')
ax.set_ylim(-0.05, 1.05)