Monin–Obukhov dispersion¶

Pasquill–Gifford dispersion driven by Monin–Obukhov scales: Golder's (1972) class boundaries map the inverse Obukhov length 1/L and the roughness length z₀ to a PG class, whose Briggs σ_y and σ_z curves then give the plume spread. GolderDispersion is built from a column of plumax.met.SurfaceLayer and passed as dispersion_params to plume_concentration or as stability_class to simulate_plume. The hard class lookup follows Golder; the optional interpolate=True blends neighbouring classes linearly in 1/L so σ is differentiable in L, which is a plumax modelling choice. PG classes given by letter remain the default. Tier I.

mo_dispersion ¶

Monin–Obukhov-driven PG dispersion via Golder's (1972) classes (implementation in plumax.gauss_plume._src.mo_dispersion).

GolderDispersion ¶

GolderDispersion(inverse_obukhov_length: Float[Array, ''] | float, z0: float, interpolate: bool = False)

Bases: Module

PG σ_y, σ_z with the class chosen from Monin–Obukhov scales (Golder 1972).

σ_i(x) = σ_i^PG(x; class(1/L, z₀)), i ∈ {y, z}, with the Briggs–PG curves of :func:~plumax.gauss_plume.calculate_briggs_dispersion. With interpolate=True the class index s = :func:golder_stability_index is continuous and σ_i(x) = (1 − w) σ_i^PG(x; ⌊s⌋) + w σ_i^PG(x; ⌊s⌋ + 1), w = s − ⌊s⌋, which is differentiable in L (a plumax modelling choice, not part of Golder). Pass it as dispersion_params to :func:~plumax.gauss_plume.plume_concentration or as stability_class to :func:~plumax.gauss_plume.simulate_plume.

Attributes:

Name	Type	Description
`inverse_obukhov_length`	`Float[Array, '']`	1/L [m⁻¹] of the column; zero is neutral.
`z0`	`float`	Roughness length z₀ > 0 [m]; clamped to :data:`GOLDER_Z0_RANGE`.
`interpolate`	`bool`	Blend neighbouring PG classes linearly in 1/L (default `False`: the hard class lookup, as in Golder and MMIF).

Raises:

Type	Description
`ParameterError`	If `z0` is not positive.

Examples:

Neutral air reproduces PG class D exactly:

>>> import jax.numpy as jnp
>>> from plumax.gauss_plume import calculate_briggs_dispersion
>>> from plumax.gauss_plume import get_dispersion_params
>>> distance = jnp.array([100.0, 1000.0])
>>> neutral = GolderDispersion(inverse_obukhov_length=0.0, z0=0.1)
>>> sigma_y, sigma_z = neutral(distance)
>>> pg_y, pg_z = calculate_briggs_dispersion(distance, get_dispersion_params("D"))
>>> bool(jnp.allclose(sigma_y, pg_y) & jnp.allclose(sigma_z, pg_z))
True

The interpolated σ moves smoothly with L inside a class:

>>> import jax
>>> def sigma_z_at_1km(inverse_length):
...     dispersion = GolderDispersion(inverse_length, z0=0.1, interpolate=True)
...     return dispersion(jnp.array(1000.0))[1]
>>> bool(jax.grad(sigma_z_at_1km)(-1.0 / 50.0) < 0.0)
True

Build it from a met column's surface-layer diagnosis:

>>> from plumax.les_fvm.grid import make_grid
>>> from plumax.met import MetField, MoninObukhov
>>> grid = make_grid(
...     (0.0, 400.0, 4), (0.0, 400.0, 4), (0.0, 200.0, 4), dtype=jnp.float64
... )
>>> shape = (1, 4, 4, 4)
>>> field = MetField(
...     plume_grid=grid, times=jnp.zeros(1), u=jnp.full(shape, 5.0),
...     v=jnp.zeros(shape), w=jnp.zeros(shape),
...     temperature=jnp.full(shape, 300.0), pressure=jnp.full(shape, 1.0e5),
...     pbl_height=jnp.full((1, 4, 4), 800.0),
... )
>>> layer = MoninObukhov()(field, z0=0.1, sensible_heat_flux=0.0)
>>> GolderDispersion.from_surface_layer(layer, z0=0.1, it=0, iy=1, ix=2).z0
0.1

References

