ge.dynamics — Soil dynamics & liquefaction¶

Three things: small-strain stiffness (\(G_{\max}\)), the modulus-reduction & damping curves used in equivalent-linear ground response analysis, and the simplified-procedure liquefaction triggering analysis (Seed-Idriss / NCEER / Idriss-Boulanger).

Small-strain shear modulus¶

Python

import geoeq as ge

# From shear-wave velocity
ge.gmax(Vs=200, gamma=20)                       # 81550 kPa = 81.5 MPa

# From Hardin & Black (1968)
ge.gmax_hardin(e=0.6, sigma_m_eff=100, OCR=2, PI=30,
                soil_type="round_grained")

Modulus reduction G/G_max and damping D¶

Darendeli (2001) is the modern empirical model:

\[ \frac{G}{G_{\max}} = \frac{1}{1 + (\gamma/\gamma_r)^{0.92}} \]

where the reference strain \(\gamma_r\) depends on mean effective stress, plasticity index, and OCR.

Python

ge.modulus_reduction(strain=1e-3, sigma_m_eff=100, PI=15, OCR=1)
ge.damping_ratio(strain=1e-3,    sigma_m_eff=100, PI=15, OCR=1)

G/Gmax curves — Modulus-reduction (solid) and damping (dashed) curves for a family of plasticity indices. Darendeli (2001) at σ'_m = 100 kPa, OCR = 1. Generated by `ge.gmax_curves_plot()`.

Liquefaction triggering — simplified procedure¶

The complete workflow in five steps. Each step is one function call.

Step 1: depth-reduction factor \(r_d\)¶

Python

ge.depth_reduction(z=6, method="idriss_1999", Mw=7.0)   # 0.971

Other methods: "liao_whitman_1986", "cetin_2004".

Step 2: cyclic stress ratio (CSR)¶

