Soil Models

Martin Vonk (2025)

Overview

Soil hydraulic models describe the relationship between:

• Pressure head (h) - the soil water suction or tension

• Water content (θ) - the volumetric water content

• Hydraulic conductivity (K) - the ease of water movement through soil

These relationships, known as the soil water retention curve (SWRC) and hydraulic conductivity function (HCF), are essential for modeling variably saturated groundwater flow.

pedon provides different soil model formulations from the scientific literature. In this notebook, we’ll demonstrate several of them and show how to compare their behavior.

import inspect
import sys

import matplotlib.pyplot as plt
import numpy as np

import pedon as pe

Available Soil Models

pedon includes the following soil hydraulic models from the literature:

# classes for soil models
soilmodels = [
    cls_name
    for cls_name, cls_obj in inspect.getmembers(sys.modules["pedon.soilmodel"])
    if inspect.isclass(cls_obj)
]
soilmodels

['Brooks',
 'Fredlund',
 'Gardner',
 'Genuchten',
 'GenuchtenGardner',
 'Haverkamp',
 'Panday',
 'Protocol',
 'Rucker',
 'SoilModel']

Setting up Common Parameters

We’ll use typical soil properties throughout this notebook. Note: pedon does not assume units, but it’s good practice to be consistent. Here we use cm for length and cm/d for conductivity.

# shared properties
k_s = 100  # saturated conductivity [cm/d]
theta_r = 0.03  # residual water content [-]
theta_s = 0.42  # saturated water content [-]

# Mualem-van Genuchten
alpha = 0.04  # shape parameter [1/cm]
n = 1.4  # shape parameter [-]

gen = pe.Genuchten(
    k_s=k_s,
    theta_r=theta_r,
    theta_s=theta_s,
    alpha=alpha,
    n=n,
)

# Brooks-Corey
h_b = 10  # bubbling pressure [cm]
# connectivity parameter [-]
l = 1.1  # noqa: E741

bro = pe.Brooks(
    k_s=k_s,
    theta_r=theta_r,
    theta_s=theta_s,
    h_b=h_b,
    l=l,
)

# haverkamp
alpha = 1.611e6
beta = 3.96
a = 1.175e6
hav = pe.Haverkamp(
    k_s=k_s,
    theta_r=theta_r,
    theta_s=theta_s,
    alpha=alpha,
    beta=beta,
    a=a,
)

# Gardner-Kozeny
c = 0.005  # shape parameter k curve
m = 0.01  # shape parameter theta curve
gar = pe.Gardner(k_s=k_s, theta_s=theta_s, c=c, m=m)

# Fredlund-Xing
a = 10.0  # shape parameter
n = 5.0  # shape parameter
m = 0.5  # shape parameter
fre = pe.Fredlund(
    k_s=k_s,
    theta_s=theta_s,
    a=a,
    n=n,
    m=m,
)

# Quick plot method
ax = gen.plot()
bro.plot(ax=ax)
gar.plot(ax=ax)
hav.plot(ax=ax)
fre.plot(ax=ax)
ax.set_xlim(0.0, 0.5)
ax.legend();

# More extensive plot method
f, axs = plt.subplots(1, 2, figsize=(7.0, 6.0), sharey=True, layout="tight")

pe.plot_swrc(gen, ax=axs[0])
pe.plot_swrc(bro, ax=axs[0])
pe.plot_swrc(gar, ax=axs[0])
pe.plot_swrc(hav, ax=axs[0])
pe.plot_swrc(fre, ax=axs[0])
axs[0].set_title("Soil Water Retention Curve")
axs[0].set_xlabel("\N{GREEK SMALL LETTER THETA} [-]")
axs[0].set_xlim(-0.01, 0.5)

pe.plot_hcf(gen, ax=axs[1])
pe.plot_hcf(bro, ax=axs[1])
pe.plot_hcf(gar, ax=axs[1])
pe.plot_hcf(hav, ax=axs[1])
pe.plot_hcf(fre, ax=axs[1])
axs[1].set_title("Hydraulic Conductivity Function")
axs[1].set_xlabel("Ks [cm/d]")
axs[1].set_xscale("log")
axs[1].set_xlim(1e-10, 1e3)

axs[0].set_yscale("log")
axs[0].set_ylabel("|\N{GREEK SMALL LETTER PSI}| [cm]")

Text(0, 0.5, '|ψ| [cm]')

Other soil models that are available are Panday (combination of Genuchten and Brooks-Corey).

Creating Custom Soil Models

