ge.site — API reference

SPT

spt_n60

Python

spt_n60(N: Union[float, ndarray], ER: float = 60.0, Cb: float = 1.0, Cs: float = 1.0, Cr: float = 1.0) -> Union[float, np.ndarray]

Compute the energy-corrected SPT blow count N₆₀.

.. math::

Text Only

N_{60} = N \times \frac{ER}{60} \times C_b \times C_s \times C_r
\qquad \text{[Skempton, 1986]}

PARAMETER	DESCRIPTION
N	Field blow count (blows per 300 mm). TYPE: float or array_like
ER	Energy ratio of the hammer (%). Safety hammer ≈ 60, donut hammer ≈ 45, automatic ≈ 80–95. TYPE: float DEFAULT: 60.0
Cb	Borehole diameter correction (1.0 for 65–115 mm, 1.05 for 150 mm, 1.15 for 200 mm). TYPE: float DEFAULT: 1.0
Cs	Sampler correction (1.0 for standard, 1.1–1.3 without liner). TYPE: float DEFAULT: 1.0
Cr	Rod length correction (0.75 for 3–4 m, 0.85 for 4–6 m, 0.95 for 6–10 m, 1.0 for 10–30 m). TYPE: float DEFAULT: 1.0

RETURNS	DESCRIPTION
float or ndarray	Energy-corrected blow count N₆₀.

References

Skempton (1986); Das (2021), Table 3.8.

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_n60
>>> spt_n60(30, ER=72, Cb=1.0, Cs=1.0, Cr=0.95)
34.2

Source code in geoeq/site/spt.py

def spt_n60(
    N: Union[float, np.ndarray],
    ER: float = 60.0,
    Cb: float = 1.0,
    Cs: float = 1.0,
    Cr: float = 1.0,
) -> Union[float, np.ndarray]:
    r"""
    Compute the energy-corrected SPT blow count N₆₀.

    .. math::

        N_{60} = N \times \frac{ER}{60} \times C_b \times C_s \times C_r
        \qquad \text{[Skempton, 1986]}

    Parameters
    ----------
    N : float or array_like
        Field blow count (blows per 300 mm).
    ER : float, default 60.0
        Energy ratio of the hammer (%). Safety hammer ≈ 60,
        donut hammer ≈ 45, automatic ≈ 80–95.
    Cb : float, default 1.0
        Borehole diameter correction (1.0 for 65–115 mm, 1.05 for 150 mm,
        1.15 for 200 mm).
    Cs : float, default 1.0
        Sampler correction (1.0 for standard, 1.1–1.3 without liner).
    Cr : float, default 1.0
        Rod length correction (0.75 for 3–4 m, 0.85 for 4–6 m,
        0.95 for 6–10 m, 1.0 for 10–30 m).

    Returns
    -------
    float or ndarray
        Energy-corrected blow count N₆₀.

    References
    ----------
    Skempton (1986); Das (2021), Table 3.8.

    Examples
    --------
    >>> from geoeq.site.spt import spt_n60
    >>> spt_n60(30, ER=72, Cb=1.0, Cs=1.0, Cr=0.95)
    34.2
    """
    N_arr = np.asarray(N, dtype=float)
    check_non_negative(N_arr, "N")
    check_positive(ER, "ER")

    result = N_arr * (ER / 60.0) * Cb * Cs * Cr
    return _scalar_or_array(result, N)

spt_n160

Python

spt_n160(N60: Union[float, ndarray], sigma_v: Union[float, ndarray], method: str = 'liao_whitman', pa: float = 100.0) -> Union[float, np.ndarray]

Compute the overburden-corrected SPT blow count (N₁)₆₀.

.. math::

Text Only

(N_1)_{60} = C_N \times N_{60}

PARAMETER	DESCRIPTION
N60	Energy-corrected blow count. TYPE: float or array_like
sigma_v	Effective vertical overburden stress (kPa). TYPE: float or array_like
method	Correction method: 'liao_whitman' : :math:C_N = \sqrt{p_a / \sigma'_v} (Liao & Whitman, 1986). 'skempton' : :math:C_N = 2 / (1 + \sigma'_v / p_a) (Skempton, 1986). 'peck' : :math:C_N = 0.77 \log_{10}(20 \, p_a / \sigma'_v) (Peck et al., 1974). TYPE: str DEFAULT: ``'liao_whitman'``
pa	Atmospheric pressure (kPa). TYPE: float DEFAULT: 100.0

RETURNS	DESCRIPTION
float or ndarray	Overburden-corrected blow count (N₁)₆₀. Capped so CN ≤ 2.0.

References

Liao & Whitman (1986); Skempton (1986); Peck et al. (1974).

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_n160
>>> round(spt_n160(20, sigma_v=100), 1)
20.0
>>> round(spt_n160(20, sigma_v=50), 1)
28.3

Source code in geoeq/site/spt.py

def spt_n160(
    N60: Union[float, np.ndarray],
    sigma_v: Union[float, np.ndarray],
    method: str = "liao_whitman",
    pa: float = 100.0,
) -> Union[float, np.ndarray]:
    r"""
    Compute the overburden-corrected SPT blow count (N₁)₆₀.

    .. math::

        (N_1)_{60} = C_N \times N_{60}

    Parameters
    ----------
    N60 : float or array_like
        Energy-corrected blow count.
    sigma_v : float or array_like
        Effective vertical overburden stress (kPa).
    method : str, default ``'liao_whitman'``
        Correction method:

        - ``'liao_whitman'`` : :math:`C_N = \sqrt{p_a / \sigma'_v}`
          (Liao & Whitman, 1986).
        - ``'skempton'`` : :math:`C_N = 2 / (1 + \sigma'_v / p_a)`
          (Skempton, 1986).
        - ``'peck'`` : :math:`C_N = 0.77 \log_{10}(20 \, p_a / \sigma'_v)`
          (Peck et al., 1974).
    pa : float, default 100.0
        Atmospheric pressure (kPa).

    Returns
    -------
    float or ndarray
        Overburden-corrected blow count (N₁)₆₀. Capped so CN ≤ 2.0.

    References
    ----------
    Liao & Whitman (1986); Skempton (1986); Peck et al. (1974).

    Examples
    --------
    >>> from geoeq.site.spt import spt_n160
    >>> round(spt_n160(20, sigma_v=100), 1)
    20.0
    >>> round(spt_n160(20, sigma_v=50), 1)
    28.3
    """
    n60 = np.asarray(N60, dtype=float)
    sv = np.asarray(sigma_v, dtype=float)
    check_non_negative(n60, "N60")
    check_positive(sv, "sigma_v")

    method = method.lower().replace("-", "_").replace(" ", "_")

    if method == "liao_whitman":
        CN = np.sqrt(pa / sv)
    elif method == "skempton":
        CN = 2.0 / (1.0 + sv / pa)
    elif method == "peck":
        CN = 0.77 * np.log10(20.0 * pa / sv)
    else:
        raise ValueError(
            f"Unknown method '{method}'. Choose from: liao_whitman, skempton, peck."
        )

    CN = np.minimum(CN, 2.0)
    result = CN * n60
    return _scalar_or_array(result, N60, sigma_v)

spt_n160cs

Python

spt_n160cs(N160: Union[float, ndarray], FC: Union[float, ndarray]) -> Union[float, np.ndarray]

Fines-content-corrected (N₁)₆₀cs for liquefaction analysis.

.. math::

Text Only

(N_1)_{60cs} = \alpha + \beta \times (N_1)_{60}

where α and β depend on fines content FC:

• FC ≤ 5%: α = 0, β = 1.0
• 5 < FC < 35%: α = exp(1.76 − 190/FC²), β = 0.99 + FC^1.5 / 1000
• FC ≥ 35%: α = 5.0, β = 1.2

PARAMETER	DESCRIPTION
N160	Overburden-corrected blow count (N₁)₆₀. TYPE: float or array_like
FC	Fines content (%). TYPE: float or array_like

RETURNS	DESCRIPTION
float or ndarray	Fines-corrected (N₁)₆₀cs.

References

Youd, T. L. et al. (2001). Liquefaction Resistance of Soils. JGGE, ASCE, 127(10), 817–833.

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_n160cs
>>> round(spt_n160cs(15, FC=3), 1)
15.0
>>> round(spt_n160cs(15, FC=20), 1)
19.6

Source code in geoeq/site/spt.py

def spt_n160cs(
    N160: Union[float, np.ndarray],
    FC: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""
    Fines-content-corrected (N₁)₆₀cs for liquefaction analysis.

    .. math::

        (N_1)_{60cs} = \alpha + \beta \times (N_1)_{60}

    where α and β depend on fines content FC:

    - FC ≤ 5%: α = 0, β = 1.0
    - 5 < FC < 35%: α = exp(1.76 − 190/FC²), β = 0.99 + FC^1.5 / 1000
    - FC ≥ 35%: α = 5.0, β = 1.2

    Parameters
    ----------
    N160 : float or array_like
        Overburden-corrected blow count (N₁)₆₀.
    FC : float or array_like
        Fines content (%).

    Returns
    -------
    float or ndarray
        Fines-corrected (N₁)₆₀cs.

    References
    ----------
    Youd, T. L. et al. (2001). Liquefaction Resistance of Soils. *JGGE*,
    ASCE, 127(10), 817–833.

    Examples
    --------
    >>> from geoeq.site.spt import spt_n160cs
    >>> round(spt_n160cs(15, FC=3), 1)
    15.0
    >>> round(spt_n160cs(15, FC=20), 1)
    19.6
    """
    n = np.asarray(N160, dtype=float)
    fc = np.asarray(FC, dtype=float)
    check_non_negative(n, "N160")
    check_non_negative(fc, "FC")

    alpha = np.where(fc <= 5, 0.0,
                     np.where(fc >= 35, 5.0,
                              np.exp(1.76 - 190.0 / fc ** 2)))
    beta = np.where(fc <= 5, 1.0,
                    np.where(fc >= 35, 1.2,
                             0.99 + fc ** 1.5 / 1000.0))

    result = alpha + beta * n
    return _scalar_or_array(result, N160, FC)

spt_friction_angle

Python

spt_friction_angle(N: Union[float, ndarray], sigma_v: Optional[Union[float, ndarray]] = None, method: str = 'hatanaka', pa: float = 100.0) -> Union[float, np.ndarray]

Estimate friction angle φ' from SPT blow count.

PARAMETER	DESCRIPTION
N	Blow count. Interpretation depends on method: - 'hatanaka': N should be (N₁)₆₀. - 'kulhawy': N should be (N₁)₆₀; σ'_v is ignored. - 'peck': N should be field N (uncorrected). TYPE: float or array_like
sigma_v	Effective overburden stress (kPa). Not used by all methods. TYPE: float or array_like DEFAULT: None
method	'hatanaka' : :math:\phi' = \sqrt{20\,(N_1)_{60}} + 20 (Hatanaka & Uchida, 1996). 'kulhawy' : :math:\phi' = \tan^{-1}\!\bigl[\frac{(N_1)_{60}}{12.2+20.3\,(\sigma'_v/p_a)}\bigr]^{0.34} (Kulhawy & Mayne, 1990). 'peck' : Peck, Hanson & Thornburn (1974) empirical table. TYPE: str DEFAULT: ``'hatanaka'``
pa	TYPE: float DEFAULT: 100.0

RETURNS	DESCRIPTION
float or ndarray	Friction angle φ' (degrees).

References

Hatanaka & Uchida (1996); Kulhawy & Mayne (1990); Peck et al. (1974).

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_friction_angle
>>> round(spt_friction_angle(20, method='hatanaka'), 1)
40.0

Source code in geoeq/site/spt.py

def spt_friction_angle(
    N: Union[float, np.ndarray],
    sigma_v: Optional[Union[float, np.ndarray]] = None,
    method: str = "hatanaka",
    pa: float = 100.0,
) -> Union[float, np.ndarray]:
    r"""
    Estimate friction angle φ' from SPT blow count.

    Parameters
    ----------
    N : float or array_like
        Blow count. Interpretation depends on *method*:
        - ``'hatanaka'``: N should be (N₁)₆₀.
        - ``'kulhawy'``: N should be (N₁)₆₀; σ'_v is ignored.
        - ``'peck'``: N should be field N (uncorrected).
    sigma_v : float or array_like, optional
        Effective overburden stress (kPa). Not used by all methods.
    method : str, default ``'hatanaka'``
        - ``'hatanaka'`` : :math:`\phi' = \sqrt{20\,(N_1)_{60}} + 20`
          (Hatanaka & Uchida, 1996).
        - ``'kulhawy'`` : :math:`\phi' = \tan^{-1}\!\bigl[\frac{(N_1)_{60}}{12.2+20.3\,(\sigma'_v/p_a)}\bigr]^{0.34}`
          (Kulhawy & Mayne, 1990).
        - ``'peck'`` : Peck, Hanson & Thornburn (1974) empirical table.
    pa : float, default 100.0

    Returns
    -------
    float or ndarray
        Friction angle φ' (degrees).

    References
    ----------
    Hatanaka & Uchida (1996); Kulhawy & Mayne (1990); Peck et al. (1974).

    Examples
    --------
    >>> from geoeq.site.spt import spt_friction_angle
    >>> round(spt_friction_angle(20, method='hatanaka'), 1)
    40.0
    """
    N_arr = np.asarray(N, dtype=float)
    check_non_negative(N_arr, "N")
    method = method.lower()

    if method == "hatanaka":
        phi = np.sqrt(20.0 * N_arr) + 20.0
    elif method == "kulhawy":
        if sigma_v is None:
            sigma_v = pa
        sv = np.asarray(sigma_v, dtype=float)
        phi = np.degrees(np.arctan((N_arr / (12.2 + 20.3 * (sv / pa))) ** 0.34))
    elif method == "peck":
        # Peck, Hanson & Thornburn (1974) empirical approximation
        phi = 26.0 + 0.3 * N_arr - 0.00054 * N_arr ** 2
        phi = np.clip(phi, 26.0, 45.0)
    else:
        raise ValueError(
            f"Unknown method '{method}'. Choose: hatanaka, kulhawy, peck."
        )

    return _scalar_or_array(phi, N)

spt_su

Python

spt_su(N60: Union[float, ndarray], method: str = 'stroud', PI: Optional[float] = None) -> Union[float, np.ndarray]

Estimate undrained shear strength Su from SPT.

PARAMETER	DESCRIPTION
N60	Energy-corrected blow count N₆₀. TYPE: float or array_like
method	'stroud' : :math:S_u = f_1 \times N_{60} (Stroud, 1974). f₁ depends on PI: PI < 20 → f₁ = 4.5; 20–30 → 5.0; 30–40 → 5.5; 40 → 6.0. 'hara' : :math:S_u = 29 \times N_{60}^{0.72} (Hara et al., 1974). TYPE: str DEFAULT: ``'stroud'``
PI	Plasticity index (%). Required for 'stroud'. TYPE: float DEFAULT: None

RETURNS	DESCRIPTION
float or ndarray	Undrained shear strength Su (kPa).

References

Stroud (1974); Hara et al. (1974).

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_su
>>> round(spt_su(10, method='hara'), 1)
152.4

Source code in geoeq/site/spt.py

def spt_su(
    N60: Union[float, np.ndarray],
    method: str = "stroud",
    PI: Optional[float] = None,
) -> Union[float, np.ndarray]:
    r"""
    Estimate undrained shear strength Su from SPT.

    Parameters
    ----------
    N60 : float or array_like
        Energy-corrected blow count N₆₀.
    method : str, default ``'stroud'``
        - ``'stroud'`` : :math:`S_u = f_1 \times N_{60}` (Stroud, 1974).
          f₁ depends on PI: PI < 20 → f₁ = 4.5; 20–30 → 5.0; 30–40 → 5.5;
          > 40 → 6.0.
        - ``'hara'`` : :math:`S_u = 29 \times N_{60}^{0.72}` (Hara et al., 1974).
    PI : float, optional
        Plasticity index (%). Required for ``'stroud'``.

    Returns
    -------
    float or ndarray
        Undrained shear strength Su (kPa).

    References
    ----------
    Stroud (1974); Hara et al. (1974).

    Examples
    --------
    >>> from geoeq.site.spt import spt_su
    >>> round(spt_su(10, method='hara'), 1)
    152.4
    """
    n60 = np.asarray(N60, dtype=float)
    check_non_negative(n60, "N60")
    method = method.lower()

    if method == "stroud":
        if PI is None:
            f1 = 5.0
        elif PI < 20:
            f1 = 4.5
        elif PI < 30:
            f1 = 5.0
        elif PI < 40:
            f1 = 5.5
        else:
            f1 = 6.0
        Su = f1 * n60
    elif method == "hara":
        Su = 29.0 * n60 ** 0.72
    else:
        raise ValueError(f"Unknown method '{method}'. Choose: stroud, hara.")

    return _scalar_or_array(Su, N60)

spt_dr

Python

spt_dr(N160: Union[float, ndarray], method: str = 'meyerhof') -> Union[float, np.ndarray]

Estimate relative density Dr from SPT.

PARAMETER	DESCRIPTION
N160	Overburden-corrected (N₁)₆₀. TYPE: float or array_like
method	'meyerhof' : :math:D_r (\%) = \sqrt{(N_1)_{60} / 0.21} × 100 (adapted from Meyerhof, 1957 — coeff for clean normally-consolidated sands). 'skempton' : :math:D_r (\%) = \sqrt{(N_1)_{60} / 0.28} × 100 (Skempton, 1986 — fine sands). 'kulhawy' : :math:D_r (\%) = \sqrt{(N_1)_{60} / C_p} × 100, where Cp = 60 + 25 log D₅₀ (Kulhawy & Mayne, 1990). Uses Cp = 60 (default for medium sand D₅₀ ≈ 1 mm). TYPE: str DEFAULT: ``'meyerhof'``

RETURNS	DESCRIPTION
float or ndarray	Relative density Dr (%).

References

Meyerhof (1957); Skempton (1986); Kulhawy & Mayne (1990).

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_dr
>>> round(spt_dr(20, method='meyerhof'), 0)
98.0

Source code in geoeq/site/spt.py

def spt_dr(
    N160: Union[float, np.ndarray],
    method: str = "meyerhof",
) -> Union[float, np.ndarray]:
    r"""
    Estimate relative density Dr from SPT.

    Parameters
    ----------
    N160 : float or array_like
        Overburden-corrected (N₁)₆₀.
    method : str, default ``'meyerhof'``
        - ``'meyerhof'`` : :math:`D_r (\%) = \sqrt{(N_1)_{60} / 0.21}`
          × 100 (adapted from Meyerhof, 1957 — coeff for clean
          normally-consolidated sands).
        - ``'skempton'`` : :math:`D_r (\%) = \sqrt{(N_1)_{60} / 0.28}`
          × 100 (Skempton, 1986 — fine sands).
        - ``'kulhawy'`` : :math:`D_r (\%) = \sqrt{(N_1)_{60} / C_p}`
          × 100, where Cp = 60 + 25 log D₅₀ (Kulhawy & Mayne, 1990).
          Uses Cp = 60 (default for medium sand D₅₀ ≈ 1 mm).

    Returns
    -------
    float or ndarray
        Relative density Dr (%).

    References
    ----------
    Meyerhof (1957); Skempton (1986); Kulhawy & Mayne (1990).

    Examples
    --------
    >>> from geoeq.site.spt import spt_dr
    >>> round(spt_dr(20, method='meyerhof'), 0)
    98.0
    """
    n = np.asarray(N160, dtype=float)
    check_non_negative(n, "N160")
    method = method.lower()

    if method == "meyerhof":
        # Dr² = N1,60 / 0.21 → Dr = √(N/0.21) as fraction, ×100 for %
        # More standard: Dr(%) = √(N1,60 / 41) × 100 from Meyerhof tables
        Dr = np.sqrt(n / 41.0) * 100.0
    elif method == "skempton":
        Dr = np.sqrt(n / 55.0) * 100.0
    elif method == "kulhawy":
        Cp = 60.0
        Dr = np.sqrt(n / Cp) * 100.0
    else:
        raise ValueError(f"Unknown method '{method}'. Choose: meyerhof, skempton, kulhawy.")

    return _scalar_or_array(np.clip(Dr, 0, 100), N160)

spt_modulus

Python

spt_modulus(N: Union[float, ndarray], soil_type: str = 'sand') -> Union[float, np.ndarray]

Estimate elastic (Young's) modulus Es from SPT.

PARAMETER	DESCRIPTION
N	SPT blow count (N or N₆₀ depending on source — typically N₆₀). TYPE: float or array_like
soil_type	'sand' : Es = 500 (N + 15) kPa (Bowles, 1996, for NC sand). 'sand_oc' : Es = 750 (N + 24) kPa (Bowles, 1996, OC sand). 'gravel' : Es = 1200 (N + 6) kPa (Bowles, 1996). 'clay_soft' : Es = 300 (N + 6) kPa (Bowles, 1996). 'clay_stiff' : Es = 500 (N + 15) kPa (Bowles, 1996). 'silt' : Es = 300 (N + 6) kPa (Bowles, 1996). TYPE: str DEFAULT: ``'sand'``

RETURNS	DESCRIPTION
float or ndarray	Elastic modulus Es (kPa).

References

Bowles (1996), Table 5.6.

Examples:

Python Console Session

>>> from geoeq.site.spt import spt_modulus
>>> spt_modulus(20, soil_type='sand')
17500.0

Source code in geoeq/site/spt.py

def spt_modulus(
    N: Union[float, np.ndarray],
    soil_type: str = "sand",
) -> Union[float, np.ndarray]:
    r"""
    Estimate elastic (Young's) modulus Es from SPT.

    Parameters
    ----------
    N : float or array_like
        SPT blow count (N or N₆₀ depending on source — typically N₆₀).
    soil_type : str, default ``'sand'``
        - ``'sand'`` : Es = 500 (N + 15) kPa (Bowles, 1996, for NC sand).
        - ``'sand_oc'`` : Es = 750 (N + 24) kPa (Bowles, 1996, OC sand).
        - ``'gravel'`` : Es = 1200 (N + 6) kPa (Bowles, 1996).
        - ``'clay_soft'`` : Es = 300 (N + 6) kPa (Bowles, 1996).
        - ``'clay_stiff'`` : Es = 500 (N + 15) kPa (Bowles, 1996).
        - ``'silt'`` : Es = 300 (N + 6) kPa (Bowles, 1996).

    Returns
    -------
    float or ndarray
        Elastic modulus Es (kPa).

    References
    ----------
    Bowles (1996), Table 5.6.

    Examples
    --------
    >>> from geoeq.site.spt import spt_modulus
    >>> spt_modulus(20, soil_type='sand')
    17500.0
    """
    n = np.asarray(N, dtype=float)
    check_non_negative(n, "N")
    st = soil_type.lower().replace(" ", "_")

    lookup = {
        "sand": (500.0, 15.0),
        "sand_nc": (500.0, 15.0),
        "sand_oc": (750.0, 24.0),
        "gravel": (1200.0, 6.0),
        "clay_soft": (300.0, 6.0),
        "clay_stiff": (500.0, 15.0),
        "silt": (300.0, 6.0),
    }

    if st not in lookup:
        raise ValueError(
            f"Unknown soil_type '{soil_type}'. Choose from: "
            + ", ".join(lookup.keys())
        )

    a, b = lookup[st]
    Es = a * (n + b)
    return _scalar_or_array(Es, N)

CPT

cpt_normalize

Python

cpt_normalize(qt: Union[float, ndarray], fs: Union[float, ndarray], sigma_v: Union[float, ndarray], sigma_v_eff: Union[float, ndarray], u2: Union[float, ndarray] = 0.0, u0: Union[float, ndarray] = 0.0) -> Dict[str, Union[float, np.ndarray]]

Compute normalized CPT parameters.

.. math::

Text Only

Q_t = \frac{q_t - \sigma_{v0}}{\sigma'_{v0}} \qquad
F_r = \frac{f_s}{q_t - \sigma_{v0}} \times 100 \qquad
B_q = \frac{u_2 - u_0}{q_t - \sigma_{v0}}

PARAMETER	DESCRIPTION
qt	Corrected cone resistance (kPa). TYPE: float or array_like
fs	Sleeve friction (kPa). TYPE: float or array_like
sigma_v	Total vertical stress σ_v0 (kPa). TYPE: float or array_like
sigma_v_eff	Effective vertical stress σ'_v0 (kPa). TYPE: float or array_like
u2	Pore pressure at shoulder (kPa). TYPE: float or array_like DEFAULT: 0.0
u0	Equilibrium pore pressure (kPa). TYPE: float or array_like DEFAULT: 0.0

RETURNS	DESCRIPTION
dict	'Qt', 'Fr', 'Bq'.

References

Robertson (1990, 2009).

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_normalize
>>> res = cpt_normalize(qt=5000, fs=50, sigma_v=100, sigma_v_eff=60)
>>> round(res['Qt'], 1)
81.7

Source code in geoeq/site/cpt.py

def cpt_normalize(
    qt: Union[float, np.ndarray],
    fs: Union[float, np.ndarray],
    sigma_v: Union[float, np.ndarray],
    sigma_v_eff: Union[float, np.ndarray],
    u2: Union[float, np.ndarray] = 0.0,
    u0: Union[float, np.ndarray] = 0.0,
) -> Dict[str, Union[float, np.ndarray]]:
    r"""
    Compute normalized CPT parameters.

    .. math::

        Q_t = \frac{q_t - \sigma_{v0}}{\sigma'_{v0}} \qquad
        F_r = \frac{f_s}{q_t - \sigma_{v0}} \times 100 \qquad
        B_q = \frac{u_2 - u_0}{q_t - \sigma_{v0}}

    Parameters
    ----------
    qt : float or array_like
        Corrected cone resistance (kPa).
    fs : float or array_like
        Sleeve friction (kPa).
    sigma_v : float or array_like
        Total vertical stress σ_v0 (kPa).
    sigma_v_eff : float or array_like
        Effective vertical stress σ'_v0 (kPa).
    u2 : float or array_like, default 0.0
        Pore pressure at shoulder (kPa).
    u0 : float or array_like, default 0.0
        Equilibrium pore pressure (kPa).

    Returns
    -------
    dict
        ``'Qt'``, ``'Fr'``, ``'Bq'``.

    References
    ----------
    Robertson (1990, 2009).

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_normalize
    >>> res = cpt_normalize(qt=5000, fs=50, sigma_v=100, sigma_v_eff=60)
    >>> round(res['Qt'], 1)
    81.7
    """
    qt_a = np.asarray(qt, dtype=float)
    fs_a = np.asarray(fs, dtype=float)
    sv = np.asarray(sigma_v, dtype=float)
    sve = np.asarray(sigma_v_eff, dtype=float)
    u2_a = np.asarray(u2, dtype=float)
    u0_a = np.asarray(u0, dtype=float)

    net = qt_a - sv
    net = np.maximum(net, 1.0)  # avoid division by zero

    Qt = net / sve
    Fr = (fs_a / net) * 100.0
    Bq = (u2_a - u0_a) / net

    return {
        "Qt": _scalar_or_array(Qt, qt, fs),
        "Fr": _scalar_or_array(Fr, qt, fs),
        "Bq": _scalar_or_array(Bq, qt, fs),
    }

cpt_ic

Python

cpt_ic(Qt: Union[float, ndarray], Fr: Union[float, ndarray]) -> Union[float, np.ndarray]

Compute the Soil Behaviour Type Index Ic.

.. math::

Text Only

I_c = \sqrt{(3.47 - \log Q_t)^2 + (\log F_r + 1.22)^2}

PARAMETER	DESCRIPTION
Qt	Normalized cone resistance. TYPE: float or array_like
Fr	Normalized friction ratio (%). TYPE: float or array_like

RETURNS	DESCRIPTION
float or ndarray	Soil behaviour type index Ic.

References

Robertson (2009), Eq. 2.

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_ic
>>> round(cpt_ic(Qt=80, Fr=1.0), 2)
1.62

Source code in geoeq/site/cpt.py

def cpt_ic(
    Qt: Union[float, np.ndarray],
    Fr: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""
    Compute the Soil Behaviour Type Index Ic.

    .. math::

        I_c = \sqrt{(3.47 - \log Q_t)^2 + (\log F_r + 1.22)^2}

    Parameters
    ----------
    Qt : float or array_like
        Normalized cone resistance.
    Fr : float or array_like
        Normalized friction ratio (%).

    Returns
    -------
    float or ndarray
        Soil behaviour type index Ic.

    References
    ----------
    Robertson (2009), Eq. 2.

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_ic
    >>> round(cpt_ic(Qt=80, Fr=1.0), 2)
    1.62
    """
    Qt_a = np.asarray(Qt, dtype=float)
    Fr_a = np.asarray(Fr, dtype=float)
    Fr_a = np.maximum(Fr_a, 0.1)  # avoid log(0)

    Ic = np.sqrt((3.47 - np.log10(Qt_a)) ** 2 + (np.log10(Fr_a) + 1.22) ** 2)
    return _scalar_or_array(Ic, Qt, Fr)

cpt_sbt

Python

cpt_sbt(Qt: Union[float, ndarray], Fr: Union[float, ndarray]) -> Dict[str, Union[int, str, float, np.ndarray]]

Robertson (2009) Soil Behaviour Type classification from Ic.

Zone boundaries (Robertson 2009):

• Ic > 3.60 → Zone 2: Organic soils — clay
• 2.95 < Ic ≤ 3.60 → Zone 3: Clay — silty clay
• 2.60 < Ic ≤ 2.95 → Zone 4: Silt mixtures
• 2.05 < Ic ≤ 2.60 → Zone 5: Sand mixtures
• 1.31 < Ic ≤ 2.05 → Zone 6: Sands — clean sand
• Ic ≤ 1.31 → Zone 7: Gravelly sand to dense sand

PARAMETER	DESCRIPTION
Qt	Normalized cone resistance. TYPE: float or array_like
Fr	Normalized friction ratio (%). TYPE: float or array_like

RETURNS	DESCRIPTION
dict	'Ic' : float or ndarray — Soil behaviour type index. 'zone' : int or list — Robertson zone number (2–7). 'description' : str or list — soil type description.

References

Robertson (2009), Table 1.

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_sbt
>>> res = cpt_sbt(Qt=80, Fr=1.0)
>>> res['zone']
6

Source code in geoeq/site/cpt.py

def cpt_sbt(
    Qt: Union[float, np.ndarray],
    Fr: Union[float, np.ndarray],
) -> Dict[str, Union[int, str, float, np.ndarray]]:
    r"""
    Robertson (2009) Soil Behaviour Type classification from Ic.

    Zone boundaries (Robertson 2009):

    - Ic > 3.60 → Zone 2: Organic soils — clay
    - 2.95 < Ic ≤ 3.60 → Zone 3: Clay — silty clay
    - 2.60 < Ic ≤ 2.95 → Zone 4: Silt mixtures
    - 2.05 < Ic ≤ 2.60 → Zone 5: Sand mixtures
    - 1.31 < Ic ≤ 2.05 → Zone 6: Sands — clean sand
    - Ic ≤ 1.31 → Zone 7: Gravelly sand to dense sand

    Parameters
    ----------
    Qt : float or array_like
        Normalized cone resistance.
    Fr : float or array_like
        Normalized friction ratio (%).

    Returns
    -------
    dict
        ``'Ic'`` : float or ndarray — Soil behaviour type index.
        ``'zone'`` : int or list — Robertson zone number (2–7).
        ``'description'`` : str or list — soil type description.

    References
    ----------
    Robertson (2009), Table 1.

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_sbt
    >>> res = cpt_sbt(Qt=80, Fr=1.0)
    >>> res['zone']
    6
    """
    Ic_val = cpt_ic(Qt, Fr)

    _ZONES = [
        (3.60, 2, "Organic soils — clay"),
        (2.95, 3, "Clay — silty clay to clay"),
        (2.60, 4, "Silt mixtures — clayey silt to silty clay"),
        (2.05, 5, "Sand mixtures — silty sand to sandy silt"),
        (1.31, 6, "Sands — clean sand to silty sand"),
        (0.00, 7, "Gravelly sand to dense sand"),
    ]

    def _classify_one(ic):
        for threshold, zone, desc in _ZONES:
            if ic > threshold:
                return zone, desc
        return 7, "Gravelly sand to dense sand"

    if np.ndim(Ic_val) == 0:
        zone, desc = _classify_one(float(Ic_val))
        return {"Ic": float(Ic_val), "zone": zone, "description": desc}
    else:
        ic_arr = np.atleast_1d(Ic_val)
        zones = []
        descs = []
        for ic in ic_arr:
            z, d = _classify_one(float(ic))
            zones.append(z)
            descs.append(d)
        return {"Ic": Ic_val, "zone": zones, "description": descs}

cpt_friction_angle

Python

cpt_friction_angle(qt: Union[float, ndarray], sigma_v_eff: Union[float, ndarray]) -> Union[float, np.ndarray]

Estimate drained friction angle from CPT.

.. math::

Text Only

\phi' = \arctan\!\bigl[0.1 + 0.38 \log(q_t / \sigma'_{v0})\bigr]
\qquad \text{[Robertson \& Campanella, 1983]}

PARAMETER	DESCRIPTION
qt	Corrected cone resistance (kPa). TYPE: float or array_like
sigma_v_eff	Effective vertical stress (kPa). TYPE: float or array_like

RETURNS	DESCRIPTION
float or ndarray	Friction angle φ' (degrees).

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_friction_angle
>>> round(cpt_friction_angle(5000, 60), 1)
35.1

Source code in geoeq/site/cpt.py

def cpt_friction_angle(
    qt: Union[float, np.ndarray],
    sigma_v_eff: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""
    Estimate drained friction angle from CPT.

    .. math::

        \phi' = \arctan\!\bigl[0.1 + 0.38 \log(q_t / \sigma'_{v0})\bigr]
        \qquad \text{[Robertson \& Campanella, 1983]}

    Parameters
    ----------
    qt : float or array_like
        Corrected cone resistance (kPa).
    sigma_v_eff : float or array_like
        Effective vertical stress (kPa).

    Returns
    -------
    float or ndarray
        Friction angle φ' (degrees).

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_friction_angle
    >>> round(cpt_friction_angle(5000, 60), 1)
    35.1
    """
    qt_a = np.asarray(qt, dtype=float)
    sve = np.asarray(sigma_v_eff, dtype=float)
    check_positive(qt_a, "qt")
    check_positive(sve, "sigma_v_eff")

    ratio = qt_a / sve
    phi = np.degrees(np.arctan(0.1 + 0.38 * np.log10(ratio)))
    return _scalar_or_array(phi, qt, sigma_v_eff)

cpt_su

Python

cpt_su(qt: Union[float, ndarray], sigma_v: Union[float, ndarray], Nkt: float = 14.0) -> Union[float, np.ndarray]

Estimate undrained shear strength from CPT.

.. math::

Text Only

S_u = \frac{q_t - \sigma_{v0}}{N_{kt}}

PARAMETER	DESCRIPTION
qt	Corrected cone resistance (kPa). TYPE: float or array_like
sigma_v	Total vertical stress (kPa). TYPE: float or array_like
Nkt	Cone factor (typically 10–20; 14 is common for medium- sensitivity clays). TYPE: float DEFAULT: 14.0

RETURNS	DESCRIPTION
float or ndarray	Undrained shear strength Su (kPa).

References

Lunne et al. (1997), Ch. 6.

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_su
>>> cpt_su(qt=1000, sigma_v=100, Nkt=14)
64.28571428571429

Source code in geoeq/site/cpt.py

def cpt_su(
    qt: Union[float, np.ndarray],
    sigma_v: Union[float, np.ndarray],
    Nkt: float = 14.0,
) -> Union[float, np.ndarray]:
    r"""
    Estimate undrained shear strength from CPT.

    .. math::

        S_u = \frac{q_t - \sigma_{v0}}{N_{kt}}

    Parameters
    ----------
    qt : float or array_like
        Corrected cone resistance (kPa).
    sigma_v : float or array_like
        Total vertical stress (kPa).
    Nkt : float, default 14.0
        Cone factor (typically 10–20; 14 is common for medium-
        sensitivity clays).

    Returns
    -------
    float or ndarray
        Undrained shear strength Su (kPa).

    References
    ----------
    Lunne et al. (1997), Ch. 6.

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_su
    >>> cpt_su(qt=1000, sigma_v=100, Nkt=14)
    64.28571428571429
    """
    qt_a = np.asarray(qt, dtype=float)
    sv = np.asarray(sigma_v, dtype=float)
    check_positive(Nkt, "Nkt")

    Su = (qt_a - sv) / Nkt
    return _scalar_or_array(np.maximum(Su, 0.0), qt, sigma_v)

cpt_dr

Python

cpt_dr(qc: Union[float, ndarray], sigma_v_eff: Union[float, ndarray], C0: float = 157.0, C1: float = 0.55, C2: float = 2.41) -> Union[float, np.ndarray]

Estimate relative density from CPT for normally consolidated uncemented quartz sands.

.. math::

Text Only

D_r = \frac{1}{C_2} \ln\!\left(\frac{q_c}{C_0 \,(\sigma'_{v0})^{C_1}}\right)
\qquad \text{[Baldi et al., 1986]}

Default coefficients C₀ = 157, C₁ = 0.55, C₂ = 2.41 are for Ticino sand (Baldi et al., 1986).

PARAMETER	DESCRIPTION
qc	Cone resistance (kPa). TYPE: float or array_like
sigma_v_eff	Effective vertical stress (kPa). TYPE: float or array_like
C0	Calibration constants. TYPE: float DEFAULT: 157.0
C1	Calibration constants. TYPE: float DEFAULT: 157.0
C2	Calibration constants. TYPE: float DEFAULT: 157.0

RETURNS	DESCRIPTION
float or ndarray	Relative density Dr (fraction, 0–1).

References

Baldi et al. (1986).

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_dr
>>> round(cpt_dr(10000, 100) * 100, 0)
55.0

Source code in geoeq/site/cpt.py

def cpt_dr(
    qc: Union[float, np.ndarray],
    sigma_v_eff: Union[float, np.ndarray],
    C0: float = 157.0,
    C1: float = 0.55,
    C2: float = 2.41,
) -> Union[float, np.ndarray]:
    r"""
    Estimate relative density from CPT for normally consolidated
    uncemented quartz sands.

    .. math::

        D_r = \frac{1}{C_2} \ln\!\left(\frac{q_c}{C_0 \,(\sigma'_{v0})^{C_1}}\right)
        \qquad \text{[Baldi et al., 1986]}

    Default coefficients C₀ = 157, C₁ = 0.55, C₂ = 2.41 are for
    Ticino sand (Baldi et al., 1986).

    Parameters
    ----------
    qc : float or array_like
        Cone resistance (kPa).
    sigma_v_eff : float or array_like
        Effective vertical stress (kPa).
    C0, C1, C2 : float
        Calibration constants.

    Returns
    -------
    float or ndarray
        Relative density Dr (fraction, 0–1).

    References
    ----------
    Baldi et al. (1986).

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_dr
    >>> round(cpt_dr(10000, 100) * 100, 0)
    55.0
    """
    qc_a = np.asarray(qc, dtype=float)
    sve = np.asarray(sigma_v_eff, dtype=float)
    check_positive(qc_a, "qc")
    check_positive(sve, "sigma_v_eff")

    Dr = (1.0 / C2) * np.log(qc_a / (C0 * sve ** C1))
    Dr = np.clip(Dr, 0.0, 1.0)
    return _scalar_or_array(Dr, qc, sigma_v_eff)

cpt_modulus

Python

cpt_modulus(qt: Union[float, ndarray], sigma_v: Union[float, ndarray], Ic: Optional[Union[float, ndarray]] = None, alpha: Optional[float] = None) -> Union[float, np.ndarray]

Estimate constrained modulus M from CPT.

.. math::

Text Only

M = \alpha_M \times (q_t - \sigma_{v0})

where α_M depends on Ic (Robertson, 2009):

• Ic > 2.2 (clays): α_M = Qt when Qt < 14, else α_M = 14
• Ic ≤ 2.2 (sands): α_M = 0.0188 × 10^(0.55 Ic + 1.68)

PARAMETER	DESCRIPTION
qt	Corrected cone resistance (kPa). TYPE: float or array_like
sigma_v	Total vertical stress (kPa). TYPE: float or array_like
Ic	Soil behaviour type index. If None, alpha must be provided. TYPE: float or array_like DEFAULT: None
alpha	Direct multiplier override. TYPE: float DEFAULT: None

RETURNS	DESCRIPTION
float or ndarray	Constrained modulus M (kPa).

References

Robertson (2009), Eq. 11.

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_modulus
>>> cpt_modulus(qt=5000, sigma_v=100, alpha=5)
24500.0

Source code in geoeq/site/cpt.py

def cpt_modulus(
    qt: Union[float, np.ndarray],
    sigma_v: Union[float, np.ndarray],
    Ic: Optional[Union[float, np.ndarray]] = None,
    alpha: Optional[float] = None,
) -> Union[float, np.ndarray]:
    r"""
    Estimate constrained modulus M from CPT.

    .. math::

        M = \alpha_M \times (q_t - \sigma_{v0})

    where α_M depends on Ic (Robertson, 2009):

    - Ic > 2.2 (clays): α_M = Qt when Qt < 14, else α_M = 14
    - Ic ≤ 2.2 (sands): α_M = 0.0188 × 10^(0.55 Ic + 1.68)

    Parameters
    ----------
    qt : float or array_like
        Corrected cone resistance (kPa).
    sigma_v : float or array_like
        Total vertical stress (kPa).
    Ic : float or array_like, optional
        Soil behaviour type index. If None, alpha must be provided.
    alpha : float, optional
        Direct multiplier override.

    Returns
    -------
    float or ndarray
        Constrained modulus M (kPa).

    References
    ----------
    Robertson (2009), Eq. 11.

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_modulus
    >>> cpt_modulus(qt=5000, sigma_v=100, alpha=5)
    24500.0
    """
    qt_a = np.asarray(qt, dtype=float)
    sv = np.asarray(sigma_v, dtype=float)
    net = qt_a - sv

    if alpha is not None:
        M = alpha * net
    elif Ic is not None:
        Ic_a = np.asarray(Ic, dtype=float)
        alpha_M = np.where(
            Ic_a > 2.2,
            np.minimum(net / np.maximum(sv, 1.0), 14.0),
            0.0188 * 10 ** (0.55 * Ic_a + 1.68),
        )
        M = alpha_M * net
    else:
        raise ValueError("Either Ic or alpha must be provided.")

    return _scalar_or_array(np.maximum(M, 0.0), qt, sigma_v)

cpt_sbt_plot

Python

cpt_sbt_plot(Qt: Union[float, list, ndarray], Fr: Union[float, list, ndarray], labels=None, ax=None, save_as: Optional[str] = None) -> Dict

Plot data on the Robertson (2009) SBT chart (Qt vs Fr).

PARAMETER	DESCRIPTION
Qt	Normalized cone resistance. TYPE: array_like
Fr	Normalized friction ratio (%). TYPE: array_like
labels	Labels for each data point. TYPE: list DEFAULT: None
ax	TYPE: Axes DEFAULT: None
save_as	TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
dict	'ax' : Axes.

Examples:

Python Console Session

>>> from geoeq.site.cpt import cpt_sbt_plot
>>> res = cpt_sbt_plot([80, 20, 5], [1.0, 3.0, 5.0])

Source code in geoeq/site/cpt.py

def cpt_sbt_plot(
    Qt: Union[float, list, np.ndarray],
    Fr: Union[float, list, np.ndarray],
    labels=None,
    ax=None,
    save_as: Optional[str] = None,
) -> Dict:
    r"""
    Plot data on the Robertson (2009) SBT chart (Qt vs Fr).

    Parameters
    ----------
    Qt : array_like
        Normalized cone resistance.
    Fr : array_like
        Normalized friction ratio (%).
    labels : list, optional
        Labels for each data point.
    ax : matplotlib.axes.Axes, optional
    save_as : str, optional

    Returns
    -------
    dict
        ``'ax'`` : Axes.

    Examples
    --------
    >>> from geoeq.site.cpt import cpt_sbt_plot
    >>> res = cpt_sbt_plot([80, 20, 5], [1.0, 3.0, 5.0])
    """
    import matplotlib.pyplot as plt

    Qt_a = np.atleast_1d(np.asarray(Qt, dtype=float))
    Fr_a = np.atleast_1d(np.asarray(Fr, dtype=float))

    if ax is None:
        fig, ax = plt.subplots(figsize=(9, 7))
    else:
        fig = ax.get_figure()

    # Draw Ic contours
    Fr_grid = np.logspace(-1, 1, 300)
    for Ic_val, ls in [(1.31, "--"), (2.05, "--"), (2.60, "--"), (2.95, "--"), (3.60, "--")]:
        Qt_contour = 10 ** (3.47 - np.sqrt(Ic_val ** 2 - (np.log10(Fr_grid) + 1.22) ** 2))
        mask = np.isfinite(Qt_contour) & (Qt_contour > 0)
        ax.plot(Fr_grid[mask], Qt_contour[mask], "k", ls=ls, alpha=0.4, linewidth=0.8)

    # Zone labels
    zone_labels = {
        (0.5, 300): "7\nGravelly\nsand",
        (0.7, 50): "6\nClean sand",
        (1.5, 15): "5\nSand\nmixtures",
        (3.0, 5): "4\nSilt\nmixtures",
        (4.0, 2): "3\nClay",
        (6.0, 1.5): "2\nOrganic",
    }
    for (fr_pos, qt_pos), label in zone_labels.items():
        ax.text(fr_pos, qt_pos, label, fontsize=7, ha="center", va="center",
                color="#666666", alpha=0.7)

    # Data points
    ax.scatter(Fr_a, Qt_a, s=60, c="steelblue", edgecolors="navy", zorder=5)
    if labels is not None:
        for i, lbl in enumerate(labels):
            ax.annotate(str(lbl), (Fr_a[i], Qt_a[i]), fontsize=8,
                        xytext=(5, 5), textcoords="offset points")

    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xlabel("Normalized Friction Ratio Fr (%)", fontweight="bold")
    ax.set_ylabel("Normalized Cone Resistance Qt", fontweight="bold")
    ax.set_title("Robertson (2009) CPT Soil Behaviour Type Chart", fontweight="bold")
    ax.set_xlim(0.1, 10)
    ax.set_ylim(1, 1000)
    ax.grid(True, which="both", alpha=0.2)

    if save_as:
        plt.savefig(save_as, bbox_inches="tight", dpi=300)

    return {"ax": ax}

Field vane shear

vane_su

Python

vane_su(T: Union[float, ndarray], D: float, H: float) -> Union[float, np.ndarray]

Undrained shear strength from field vane torque.

For a rectangular vane with height H and diameter D, assuming uniform shear stress on the cylindrical surface and both end caps:

.. math::

Text Only

S_u = \frac{T}{\pi \left(\frac{D^2 H}{2} + \frac{D^3}{6}\right)}
\qquad \text{[ASTM D2573]}

For the standard vane with H/D = 2:

.. math::

Text Only

S_u = \frac{6\,T}{7\,\pi\,D^3}

PARAMETER	DESCRIPTION
T	Maximum torque (kN·m). TYPE: float or array_like
D	Vane diameter (m). Standard sizes: 0.0508 m (2 in.), 0.0635 m (2.5 in.), 0.0925 m (3.64 in.). TYPE: float
H	Vane height (m). Standard H/D = 2. TYPE: float

RETURNS	DESCRIPTION
float or ndarray	Undrained shear strength Su (kPa).

References

ASTM D2573 (2018); Das (2021), Section 3.15.

Examples:

Python Console Session

>>> from geoeq.site.vane import vane_su
>>> round(vane_su(T=0.012, D=0.065, H=0.130), 1)
20.9

Source code in geoeq/site/vane.py

def vane_su(
    T: Union[float, np.ndarray],
    D: float,
    H: float,
) -> Union[float, np.ndarray]:
    r"""
    Undrained shear strength from field vane torque.

    For a rectangular vane with height H and diameter D, assuming
    uniform shear stress on the cylindrical surface and both end caps:

    .. math::

        S_u = \frac{T}{\pi \left(\frac{D^2 H}{2} + \frac{D^3}{6}\right)}
        \qquad \text{[ASTM D2573]}

    For the standard vane with H/D = 2:

    .. math::

        S_u = \frac{6\,T}{7\,\pi\,D^3}

    Parameters
    ----------
    T : float or array_like
        Maximum torque (kN·m).
    D : float
        Vane diameter (m).  Standard sizes: 0.0508 m (2 in.),
        0.0635 m (2.5 in.), 0.0925 m (3.64 in.).
    H : float
        Vane height (m).  Standard H/D = 2.

    Returns
    -------
    float or ndarray
        Undrained shear strength Su (kPa).

    References
    ----------
    ASTM D2573 (2018); Das (2021), Section 3.15.

    Examples
    --------
    >>> from geoeq.site.vane import vane_su
    >>> round(vane_su(T=0.012, D=0.065, H=0.130), 1)
    20.9
    """
    T_arr = np.asarray(T, dtype=float)
    check_non_negative(T_arr, "T")
    check_positive(D, "D")
    check_positive(H, "H")

    K = np.pi * (D**2 * H / 2.0 + D**3 / 6.0)
    Su = T_arr / K
    return _scalar_or_array(Su, T)

vane_correction

Python

vane_correction(Su_fvt: Union[float, ndarray], PI: Union[float, ndarray]) -> Union[float, np.ndarray]

Bjerrum (1972) correction for field vane undrained shear strength.

The field vane overestimates the mobilised Su on the failure surface. The correction factor μ depends on the soil's plasticity index:

.. math::

Text Only

S_{u,\text{design}} = \mu \times S_{u,\text{FVT}}

Bjerrum (1972):

.. math::

Text Only

\mu = 1.7 - 0.54 \log_{10}(\text{PI})
\qquad \text{for PI in \%}

For very low PI (< 20 %), μ is capped at 1.0.

PARAMETER	DESCRIPTION
Su_fvt	Field vane undrained shear strength (kPa). TYPE: float or array_like
PI	Plasticity index (%). TYPE: float or array_like

RETURNS	DESCRIPTION
float or ndarray	Corrected undrained shear strength Su,design (kPa).

Notes

• Bjerrum's factor typically ranges from 0.5 to 1.0.
• Aas et al. (1986) proposed an alternative correction based on Su/σ'v ratio; use the Bjerrum factor for routine design.

References

Bjerrum (1972); Das (2021), Eq. 3.39.

Examples:

Python Console Session

>>> from geoeq.site.vane import vane_correction
>>> round(vane_correction(Su_fvt=50, PI=40), 1)
38.2

Source code in geoeq/site/vane.py

def vane_correction(
    Su_fvt: Union[float, np.ndarray],
    PI: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    r"""
    Bjerrum (1972) correction for field vane undrained shear strength.

    The field vane overestimates the mobilised Su on the failure
    surface.  The correction factor μ depends on the soil's plasticity
    index:

    .. math::

        S_{u,\text{design}} = \mu \times S_{u,\text{FVT}}

    Bjerrum (1972):

    .. math::

        \mu = 1.7 - 0.54 \log_{10}(\text{PI})
        \qquad \text{for PI in \%}

    For very low PI (< 20 %), μ is capped at 1.0.

    Parameters
    ----------
    Su_fvt : float or array_like
        Field vane undrained shear strength (kPa).
    PI : float or array_like
        Plasticity index (%).

    Returns
    -------
    float or ndarray
        Corrected undrained shear strength Su,design (kPa).

    Notes
    -----
    * Bjerrum's factor typically ranges from 0.5 to 1.0.
    * Aas et al. (1986) proposed an alternative correction based on
      Su/σ'v ratio; use the Bjerrum factor for routine design.

    References
    ----------
    Bjerrum (1972); Das (2021), Eq. 3.39.

    Examples
    --------
    >>> from geoeq.site.vane import vane_correction
    >>> round(vane_correction(Su_fvt=50, PI=40), 1)
    38.2
    """
    Su = np.asarray(Su_fvt, dtype=float)
    pi = np.asarray(PI, dtype=float)
    check_non_negative(Su, "Su_fvt")
    check_positive(pi, "PI")

    mu = 1.7 - 0.54 * np.log10(pi)
    mu = np.clip(mu, 0.5, 1.0)

    result = mu * Su
    return _scalar_or_array(result, Su_fvt, PI)

vane_remolded

Python

vane_remolded(T_peak: float, T_remolded: float, D: float, H: float) -> Dict[str, float]

Undrained shear strength (peak and remolded) and sensitivity.

After the peak torque test, the vane is rotated rapidly through several full revolutions to remold the soil, then retested.

.. math::

Text Only

S_t = \frac{S_{u,\text{peak}}}{S_{u,\text{remolded}}}

PARAMETER	DESCRIPTION
T_peak	Peak torque (kN·m). TYPE: float
T_remolded	Remolded torque (kN·m). TYPE: float
D	Vane diameter (m). TYPE: float
H	Vane height (m). TYPE: float

RETURNS	DESCRIPTION
dict	'Su_peak' : float — Peak Su (kPa). 'Su_remolded' : float — Remolded Su (kPa). 'sensitivity' : float — Sensitivity St. 'classification' : str — Sensitivity classification.

References

Das (2021), Table 3.7; ASTM D2573.

Examples:

Python Console Session

>>> from geoeq.site.vane import vane_remolded
>>> res = vane_remolded(T_peak=0.015, T_remolded=0.003, D=0.065, H=0.130)
>>> res['sensitivity']
5.0

Source code in geoeq/site/vane.py

def vane_remolded(
    T_peak: float,
    T_remolded: float,
    D: float,
    H: float,
) -> Dict[str, float]:
    r"""
    Undrained shear strength (peak and remolded) and sensitivity.

    After the peak torque test, the vane is rotated rapidly through
    several full revolutions to remold the soil, then retested.

    .. math::

        S_t = \frac{S_{u,\text{peak}}}{S_{u,\text{remolded}}}

    Parameters
    ----------
    T_peak : float
        Peak torque (kN·m).
    T_remolded : float
        Remolded torque (kN·m).
    D : float
        Vane diameter (m).
    H : float
        Vane height (m).

    Returns
    -------
    dict
        ``'Su_peak'`` : float — Peak Su (kPa).
        ``'Su_remolded'`` : float — Remolded Su (kPa).
        ``'sensitivity'`` : float — Sensitivity St.
        ``'classification'`` : str — Sensitivity classification.

    References
    ----------
    Das (2021), Table 3.7; ASTM D2573.

    Examples
    --------
    >>> from geoeq.site.vane import vane_remolded
    >>> res = vane_remolded(T_peak=0.015, T_remolded=0.003, D=0.065, H=0.130)
    >>> res['sensitivity']
    5.0
    """
    check_non_negative(T_peak, "T_peak")
    check_non_negative(T_remolded, "T_remolded")
    check_positive(D, "D")
    check_positive(H, "H")

    Su_peak = float(vane_su(T_peak, D, H))
    Su_rem = float(vane_su(T_remolded, D, H))

    if Su_rem > 0:
        St = Su_peak / Su_rem
    else:
        St = float("inf")

    if St < 2:
        cls_ = "Insensitive"
    elif St < 4:
        cls_ = "Low sensitivity"
    elif St < 8:
        cls_ = "Medium sensitivity"
    elif St < 16:
        cls_ = "Sensitive"
    else:
        cls_ = "Quick clay"

    return {
        "Su_peak": Su_peak,
        "Su_remolded": Su_rem,
        "sensitivity": float(St),
        "classification": cls_,
    }

Pressuremeter

pmt_parameters

Python

pmt_parameters(pressure: Union[list, ndarray], volume: Union[list, ndarray]) -> Dict[str, float]

Extract pressuremeter modulus Em and limit pressure pL from the pressure–volume curve.

The pseudo-elastic phase is identified as the steepest linear portion of the p–ΔV curve. The pressuremeter modulus is:

.. math::

Text Only

E_m = 2\,(1+\nu)\,(V_0 + V_m)\,\frac{\Delta p}{\Delta V}

where :math:V_0 is the initial probe volume, :math:V_m is the mid-volume of the linear phase, :math:\Delta p / \Delta V is the slope of the linear portion, and ν = 0.33.

The limit pressure :math:p_L is the pressure at which the cavity volume has doubled from its initial value.

PARAMETER	DESCRIPTION
pressure	Applied pressures (kPa), monotonically increasing. TYPE: array_like
volume	Corresponding injected volumes (cm³). TYPE: array_like

RETURNS	DESCRIPTION
dict	'Em' : float — Pressuremeter modulus (kPa). 'pL' : float — Limit pressure (kPa). 'p0' : float — Initial pressure (kPa) at start of elastic phase. 'pf' : float — Creep pressure (kPa) at end of elastic phase.

Notes

• The probe volume V₀ is typically 535 cm³ for the Ménard G-probe.
• This function uses an automated algorithm to detect the linear phase; for critical projects, manual selection is recommended.

References

Ménard (1956); Briaud (1992), Ch. 4.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_parameters
>>> p = [0, 50, 100, 200, 300, 400, 500, 600, 700, 800]
>>> v = [0, 20, 40, 60, 80, 110, 160, 240, 370, 600]
>>> res = pmt_parameters(p, v)
>>> res['Em'] > 0
True

Source code in geoeq/site/pmt.py

def pmt_parameters(
    pressure: Union[list, np.ndarray],
    volume: Union[list, np.ndarray],
) -> Dict[str, float]:
    r"""
    Extract pressuremeter modulus Em and limit pressure pL from the
    pressure–volume curve.

    The pseudo-elastic phase is identified as the steepest linear
    portion of the p–ΔV curve.  The pressuremeter modulus is:

    .. math::

        E_m = 2\,(1+\nu)\,(V_0 + V_m)\,\frac{\Delta p}{\Delta V}

    where :math:`V_0` is the initial probe volume, :math:`V_m` is the
    mid-volume of the linear phase, :math:`\Delta p / \Delta V` is
    the slope of the linear portion, and ν = 0.33.

    The limit pressure :math:`p_L` is the pressure at which the
    cavity volume has doubled from its initial value.

    Parameters
    ----------
    pressure : array_like
        Applied pressures (kPa), monotonically increasing.
    volume : array_like
        Corresponding injected volumes (cm³).

    Returns
    -------
    dict
        ``'Em'`` : float — Pressuremeter modulus (kPa).
        ``'pL'`` : float — Limit pressure (kPa).
        ``'p0'`` : float — Initial pressure (kPa) at start of elastic phase.
        ``'pf'`` : float — Creep pressure (kPa) at end of elastic phase.

    Notes
    -----
    * The probe volume V₀ is typically 535 cm³ for the Ménard G-probe.
    * This function uses an automated algorithm to detect the linear
      phase; for critical projects, manual selection is recommended.

    References
    ----------
    Ménard (1956); Briaud (1992), Ch. 4.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_parameters
    >>> p = [0, 50, 100, 200, 300, 400, 500, 600, 700, 800]
    >>> v = [0, 20, 40, 60, 80, 110, 160, 240, 370, 600]
    >>> res = pmt_parameters(p, v)
    >>> res['Em'] > 0
    True
    """
    p = np.asarray(pressure, dtype=float)
    v = np.asarray(volume, dtype=float)
    if len(p) != len(v):
        raise ValueError("pressure and volume must have the same length.")
    if len(p) < 4:
        raise ValueError("Need at least 4 data points.")

    idx = np.argsort(p)
    p = p[idx]
    v = v[idx]

    V0 = 535.0  # Ménard G-probe standard volume (cm³)
    nu = 0.33

    n = len(p)
    best_r2 = -1.0
    best_i0 = 0
    best_i1 = n - 1

    for i0 in range(n - 2):
        for i1 in range(i0 + 2, min(i0 + n, n)):
            seg_p = p[i0:i1 + 1]
            seg_v = v[i0:i1 + 1]
            if len(seg_p) < 3:
                continue
            coeffs = np.polyfit(seg_v, seg_p, 1)
            pred = np.polyval(coeffs, seg_v)
            ss_res = np.sum((seg_p - pred) ** 2)
            ss_tot = np.sum((seg_p - np.mean(seg_p)) ** 2)
            if ss_tot == 0:
                continue
            r2 = 1.0 - ss_res / ss_tot
            span = i1 - i0
            if r2 > 0.95 and span > best_i1 - best_i0:
                best_r2 = r2
                best_i0 = i0
                best_i1 = i1

    seg_v = v[best_i0:best_i1 + 1]
    seg_p = p[best_i0:best_i1 + 1]
    dp_dv = np.polyfit(seg_v, seg_p, 1)[0]

    Vm = (seg_v[0] + seg_v[-1]) / 2.0
    Em = 2.0 * (1.0 + nu) * (V0 + Vm) * dp_dv

    p0 = float(p[best_i0])
    pf = float(p[best_i1])

    # Limit pressure: pressure at which volume doubles from V0
    V_limit = V0
    if v[-1] >= V_limit:
        pL = float(np.interp(V_limit, v, p))
    else:
        pL = float(p[-1])

    return {
        "Em": float(Em),
        "pL": pL,
        "p0": p0,
        "pf": pf,
    }

pmt_modulus

Python

pmt_modulus(Em: float, alpha: float = 0.5) -> float

Young's modulus from pressuremeter modulus.

.. math::

Text Only

E = \frac{E_m}{\alpha}

PARAMETER	DESCRIPTION
Em	Pressuremeter modulus (kPa). TYPE: float
alpha	Rheological coefficient (Ménard, 1975). ============= ============= Soil type α ============= ============= Sand 0.33 -- 0.50 Silt 0.50 Clay 0.50 -- 0.67 Rock 0.33 ============= ============= TYPE: float DEFAULT: 0.5

RETURNS	DESCRIPTION
float	Young's modulus E (kPa).

References

Ménard (1975); Briaud (1992), Table 4.2.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_modulus
>>> pmt_modulus(Em=15000, alpha=0.5)
30000.0

Source code in geoeq/site/pmt.py

def pmt_modulus(
    Em: float,
    alpha: float = 0.5,
) -> float:
    r"""
    Young's modulus from pressuremeter modulus.

    .. math::

        E = \frac{E_m}{\alpha}

    Parameters
    ----------
    Em : float
        Pressuremeter modulus (kPa).
    alpha : float, default 0.5
        Rheological coefficient (Ménard, 1975).

        ============= =============
        Soil type     α
        ============= =============
        Sand          0.33 -- 0.50
        Silt          0.50
        Clay          0.50 -- 0.67
        Rock          0.33
        ============= =============

    Returns
    -------
    float
        Young's modulus E (kPa).

    References
    ----------
    Ménard (1975); Briaud (1992), Table 4.2.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_modulus
    >>> pmt_modulus(Em=15000, alpha=0.5)
    30000.0
    """
    check_positive(Em, "Em")
    check_positive(alpha, "alpha")
    return float(Em / alpha)

pmt_su

Python

pmt_su(pL: float, p0: float, sigma_h0: float = 0.0) -> float

Undrained shear strength from pressuremeter limit pressure.

.. math::

Text Only

S_u = \frac{p_L - p_0}{N_p}

where :math:N_p \approx 5.5 for the Ménard pressuremeter (Baguelin et al., 1978), and :math:p_0 can be approximated as the in-situ total horizontal stress.

Alternatively, using the net limit pressure:

.. math::

Text Only

p_L^* = p_L - \sigma_{h0}
\qquad
S_u = \frac{p_L^*}{N_p}

PARAMETER	DESCRIPTION
pL	Limit pressure (kPa). TYPE: float
p0	Initial pressure / in-situ horizontal stress (kPa). TYPE: float
sigma_h0	In-situ total horizontal stress (kPa). If provided and non-zero, the net limit pressure method is used. TYPE: float DEFAULT: 0.0

RETURNS	DESCRIPTION
float	Undrained shear strength Su (kPa).

References

Baguelin et al. (1978); Briaud (1992), Eq. 4.12.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_su
>>> round(pmt_su(pL=600, p0=100), 1)
90.9

Source code in geoeq/site/pmt.py

def pmt_su(
    pL: float,
    p0: float,
    sigma_h0: float = 0.0,
) -> float:
    r"""
    Undrained shear strength from pressuremeter limit pressure.

    .. math::

        S_u = \frac{p_L - p_0}{N_p}

    where :math:`N_p \approx 5.5` for the Ménard pressuremeter
    (Baguelin et al., 1978), and :math:`p_0` can be approximated
    as the in-situ total horizontal stress.

    Alternatively, using the net limit pressure:

    .. math::

        p_L^* = p_L - \sigma_{h0}
        \qquad
        S_u = \frac{p_L^*}{N_p}

    Parameters
    ----------
    pL : float
        Limit pressure (kPa).
    p0 : float
        Initial pressure / in-situ horizontal stress (kPa).
    sigma_h0 : float, default 0.0
        In-situ total horizontal stress (kPa).  If provided and
        non-zero, the net limit pressure method is used.

    Returns
    -------
    float
        Undrained shear strength Su (kPa).

    References
    ----------
    Baguelin et al. (1978); Briaud (1992), Eq. 4.12.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_su
    >>> round(pmt_su(pL=600, p0=100), 1)
    90.9
    """
    check_positive(pL, "pL")
    check_non_negative(p0, "p0")

    Np = 5.5
    if sigma_h0 > 0:
        Su = (pL - sigma_h0) / Np
    else:
        Su = (pL - p0) / Np

    return float(max(Su, 0.0))

pmt_bearing

Python

pmt_bearing(pL: float, p0: float, sigma_v: float, shape: str = 'strip') -> Dict[str, float]

Ultimate bearing capacity from pressuremeter data (Ménard method).

.. math::

Text Only

q_u = \sigma_{v0} + k \,(p_L - p_0)

where k is a bearing factor depending on foundation shape and embedment:

============= ==== Shape k ============= ==== Strip 0.8 Square/Circle 1.2 ============= ====

PARAMETER	DESCRIPTION
pL	Limit pressure (kPa). TYPE: float
p0	Initial pressure (kPa). TYPE: float
sigma_v	Overburden stress at foundation level (kPa). TYPE: float
shape	'strip' or 'square'/'circle'. TYPE: str DEFAULT: ``'strip'``

RETURNS	DESCRIPTION
dict	'qu' : float — Ultimate bearing capacity (kPa). 'net_limit_pressure' : float — pL* = pL - p0 (kPa). 'k' : float — Bearing factor used.

References

Ménard (1963); Briaud (1992), Section 5.3.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_bearing
>>> res = pmt_bearing(pL=800, p0=100, sigma_v=50, shape='square')
>>> res['qu']
890.0

Source code in geoeq/site/pmt.py

def pmt_bearing(
    pL: float,
    p0: float,
    sigma_v: float,
    shape: str = "strip",
) -> Dict[str, float]:
    r"""
    Ultimate bearing capacity from pressuremeter data (Ménard method).

    .. math::

        q_u = \sigma_{v0} + k \,(p_L - p_0)

    where k is a bearing factor depending on foundation shape and
    embedment:

    ============= ====
    Shape         k
    ============= ====
    Strip         0.8
    Square/Circle 1.2
    ============= ====

    Parameters
    ----------
    pL : float
        Limit pressure (kPa).
    p0 : float
        Initial pressure (kPa).
    sigma_v : float
        Overburden stress at foundation level (kPa).
    shape : str, default ``'strip'``
        ``'strip'`` or ``'square'``/``'circle'``.

    Returns
    -------
    dict
        ``'qu'`` : float — Ultimate bearing capacity (kPa).
        ``'net_limit_pressure'`` : float — pL* = pL - p0 (kPa).
        ``'k'`` : float — Bearing factor used.

    References
    ----------
    Ménard (1963); Briaud (1992), Section 5.3.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_bearing
    >>> res = pmt_bearing(pL=800, p0=100, sigma_v=50, shape='square')
    >>> res['qu']
    890.0
    """
    check_positive(pL, "pL")
    check_non_negative(p0, "p0")
    check_non_negative(sigma_v, "sigma_v")

    shape_l = shape.lower()
    if shape_l in ("strip", "continuous"):
        k = 0.8
    elif shape_l in ("square", "circle", "circular", "rectangular"):
        k = 1.2
    else:
        raise ValueError(f"Unknown shape '{shape}'. Choose: strip, square, circle.")

    pL_star = pL - p0
    qu = sigma_v + k * pL_star

    return {
        "qu": float(qu),
        "net_limit_pressure": float(pL_star),
        "k": float(k),
    }

pmt_settlement

Python

pmt_settlement(q: float, B: float, Em: float, alpha: float = 0.5, shape: str = 'circle', B0: float = 0.6) -> float

Settlement from Ménard pressuremeter method.

.. math::

Text Only

s = \frac{2}{9\,E_m}\,q\,B_0\,
    \left(\frac{\alpha \, B}{B_0}\right)^{\!\beta}

where β = shape factor (α from Ménard, here separate from the rheological factor). The simplified Ménard (1963) approach:

.. math::

Text Only

s = \frac{q \, B}{E / \alpha_s}

Using the deviatoric and spherical components (full Ménard approach):

.. math::

Text Only

s = \frac{q}{9\,E_m}\bigl[2\,B_0\,\alpha_d\,(B/B_0)^{\alpha_s}
    + \alpha_s\,B\bigr]

This function uses the simplified version for routine estimates.

PARAMETER	DESCRIPTION
q	Net applied pressure (kPa). TYPE: float
B	Foundation width or diameter (m). TYPE: float
Em	Pressuremeter modulus (kPa). TYPE: float
alpha	Rheological coefficient. TYPE: float DEFAULT: 0.5
shape	'circle' → λ_d = 1.0, 'square' → λ_d = 1.12, 'strip' → λ_d = 1.53. TYPE: str DEFAULT: ``'circle'``
B0	Reference width (m) = 0.60 m = 60 cm. TYPE: float DEFAULT: 0.6

RETURNS	DESCRIPTION
float	Settlement (m).

References

Ménard (1963); Briaud (1992), Section 6.2.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_settlement
>>> round(pmt_settlement(q=150, B=2.0, Em=15000) * 1000, 1)
10.7

Source code in geoeq/site/pmt.py

def pmt_settlement(
    q: float,
    B: float,
    Em: float,
    alpha: float = 0.5,
    shape: str = "circle",
    B0: float = 0.6,
) -> float:
    r"""
    Settlement from Ménard pressuremeter method.

    .. math::

        s = \frac{2}{9\,E_m}\,q\,B_0\,
            \left(\frac{\alpha \, B}{B_0}\right)^{\!\beta}

    where β = shape factor (α from Ménard, here separate from the
    rheological factor).  The simplified Ménard (1963) approach:

    .. math::

        s = \frac{q \, B}{E / \alpha_s}

    Using the deviatoric and spherical components
    (full Ménard approach):

    .. math::

        s = \frac{q}{9\,E_m}\bigl[2\,B_0\,\alpha_d\,(B/B_0)^{\alpha_s}
            + \alpha_s\,B\bigr]

    This function uses the simplified version for routine estimates.

    Parameters
    ----------
    q : float
        Net applied pressure (kPa).
    B : float
        Foundation width or diameter (m).
    Em : float
        Pressuremeter modulus (kPa).
    alpha : float, default 0.5
        Rheological coefficient.
    shape : str, default ``'circle'``
        ``'circle'`` → λ_d = 1.0, ``'square'`` → λ_d = 1.12,
        ``'strip'`` → λ_d = 1.53.
    B0 : float, default 0.6
        Reference width (m) = 0.60 m = 60 cm.

    Returns
    -------
    float
        Settlement (m).

    References
    ----------
    Ménard (1963); Briaud (1992), Section 6.2.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_settlement
    >>> round(pmt_settlement(q=150, B=2.0, Em=15000) * 1000, 1)
    10.7
    """
    check_positive(q, "q")
    check_positive(B, "B")
    check_positive(Em, "Em")
    check_positive(alpha, "alpha")

    shape_l = shape.lower()
    lam = {"circle": 1.0, "square": 1.12, "strip": 1.53}.get(shape_l, 1.0)

    # Deviatoric component
    s_d = (2.0 * q * B0) / (9.0 * Em) * lam * (B / B0) ** alpha
    # Spherical component
    s_s = (q * alpha * B) / (9.0 * Em)

    return float(s_d + s_s)

pmt_ko

Python

pmt_ko(p0: float, sigma_v: float, u0: float = 0.0) -> float

At-rest earth pressure coefficient from PMT initial pressure.

.. math::

Text Only

K_0 = \frac{p_0 - u_0}{\sigma'_{v0}}

PARAMETER	DESCRIPTION
p0	In-situ horizontal stress from PMT (kPa). TYPE: float
sigma_v	Total vertical stress (kPa). TYPE: float
u0	Pore water pressure (kPa). TYPE: float DEFAULT: 0.0

RETURNS	DESCRIPTION
float	Coefficient of earth pressure at rest K₀.

Examples:

Python Console Session

>>> from geoeq.site.pmt import pmt_ko
>>> pmt_ko(p0=80, sigma_v=150, u0=30)
0.4166666666666667

Source code in geoeq/site/pmt.py

def pmt_ko(
    p0: float,
    sigma_v: float,
    u0: float = 0.0,
) -> float:
    r"""
    At-rest earth pressure coefficient from PMT initial pressure.

    .. math::

        K_0 = \frac{p_0 - u_0}{\sigma'_{v0}}

    Parameters
    ----------
    p0 : float
        In-situ horizontal stress from PMT (kPa).
    sigma_v : float
        Total vertical stress (kPa).
    u0 : float, default 0.0
        Pore water pressure (kPa).

    Returns
    -------
    float
        Coefficient of earth pressure at rest K₀.

    Examples
    --------
    >>> from geoeq.site.pmt import pmt_ko
    >>> pmt_ko(p0=80, sigma_v=150, u0=30)
    0.4166666666666667
    """
    check_non_negative(p0, "p0")
    check_positive(sigma_v, "sigma_v")
    check_non_negative(u0, "u0")

    sigma_v_eff = sigma_v - u0
    if sigma_v_eff <= 0:
        raise ValueError("Effective vertical stress must be positive.")

    return float((p0 - u0) / sigma_v_eff)

Plate load test

plt_bearing

Python

plt_bearing(pressure: Union[list, ndarray], settlement: Union[list, ndarray], failure_criterion: str = 'log_log') -> Dict[str, float]

Determine ultimate bearing capacity from plate load test data.

PARAMETER	DESCRIPTION
pressure	Applied pressures (kPa). TYPE: array_like
settlement	Corresponding settlements (mm). TYPE: array_like
failure_criterion	'log_log' : De Beer's double-log method — intersection of two straight-line segments in log–log space. 'tangent' : Tangent intersection method — failure at the intersection of the initial tangent and the final tangent. 'settlement' : Failure at settlement = B/10 (where B = 0.30 m standard plate → 30 mm). TYPE: str DEFAULT: ``'log_log'``

RETURNS	DESCRIPTION
dict	'qu' : float — Ultimate bearing capacity (kPa). 'settlement_at_failure' : float — Settlement at failure (mm).

Notes

• Standard plate diameter is 300 mm (12 in.), but 450 mm and 600 mm (18/24 in.) plates are also used.
• For clay soils, the plate test gives the bearing capacity of the plate, which equals that of the actual foundation (independent of size for φ = 0).
• For sandy soils, scale correction is required — see :func:plt_settlement_correction.

References

Das (2021), Section 4.14; IS 1888 (1982).

Examples:

Python Console Session

>>> from geoeq.site.plt_test import plt_bearing
>>> p = [0, 50, 100, 150, 200, 250, 300, 350, 400]
>>> s = [0, 0.5, 1.2, 2.5, 5.0, 10.0, 20.0, 35.0, 55.0]
>>> res = plt_bearing(p, s, failure_criterion='settlement')
>>> res['qu'] > 0
True

Source code in geoeq/site/plt_test.py

def plt_bearing(
    pressure: Union[list, np.ndarray],
    settlement: Union[list, np.ndarray],
    failure_criterion: str = "log_log",
) -> Dict[str, float]:
    r"""
    Determine ultimate bearing capacity from plate load test data.

    Parameters
    ----------
    pressure : array_like
        Applied pressures (kPa).
    settlement : array_like
        Corresponding settlements (mm).
    failure_criterion : str, default ``'log_log'``
        - ``'log_log'`` : De Beer's double-log method — intersection
          of two straight-line segments in log–log space.
        - ``'tangent'`` : Tangent intersection method — failure at the
          intersection of the initial tangent and the final tangent.
        - ``'settlement'`` : Failure at settlement = B/10
          (where B = 0.30 m standard plate → 30 mm).

    Returns
    -------
    dict
        ``'qu'`` : float — Ultimate bearing capacity (kPa).
        ``'settlement_at_failure'`` : float — Settlement at failure (mm).

    Notes
    -----
    * Standard plate diameter is 300 mm (12 in.), but 450 mm and
      600 mm (18/24 in.) plates are also used.
    * For clay soils, the plate test gives the **bearing capacity of
      the plate**, which equals that of the actual foundation
      (independent of size for φ = 0).
    * For sandy soils, scale correction is required — see
      :func:`plt_settlement_correction`.

    References
    ----------
    Das (2021), Section 4.14; IS 1888 (1982).

    Examples
    --------
    >>> from geoeq.site.plt_test import plt_bearing
    >>> p = [0, 50, 100, 150, 200, 250, 300, 350, 400]
    >>> s = [0, 0.5, 1.2, 2.5, 5.0, 10.0, 20.0, 35.0, 55.0]
    >>> res = plt_bearing(p, s, failure_criterion='settlement')
    >>> res['qu'] > 0
    True
    """
    p = np.asarray(pressure, dtype=float)
    s = np.asarray(settlement, dtype=float)
    if len(p) != len(s):
        raise ValueError("pressure and settlement must have the same length.")
    if len(p) < 3:
        raise ValueError("Need at least 3 data points.")

    criterion = failure_criterion.lower().replace("-", "_").replace(" ", "_")

    if criterion == "settlement":
        # Failure at 30 mm (B/10 for standard 300 mm plate)
        s_fail = 30.0
        mask = s > 0
        if np.any(s >= s_fail):
            qu = float(np.interp(s_fail, s, p))
        else:
            qu = float("nan")
        return {"qu": qu, "settlement_at_failure": s_fail}

    elif criterion == "log_log":
        mask = (p > 0) & (s > 0)
        p_m = p[mask]
        s_m = s[mask]
        if len(p_m) < 4:
            raise ValueError("Need at least 4 positive data points for log-log method.")

        logp = np.log10(p_m)
        logs = np.log10(s_m)

        n = len(p_m)
        best_ssr = float("inf")
        best_k = 2
        for k in range(2, n - 1):
            c1 = np.polyfit(logp[:k + 1], logs[:k + 1], 1)
            c2 = np.polyfit(logp[k:], logs[k:], 1)
            r1 = logs[:k + 1] - np.polyval(c1, logp[:k + 1])
            r2 = logs[k:] - np.polyval(c2, logp[k:])
            ssr = np.sum(r1**2) + np.sum(r2**2)
            if ssr < best_ssr:
                best_ssr = ssr
                best_k = k

        c1 = np.polyfit(logp[:best_k + 1], logs[:best_k + 1], 1)
        c2 = np.polyfit(logp[best_k:], logs[best_k:], 1)

        if abs(c1[0] - c2[0]) < 1e-12:
            return {"qu": float("nan"), "settlement_at_failure": float("nan")}

        logp_int = (c2[1] - c1[1]) / (c1[0] - c2[0])
        logs_int = c1[0] * logp_int + c1[1]

        qu = float(10**logp_int)
        s_fail = float(10**logs_int)
        return {"qu": qu, "settlement_at_failure": s_fail}

    elif criterion == "tangent":
        # Initial and final tangent intersection
        n = len(p)
        # Initial tangent from first few points
        n_init = min(3, n // 2)
        c_init = np.polyfit(s[:n_init], p[:n_init], 1)
        # Final tangent from last few points
        n_fin = min(3, n // 2)
        c_fin = np.polyfit(s[-n_fin:], p[-n_fin:], 1)

        if abs(c_init[0] - c_fin[0]) < 1e-12:
            return {"qu": float("nan"), "settlement_at_failure": float("nan")}

        s_int = (c_fin[1] - c_init[1]) / (c_init[0] - c_fin[0])
        qu = float(c_init[0] * s_int + c_init[1])
        return {"qu": max(qu, 0.0), "settlement_at_failure": float(max(s_int, 0.0))}

    else:
        raise ValueError(
            f"Unknown criterion '{failure_criterion}'. "
            "Choose: log_log, tangent, settlement."
        )

plt_subgrade_modulus

Python

plt_subgrade_modulus(pressure: Union[float, ndarray], settlement: Union[float, ndarray], plate_size: float = 0.3) -> Union[float, Dict[str, float]]

Modulus of subgrade reaction from plate load test.

.. math::

Text Only

k_s = \frac{q}{\delta}

where q is the applied pressure (kPa) and δ is the settlement (m).

For a standard 300 mm plate, correction to foundation size:

.. math::

Text Only

k_{s,\text{foundation}} = k_{s,\text{plate}} \times
\left(\frac{B_p}{B_f}\right)
\qquad \text{(clay, φ = 0)}

.. math::

Text Only

k_{s,\text{foundation}} = k_{s,\text{plate}} \times
\left(\frac{B_f + B_p}{2 \, B_f}\right)^{\!2}
\qquad \text{(sand)}

PARAMETER	DESCRIPTION
pressure	Applied pressure(s) (kPa). TYPE: float or array_like
settlement	Corresponding settlement(s) (mm). TYPE: float or array_like
plate_size	Plate diameter or width (m). TYPE: float DEFAULT: 0.3

RETURNS	DESCRIPTION
float or dict	If scalar inputs: float — ks (kN/m³). If array inputs: dict with 'ks' array and 'mean_ks'.

Notes

Settlement is converted from mm to m internally.

References

Das (2021), Eq. 4.35; Terzaghi (1955).

Examples:

Python Console Session

>>> from geoeq.site.plt_test import plt_subgrade_modulus
>>> plt_subgrade_modulus(pressure=200, settlement=5.0)
40000.0

Source code in geoeq/site/plt_test.py

def plt_subgrade_modulus(
    pressure: Union[float, np.ndarray],
    settlement: Union[float, np.ndarray],
    plate_size: float = 0.3,
) -> Union[float, Dict[str, float]]:
    r"""
    Modulus of subgrade reaction from plate load test.

    .. math::

        k_s = \frac{q}{\delta}

    where q is the applied pressure (kPa) and δ is the settlement (m).

    For a standard 300 mm plate, correction to foundation size:

    .. math::

        k_{s,\text{foundation}} = k_{s,\text{plate}} \times
        \left(\frac{B_p}{B_f}\right)
        \qquad \text{(clay, φ = 0)}

    .. math::

        k_{s,\text{foundation}} = k_{s,\text{plate}} \times
        \left(\frac{B_f + B_p}{2 \, B_f}\right)^{\!2}
        \qquad \text{(sand)}

    Parameters
    ----------
    pressure : float or array_like
        Applied pressure(s) (kPa).
    settlement : float or array_like
        Corresponding settlement(s) (mm).
    plate_size : float, default 0.3
        Plate diameter or width (m).

    Returns
    -------
    float or dict
        If scalar inputs: float — ks (kN/m³).
        If array inputs: dict with ``'ks'`` array and ``'mean_ks'``.

    Notes
    -----
    Settlement is converted from mm to m internally.

    References
    ----------
    Das (2021), Eq. 4.35; Terzaghi (1955).

    Examples
    --------
    >>> from geoeq.site.plt_test import plt_subgrade_modulus
    >>> plt_subgrade_modulus(pressure=200, settlement=5.0)
    40000.0
    """
    check_positive(plate_size, "plate_size")

    if np.ndim(pressure) == 0 and np.ndim(settlement) == 0:
        check_positive(pressure, "pressure")
        check_positive(settlement, "settlement")
        s_m = settlement / 1000.0  # mm → m
        return float(pressure / s_m)

    p = np.asarray(pressure, dtype=float)
    s = np.asarray(settlement, dtype=float)
    check_positive(p, "pressure")
    check_positive(s, "settlement")

    s_m = s / 1000.0
    ks = p / s_m
    return {"ks": ks, "mean_ks": float(np.mean(ks))}

plt_settlement_correction

Python

plt_settlement_correction(s_plate: float, B_plate: float, B_foundation: float, soil_type: str = 'sand') -> float

Terzaghi & Peck settlement correction from plate to foundation.

For sand (Terzaghi & Peck, 1967):

.. math::

Text Only

\frac{S_f}{S_p} = \left(\frac{B_f \,(B_p + 0.3)}
                               {B_p \,(B_f + 0.3)}\right)^{\!2}

For clay (φ = 0, bearing capacity independent of size):

.. math::

Text Only

\frac{S_f}{S_p} = \frac{B_f}{B_p}

PARAMETER	DESCRIPTION
s_plate	Settlement of the test plate (mm). TYPE: float
B_plate	Plate size (m). TYPE: float
B_foundation	Foundation size (m). TYPE: float
soil_type	'sand' or 'clay'. TYPE: str DEFAULT: ``'sand'``

RETURNS	DESCRIPTION
float	Estimated foundation settlement (mm).

References

Terzaghi & Peck (1967); Das (2021), Eq. 4.34.

Examples:

Python Console Session

>>> from geoeq.site.plt_test import plt_settlement_correction
>>> round(plt_settlement_correction(s_plate=5.0, B_plate=0.3,
...       B_foundation=2.0, soil_type='sand'), 1)
17.7

Source code in geoeq/site/plt_test.py

def plt_settlement_correction(
    s_plate: float,
    B_plate: float,
    B_foundation: float,
    soil_type: str = "sand",
) -> float:
    r"""
    Terzaghi & Peck settlement correction from plate to foundation.

    For **sand** (Terzaghi & Peck, 1967):

    .. math::

        \frac{S_f}{S_p} = \left(\frac{B_f \,(B_p + 0.3)}
                                       {B_p \,(B_f + 0.3)}\right)^{\!2}

    For **clay** (φ = 0, bearing capacity independent of size):

    .. math::

        \frac{S_f}{S_p} = \frac{B_f}{B_p}

    Parameters
    ----------
    s_plate : float
        Settlement of the test plate (mm).
    B_plate : float
        Plate size (m).
    B_foundation : float
        Foundation size (m).
    soil_type : str, default ``'sand'``
        ``'sand'`` or ``'clay'``.

    Returns
    -------
    float
        Estimated foundation settlement (mm).

    References
    ----------
    Terzaghi & Peck (1967); Das (2021), Eq. 4.34.

    Examples
    --------
    >>> from geoeq.site.plt_test import plt_settlement_correction
    >>> round(plt_settlement_correction(s_plate=5.0, B_plate=0.3,
    ...       B_foundation=2.0, soil_type='sand'), 1)
    17.7
    """
    check_non_negative(s_plate, "s_plate")
    check_positive(B_plate, "B_plate")
    check_positive(B_foundation, "B_foundation")

    Bp = B_plate
    Bf = B_foundation

    soil = soil_type.lower()
    if soil == "sand":
        ratio = (Bf * (Bp + 0.3) / (Bp * (Bf + 0.3))) ** 2
    elif soil == "clay":
        ratio = Bf / Bp
    else:
        raise ValueError(f"Unknown soil_type '{soil_type}'. Choose: sand, clay.")

    return float(s_plate * ratio)

plt_elastic_modulus

Python

plt_elastic_modulus(pressure: float, settlement: float, plate_diameter: float = 0.3, poisson_ratio: float = 0.3, Is: float = 0.79) -> float

Elastic modulus from plate load test (Timoshenko & Goodier).

.. math::

Text Only

E = \frac{q \, B \, (1 - \nu^2) \, I_s}{\delta}

PARAMETER	DESCRIPTION
pressure	Applied pressure (kPa). TYPE: float
settlement	Settlement (mm). TYPE: float
plate_diameter	Plate diameter (m). TYPE: float DEFAULT: 0.3
poisson_ratio	Poisson's ratio. TYPE: float DEFAULT: 0.3
Is	Influence factor. 0.79 for rigid circular plate, π/4 for flexible circular plate. TYPE: float DEFAULT: 0.79

RETURNS	DESCRIPTION
float	Elastic modulus E (kPa).

References

Timoshenko & Goodier (1970); Das (2021), Eq. 4.33.

Examples:

Python Console Session

>>> from geoeq.site.plt_test import plt_elastic_modulus
>>> round(plt_elastic_modulus(200, 5.0, plate_diameter=0.3), 0)
8658.0

Source code in geoeq/site/plt_test.py

def plt_elastic_modulus(
    pressure: float,
    settlement: float,
    plate_diameter: float = 0.3,
    poisson_ratio: float = 0.3,
    Is: float = 0.79,
) -> float:
    r"""
    Elastic modulus from plate load test (Timoshenko & Goodier).

    .. math::

        E = \frac{q \, B \, (1 - \nu^2) \, I_s}{\delta}

    Parameters
    ----------
    pressure : float
        Applied pressure (kPa).
    settlement : float
        Settlement (mm).
    plate_diameter : float, default 0.3
        Plate diameter (m).
    poisson_ratio : float, default 0.3
        Poisson's ratio.
    Is : float, default 0.79
        Influence factor. 0.79 for rigid circular plate,
        π/4 for flexible circular plate.

    Returns
    -------
    float
        Elastic modulus E (kPa).

    References
    ----------
    Timoshenko & Goodier (1970); Das (2021), Eq. 4.33.

    Examples
    --------
    >>> from geoeq.site.plt_test import plt_elastic_modulus
    >>> round(plt_elastic_modulus(200, 5.0, plate_diameter=0.3), 0)
    8658.0
    """
    check_positive(pressure, "pressure")
    check_positive(settlement, "settlement")
    check_positive(plate_diameter, "plate_diameter")

    s_m = settlement / 1000.0  # mm → m
    E = pressure * plate_diameter * (1.0 - poisson_ratio**2) * Is / s_m
    return float(E)

plt_plot

Python

plt_plot(pressure: Union[list, ndarray], settlement: Union[list, ndarray], ax=None, save_as: Optional[str] = None) -> Dict

Plot plate load test pressure–settlement curve.

PARAMETER	DESCRIPTION
pressure	Applied pressures (kPa). TYPE: array_like
settlement	Corresponding settlements (mm). TYPE: array_like
ax	TYPE: Axes DEFAULT: None
save_as	TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
dict	'ax' : Axes.

Examples:

Python Console Session

>>> from geoeq.site.plt_test import plt_plot
>>> res = plt_plot([0, 100, 200, 300], [0, 2, 5, 12])

Source code in geoeq/site/plt_test.py

def plt_plot(
    pressure: Union[list, np.ndarray],
    settlement: Union[list, np.ndarray],
    ax=None,
    save_as: Optional[str] = None,
) -> Dict:
    r"""
    Plot plate load test pressure–settlement curve.

    Parameters
    ----------
    pressure : array_like
        Applied pressures (kPa).
    settlement : array_like
        Corresponding settlements (mm).
    ax : matplotlib.axes.Axes, optional
    save_as : str, optional

    Returns
    -------
    dict
        ``'ax'`` : Axes.

    Examples
    --------
    >>> from geoeq.site.plt_test import plt_plot
    >>> res = plt_plot([0, 100, 200, 300], [0, 2, 5, 12])
    """
    import matplotlib.pyplot as plt

    p = np.asarray(pressure, dtype=float)
    s = np.asarray(settlement, dtype=float)

    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 6))
    else:
        fig = ax.get_figure()

    ax.plot(p, s, "o-", color="steelblue", linewidth=1.5, markersize=6,
            markerfacecolor="white", markeredgewidth=1.5)
    ax.set_xlabel("Applied Pressure (kPa)", fontweight="bold")
    ax.set_ylabel("Settlement (mm)", fontweight="bold")
    ax.set_title("Plate Load Test — Pressure vs Settlement", fontweight="bold")
    ax.invert_yaxis()
    ax.grid(True, alpha=0.3)

    if save_as:
        plt.savefig(save_as, bbox_inches="tight", dpi=300)

    return {"ax": ax}

Pile load tests

davisson

Python

davisson(load: Union[List[float], ndarray], settlement: Union[List[float], ndarray], diameter: float, length: float, area: float, elastic_modulus: float) -> Dict[str, float]

Davisson's Offset Limit Method (1972).

The failure load is defined as the load at which the measured load-settlement curve intersects an offset elastic compression line.

The offset line is:

.. math::

Text Only

\Delta_{\text{offset}} = \frac{P \, L}{A \, E}
    + \left(0.15 + \frac{D}{120}\right) \times 25.4 \; \text{[mm]}

where the term in parentheses is in inches and :math:D is the pile diameter in inches. The constant 0.15 in. (3.81 mm) is the quake offset, and :math:D/120 accounts for pile-tip displacement.

Equivalently in SI (all mm):

.. math::

Text Only

\Delta_{\text{offset}} = \frac{P \, L}{A \, E}
    + 3.81 + \frac{D}{4.724}

PARAMETER	DESCRIPTION
load	Applied loads (kN), monotonically increasing. TYPE: array_like
settlement	Corresponding pile-head settlements (mm). TYPE: array_like
diameter	Pile diameter or equivalent diameter (mm). TYPE: float
length	Pile length (mm). TYPE: float
area	Cross-sectional area of the pile (mm^2). TYPE: float
elastic_modulus	Elastic modulus of the pile material (kN/mm^2 = GPa). Typical: steel 200 GPa, concrete 25-35 GPa. TYPE: float

RETURNS	DESCRIPTION
dict	'Qu_davisson' : float Davisson failure load (kN), or nan if no intersection. 'settlement_at_failure' : float Settlement at the failure load (mm). 'elastic_compression_offset' : float The fixed offset :math:\Delta_0 (mm) = 3.81 + D/4.724.

Notes

• Originally formulated for driven piles up to 600 mm (24 in.) diameter.
• For larger piles, consider the FHWA 5 %% criterion (:func:fhwa_5_percent).
• The elastic compression slope is :math:L / (A \, E).

References

Davisson (1972); Fellenius (2020), Ch. 8.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import davisson
>>> Q = [0, 200, 400, 600, 800, 1000, 1200, 1400]
>>> s = [0, 1.0, 2.5, 5.0, 8.5, 14.0, 22.0, 35.0]
>>> res = davisson(Q, s, diameter=300, length=12000, area=70686, elastic_modulus=30)
>>> 0 < res['Qu_davisson'] < 1500
True

Source code in geoeq/site/pile_load.py

def davisson(
    load: Union[List[float], np.ndarray],
    settlement: Union[List[float], np.ndarray],
    diameter: float,
    length: float,
    area: float,
    elastic_modulus: float,
) -> Dict[str, float]:
    r"""
    Davisson's Offset Limit Method (1972).

    The failure load is defined as the load at which the measured
    load-settlement curve intersects an *offset elastic compression line*.

    The offset line is:

    .. math::

        \Delta_{\text{offset}} = \frac{P \, L}{A \, E}
            + \left(0.15 + \frac{D}{120}\right) \times 25.4 \; \text{[mm]}

    where the term in parentheses is in **inches** and :math:`D` is the
    pile diameter in **inches**.  The constant 0.15 in. (3.81 mm) is the
    quake offset, and :math:`D/120` accounts for pile-tip displacement.

    Equivalently in SI (all mm):

    .. math::

        \Delta_{\text{offset}} = \frac{P \, L}{A \, E}
            + 3.81 + \frac{D}{4.724}

    Parameters
    ----------
    load : array_like
        Applied loads (kN), monotonically increasing.
    settlement : array_like
        Corresponding pile-head settlements (mm).
    diameter : float
        Pile diameter or equivalent diameter (mm).
    length : float
        Pile length (mm).
    area : float
        Cross-sectional area of the pile (mm^2).
    elastic_modulus : float
        Elastic modulus of the pile material (kN/mm^2 = GPa).
        Typical: steel 200 GPa, concrete 25-35 GPa.

    Returns
    -------
    dict
        ``'Qu_davisson'`` : float
            Davisson failure load (kN), or ``nan`` if no intersection.
        ``'settlement_at_failure'`` : float
            Settlement at the failure load (mm).
        ``'elastic_compression_offset'`` : float
            The fixed offset :math:`\Delta_0` (mm) = 3.81 + D/4.724.

    Notes
    -----
    * Originally formulated for driven piles up to 600 mm (24 in.) diameter.
    * For larger piles, consider the FHWA 5 %% criterion
      (:func:`fhwa_5_percent`).
    * The elastic compression slope is :math:`L / (A \, E)`.

    References
    ----------
    Davisson (1972); Fellenius (2020), Ch. 8.

    Examples
    --------
    >>> from geoeq.site.pile_load import davisson
    >>> Q = [0, 200, 400, 600, 800, 1000, 1200, 1400]
    >>> s = [0, 1.0, 2.5, 5.0, 8.5, 14.0, 22.0, 35.0]
    >>> res = davisson(Q, s, diameter=300, length=12000, area=70686, elastic_modulus=30)
    >>> 0 < res['Qu_davisson'] < 1500
    True
    """
    Q = np.asarray(load, dtype=float)
    s = np.asarray(settlement, dtype=float)
    if len(Q) != len(s):
        raise ValueError("load and settlement must have the same length.")
    if len(Q) < 3:
        raise ValueError("Need at least 3 data points.")
    check_positive(diameter, "diameter")
    check_positive(length, "length")
    check_positive(area, "area")
    check_positive(elastic_modulus, "elastic_modulus")

    # Sort by load
    idx = np.argsort(Q)
    Q = Q[idx]
    s = s[idx]

    # Convert diameter mm -> inches for offset formula, then back to mm
    D_in = diameter / 25.4
    offset_in = 0.15 + D_in / 120.0  # inches
    offset_mm = offset_in * 25.4      # mm

    # Elastic compression slope: delta_elastic = P * L / (A * E)
    # slope_elastic has units mm/kN
    slope_elastic = length / (area * elastic_modulus)

    # Offset line: s_offset(P) = P * slope_elastic + offset_mm
    s_offset = Q * slope_elastic + offset_mm

    # Find intersection: where measured settlement crosses the offset line
    # diff = s_measured - s_offset; intersection where diff changes sign
    diff = s - s_offset

    Qu = float("nan")
    s_fail = float("nan")

    for i in range(1, len(diff)):
        if diff[i - 1] <= 0 and diff[i] > 0:
            # Linear interpolation
            frac = diff[i - 1] / (diff[i - 1] - diff[i])
            Qu = float(Q[i - 1] + frac * (Q[i] - Q[i - 1]))
            s_fail = float(s[i - 1] + frac * (s[i] - s[i - 1]))
            break

    return {
        "Qu_davisson": Qu,
        "settlement_at_failure": s_fail,
        "elastic_compression_offset": offset_mm,
    }

chin

Python

chin(load: Union[List[float], ndarray], settlement: Union[List[float], ndarray], start_index: Optional[int] = None, end_index: Optional[int] = None) -> Dict[str, float]

Chin's Hyperbolic Method (1970).

Plot :math:\Delta / Q (settlement / load) versus :math:\Delta (settlement). The data approach a straight line at larger settlements. The ultimate capacity is:

.. math::

Text Only

Q_u = \frac{1}{C_1}

where :math:C_1 is the slope of the linear portion of :math:\Delta / Q versus :math:\Delta.

PARAMETER	DESCRIPTION
load	Applied loads (kN). Must be > 0 for the data points used. TYPE: array_like
settlement	Corresponding settlements (mm). TYPE: array_like
start_index	Index at which the linear portion begins (0-based). If None, uses the last 50 %% of data points. TYPE: int DEFAULT: None
end_index	Index at which the linear portion ends (exclusive). TYPE: int DEFAULT: None

RETURNS	DESCRIPTION
dict	'Qu_chin' : float Chin ultimate load (kN). 'slope' : float Slope C1 of the linear regression (mm/kN per mm = 1/kN). 'intercept' : float Intercept C2 of the linear regression. 'r_squared' : float Coefficient of determination of the linear fit.

Notes

• Chin's method typically overestimates the actual failure load by 20--40 %%. Apply a reduction factor (commonly 0.80) for design.
• Only data points where Q > 0 are used.

References

Chin (1970); Fellenius (2020), Section 8.9.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import chin
>>> Q = [100, 200, 400, 600, 800, 1000, 1200]
>>> s = [0.5, 1.2, 3.0, 5.5, 9.0, 14.0, 21.0]
>>> res = chin(Q, s)
>>> res['Qu_chin'] > 0
True

Source code in geoeq/site/pile_load.py

def chin(
    load: Union[List[float], np.ndarray],
    settlement: Union[List[float], np.ndarray],
    start_index: Optional[int] = None,
    end_index: Optional[int] = None,
) -> Dict[str, float]:
    r"""
    Chin's Hyperbolic Method (1970).

    Plot :math:`\Delta / Q` (settlement / load) versus :math:`\Delta`
    (settlement).  The data approach a straight line at larger
    settlements.  The ultimate capacity is:

    .. math::

        Q_u = \frac{1}{C_1}

    where :math:`C_1` is the **slope** of the linear portion of
    :math:`\Delta / Q` versus :math:`\Delta`.

    Parameters
    ----------
    load : array_like
        Applied loads (kN).  Must be > 0 for the data points used.
    settlement : array_like
        Corresponding settlements (mm).
    start_index : int, optional
        Index at which the linear portion begins (0-based).
        If ``None``, uses the last 50 %% of data points.
    end_index : int, optional
        Index at which the linear portion ends (exclusive).

    Returns
    -------
    dict
        ``'Qu_chin'`` : float
            Chin ultimate load (kN).
        ``'slope'`` : float
            Slope C1 of the linear regression (mm/kN per mm = 1/kN).
        ``'intercept'`` : float
            Intercept C2 of the linear regression.
        ``'r_squared'`` : float
            Coefficient of determination of the linear fit.

    Notes
    -----
    * Chin's method typically **overestimates** the actual failure load
      by 20--40 %%.  Apply a reduction factor (commonly 0.80) for
      design.
    * Only data points where Q > 0 are used.

    References
    ----------
    Chin (1970); Fellenius (2020), Section 8.9.

    Examples
    --------
    >>> from geoeq.site.pile_load import chin
    >>> Q = [100, 200, 400, 600, 800, 1000, 1200]
    >>> s = [0.5, 1.2, 3.0, 5.5, 9.0, 14.0, 21.0]
    >>> res = chin(Q, s)
    >>> res['Qu_chin'] > 0
    True
    """
    Q = np.asarray(load, dtype=float)
    s = np.asarray(settlement, dtype=float)
    if len(Q) != len(s):
        raise ValueError("load and settlement must have the same length.")

    # Filter Q > 0
    mask = Q > 0
    Q = Q[mask]
    s = s[mask]
    if len(Q) < 3:
        raise ValueError("Need at least 3 data points with Q > 0.")

    # Select the linear portion
    n = len(Q)
    i0 = start_index if start_index is not None else n // 2
    i1 = end_index if end_index is not None else n
    if i0 < 0 or i1 > n or i0 >= i1:
        raise ValueError("Invalid start_index / end_index range.")

    delta_over_Q = s / Q  # mm / kN  →  units of 1/kN * mm

    # Linear regression: delta/Q = C1 * delta + C2
    x = s[i0:i1]
    y = delta_over_Q[i0:i1]

    coeffs = np.polyfit(x, y, 1)
    C1 = coeffs[0]  # slope  (1/kN)
    C2 = coeffs[1]  # intercept

    # R-squared
    y_pred = np.polyval(coeffs, x)
    ss_res = np.sum((y - y_pred) ** 2)
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r_sq = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0

    Qu = 1.0 / C1 if C1 > 0 else float("nan")

    return {
        "Qu_chin": float(Qu),
        "slope": float(C1),
        "intercept": float(C2),
        "r_squared": float(r_sq),
    }

de_beer

Python

de_beer(load: Union[List[float], ndarray], settlement: Union[List[float], ndarray]) -> Dict[str, float]

De Beer's Double-Logarithmic Method (1967).

Plot :math:\log_{10}(Q) versus :math:\log_{10}(\Delta). The data form two approximately straight-line segments. The failure load is the load at the intersection of these two segments.

The method fits two straight lines to the first and second halves of the double-log data and finds their intersection.

PARAMETER	DESCRIPTION
load	Applied loads (kN). Must be > 0. TYPE: array_like
settlement	Corresponding settlements (mm). Must be > 0. TYPE: array_like

RETURNS	DESCRIPTION
dict	'Qu_de_beer' : float Failure load (kN) at the intersection. 'settlement_at_failure' : float Settlement at failure (mm). 'log_load_at_failure' : float log10(Qu). 'log_settlement_at_failure' : float log10(settlement) at failure.

Notes

• All data points must have Q > 0 and settlement > 0.
• The method uses an automated split: it searches for the break point that minimises total sum-of-squared residuals from two piecewise-linear fits in log-log space.

References

De Beer (1967); Fellenius (2020), Section 8.8.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import de_beer
>>> Q = [50, 100, 200, 400, 600, 800, 1000, 1200]
>>> s = [0.3, 0.7, 1.5, 3.5, 7.0, 14.0, 25.0, 42.0]
>>> res = de_beer(Q, s)
>>> res['Qu_de_beer'] > 0
True

Source code in geoeq/site/pile_load.py

def de_beer(
    load: Union[List[float], np.ndarray],
    settlement: Union[List[float], np.ndarray],
) -> Dict[str, float]:
    r"""
    De Beer's Double-Logarithmic Method (1967).

    Plot :math:`\log_{10}(Q)` versus :math:`\log_{10}(\Delta)`.  The
    data form two approximately straight-line segments.  The failure
    load is the load at the **intersection** of these two segments.

    The method fits two straight lines to the first and second halves
    of the double-log data and finds their intersection.

    Parameters
    ----------
    load : array_like
        Applied loads (kN).  Must be > 0.
    settlement : array_like
        Corresponding settlements (mm).  Must be > 0.

    Returns
    -------
    dict
        ``'Qu_de_beer'`` : float
            Failure load (kN) at the intersection.
        ``'settlement_at_failure'`` : float
            Settlement at failure (mm).
        ``'log_load_at_failure'`` : float
            log10(Qu).
        ``'log_settlement_at_failure'`` : float
            log10(settlement) at failure.

    Notes
    -----
    * All data points must have Q > 0 and settlement > 0.
    * The method uses an automated split: it searches for the break
      point that minimises total sum-of-squared residuals from two
      piecewise-linear fits in log-log space.

    References
    ----------
    De Beer (1967); Fellenius (2020), Section 8.8.

    Examples
    --------
    >>> from geoeq.site.pile_load import de_beer
    >>> Q = [50, 100, 200, 400, 600, 800, 1000, 1200]
    >>> s = [0.3, 0.7, 1.5, 3.5, 7.0, 14.0, 25.0, 42.0]
    >>> res = de_beer(Q, s)
    >>> res['Qu_de_beer'] > 0
    True
    """
    Q = np.asarray(load, dtype=float)
    s = np.asarray(settlement, dtype=float)
    if len(Q) != len(s):
        raise ValueError("load and settlement must have the same length.")
    if np.any(Q <= 0) or np.any(s <= 0):
        raise ValueError("All load and settlement values must be > 0.")
    if len(Q) < 4:
        raise ValueError("Need at least 4 data points.")

    logQ = np.log10(Q)
    logS = np.log10(s)

    # Find optimal break point by minimising combined residuals
    n = len(Q)
    best_ssr = float("inf")
    best_k = 2

    for k in range(2, n - 1):
        # Fit line to first segment [0..k] and second segment [k..n]
        c1 = np.polyfit(logS[:k + 1], logQ[:k + 1], 1)
        c2 = np.polyfit(logS[k:], logQ[k:], 1)

        r1 = logQ[:k + 1] - np.polyval(c1, logS[:k + 1])
        r2 = logQ[k:] - np.polyval(c2, logS[k:])
        ssr = np.sum(r1 ** 2) + np.sum(r2 ** 2)

        if ssr < best_ssr:
            best_ssr = ssr
            best_k = k

    # Final fits at best break point
    c1 = np.polyfit(logS[:best_k + 1], logQ[:best_k + 1], 1)
    c2 = np.polyfit(logS[best_k:], logQ[best_k:], 1)

    # Intersection: c1[0]*x + c1[1] = c2[0]*x + c2[1]
    if abs(c1[0] - c2[0]) < 1e-12:
        # Parallel lines — no intersection
        return {
            "Qu_de_beer": float("nan"),
            "settlement_at_failure": float("nan"),
            "log_load_at_failure": float("nan"),
            "log_settlement_at_failure": float("nan"),
        }

    logS_int = (c2[1] - c1[1]) / (c1[0] - c2[0])
    logQ_int = c1[0] * logS_int + c1[1]

    return {
        "Qu_de_beer": float(10 ** logQ_int),
        "settlement_at_failure": float(10 ** logS_int),
        "log_load_at_failure": float(logQ_int),
        "log_settlement_at_failure": float(logS_int),
    }

hansen_80

Python

hansen_80(load: Union[List[float], ndarray], settlement: Union[List[float], ndarray], start_index: Optional[int] = None, end_index: Optional[int] = None) -> Dict[str, float]

Hansen's 80 %% Method (1963).

Plot :math:\sqrt{\Delta / Q} versus :math:\Delta. At large settlements the relationship becomes linear:

.. math::

Text Only

\sqrt{\frac{\Delta}{Q}} = a \, \Delta + b

The ultimate capacity is:

.. math::

Text Only

Q_u = \frac{1}{2 \sqrt{a \, b}}

and the settlement at failure is:

.. math::

Text Only

\Delta_u = \frac{b}{a}

The name "80 %% method" comes from the fact that at the ultimate load the pile mobilises 80 %% of its capacity for any settlement beyond :math:\Delta_u.

PARAMETER	DESCRIPTION
load	Applied loads (kN). Must be > 0. TYPE: array_like
settlement	Corresponding settlements (mm). Must be > 0. TYPE: array_like
start_index	Start of linear portion (0-based). Default: last 50 %%. TYPE: int DEFAULT: None
end_index	End of linear portion (exclusive). TYPE: int DEFAULT: None

RETURNS	DESCRIPTION
dict	'Qu_hansen' : float Hansen ultimate load (kN). 'settlement_at_failure' : float Settlement at failure (mm) = b/a. 'a' : float Slope of the linear fit. 'b' : float Intercept of the linear fit. 'r_squared' : float Coefficient of determination.

Notes

• Like Chin's method, Hansen's method also tends to overestimate the actual failure load.
• If a or b is negative the result is physically meaningless and nan is returned.

References

Hansen (1963); Fellenius (2020), Section 8.10.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import hansen_80
>>> Q = [100, 200, 400, 600, 800, 1000, 1200]
>>> s = [0.5, 1.2, 3.0, 5.5, 9.0, 14.0, 21.0]
>>> res = hansen_80(Q, s)
>>> res['Qu_hansen'] > 0
True

Source code in geoeq/site/pile_load.py

def hansen_80(
    load: Union[List[float], np.ndarray],
    settlement: Union[List[float], np.ndarray],
    start_index: Optional[int] = None,
    end_index: Optional[int] = None,
) -> Dict[str, float]:
    r"""
    Hansen's 80 %% Method (1963).

    Plot :math:`\sqrt{\Delta / Q}` versus :math:`\Delta`.  At large
    settlements the relationship becomes linear:

    .. math::

        \sqrt{\frac{\Delta}{Q}} = a \, \Delta + b

    The ultimate capacity is:

    .. math::

        Q_u = \frac{1}{2 \sqrt{a \, b}}

    and the settlement at failure is:

    .. math::

        \Delta_u = \frac{b}{a}

    The name "80 %% method" comes from the fact that at the ultimate
    load the pile mobilises 80 %% of its capacity for any
    settlement beyond :math:`\Delta_u`.

    Parameters
    ----------
    load : array_like
        Applied loads (kN).  Must be > 0.
    settlement : array_like
        Corresponding settlements (mm).  Must be > 0.
    start_index : int, optional
        Start of linear portion (0-based).  Default: last 50 %%.
    end_index : int, optional
        End of linear portion (exclusive).

    Returns
    -------
    dict
        ``'Qu_hansen'`` : float
            Hansen ultimate load (kN).
        ``'settlement_at_failure'`` : float
            Settlement at failure (mm) = b/a.
        ``'a'`` : float
            Slope of the linear fit.
        ``'b'`` : float
            Intercept of the linear fit.
        ``'r_squared'`` : float
            Coefficient of determination.

    Notes
    -----
    * Like Chin's method, Hansen's method also tends to **overestimate**
      the actual failure load.
    * If a or b is negative the result is physically meaningless and
      ``nan`` is returned.

    References
    ----------
    Hansen (1963); Fellenius (2020), Section 8.10.

    Examples
    --------
    >>> from geoeq.site.pile_load import hansen_80
    >>> Q = [100, 200, 400, 600, 800, 1000, 1200]
    >>> s = [0.5, 1.2, 3.0, 5.5, 9.0, 14.0, 21.0]
    >>> res = hansen_80(Q, s)
    >>> res['Qu_hansen'] > 0
    True
    """
    Q = np.asarray(load, dtype=float)
    s = np.asarray(settlement, dtype=float)
    if len(Q) != len(s):
        raise ValueError("load and settlement must have the same length.")

    mask = (Q > 0) & (s > 0)
    Q = Q[mask]
    s = s[mask]
    if len(Q) < 3:
        raise ValueError("Need at least 3 data points with Q > 0 and s > 0.")

    n = len(Q)
    i0 = start_index if start_index is not None else n // 2
    i1 = end_index if end_index is not None else n
    if i0 < 0 or i1 > n or i0 >= i1:
        raise ValueError("Invalid start_index / end_index range.")

    y = np.sqrt(s / Q)  # sqrt(mm / kN)

    x_fit = s[i0:i1]
    y_fit = y[i0:i1]

    coeffs = np.polyfit(x_fit, y_fit, 1)
    a = coeffs[0]  # slope
    b = coeffs[1]  # intercept

    # R-squared
    y_pred = np.polyval(coeffs, x_fit)
    ss_res = np.sum((y_fit - y_pred) ** 2)
    ss_tot = np.sum((y_fit - np.mean(y_fit)) ** 2)
    r_sq = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0

    if a > 0 and b > 0:
        Qu = 1.0 / (2.0 * np.sqrt(a * b))
        s_fail = b / a
    else:
        Qu = float("nan")
        s_fail = float("nan")

    return {
        "Qu_hansen": float(Qu),
        "settlement_at_failure": float(s_fail),
        "a": float(a),
        "b": float(b),
        "r_squared": float(r_sq),
    }

fhwa_5_percent

Python

fhwa_5_percent(load: Union[List[float], ndarray], settlement: Union[List[float], ndarray], diameter: float) -> Dict[str, float]

FHWA / Modified Davisson 5 %% Diameter Criterion for drilled shafts.

The failure load is the load corresponding to a settlement equal to 5 %% of the pile (shaft) diameter:

.. math::

Text Only

\Delta_{\text{fail}} = 0.05 \, D

This criterion is recommended by FHWA for drilled shafts and large-diameter bored piles (> 600 mm / 24 in.).

PARAMETER	DESCRIPTION
load	Applied loads (kN). TYPE: array_like
settlement	Corresponding settlements (mm). TYPE: array_like
diameter	Shaft diameter (mm). TYPE: float

RETURNS	DESCRIPTION
dict	'Qu_fhwa' : float Failure load (kN) at 5 %% diameter settlement. 'settlement_criterion' : float 5 %% of diameter (mm).

References

FHWA (2010). Drilled Shafts: Construction Procedures and LRFD Design Methods, FHWA-NHI-10-016, Section 13.3.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import fhwa_5_percent
>>> Q = [0, 500, 1000, 1500, 2000, 2500, 3000]
>>> s = [0, 2.0, 5.0, 12.0, 25.0, 45.0, 75.0]
>>> res = fhwa_5_percent(Q, s, diameter=600)
>>> res['settlement_criterion']
30.0

Source code in geoeq/site/pile_load.py

def fhwa_5_percent(
    load: Union[List[float], np.ndarray],
    settlement: Union[List[float], np.ndarray],
    diameter: float,
) -> Dict[str, float]:
    r"""
    FHWA / Modified Davisson 5 %% Diameter Criterion for drilled shafts.

    The failure load is the load corresponding to a settlement equal to
    5 %% of the pile (shaft) diameter:

    .. math::

        \Delta_{\text{fail}} = 0.05 \, D

    This criterion is recommended by FHWA for drilled shafts and
    large-diameter bored piles (> 600 mm / 24 in.).

    Parameters
    ----------
    load : array_like
        Applied loads (kN).
    settlement : array_like
        Corresponding settlements (mm).
    diameter : float
        Shaft diameter (mm).

    Returns
    -------
    dict
        ``'Qu_fhwa'`` : float
            Failure load (kN) at 5 %% diameter settlement.
        ``'settlement_criterion'`` : float
            5 %% of diameter (mm).

    References
    ----------
    FHWA (2010). *Drilled Shafts: Construction Procedures and LRFD
    Design Methods*, FHWA-NHI-10-016, Section 13.3.

    Examples
    --------
    >>> from geoeq.site.pile_load import fhwa_5_percent
    >>> Q = [0, 500, 1000, 1500, 2000, 2500, 3000]
    >>> s = [0, 2.0, 5.0, 12.0, 25.0, 45.0, 75.0]
    >>> res = fhwa_5_percent(Q, s, diameter=600)
    >>> res['settlement_criterion']
    30.0
    """
    Q = np.asarray(load, dtype=float)
    s = np.asarray(settlement, dtype=float)
    check_positive(diameter, "diameter")
    if len(Q) != len(s):
        raise ValueError("load and settlement must have the same length.")

    s_crit = 0.05 * diameter  # mm

    idx = np.argsort(s)
    Q_sorted = Q[idx]
    s_sorted = s[idx]

    if s_crit > s_sorted[-1]:
        Qu = float("nan")
    elif s_crit <= s_sorted[0]:
        Qu = float(Q_sorted[0])
    else:
        from scipy.interpolate import interp1d
        interp = interp1d(s_sorted, Q_sorted, kind="linear")
        Qu = float(interp(s_crit))

    return {
        "Qu_fhwa": Qu,
        "settlement_criterion": float(s_crit),
    }

case_method

Python

case_method(F1: float, F2: float, v1: float, v2: float, impedance: float, Jc: float = 0.5) -> Dict[str, float]

Case Method for dynamic pile capacity (Goble et al. 1975).

From force and velocity measurements at the pile head at times :math:t_1 (impact peak) and :math:t_2 = t_1 + 2L/c:

.. math::

Text Only

R_{\text{total}} = \frac{1}{2}\bigl[(F_1 + F_2) +
    Z\,(v_1 - v_2)\bigr]

.. math::

Text Only

R_{\text{static}} = R_{\text{total}}\,(1 - J_c) +
    \frac{1}{2}\bigl[(F_2 + Z\,v_2)\bigr]\,J_c

Or equivalently (closed-form used in practice):

.. math::

Text Only

R_{\text{static}} = \frac{(1-J_c)}{2}(F_1 + Z\,v_1)
    + \frac{(1+J_c)}{2}(F_2 - Z\,v_2) \; \times (-1)

The simplified FHWA form used here:

.. math::

Text Only

R_{\text{static}} = \frac{1}{2}\bigl[(1-J_c)(F_1+Z\,v_1)
    + (1+J_c)(F_2-Z\,v_2)\bigr] \times \frac{1}{1}

Note: the sign convention uses compressive force positive, velocity positive downward. See Rausche et al. (1985) for derivation.

PARAMETER	DESCRIPTION
F1	Force at pile head at time t1 — first peak (kN). TYPE: float
F2	Force at pile head at time t2 = t1 + 2L/c (kN). TYPE: float
v1	Velocity at pile head at time t1 (m/s). TYPE: float
v2	Velocity at pile head at time t2 (m/s). TYPE: float
impedance	Pile impedance Z = EA/c (kNs/m). Compute with :func:pile_impedance. TYPE: float
Jc	Case damping factor (dimensionless). Typical ranges by soil type: ======== ============ Soil Jc ======== ============ Sand 0.10 -- 0.20 Silt 0.20 -- 0.40 Clay 0.40 -- 0.70 Peat 0.70 -- 1.00 ======== ============ TYPE: float DEFAULT: 0.5

RETURNS	DESCRIPTION
dict	'R_total' : float — Total dynamic resistance (kN). 'R_static' : float — Estimated static resistance (kN). 'Jc' : float — Damping factor used.

References

Goble, Rausche & Likins (1975); Rausche, F., Goble, G. G., & Likins, G. E. (1985). "Dynamic determination of pile capacity." ASCE J. Geotech. Eng., 111(3).

Examples:

Python Console Session

>>> from geoeq.site.pile_load import case_method
>>> res = case_method(F1=2500, F2=800, v1=3.0, v2=0.5,
...                   impedance=4000, Jc=0.4)
>>> res['R_static'] > 0
True

Source code in geoeq/site/pile_load.py

def case_method(
    F1: float,
    F2: float,
    v1: float,
    v2: float,
    impedance: float,
    Jc: float = 0.5,
) -> Dict[str, float]:
    r"""
    Case Method for dynamic pile capacity (Goble et al. 1975).

    From force and velocity measurements at the pile head at times
    :math:`t_1` (impact peak) and :math:`t_2 = t_1 + 2L/c`:

    .. math::

        R_{\text{total}} = \frac{1}{2}\bigl[(F_1 + F_2) +
            Z\,(v_1 - v_2)\bigr]

    .. math::

        R_{\text{static}} = R_{\text{total}}\,(1 - J_c) +
            \frac{1}{2}\bigl[(F_2 + Z\,v_2)\bigr]\,J_c

    Or equivalently (closed-form used in practice):

    .. math::

        R_{\text{static}} = \frac{(1-J_c)}{2}(F_1 + Z\,v_1)
            + \frac{(1+J_c)}{2}(F_2 - Z\,v_2) \; \times (-1)

    The simplified FHWA form used here:

    .. math::

        R_{\text{static}} = \frac{1}{2}\bigl[(1-J_c)(F_1+Z\,v_1)
            + (1+J_c)(F_2-Z\,v_2)\bigr] \times \frac{1}{1}

    Note: the sign convention uses compressive force positive, velocity
    positive downward.  See Rausche et al. (1985) for derivation.

    Parameters
    ----------
    F1 : float
        Force at pile head at time t1 — first peak (kN).
    F2 : float
        Force at pile head at time t2 = t1 + 2L/c (kN).
    v1 : float
        Velocity at pile head at time t1 (m/s).
    v2 : float
        Velocity at pile head at time t2 (m/s).
    impedance : float
        Pile impedance Z = E*A/c (kN*s/m).
        Compute with :func:`pile_impedance`.
    Jc : float, default 0.5
        Case damping factor (dimensionless).

        Typical ranges by soil type:

        ======== ============
        Soil     Jc
        ======== ============
        Sand     0.10 -- 0.20
        Silt     0.20 -- 0.40
        Clay     0.40 -- 0.70
        Peat     0.70 -- 1.00
        ======== ============

    Returns
    -------
    dict
        ``'R_total'`` : float — Total dynamic resistance (kN).
        ``'R_static'`` : float — Estimated static resistance (kN).
        ``'Jc'`` : float — Damping factor used.

    References
    ----------
    Goble, Rausche & Likins (1975);
    Rausche, F., Goble, G. G., & Likins, G. E. (1985). "Dynamic
    determination of pile capacity." *ASCE J. Geotech. Eng.*, 111(3).

    Examples
    --------
    >>> from geoeq.site.pile_load import case_method
    >>> res = case_method(F1=2500, F2=800, v1=3.0, v2=0.5,
    ...                   impedance=4000, Jc=0.4)
    >>> res['R_static'] > 0
    True
    """
    check_positive(impedance, "impedance")
    check_range(Jc, "Jc", 0.0, 1.0)

    Z = impedance

    R_total = 0.5 * ((F1 + F2) + Z * (v1 - v2))

    # Standard FHWA Case Method (GRL formulation)
    R_static = 0.5 * ((1 - Jc) * (F1 + Z * v1) + (1 + Jc) * (F2 - Z * v2))

    return {
        "R_total": float(R_total),
        "R_static": float(R_static),
        "Jc": float(Jc),
    }

hiley

Python

hiley(hammer_weight: float, drop_height: float, efficiency: float, set_per_blow: float, elastic_compression: float) -> Dict[str, float]

Hiley Driving Formula — ultimate pile capacity from driving data.

.. math::

Text Only

Q_u = \frac{W \, h \, \eta}{s + c/2}

PARAMETER	DESCRIPTION
hammer_weight	Weight of the hammer, W (kN). TYPE: float
drop_height	Drop height, h (m). TYPE: float
efficiency	Hammer efficiency, :math:\eta (dimensionless, 0 to 1). Typical values: ======================== =========== Hammer type Efficiency ======================== =========== Drop hammer (free fall) 0.75 -- 1.0 Single-acting air/steam 0.75 -- 0.85 Double-acting air/steam 0.70 -- 0.80 Diesel hammer 0.70 -- 0.90 Hydraulic hammer 0.80 -- 0.95 ======================== =========== TYPE: float
set_per_blow	Permanent set per blow, s (m). Typically measured as the average of the last 10 blows. s > 0. TYPE: float
elastic_compression	Total temporary (elastic) compression, c (m). Sum of pile, soil, and driving cap compressions. Typical range: 0.002 -- 0.025 m. TYPE: float

RETURNS	DESCRIPTION
dict	'Qu_hiley' : float — Ultimate pile capacity (kN). 'energy_input' : float — Hammer energy Wheta (kN*m).

References

Hiley, A. (1925). "A rational pile driving formula and its application in practice explained." Engineering (London), 119.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import hiley
>>> res = hiley(hammer_weight=50, drop_height=1.5, efficiency=0.85,
...            set_per_blow=0.010, elastic_compression=0.005)
>>> res['Qu_hiley'] > 0
True

Source code in geoeq/site/pile_load.py

def hiley(
    hammer_weight: float,
    drop_height: float,
    efficiency: float,
    set_per_blow: float,
    elastic_compression: float,
) -> Dict[str, float]:
    r"""
    Hiley Driving Formula — ultimate pile capacity from driving data.

    .. math::

        Q_u = \frac{W \, h \, \eta}{s + c/2}

    Parameters
    ----------
    hammer_weight : float
        Weight of the hammer, W (kN).
    drop_height : float
        Drop height, h (m).
    efficiency : float
        Hammer efficiency, :math:`\eta` (dimensionless, 0 to 1).

        Typical values:

        ======================== ===========
        Hammer type              Efficiency
        ======================== ===========
        Drop hammer (free fall)  0.75 -- 1.0
        Single-acting air/steam  0.75 -- 0.85
        Double-acting air/steam  0.70 -- 0.80
        Diesel hammer            0.70 -- 0.90
        Hydraulic hammer         0.80 -- 0.95
        ======================== ===========

    set_per_blow : float
        Permanent set per blow, s (m).  Typically measured as the
        average of the last 10 blows.  s > 0.
    elastic_compression : float
        Total temporary (elastic) compression, c (m).
        Sum of pile, soil, and driving cap compressions.
        Typical range: 0.002 -- 0.025 m.

    Returns
    -------
    dict
        ``'Qu_hiley'`` : float — Ultimate pile capacity (kN).
        ``'energy_input'`` : float — Hammer energy W*h*eta (kN*m).

    References
    ----------
    Hiley, A. (1925). "A rational pile driving formula and its
    application in practice explained." *Engineering (London)*, 119.

    Examples
    --------
    >>> from geoeq.site.pile_load import hiley
    >>> res = hiley(hammer_weight=50, drop_height=1.5, efficiency=0.85,
    ...            set_per_blow=0.010, elastic_compression=0.005)
    >>> res['Qu_hiley'] > 0
    True
    """
    check_positive(hammer_weight, "hammer_weight")
    check_positive(drop_height, "drop_height")
    check_range(efficiency, "efficiency", 0.0, 1.0)
    check_positive(set_per_blow, "set_per_blow")
    check_non_negative(elastic_compression, "elastic_compression")

    W = hammer_weight
    h = drop_height
    eta = efficiency
    s = set_per_blow
    c = elastic_compression

    energy = W * h * eta
    Qu = energy / (s + c / 2.0)

    return {
        "Qu_hiley": float(Qu),
        "energy_input": float(energy),
    }

danish_formula

Python

danish_formula(hammer_weight: float, drop_height: float, efficiency: float, set_per_blow: float, pile_length: float, pile_area: float, pile_modulus: float) -> Dict[str, float]

Danish Driving Formula — ultimate pile capacity.

.. math::

Text Only

Q_u = \frac{e_h \, W \, h}
{s + \sqrt{\frac{e_h \, W \, h \, L}{2 \, A \, E}}}

where :math:e_h = \eta is the hammer efficiency.

PARAMETER	DESCRIPTION
hammer_weight	Weight of the hammer, W (kN). TYPE: float
drop_height	Drop height, h (m). TYPE: float
efficiency	Hammer efficiency, :math:e_h (dimensionless, 0 to 1). TYPE: float
set_per_blow	Permanent set per blow, s (m). s > 0. TYPE: float
pile_length	Pile length, L (m). TYPE: float
pile_area	Cross-sectional area, A (m^2). TYPE: float
pile_modulus	Elastic modulus, E (kN/m^2 = kPa). Typical: steel ~200e6 kPa, concrete ~25e6 kPa. TYPE: float

RETURNS	DESCRIPTION
dict	'Qu_danish' : float — Ultimate pile capacity (kN). 'elastic_compression_term' : float The square-root term (m).

References

Sorensen, T. & Hansen, B. (1957). "Pile driving formulae — an investigation based on dimensional analysis and statistical analysis." Proc., 4^th ICSMFE, London, Vol. 2, pp. 61--65.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import danish_formula
>>> res = danish_formula(hammer_weight=50, drop_height=1.5,
...     efficiency=0.85, set_per_blow=0.010,
...     pile_length=15, pile_area=0.07, pile_modulus=25e6)
>>> res['Qu_danish'] > 0
True

Source code in geoeq/site/pile_load.py

def danish_formula(
    hammer_weight: float,
    drop_height: float,
    efficiency: float,
    set_per_blow: float,
    pile_length: float,
    pile_area: float,
    pile_modulus: float,
) -> Dict[str, float]:
    r"""
    Danish Driving Formula — ultimate pile capacity.

    .. math::

        Q_u = \frac{e_h \, W \, h}
        {s + \sqrt{\frac{e_h \, W \, h \, L}{2 \, A \, E}}}

    where :math:`e_h = \eta` is the hammer efficiency.

    Parameters
    ----------
    hammer_weight : float
        Weight of the hammer, W (kN).
    drop_height : float
        Drop height, h (m).
    efficiency : float
        Hammer efficiency, :math:`e_h` (dimensionless, 0 to 1).
    set_per_blow : float
        Permanent set per blow, s (m).  s > 0.
    pile_length : float
        Pile length, L (m).
    pile_area : float
        Cross-sectional area, A (m^2).
    pile_modulus : float
        Elastic modulus, E (kN/m^2 = kPa).
        Typical: steel ~200e6 kPa, concrete ~25e6 kPa.

    Returns
    -------
    dict
        ``'Qu_danish'`` : float — Ultimate pile capacity (kN).
        ``'elastic_compression_term'`` : float
            The square-root term (m).

    References
    ----------
    Sorensen, T. & Hansen, B. (1957). "Pile driving formulae — an
    investigation based on dimensional analysis and statistical
    analysis." *Proc., 4th ICSMFE*, London, Vol. 2, pp. 61--65.

    Examples
    --------
    >>> from geoeq.site.pile_load import danish_formula
    >>> res = danish_formula(hammer_weight=50, drop_height=1.5,
    ...     efficiency=0.85, set_per_blow=0.010,
    ...     pile_length=15, pile_area=0.07, pile_modulus=25e6)
    >>> res['Qu_danish'] > 0
    True
    """
    check_positive(hammer_weight, "hammer_weight")
    check_positive(drop_height, "drop_height")
    check_range(efficiency, "efficiency", 0.0, 1.0)
    check_positive(set_per_blow, "set_per_blow")
    check_positive(pile_length, "pile_length")
    check_positive(pile_area, "pile_area")
    check_positive(pile_modulus, "pile_modulus")

    W = hammer_weight
    h = drop_height
    eh = efficiency
    s = set_per_blow
    L = pile_length
    A = pile_area
    E = pile_modulus

    energy = eh * W * h
    elastic_term = np.sqrt(energy * L / (2.0 * A * E))
    Qu = energy / (s + elastic_term)

    return {
        "Qu_danish": float(Qu),
        "elastic_compression_term": float(elastic_term),
    }

enr

Python

enr(hammer_weight: float, drop_height: float, set_per_blow: float, rebound: float = 0.0254, factor_of_safety: float = 6.0) -> Dict[str, float]

Engineering News Record (ENR) Formula — allowable pile capacity.

.. math::

Text Only

Q_a = \frac{W \, h}{F_s \, (s + c)}

PARAMETER	DESCRIPTION
hammer_weight	Weight of the hammer, W (kN). TYPE: float
drop_height	Drop height, h (m). TYPE: float
set_per_blow	Permanent set per blow, s (m). TYPE: float
rebound	Elastic rebound / temporary compression, c (m). Original ENR uses c = 25.4 mm (1 inch) for drop hammers, c = 2.54 mm (0.1 inch) for steam/air hammers. TYPE: float DEFAULT: 0.0254
factor_of_safety	Built-in factor of safety. The original ENR formula uses FS = 6. The Modified ENR uses a variable FS. TYPE: float DEFAULT: 6.0

RETURNS	DESCRIPTION
dict	'Qa_enr' : float — Allowable pile capacity (kN). 'Qu_enr' : float — Implied ultimate capacity = Qa * FS (kN).

Notes

• The ENR formula is considered unreliable by modern standards. It has a very high coefficient of variation (> 0.5) and should only be used for preliminary estimates or where required by local building codes.
• ASCE and FHWA discourage its use for design.

References

Wellington, A. M. (1893). "Piles and pile driving." Engineering News Record. Das (2021), Section 11.16.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import enr
>>> res = enr(hammer_weight=30, drop_height=1.0, set_per_blow=0.010)
>>> res['Qa_enr'] > 0
True

Source code in geoeq/site/pile_load.py

def enr(
    hammer_weight: float,
    drop_height: float,
    set_per_blow: float,
    rebound: float = 0.0254,
    factor_of_safety: float = 6.0,
) -> Dict[str, float]:
    r"""
    Engineering News Record (ENR) Formula — allowable pile capacity.

    .. math::

        Q_a = \frac{W \, h}{F_s \, (s + c)}

    Parameters
    ----------
    hammer_weight : float
        Weight of the hammer, W (kN).
    drop_height : float
        Drop height, h (m).
    set_per_blow : float
        Permanent set per blow, s (m).
    rebound : float, default 0.0254
        Elastic rebound / temporary compression, c (m).
        Original ENR uses c = 25.4 mm (1 inch) for drop hammers,
        c = 2.54 mm (0.1 inch) for steam/air hammers.
    factor_of_safety : float, default 6.0
        Built-in factor of safety.  The original ENR formula uses
        FS = 6.  The **Modified ENR** uses a variable FS.

    Returns
    -------
    dict
        ``'Qa_enr'`` : float — Allowable pile capacity (kN).
        ``'Qu_enr'`` : float — Implied ultimate capacity = Qa * FS (kN).

    Notes
    -----
    * The ENR formula is considered **unreliable** by modern standards.
      It has a very high coefficient of variation (> 0.5) and should
      only be used for preliminary estimates or where required by
      local building codes.
    * ASCE and FHWA discourage its use for design.

    References
    ----------
    Wellington, A. M. (1893). "Piles and pile driving."
    *Engineering News Record*.
    Das (2021), Section 11.16.

    Examples
    --------
    >>> from geoeq.site.pile_load import enr
    >>> res = enr(hammer_weight=30, drop_height=1.0, set_per_blow=0.010)
    >>> res['Qa_enr'] > 0
    True
    """
    check_positive(hammer_weight, "hammer_weight")
    check_positive(drop_height, "drop_height")
    check_positive(set_per_blow, "set_per_blow")
    check_non_negative(rebound, "rebound")
    check_positive(factor_of_safety, "factor_of_safety")

    W = hammer_weight
    h = drop_height
    s = set_per_blow
    c = rebound
    FS = factor_of_safety

    Qa = (W * h) / (FS * (s + c))
    Qu = Qa * FS

    return {
        "Qa_enr": float(Qa),
        "Qu_enr": float(Qu),
    }

pile_impedance

Python

pile_impedance(elastic_modulus: float, area: float, wave_speed: float) -> Dict[str, float]

Pile impedance — fundamental property for wave-equation analysis.

.. math::

Text Only

Z = \frac{E \, A}{c}

PARAMETER	DESCRIPTION
elastic_modulus	Elastic modulus, E (kN/m^2 = kPa). Typical: steel ~200e6 kPa, concrete ~30e6 kPa. TYPE: float
area	Cross-sectional area, A (m^2). TYPE: float
wave_speed	Stress-wave speed, c (m/s). Typical: steel ~5120 m/s, concrete ~3600--4000 m/s. TYPE: float

RETURNS	DESCRIPTION
dict	'impedance' : float — Impedance Z (kN*s/m).

Notes

• Impedance changes along the pile indicate defects:

• Decrease in Z → necking, crack, or void.

• Increase in Z → bulging or inclusion of harder material.

• Related: :math:c = \sqrt{E / \rho}, where :math:\rho is the mass density (kg/m^3).

References

Rausche et al. (1985); ASTM D4945 (PDA testing).

Examples:

Python Console Session

>>> from geoeq.site.pile_load import pile_impedance
>>> res = pile_impedance(elastic_modulus=30e6, area=0.07,
...                      wave_speed=3800)
>>> round(res['impedance'], 1)
552.6

Source code in geoeq/site/pile_load.py

def pile_impedance(
    elastic_modulus: float,
    area: float,
    wave_speed: float,
) -> Dict[str, float]:
    r"""
    Pile impedance — fundamental property for wave-equation analysis.

    .. math::

        Z = \frac{E \, A}{c}

    Parameters
    ----------
    elastic_modulus : float
        Elastic modulus, E (kN/m^2 = kPa).
        Typical: steel ~200e6 kPa, concrete ~30e6 kPa.
    area : float
        Cross-sectional area, A (m^2).
    wave_speed : float
        Stress-wave speed, c (m/s).
        Typical: steel ~5120 m/s, concrete ~3600--4000 m/s.

    Returns
    -------
    dict
        ``'impedance'`` : float — Impedance Z (kN*s/m).

    Notes
    -----
    * Impedance changes along the pile indicate defects:

      - **Decrease** in Z → necking, crack, or void.
      - **Increase** in Z → bulging or inclusion of harder material.

    * Related: :math:`c = \sqrt{E / \rho}`, where :math:`\rho` is
      the mass density (kg/m^3).

    References
    ----------
    Rausche et al. (1985); ASTM D4945 (PDA testing).

    Examples
    --------
    >>> from geoeq.site.pile_load import pile_impedance
    >>> res = pile_impedance(elastic_modulus=30e6, area=0.07,
    ...                      wave_speed=3800)
    >>> round(res['impedance'], 1)
    552.6
    """
    check_positive(elastic_modulus, "elastic_modulus")
    check_positive(area, "area")
    check_positive(wave_speed, "wave_speed")

    Z = elastic_modulus * area / wave_speed

    return {"impedance": float(Z)}

beta_integrity

Python

beta_integrity(measured_length: float, known_length: float) -> Dict[str, float]

Beta integrity factor — pile length verification from sonic test.

.. math::

Text Only

\beta = \frac{L_{\text{measured}}}{L_{\text{known}}}

PARAMETER	DESCRIPTION
measured_length	Length obtained from low-strain integrity test (m). TYPE: float
known_length	As-built or design pile length (m). TYPE: float

RETURNS	DESCRIPTION
dict	'beta' : float — Ratio of measured to known length. 'interpretation' : str — Qualitative assessment.

Notes

Interpretation guide:

============= ========================================= Beta range Interpretation ============= ========================================= 0.95 -- 1.05 Intact — no significant defect detected 0.80 -- 0.95 Possible defect — further investigation < 0.80 Likely severe defect or shortened pile

1.05 Possible bulging or wave speed error ============= =========================================

References

ASTM D5882 (2016). Standard Test Method for Low Strain Integrity Testing of Deep Foundations.

Examples:

Python Console Session

>>> from geoeq.site.pile_load import beta_integrity
>>> res = beta_integrity(measured_length=14.5, known_length=15.0)
>>> round(res['beta'], 3)
0.967

Source code in geoeq/site/pile_load.py

def beta_integrity(
    measured_length: float,
    known_length: float,
) -> Dict[str, float]:
    r"""
    Beta integrity factor — pile length verification from sonic test.

    .. math::

        \beta = \frac{L_{\text{measured}}}{L_{\text{known}}}

    Parameters
    ----------
    measured_length : float
        Length obtained from low-strain integrity test (m).
    known_length : float
        As-built or design pile length (m).

    Returns
    -------
    dict
        ``'beta'`` : float — Ratio of measured to known length.
        ``'interpretation'`` : str — Qualitative assessment.

    Notes
    -----
    Interpretation guide:

    ============= =========================================
    Beta range    Interpretation
    ============= =========================================
    0.95 -- 1.05  Intact — no significant defect detected
    0.80 -- 0.95  Possible defect — further investigation
    < 0.80        Likely severe defect or shortened pile
    > 1.05        Possible bulging or wave speed error
    ============= =========================================

    References
    ----------
    ASTM D5882 (2016). Standard Test Method for Low Strain Integrity
    Testing of Deep Foundations.

    Examples
    --------
    >>> from geoeq.site.pile_load import beta_integrity
    >>> res = beta_integrity(measured_length=14.5, known_length=15.0)
    >>> round(res['beta'], 3)
    0.967
    """
    check_positive(measured_length, "measured_length")
    check_positive(known_length, "known_length")

    beta = measured_length / known_length

    if 0.95 <= beta <= 1.05:
        interp = "Intact — no significant defect detected."
    elif 0.80 <= beta < 0.95:
        interp = "Possible defect — further investigation recommended."
    elif beta < 0.80:
        interp = "Likely severe defect or shortened pile."
    else:
        interp = "Possible bulging or wave speed calibration error."

    return {
        "beta": float(beta),
        "interpretation": interp,
    }

Field permeability

slug_test

Python

slug_test(r: float, R: float, Le: float, T0: float, method: str = 'hvorslev') -> float

Hydraulic conductivity from slug test (Hvorslev method).

The Hvorslev (1951) basic time lag method for a well point or piezometer with length/diameter ratio > 8:

.. math::

Text Only

k = \frac{r^2 \ln(L_e / R)}{2\,L_e\,T_0}

PARAMETER	DESCRIPTION
r	Standpipe (casing) radius (m). TYPE: float
R	Well screen or intake radius (m). TYPE: float
Le	Length of the well screen or intake (m). TYPE: float
T0	Basic time lag (s) — time for head to reach 37% (1/e) of initial displacement. Obtained from the slope of ln(h/h₀) vs t plot. TYPE: float
method	Currently only 'hvorslev' is implemented. TYPE: str DEFAULT: ``'hvorslev'``

RETURNS	DESCRIPTION
float	Hydraulic conductivity k (m/s).

Notes

• Valid for Le/R > 8.
• For Le/R < 8, use the Bouwer & Rice (1976) method.
• T₀ is found from the -1/slope of the ln(h/h₀) vs time plot (the time at which ln(h/h₀) = -1).

References

Hvorslev (1951); Das (2021), Section 5.7.

Examples:

Python Console Session

>>> from geoeq.site.field_perm import slug_test
>>> k = slug_test(r=0.025, R=0.05, Le=1.5, T0=300)
>>> f'{k:.2e}'
'3.05e-07'

Source code in geoeq/site/field_perm.py

def slug_test(
    r: float,
    R: float,
    Le: float,
    T0: float,
    method: str = "hvorslev",
) -> float:
    r"""
    Hydraulic conductivity from slug test (Hvorslev method).

    The Hvorslev (1951) basic time lag method for a well point or
    piezometer with length/diameter ratio > 8:

    .. math::

        k = \frac{r^2 \ln(L_e / R)}{2\,L_e\,T_0}

    Parameters
    ----------
    r : float
        Standpipe (casing) radius (m).
    R : float
        Well screen or intake radius (m).
    Le : float
        Length of the well screen or intake (m).
    T0 : float
        Basic time lag (s) — time for head to reach 37% (1/e) of
        initial displacement.  Obtained from the slope of
        ln(h/h₀) vs t plot.
    method : str, default ``'hvorslev'``
        Currently only ``'hvorslev'`` is implemented.

    Returns
    -------
    float
        Hydraulic conductivity k (m/s).

    Notes
    -----
    * Valid for Le/R > 8.
    * For Le/R < 8, use the Bouwer & Rice (1976) method.
    * T₀ is found from the -1/slope of the ln(h/h₀) vs time plot
      (the time at which ln(h/h₀) = -1).

    References
    ----------
    Hvorslev (1951); Das (2021), Section 5.7.

    Examples
    --------
    >>> from geoeq.site.field_perm import slug_test
    >>> k = slug_test(r=0.025, R=0.05, Le=1.5, T0=300)
    >>> f'{k:.2e}'
    '3.05e-07'
    """
    check_positive(r, "r")
    check_positive(R, "R")
    check_positive(Le, "Le")
    check_positive(T0, "T0")

    method_l = method.lower()
    if method_l != "hvorslev":
        raise ValueError(f"Unknown method '{method}'. Currently only 'hvorslev' is supported.")

    k = r**2 * np.log(Le / R) / (2.0 * Le * T0)
    return float(k)

pumping_test_confined

Python

pumping_test_confined(Q: float, h1: float, h2: float, r1: float, r2: float, H: float) -> Dict[str, float]

Hydraulic conductivity from pumping test — confined aquifer.

Thiem (1906) equilibrium equation for steady-state flow to a well in a confined aquifer:

.. math::

Text Only

k = \frac{Q \ln(r_2 / r_1)}{2\,\pi\,H\,(h_2 - h_1)}

PARAMETER	DESCRIPTION
Q	Steady-state pumping rate (m³/s). TYPE: float
h1	Piezometric head at observation well 1 (m). TYPE: float
h2	Piezometric head at observation well 2 (m), h₂ > h₁. TYPE: float
r1	Radial distance of observation well 1 from pumped well (m). TYPE: float
r2	Radial distance of observation well 2 from pumped well (m), r₂ > r₁. TYPE: float
H	Aquifer thickness (m). TYPE: float

RETURNS	DESCRIPTION
dict	'k' : float — Hydraulic conductivity (m/s). 'T' : float — Transmissivity T = kH (m²/s).

Notes

• h₂ > h₁ since the head increases with distance from the well.
• r₂ > r₁ for the outer observation well.

References

Thiem (1906); Das (2021), Eq. 5.16.

Examples:

Python Console Session

>>> from geoeq.site.field_perm import pumping_test_confined
>>> res = pumping_test_confined(Q=0.004, h1=18.0, h2=19.5,
...                             r1=10, r2=50, H=20)
>>> f"{res['k']:.2e}"
'3.42e-05'

Source code in geoeq/site/field_perm.py

def pumping_test_confined(
    Q: float,
    h1: float,
    h2: float,
    r1: float,
    r2: float,
    H: float,
) -> Dict[str, float]:
    r"""
    Hydraulic conductivity from pumping test — confined aquifer.

    Thiem (1906) equilibrium equation for steady-state flow to a well
    in a confined aquifer:

    .. math::

        k = \frac{Q \ln(r_2 / r_1)}{2\,\pi\,H\,(h_2 - h_1)}

    Parameters
    ----------
    Q : float
        Steady-state pumping rate (m³/s).
    h1 : float
        Piezometric head at observation well 1 (m).
    h2 : float
        Piezometric head at observation well 2 (m), h₂ > h₁.
    r1 : float
        Radial distance of observation well 1 from pumped well (m).
    r2 : float
        Radial distance of observation well 2 from pumped well (m),
        r₂ > r₁.
    H : float
        Aquifer thickness (m).

    Returns
    -------
    dict
        ``'k'`` : float — Hydraulic conductivity (m/s).
        ``'T'`` : float — Transmissivity T = kH (m²/s).

    Notes
    -----
    * h₂ > h₁ since the head increases with distance from the well.
    * r₂ > r₁ for the outer observation well.

    References
    ----------
    Thiem (1906); Das (2021), Eq. 5.16.

    Examples
    --------
    >>> from geoeq.site.field_perm import pumping_test_confined
    >>> res = pumping_test_confined(Q=0.004, h1=18.0, h2=19.5,
    ...                             r1=10, r2=50, H=20)
    >>> f"{res['k']:.2e}"
    '3.42e-05'
    """
    check_positive(Q, "Q")
    check_positive(r1, "r1")
    check_positive(r2, "r2")
    check_positive(H, "H")
    if r2 <= r1:
        raise ValueError("r2 must be greater than r1.")
    if h2 <= h1:
        raise ValueError("h2 must be greater than h1 (head increases with distance).")

    k = Q * np.log(r2 / r1) / (2.0 * np.pi * H * (h2 - h1))
    T = k * H

    return {"k": float(k), "T": float(T)}

pumping_test_unconfined

Python

pumping_test_unconfined(Q: float, h1: float, h2: float, r1: float, r2: float) -> Dict[str, float]

Hydraulic conductivity from pumping test — unconfined aquifer.

Thiem equation for steady-state flow in an unconfined (phreatic/water-table) aquifer:

.. math::

Text Only

k = \frac{Q \ln(r_2 / r_1)}{\pi\,(h_2^2 - h_1^2)}

PARAMETER	DESCRIPTION
Q	Steady-state pumping rate (m³/s). TYPE: float
h1	Water table height at observation well 1 (m from base). TYPE: float
h2	Water table height at observation well 2 (m from base), h₂ > h₁. TYPE: float
r1	Radial distance to observation well 1 (m). TYPE: float
r2	Radial distance to observation well 2 (m), r₂ > r₁. TYPE: float

RETURNS	DESCRIPTION
dict	'k' : float — Hydraulic conductivity (m/s).

References

Thiem (1906); Das (2021), Eq. 5.17.

Examples:

Python Console Session

>>> from geoeq.site.field_perm import pumping_test_unconfined
>>> res = pumping_test_unconfined(Q=0.003, h1=15.0, h2=18.0,
...                               r1=10, r2=50)
>>> f"{res['k']:.2e}"
'1.55e-05'

Source code in geoeq/site/field_perm.py

def pumping_test_unconfined(
    Q: float,
    h1: float,
    h2: float,
    r1: float,
    r2: float,
) -> Dict[str, float]:
    r"""
    Hydraulic conductivity from pumping test — unconfined aquifer.

    Thiem equation for steady-state flow in an unconfined
    (phreatic/water-table) aquifer:

    .. math::

        k = \frac{Q \ln(r_2 / r_1)}{\pi\,(h_2^2 - h_1^2)}

    Parameters
    ----------
    Q : float
        Steady-state pumping rate (m³/s).
    h1 : float
        Water table height at observation well 1 (m from base).
    h2 : float
        Water table height at observation well 2 (m from base),
        h₂ > h₁.
    r1 : float
        Radial distance to observation well 1 (m).
    r2 : float
        Radial distance to observation well 2 (m), r₂ > r₁.

    Returns
    -------
    dict
        ``'k'`` : float — Hydraulic conductivity (m/s).

    References
    ----------
    Thiem (1906); Das (2021), Eq. 5.17.

    Examples
    --------
    >>> from geoeq.site.field_perm import pumping_test_unconfined
    >>> res = pumping_test_unconfined(Q=0.003, h1=15.0, h2=18.0,
    ...                               r1=10, r2=50)
    >>> f"{res['k']:.2e}"
    '1.55e-05'
    """
    check_positive(Q, "Q")
    check_positive(r1, "r1")
    check_positive(r2, "r2")
    check_positive(h1, "h1")
    check_positive(h2, "h2")
    if r2 <= r1:
        raise ValueError("r2 must be greater than r1.")
    if h2 <= h1:
        raise ValueError("h2 must be greater than h1.")

    k = Q * np.log(r2 / r1) / (np.pi * (h2**2 - h1**2))
    return {"k": float(k)}

lefranc_test

Python

lefranc_test(Q: float, H: float, D: float, method: str = 'constant_head') -> float

Lefranc permeability test — in-situ hydraulic conductivity.

For a constant-head test in a borehole with an open-ended cylindrical cavity of diameter D and length H:

.. math::

Text Only

k = \frac{Q}{F \, H_w}

where :math:H_w is the constant head difference and F is a shape factor. For a cylindrical cavity with L/D > 4:

.. math::

Text Only

F = \frac{2\,\pi\,H}{\ln(2\,H / D)}

PARAMETER	DESCRIPTION
Q	Flow rate (m³/s). TYPE: float
H	Head difference (m) — for constant-head test, this is the applied head; flow rate must correspond. TYPE: float
D	Borehole diameter (m). TYPE: float
method	Currently 'constant_head' is implemented. TYPE: str DEFAULT: ``'constant_head'``

RETURNS	DESCRIPTION
float	Hydraulic conductivity k (m/s).

Notes

• The Lefranc test is common in European practice (NF P 94-132).
• The shape factor assumes L/D > 4 for the simplified formula.

References

Lefranc (1936); Cassan, M. (2005). Les essais d'eau in situ.

Examples:

Python Console Session

>>> from geoeq.site.field_perm import lefranc_test
>>> k = lefranc_test(Q=1e-5, H=2.0, D=0.076)
>>> f'{k:.2e}'
'2.78e-06'

Source code in geoeq/site/field_perm.py

def lefranc_test(
    Q: float,
    H: float,
    D: float,
    method: str = "constant_head",
) -> float:
    r"""
    Lefranc permeability test — in-situ hydraulic conductivity.

    For a **constant-head** test in a borehole with an open-ended
    cylindrical cavity of diameter D and length H:

    .. math::

        k = \frac{Q}{F \, H_w}

    where :math:`H_w` is the constant head difference and F is a
    shape factor.  For a cylindrical cavity with L/D > 4:

    .. math::

        F = \frac{2\,\pi\,H}{\ln(2\,H / D)}

    Parameters
    ----------
    Q : float
        Flow rate (m³/s).
    H : float
        Head difference (m) — for constant-head test, this is the
        applied head; flow rate must correspond.
    D : float
        Borehole diameter (m).
    method : str, default ``'constant_head'``
        Currently ``'constant_head'`` is implemented.

    Returns
    -------
    float
        Hydraulic conductivity k (m/s).

    Notes
    -----
    * The Lefranc test is common in European practice (NF P 94-132).
    * The shape factor assumes L/D > 4 for the simplified formula.

    References
    ----------
    Lefranc (1936); Cassan, M. (2005). *Les essais d'eau in situ*.

    Examples
    --------
    >>> from geoeq.site.field_perm import lefranc_test
    >>> k = lefranc_test(Q=1e-5, H=2.0, D=0.076)
    >>> f'{k:.2e}'
    '2.78e-06'
    """
    check_positive(Q, "Q")
    check_positive(H, "H")
    check_positive(D, "D")

    # Shape factor for cylindrical cavity, L=H, diameter D
    F = 2.0 * np.pi * H / np.log(2.0 * H / D)
    k = Q / (F * H)
    return float(k)

Field CBR / DCP

dcp_cbr

Python

dcp_cbr(dcpi: Union[float, ndarray], method: str = 'webster') -> Union[float, np.ndarray]

Estimate CBR from Dynamic Cone Penetrometer Index (DCPI).

The DCP measures penetration per blow (mm/blow). Several empirical correlations relate DCPI to CBR:

Webster et al. (1992) — US Army Corps of Engineers:

.. math::

Text Only

\log_{10}(\text{CBR}) = 1.12 - 0.39 \,\log_{10}(\text{DCPI})

i.e. :math:\text{CBR} = 10^{1.12 - 0.39 \log(\text{DCPI})}

TRL (1993) — Transport Research Laboratory:

.. math::

Text Only

\log_{10}(\text{CBR}) = 2.48 - 1.057 \,\log_{10}(\text{DCPI})

Kleyn (1975) — South African method:

.. math::

Text Only

\log_{10}(\text{CBR}) = 2.62 - 1.27 \,\log_{10}(\text{DCPI})

PARAMETER	DESCRIPTION
dcpi	DCP Index — penetration per blow (mm/blow). TYPE: float or array_like
method	Correlation method: 'webster', 'trl', or 'kleyn'. TYPE: str DEFAULT: ``'webster'``

RETURNS	DESCRIPTION
float or ndarray	Estimated CBR (%).

References

Webster et al. (1992); TRL (1993); Kleyn (1975).

Examples:

Python Console Session

>>> from geoeq.site.field_cbr import dcp_cbr
>>> round(dcp_cbr(dcpi=10, method='webster'), 1)
5.4
>>> round(dcp_cbr(dcpi=10, method='trl'), 1)
26.5

Source code in geoeq/site/field_cbr.py

def dcp_cbr(
    dcpi: Union[float, np.ndarray],
    method: str = "webster",
) -> Union[float, np.ndarray]:
    r"""
    Estimate CBR from Dynamic Cone Penetrometer Index (DCPI).

    The DCP measures penetration per blow (mm/blow).  Several
    empirical correlations relate DCPI to CBR:

    **Webster et al. (1992)** — US Army Corps of Engineers:

    .. math::

        \log_{10}(\text{CBR}) = 1.12 - 0.39 \,\log_{10}(\text{DCPI})

    i.e.  :math:`\text{CBR} = 10^{1.12 - 0.39 \log(\text{DCPI})}`

    **TRL (1993)** — Transport Research Laboratory:

    .. math::

        \log_{10}(\text{CBR}) = 2.48 - 1.057 \,\log_{10}(\text{DCPI})

    **Kleyn (1975)** — South African method:

    .. math::

        \log_{10}(\text{CBR}) = 2.62 - 1.27 \,\log_{10}(\text{DCPI})

    Parameters
    ----------
    dcpi : float or array_like
        DCP Index — penetration per blow (mm/blow).
    method : str, default ``'webster'``
        Correlation method: ``'webster'``, ``'trl'``, or ``'kleyn'``.

    Returns
    -------
    float or ndarray
        Estimated CBR (%).

    References
    ----------
    Webster et al. (1992); TRL (1993); Kleyn (1975).

    Examples
    --------
    >>> from geoeq.site.field_cbr import dcp_cbr
    >>> round(dcp_cbr(dcpi=10, method='webster'), 1)
    5.4
    >>> round(dcp_cbr(dcpi=10, method='trl'), 1)
    26.5
    """
    d = np.asarray(dcpi, dtype=float)
    check_positive(d, "dcpi")
    method_l = method.lower()

    if method_l == "webster":
        cbr = 10 ** (1.12 - 0.39 * np.log10(d))
    elif method_l == "trl":
        cbr = 10 ** (2.48 - 1.057 * np.log10(d))
    elif method_l == "kleyn":
        cbr = 10 ** (2.62 - 1.27 * np.log10(d))
    else:
        raise ValueError(
            f"Unknown method '{method}'. Choose: webster, trl, kleyn."
        )

    if all(np.ndim(x) == 0 for x in [dcpi]):
        return float(cbr)
    return np.asarray(cbr)

field_cbr_test

Python

field_cbr_test(penetration: Union[list, ndarray], load: Union[list, ndarray], area: float = 1935.5) -> Dict[str, float]

Field (in-situ) CBR from load–penetration data.

The procedure is identical to the laboratory CBR (ASTM D1883) but performed on the in-situ subgrade.

.. math::

Text Only

\text{CBR} = \frac{\text{Test load at standard penetration}}
                  {\text{Standard load}} \times 100

Standard reference loads (for 1935.5 mm² = 3 in² piston):

• At 2.54 mm (0.1 in.): 13.24 kN
• At 5.08 mm (0.2 in.): 19.96 kN

PARAMETER	DESCRIPTION
penetration	Penetration values (mm). TYPE: array_like
load	Corresponding loads (kN). TYPE: array_like
area	Piston area (mm²). Standard 3 in² = 1935.5 mm². TYPE: float DEFAULT: 1935.5

RETURNS	DESCRIPTION
dict	'CBR' : float — Governing CBR (%). 'CBR_2_54' : float — CBR at 2.54 mm. 'CBR_5_08' : float — CBR at 5.08 mm. 'load_2_54' : float — Load at 2.54 mm (kN). 'load_5_08' : float — Load at 5.08 mm (kN).

References

ASTM D4429 (2009). Standard Test Method for CBR (California Bearing Ratio) of Soils in Place.

Examples:

Python Console Session

>>> from geoeq.site.field_cbr import field_cbr_test
>>> pen = [0, 0.64, 1.27, 1.91, 2.54, 3.81, 5.08, 7.62, 10.16, 12.70]
>>> load = [0, 0.5, 1.2, 2.5, 4.0, 6.5, 9.0, 13.0, 16.0, 18.0]
>>> res = field_cbr_test(pen, load)
>>> res['CBR'] > 0
True

Source code in geoeq/site/field_cbr.py

def field_cbr_test(
    penetration: Union[list, np.ndarray],
    load: Union[list, np.ndarray],
    area: float = 1935.5,
) -> Dict[str, float]:
    r"""
    Field (in-situ) CBR from load–penetration data.

    The procedure is identical to the laboratory CBR (ASTM D1883)
    but performed on the in-situ subgrade.

    .. math::

        \text{CBR} = \frac{\text{Test load at standard penetration}}
                          {\text{Standard load}} \times 100

    Standard reference loads (for 1935.5 mm² = 3 in² piston):

    - At 2.54 mm (0.1 in.): 13.24 kN
    - At 5.08 mm (0.2 in.): 19.96 kN

    Parameters
    ----------
    penetration : array_like
        Penetration values (mm).
    load : array_like
        Corresponding loads (kN).
    area : float, default 1935.5
        Piston area (mm²).  Standard 3 in² = 1935.5 mm².

    Returns
    -------
    dict
        ``'CBR'`` : float — Governing CBR (%).
        ``'CBR_2_54'`` : float — CBR at 2.54 mm.
        ``'CBR_5_08'`` : float — CBR at 5.08 mm.
        ``'load_2_54'`` : float — Load at 2.54 mm (kN).
        ``'load_5_08'`` : float — Load at 5.08 mm (kN).

    References
    ----------
    ASTM D4429 (2009). Standard Test Method for CBR (California
    Bearing Ratio) of Soils in Place.

    Examples
    --------
    >>> from geoeq.site.field_cbr import field_cbr_test
    >>> pen = [0, 0.64, 1.27, 1.91, 2.54, 3.81, 5.08, 7.62, 10.16, 12.70]
    >>> load = [0, 0.5, 1.2, 2.5, 4.0, 6.5, 9.0, 13.0, 16.0, 18.0]
    >>> res = field_cbr_test(pen, load)
    >>> res['CBR'] > 0
    True
    """
    pen = np.asarray(penetration, dtype=float)
    ld = np.asarray(load, dtype=float)
    if len(pen) != len(ld):
        raise ValueError("penetration and load must have the same length.")
    check_positive(area, "area")

    # Standard reference loads (kN) for 1935.5 mm² piston
    std_load_2_54 = 13.24
    std_load_5_08 = 19.96

    # Interpolate loads at standard penetrations
    load_2_54 = float(np.interp(2.54, pen, ld))
    load_5_08 = float(np.interp(5.08, pen, ld))

    CBR_2_54 = (load_2_54 / std_load_2_54) * 100.0
    CBR_5_08 = (load_5_08 / std_load_5_08) * 100.0

    CBR = max(CBR_2_54, CBR_5_08)

    return {
        "CBR": float(CBR),
        "CBR_2_54": float(CBR_2_54),
        "CBR_5_08": float(CBR_5_08),
        "load_2_54": load_2_54,
        "load_5_08": load_5_08,
    }