Golder, D. (1972). Relations among stability parameters in the surface layer. Boundary-Layer Meteorology 3, 47–58. US EPA (2024). MMIF v4.1.1, pasquill_gifford.f90 (pg_golder).

Source code in src/plumax/gauss_plume/_src/mo_dispersion.py

def __init__(
    self,
    inverse_obukhov_length: Float[Array, ""] | float,
    z0: float,
    interpolate: bool = False,
) -> None:
    self.inverse_obukhov_length = jnp.asarray(inverse_obukhov_length)
    self.z0 = float(z0)
    self.interpolate = bool(interpolate)

from_surface_layer `classmethod` ¶

from_surface_layer(layer: SurfaceLayer, *, z0: float, it: int, iy: int, ix: int, interpolate: bool = False) -> GolderDispersion

Golder dispersion at one column of a :class:plumax.met.SurfaceLayer.

Parameters:

Name	Type	Description	Default
`layer`	`SurfaceLayer`	Surface-layer scales from :class:`plumax.met.MoninObukhov`; only its 1/L is used.	required
`z0`	`float`	The roughness length z₀ [m] the layer was diagnosed with.	required
`it`	`int`	Met time index.	required
`iy`	`int`	Analysis y index.	required
`ix`	`int`	Analysis x index.	required
`interpolate`	`bool`	See the class attribute.	`False`

Returns:

Type	Description
`GolderDispersion`	The dispersion model for that column.

Source code in src/plumax/gauss_plume/_src/mo_dispersion.py

@classmethod
def from_surface_layer(
    cls,
    layer: SurfaceLayer,
    *,
    z0: float,
    it: int,
    iy: int,
    ix: int,
    interpolate: bool = False,
) -> GolderDispersion:
    """Golder dispersion at one column of a :class:`plumax.met.SurfaceLayer`.

    Args:
        layer: Surface-layer scales from :class:`plumax.met.MoninObukhov`;
            only its 1/L is used.
        z0: The roughness length z₀ [m] the layer was diagnosed with.
        it: Met time index.
        iy: Analysis y index.
        ix: Analysis x index.
        interpolate: See the class attribute.

    Returns:
        The dispersion model for that column.
    """
    inverse_length = layer.inverse_obukhov_length[it, iy, ix]
    return cls(inverse_length, z0=z0, interpolate=interpolate)

golder_class_boundaries ¶

golder_class_boundaries(z0: float) -> Float[Array, ' boundary']

Inverse Obukhov lengths at Golder's five PG class boundaries.

(1/L)_k = (ln z₀ − d_k) / c_k for the A/B, B/C, C/D, D/E and E/F boundaries, with z₀ clamped to :data:GOLDER_Z0_RANGE.

Parameters:

Name	Type	Description	Default
`z0`	`float`	Roughness length z₀ > 0 [m].	required

Returns:

Type	Description
`Float[Array, ' boundary']`	The five boundaries [m⁻¹], increasing from A/B (most negative) to E/F.

Examples:

At z₀ = 0.1 m the unstable boundaries sit at L ≈ −10.5, −30.4 and −104 m, the stable ones at L ≈ 112 and 26.9 m:

>>> import jax.numpy as jnp
>>> boundaries = golder_class_boundaries(0.1)
>>> [round(float(1.0 / b), 1) for b in boundaries]
[-10.5, -30.4, -104.1, 111.6, 26.9]

z₀ beyond 0.5 m is clamped:

>>> bool(jnp.allclose(golder_class_boundaries(2.0), golder_class_boundaries(0.5)))
True

Source code in src/plumax/gauss_plume/_src/mo_dispersion.py