The built-in models cover most common use cases, but pedon also supports creating your own custom soil models.

Custom models are useful when:

• You need to implement a proprietary or research-specific formulation

• You want to combine features from different models

• You’re testing new hydraulic function formulations

The key requirement is implementing the SoilModel protocol, which defines six methods: theta(), s(), k_r(), k(), h(), and plot().

pe.SoilModel?

Init signature: pe.SoilModel(*args, **kwargs)
Docstring:
Protocol for soil models

This protocol defines the interface for custom soil models. To create a custom
soil model, implement all methods defined below. Your class will automatically
conform to this protocol.

Example:
    >>> @dataclass
    ... class CustomModel:
    ...     k_s: float
    ...     theta_s: float
    ...
    ...     def theta(self, h: FloatArray) -> FloatArray:
    ...         # Calculate soil moisture content from pressure head
    ...         return ...
    ...
    ...     def s(self, h: FloatArray) -> FloatArray:
    ...         # Calculate effective saturation from pressure head
    ...         return (self.theta(h) - theta_r) / (theta_s - theta_r)
    ...
    ...     def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
    ...         # Calculate relative permeability from pressure head or saturation
    ...         return ...
    ...
    ...     def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
    ...         # Calculate hydraulic conductivity from pressure head or saturation
    ...         return self.k_s * self.k_r(h=h, s=s)
    ...
    ...     def h(self, theta: FloatArray) -> FloatArray:
    ...         # Inverse of theta method
    ...         return ...
    ...
    ...     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
    ...         # Plot the soil water retention curve by calling `plot_swrc`
    ...         return plot_swrc(self, ax=ax)
File:           ~/repos/pedon/src/pedon/soilmodel.py
Type:           _ProtocolMeta
Subclasses:

from dataclasses import dataclass


@dataclass
class CustomSoilModel:
    """Power-law soil model with custom permeability

    A more flexible model that uses power-law relationships for both
    water retention and relative permeability.
    """

    k_s: float  # Saturated hydraulic conductivity [cm/d]
    theta_r: float  # Residual water content [-]
    theta_s: float  # Saturated water content [-]
    h_ref: float  # Reference pressure head [cm]
    power_retention: float  # Power exponent for retention curve
    power_conductivity: float  # Power exponent for conductivity

    def theta(self, h):
        """Water content using power law"""
        # Normalized pressure head
        h_norm = np.abs(h) / self.h_ref

        s = 1.0 / (1.0 + h_norm**self.power_retention)
        return self.theta_r + (self.theta_s - self.theta_r) * s

    def s(self, h):
        """Effective saturation"""
        return (self.theta(h) - self.theta_r) / (self.theta_s - self.theta_r)

    def k_r(self, h, s=None):
        """Relative permeability"""
        if s is None:
            s = self.s(h)
        # Custom power law for conductivity
        return s**self.power_conductivity

    def k(self, h, s=None):
        """Hydraulic conductivity"""
        return self.k_s * self.k_r(h=h, s=s)

    def h(self, theta):
        """Inverse function"""
        s = (theta - self.theta_r) / (self.theta_s - self.theta_r)
        # From: s = 1/(1 + h_norm^power) => h_norm = (1/s - 1)^(1/power)
        h_norm = (1.0 / s - 1.0) ** (1.0 / self.power_retention)
        return h_norm * self.h_ref

    def plot(self, ax=None):
        """Plot the soil water retention curve"""
        return pe.plot_swrc(self, ax=ax)

# Create instances with different exponents
soil_one = CustomSoilModel(
    k_s=10,
    theta_r=0.02,
    theta_s=0.4,
    h_ref=10.0,
    power_retention=2.5,
    power_conductivity=4.0,
)

soil_two = CustomSoilModel(
    k_s=200,
    theta_r=0.08,
    theta_s=0.45,
    h_ref=5.0,
    power_retention=1.0,
    power_conductivity=2.0,
)

Due to the duck typing the CustomSoilModel is checkable as an instance of the SoilModel protocol. This ensures that the model can also be used for other methods in pedon.

print(f"CustomSoilModel (clay) is a SoilModel: {isinstance(soil_one, pe.SoilModel)}")
print(f"CustomSoilModel (sand) is a SoilModel: {isinstance(soil_two, pe.SoilModel)}")

CustomSoilModel (clay) is a SoilModel: True
CustomSoilModel (sand) is a SoilModel: True

Quick plot to check if it works.

ax = soil_one.plot()
soil_two.plot(ax=ax)

<Axes: >