soil_heat.gao_et_al
Functions
|
Compute the soil thermal conductivity \(\lambda_s\) |
|
Calculate the soil thermal diffusivity \(k_s\) |
|
Compute the volumetric heat capacity \(C_v\) of soil: |
|
Calculate the normalized mean error (NME) between calculated |
|
Compute the root-mean-square error (RMSE) between calculated |
|
Estimate the soil heat flux at a target depth z |
|
Estimate soil heat flux \(G(z)\) at a shallow depth |
|
Estimate soil heat flux \(G(z,t)\) at depth z based on a |
|
Compute soil heat flux \(G(z,t)\) from the **H04 harmonic |
|
Compute soil heat flux at the end of a temperature record |
|
Estimate the thermal damping depth d (m) from the exponential |
|
Extrapolate soil heat flux from a reference depth z_r to a |
|
Ground-heat-flux estimate at a target depth z using the |
|
Estimate soil heat flux at a shallow depth z from a flux-plate |
|
WBZ12-S analytic–numeric solution for soil heat flux G(z, t) at |
|
Compute the exact analytical soil‐temperature profile |
|
Exact analytical soil-heat-flux solution for sinusoidal surface forcing. |
Module Contents
- soil_heat.gao_et_al.lambda_s(theta: numpy.ndarray | float) numpy.ndarray | float[source]
Compute the soil thermal conductivity \(\lambda_s\) as a function of volumetric water content (θ).
The relationship is taken from Gao et al. (2017, Eq. 12):
\[\lambda_s(\theta) = 0.20 + \exp\bigl[\,1.46\,(\theta - 0.34)\bigr]\]- Parameters:
theta (float or array_like) – Volumetric water content (m³ m⁻³). Accepts a scalar value or any array-like object that can be converted to a
numpy.ndarray. Values should lie in the closed interval[0, 1].- Returns:
Soil thermal conductivity λ_s in W m⁻¹ K⁻¹. The returned type matches the input: a scalar is returned for a scalar θ, and a NumPy array for array-like θ.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element of θ is outside the physically meaningful range
[0, 1].
Notes
Vectorization – The function is fully vectorized; it operates element-wise on NumPy arrays and broadcasts according to NumPy broadcasting rules.
Empirical limits – Gao et al. recommend the equation for mineral soils where 0 ≤ θ ≤ 0.5. Extrapolation beyond this range may introduce error.
Units – The output uses SI units (W m⁻¹ K⁻¹).
References
Gao, Z., Niu, G.-Y., & Hedquist, B. C. (2017). Soil thermal conductivity parameterization: A closed-form equation and its evaluation. Journal of Geophysical Research: Atmospheres, 122(6), 3466–3478. DOI: 10.1002/2016JD025992
Examples
>>> lambda_s(0.25) 0.463... >>> import numpy as np >>> theta = np.linspace(0, 0.5, 5) >>> lambda_s(theta) array([0.20 , 0.31243063, 0.41172297, 0.51210322, 0.61861198])
- soil_heat.gao_et_al.k_s(theta: numpy.ndarray | float) numpy.ndarray | float[source]
Calculate the soil thermal diffusivity \(k_s\) as a function of volumetric water content (θ).
The empirical relationship follows Gao et al. (2017, Eq. 13):
\[k_s(\theta) = \bigl[0.69 + \exp\bigl(3.06\,(\theta - 0.26)\bigr)\bigr] \times 10^{-7}\]- Parameters:
theta (float or array_like) – Volumetric water content (m³ m⁻³). Accepts a scalar or any array-like sequence convertible to a
numpy.ndarray. Values must lie in the closed interval[0, 1].- Returns:
Soil thermal diffusivity k_s in m² s⁻¹. A scalar is returned if θ is a scalar; otherwise a NumPy array of matching shape is returned.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element of θ is outside the physically meaningful range
[0, 1].
Notes
Vectorization – The calculation is fully vectorized and broadcasts according to NumPy rules.
Applicability – Gao et al. derived this expression for mineral soils under typical field conditions. Extrapolating beyond 0 ≤ θ ≤ 0.5 may reduce accuracy.
Units – Output is in SI units (m² s⁻¹).
References
Gao, Z., Niu, G.-Y., & Hedquist, B. C. (2017). Soil thermal conductivity parameterization: A closed-form equation and its evaluation. Journal of Geophysical Research: Atmospheres, **122**(6), 3466–3478. https://doi.org/10.1002/2016JD025992
Examples
>>> k_s(0.25) 1.321...e-07 >>> thetas = np.linspace(0, 0.5, 6) >>> k_s(thetas) array([1.14...e-07, 1.21...e-07, 1.30...e-07, 1.41...e-07, 1.54...e-07, 1.69...e-07])
- soil_heat.gao_et_al.volumetric_heat_capacity(lambda_s_val: numpy.ndarray | float, k_s_val: numpy.ndarray | float) numpy.ndarray | float[source]
Compute the volumetric heat capacity \(C_v\) of soil:
\[C_v \;=\; \frac{\lambda_s}{k_s}\]where \(\lambda_s\) is the thermal conductivity (W m⁻¹ K⁻¹) and \(k_s\) is the thermal diffusivity (m² s⁻¹).
The resulting heat capacity has units of J m⁻³ K⁻¹.
- Parameters:
- Returns:
Volumetric heat capacity C_v (J m⁻³ K⁻¹). The return type matches the input: a scalar for scalar inputs, or a NumPy array for array-like inputs.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element in k_s_val is zero or negative, which would lead to division by zero or non-physical results.
Notes
Vectorization – The function is fully vectorized; both inputs are converted to
numpy.ndarrayand follow NumPy broadcasting rules.Physical meaning – Volumetric heat capacity represents the energy required to raise the temperature of a unit volume of soil by one kelvin. High values correspond to moist or water-logged soils; dry mineral soils have lower C_v.
Examples
>>> #Scalar inputs >>> volumetric_heat_capacity(0.8, 1.2e-6) 666666.666...
>>> #Vectorized inputs >>> lambda_vals = np.array([0.5, 0.6, 0.7]) >>> diffusivities = np.array([1.1e-6, 1.2e-6, 1.3e-6]) >>> volumetric_heat_capacity(lambda_vals, diffusivities) array([454545.4545..., 500000. ..., 538461.5384...])
- soil_heat.gao_et_al.nme(calc: numpy.ndarray | float, meas: numpy.ndarray | float) float[source]
Calculate the normalized mean error (NME) between calculated and measured values, expressed as a percentage.
The formulation follows Gao et al. (2017, Eq. 14):
\[\text{NME} \;=\; 100\,\frac{\sum\limits_{i}\left|\hat{y}_i - y_i\right|} {\sum\limits_{i}\left|y_i\right|}\]where \(\hat{y}_i\) are the calculated (model) values and \(y_i\) are the measured (reference) values.
- Parameters:
- Returns:
Normalized mean error (percent). A value of 0 % indicates perfect agreement; larger values indicate greater error.
- Return type:
- Raises:
If calc and meas cannot be broadcast to a common shape.
If the denominator
Σ|meas|equals zero (e.g., all measured values are zero), which would make NME undefined.
Notes
Range – NME is non-negative and unbounded above.
Interpretation – Because both numerator and denominator use absolute values, NME is insensitive to the direction of the error (over- vs under-prediction) and therefore complements signed error metrics such as mean bias.
Vectorization – The implementation is fully vectorized and adheres to NumPy broadcasting rules.
References
Gao, Z., Niu, G.-Y., & Hedquist, B. C. (2017). Soil thermal conductivity parameterization: A closed-form equation and its evaluation. Journal of Geophysical Research: Atmospheres, 122(6), 3466–3478. https://doi.org/10.1002/2016JD025992
Examples
>>> # Single values >>> nme(4.5, 5.0) 10.0
>>> # Vectors >>> calc = np.array([1.0, 2.1, 3.2]) >>> meas = np.array([1.2, 2.0, 3.0]) >>> nme(calc, meas) 4.7619...
>>> # Broadcasting (scalar vs array) >>> nme(2.0, np.array([1.5, 2.5, 2.0])) 16.6666...
- soil_heat.gao_et_al.rmse(calc: numpy.ndarray | float, meas: numpy.ndarray | float) float[source]
Compute the root-mean-square error (RMSE) between calculated (model) and measured (reference) values.
Following Gao et al. (2017, Eq. 15):
\[\text{RMSE} \;=\; \sqrt{\frac{1}{N}\sum_{i=1}^{N}\bigl(\hat{y}_i - y_i\bigr)^2}\]where \(\hat{y}_i\) are calculated values and \(y_i\) are measured values.
- Parameters:
- Returns:
Root-mean-square error (same units as the inputs).
- Return type:
- Raises:
If calc and meas cannot be broadcast to a common shape.
If the input arrays are empty (
N = 0).
Notes
Interpretation – RMSE represents the sample-standard-deviation of the differences between two datasets; lower values indicate better agreement.
Vectorization – The function is fully vectorized and respects NumPy broadcasting rules.
Units – RMSE preserves the units of the input variables.
References
Gao, Z., Niu, G.-Y., & Hedquist, B. C. (2017). Soil thermal conductivity parameterization: A closed-form equation and its evaluation. Journal of Geophysical Research: Atmospheres, **122**(6), 3466–3478. https://doi.org/10.1002/2016JD025992
Examples
>>> # Scalar inputs >>> rmse(4.5, 5.0) 0.5
>>> # Vector inputs >>> calc = np.array([2.1, 3.0, 4.2]) >>> meas = np.array([2.0, 3.5, 4.0]) >>> rmse(calc, meas) 0.2645...
>>> # Broadcasting >>> rmse(3.0, np.array([2.5, 3.5, 3.0])) 0.4082...
- soil_heat.gao_et_al.calorimetric_gz(g_zr: numpy.ndarray | float, cv_layers: numpy.ndarray | float, dT_dt_layers: numpy.ndarray | float, dz_layers: numpy.ndarray | float) numpy.ndarray | float[source]
Estimate the soil heat flux at a target depth z (typically 5 cm) using the calorimetric method.
The method corrects an in-situ heat-flux plate reading made at a reference depth z_r by adding the change in heat storage of the soil column located between z and z_r. Mathematically (Liebethal & Foken 2007, Eq. 1):
\[G(z,t) \;=\; G(z_r,t) \;+\; \sum_{l=1}^{N} C_{v,l}\,\frac{\partial \bar{T}_l}{\partial t}\,\Delta z_l\]where
\(G(z,t)\) … heat flux at depth z (W m⁻²)
\(G(z_r,t)\) … measured plate flux at reference depth z_r
\(C_{v,l}\) … volumetric heat capacity of layer l (J m⁻³ K⁻¹)
\(\partial \bar{T}_l / \partial t\) … time derivative of the layer-averaged temperature (K s⁻¹)
\(\Delta z_l\) … thickness of layer l (m)
- Parameters:
g_zr (float or array_like) – Heat-flux plate measurement at depth z_r (W m⁻²). Can be a scalar or time series.
cv_layers (array_like) – Volumetric heat capacity \(C_{v,l}\) for each of N soil sub-layers (J m⁻³ K⁻¹). Shape
(N, …)where the trailing dimensions (…) must be broadcast-compatible with g_zr.dT_dt_layers (array_like) – Time derivative of layer-mean temperature \(\partial \bar{T}_l / \partial t\) (K s⁻¹); same shape as cv_layers.
dz_layers (array_like) – Thickness of each layer \(\Delta z_l\) (m); shape
(N,)or broadcast-compatible with the first axis of cv_layers.
- Returns:
Calorimetrically corrected soil heat flux G(z) (W m⁻²). Scalar if all inputs are scalar; otherwise a NumPy array matching the broadcast shape of g_zr.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If the inputs cannot be broadcast to a common shape, or if the number of layers inferred from cv_layers, dT_dt_layers, and dz_layers are inconsistent.
Notes
Vectorization – All operations are fully vectorized. Inputs are converted to
numpy.ndarrayand follow NumPy broadcasting rules.Layer axis – The first dimension of cv_layers and dT_dt_layers (axis 0) is interpreted as the layer index l. Layer thicknesses dz_layers are broadcast across any additional dimensions.
Units – Ensure consistent SI units: W m⁻², J m⁻³ K⁻¹, K s⁻¹, and m.
References
Liebethal, C., & Foken, T. (2007). Evaluation of six parameterization approaches for the ground heat flux. Agricultural and Forest Meteorology, **143**(1–2), 65-80. https://doi.org/10.1016/j.agrformet.2006.11.001
Examples
>>> # Three-layer example, single time step >>> g_plate = -15.2 # W m-2 at z_r = −0.05 m >>> Cv = np.array([2.5e6, 2.3e6, 2.1e6]) # J m-3 K-1 >>> dTdt = np.array([1.2e-4, 0.9e-4, 0.6e-4]) # K s-1 >>> dz = np.array([0.02, 0.02, 0.01]) # m >>> calorimetric_gz(g_plate, Cv, dTdt, dz) -4.06...
>>> # Vectorized daily time series with two layers >>> g_plate = np.random.normal(-10, 2, 1440) # per minute >>> Cv = np.array([[2.4e6], [2.2e6]]) # (2,1) >>> dTdt = np.random.normal(5e-5, 2e-5, (2,1440)) >>> dz = np.array([0.03, 0.02]) >>> Gz = calorimetric_gz(g_plate, Cv, dTdt, dz) # shape (1440,)
- soil_heat.gao_et_al.force_restore_gz(cv: numpy.ndarray | float, dTg_dt: numpy.ndarray | float, Tg: numpy.ndarray | float, Tg_bar: numpy.ndarray | float, delta_z: float = 0.05, omega: float = OMEGA_DAY) numpy.ndarray | float[source]
Estimate soil heat flux \(G(z)\) at a shallow depth (\(z = \delta z\), default 5 cm) with the force–restore method.
The formulation (Liebethal & Foken 2007, Eq. 2) corrects the calorimetric storage term with a restore component that accounts for the difference between the instantaneous ground temperature T_g and its running mean bar{T}_g:
\[G(z,t) \;=\; C_v \, \delta z \, \frac{\partial T_g}{\partial t} \;+ \sqrt{\omega \, C_v \, \lambda_s(C_v)} \left[ \frac{1}{\omega}\,\frac{\partial T_g}{\partial t} + \bigl(T_g - \bar{T}_g\bigr) \right]\]where
\(C_v\) … volumetric heat capacity (J m⁻³ K⁻¹)
\(\partial T_g / \partial t\) … ground-temperature tendency (K s⁻¹)
\(\lambda_s(C_v)\) … soil thermal conductivity derived from C_v via
lambda_s_from_cv()(W m⁻¹ K⁻¹)\(\omega\) … angular frequency of the diurnal cycle (rad s⁻¹)
\(T_g\) / \(\bar{T}_g\) … instantaneous and running-mean ground temperature (K)
\(\delta z\) … sensor depth below the surface (m)
- Parameters:
cv (float or array_like) – Volumetric heat capacity \(C_v\) (J m⁻³ K⁻¹). Must be positive. Accepts scalars or NumPy-broadcastable arrays.
dTg_dt (float or array_like) – Time derivative \(\partial T_g/\partial t\) (K s⁻¹).
Tg (float or array_like) – Instantaneous ground (surface) temperature \(T_g\) (K or °C) at depth δz.
Tg_bar (float or array_like) – Running mean ground temperature \(\bar{T}_g\) over the diurnal cycle (same units as Tg).
delta_z (float, optional) – Depth \(\delta z\) in metres; default is 0.05 m (5 cm).
omega (float, optional) – Angular frequency \(\omega\) in rad s⁻¹. Defaults to :pydata:`OMEGA_DAY` (2π / 86 400 s).
- Returns:
Soil heat flux G(z) at depth δz (W m⁻²). The output shape follows NumPy broadcasting rules applied to the inputs.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any input arrays cannot be broadcast to a common shape, or if cv contains non-positive values.
Notes
λ_s(C_v) mapping – The function relies on a helper
lambda_s_from_cv()that converts volumetric heat capacity to thermal conductivity. Ensure that this helper is present in the import path.Units – Keep units internally consistent (SI).
Vectorization – All operations are vectorized; the mathematical expression is evaluated element-wise for array inputs.
Interpretation – The first term represents storage (calorimetric), while the second tempers short-term fluctuations, “restoring” T_g toward its mean.
References
Liebethal, C., & Foken, T. (2007). Evaluation of six parameterization approaches for the ground heat flux. Agricultural and Forest Meteorology, 143(1–2), 65-80. https://doi.org/10.1016/j.agrformet.2006.11.001
Examples
>>> Cv = 2.4e6 # J m-3 K-1 >>> dTgdt = 1.0e-4 # K s-1 >>> Tg = 293.5 # K >>> Tg_bar = 291.7 # K (running mean) >>> force_restore_gz(Cv, dTgdt, Tg, Tg_bar) 19.3... # W m-2
>>> #Vectorized daily record >>> cv_arr = np.full(1440, 2.2e6) >>> dT_arr = np.gradient(np.sin(np.linspace(0, 2*np.pi, 1440))) / 60 >>> Gz_ts = force_restore_gz(cv_arr, dT_arr, 298+2*np.sin(...), ... 298*np.ones_like(dT_arr))
- soil_heat.gao_et_al.gao2010_gz(AT: numpy.ndarray | float, lambda_s_val: numpy.ndarray | float, k_s_val: numpy.ndarray | float, t: numpy.ndarray | float, omega: float = OMEGA_DAY) numpy.ndarray | float[source]
Estimate soil heat flux \(G(z,t)\) at depth z based on a sinusoidal ground-temperature forcing (Gao et al. 2010, Eq. 3).
The analytical solution assumes that the surface (or ground-contact) temperature varies sinusoidally with amplitude A_T and angular frequency ω. Under these conditions the heat flux at any depth z can be written
\[G(z,t) \;=\; \sqrt{2}\, \frac{\lambda_s A_T}{d}\, \sin\bigl(\omega t + \tfrac{\pi}{4}\bigr),\]where the thermal damping depth
\[d \;=\; \sqrt{\frac{2 k_s}{\omega}}\]is determined by the soil thermal diffusivity k_s and the forcing frequency ω.
- Parameters:
AT (float or array_like) – Amplitude of the sinusoidal ground-surface temperature (K).
lambda_s_val (float or array_like) – Soil thermal conductivity \(\lambda_s\) (W m⁻¹ K⁻¹).
k_s_val (float or array_like) – Soil thermal diffusivity \(k_s\) (m² s⁻¹).
t (float or array_like) – Time variable (s). Can be absolute time since epoch or simply seconds since the start of the cycle—as long as ω t is dimensionless.
omega (float, optional) – Angular frequency ω (rad s⁻¹). Defaults to OMEGA_DAY (≈ 7.272 × 10⁻⁵ s⁻¹, i.e. 2π / 86 400 s).
- Returns:
Heat flux G(z,t) (W m⁻²). The output follows NumPy’s broadcasting rules applied to the inputs.
- Return type:
float or numpy.ndarray
- Raises:
If lambda_s_val, k_s_val, and t cannot be broadcast to a common shape.
If any element of k_s_val or omega is non-positive.
Notes
Vectorization – All inputs are internally converted to
numpy.ndarray; the formula is evaluated element-wise and fully supports broadcasting.Units – Ensure consistent SI units: W m⁻¹ K⁻¹ (λ_s), m² s⁻¹ (k_s), s (t), and rad s⁻¹ (ω).
Scope – The solution presumes purely sinusoidal boundary forcing and homogeneous soil properties; real-world deviations (e.g., non-sine forcing, stratified soils, moisture variation) will introduce error.
References
Gao, Z., Horton, R., Luo, L., & Kucharik, C. J. (2010). A simple method to measure soil temperature dynamics: Theory and application. Soil Science Society of America Journal, 74(2), 580-588. https://doi.org/10.2136/sssaj2009.0169
Examples
>>> # Daily cycle at 5 cm depth >>> AT = 8.0 # K >>> lambda_ = 1.2 # W m-1 K-1 >>> kappa = 1.0e-6 # m2 s-1 >>> t_day = np.linspace(0, 86400, 97) # 15-min steps >>> Gz = gao2010_gz(AT, lambda_, kappa, t_day) >>> Gz.shape (97,)
- soil_heat.gao_et_al.heusinkveld_gz(A_n: numpy.ndarray | float, Phi_n: numpy.ndarray | float, n_max: int, k_s_val: numpy.ndarray | float, lambda_s_val: numpy.ndarray | float, w: float) numpy.ndarray | float[source]
Compute soil heat flux \(G(z,t)\) from the H04 harmonic series solution proposed by Heusinkveld et al. (2004, Eq. 4).
The approach represents the surface (ground) temperature as a Fourier series with harmonics up to order n_max. For each harmonic \(n\) the heat-flux contribution at depth z can be written
\[G_n(z,t) \;=\; \frac{\lambda_s}{10\,\pi}\, A_n\; \sqrt{\frac{1}{k_s\,n\,\omega\,k_s}}\; \sin\bigl(n\,\omega\,t + \Phi_n + \tfrac{\pi}{4}\bigr),\]and the full signal is obtained by summing over n = 1…n_max.
Note The implementation below follows the algebraic form that appeared in H04; consult the original paper for derivation details and recommended parameter ranges.
- Parameters:
A_n (float or array_like) – Amplitudes \(A_n\) of the n-th harmonic of surface-temperature forcing (K). Provide either * a single scalar applied to every harmonic, or * an array of length ≥ n_max giving amplitude for each harmonic.
Phi_n (float or array_like) – Phase shifts \(\Phi_n\) (rad) corresponding to each harmonic order. Same broadcasting rules as A_n.
n_max (int) – Highest harmonic order to include in the summation.
k_s_val (float or array_like) – Soil thermal diffusivity \(k_s\) (m² s⁻¹).
lambda_s_val (float or array_like) – Soil thermal conductivity \(\lambda_s\) (W m⁻¹ K⁻¹).
w (float) – Fundamental angular frequency \(\omega\) (rad s⁻¹), e.g. \(2\pi/86400\) for a 24-h cycle.
- Returns:
Soil heat flux G(z,t) (W m⁻²). The return shape is the broadcast shape of the input arrays (excluding the harmonic axis).
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If n_max is less than 1, or if k_s_val or w are non-positive, or if the amplitudes/phases cannot be broadcast to length n_max.
Notes
Vectorization – Internally, the harmonic index axis has length n_max (
n = np.arange(1, n_max+1)). All other dimensions come from broadcasting A_n, Phi_n, k_s_val, lambda_s_val, and t (if vectorized in the caller).Units – Consistency with SI units is assumed.
Interpretation – The prefactor
λ_s / (10 π)appears in H04’s original derivation. If you adopt a different convention (e.g., depth-explicit damping), modify accordingly.
References
Heusinkveld, B. G., Jacobs, A. F. G., Holtslag, A. A. M., & et al. (2004). Surface energy balance closure in an arid region: The role of soil heat flux. Agricultural and Forest Meteorology, **122**(1), 21-37. https://doi.org/10.1016/j.agrformet.2003.09.005
Examples
>>> # Daily cycle with three harmonics >>> A = np.array([10, 4, 1.5]) # K >>> Phi = np.array([0, -np.pi/6, np.pi/8]) >>> Gz = heusinkveld_gz(A, Phi, n_max=3, ... k_s_val=1e-6, ... lambda_s_val=1.0, ... w=2*np.pi/86400) >>> Gz.shape ()
- soil_heat.gao_et_al.hsieh2009_gz(tz_series: numpy.ndarray | list | tuple, time_series: numpy.ndarray | list | tuple, cv_series: numpy.ndarray | list | tuple, ks_series: numpy.ndarray | list | tuple) float[source]
Compute soil heat flux at the end of a temperature record using the half-order integral (Hsieh et al., 2009, Eq. 5).
The method exploits the analytical solution of the one-dimensional heat-conduction equation for a semi-infinite medium, leading to a convolution integral of order ½ that relates the time series of near-surface soil temperature to the downward heat flux:
\[G(t_N) \;=\; 2\sqrt{\frac{k_s(t_N)\,C_v(t_N)}{\pi}}\; \int_{t_0}^{t_N} \frac{\partial T(z,t')}{\partial t'} \left(t_N - t'\right)^{-1/2}\,dt'\]In discrete form with linear interpolation between measurements \(t_i \;(i = 0,\dots,N)\),
\[\int_{t_0}^{t_N}\! \frac{dT}{dt'}\,(t_N-t')^{-1/2}dt' \approx \sum_{i=0}^{N-1} \frac{\Delta T_i}{\Delta t_i}\! \left[(t_N-t_i)^{1/2} - (t_N-t_{i+1})^{1/2}\right],\]where \(\Delta T_i = T_{i+1}-T_i\) and \(\Delta t_i = t_{i+1}-t_i\).
The final scalar result corresponds to the flux at \(t_N\;(=\text{time_series}[-1])\).
- Parameters:
tz_series (array_like) – Near-surface (or shallow-depth) soil temperature observations T(z,t) (K). Length N ≥ 2.
time_series (array_like) – Strictly monotonically increasing time stamps (s) matching tz_series. Same length N.
cv_series (array_like) – Volumetric heat capacity C_v (J m⁻³ K⁻¹) at each time stamp. Same length N.
ks_series (array_like) – Thermal diffusivity k_s (m² s⁻¹) at each time stamp. Same length N.
- Returns:
Soil heat flux G(t_N) at the final time point (W m⁻²).
- Return type:
- Raises:
If the four series differ in length or have fewer than two samples.
If time_series is not strictly increasing.
If any element of cv_series or ks_series is non-positive.
Notes
The algorithm uses a forward finite-difference for \(dT/dt'\) and trapezoidal integration over each interval [t_i, t_{i+1}].
Only the latest values of k_s and C_v are used, consistent with the Hsieh et al. derivation. Replace with a time-variable kernel if property changes are large over the record.
All inputs are cast to
numpy.ndarraywithdtype=float.
References
Hsieh, C.-I., Katul, G., & Chi, T.-C. (2009). Retrieval of soil heat flux from soil temperature data by the continuous-time heat equation model. Water Resources Research, 45, W08433. https://doi.org/10.1029/2009WR007891
Examples
>>> times = np.array([0, 600, 1200, 1800]) # every 10 min >>> temps = np.array([292.5, 293.0, 293.6, 294.0]) # K >>> Cv = np.full_like(temps, 2.3e6) # J m-3 K-1 >>> ks = np.full_like(temps, 1.1e-6) # m2 s-1 >>> hsieh2009_gz(temps, times, Cv, ks) -2.13...
- soil_heat.gao_et_al.leuning_damping_depth(z: numpy.ndarray | float, zr: numpy.ndarray | float, AT_z: numpy.ndarray | float, AT_zr: numpy.ndarray | float) numpy.ndarray | float[source]
Estimate the thermal damping depth d (m) from the exponential decay of temperature–wave amplitude with depth.
The formulation is based on Eq. (6) of Leuning et al. (1985) for a one-dimensional, sinusoidally forced soil column:
\[d \;=\; \frac{z_r - z}{\ln\bigl(A_T(z) / A_T(z_r)\bigr)}\]where
\(A_T(z)\) — amplitude of the temperature wave at depth z
\(z_r\) — reference depth (m)
\(A_T(z_r)\) — amplitude at reference depth
The equation assumes amplitudes decrease monotonically with depth following \(A_T(z) = A_T(0)\,\exp(-z/d)\).
- Parameters:
z (float or array_like) – Depth(s) (m) at which the amplitude A_T(z) was measured. Positive values denote depth below the surface.
zr (float or array_like) – Reference depth(s) z_r (m). Must be broadcast-compatible with z. If zr < z the denominator switches sign, yielding a negative d that should be interpreted carefully.
AT_z (float or array_like) – Temperature-wave amplitude at z (K).
AT_zr (float or array_like) – Temperature-wave amplitude at z_r (K).
- Returns:
Thermal damping depth d (m). The output shape follows NumPy broadcasting rules applied to the four inputs.
- Return type:
float or numpy.ndarray
- Raises:
If any of the following conditions occur:
AT_z or AT_zr contain non-positive values (log undefined).
AT_z and AT_zr are equal everywhere, leading to
ln(1) = 0and division by zero.The inputs cannot be broadcast to a common shape.
Notes
Vectorisation — Internally, all inputs are converted to
ndarray; the formula is applied element-wise and fully supports broadcasting.Physical meaning — Larger d indicates slower attenuation of temperature fluctuations with depth (i.e. higher thermal diffusivity or conductivity). Very small logarithmic denominators (AT_z ≈ AT_zr) imply an extremely deep damping depth that may fall outside the soil column considered.
Units — Depths must share the same units (metres). Amplitudes may be in kelvins or degrees Celsius provided they use identical scaling.
References
Leuning, R., Barlow, E. W. R., & Paltridge, G. (1985). Temperature gradients in a soil: Theory and measurement. Agricultural and Forest Meteorology, 35(1–4), 127–143. https://doi.org/10.1016/0168-1923(85)90067-3
Examples
>>> d = leuning_damping_depth(z=0.10, zr=0.05, AT_z=4.0, AT_zr=8.0) >>> round(d, 4) 0.0721
>>> #Broadcasting with arrays >>> depths = np.array([0.05, 0.10, 0.20]) >>> amps = np.array([8.0, 4.0, 1.5]) >>> d_arr = leuning_damping_depth(depths, 0.02, amps, 10.0) >>> d_arr array([0.0380..., 0.0495..., 0.0854...])
- soil_heat.gao_et_al.leuning_gz(g_zr: numpy.ndarray | float, z: numpy.ndarray | float, zr: numpy.ndarray | float, d: numpy.ndarray | float) numpy.ndarray | float[source]
Extrapolate soil heat flux from a reference depth z_r to a shallower (or deeper) target depth z using an exponential attenuation model (Leuning et al. 1985, Eq. 7).
The model assumes the amplitude of the heat-flux wave decays exponentially with depth at the same damping depth d that governs temperature attenuation:
\[G(z,t) \;=\; G(z_r,t)\;\exp\!\left(\frac{z_r - z}{d}\right)\]where \(G(z_r,t)\) is the measured flux at z_r (W m⁻²).
- Parameters:
g_zr (float or array_like) – Heat flux measured at reference depth z_r (W m⁻²). Can be a scalar or a NumPy-broadcastable array (e.g. a time series).
z (float or array_like) – Target depth z (m, positive downward). Must be broadcast- compatible with g_zr.
zr (float or array_like) – Reference depth z_r (m, positive downward). Broadcast compatible with z.
d (float or array_like) – Thermal damping depth (m). Positive; same shape rules as z.
- Returns:
Estimated heat flux at depth z (W m⁻²). Follows NumPy broadcasting rules applied to the inputs.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element of d is non-positive, or if inputs cannot be broadcast to a common shape.
Notes
Direction of extrapolation – z < z_r (shallower) ⇒ magnitude increases toward the surface; z > z_r (deeper) ⇒ magnitude decreases with depth.
Vectorisation – Inputs are converted to
numpy.ndarrayand the formula is evaluated element-wise.Limitations – The simple exponential form ignores soil layering and moisture variability; best suited to homogeneous profiles over the depth interval considered.
References
Leuning, R., Barlow, E. W. R., & Paltridge, G. (1985). Temperature gradients in a soil: Theory and measurement. Agricultural and Forest Meteorology, 35(1–4), 127-143. https://doi.org/10.1016/0168-1923(85)90067-3
Examples
>>> # Single depth conversion >>> leuning_gz(g_zr=-12.0, z=0.05, zr=0.08, d=0.07) -18.841...
>>> # Vectorized daily time series >>> g_plate = np.random.normal(-10, 3, 1440) # W m-2 @ 8 cm >>> d = 0.07 # m >>> Gz_5cm = leuning_gz(g_plate, z=0.05, zr=0.08, d=d) >>> Gz_5cm.shape (1440,)
- soil_heat.gao_et_al.simple_measurement_gz(g_zr: numpy.ndarray | float, cv_layers: numpy.ndarray | float, tz_layers: numpy.ndarray | float, dt: float, dz_layers: numpy.ndarray | float) numpy.ndarray | float[source]
Ground-heat-flux estimate at a target depth z using the simple-measurement variant of the calorimetric method (Liebethal & Foken 2007, Eq. 8).
The algorithm corrects a heat-flux‐plate measurement taken at a reference depth z_r by adding an approximation of the heat storage change ΔS in the soil slab between the plate and the surface, computed from two successive soil-temperature profiles:
\[G(z,t_j) \;=\; G(z_r,t_j) \;+\; \sum_{l=1}^{N} C_{v,l}\,\Delta z_l\,\]rac{Delta T_{l,j} + rac{1}{2}igl(Delta T_{l,j} - Delta T_{l,j-1}igr)}
{Delta t}
where
\(C_{v,l}\) — volumetric heat capacity of layer l (J m⁻³ K⁻¹)
\(\Delta z_l\) — layer thickness (m)
\(\Delta T_{l,j}\) — temperature change in layer l between t_{j-1} and t_j (K)
\(\Delta t\) — measurement interval (s)
The centred finite-difference term in brackets approximates the mid-interval temperature tendency, providing second-order accuracy from simple time-series measurements (tz_layers).
- g_zrfloat or array_like
Plate-measured soil heat flux at reference depth z_r (W m⁻²). Length M time steps.
- cv_layersarray_like
Volumetric heat capacity for each of N soil layers (J m⁻³ K⁻¹). Shape
(N,)or broadcast-compatible with the first axis of tz_layers.- tz_layersarray_like
Layer-average soil temperatures. Shape
(N, M)where N is the number of layers (matching cv_layers/dz_layers) and M is the number of time stamps (matching g_zr).- dtfloat
Constant time step between consecutive temperature samples (s). Must be positive.
- dz_layersarray_like
Thickness of each layer (m). Shape
(N,)broadcast-compatible with cv_layers.
- float or numpy.ndarray
Calorimetrically corrected soil heat flux at depth z (W m⁻²). Length M − 1 because the centred difference requires two successive profiles.
- ValueError
If input dimensions are inconsistent, dt ≤ 0, or any layer arrays contain non-finite values.
Temporal alignment – The first output value corresponds to the mid-point of the first two temperature samples
t[0] … t[1]; therefore the result is shifted ½ Δt relative to the plate-flux time stamps.Vectorisation – All calculations are fully vectorized using NumPy broadcasting. Inputs are internally cast to
numpy.ndarray.Applicability – Best suited to homogeneous layers with high-quality temperature measurements at identical timestamps.
Liebethal, C., & Foken, T. (2007). Evaluation of six parameterization approaches for the ground heat flux. Agricultural and Forest Meteorology, 143(1–2), 65–80. https://doi.org/10.1016/j.agrformet.2006.11.001
>>> g_plate = np.array([-12.3, -10.9, -8.5]) # W m-2 at z_r >>> Cv = np.array([2.3e6, 2.1e6]) # two layers >>> T = np.array([[15.0, 15.5, 16.0], # °C layer 1 ... [14.0, 14.2, 14.4]]) # °C layer 2 >>> dz = np.array([0.03, 0.02]) # m >>> Gz = simple_measurement_gz(g_plate, Cv, T, dt=1800, dz_layers=dz) >>> Gz array([-10.700..., -8.317...])
- soil_heat.gao_et_al.wbz12_g_gz(g_zr_series: numpy.ndarray | list | tuple, time_series: numpy.ndarray | list | tuple, z: float, zr: float, k_s_val: float) numpy.ndarray[source]
Estimate soil heat flux at a shallow depth z from a flux-plate record at reference depth z_r using the WBZ12-G convolution method (Wang, Bou-Zeid & Zhang 2012, their Eq. 9–10).
The algorithm inverts the one-dimensional heat-conduction equation for a semi-infinite homogeneous soil, assuming a step-response kernel based on the complementary error function erfc:
\[G(z,t) \;=\; 2\,G(z_r,t) \;-\; \frac{J(t)}{\Delta F_z(t)}\]with the discrete convolution integral
\[J(t_n) \;=\; \sum_{j=0}^{n-1} \tfrac{\bigl[G(z_r,t_{n-j-1}) + G(z_r,t_{n-j})\bigr]}{2}\; \Delta F_z(t_{j+1})\]and the transfer function increment
\[\Delta F_z(t) = \operatorname{erfc}\! \left[\frac{z_r - z}{2\sqrt{k_s (t - t_0)}}\right] \;-\; \operatorname{erfc}\! \left[\frac{z_r - z}{2\sqrt{k_s (t - t_0 - \Delta t)}}\right].\]- Parameters:
g_zr_series (array_like) – Time series of heat-flux-plate measurements at z_r (W m⁻²); length N.
time_series (array_like) – Monotonically increasing time stamps (s) corresponding to g_zr_series. Must be the same length N.
z (float) – Target depth (m, positive downward).
zr (float) – Reference depth of the heat-flux plate (m, positive downward).
k_s_val (float) – Soil thermal diffusivity k_s (m² s⁻¹). Must be positive.
- Returns:
Estimated heat-flux series at depth z, same length N as the input series (W m⁻²).
- Return type:
numpy.ndarray
- Raises:
ValueError – If k_s_val ≤ 0, the two series differ in length, time stamps are non-monotonic, or contain fewer than two samples.
Notes
The first element of the output corresponds to t₀ (no storage correction possible yet); subsequent points incorporate the convolution up to that instant.
A very small term
+ 1e-12is added under the square root to avoid division by zero at t = t₀.The complementary error function is evaluated with
numpy.erfc(), which is vectorized and avoids the SciPy dependency.
References
Wang, W., Bou-Zeid, E., & Zhang, Y. (2012). Estimating surface heat fluxes using the surface renewal method. Boundary-Layer Meteorology, 144(2), 407-422. https://doi.org/10.1007/s10546-012-9730-9
Examples
>>> # 10-min sampled day-long plate record at zr = 0.08 m >>> t = np.arange(0, 24*3600 + 1, 600) # s >>> Gp = -10 + 5*np.sin(2*np.pi*t/86400) >>> Gz = wbz12_g_gz(Gp, t, z=0.05, zr=0.08, k_s_val=1.0e-6) >>> Gz.shape (145,)
- soil_heat.gao_et_al.wbz12_s_gz(Ag: numpy.ndarray | float, ks_val: numpy.ndarray | float, zr: float, z: float, t: numpy.ndarray | float, eps: numpy.ndarray | float, omega: float = OMEGA_DAY) numpy.ndarray | float[source]
WBZ12-S analytic–numeric solution for soil heat flux G(z, t) at depth z (Wang, Bou-Zeid & Zhang 2012, Eq. 11).
WBZ12-S combines a closed-form forcing term (sinusoidal surface temperature) with a storage term that is expressed as a Fourier‐ Laplace integral and must be evaluated numerically. The governing equation assumes homogeneous soil properties and a sinusoidal surface temperature of amplitude A_g and phase eps:
\[G(z,t) \;=\; A_g\,e^{-\eta}\,\sin\bigl(\omega t +\]- arepsilon - etabigr)
;-; frac{2 A_g k_s}{pi} int_{0}^{infty}
- frac{
k_s zeta^{2} igl[sinvarepsilon - omega cosvarepsilonbigr]
- }{
omega^{2} + k_{s}^{2} zeta^{4}
} sinbigl[zeta,(z_r - z)bigr], e^{-k_s zeta^{2} t},dzeta
with
\[\eta \;=\; (z_r - z)\,\sqrt{\omega / (2 k_s)}.\]The first term (prefactor × sin_term) represents the periodic steady-state component, while the integral accounts for the transient adjustment of the soil profile.
- Agfloat or array_like
Amplitude of the sinusoidal surface/ground temperature (K). Can be broadcast with t and eps.
- ks_valfloat or array_like
Soil thermal diffusivity \(k_s\) (m² s⁻¹). Must be positive.
- zrfloat
Reference depth of the heat-flux plate (m, positive downward).
- zfloat
Target depth for the output heat flux (m, positive downward).
- tfloat or array_like
Time (s) since the start of the forcing cycle. Broadcast- compatible with Ag and eps.
- epsfloat or array_like
Phase shift \(\varepsilon\) (rad) of the surface temperature. Broadcast-compatible with t.
- omegafloat, optional
Angular frequency \(\omega\) (rad s⁻¹) of the sinusoidal forcing (default is :pydata:`OMEGA_DAY` ≈ 7.272 × 10⁻⁵ s⁻¹).
- float or numpy.ndarray
Soil heat flux G(z, t) (W m⁻²). The shape matches NumPy broadcasting over (Ag, ks_val, t, eps).
- ValueError
If ks_val or omega are non-positive, or if the inputs cannot be broadcast to a common shape.
Numerical integral – The transient term is evaluated using
scipy.integrate.quad()over \(\zeta \in [0, \infty)\). The integral can be expensive for long time-series; cache results or vectorise with more sophisticated quadrature if performance is critical.Vectorisation – quad is called individually for every element in the broadcasted inputs; large arrays may thus incur heavy computational cost.
Units – Use consistent SI units: W m⁻² (output), K (Ag), m² s⁻¹ (ks_val), rad s⁻¹ (omega), seconds (t), metres (depths).
Wang, W., Bou-Zeid, E., & Zhang, Y. (2012). Estimating surface heat fluxes using the surface renewal method. Boundary-Layer Meteorology, 144(2), 407–422. https://doi.org/10.1007/s10546-012-9730-9
>>> # Daily sinusoid, single point >>> G = wbz12_s_gz( ... Ag=8.0, ks_val=1.0e-6, ... zr=0.08, z=0.05, ... t=np.linspace(0, 86400, 97), ... eps=0.0 ... ) >>> G.shape (97,)
- soil_heat.gao_et_al.exact_temperature_gz(z: numpy.ndarray | float, AT: numpy.ndarray | float, t: numpy.ndarray | float, d: numpy.ndarray | float, omega: float = OMEGA_DAY, T_i: float = 298.15) numpy.ndarray | float[source]
Compute the exact analytical soil‐temperature profile under sinusoidal surface forcing.
A sinusoidal temperature wave propagating downward through a homogeneous soil decays exponentially with depth and lags in phase. The closed‐form solution at depth z (m) is
\[T(z, t) \;=\; T_i \;+\; A_T \, e^{-z/d}\, \sin\bigl(\omega t \;-\; z/d\bigr),\]where
\(T_i\) — initial (mean) temperature of the soil column (K)
\(A_T\) — amplitude of the surface-temperature oscillation (K)
\(d\) — thermal damping depth (m)
\(\omega\) — angular frequency of the forcing (rad s⁻¹)
\(t\) — time since the start of oscillation (s)
- Parameters:
z (float or array_like) – Depth below the surface (m). Positive downward. Can be a scalar or an array broadcastable with t and AT.
AT (float or array_like) – Amplitude A_T of the surface temperature wave (K).
t (float or array_like) – Time variable (s). Supports vectorisation.
d (float or array_like) – Thermal damping depth (m). Must be positive. Can be scalar or broadcastable with z.
omega (float, optional) – Angular frequency ω (rad s⁻¹). Defaults to :pydata:`OMEGA_DAY` (≈ 7.272 × 10⁻⁵ s⁻¹ for 24 h).
T_i (float, optional) – Mean (initial) soil temperature (K). Defaults to 298.15 K (25 °C).
- Returns:
Soil temperature T(z, t) (same units as T_i). The result shape matches NumPy broadcasting over the inputs.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element of d or omega is non-positive, or if inputs cannot be broadcast to a common shape.
Notes
Phase lag – Each e-folding depth d introduces a phase delay of 1 radian (~57.3 °) relative to the surface signal.
Vectorisation – All inputs are converted to
numpy.ndarray; the formula is evaluated element-wise and fully supports broadcasting.Units – Ensure consistent SI units (metres, seconds, kelvins).
References
Adapted from Gao, Z., Horton, R., Luo, L., & Kucharik, C. J. (2010). A simple method to measure soil temperature dynamics: Theory and application. Soil Science Society of America Journal, 74(2), 580-588. https://doi.org/10.2136/sssaj2009.0169
Examples
>>> # Daily cycle at 10 cm depth >>> z = 0.10 # m >>> AT = 8.0 # K >>> d = 0.12 # m >>> t_day = np.linspace(0, 86400, 97) # 15-min resolution >>> Tz = exact_temperature_gz(z, AT, t_day, d) >>> Tz.shape (97,)
- soil_heat.gao_et_al.exact_gz(z: numpy.ndarray | float, AT: numpy.ndarray | float, lambda_s_val: numpy.ndarray | float, d: numpy.ndarray | float, t: numpy.ndarray | float, omega: float = OMEGA_DAY) numpy.ndarray | float[source]
Exact analytical soil-heat-flux solution for sinusoidal surface forcing.
A sinusoidal surface–temperature wave of amplitude A_T and angular frequency ω generates a heat-flux wave that decays exponentially with depth and lags the temperature signal by π / 4 radians (Gao et al., 2010, Eq. 17):
\[G(z,t) \;=\; \sqrt{2}\; \frac{\lambda_s A_T}{d}\; e^{-z/d}\; \sin\bigl(\omega t - z/d + \tfrac{\pi}{4}\bigr)\]where
\(\lambda_s\) — soil thermal conductivity (W m⁻¹ K⁻¹)
\(d\) — thermal damping depth (m)
\(z\) — depth below the surface (m, positive downward)
\(t\) — time (s) since the start of oscillation
\(\omega\) — angular frequency (rad s⁻¹)
- Parameters:
z (float or array_like) – Depth(s) below the surface (m), positive downward.
AT (float or array_like) – Amplitude A_T of the surface-temperature oscillation (K).
lambda_s_val (float or array_like) – Soil thermal conductivity λ_s (W m⁻¹ K⁻¹).
d (float or array_like) – Thermal damping depth (m). Must be positive.
t (float or array_like) – Time variable (s).
omega (float, optional) – Angular frequency ω (rad s⁻¹). Defaults to :pydata:`OMEGA_DAY` (≈ 7.272 × 10⁻⁵ s⁻¹, i.e. 2π / 86 400 s).
- Returns:
Soil heat flux G(z, t) (W m⁻²). The return shape follows NumPy broadcasting rules applied to the inputs.
- Return type:
float or numpy.ndarray
- Raises:
ValueError – If any element of d or omega is non-positive, or if inputs cannot be broadcast to a common shape.
Notes
Phase shift – At any given depth, the heat-flux wave lags the temperature wave by 45 ° (π / 4 rad).
Vectorisation – All inputs are converted to
numpy.ndarray; the expression is evaluated element-wise and fully supports broadcasting.Units – Ensure all quantities use consistent SI units.
References
Gao, Z., Horton, R., Luo, L., & Kucharik, C. J. (2010). A simple method to measure soil temperature dynamics: Theory and application. Soil Science Society of America Journal, 74(2), 580–588. https://doi.org/10.2136/sssaj2009.0169
Examples
>>> # Scalar example >>> exact_gz( ... z=0.05, ... AT=8.0, ... lambda_s_val=1.2, ... d=0.12, ... t=3600, ... ) 21.52...
>>> # Vectorized daily cycle >>> t_day = np.linspace(0, 86400, 97) # 15-min resolution >>> Gz = exact_gz( ... z=0.10, ... AT=6.0, ... lambda_s_val=1.1, ... d=0.11, ... t=t_day, ... ) >>> Gz.shape (97,)