def golder_class_boundaries(z0: float) -> Float[Array, " boundary"]:
    """Inverse Obukhov lengths at Golder's five PG class boundaries.

    (1/L)_k = (ln z₀ − d_k) / c_k for the A/B, B/C, C/D, D/E and E/F
    boundaries, with z₀ clamped to :data:`GOLDER_Z0_RANGE`.

    Args:
        z0: Roughness length z₀ > 0 [m].

    Returns:
        The five boundaries [m⁻¹], increasing from A/B (most negative) to E/F.

    Examples:
        At z₀ = 0.1 m the unstable boundaries sit at L ≈ −10.5, −30.4 and
        −104 m, the stable ones at L ≈ 112 and 26.9 m:

        >>> import jax.numpy as jnp
        >>> boundaries = golder_class_boundaries(0.1)
        >>> [round(float(1.0 / b), 1) for b in boundaries]
        [-10.5, -30.4, -104.1, 111.6, 26.9]

        z₀ beyond 0.5 m is clamped:

        >>> bool(jnp.allclose(golder_class_boundaries(2.0), golder_class_boundaries(0.5)))
        True
    """
    roughness = jnp.asarray(z0)
    roughness = eqx.error_if(roughness, roughness <= 0.0, "z0 must be positive")
    clamped = jnp.clip(roughness, *GOLDER_Z0_RANGE)
    log_z0 = jnp.log(clamped)
    scale = jnp.asarray(GOLDER_BOUNDARY_SCALE, dtype=log_z0.dtype)
    offset = jnp.asarray(GOLDER_BOUNDARY_OFFSET, dtype=log_z0.dtype)
    # (1/L)_k = (ln z₀ − d_k) / c_k
    return (log_z0 - offset) / scale

golder_stability_class ¶

golder_stability_class(inverse_obukhov_length: Float[Array, '*shape'] | float, z0: float) -> Int[Array, '*shape']

Golder's (1972) PG class for 1/L and z₀, as an index 0 (A) … 5 (F).

The class is the number of boundaries (1/L)_k strictly below 1/L, so a value on a boundary goes to the less stable class and 1/L = 0 is D.

Parameters:

Name	Type	Description	Default
`inverse_obukhov_length`	`Float[Array, '*shape'] \| float`	1/L [m⁻¹]; zero is neutral, negative unstable.	required
`z0`	`float`	Roughness length z₀ > 0 [m].	required

Returns:

Type	Description
`Int[Array, '*shape']`	The class index into :data:`~plumax.gauss_plume.STABILITY_CLASSES`.

Examples:

Over z₀ = 0.1 m, L = −20 m is B, neutral is D and L = 50 m is E:

>>> import jax.numpy as jnp
>>> inverse = jnp.array([-1.0 / 20.0, 0.0, 1.0 / 50.0])
>>> golder_stability_class(inverse, 0.1).tolist()
[1, 3, 4]

Source code in src/plumax/gauss_plume/_src/mo_dispersion.py

def golder_stability_class(
    inverse_obukhov_length: Float[Array, "*shape"] | float, z0: float
) -> Int[Array, "*shape"]:
    """Golder's (1972) PG class for 1/L and z₀, as an index 0 (A) … 5 (F).

    The class is the number of boundaries (1/L)_k strictly below 1/L, so a
    value on a boundary goes to the less stable class and 1/L = 0 is D.

    Args:
        inverse_obukhov_length: 1/L [m⁻¹]; zero is neutral, negative unstable.
        z0: Roughness length z₀ > 0 [m].

    Returns:
        The class index into :data:`~plumax.gauss_plume.STABILITY_CLASSES`.

    Examples:
        Over z₀ = 0.1 m, L = −20 m is B, neutral is D and L = 50 m is E:

        >>> import jax.numpy as jnp
        >>> inverse = jnp.array([-1.0 / 20.0, 0.0, 1.0 / 50.0])
        >>> golder_stability_class(inverse, 0.1).tolist()
        [1, 3, 4]
    """
    inverse_obukhov_length = jnp.asarray(inverse_obukhov_length)
    boundaries = golder_class_boundaries(z0)
    above = inverse_obukhov_length[..., None] > boundaries
    return jnp.sum(above, axis=-1)

golder_stability_index ¶

golder_stability_index(inverse_obukhov_length: Float[Array, '*shape'] | float, z0: float) -> Float[Array, '*shape']

Continuous PG class index from 1/L and z₀, piecewise linear in 1/L.

s(1/L) interpolates linearly through the knots (1/L, s) = ((1/L)_AB − Δ_A, 0), ((1/L)_AB, 0.5), ((1/L)_BC, 1.5), ((1/L)_CD, 2.5), (0, 3), ((1/L)_DE, 3.5), ((1/L)_EF, 4.5), ((1/L)_EF + Δ_F, 5), and is clamped to [0, 5] outside them, with Δ_A = ((1/L)_BC − (1/L)_AB)/2 and Δ_F = ((1/L)_EF − (1/L)_DE)/2 (the half-width of the neighbouring class).

