Skip to content

Seepage & flow

Five functions covering Darcy's law, the conditions for piping, and flow-net based seepage estimation.

Quick example

Python
import geoeq as ge

# Darcy's law: Q = k * i * A
ge.darcy_flow(k=1e-4, i=0.1, A=1.0)              # 1e-5 m^3/s

# Hydraulic gradient
ge.hydraulic_gradient(dh=2.0, L=10.0)            # 0.20

# Critical (boiling/heave) gradient
ge.critical_gradient(Gs=2.65, e=0.65)            # 1.0

# Equivalent permeability of a layered system
ge.equivalent_k([1e-3, 1e-5], [5, 5], direction="horizontal")
ge.equivalent_k([1e-3, 1e-5], [5, 5], direction="vertical")

# Seepage from a flow net
ge.flow_net(Nf=4, Nd=8, k=1e-5, dh=10)           # 5e-5 m^3/s per m

Critical gradient — when quicksand happens

\[ i_{\text{cr}} = \frac{G_s - 1}{1 + e} \]

When the upward seepage gradient reaches \(i_{\text{cr}}\), effective stress vanishes and the soil "boils". For a typical sand (\(G_s = 2.65\), \(e = 0.65\)), \(i_{\text{cr}} = 1.0\) — the classic textbook answer.

Anisotropy in layered ground

Horizontal and vertical permeabilities can differ by orders of magnitude in a stratified profile. With \(k_1 = 10^{-3}\) and \(k_2 = 10^{-5}\) in equal-thickness layers:

  • \(k_{\text{eq}}^{H} \approx 5 \times 10^{-4}\) — the high-k layer dominates
  • \(k_{\text{eq}}^{V} \approx 2 \times 10^{-5}\) — the low-k layer dominates

This 25× ratio is the engineering argument for filter drains (horizontal flow allowed) and impermeable cores (vertical flow blocked).

API reference

darcy_flow

Python
darcy_flow(k: float, i: float, A: float) -> float

Darcy's law: Q = k * i * A (m^3/s).

PARAMETER DESCRIPTION
k

Hydraulic conductivity (m/s).

TYPE: float

i

Hydraulic gradient (-, dimensionless).

TYPE: float

A

Cross-sectional area (m^2).

TYPE: float

Reference

Darcy (1856); Das (2010) Eq. 5.11.

Source code in geoeq/design/seepage.py
Python
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def darcy_flow(k: float, i: float, A: float) -> float:
    """Darcy's law: Q = k * i * A (m^3/s).

    Parameters
    ----------
    k : float
        Hydraulic conductivity (m/s).
    i : float
        Hydraulic gradient (-, dimensionless).
    A : float
        Cross-sectional area (m^2).

    Reference
    ---------
    Darcy (1856); Das (2010) Eq. 5.11.
    """
    check_positive(k, "k")
    check_positive(A, "A")
    return k * i * A

hydraulic_gradient

Python
hydraulic_gradient(dh: float, L: float) -> float

Hydraulic gradient i = delta h / L.

PARAMETER DESCRIPTION
dh

Head loss (m).

TYPE: float

L

Flow-path length (m).

TYPE: float

Reference

Das (2010) Eq. 5.10.

Source code in geoeq/design/seepage.py
Python
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def hydraulic_gradient(dh: float, L: float) -> float:
    """Hydraulic gradient i = delta h / L.

    Parameters
    ----------
    dh : float
        Head loss (m).
    L : float
        Flow-path length (m).

    Reference
    ---------
    Das (2010) Eq. 5.10.
    """
    check_positive(L, "L")
    return dh / L

critical_gradient

Python
critical_gradient(Gs: float, e: float) -> float

Critical (boiling) hydraulic gradient i_cr = (Gs - 1)/(1 + e).

When i >= i_cr in upward flow, effective stress vanishes and quicksand (heave) develops.

Reference

Das (2010) Eq. 5.18; Terzaghi & Peck (1948).

Source code in geoeq/design/seepage.py
Python
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def critical_gradient(Gs: float, e: float) -> float:
    """Critical (boiling) hydraulic gradient i_cr = (Gs - 1)/(1 + e).

    When i >= i_cr in upward flow, effective stress vanishes and quicksand
    (heave) develops.

    Reference
    ---------
    Das (2010) Eq. 5.18; Terzaghi & Peck (1948).
    """
    check_positive(Gs, "Gs")
    check_positive(e, "e")
    if Gs <= 1:
        raise ValueError("Gs must be > 1.")
    return (Gs - 1) / (1 + e)

equivalent_k

Python
equivalent_k(k_layers: Sequence[float], H_layers: Sequence[float], direction: str = 'horizontal') -> float

Equivalent permeability of a layered system.

horizontal: k_eq = sum(k_i * H_i) / sum(H_i) (parallel flow) vertical: k_eq = sum(H_i) / sum(H_i / k_i) (series flow)

Reference

Das (2010) Eq. 5.14 (horizontal), 5.15 (vertical).

Source code in geoeq/design/seepage.py
Python
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def equivalent_k(
    k_layers: Sequence[float],
    H_layers: Sequence[float],
    direction: str = "horizontal",
) -> float:
    """Equivalent permeability of a layered system.

    horizontal: k_eq = sum(k_i * H_i) / sum(H_i)        (parallel flow)
    vertical:   k_eq = sum(H_i) / sum(H_i / k_i)        (series flow)

    Reference
    ---------
    Das (2010) Eq. 5.14 (horizontal), 5.15 (vertical).
    """
    k = np.asarray(k_layers, dtype=float)
    H = np.asarray(H_layers, dtype=float)
    if k.shape != H.shape:
        raise ValueError("k_layers and H_layers must have the same length.")
    check_positive(k, "k_layers")
    check_positive(H, "H_layers")
    direction = direction.lower()
    if direction in ("horizontal", "h", "parallel"):
        return float(np.sum(k * H) / np.sum(H))
    elif direction in ("vertical", "v", "series"):
        return float(np.sum(H) / np.sum(H / k))
    raise ValueError("direction must be 'horizontal' or 'vertical'.")

flow_net

Python
flow_net(Nf: float, Nd: float, k: float, dh: float, L: float = 1.0) -> float

Seepage discharge from a flow net.

Q = k * dh * (Nf / Nd) * L (m^3/s per metre length normal to plane)

PARAMETER DESCRIPTION
Nf

Number of flow channels.

TYPE: float

Nd

Number of equipotential drops.

TYPE: float

k

Hydraulic conductivity (m/s).

TYPE: float

dh

Total head loss across the system (m).

TYPE: float

L

Length perpendicular to the flow plane (m). Default 1 m.

TYPE: float DEFAULT: 1.0

Reference

Das (2010) Eq. 5.20.

Source code in geoeq/design/seepage.py
Python
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def flow_net(Nf: float, Nd: float, k: float, dh: float, L: float = 1.0) -> float:
    """Seepage discharge from a flow net.

    Q = k * dh * (Nf / Nd) * L   (m^3/s per metre length normal to plane)

    Parameters
    ----------
    Nf : float
        Number of flow channels.
    Nd : float
        Number of equipotential drops.
    k : float
        Hydraulic conductivity (m/s).
    dh : float
        Total head loss across the system (m).
    L : float
        Length perpendicular to the flow plane (m). Default 1 m.

    Reference
    ---------
    Das (2010) Eq. 5.20.
    """
    check_positive(Nf, "Nf")
    check_positive(Nd, "Nd")
    check_positive(k, "k")
    return k * dh * (Nf / Nd) * L