\[ CSR = 0.65\,(a_{\max}/g)\,(\sigma_v/\sigma'_v)\,r_d \]

Python

csr = ge.liquefaction_csr(amax=0.25, sigma_v=120,
                          sigma_v_eff=70, z=6, Mw=7.0)
# {'CSR': 0.226, 'rd': 0.971, 'amax_g': 0.25}

Step 3: cyclic resistance ratio (CRR)¶

Four families. Pick by what you have:

Python

# NCEER, SPT (Youd et al. 2001) — most widely cited
ge.liquefaction_crr(N160cs=12, method="youd_2001")

# Updated SPT (Idriss & Boulanger 2008)
ge.liquefaction_crr(N160cs=15, method="idriss_boulanger_2008")

# Updated CPT (Idriss & Boulanger 2008)
ge.liquefaction_crr(qc1Ncs=100, method="idriss_boulanger_2008_cpt")

# Vs-based (Andrus & Stokoe 2000)
ge.liquefaction_crr(Vs1=150, FC=10, method="andrus_stokoe_2000")

Step 4: magnitude scaling factor (MSF)¶

Python

ge.magnitude_scaling_factor(Mw=7.0, method="idriss_1999")
ge.magnitude_scaling_factor(Mw=7.0, method="nceer")
ge.magnitude_scaling_factor(Mw=7.0, method="boulanger_idriss_2014")

Step 5: factor of safety¶

\[ FS_L = \frac{CRR \cdot MSF \cdot K_\sigma \cdot K_\alpha}{CSR} \]

Python

fs = ge.liquefaction_fos(
    CSR=csr["CSR"], CRR=0.12, Mw=7.0,
    K_sigma=1.0, K_alpha=1.0,
)
# {'FS': 0.69, 'liquefies': True, 'MSF': 1.19, ...}

The triggering chart¶

Liquefaction triggering chart. CRR curves from Youd et al. (2001) and Idriss & Boulanger (2008). Field data points coloured by their computed FS: red = liquefying, green = safe. Generated by `ge.liquefaction_chart()`.

Response spectrum¶

Python

ge.response_spectrum(T=[0.01, 0.1, 0.5, 1, 2, 4],
                      Sa=[0.40, 0.62, 0.55, 0.32, 0.18, 0.10],
                      site_class="D")

API reference¶

gmax ¶

Python

gmax(Vs: float, gamma: float = 18.0) -> float

Small-strain shear modulus Gmax from shear-wave velocity Vs.

Text Only

Gmax = rho * Vs^2    (kPa)

PARAMETER	DESCRIPTION
`Vs`	Shear-wave velocity (m/s). TYPE: `float`
`gamma`	Total unit weight (kN/m^3). Default 18. TYPE: `float` DEFAULT: `18.0`

RETURNS	DESCRIPTION
`Gmax`	Small-strain shear modulus (kPa). TYPE: `float`

Reference

Stokoe et al. (1999); Das (2010) Eq. 11.18.

Source code in geoeq/dynamics/modulus.py

Python
def gmax(Vs: float, gamma: float = 18.0) -> float:
    """Small-strain shear modulus Gmax from shear-wave velocity Vs.

        Gmax = rho * Vs^2    (kPa)

    Parameters
    ----------
    Vs : float
        Shear-wave velocity (m/s).
    gamma : float
        Total unit weight (kN/m^3). Default 18.

    Returns
    -------
    Gmax : float
        Small-strain shear modulus (kPa).

    Reference
    ---------
    Stokoe et al. (1999); Das (2010) Eq. 11.18.
    """
    check_positive(Vs, "Vs")
    check_positive(gamma, "gamma")
    rho = gamma * 1000 / 9.81  # kN/m^3 -> kg/m^3
    return float(rho * Vs ** 2 / 1000)  # Pa -> kPa

gmax_hardin ¶

Python

gmax_hardin(e: float, sigma_m_eff: float, OCR: float = 1.0, PI: float = 0, soil_type: str = 'round_grained') -> float

Gmax from Hardin & Black (1968) / Hardin & Drnevich (1972).

Text Only

Gmax = A * F(e) * OCR^k * (sigma'_m / pa)^n * pa  (kPa)

For most soils n ~ 0.5 and:

Round-grained sands: A = 6900, F(e) = (2.97 - e)^2 / (1 + e)
Angular-grained sands: A = 3300, F(e) = (2.17 - e)^2 / (1 + e)
Clays: A = 3230, F(e) = (2.97 - e)^2 / (1 + e)

PARAMETER	DESCRIPTION
`e`	Void ratio (-). TYPE: `float`
`sigma_m_eff`	Mean effective confining stress (kPa). TYPE: `float`
`OCR`	Over-consolidation ratio. TYPE: `float` DEFAULT: `1.0`
`PI`	Plasticity index. Sets k via Hardin-Drnevich table: PI=0: k=0 PI=20: k=0.18 PI=40: k=0.30 PI=60: k=0.41 PI=80: k=0.48 PI>=100: k=0.50 TYPE: `float` DEFAULT: `0`
`soil_type`	'round_grained' \| 'angular_grained' \| 'clay'. TYPE: `str` DEFAULT: `'round_grained'`

RETURNS	DESCRIPTION
`Gmax`	Small-strain shear modulus (kPa). TYPE: `float`

Reference

Hardin & Black (1968); Hardin & Drnevich (1972); Kramer (1996) Eq. 6.9.

Source code in geoeq/dynamics/modulus.py

Python
def gmax_hardin(e: float, sigma_m_eff: float, OCR: float = 1.0,
                PI: float = 0, soil_type: str = "round_grained") -> float:
    """Gmax from Hardin & Black (1968) / Hardin & Drnevich (1972).

        Gmax = A * F(e) * OCR^k * (sigma'_m / pa)^n * pa  (kPa)

    For most soils n ~ 0.5 and:

    * Round-grained sands:    A = 6900,  F(e) = (2.97 - e)^2 / (1 + e)
    * Angular-grained sands:  A = 3300,  F(e) = (2.17 - e)^2 / (1 + e)
    * Clays:                  A = 3230,  F(e) = (2.97 - e)^2 / (1 + e)

    Parameters
    ----------
    e : float
        Void ratio (-).
    sigma_m_eff : float
        Mean effective confining stress (kPa).
    OCR : float
        Over-consolidation ratio.
    PI : float
        Plasticity index. Sets k via Hardin-Drnevich table:
            PI=0:   k=0
            PI=20:  k=0.18
            PI=40:  k=0.30
            PI=60:  k=0.41
            PI=80:  k=0.48
            PI>=100: k=0.50
    soil_type : str
        'round_grained' | 'angular_grained' | 'clay'.

    Returns
    -------
    Gmax : float
        Small-strain shear modulus (kPa).

    Reference
    ---------
    Hardin & Black (1968); Hardin & Drnevich (1972); Kramer (1996) Eq. 6.9.
    """
    check_positive(e, "e")
    check_positive(sigma_m_eff, "sigma_m_eff")
    check_positive(OCR, "OCR")
    soil_type = soil_type.lower()
    if soil_type == "round_grained":
        A, F = 6900, (2.97 - e) ** 2 / (1 + e)
    elif soil_type in ("angular_grained", "angular"):
        A, F = 3300, (2.17 - e) ** 2 / (1 + e)
    elif soil_type == "clay":
        A, F = 3230, (2.97 - e) ** 2 / (1 + e)
    else:
        raise ValueError(
            "soil_type must be 'round_grained', 'angular_grained', or 'clay'.")
    # k from PI (Kramer 1996 Table 6.4)
    k_table = [(0, 0), (20, 0.18), (40, 0.30),
               (60, 0.41), (80, 0.48), (100, 0.50)]
    PI_clip = min(max(PI, 0), 100)
    # piecewise linear
    for (pi1, k1), (pi2, k2) in zip(k_table[:-1], k_table[1:]):
        if pi1 <= PI_clip <= pi2:
            k = k1 + (k2 - k1) * (PI_clip - pi1) / (pi2 - pi1)
            break
    else:
        k = 0.5
    pa = 101.325  # atmospheric pressure, kPa
    n = 0.5
    return float(A * F * OCR ** k * (sigma_m_eff / pa) ** n * pa / 100)

modulus_reduction ¶

Python

modulus_reduction(strain: float, sigma_m_eff: float = 100.0, PI: float = 15, OCR: float = 1.0, method: str = 'darendeli') -> float

Modulus-reduction ratio G/Gmax at shear strain strain (%).

Text Only

G / Gmax = 1 / (1 + (gamma / gamma_r)^a)

where gamma_r is the reference strain (Darendeli 2001) and a ~ 0.92 for typical soils.

PARAMETER	DESCRIPTION
`strain`	Shear strain (decimal, e.g. 0.001 = 0.1%). TYPE: `float`
`sigma_m_eff`	Mean effective confining stress (kPa). Default 100. TYPE: `float` DEFAULT: `100.0`
`PI`	Plasticity index. Default 15. TYPE: `float` DEFAULT: `15`
`OCR`	Over-consolidation ratio. Default 1. TYPE: `float` DEFAULT: `1.0`
`method`	'darendeli' (2001) or 'vucetic_dobry' (1991 chart values). TYPE: `str` DEFAULT: `'darendeli'`

RETURNS	DESCRIPTION
`G_over_Gmax`	Modulus reduction ratio (-). TYPE: `float`

Reference

Darendeli (2001) Eq. 8.10-8.12; Vucetic & Dobry (1991).

Source code in geoeq/dynamics/modulus.py

Python
def modulus_reduction(
    strain: float, sigma_m_eff: float = 100.0, PI: float = 15,
    OCR: float = 1.0, method: str = "darendeli",
) -> float:
    """Modulus-reduction ratio G/Gmax at shear strain ``strain`` (%).

        G / Gmax = 1 / (1 + (gamma / gamma_r)^a)

    where gamma_r is the reference strain (Darendeli 2001) and a ~ 0.92
    for typical soils.

    Parameters
    ----------
    strain : float
        Shear strain (decimal, e.g. 0.001 = 0.1%).
    sigma_m_eff : float
        Mean effective confining stress (kPa). Default 100.
    PI : float
        Plasticity index. Default 15.
    OCR : float
        Over-consolidation ratio. Default 1.
    method : str
        'darendeli' (2001) or 'vucetic_dobry' (1991 chart values).

    Returns
    -------
    G_over_Gmax : float
        Modulus reduction ratio (-).

    Reference
    ---------
    Darendeli (2001) Eq. 8.10-8.12; Vucetic & Dobry (1991).
    """
    check_non_negative(strain, "strain")
    method = method.lower()
    gamma_pct = strain * 100  # decimal -> %
    if method == "darendeli":
        gamma_r = _darendeli_gamma_r(sigma_m_eff, PI, OCR)
        a = 0.92
        return float(1.0 / (1.0 + (gamma_pct / gamma_r) ** a))
    if method in ("vucetic_dobry", "vd"):
        # Vucetic & Dobry (1991) for PI=0 sand: digitized fit.
        # G/Gmax ~ 1/(1 + (gamma/gamma_05)^1) with gamma_05 from PI:
        gamma_05 = 0.01 + 0.001 * PI  # very rough fit
        return float(1.0 / (1.0 + (gamma_pct / gamma_05)))
    raise ValueError("method must be 'darendeli' or 'vucetic_dobry'.")

damping_ratio ¶

Python

damping_ratio(strain: float, sigma_m_eff: float = 100.0, PI: float = 15, OCR: float = 1.0, freq_Hz: float = 1.0, N_cycles: float = 10, method: str = 'darendeli') -> float

Equivalent-linear damping ratio D (decimal, e.g. 0.05 = 5%).

Uses Darendeli (2001) full model: D = D_min + b * (G/Gmax)^0.1 * (D_masing - D_min)/100 [%] D_min = (phi6 + phi7 * PI * OCR^phi8) * sigma_norm^phi9 * (1 + phi10 * ln(freq))

Returns the small-strain D for very small strains, with the Masing correction at larger strains.

Reference

Darendeli (2001) Eq. 8.13-8.18.

Source code in geoeq/dynamics/modulus.py

Python
def damping_ratio(
    strain: float, sigma_m_eff: float = 100.0, PI: float = 15,
    OCR: float = 1.0, freq_Hz: float = 1.0, N_cycles: float = 10,
    method: str = "darendeli",
) -> float:
    """Equivalent-linear damping ratio D (decimal, e.g. 0.05 = 5%).

    Uses Darendeli (2001) full model:
        D = D_min + b * (G/Gmax)^0.1 * (D_masing - D_min)/100   [%]
        D_min = (phi6 + phi7 * PI * OCR^phi8) * sigma_norm^phi9
                * (1 + phi10 * ln(freq))

    Returns the small-strain D for very small strains, with the Masing
    correction at larger strains.

    Reference
    ---------
    Darendeli (2001) Eq. 8.13-8.18.
    """
    check_non_negative(strain, "strain")
    pa = 101.325
    sigma_norm = sigma_m_eff / pa
    # Coefficients from Darendeli (2001) Table 8.12 (selected):
    phi6, phi7, phi8, phi9, phi10 = 0.8005, 0.0129, -0.1069, -0.2889, 0.2919
    phi11 = 0.6329
    phi12 = -0.0057

    method = method.lower()
    # Small-strain damping (%), Darendeli Eq. 8.18.
    D_min = ((phi6 + phi7 * PI * OCR ** phi8)
             * sigma_norm ** phi9 * (1 + phi10 * np.log(freq_Hz)))
    G_Gmax = modulus_reduction(strain, sigma_m_eff, PI, OCR,
                                method=method if method != "darendeli"
                                else "darendeli")
    # Use the Hardin & Drnevich (1972) hyperbolic relation:
    #     D / D_max = 1 - G/Gmax
    # to avoid the well-known singularity of the closed-form Masing
    # damping at gamma -> 0. D_max ~ 25 % for typical sands, lower for
    # high-PI clays.
    D_max_pct = max(20.0, 25.0 - 0.05 * PI)  # %
    D_pct = D_min + (D_max_pct - D_min) * (1.0 - G_Gmax)
    return float(D_pct / 100.0)  # % -> decimal

depth_reduction ¶

Python

depth_reduction(z: float, method: str = 'idriss_1999', Mw: float = 7.5) -> float

Stress-reduction factor rd that accounts for the flexibility of the soil column under earthquake shear stress.

PARAMETER	DESCRIPTION
`z`	Depth (m). TYPE: `float`
`method`	'liao_whitman_1986' -- simple bilinear (deprecated by NCEER). 'idriss_1999' -- magnitude-dependent (NCEER-recommended). 'cetin_2004' -- empirical, magnitude-dependent. TYPE: `str` DEFAULT: `'idriss_1999'`
`Mw`	Moment magnitude (only used by idriss_1999/cetin). TYPE: `float` DEFAULT: `7.5`

RETURNS	DESCRIPTION
`rd`	TYPE: `float`

Reference

Idriss (1999); Youd et al. (2001) Eq. 4.

Source code in geoeq/dynamics/liquefaction.py

Python
def depth_reduction(z: float, method: str = "idriss_1999",
                     Mw: float = 7.5) -> float:
    """Stress-reduction factor rd that accounts for the flexibility of the
    soil column under earthquake shear stress.

    Parameters
    ----------
    z : float
        Depth (m).
    method : str
        'liao_whitman_1986' -- simple bilinear (deprecated by NCEER).
        'idriss_1999'       -- magnitude-dependent (NCEER-recommended).
        'cetin_2004'        -- empirical, magnitude-dependent.
    Mw : float
        Moment magnitude (only used by idriss_1999/cetin).

    Returns
    -------
    rd : float

    Reference
    ---------
    Idriss (1999); Youd et al. (2001) Eq. 4.
    """
    check_non_negative(z, "z")
    method = method.lower()
    if method == "liao_whitman_1986":
        if z <= 9.15:
            return float(1 - 0.00765 * z)
        if z <= 23:
            return float(1.174 - 0.0267 * z)
        return float(0.744 - 0.008 * z)
    if method == "idriss_1999":
        # Idriss (1999) magnitude-dependent rd (Youd et al. 2001 Eq. 4)
        alpha = -1.012 - 1.126 * np.sin(z / 11.73 + 5.133)
        beta = 0.106 + 0.118 * np.sin(z / 11.28 + 5.142)
        rd = np.exp(alpha + beta * Mw)
        return float(rd)
    if method == "cetin_2004":
        # Cetin et al. (2004) -- depth- and magnitude-dependent.
        num = (1 + (-23.013 - 2.949 * z + 0.999 * Mw + 0.0525 * z * Mw)
               / (16.258 + 0.201 * np.exp(0.341 * (-z + 0.0785 * Mw + 7.586))))
        den = (1 + (-23.013 + 0.999 * Mw)
               / (16.258 + 0.201 * np.exp(0.341 * (0.0785 * Mw + 7.586))))
        return float(num / den)
    raise ValueError("method must be one of 'liao_whitman_1986', "
                     "'idriss_1999', or 'cetin_2004'.")

liquefaction_csr ¶

Python

liquefaction_csr(amax: float, sigma_v: float, sigma_v_eff: float, Mw: float = 7.5, z: float = None, rd: float = None, g: float = 9.81) -> dict

Cyclic stress ratio induced by the earthquake.

Text Only

CSR = 0.65 * (amax / g) * (sigma_v / sigma_v_eff) * rd

PARAMETER	DESCRIPTION
`amax`	Peak horizontal ground acceleration as a fraction of g (e.g. 0.25 for 0.25g). TYPE: `float`
`sigma_v`	Total vertical stress at depth (kPa). TYPE: `float`
`sigma_v_eff`	Effective vertical stress at depth (kPa). TYPE: `float`
`Mw`	Moment magnitude (used only for rd if rd is None). TYPE: `float` DEFAULT: `7.5`
`z`	Depth (m). Required if rd is None. TYPE: `float` DEFAULT: `None`
`rd`	Pre-computed depth-reduction factor. If None, computed from z. TYPE: `float` DEFAULT: `None`

RETURNS	DESCRIPTION
`dict`	`{'CSR': ..., 'rd': ...}`.

Reference

Seed & Idriss (1971); Youd et al. (2001) Eq. 1.

Source code in geoeq/dynamics/liquefaction.py

Python
def liquefaction_csr(
    amax: float, sigma_v: float, sigma_v_eff: float,
    Mw: float = 7.5, z: float = None, rd: float = None,
    g: float = 9.81,
) -> dict:
    """Cyclic stress ratio induced by the earthquake.

        CSR = 0.65 * (amax / g) * (sigma_v / sigma_v_eff) * rd

    Parameters
    ----------
    amax : float
        Peak horizontal ground acceleration as a fraction of g
        (e.g. 0.25 for 0.25g).
    sigma_v : float
        Total vertical stress at depth (kPa).
    sigma_v_eff : float
        Effective vertical stress at depth (kPa).
    Mw : float
        Moment magnitude (used only for rd if rd is None).
    z : float
        Depth (m). Required if rd is None.
    rd : float
        Pre-computed depth-reduction factor. If None, computed from z.

    Returns
    -------
    dict
        ``{'CSR': ..., 'rd': ...}``.

    Reference
    ---------
    Seed & Idriss (1971); Youd et al. (2001) Eq. 1.
    """
    check_positive(sigma_v, "sigma_v")
    check_positive(sigma_v_eff, "sigma_v_eff")
    check_positive(amax, "amax")
    amax_over_g = float(amax)
    if rd is None:
        if z is None:
            raise ValueError("Provide rd or z (and Mw).")
        rd = depth_reduction(z, method="idriss_1999", Mw=Mw)
    CSR = 0.65 * amax_over_g * (sigma_v / sigma_v_eff) * rd
    return {"CSR": float(CSR), "rd": float(rd), "amax_g": float(amax_over_g)}

liquefaction_crr ¶

Python

liquefaction_crr(N160cs: float = None, qc1Ncs: float = None, Vs1: float = None, FC: float = 0, method: str = 'youd_2001', Mw: float = 7.5) -> dict

Cyclic resistance ratio at Mw=7.5 from SPT, CPT, or Vs.

Methods: * youd_2001 -- NCEER SPT (Youd et al. 2001 Eq. 5): CRR7.5 = 1/(34 - N160cs) + N160cs/135 + 50/(10N160cs+45)^2 - 1/200 * idriss_boulanger_2008 -- updated SPT (Eq. 4): CRR7.5 = exp[ N160cs/14.1 + (N160cs/126)^2 - (N160cs/23.6)^3 + (N160cs/25.4)^4 - 2.8 ] * idriss_boulanger_2008_cpt -- CPT-based (Eq. 6): CRR7.5 = exp[ qc1Ncs/540 + (qc1Ncs/67)^2 - (qc1Ncs/80)^3 + (qc1Ncs/114)^4 - 3 ] * andrus_stokoe_2000 -- Vs-based (Eq. 9): CRR7.5 = 0.022(Vs1/100)^2 + 2.8(1/(215-Vs1) - 1/215)

RETURNS	DESCRIPTION
`dict with 'CRR', 'method'.`

Reference

Youd et al. (2001); Idriss & Boulanger (2008); Andrus & Stokoe (2000).

Source code in geoeq/dynamics/liquefaction.py

Python
def liquefaction_crr(
    N160cs: float = None, qc1Ncs: float = None,
    Vs1: float = None, FC: float = 0,
    method: str = "youd_2001", Mw: float = 7.5,
) -> dict:
    """Cyclic resistance ratio at Mw=7.5 from SPT, CPT, or Vs.

    Methods:
    * **youd_2001** -- NCEER SPT (Youd et al. 2001 Eq. 5):
        CRR7.5 = 1/(34 - N160cs) + N160cs/135 + 50/(10*N160cs+45)^2 - 1/200
    * **idriss_boulanger_2008** -- updated SPT (Eq. 4):
        CRR7.5 = exp[ N160cs/14.1 + (N160cs/126)^2 - (N160cs/23.6)^3
                    + (N160cs/25.4)^4 - 2.8 ]
    * **idriss_boulanger_2008_cpt** -- CPT-based (Eq. 6):
        CRR7.5 = exp[ qc1Ncs/540 + (qc1Ncs/67)^2 - (qc1Ncs/80)^3
                    + (qc1Ncs/114)^4 - 3 ]
    * **andrus_stokoe_2000** -- Vs-based (Eq. 9):
        CRR7.5 = 0.022(Vs1/100)^2 + 2.8*(1/(215-Vs1) - 1/215)

    Returns
    -------
    dict with 'CRR', 'method'.

    Reference
    ---------
    Youd et al. (2001); Idriss & Boulanger (2008); Andrus & Stokoe (2000).
    """
    method = method.lower()
    if method in ("youd_2001", "nceer", "youd"):
        if N160cs is None:
            raise ValueError("N160cs required for Youd 2001.")
        if N160cs >= 30:
            CRR = 2.0  # capped (non-liquefiable)
        else:
            CRR = (1 / (34 - N160cs)
                   + N160cs / 135
                   + 50 / (10 * N160cs + 45) ** 2
                   - 1 / 200)
        out = {"CRR": float(CRR), "method": "youd_2001"}
    elif method in ("idriss_boulanger_2008", "ib2008", "ib"):
        if N160cs is None:
            raise ValueError("N160cs required for IB2008 SPT.")
        x = N160cs
        if x > 37.5:  # cap at upper bound
            x = 37.5
        CRR = np.exp(x / 14.1
                     + (x / 126) ** 2
                     - (x / 23.6) ** 3
                     + (x / 25.4) ** 4 - 2.8)
        out = {"CRR": float(CRR), "method": "idriss_boulanger_2008"}
    elif method in ("idriss_boulanger_2008_cpt", "ib2008_cpt", "cpt"):
        if qc1Ncs is None:
            raise ValueError("qc1Ncs required for IB2008 CPT.")
        x = qc1Ncs
        CRR = np.exp(x / 540
                     + (x / 67) ** 2
                     - (x / 80) ** 3
                     + (x / 114) ** 4 - 3.0)
        out = {"CRR": float(CRR), "method": "idriss_boulanger_2008_cpt"}
    elif method in ("andrus_stokoe_2000", "as2000", "vs"):
        if Vs1 is None:
            raise ValueError("Vs1 required for Andrus & Stokoe 2000.")
        if FC <= 5:
            Vs1_star = 215
        elif FC < 35:
            Vs1_star = 215 - 0.5 * (FC - 5)
        else:
            Vs1_star = 200
        if Vs1 >= Vs1_star:
            CRR = 2.0
        else:
            CRR = 0.022 * (Vs1 / 100) ** 2 + 2.8 * (1 / (Vs1_star - Vs1) - 1 / Vs1_star)
        out = {"CRR": float(CRR), "method": "andrus_stokoe_2000",
               "Vs1_star": float(Vs1_star)}
    else:
        raise ValueError(
            "method must be 'youd_2001', 'idriss_boulanger_2008', "
            "'idriss_boulanger_2008_cpt', or 'andrus_stokoe_2000'.")
    return out

magnitude_scaling_factor ¶

Python

magnitude_scaling_factor(Mw: float, method: str = 'idriss_1999') -> float

Magnitude scaling factor MSF to adjust CRR_7.5 to other magnitudes.

Methods: * idriss_1999: MSF = 6.9 * exp(-Mw/4) - 0.058 (capped at 1.8) * nceer: MSF = (10^2.24)/Mw2.56 (Youd et al. 2001 Eq. 24) * boulanger_idriss_2014: MSF = 1 + (MSFmax - 1)*(8.64 * exp(-Mw/4) - 1.325)

Reference

Idriss (1999); Youd et al. (2001); Boulanger & Idriss (2014).

Source code in geoeq/dynamics/liquefaction.py

Python
def magnitude_scaling_factor(Mw: float,
                              method: str = "idriss_1999") -> float:
    """Magnitude scaling factor MSF to adjust CRR_7.5 to other magnitudes.

    Methods:
    * **idriss_1999**:  MSF = 6.9 * exp(-Mw/4) - 0.058  (capped at 1.8)
    * **nceer**:        MSF = (10^2.24)/Mw^2.56 (Youd et al. 2001 Eq. 24)
    * **boulanger_idriss_2014**:
        MSF = 1 + (MSFmax - 1)*(8.64 * exp(-Mw/4) - 1.325)

    Reference
    ---------
    Idriss (1999); Youd et al. (2001); Boulanger & Idriss (2014).
    """
    check_positive(Mw, "Mw")
    method = method.lower()
    if method in ("idriss_1999", "idriss"):
        return float(min(1.8, 6.9 * np.exp(-Mw / 4) - 0.058))
    if method in ("nceer", "youd_2001"):
        return float(10 ** 2.24 / Mw ** 2.56)
    if method in ("boulanger_idriss_2014", "bi2014"):
        MSF_max = 1.8  # for sands
        return float(1 + (MSF_max - 1) * (8.64 * np.exp(-Mw / 4) - 1.325))
    raise ValueError(
        "method must be 'idriss_1999', 'nceer', or 'boulanger_idriss_2014'.")

liquefaction_fos ¶

Python

liquefaction_fos(CSR: float, CRR: float, Mw: float = 7.5, MSF: float = None, K_sigma: float = 1.0, K_alpha: float = 1.0, MSF_method: str = 'idriss_1999') -> dict

Factor of safety against liquefaction triggering.

Text Only

FS_L = (CRR * MSF * K_sigma * K_alpha) / CSR

PARAMETER	DESCRIPTION
`CSR`	Cyclic stress ratio from earthquake loading. TYPE: `float`
`CRR`	Cyclic resistance ratio at Mw=7.5. TYPE: `float`
`Mw`	Moment magnitude (for MSF if not provided). TYPE: `float` DEFAULT: `7.5`
`MSF`	Pre-computed magnitude scaling factor. TYPE: `float` DEFAULT: `None`
`K_sigma`	Overburden correction (default 1). TYPE: `float` DEFAULT: `1.0`
`K_alpha`	Static shear stress correction (default 1). TYPE: `float` DEFAULT: `1.0`

RETURNS	DESCRIPTION
`dict`	`{'FS': ..., 'MSF': ..., 'liquefies': bool}`.

Reference

Youd et al. (2001) Eq. 26; Idriss & Boulanger (2008).

Source code in geoeq/dynamics/liquefaction.py

Python
def liquefaction_fos(
    CSR: float, CRR: float, Mw: float = 7.5,
    MSF: float = None, K_sigma: float = 1.0, K_alpha: float = 1.0,
    MSF_method: str = "idriss_1999",
) -> dict:
    """Factor of safety against liquefaction triggering.

        FS_L = (CRR * MSF * K_sigma * K_alpha) / CSR

    Parameters
    ----------
    CSR : float
        Cyclic stress ratio from earthquake loading.
    CRR : float
        Cyclic resistance ratio at Mw=7.5.
    Mw : float
        Moment magnitude (for MSF if not provided).
    MSF : float, optional
        Pre-computed magnitude scaling factor.
    K_sigma : float
        Overburden correction (default 1).
    K_alpha : float
        Static shear stress correction (default 1).

    Returns
    -------
    dict
        ``{'FS': ..., 'MSF': ..., 'liquefies': bool}``.

    Reference
    ---------
    Youd et al. (2001) Eq. 26; Idriss & Boulanger (2008).
    """
    check_positive(CSR, "CSR")
    check_positive(CRR, "CRR")
    if MSF is None:
        MSF = magnitude_scaling_factor(Mw, method=MSF_method)
    FS = CRR * MSF * K_sigma * K_alpha / CSR
    return {
        "FS": float(FS), "MSF": float(MSF),
        "K_sigma": float(K_sigma), "K_alpha": float(K_alpha),
        "liquefies": bool(FS < 1.0),
    }

response_spectrum ¶

Python

response_spectrum(T: Sequence[float], Sa: Sequence[float], site_class: str = None, ax=None, save_as: str = None, **kwargs)

Plot a response spectrum.

PARAMETER	DESCRIPTION
`T`	Periods (s). TYPE: `sequence`
`Sa`	Spectral accelerations (g). TYPE: `sequence`
`site_class`	Label printed on the figure (e.g. 'Site Class D'). TYPE: `str` DEFAULT: `None`
`ax`	Embed in an existing axis. TYPE: `matplotlib axis` DEFAULT: `None`
`save_as`	Path to save the figure. TYPE: `str` DEFAULT: `None`
`**kwargs`	Passed to Matplotlib `plot`. DEFAULT: `{}`

RETURNS	DESCRIPTION
`Figure`

Source code in geoeq/dynamics/response.py

Python
def response_spectrum(T: Sequence[float], Sa: Sequence[float],
                       site_class: str = None, ax=None,
                       save_as: str = None, **kwargs):
    """Plot a response spectrum.

    Parameters
    ----------
    T : sequence
        Periods (s).
    Sa : sequence
        Spectral accelerations (g).
    site_class : str, optional
        Label printed on the figure (e.g. 'Site Class D').
    ax : matplotlib axis, optional
        Embed in an existing axis.
    save_as : str, optional
        Path to save the figure.
    **kwargs
        Passed to Matplotlib ``plot``.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt
    T = np.asarray(T)
    Sa = np.asarray(Sa)
    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 5))
    else:
        fig = ax.figure
    style = {"color": "#1f3a93", "linewidth": 1.8}
    style.update(kwargs)
    ax.plot(T, Sa, **style)
    ax.set_xlabel("Period $T$ (s)")
    ax.set_ylabel("Spectral acceleration $S_a$ (g)")
    title = "Response spectrum"
    if site_class:
        title += f" ({site_class})"
    ax.set_title(title)
    ax.set_xscale("log")
    ax.grid(True, which="both", alpha=0.3)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

gmax_curves_plot ¶

Python

gmax_curves_plot(PI_values=(0, 15, 30, 50, 100), sigma_m_eff: float = 100.0, OCR: float = 1.0, strain_range=None, method: str = 'darendeli', ax=None, save_as: str = None)

Plot G/Gmax (left axis) and damping D (right axis) versus shear strain for a family of plasticity indices.

Reproduces the Vucetic-Dobry style family of curves using Darendeli (2001)'s closed-form expressions.

Source code in geoeq/dynamics/plots.py

Python
def gmax_curves_plot(
    PI_values=(0, 15, 30, 50, 100),
    sigma_m_eff: float = 100.0, OCR: float = 1.0,
    strain_range=None, method: str = "darendeli",
    ax=None, save_as: str = None,
):
    """Plot G/Gmax (left axis) and damping D (right axis) versus
    shear strain for a family of plasticity indices.

    Reproduces the Vucetic-Dobry style family of curves using
    Darendeli (2001)'s closed-form expressions.
    """
    import matplotlib.pyplot as plt
    if strain_range is None:
        strain_range = np.logspace(-6, -2, 80)  # 0.0001 % .. 1 %
    if ax is None:
        fig, ax1 = plt.subplots(figsize=(8, 5))
    else:
        ax1 = ax
        fig = ax1.figure
    ax2 = ax1.twinx()

    cmap = plt.cm.viridis(np.linspace(0.2, 0.9, len(PI_values)))
    for PI, colour in zip(PI_values, cmap):
        G_ratio = [modulus_reduction(s, sigma_m_eff, PI, OCR, method)
                   for s in strain_range]
        D = [damping_ratio(s, sigma_m_eff, PI, OCR, method=method) * 100
             for s in strain_range]
        ax1.plot(strain_range * 100, G_ratio, color=colour, linewidth=2,
                 label=f"PI = {PI}")
        ax2.plot(strain_range * 100, D, color=colour, linewidth=1.2,
                 linestyle="--", alpha=0.7)

    ax1.set_xscale("log")
    ax1.set_xlim(1e-4, 1.0)
    ax1.set_xlabel(r"Cyclic shear strain $\gamma$ (%)")
    ax1.set_ylabel(r"$G / G_{\max}$ (solid)")
    ax2.set_ylabel(r"Damping ratio $D$ (\%), dashed")
    ax1.set_ylim(0, 1.05)
    ax2.set_ylim(0, 30)
    ax1.set_title(
        f"Modulus reduction + damping — Darendeli (2001), "
        f"$\\sigma'_m={sigma_m_eff}$ kPa, OCR$={OCR}$")
    ax1.grid(True, which="both", alpha=0.3)
    ax1.legend(loc="center left", fontsize=9)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

liquefaction_chart ¶

Python

liquefaction_chart(data_N: Sequence[float] = None, data_CSR: Sequence[float] = None, data_FS: Sequence[float] = None, Mw: float = 7.5, methods=('youd_2001', 'idriss_boulanger_2008'), ax=None, save_as: str = None)

Liquefaction-triggering chart: CRR curve(s) plus optional CSR data.

PARAMETER	DESCRIPTION
`data_N`	Field data points to overlay (corrected blow counts and CSR). TYPE: `sequences` DEFAULT: `None`
`data_CSR`	Field data points to overlay (corrected blow counts and CSR). TYPE: `sequences` DEFAULT: `None`
`data_FS`	Per-point factor of safety; if given, points are coloured red (FS < 1) / green (FS >= 1). TYPE: `sequence` DEFAULT: `None`
`methods`	CRR curves to draw. Default both Youd-2001 and IB-2008. TYPE: `sequence of str` DEFAULT: `('youd_2001', 'idriss_boulanger_2008')`

RETURNS	DESCRIPTION
`Figure`

Source code in geoeq/dynamics/plots.py

Python
def liquefaction_chart(
    data_N: Sequence[float] = None,
    data_CSR: Sequence[float] = None,
    data_FS: Sequence[float] = None,
    Mw: float = 7.5,
    methods=("youd_2001", "idriss_boulanger_2008"),
    ax=None, save_as: str = None,
):
    """Liquefaction-triggering chart: CRR curve(s) plus optional CSR data.

    Parameters
    ----------
    data_N, data_CSR : sequences
        Field data points to overlay (corrected blow counts and CSR).
    data_FS : sequence, optional
        Per-point factor of safety; if given, points are coloured
        red (FS < 1) / green (FS >= 1).
    methods : sequence of str
        CRR curves to draw. Default both Youd-2001 and IB-2008.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt
    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 6))
    else:
        fig = ax.figure

    N_arr = np.linspace(1, 30, 200)
    colours = {"youd_2001": "#1f3a93",
               "idriss_boulanger_2008": "#c0392b"}
    labels = {"youd_2001": "NCEER (Youd et al. 2001)",
              "idriss_boulanger_2008": "Idriss \\& Boulanger (2008)"}
    for m in methods:
        crr = []
        for N in N_arr:
            res = liquefaction_crr(N160cs=N, method=m)
            crr.append(res["CRR"])
        ax.plot(N_arr, crr, color=colours.get(m, "k"),
                linewidth=2.0, label=labels.get(m, m))

    # Shaded liquefaction zone
    ax.fill_between(N_arr, 0, [liquefaction_crr(N160cs=n, method="youd_2001")["CRR"]
                                for n in N_arr],
                    color="#c0392b", alpha=0.08, label="Liquefaction zone")

    # Field data overlay
    if data_N is not None and data_CSR is not None:
        if data_FS is not None:
            colors_pts = ["#c0392b" if fs < 1 else "#27ae60" for fs in data_FS]
        else:
            colors_pts = "#000000"
        ax.scatter(data_N, data_CSR, c=colors_pts, s=70,
                    edgecolors="black", linewidths=0.8, zorder=5)

    ax.set_xlim(0, 30)
    ax.set_ylim(0, 0.6)
    ax.set_xlabel(r"Corrected SPT blow count $(N_1)_{60,cs}$")
    ax.set_ylabel(r"CSR or CRR$_{7.5}$")
    ax.set_title(f"Liquefaction triggering chart ($M_w = {Mw}$)")
    ax.grid(alpha=0.3)
    ax.legend(loc="upper left", fontsize=9)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig