"""soil_ground_heat_flux.py
====================================================
Python utilities implementing the equations from:
Wang, Z.-H., & Bou-Zeid, E. (2012). *A novel approach for the estimation of
soil ground heat flux*. *Agricultural and Forest Meteorology*, 154-155,
214-221.
All equations appearing in that paper are reproduced as vectorised Python
functions, together with a few convenience helpers. The library is fully
NumPy-aware and can be used on scalars or on time-series arrays.
Units
-----
All functions assume SI units **throughout**:
* depth ``z`` ― m (positive downward)
* time ``t`` ― s
* temperature ``T`` ― K (or °C provided the 0-offset is handled
consistently)
* heat flux ``G, H, LE`` ― W m-2
* radiative flux ``Rn`` ― W m-2
* soil properties
* thermal conductivity ``k`` ― W m-1 K-1
* heat capacity ``rho_c`` ― J m-3 K-1
* thermal diffusivity ``kappa = k / rho_c`` ― m² s-1
Dependencies
------------
>>> pip install numpy scipy
Example
-------
>>> import numpy as np, wang_and_bouzeid as sghf
>>> # 30-minute series (dt = 1800 s) of flux-plate measurements at z = 0.08 m
>>> Gz = np.loadtxt('Gz_8cm.txt')
>>> G0 = sghf.estimate_G0_from_Gz(Gz, z_r=0.08, kappa=0.7e-6, dt=1800)
"""
from __future__ import annotations
import math
from typing import Callable, Sequence
import numpy as np
from scipy.special import erfc, gammaincc
_trapz = getattr(np, "trapezoid", None) or getattr(np, "trapz")
__all__ = [
# Energy balance & residuals
"energy_balance_residual",
"surface_energy_residual",
# Conventional ground-heat-flux estimator (Eq. 2)
"ground_heat_flux_conventional",
# Heat-conduction fundamentals
"green_function_temperature",
"temperature_convolution_solution",
"soil_heat_flux_from_G0",
"estimate_G0_from_Gz",
# Sinusoidal analytical solutions (Eqs. 13–15)
"sinusoidal_boundary_flux",
"soil_temperature_sinusoidal",
"soil_heat_flux_sinusoidal",
# Soil-property parameterisations (Eqs. 16–18)
"heat_capacity_moist_soil",
"pf_from_theta",
"thermal_conductivity_moist_soil",
"thermal_diffusivity",
]
# -----------------------------------------------------------------------------
# 1. Energy-balance bookkeeping (Eqs. 1 & 19)
# -----------------------------------------------------------------------------
[docs]
def energy_balance_residual(
Rn: float | np.ndarray,
H: float | np.ndarray,
LE: float | np.ndarray,
G0: float | np.ndarray,
) -> float | np.ndarray:
"""
Compute the **closure residual** of the surface energy balance (SEB).
The classical SEB for land–atmosphere exchange is
.. math::
R_n - G_0 \;=\; H + LE + \varepsilon ,
where the residual term :math:`\varepsilon` quantifies the lack of
closure. Rearranging gives
.. math::
\varepsilon \;=\; R_n - G_0 - H - LE ,
which is what this helper returns.
Parameters
----------
Rn : float or array_like
Net radiation *Rₙ* (W m⁻²). Positive downward.
H : float or array_like
Sensible heat flux *H* (W m⁻²). Positive upward (atmosphere ← surface).
LE : float or array_like
Latent heat flux *LE* (W m⁻²). Positive upward.
G0 : float or array_like
Ground (soil) heat flux *G₀* (W m⁻²). Positive downward
(into the soil). Some authors use the opposite sign convention;
ensure consistency with *Rn*.
Returns
-------
float or ndarray
Energy‐balance residual :math:`\varepsilon` (W m⁻²) with the
broadcast shape of the inputs.
**Positive** values indicate missing energy
(surface gains > turbulent + ground fluxes),
whereas **negative** values mean an apparent energy surplus.
Notes
-----
* **Broadcasting** – All inputs are treated with NumPy broadcasting,
allowing scalars, 1-D arrays, or DataFrame columns.
* **Closure diagnostics** – The residual can be summarised as a mean
bias or expressed as a relative closure fraction
``1 − ε / (Rn − G0)``.
* **Typical magnitudes** – Eddy‐covariance towers often report
10–30 % non-closure. Persistently large |ε| values may indicate
advective transport, sensor tilt, or data‐processing issues.
Examples
--------
>>> ε = energy_balance_residual(
... Rn=df["Rn"], H=df["H"], LE=df["LE"], G0=df["G"]
... )
>>> ε.describe(percentiles=[0.05, 0.5, 0.95])
count 17280.000000
mean -9.73
std 34.51
5% -58.42
50% -8.11
95% 39.12
Name: residual, dtype: float64
Plot daily average residuals:
>>> ε.resample("D").mean().plot(marker="o")
>>> plt.axhline(0, color="k", lw=0.8)
>>> plt.ylabel("Energy balance residual (W m$^{-2}$)")
>>> plt.title("Daily SEB closure")
"""
return Rn - G0 - H - LE
# Alias used later in the module.
[docs]
surface_energy_residual = energy_balance_residual
# -----------------------------------------------------------------------------
# 2. Conventional ground-heat-flux estimator (gradient + calorimetry) – Eq. 2
# -----------------------------------------------------------------------------
[docs]
def ground_heat_flux_conventional(
k: float,
dT_dz_at_zr: float,
rho_c: float,
dT_dt_profile: Sequence[float],
z_profile: Sequence[float],
) -> float:
"""Conventional estimate of *surface* ground heat flux, Eq. (2).
Parameters
----------
k
Effective soil thermal conductivity **(W m-1 K-1)**.
dT_dz_at_zr
Vertical temperature gradient evaluated at the flux-plate depth ``z_r``
**(K m-1)**. A *negative* gradient means temperature decreases with
depth.
rho_c
Volumetric heat capacity of the soil **(J m-3 K-1)**.
dT_dt_profile
Time derivatives ∂T/∂t for *each* node between the surface and ``z_r``
**(K s-1)**. Any iterable (list, ndarray, …). Must align with
``z_profile``.
z_profile
Depth of each node in the temperature profile **(m)**. Increasing,
positive downward, **excluding** the surface (z = 0) but *including*
``z_r`` (last element).
Returns
-------
G0 : float
Ground heat flux at the surface **(W m-2)**. Positive *into* the soil.
"""
# Fourier conduction term (gradient method)
G_conduction = -k * dT_dz_at_zr
# Heat-storage (calorimetry) – numerical integration by trapezoid
z = np.asarray(z_profile)
dT_dt = np.asarray(dT_dt_profile)
if z.shape != dT_dt.shape:
raise ValueError("z_profile and dT_dt_profile must have same length")
# Integration bounds: surface (0) to z_r (last node) – prepend surface
z_nodes = np.concatenate(([0.0], z))
dT_dt_nodes = np.concatenate(([dT_dt[0]], dT_dt))
storage = _trapz(dT_dt_nodes, x=z_nodes)
G_storage = rho_c * storage
return G_conduction + G_storage # type: ignore
# -----------------------------------------------------------------------------
# 3. Heat-conduction fundamentals (Eqs. 3–12)
# -----------------------------------------------------------------------------
[docs]
def green_function_temperature(
z: float,
t: float,
kappa: float,
) -> float:
"""
Green-function solution :math:`g_z(t)` for the **one‐dimensional,
semi-infinite heat equation** with a unit surface‐flux impulse at
:math:`t=0,\,z=0`.
The governing partial differential equation is
.. math::
\\frac{\\partial T}{\\partial t}
\;=\;
\\kappa\\,\\frac{\\partial^{2} T}{\\partial z^{2}},
\\qquad z \\ge 0,\\; t > 0,
with initial condition :math:`T(z,0)=0` and boundary condition
corresponding to a Dirac δ–heat pulse applied at the surface
(:math:`z=0`). The resulting Green function (Carslaw & Jaeger,
1959, Eq. 7) is
.. math::
g_z(t)
\;=\;
\\frac{2}{\\sqrt{\\pi}}
\\,\\sqrt{\\kappa t}\;
\\exp\\!\Bigl[-\\frac{z^{2}}{4\\kappa t}\\Bigr]
\;-\;
z\,\\operatorname{erfc}\\!
\\Bigl[\\frac{z}{2\\sqrt{\\kappa t}}\\Bigr],
\\qquad t>0,
and :math:`g_z(t)=0` for :math:`t\\le 0` (causality).
Parameters
----------
z : float
Depth below the surface (m, positive downward). Must be
non-negative.
t : float
Time since the surface impulse (s). Values :math:`t \\le 0`
return 0 by definition.
kappa : float
Thermal diffusivity :math:`\\kappa` of the half-space
(m² s⁻¹).
Returns
-------
float
Green-function value :math:`g_z(t)` (units **m**, because the
solution integrates heat-flux density with respect to depth to
yield temperature).
Notes
-----
* **Causality check** – If ``t`` ≤ 0 the function short-circuits and
returns 0.0.
* **Vectorisation** – For vector or array input use
:func:`numpy.vectorize` or wrap the function in
:func:`numpy.frompyfunc`; the core implementation is scalar for
numerical clarity.
* **Usage** – ``g_z(t)`` can be convolved with an arbitrary surface
heat-flux time series ``q₀(t)`` to obtain temperature at depth
via
.. math::
T(z,t) \;=\; \\int_{0}^{t} g_z(t-τ)\,q_0(τ)\;\\mathrm dτ .
References
----------
Carslaw, H. S., & Jaeger, J. C. (1959).
*Conduction of Heat in Solids* (2nd ed., p. 100).
Oxford University Press.
Examples
--------
>>> g = green_function_temperature(z=0.05, t=3_600, kappa=1.4e-7)
>>> print(f"g(5 cm, 1 h) = {g:.3e} m")
g(5 cm, 1 h) = 7.42e-04 m
"""
if t <= 0:
return 0.0
return 2.0 / math.sqrt(math.pi) * math.sqrt(kappa * t) * math.exp(
-(z**2) / (4.0 * kappa * t)
) - z * erfc(z / (2.0 * math.sqrt(kappa * t)))
[docs]
def temperature_convolution_solution(
z: float, t_series: np.ndarray, f_series: np.ndarray, kappa: float, Ti: float = 0.0
) -> np.ndarray:
"""Temperature time-series at depth *z* via Duhamel convolution (Eq. 6).
``T(z,t) = Ti + ∫ f(t-τ) d g_z(τ)``
The integral becomes a discrete convolution where *f* is the boundary
heat-flux series (W m-2 → ∂T/∂z via Fourier).
"""
if t_series.ndim != 1 or f_series.ndim != 1:
raise ValueError("t_series and f_series must be 1-D arrays of equal length")
if t_series.size != f_series.size:
raise ValueError("t_series and f_series must be the same length")
dt = np.diff(t_series)
if not np.allclose(dt, dt[0]):
raise ValueError("Time vector must be uniformly spaced")
dt = dt[0]
g = np.array([green_function_temperature(z, t, kappa) for t in t_series])
dg = np.diff(g, prepend=0.0) # discrete derivative → Stieltjes measure
# Convolution implementation: cumulative sum of f * dg (causal)
T = Ti + np.cumsum(f_series[::-1] * dg)[::-1] # reversed for (t-τ)
return T
[docs]
def soil_heat_flux_from_G0(
z: float, t_series: np.ndarray, G0_series: np.ndarray, kappa: float
) -> np.ndarray:
"""Compute *G(z,t)* from a known surface flux series *G0* (Eq. 9)."""
if t_series.ndim != 1 or G0_series.ndim != 1:
raise ValueError("t_series and G0_series must be 1-D")
if t_series.size != G0_series.size:
raise ValueError("t_series and G0_series must align")
dt = np.diff(t_series)
if not np.allclose(dt, dt[0]):
raise ValueError("Time vector must be uniformly spaced")
dt = dt[0]
# Build F_z(t) = erfc(z / 2√(κ t)) (with F_z(0) = 0 by limit)
with np.errstate(divide="ignore", invalid="ignore"):
Fz = erfc(z / (2.0 * np.sqrt(kappa * t_series)))
Fz[0] = 0.0
dF = np.diff(Fz, prepend=0.0)
# Convolution similar to temperature_convolution_solution
Gz = np.cumsum(G0_series[::-1] * dF)[::-1]
return Gz
[docs]
def estimate_G0_from_Gz(
Gz_series: np.ndarray, z_r: float, kappa: float, dt: float
) -> np.ndarray:
"""Estimate *surface* ground heat flux *G0* from plate measurements *Gz*.
Implements discretised Eq. (11) – the recursion proposed by Wang & Bou-Zeid
(2012). Time-series must be *regularly* sampled.
Parameters
----------
Gz_series : np.ndarray
Soil heat-flux measurements at depth *z_r* **(W m-2)**.
z_r : float
Plate depth **(m)**.
kappa : float
Thermal diffusivity **(m² s-1)**.
dt : float
Sampling interval **(s)**.
Returns
-------
G0 : np.ndarray
Estimated surface heat-flux series **(W m-2)**.
"""
Gz_series = np.asarray(Gz_series, dtype=float)
n_steps = Gz_series.size
# Pre-compute ΔF_z(j) for j = 1 … n-1 (Eq. 10)
j = np.arange(n_steps) # 0 … n-1
t_j = j * dt
with np.errstate(divide="ignore", invalid="ignore"):
Fz = erfc(z_r / (2.0 * np.sqrt(kappa * t_j)))
Fz[0] = 0.0
dF = np.diff(Fz, prepend=0.0)
G0 = np.zeros_like(Gz_series)
for n in range(1, n_steps):
# J_{n-1} term (Eq. 12)
J = 0.0
for j in range(1, n):
J += 0.5 * (G0[n - j] + G0[n - j - 1]) * dF[j]
G0[n] = (2.0 * Gz_series[n] - J) / dF[1]
# By construction dF[1] > 0 (t = dt)
G0[0] = Gz_series[0] # first guess – no history available
return G0
# -----------------------------------------------------------------------------
# 4. Sinusoidal analytical solutions (Eqs. 13–15)
# -----------------------------------------------------------------------------
[docs]
def sinusoidal_boundary_flux(
t: float | np.ndarray,
A: float,
omega: float,
epsilon: float,
) -> float | np.ndarray:
"""
Evaluate a **sinusoidal surface heat-flux forcing**
.. math::
q_0(t) \;=\; A \,\sin\!\bigl(\omega t + \varepsilon\bigr),
which is commonly used as a boundary condition for analytical soil-
heat-flux solutions (see Eq. 13 of the companion text).
Parameters
----------
t : float or array_like
Time since the start of the simulation (s).
May be scalar or any NumPy-broadcastable shape; units must be
consistent with ``omega``.
A : float
Amplitude of the surface heat flux (W m⁻²). Positive **downward**
into the soil.
omega : float
Angular frequency (rad s⁻¹). For a diurnal cycle
``omega = 2 π / 86 400`` ≈ 7.272 × 10⁻⁵ rad s⁻¹.
epsilon : float
Phase shift **ε** (rad). Positive values delay the flux peak,
negative values advance it.
Returns
-------
ndarray or float
Instantaneous surface heat flux *q₀(t)* (W m⁻²) with shape given
by NumPy broadcasting of the inputs.
Notes
-----
* **Sign convention** — Positive *q₀* adds energy to the soil column; be
sure it matches the sign convention of your governing equation.
* **Vectorisation** — The implementation is a single call to
``numpy.sin`` and therefore fully vectorised.
* **Period** — The period *P* (s) of the forcing is related to
``omega`` by *P = 2 π / ω*.
Examples
--------
>>> import numpy as np, matplotlib.pyplot as plt
>>> t = np.linspace(0, 2*86400, 1_000) # 2 days
>>> q0 = sinusoidal_boundary_flux(
... t, A=120, omega=2*np.pi/86400, epsilon=0)
>>> plt.plot(t/3600, q0)
>>> plt.xlabel("Time (h)")
>>> plt.ylabel("Surface heat flux $q_0$ (W m$^{-2}$)")
>>> plt.title("Sinusoidal surface forcing (A = 120 W m⁻²)")
>>> plt.show()
"""
return A * np.sin(omega * t + epsilon)
[docs]
def soil_temperature_sinusoidal(
z: float,
t: float | np.ndarray,
A: float,
omega: float,
epsilon: float,
Ti: float,
kappa: float,
) -> float | np.ndarray:
"""
Analytical solution for **soil temperature** beneath a sinusoidally
forced surface heat‐flux boundary.
The temperature response of a semi-infinite, homogeneous soil column
to a surface heat flux
.. math::
q_0(t) \;=\; A \sin(\omega t + \varepsilon)
is (Carslaw & Jaeger 1959, Eq. 14)
.. math::
T(z,t)
\;=\;
T_i
\;+\;
\\underbrace{\\frac{A}{\\kappa\\sqrt{\\omega}}
e^{-z r}\,
\\sin\!\bigl(\\omega t + \varepsilon - z r - π/4\\bigr)}
_{\\text{steady harmonic}}
\;+\;
\\underbrace{T_{\\text{trans}}(z,t)}
_{\\text{transient integral}} ,
where :math:`r = \\sqrt{\\omega / 2\\kappa}`.
The first term is the *steady* periodic component that
propagates downward with exponentially damped amplitude and a
depth-dependent phase lag. The second term accounts for the
*transient* adjustment from the initial uniform temperature *Tᵢ* and
is evaluated numerically here by vectorised trapezoidal quadrature.
Parameters
----------
z : float
Depth below the soil surface (m, positive downward).
t : float or array_like
Time since the start of the forcing (s). Scalar or NumPy array.
A : float
Amplitude of the sinusoidal **surface heat flux** (W m⁻², positive
downward). Consistent with
:func:`sinusoidal_boundary_flux`.
omega : float
Angular frequency of the forcing (rad s⁻¹).
For a diurnal wave ``omega = 2 * π / 86_400``.
epsilon : float
Phase shift ε (rad) of the surface heat-flux wave.
Ti : float
Initial uniform soil temperature *Tᵢ* (°C or K).
kappa : float
Thermal diffusivity κ of the soil (m² s⁻¹).
Returns
-------
float or ndarray
Soil temperature *T(z, t)* in the same units as *Ti*.
If *t* is an array the returned array has the same shape.
Notes
-----
* **Steady component** –
The exponential term ``exp(-z*r)`` dampens amplitude with depth
while the phase lag is ``z*r + π/4``.
* **Transient component** –
The integral is truncated at ``x = 50 / z`` (or 50 if *z = 0*) and
evaluated with 2 000 panels, which empirically yields < 0.1 %
relative error for typical soil parameters. Adjust the limits or
panel count for higher precision.
* **Vectorisation** –
For array *t* the quadrature is performed in parallel using
broadcasting; memory usage scales with ``len(t)`` × 2 000.
* **Units** –
Ensure *A* is W m⁻² **heat flux**. If you have a surface
temperature wave instead, transform to an equivalent heat-flux
boundary or modify the formulation.
Examples
--------
>>> import numpy as np, matplotlib.pyplot as plt
>>> z = 0.05 # 5 cm
>>> k = 1.4e-7 # m² s⁻¹
>>> A = 120 # W m⁻²
>>> ω = 2*np.pi/86400 # diurnal
>>> ε = 0 # no phase shift
>>> Ti = 15.0 # °C
>>> t = np.linspace(0, 172800, 2000) # 2 days
>>> Tz = soil_temperature_sinusoidal(z, t, A, ω, ε, Ti, k)
>>> plt.plot(t/3600, Tz)
>>> plt.xlabel("Time (h)")
>>> plt.ylabel("Temperature (°C)")
>>> plt.title("Analytical diurnal soil temperature at 5 cm")
>>> plt.show()
"""
r = np.sqrt(omega / (2.0 * kappa))
exp_term = np.exp(-z * r)
phase = omega * t + epsilon - z * r - math.pi / 4.0
steady = A / (kappa * np.sqrt(omega)) * exp_term * np.sin(phase)
# Transient integral (third term) – numeric quadrature (vectorised)
def _integrand(x):
return (
(kappa * x**2 * math.sin(epsilon) - omega * math.cos(epsilon))
* np.cos(x * z)
/ (omega**2 + (kappa**2) * x**4)
)
if np.isscalar(t):
# Scalar: quad via np.trapz on a finite domain
xi = np.linspace(0.0, 50.0 / z if z else 50.0, 2000)
transient = (
-2
* A
* kappa
/ math.pi
* _trapz(_integrand(xi) * np.exp(-kappa * xi**2 * t), xi)
)
else:
transient = np.zeros_like(t, dtype=float)
xi = np.linspace(0.0, 50.0 / z if z else 50.0, 2000)
integ = _integrand(xi)[:, None] * np.exp(-kappa * xi[:, None] ** 2 * t[None, :])
transient = -2 * A * kappa / math.pi * _trapz(integ, xi, axis=0)
return Ti + steady + transient
[docs]
def soil_heat_flux_sinusoidal(
z: float,
t: float | np.ndarray,
A: float,
omega: float,
epsilon: float,
kappa: float,
) -> float | np.ndarray:
"""
Analytical **soil heat–flux** response *G(z,t)* to a sinusoidal
surface–flux boundary condition.
A semi-infinite, homogeneous soil column forced at the surface by
.. math::
q_0(t) \;=\; A \,\sin\!\bigl(\omega t + \varepsilon\bigr)
admits the Green-function solution (Carslaw & Jaeger, 1959, Eq. 15)
.. math::
G(z,t)
\;=\;
A\,e^{-z r}\,\sin(\omega t + \varepsilon - z r)
\;+\;
G_{\text{trans}}(z,t),
where :math:`r = \sqrt{ \omega / 2\kappa }` and the transient term
.. math::
G_{\text{trans}}(z,t)
\;=\;
-\frac{2 A \kappa}{\pi}
\int_{0}^{\infty}
\\frac{( \kappa x^{2}\sin\varepsilon - \omega \cos\varepsilon )
\;x \sin(x z)}
{\omega^{2} + \kappa^{2} x^{4}}
\,e^{-\kappa x^{2} t}\; \mathrm d x .
The first term is the steady, exponentially damped harmonic that
propagates downward with a depth-dependent phase lag; the second
term describes the transient adjustment from an initially unheated
half-space. The integral is evaluated here by vectorised
trapezoidal quadrature on a finite domain (*x* ≤ 50 / *z*).
Parameters
----------
z : float
Depth below the soil surface (m, positive downward).
t : float or array_like
Time since the onset of forcing (s). Scalar or NumPy array.
A : float
Amplitude of the *surface* heat flux (W m⁻², positive downward).
omega : float
Angular frequency ω of the forcing (rad s⁻¹).
For a diurnal wave ``omega = 2 π / 86 400``.
epsilon : float
Phase shift ε of the surface flux (radians).
kappa : float
Thermal diffusivity κ of the soil (m² s⁻¹).
Returns
-------
float or ndarray
Heat flux *G(z,t)* at depth *z* (W m⁻²) with shape equal to
``t`` after NumPy broadcasting.
Positive values **downward**, matching the sign convention of *A*.
Notes
-----
* **Steady component** –
Amplitude decays as ``exp(-z*r)``; phase lag increases linearly
with depth by *z r*.
* **Transient component** –
The integration limit ``50 / z`` (or 50 for *z = 0*) yields < 0.1 %
relative error for common parameter ranges. Increase the upper
bound and/or panel count (2 000) for stricter accuracy.
* **Vectorisation** –
For array *t* the quadrature is evaluated for all time steps
simultaneously via broadcasting; memory usage scales with
``len(t) × 2000``.
* **Coupling with temperature** –
The companion function
:func:`soil_temperature_sinusoidal` gives *T(z,t)* under the
same boundary forcing; *G* and *T* satisfy Fourier’s law
``G = -k ∂T/∂z`` once thermal conductivity *k* is specified.
References
----------
Carslaw, H. S., & Jaeger, J. C. (1959).
*Conduction of Heat in Solids* (2nd ed., pp. 100–102).
Oxford University Press.
Examples
--------
>>> import numpy as np, matplotlib.pyplot as plt
>>> z = 0.05 # 5 cm
>>> kappa = 1.4e-7 # m² s⁻¹
>>> A = 120 # W m⁻² downward
>>> omega = 2*np.pi/86400 # diurnal frequency
>>> eps = 0
>>> t = np.linspace(0, 172800, 2000) # 2 days
>>> G = soil_heat_flux_sinusoidal(z, t, A, omega, eps, kappa)
>>> plt.plot(t/3600, G)
>>> plt.xlabel("Time (h)")
>>> plt.ylabel("Heat flux $G$ (W m$^{-2}$)")
>>> plt.title("Analytical diurnal soil heat flux at 5 cm")
>>> plt.show()
"""
r = np.sqrt(omega / (2.0 * kappa))
exp_term = np.exp(-z * r)
phase = omega * t + epsilon - z * r
steady = A * exp_term * np.sin(phase)
# Transient integral similar to temperature – numeric quadrature
def _integrand(x):
return (
(kappa * x**2 * math.sin(epsilon) - omega * math.cos(epsilon))
* x
* np.sin(x * z)
/ (omega**2 + (kappa**2) * x**4)
)
if np.isscalar(t):
xi = np.linspace(0.0, 50.0 / z if z else 50.0, 2000)
transient = (
-2
* A
* kappa
/ math.pi
* _trapz(_integrand(xi) * np.exp(-kappa * xi**2 * t), xi)
)
else:
xi = np.linspace(0.0, 50.0 / z if z else 50.0, 2000)
integ = _integrand(xi)[:, None] * np.exp(-kappa * xi[:, None] ** 2 * t[None, :])
transient = -2 * A * kappa / math.pi * _trapz(integ, xi, axis=0)
return steady + transient
# -----------------------------------------------------------------------------
# 5. Soil-property parameterisations (Eqs. 16–18)
# -----------------------------------------------------------------------------
[docs]
def heat_capacity_moist_soil(
theta_v: float | np.ndarray,
theta_s: float,
rho_c_s: float = 1.26e6,
rho_c_w: float = 4.20e6,
) -> float | np.ndarray:
"""Volumetric heat capacity of moist soil, Eq. (16).
Parameters
----------
theta_v
Volumetric water content **(m³ m-3)**.
theta_s
Porosity (saturated volumetric water content) **(m³ m-3)**.
rho_c_s, rho_c_w
Heat capacity of dry soil / water **(J m-3 K-1)**.
"""
return theta_v * rho_c_w + (1.0 - theta_s) * rho_c_s
[docs]
def pf_from_theta(
theta_v: float | np.ndarray, theta_s: float, psi_s: float, b: float
) -> float | np.ndarray:
"""Return Pf (Eq. 18) from volumetric water content."""
return np.log10(100.0 * psi_s * (theta_s / theta_v) ** b)
[docs]
def thermal_conductivity_moist_soil(
theta_v: float | np.ndarray, theta_s: float, psi_s: float, b: float
) -> float | np.ndarray:
"""Thermal conductivity parameterisation, Eq. (17)."""
Pf = pf_from_theta(theta_v, theta_s, psi_s, b)
k = np.where(Pf <= 5.1, 0.420 * np.exp(-Pf - 2.7), 0.1744)
return k
[docs]
def thermal_diffusivity(
k: float | np.ndarray, rho_c: float | np.ndarray
) -> float | np.ndarray:
"""Return κ = k / (ρ c)."""
return k / rho_c