This is a plumax modelling choice, not part of Golder (1972): Golder gives only the class boundaries. Each boundary maps to the half-integer between its two classes, so round of the index (ties to the less stable class) is the hard class of :func:golder_stability_class, and neutral maps to exactly 3 (class D). Its derivative in 1/L is non-zero everywhere except deep inside A (1/L < (1/L)_AB − Δ_A) and F (1/L > (1/L)_EF + Δ_F).

Parameters:

Name	Type	Description	Default
`inverse_obukhov_length`	`Float[Array, '*shape'] \| float`	1/L [m⁻¹]; zero is neutral, negative unstable.	required
`z0`	`float`	Roughness length z₀ > 0 [m] (a scalar).	required

Returns:

Type	Description
`Float[Array, '*shape']`	The class index s ∈ [0, 5], 0 for A and 5 for F.

Examples:

Boundaries land on half-integers and neutral on D:

>>> import jax.numpy as jnp
>>> boundaries = golder_class_boundaries(0.1)
>>> golder_stability_index(boundaries, 0.1).tolist()
[0.5, 1.5, 2.5, 3.5, 4.5]
>>> float(golder_stability_index(0.0, 0.1))
3.0

Source code in src/plumax/gauss_plume/_src/mo_dispersion.py

def golder_stability_index(
    inverse_obukhov_length: Float[Array, "*shape"] | float, z0: float
) -> Float[Array, "*shape"]:
    """Continuous PG class index from 1/L and z₀, piecewise linear in 1/L.

    s(1/L) interpolates linearly through the knots (1/L, s) =
    ((1/L)_AB − Δ_A, 0), ((1/L)_AB, 0.5), ((1/L)_BC, 1.5), ((1/L)_CD, 2.5),
    (0, 3), ((1/L)_DE, 3.5), ((1/L)_EF, 4.5), ((1/L)_EF + Δ_F, 5), and is
    clamped to [0, 5] outside them, with Δ_A = ((1/L)_BC − (1/L)_AB)/2 and
    Δ_F = ((1/L)_EF − (1/L)_DE)/2 (the half-width of the neighbouring class).

    This is a plumax modelling choice, not part of Golder (1972): Golder gives
    only the class boundaries. Each boundary maps to the half-integer between
    its two classes, so ``round`` of the index (ties to the less stable class)
    is the hard class of :func:`golder_stability_class`, and neutral maps to
    exactly 3 (class D). Its derivative in 1/L is non-zero everywhere except
    deep inside A (1/L < (1/L)_AB − Δ_A) and F (1/L > (1/L)_EF + Δ_F).

    Args:
        inverse_obukhov_length: 1/L [m⁻¹]; zero is neutral, negative unstable.
        z0: Roughness length z₀ > 0 [m] (a scalar).

    Returns:
        The class index s ∈ [0, 5], 0 for A and 5 for F.

    Examples:
        Boundaries land on half-integers and neutral on D:

        >>> import jax.numpy as jnp
        >>> boundaries = golder_class_boundaries(0.1)
        >>> golder_stability_index(boundaries, 0.1).tolist()
        [0.5, 1.5, 2.5, 3.5, 4.5]
        >>> float(golder_stability_index(0.0, 0.1))
        3.0
    """
    inverse_obukhov_length = jnp.asarray(inverse_obukhov_length)
    boundaries = golder_class_boundaries(z0)
    # Δ_A = ((1/L)_BC − (1/L)_AB) / 2
    half_width_a = 0.5 * (boundaries[1] - boundaries[0])
    # Δ_F = ((1/L)_EF − (1/L)_DE) / 2
    half_width_f = 0.5 * (boundaries[4] - boundaries[3])
    neutral = jnp.zeros((1,), dtype=boundaries.dtype)
    knots = jnp.concatenate(
        [
            (boundaries[0] - half_width_a)[None],
            boundaries[:3],
            neutral,
            boundaries[3:],
            (boundaries[4] + half_width_f)[None],
        ]
    )
    index_knots = jnp.asarray(_INDEX_KNOTS, dtype=boundaries.dtype)
    return jnp.interp(inverse_obukhov_length, knots, index_knots)

Monin–Obukhov dispersion¶

mo_dispersion ¶

GolderDispersion ¶

from_surface_layer classmethod ¶

golder_class_boundaries ¶

golder_stability_class ¶

golder_stability_index ¶

from_surface_layer `classmethod` ¶