"""Atmospheric stability as the inverse Obukhov length ``1/L`` in ``m^-1``.
Positive values are stable. :class:`KellyGryning` generates it for every
analytical wind resource, and :class:`Charnock` gives an offshore roughness that
grows with the wind.
References
----------
Kelly, M. and Gryning, S.-E. (2010), Boundary-Layer Meteorol. 136, 377-390.
https://doi.org/10.1007/s10546-010-9509-9
Högström, U. (1988), Boundary-Layer Meteorol. 42, 55-78.
https://doi.org/10.1007/BF00119875
Charnock, H. (1955), Q. J. R. Meteorol. Soc. 81, 639-644.
https://doi.org/10.1002/qj.49708135027
Garratt, J. R. (1977), Mon. Weather Rev. 105, 915-929.
https://doi.org/10.1175/1520-0493(1977)105<0915:RODCOO>2.0.CO;2
Kelly, M. (2021), Boundary-Layer Meteorol. 179, 163-168.
https://doi.org/10.1007/s10546-020-00592-z
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import TYPE_CHECKING
import numpy as np
from .._repr import _FieldRepr
from .._validation import (
_validate_finite,
_validate_positive,
_validate_real,
_validate_real_array,
)
if TYPE_CHECKING:
from numpy.typing import ArrayLike
# Högström (1988) gradients, phi_m = (1 - 19.3 z/L)^(-1/4) unstable and
# 1 + 4.8 z/L stable, with phi_h = 0.95 + 7.8 z/L stable.
_PSI_GAMMA = 19.3
_PSI_STABLE = 4.8
_PHI_HEAT = (0.95, 7.8)
_CRITICAL_RICHARDSON = 0.25
def _stable_limit() -> float:
"""Return ``z/L`` where the gradient Richardson number reaches 0.25.
``Ri = zeta phi_h / phi_m^2`` with Högström's gradients is a quadratic in
``zeta`` at the critical value, and its positive root is the limit.
"""
a, b = _PHI_HEAT
quadratic = b - _CRITICAL_RICHARDSON * _PSI_STABLE**2
linear = a - 2.0 * _CRITICAL_RICHARDSON * _PSI_STABLE
constant = -_CRITICAL_RICHARDSON
return (-linear + np.sqrt(linear**2 - 4.0 * quadratic * constant)) / (
2.0 * quadratic
)
#: ``z/L`` beyond which the stable wind gradient is held. Here the gradient
#: Richardson number with Högström's gradients reaches the critical 0.25, so
#: the local shear stops growing and turbulence falls only slowly beyond it.
STABLE_LIMIT = float(_stable_limit())
# Gradient phi_m held beyond the limit, so 1 - phi_m is the slope of Psi.
_HELD_SLOPE = _PSI_STABLE * STABLE_LIMIT
[docs]
def stability_correction(zeta: ArrayLike) -> np.ndarray:
"""Return the diabatic correction ``Psi(z / L)`` to the log wind profile.
``zeta`` is ``z / L``, positive stable. ``Psi`` is zero at neutral, positive
when unstable and negative when stable. It is the integral
``Psi = int (1 - phi_m) d ln(z/L)`` of :func:`stability_gradient`, so
``Psi = -4.8 z/L`` up to :data:`~oowl.resources.profile.STABLE_LIMIT` and falls logarithmically
beyond it, where the gradient is held.
"""
z = _validate_real_array("zeta", zeta)
if not np.all(np.isfinite(z)):
raise ValueError("zeta must be finite")
result = np.empty_like(z)
unstable = z < 0.0
if np.any(unstable):
x = (1.0 - _PSI_GAMMA * z[unstable]) ** 0.25
result[unstable] = (
2.0 * np.log(0.5 * (1.0 + x))
+ np.log(0.5 * (1.0 + x**2))
- 2.0 * np.arctan(x)
+ 0.5 * np.pi
)
stable = z[~unstable]
result_stable = -_PSI_STABLE * stable
beyond = stable > STABLE_LIMIT
result_stable[beyond] = -_HELD_SLOPE * (1.0 + np.log(stable[beyond] / STABLE_LIMIT))
result[~unstable] = result_stable
return result
[docs]
def stability_gradient(zeta: ArrayLike) -> np.ndarray:
"""Return the dimensionless wind gradient ``phi_m(z / L)``.
Högström's ``(1 - 19.3 z/L)^(-1/4)`` when unstable and ``1 + 4.8 z/L`` when
stable, held beyond :data:`~oowl.resources.profile.STABLE_LIMIT`. It is one at neutral, below one
when unstable and above one when stable.
"""
z = _validate_real_array("zeta", zeta)
if not np.all(np.isfinite(z)):
raise ValueError("zeta must be finite")
unstable = (1.0 - _PSI_GAMMA * np.minimum(z, 0.0)) ** -0.25
return np.where(z < 0.0, unstable, 1.0 + _PSI_STABLE * np.minimum(z, STABLE_LIMIT))
def _profile(
height: float, roughness_length: ArrayLike, inverse_obukhov_length: ArrayLike
) -> np.ndarray:
"""Return ``ln(z/z0) - Psi(z/L) + Psi(z0/L)`` at ``height``.
It is the integral of ``phi_m`` over ``ln z`` from ``z0`` to ``z``, so it
stays positive however unstable the case. The surface term ``Psi(z0/L)``
is usually negligible and matters over rough ground in strong convection.
"""
inverse_length = np.asarray(inverse_obukhov_length, dtype=float)
if not np.all(np.isfinite(inverse_length)):
raise ValueError("inverse_obukhov_length must be finite")
ratio = _checked_surface(height, roughness_length)
zeta = float(height) * inverse_length
return (
np.log(ratio) - stability_correction(zeta) + stability_correction(zeta / ratio)
)
def _checked_wind_speed(wind_speed) -> np.ndarray:
u = _validate_real_array("wind_speed", wind_speed)
if not np.all(np.isfinite(u)) or np.any(u < 0):
raise ValueError("wind_speed must contain finite values >= 0")
return u
def _checked_positive_speed(wind_speed) -> np.ndarray:
u = _validate_real_array("wind_speed", wind_speed)
if not np.all(np.isfinite(u)) or np.any(u <= 0):
raise ValueError("wind_speed must contain only positive finite values")
return u
def _checked_height(height) -> float:
height = _validate_real("height", height)
if not np.isfinite(height) or height <= 0.0:
raise ValueError("height must be positive and finite")
return height
def _checked_quantiles(values, shape: tuple[int, ...]) -> np.ndarray:
q = _validate_real_array("quantiles", values)
if q.shape != shape:
raise ValueError("quantiles must have the same shape as wind_speed")
if not np.all(np.isfinite(q)) or np.any((q < 0) | (q > 1)):
raise ValueError("quantiles must contain finite values in [0, 1]")
return q
_GRAVITY = 9.81
_VON_KARMAN = 0.4
[docs]
def friction_velocity(
wind_speed: ArrayLike, height: float, roughness_length: ArrayLike
) -> np.ndarray:
"""Return the neutral friction velocity ``u_star`` in ``m/s``.
``u_star = kappa u(z) / ln(z / z0)`` with no stability correction.
``roughness_length`` is one value or one per wind speed.
"""
u = _checked_wind_speed(wind_speed)
return _VON_KARMAN * u / np.log(_checked_surface(height, roughness_length))
def _checked_surface(height, roughness_length) -> np.ndarray:
"""Return ``height / roughness_length`` for a usable surface."""
height = _validate_real("height", height)
roughness_length = _validate_real_array("roughness_length", roughness_length)
if not np.all(np.isfinite(roughness_length)) or np.any(roughness_length <= 0):
raise ValueError("roughness_length must be positive and finite")
if not np.isfinite(height) or np.any(height <= roughness_length):
raise ValueError("height must be finite and above roughness_length")
return height / roughness_length
[docs]
@dataclass(frozen=True, repr=False)
class Charnock(_FieldRepr):
"""Offshore roughness that grows with the wind, ``z0 = A_c u_star^2 / g``.
Waves roughen the water surface as the wind rises (Charnock, 1955). With
the neutral log law this is solved exactly for the wind speed ``U`` at
height ``z`` through the lower branch of the Lambert W function (Kelly,
2021),
.. math::
z_0 = \\frac{\\kappa^2 A_c}{4} \\frac{U^2 / g}
{\\left[W_{-1}\\!\\left(-\\frac{\\kappa U}{2}
\\sqrt{\\frac{A_c}{g z}}\\right)\\right]^2}.
``constant`` is the Charnock constant ``A_c``. The default ``0.0144`` is
the ocean fit of Garratt (1977), and ESDU International (1982) recommends
``0.0162``. Pass an instance as a resource's ``roughness_length``.
"""
#: Dimensionless Charnock constant ``A_c``.
constant: float = 0.0144
def __post_init__(self):
constant = _validate_positive("constant", self.constant)
object.__setattr__(self, "constant", constant)
[docs]
def roughness_length(self, wind_speed: ArrayLike, height: float) -> np.ndarray:
"""Return ``z0`` in meters for each positive wind speed at ``height``."""
from scipy.special import lambertw
u = _checked_wind_speed(wind_speed)
if np.any(u <= 0):
raise ValueError("wind_speed must be positive to define a roughness")
height = _checked_height(height)
argument = -0.5 * _VON_KARMAN * u * np.sqrt(self.constant / (_GRAVITY * height))
# The log law and Charnock meet only while the argument stays above -1/e.
if np.any(argument < -np.exp(-1.0)):
raise ValueError(
"wind_speed is too strong for a Charnock roughness at height"
)
branch = lambertw(argument, -1).real
return 0.25 * _VON_KARMAN**2 * self.constant * u**2 / (_GRAVITY * branch**2)
def _resolve_roughness(roughness_length):
"""Return a validated roughness, a length or a Charnock model."""
if isinstance(roughness_length, Charnock):
return roughness_length
if isinstance(roughness_length, (str, bytes, bool, np.bool_)):
value = None
else:
try:
value = float(roughness_length)
except (TypeError, ValueError):
value = None
if value is None:
raise TypeError(
"roughness_length must be a positive length or a Charnock model, "
f"not {type(roughness_length).__name__}"
)
if not np.isfinite(value) or value <= 0:
raise ValueError("roughness_length must be positive and finite")
return value
# Kelly and Gryning (2010) universal constants C+ and C- of eq. 5.
_KG_STABLE_CONSTANT = 1.0
_KG_UNSTABLE_CONSTANT = 3.0
# Surface-layer scaling of 1/L with u_star at fixed thermal forcing.
_KG_STABLE_EXPONENT = 2.0
_KG_UNSTABLE_EXPONENT = 3.0
# Range searched when fitting the exponents to a record.
_EXPONENT_RANGE = (0.0, 6.0)
# The largest probability below one, so the tails stay finite.
_LAST_PROBABILITY = float(np.nextafter(1.0, 0.0))
[docs]
@dataclass(frozen=True, repr=False, kw_only=True)
class KellyGryning(_FieldRepr):
"""Stability climate of Kelly and Gryning (2010), conditioned on the wind.
Each side of the distribution of ``1/L`` is a stretched exponential,
.. math::
P(L^{-1}) = n_\\pm \\frac{C_\\pm}{\\sigma_\\pm \\Gamma(5/2)}
\\exp\\left[-\\left(\\frac{C_\\pm |L^{-1}|}{\\sigma_\\pm}\\right)^{2/3}\\right],
with the universal ``C+ = 1`` and ``C- = 3``. ``n+`` and ``n-`` are the
stable and unstable fractions and ``sigma+`` and ``sigma-`` the scales of
``1/L`` on each side. Kelly and Gryning derive the scales from heat-flux
statistics as ``sigma = (g / T0) std(H) / <u_star>^3``. :meth:`fit`
estimates them from a record of ``1/L`` by maximum likelihood, which at
Høvsøre gives an unstable scale a third of the heat-flux estimate.
Conditioning on the wind speed is an OOWL extension. The scales hold at the
site's mean friction velocity ``reference_friction_velocity``. A case whose
neutral friction velocity is ``u_star`` uses
``sigma+ (u_ref / u_star)^p+`` and ``sigma- (u_ref / u_star)^p-``. A stable
surface layer keeps its temperature scale, so ``1/L`` falls as
``u_star^-2``, while convection keeps its heat flux, so ``1/L`` falls as
``u_star^-3``. These are the default exponents, and :meth:`fit` can
estimate them instead.
``unstable_fraction`` is one number, or one per friction velocity in
``fraction_friction_velocity``, interpolated linearly in the neutral
``u_star`` and held beyond the ends. Sites often turn unstable in strong
wind, when cold air crosses warmer water, and at low wind by day over land.
The defaults are round offshore magnitudes, not a site calibration.
"""
#: Fraction ``n-`` of unstable cases, one or one per friction velocity.
unstable_fraction: float | tuple[float, ...] = 0.6
#: Stable scale ``sigma+`` in ``m^-1`` at ``reference_friction_velocity``.
stable_scale: float = 0.015
#: Unstable scale ``sigma-`` in ``m^-1`` at ``reference_friction_velocity``.
unstable_scale: float = 0.03
#: Mean friction velocity in ``m/s`` at which the scales hold.
reference_friction_velocity: float = 0.3
#: Exponent ``p+`` of ``u_ref / u_star`` scaling the stable side.
stable_exponent: float = _KG_STABLE_EXPONENT
#: Exponent ``p-`` of ``u_ref / u_star`` scaling the unstable side.
unstable_exponent: float = _KG_UNSTABLE_EXPONENT
#: Neutral friction velocities in ``m/s`` of a table of unstable fractions.
fraction_friction_velocity: tuple[float, ...] | None = None
def __post_init__(self):
self._resolve_fractions()
for name in (
"stable_scale",
"unstable_scale",
"reference_friction_velocity",
"stable_exponent",
"unstable_exponent",
):
value = _validate_finite(name, getattr(self, name))
object.__setattr__(self, name, value)
for name in ("stable_scale", "unstable_scale", "reference_friction_velocity"):
if getattr(self, name) <= 0.0:
raise ValueError(f"{name} must be positive")
for name in ("stable_exponent", "unstable_exponent"):
if getattr(self, name) < 0.0:
raise ValueError(f"{name} must not be negative")
def _resolve_fractions(self) -> None:
"""Store the fraction as a float, or as a table over friction velocity."""
fraction = _validate_real_array("unstable_fraction", self.unstable_fraction)
nodes = self.fraction_friction_velocity
if not np.all(np.isfinite(fraction)) or np.any(
(fraction < 0.0) | (fraction > 1.0)
):
raise ValueError("unstable_fraction must lie in [0, 1]")
if fraction.ndim == 0:
if nodes is not None:
raise ValueError(
"fraction_friction_velocity needs one unstable_fraction each"
)
object.__setattr__(self, "unstable_fraction", float(fraction))
return
nodes = _validate_real_array(
"fraction_friction_velocity", [] if nodes is None else nodes
)
if fraction.ndim != 1 or nodes.shape != fraction.shape or len(nodes) < 2:
raise ValueError(
"a table of unstable_fraction needs one fraction_friction_velocity "
"each, at least two"
)
if (
not np.all(np.isfinite(nodes))
or nodes[0] <= 0
or np.any(np.diff(nodes) <= 0)
):
raise ValueError(
"fraction_friction_velocity must be positive and increasing"
)
object.__setattr__(self, "unstable_fraction", tuple(fraction.tolist()))
object.__setattr__(self, "fraction_friction_velocity", tuple(nodes.tolist()))
[docs]
def fractions(
self, wind_speed: ArrayLike, *, height: float, roughness_length: ArrayLike
) -> np.ndarray:
"""Return the unstable fraction ``n-`` at each wind speed."""
u = _checked_wind_speed(wind_speed)
if self.fraction_friction_velocity is None:
return np.full(u.shape, self.unstable_fraction)
return np.interp(
friction_velocity(u, height, roughness_length),
self.fraction_friction_velocity,
self.unstable_fraction,
)
[docs]
@classmethod
def fit(
cls,
inverse_obukhov_length: ArrayLike,
wind_speed: ArrayLike,
*,
height: float,
roughness_length: float,
fit_exponents: bool = False,
fraction_bins: int = 1,
) -> KellyGryning:
"""Return the climate that best explains a record of ``1/L``.
``inverse_obukhov_length`` in ``m^-1`` and ``wind_speed`` at ``height``
are paired records, for example ten-minute means, over a surface of
``roughness_length``. The unstable fraction is the share of negative
``1/L`` and ``reference_friction_velocity`` the mean neutral friction
velocity. Given the exponents, each scale maximizes the likelihood in
closed form,
.. math::
\\sigma_\\pm = \\left[\\frac{2}{3N_\\pm}\\sum_i
\\left(\\frac{C_\\pm |L_i^{-1}|}{r_i^{p_\\pm}}\\right)^{2/3}
\\right]^{3/2},
\\qquad r_i = \\frac{u_{\\mathrm{ref}}}{u_{*,i}}.
With ``fit_exponents`` the exponents also maximize the likelihood,
searched between 0 and 6. Both sides need records.
``fraction_bins`` above one splits the records into that many
equal-count bins of friction velocity and tabulates the unstable share
of each at its mean ``u_star``. Records of one friction velocity share
the bin of its first record in sorted order, so a record of few
distinct wind speeds, such as speeds stored at a fixed resolution, gives
fewer bins than requested, and a single bin gives one fraction. The
scales describe each side given its sign, so they do not change.
"""
from scipy.optimize import minimize_scalar
inverse_length = np.asarray(inverse_obukhov_length, dtype=float)
speed = _checked_wind_speed(wind_speed)
if inverse_length.ndim != 1 or inverse_length.shape != speed.shape:
raise ValueError("records must be 1-D and of equal length")
if not np.all(np.isfinite(inverse_length)):
raise ValueError("inverse_obukhov_length must be finite")
if np.any(speed <= 0):
raise ValueError("wind_speed must be positive")
u_star = friction_velocity(speed, height, roughness_length)
reference = float(np.mean(u_star))
log_ratio = np.log(reference / u_star)
unstable = inverse_length < 0.0
if unstable.all() or not unstable.any():
raise ValueError("fitting needs both stable and unstable records")
if isinstance(fraction_bins, (bool, np.bool_)) or not isinstance(
fraction_bins, (int, float, np.integer, np.floating)
):
raise TypeError(
f"fraction_bins must be an integer, not {type(fraction_bins).__name__}"
)
if not isinstance(fraction_bins, (int, np.integer)) or not (
1 <= fraction_bins <= len(speed)
):
raise ValueError("fraction_bins must be an integer from 1 to the records")
fraction: float | tuple[float, ...] = float(unstable.mean())
nodes = None
if fraction_bins > 1:
order = np.argsort(u_star)
sorted_speed = u_star[order]
# Each record joins the equal-count bin of its value's first
# occurrence, so ties share a bin and empty bins drop out.
group_starts = np.flatnonzero(
np.r_[True, sorted_speed[1:] != sorted_speed[:-1]]
)
size, extra = divmod(len(order), fraction_bins)
counts = np.full(fraction_bins - 1, size)
counts[:extra] += 1
starts = np.r_[group_starts, len(order)][
np.searchsorted(group_starts, np.cumsum(counts))
]
bins = np.split(order, np.unique(starts[starts < len(order)]))
if len(bins) > 1:
fraction = tuple(float(unstable[i].mean()) for i in bins)
nodes = tuple(float(u_star[i].mean()) for i in bins)
def side(mask, constant, exponent):
magnitude = constant * np.abs(inverse_length[mask])
reduced = magnitude * np.exp(-exponent * log_ratio[mask])
return (2.0 / (3.0 * mask.sum()) * np.sum(reduced ** (2.0 / 3.0))) ** 1.5
def exponent_for(mask, constant, default):
if not fit_exponents:
return default
# Negative log-likelihood with the scale at its optimum for p.
def objective(exponent):
scale = side(mask, constant, exponent)
return mask.sum() * np.log(scale) + exponent * log_ratio[mask].sum()
return float(
minimize_scalar(objective, bounds=_EXPONENT_RANGE, method="bounded").x
)
stable_exponent = exponent_for(
~unstable, _KG_STABLE_CONSTANT, _KG_STABLE_EXPONENT
)
unstable_exponent = exponent_for(
unstable, _KG_UNSTABLE_CONSTANT, _KG_UNSTABLE_EXPONENT
)
return cls(
unstable_fraction=fraction,
fraction_friction_velocity=nodes,
stable_scale=float(side(~unstable, _KG_STABLE_CONSTANT, stable_exponent)),
unstable_scale=float(
side(unstable, _KG_UNSTABLE_CONSTANT, unstable_exponent)
),
reference_friction_velocity=reference,
stable_exponent=stable_exponent,
unstable_exponent=unstable_exponent,
)
[docs]
def scales(
self, wind_speed: ArrayLike, *, height: float, roughness_length: ArrayLike
) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(sigma+, sigma-)`` in ``m^-1`` at each wind speed."""
u = _checked_wind_speed(wind_speed)
if np.any(u <= 0):
raise ValueError("wind_speed must be positive to define a stability")
ratio = self.reference_friction_velocity / friction_velocity(
u, height, roughness_length
)
return (
self.stable_scale * ratio**self.stable_exponent,
self.unstable_scale * ratio**self.unstable_exponent,
)
[docs]
def sample(
self,
wind_speed: ArrayLike,
quantiles: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
) -> np.ndarray:
"""Return ``1/L`` in ``m^-1`` at each wind speed and quantile, exactly.
Quantiles below the unstable fraction are unstable and the rest stable.
On each side ``(C |1/L| / sigma)^(2/3)`` follows a gamma law of shape
``3/2``, whose inverse gives the tail.
"""
surface = {"height": height, "roughness_length": roughness_length}
stable_scale, unstable_scale = self.scales(wind_speed, **surface)
q = _checked_quantiles(quantiles, np.shape(stable_scale))
fraction = self.fractions(wind_speed, **surface)
return self._inverse(q, fraction, stable_scale, unstable_scale)
def _quantile_table(
self,
wind_speed: np.ndarray,
quantiles: np.ndarray,
*,
height: float,
roughness_length: np.ndarray,
) -> np.ndarray:
"""Return ``1/L`` with one row per wind speed and one column per quantile.
As :meth:`sample` at every pair. A single unstable fraction gives every
speed the same depths, so they are inverted only once.
"""
surface = {"height": height, "roughness_length": roughness_length}
stable_scale, unstable_scale = self.scales(wind_speed, **surface)
fraction = (
self.unstable_fraction
if self.fraction_friction_velocity is None
else self.fractions(wind_speed, **surface)[:, None]
)
return self._inverse(
quantiles, fraction, stable_scale[:, None], unstable_scale[:, None]
)
@staticmethod
def _inverse(q, fraction, stable_scale, unstable_scale) -> np.ndarray:
"""Return ``1/L`` at quantiles ``q``, all arguments broadcast together."""
from scipy.special import gammaincinv
unstable = q < fraction
# Depth into each side, from zero at neutral to one in the far tail.
depth = np.zeros(unstable.shape)
np.divide(fraction - q, fraction, out=depth, where=unstable)
np.divide(
q - fraction, 1.0 - fraction, out=depth, where=~unstable & (fraction < 1.0)
)
# Quadratures repeat the same depths at every speed, so each distinct
# depth is inverted once.
levels, index = np.unique(
np.clip(depth, 0.0, _LAST_PROBABILITY).ravel(), return_inverse=True
)
magnitude = (gammaincinv(1.5, levels) ** 1.5)[index].reshape(depth.shape)
return np.where(
unstable,
-unstable_scale / _KG_UNSTABLE_CONSTANT * magnitude,
stable_scale / _KG_STABLE_CONSTANT * magnitude,
)
[docs]
def bounds(
self,
wind_speed: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
coverage: float = 0.98,
) -> tuple[np.ndarray, np.ndarray]:
"""Return the central ``coverage`` interval of ``1/L``, exactly."""
if not 0.0 < coverage < 1.0:
raise ValueError("coverage must lie in (0, 1)")
shape = np.shape(wind_speed)
surface = {"height": height, "roughness_length": roughness_length}
tail = 0.5 * (1.0 - coverage)
low = self.sample(wind_speed, np.full(shape, tail), **surface)
high = self.sample(wind_speed, np.full(shape, 1.0 - tail), **surface)
return low, high
def _resolve_model(stability):
"""Return the stability model accepted by a wind resource."""
if stability is None:
return KellyGryning()
if isinstance(stability, KellyGryning):
return stability
raise TypeError("stability must be a KellyGryning model")
