soil_heat
This package provides a collection of Python functions for calculating soil heat flux and related soil thermal properties. The implementations are based on various models and parameterizations published in the peer-reviewed scientific literature.
The library is organized into modules, each corresponding to a specific publication or a set of related equations.
Available submodules:
soil_heat: Core functions and utilities.
gao_et_al: Functions from Gao et al. (2017).
liebethal_and_folken: Functions from Liebethal & Foken (2006).
wang_and_bouzeid: Functions from Wang & Bou-Zeid (2012).
wang_and_yang: Functions from Yang & Wang (2008).
fao56_soil_heat_flux: FAO-56 and ASCE soil heat flux methods.
storage_calculations: Soil and canopy heat storage calculations.
johansen: Soil thermal property models (Johansen, Lu et al.) and advanced heat flux methods.
All functions are designed to work with NumPy arrays for efficient, vectorized computations.
Submodules
Attributes
|
Calculate the residual of the surface energy budget (Yang & Wang 2008, Eq. 12). |
Functions
|
Estimate near-surface soil heat flux using Fourier’s law. |
|
Calculate surface soil heat flux via the calorimetric (heat-storage) method. |
|
Compute the vertical temperature gradient between two sensors. |
|
Compute heat flux G at cell interfaces using Fourier’s law. |
|
Estimate soil volumetric heat capacity using the de Vries (1963) linear |
|
Convert thermal diffusivity ( |
|
Compute the daily diurnal amplitude of a time-series. |
|
Compute the daily peak‐time lag (Δt) between two diurnal signals. |
|
Fit a sinusoidal model to time–series data using non-linear least |
|
Evaluate a sinusoidal wave of the form |
|
Estimate soil thermal diffusivity ( |
|
Estimate soil thermal diffusivity (α) from the phase lag of a temperature wave. |
|
Estimate soil thermal diffusivity ( |
|
Estimate soil thermal diffusivity ( |
|
Estimate soil thermal diffusivity ( |
|
Compute thermal diffusivity, thermal conductivity, and **soil heat |
|
Estimate the volumetric heat capacity of dry soil ( |
|
Calculate soil heat storage and heat flux. |
|
Compute the soil thermal conductivity \(\lambda_s\) |
|
Calculate the soil thermal diffusivity \(k_s\) |
|
Estimate soil volumetric heat capacity using the de Vries (1963) linear |
|
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. |
|
Compute the reference ground‑heat flux G₀,M using the |
|
Estimate ground heat flux as a fixed fraction of net radiation |
|
Estimate ground heat flux using a linear regression against net |
|
Compute the parameters A and B for the universal net-radiation |
|
Estimate ground heat flux using the universal net-radiation |
|
Estimate the diurnal surface-temperature amplitude (ΔT_s) from |
|
Calculate the empirical parameter φ based on soil moisture content |
|
Estimate ground heat flux from the sensible heat flux (H) |
|
Estimate surface ground heat flux using the "simple measurement" |
|
Calculate the thickness of the active soil layer (δz) for the |
|
Estimate ground heat flux using the two-layer force-restore method |
|
Compute the closure residual of the surface energy balance (SEB). |
|
Conventional estimate of surface ground heat flux, Eq. (2). |
|
Green-function solution \(g_z(t)\) for the **one‐dimensional, |
|
Temperature time-series at depth z via Duhamel convolution (Eq. 6). |
|
Compute G(z,t) from a known surface flux series G0 (Eq. 9). |
|
Estimate surface ground heat flux G0 from plate measurements Gz. |
|
Evaluate a sinusoidal surface heat-flux forcing |
|
Analytical solution for soil temperature beneath a sinusoidally |
|
Analytical soil heat–flux response G(z,t) to a sinusoidal |
|
Volumetric heat capacity of moist soil, Eq. (16). |
|
Return Pf (Eq. 18) from volumetric water content. |
|
Thermal conductivity parameterisation, Eq. (17). |
|
Return κ = k / (ρ c). |
|
Compute heat flux G at cell interfaces using Fourier’s law. |
|
Calculate the soil heat flux profile by integrating the change in |
|
Estimate soil volumetric heat capacity using the de Vries (1963) linear |
|
Generate layer thicknesses for a vertically stretched grid. |
|
Solve the 1D thermal diffusion equation for one time step using an |
|
Correct a modeled temperature profile using observed temperatures. |
|
Calculate surface temperature from upward and downward long-wave |
|
Estimate soil thermal conductivity based on soil moisture and dry |
|
Calculate the diagnostic error in heat flux that arises from assuming |
|
Calculate the residual of the surface energy budget (Yang & Wang 2008, Eq. 12). |
|
Perform a single integration step of the TDEC (Temperature Diffusion |
|
Soil heat flux from the general calorimetric equation (FAO-56 Eq. 41). |
|
Daily soil heat flux for a grass reference surface (FAO-56 Eq. 42). |
|
Monthly soil heat flux using adjacent months (FAO-56 Eq. 43). |
|
Monthly soil heat flux — next month unavailable (FAO-56 Eq. 44). |
|
Hourly (sub-daily) soil heat flux as a fraction of net radiation. |
|
Hourly soil heat flux, auto-detecting day/night from Rn sign. |
|
Compute soil heat flux for each month in a temperature time series. |
|
Monthly G for a complete 12-month annual cycle using Eq. 43. |
|
Convert energy flux to equivalent evaporation depth (FAO-56 Eq. 20). |
|
Convert equivalent evaporation depth back to energy flux. |
|
Convert flux from MJ m⁻² day⁻¹ to W m⁻². |
|
Convert flux from W m⁻² to MJ m⁻² day⁻¹. |
|
Reproduce FAO-56 Example 13 and other verification cases. |
|
Calculate the soil heat storage flux (S) in W/m² using a change-in-storage |
|
Calculate the energy storage flux of a plant canopy (S_canopy) in W/m². |
|
Estimate soil volumetric heat capacity using the de Vries (1963) linear |
|
Estimate soil thermal conductivity using the Johansen (1975) model with the |
|
Estimate soil thermal conductivity using the Lu et al. (2007) empirical model. |
|
Estimate surface ground heat flux using a temperature gradient from the canopy |
|
Compute soil heat storage in the layer from the surface to the heat flux plate depth. |
|
Estimate surface ground heat flux using a conductive gradient at a reference depth |
Package Contents
- soil_heat.compute_heat_flux_conduction(df: pandas.DataFrame, depth1: float = 0.05, depth2: float = 0.1, col_T1: str = 'T5cm', col_T2: str = 'T10cm', col_theta1: str = 'VWC5cm', col_theta2: str = 'VWC10cm', porosity: float = 0.4, k_dry: float = 0.25, k_sat: float = 1.5) pandas.Series[source]
Estimate near-surface soil heat flux using Fourier’s law.
This “gradient” approach computes conductive ground-heat flux \(G\) between two depths by multiplying the vertical temperature gradient with an effective thermal conductivity that varies with volumetric water content (VWC).
- Parameters:
df (pandas.DataFrame) – Time-indexed data containing at least the four columns specified by col_T1, col_T2, col_theta1, and col_theta2. The index spacing defines the temporal resolution of the output.
depth1 (float, default (0.05, 0.10)) – Sensor depths (m). depth2 must be greater (deeper) than depth1.
depth2 (float, default (0.05, 0.10)) – Sensor depths (m). depth2 must be greater (deeper) than depth1.
col_T1 (str, default ("T5cm", "T10cm")) – Column names for temperature (°C or K) at depth1 and depth2.
col_T2 (str, default ("T5cm", "T10cm")) – Column names for temperature (°C or K) at depth1 and depth2.
col_theta1 (str, default ("VWC5cm", "VWC10cm")) – Column names for volumetric water content (m³ m⁻³) at depth1 and depth2.
col_theta2 (str, default ("VWC5cm", "VWC10cm")) – Column names for volumetric water content (m³ m⁻³) at depth1 and depth2.
porosity (float, default 0.40) – Soil total porosity (saturated VWC, m³ m⁻³).
k_dry (float, default 0.25) – Dry-soil thermal conductivity (W m⁻¹ K⁻¹).
k_sat (float, default 1.50) – Saturated-soil thermal conductivity (W m⁻¹ K⁻¹).
- Returns:
Half-hourly (or whatever the index step is) ground-heat-flux series with name
"G_conduction". Units are W m⁻². Positive values indicate downward flux.- Return type:
pandas.Series
Notes
The effective thermal conductivity is computed by a simple linear mixing model:
\[\lambda_\text{eff} = k_\text{dry} + \frac{\bar{\theta}}{\phi} \bigl(k_\text{sat} - k_\text{dry}\bigr),\]where \(\bar{\theta}\) is the mean VWC of the two depths and \(\phi\) is porosity. More sophisticated models (e.g. Johansen, de Vries) can be substituted if site-specific calibration is available.
References
Campbell & Norman (2012) An Introduction to Environmental Biophysics, ch. 7.
Gao et al. (2017) Agricultural and Forest Meteorology, 240 – 241, 194–204.
Examples
>>> G = compute_heat_flux_conduction(df_site, ... depth1=0.05, depth2=0.10, ... col_T1="T_05", ... col_T2="T_10", ... col_theta1="VWC_05", ... col_theta2="VWC_10") >>> G.plot(title="Soil heat flux (gradient method)")
- soil_heat.compute_heat_flux_calorimetric(df: pandas.DataFrame, depth_levels: list[float], T_cols: list[str], theta_cols: list[str], C_dry: float = 2100000.0, C_w: float = 4200000.0) pandas.Series[source]
Calculate surface soil heat flux via the calorimetric (heat-storage) method.
The calorimetric method integrates the transient change in heat storage within a multilayer soil column. For a surface-to-depth layer of thickness \(z_{\text{ref}}\), the surface flux \(G_0\) is approximated by
\[G_0 \;\approx\; \frac{\Delta Q}{\Delta t} \;=\; \frac{1}{\Delta t} \sum_{i=1}^{N_\text{layers}} C_i \, \Delta T_i \, \Delta z_i,\]where \(C_i\) is volumetric heat capacity (J m⁻³ K⁻¹), \(\Delta T_i\) is the average temperature change (K) in layer i, and \(\Delta z_i\) is layer thickness (m). No heat-flux-plate reading is required if the deepest measurement depth lies below the diurnal damping depth such that \(G(z_{\text{ref}}) \approx 0\).
- Parameters:
df (pandas.DataFrame) – Time-indexed data containing temperature and VWC columns for all depths specified in T_cols and theta_cols. Index spacing sets the output time step.
depth_levels (list of float) – Depths (m) corresponding in order to the entries in T_cols and theta_cols. Must be strictly increasing.
T_cols (list of str) – Column names for soil temperatures (°C or K) at depth_levels.
theta_cols (list of str) – Column names for volumetric water content (m³ m⁻³) at depth_levels.
C_dry (float, default 2.1e6) – Volumetric heat capacity of dry soil matrix (J m⁻³ K⁻¹).
C_w (float, default 4.2e6) – Volumetric heat capacity of liquid water (J m⁻³ K⁻¹).
- Returns:
Surface ground-heat-flux series,
"G_calorimetric"(W m⁻²). Positive values denote downward flux. The first time step is set to NaN because a preceding interval is required.- Return type:
pandas.Series
Notes
Heat capacity model
A simple two-component mixture is assumed:
\[C = (1 - \theta)\,C_{\text{dry}} + \theta\,C_w.\]If bulk density or mineral fraction data are available, replace this linear approximation with a mass-weighted formulation.
Boundary assumption
The deepest temperature is treated as a “no-flux” boundary (storage only). If diurnal waves penetrate deeper at your site, include an additional flux-plate term or extend depth_levels downward.
References
Mayocchi & Bristow (1995) Agricultural and Forest Meteorology 75, 93–109.
Oke (2002) Boundary-Layer Climates, 2nd ed., §2.3.
Fluxnet2015 “G” best-practice guide (https://fluxnet.org/sites/default/files/soil_heat_flux_guide.pdf).
Examples
>>> depths = [0.05, 0.10, 0.20, 0.50] # m >>> Tcols = ["T5", "T10", "T20", "T50"] # °C >>> Vcols = ["VWC5", "VWC10", "VWC20", "VWC50"] >>> G0 = compute_heat_flux_calorimetric(df_site, ... depths, Tcols, Vcols) >>> G0.resample("D").mean().plot() >>> plt.ylabel("Daily mean G₀ (W m$^{-2}$)")
- soil_heat.temperature_gradient(T_upper: numpy.ndarray | float, T_lower: numpy.ndarray | float, depth_upper: float, depth_lower: float) numpy.ndarray | float[source]
Compute the vertical temperature gradient between two sensors.
The gradient is defined as the change in temperature divided by the change in depth (positive downward):
\[\frac{∂T}{∂z} \;=\; \frac{T_{\text{lower}} - T_{\text{upper}}} {z_{\text{lower}} - z_{\text{upper}}} \;\;[^{\circ}\text{C m}^{-1}]\]- Parameters:
T_upper (float or array_like) – Temperature at the shallower depth
depth_upper(°C).T_lower (float or array_like) – Temperature at the deeper depth
depth_lower(°C). Must be broadcast-compatible withT_upper.depth_upper (float) – Depth of the upper sensor (m, positive downward).
depth_lower (float) – Depth of the lower sensor (m, positive downward). Must satisfy
depth_lower > depth_upperfor a meaningful gradient.
- Returns:
Temperature gradient ∂T/∂z (°C m⁻¹). Shape follows NumPy broadcasting of
T_upperandT_lower.- Return type:
ndarray or float
- Raises:
ValueError – If
depth_lower≤depth_upper.
Notes
Sign convention – A positive gradient indicates temperatures increase with depth (warmer below).
Vectorised – The arithmetic is fully NumPy-broadcasted; use it on scalar values, 1-D arrays, or entire DataFrames’ columns.
Units – Because depth is in metres and temperature in degrees Celsius, the result is °C m⁻¹ (identical to K m⁻¹).
Examples
>>> grad = temperature_gradient( ... T_upper=18.6, T_lower=20.1, ... depth_upper=0.05, depth_lower=0.10, ... ) >>> print(f"Gradient = {grad:.2f} °C/m") Gradient = 30.00 °C/m
Array input:
>>> T_up = np.array([15.0, 16.2, 17.1]) >>> T_low = np.array([14.0, 15.8, 16.9]) >>> temperature_gradient(T_up, T_low, 0.02, 0.08) array([-16.66666667, -6.66666667, -3.33333333]) # °C/m
- soil_heat.soil_heat_flux(T_upper, T_lower, depth_upper, depth_lower, k)[source]
Calculate soil heat flux (G) using Fourier’s law of heat conduction.
This function computes the one-dimensional heat flux in the soil based on the temperature gradient between two depths and the soil’s thermal conductivity. The formula used is:
\[G = -k \frac{\Delta T}{\Delta z}\]where \(G\) is the soil heat flux, \(k\) is the thermal conductivity, \(\Delta T\) is the temperature difference, and \(\Delta z\) is the distance between the two measurement depths.
- Parameters:
T_upper (float or array_like) – Temperature at the upper depth (°C or K).
T_lower (float or array_like) – Temperature at the lower depth (°C or K).
depth_upper (float) – Depth of the upper sensor (m, positive downward).
depth_lower (float) – Depth of the lower sensor (m, positive downward).
k (float or array_like) – Thermal conductivity of the soil layer between the sensors (W m⁻¹ K⁻¹).
- Returns:
The calculated soil heat flux (W m⁻²). A positive value indicates a downward flux (from the surface into the soil).
- Return type:
float or array_like
- soil_heat.volumetric_heat_capacity(theta_v)[source]
Estimate volumetric heat capacity (Cv) from soil moisture content.
This function uses a simple mixing model to estimate the volumetric heat capacity of moist soil. It assumes the soil is a two-component mixture of dry soil and water. The formula is:
\[C_v = (1 - \theta_v) C_{soil} + \theta_v C_{water}\]where \(\theta_v\) is the volumetric water content, and \(C_{soil}\) and \(C_{water}\) are the volumetric heat capacities of dry soil and water, respectively.
- soil_heat.thermal_conductivity(alpha: numpy.ndarray | float, theta_v: numpy.ndarray | float) numpy.ndarray | float[source]
Convert thermal diffusivity (
α) to thermal conductivity (k) using the bulk volumetric heat capacity of moist soil.The relationship is
\[k \;=\; α \, C_v(θ_v),\]where
α – thermal diffusivity (m² s⁻¹),
C_v(θ_v) – volumetric heat capacity (J m⁻³ K⁻¹) as a function of volumetric water content θ_v (m³ m⁻³). It is obtained from
volumetric_heat_capacity().
- Parameters:
- Returns:
Thermal conductivity k (W m⁻¹ K⁻¹) with the broadcast shape of the inputs.
- Return type:
ndarray or float
Notes
Volumetric heat capacity model –
volumetric_heat_capacity()typically assumes a two‐phase mixture of mineral soil and water:\[C_v(θ_v) \;=\; (1-θ_v)\,ρc_ ext{dry} \;+\; θ_v\,ρc_ ext{w} ,\]where
ρc_dry(≈ 2.0 MJ m⁻³ K⁻¹) andρc_w(4.18 MJ m⁻³ K⁻¹) are the volumetric heat capacities of dry soil and liquid water, respectively. Ensure these defaults suit your substrate.Vectorisation – The function is a one‐liner,
alpha * Cv, and thus inherits full NumPy broadcasting rules.Temperature units – Because heat capacity is per kelvin, k is returned in W m⁻¹ K⁻¹ (equivalent to W m⁻¹ °C⁻¹).
Examples
>>> α = np.array([1.4e-7, 1.6e-7]) # m² s⁻¹ >>> θ = np.array([0.10, 0.25]) # m³ m⁻³ >>> k = thermal_conductivity(α, θ) >>> k array([0.29, 0.54]) # W m⁻¹ K⁻¹
Plot conductivity versus moisture:
>>> θ_range = np.linspace(0, 0.45, 100) >>> k_vals = thermal_conductivity(1.5e-7, θ_range) >>> plt.plot(θ_range, k_vals) >>> plt.xlabel("Volumetric water content (m³ m⁻³)") >>> plt.ylabel("Thermal conductivity (W m⁻¹ K⁻¹)")
- soil_heat.diurnal_amplitude(series: pandas.Series) pandas.Series[source]
Compute the daily diurnal amplitude of a time-series.
The diurnal amplitude for a given calendar day is defined as the difference between that day’s maximum and minimum values:
\[A_d \;=\; \max\_{t \in d} x(t) \;-\; \min\_{t \in d} x(t)\]This metric is frequently used for temperature, soil-heat, or other environmental data to characterise the strength of the diurnal cycle.
- Parameters:
series (pandas.Series) – Time-indexed observations with a DatetimeIndex. Any frequency is accepted, but the index must be sorted and monotonic. Missing values (NaN) are ignored within each daily window.
- Returns:
Daily diurnal amplitude, indexed by date (midnight
00:00of each day). Units are the same as those of the inputseries.- Return type:
pandas.Series
Notes
Resampling rule – The computation uses
>>> daily_max = series.resample("D").max() >>> daily_min = series.resample("D").min()
which bins data by calendar day in the series’ timezone. Incomplete trailing days yield NaN.
Timezone safety – If the series’ index spans daylight-saving transitions, consider converting to UTC prior to analysis to avoid artificial jumps in daily windows.
Robustness – For noisy signals, you may wish to smooth
series(e.g. rolling median) before calling this function.
Examples
>>> amp = diurnal_amplitude(df["air_temperature"]) >>> amp.plot(title="Daily Temperature Amplitude") >>> amp.describe().loc[["min", "mean", "max"]] min 4.3 mean 9.7 max 15.2 Name: air_temperature, dtype: float64
- soil_heat.diurnal_peak_lag(series1: pandas.Series, series2: pandas.Series) pandas.Series[source]
Compute the daily peak‐time lag (Δt) between two diurnal signals.
For each calendar day the function identifies the clock time at which each series reaches its maximum value and returns the signed time difference in hours (
series1minusseries2). A modular correction confines the result to the interval[-12, 12]h so that, for example, a raw lag of –23 h becomes +1 h.- Parameters:
series1 (pandas.Series) – Time-indexed observations of equal length, preferably temperature or some other quantity exhibiting a clear diurnal cycle. The index must be DatetimeIndex and should be timezone-aware and aligned in frequency. Missing values are ignored within each daily resampling window.
series2 (pandas.Series) – Time-indexed observations of equal length, preferably temperature or some other quantity exhibiting a clear diurnal cycle. The index must be DatetimeIndex and should be timezone-aware and aligned in frequency. Missing values are ignored within each daily resampling window.
- Returns:
Daily peak-lag values (float, hours) indexed by the date of the peak (00:00 of each day). Positive lags mean the peak of
series1occurs later than the peak ofseries2on that day; negative lags indicate the opposite.- Return type:
pandas.Series
Notes
Resampling rule – Peaks are detected with
series.resample('D').apply(lambda x: x.idxmax()). Ensure the input data span whole days; incomplete trailing days yield NaN.Wrap-around correction – The transformation
(lag + 12) % 24 − 12folds lags so that the maximum absolute value is always < 12 h, which prevents a late-evening peak at 23:30 and an early-morning peak at 00:30 from being reported as –23 h.Daylight-saving – If the index carries a timezone subject to DST transitions, consider converting to UTC prior to analysis to avoid spurious 1-h jumps.
Examples
>>> peak_lag = diurnal_peak_lag(df['ts_05cm'], df['ts_10cm']) >>> peak_lag.describe() count 90.000000 mean 1.42 std 0.53 min -0.73 25% 1.11 50% 1.42 75% 1.74 max 2.33 Name: ts_05cm, dtype: float64
Plot the distribution:
>>> peak_lag.plot(kind='hist', bins=24) >>> plt.xlabel('Peak lag (h)') >>> plt.title('Daily phase lag: 5 cm vs 10 cm temperature')
- soil_heat.fit_sinusoid(t: numpy.ndarray, data: numpy.ndarray) tuple[numpy.ndarray, numpy.ndarray][source]
Fit a sinusoidal model to time–series data using non-linear least squares.
The model is
\[y(t)\;=\;A \sin( \omega t +\]arphi ) + C ,
where
Ais the amplitude,ωthe angular frequency,φthe phase shift, andCthe vertical offset. Initial parameter guesses are derived from the sample statistics of data and an assumed daily frequency.- tndarray
1-D array of time stamps (s). Must be the same length as
data.- datandarray
Observed values corresponding to
t(e.g. temperature, °C).
- poptndarray, shape (4,)
Optimal parameters
[A, ω, φ, C]that minimise the sum-of-squares error.- pcovndarray, shape (4, 4)
Covariance matrix of the parameter estimates returned by
scipy.optimize.curve_fit().
Initial guess – Amplitude is set to the sample standard deviation of data, frequency to a 24-h cycle (
ω = 2π / 86 400s⁻¹), phase to 0, and offset to the sample mean. Adjust these if fitting to non-diurnal signals.Robustness – If convergence issues arise, provide a closer initial guess or bound parameters via
curve_fit’sbounds=keyword.
>>> import numpy as np >>> from scipy.optimize import curve_fit >>> t = np.arange(0, 3*86400, 1800) # 3 days, 30-min Δt >>> true = sinusoid(t, 7, 2*np.pi/86400, 0.3, 15) >>> rng = np.random.default_rng(0) >>> y = true + rng.normal(0, 0.5, t.size) # add noise >>> params, _ = fit_sinusoid(t, y) >>> A, ω, φ, C = params >>> print(f"Amplitude={A:.2f}, Period={2*np.pi/ω/3600:.2f} h") Amplitude=7.01, Period=24.00 h
- soil_heat.sinusoid(t: numpy.ndarray | float, A: float, omega: float, phase: float, offset: float) numpy.ndarray | float[source]
Evaluate a sinusoidal wave of the form
\[f(t) \;=\; A \sin(\omega\, t +\]arphi) + C ,
where \(A\) is the amplitude, \(\omega\) the angular frequency, :math:`
- arphi` the phase shift, and \(C\)
a constant vertical offset.
The function is vectorised with NumPy broadcasting, so
tmay be a scalar, 1-D array, or any shape compatible with the parameters.- tfloat or array_like
Independent variable (time, angle, etc.). Units are arbitrary, but must be consistent with
omega(e.g. seconds ifomegais rad s⁻¹).- Afloat
Wave amplitude. Sets the peak deviation from
offset.- omegafloat
Angular frequency (rad ×
t``⁻¹). For a *temporal* signal ``ω = 2π / Pwhere P is the period.- phasefloat
Phase shift :math:`
- arphi` in radians.
Positive values delay the wave (right shift), negative values advance it.
- offsetfloat
Constant vertical shift \(C\). Often the long-term mean or base-line value of the signal.
- ndarray or float
Value(s) of the sinusoid at
twith the same shape as the broadcast result of the inputs.
Vectorisation – Internally relies on
numpy.sin; all standard broadcasting rules apply.Period – The fundamental period P is related to
omegaby P = 2π / ω. Specify ω rather than P to avoid repeated division operations when fitting.
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> t = np.linspace(0, 24, 1000) # hours >>> temp = sinusoid(t, A=6, omega=2*np.pi/24, phase=0, offset=15) >>> plt.plot(t, temp) >>> plt.xlabel("Time (h)") >>> plt.ylabel("Temperature (°C)") >>> plt.title("Idealised diurnal temperature wave") >>> plt.show()
Fit a sinusoid to noisy data with
scipy.optimize.curve_fit():>>> from scipy.optimize import curve_fit >>> rng = np.random.default_rng(42) >>> y_obs = sinusoid(t, 6, 2*np.pi/24, 0.2, 15) + rng.normal(0, 0.5, t.size) >>> p0 = (5, 2*np.pi/24, 0, 15) # initial guess >>> popt, _ = curve_fit(sinusoid, t, y_obs, p0=p0) >>> amp, omg, ph, off = popt >>> print(f"Amplitude = {amp:.2f}, Phase = {ph:.2f} rad")
- soil_heat.thermal_diffusivity_amplitude(A1: float, A2: float, z1: float, z2: float, period: int = 86400) float[source]
Estimate soil thermal diffusivity (
α) from the damping of harmonic amplitude between two depths.A one–dimensional soil column subject to a sinusoidal surface temperature oscillation exhibits an exponential decay of amplitude with depth (Carslaw & Jaeger, 1959). For a single angular frequency \(ω = 2π/P\), the analytical solution yields
\[α \;=\; \frac{π\, (z_2 - z_1)^2} {P \;\bigl[\,\ln(A_1/A_2)\bigr]^2} ,\]where
A₁ and A₂ are the harmonic amplitudes at depths z₁ and z₂, respectively (A₁ > A₂),
P is the forcing period, and
z₂ – z₁ is the vertical separation of the two sensors.
- Parameters:
A1 (float) – Diurnal (or other fundamental) temperature amplitudes at the shallow depth
z1and deeper depthz2. Units °C or K (identical for both).A2 (float) – Diurnal (or other fundamental) temperature amplitudes at the shallow depth
z1and deeper depthz2. Units °C or K (identical for both).z1 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.z2 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.period (int, default
86_400) – Fundamental period P of the temperature wave in seconds.86 400s corresponds to a 24-hour diurnal cycle.
- Returns:
Thermal diffusivity α in m² s⁻¹.
- Return type:
- Raises:
ValueError – If
A1 <= A2(violates physical damping assumption) or ifz2 <= z1.
Notes
Amplitude extraction –
A1andA2should be obtained from a harmonic fit or spectral decomposition that isolates the target frequency; raw peak–trough differences are less robust.Logarithmic sensitivity – Because the formula involves
ln(A1/A2), small uncertainties in amplitudes propagate non-linearly; ensure adequate signal-to-noise ratio.Once
αis known, thermal conductivitykfollows fromk = ρc αgiven an independent estimate of volumetric heat capacityρc.
References
Carslaw, H. S., & Jaeger, J. C. (1959). Conduction of Heat in Solids (2nd ed., pp. 501–502). Oxford University Press.
Examples
>>> # Amplitudes from harmonic regression at 5 cm and 10 cm depths >>> alpha = thermal_diffusivity_amplitude( ... A1=6.3, A2=4.1, z1=0.05, z2=0.10 ... ) >>> print(f"α = {alpha:.2e} m² s⁻¹") α = 1.38e-07 m² s⁻¹
- soil_heat.thermal_diffusivity_lag(delta_t, z1, z2, period=86400)[source]
Estimate soil thermal diffusivity (α) from the phase lag of a temperature wave.
This method calculates thermal diffusivity based on the time it takes for a temperature wave (e.g., the diurnal cycle) to travel between two depths in the soil. The formula is derived from the analytical solution to the one-dimensional heat conduction equation for a periodic boundary condition:
\[\alpha = \frac{P (z_2 - z_1)^2}{4 \pi (\Delta t)^2}\]where \(P\) is the period of the wave, \((z_2 - z_1)\) is the distance between the sensors, and \(\Delta t\) is the time lag of the temperature peak between the two depths.
- Parameters:
delta_t (float or array_like) – Time lag between the temperature peaks at two depths (seconds).
z1 (float) – Depth of the upper sensor (m, positive downward).
z2 (float) – Depth of the lower sensor (m, positive downward).
period (int, optional) – The period of the temperature wave (seconds). The default is 86400, which corresponds to the daily (diurnal) cycle.
- Returns:
The calculated thermal diffusivity (α) in m² s⁻¹.
- Return type:
float or array_like
References
Nerpin, S.V., and Chudnovskii, A.F. (1967). Soil physics. (Moscow: Nauka) p 584. (In Russian)
- soil_heat.thermal_diffusivity_logrithmic(t1z1: float, t2z1: float, t3z1: float, t4z1: float, t1z2: float, t2z2: float, t3z2: float, t4z2: float, z1: float, z2: float, period: int = 86400) float[source]
Estimate soil thermal diffusivity (
α) between two depths using the Seemann four–temperature logarithmic method (also known as the Kolmogorov–Seemann method).The approach utilises two consecutive half‐period pairs of temperature measurements at a shallow depth
z1and a deeper depthz2. LetT₁–T₄denote the temperatures sampled at equal time steps (¼ P) apart, where P is the fundamental period of the harmonic forcing. The solution of the 1-D heat conduction equation for a sinusoidal boundary yields\[α \;=\; \frac{4 \, π \, (z_2 - z_1)^2} {P \;\bigl[\, \ln\bigl( ΔT_{z1} / ΔT_{z2} \bigr)\bigr]^2}\]with amplitude decrements
\[ΔT_{zij} = \sqrt{(T_1 - T_3)^2 + (T_2 - T_4)^2}\;.\]The formulation is advantageous when only a short record is available (four points suffice) but is sensitive to sensor noise and non-sinusoidal disturbances.
- Parameters:
t1z1 (float) – Temperatures (°C) at depth
z1sampled at four successive quarter-period intervals.t2z1 (float) – Temperatures (°C) at depth
z1sampled at four successive quarter-period intervals.t3z1 (float) – Temperatures (°C) at depth
z1sampled at four successive quarter-period intervals.t4z1 (float) – Temperatures (°C) at depth
z1sampled at four successive quarter-period intervals.t1z2 (float) – Temperatures (°C) at depth
z2sampled at the same times as the readings atz1.t2z2 (float) – Temperatures (°C) at depth
z2sampled at the same times as the readings atz1.t3z2 (float) – Temperatures (°C) at depth
z2sampled at the same times as the readings atz1.t4z2 (float) – Temperatures (°C) at depth
z2sampled at the same times as the readings atz1.z1 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1for a meaningful diffusivity.z2 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1for a meaningful diffusivity.period (int, default
86_400) – Fundamental period P of the temperature oscillation in seconds.86 400s corresponds to a 24-hour diurnal wave.
- Returns:
Thermal diffusivity α in m² s⁻¹.
- Return type:
Notes
Sampling interval – The four readings should be equidistant in time and span a full period P. A common practice is to use the peak, trough, and two mid-slope points of the diurnal cycle.
Noise sensitivity – Because the method involves logarithms of amplitude ratios, small errors in temperature can propagate strongly; consider pre-smoothing or repeating the calculation on multiple windows and averaging.
Relation to conductivity – Once
αis known, bulk thermal conductivitykfollows fromk = ρc αwith an independent estimate of volumetric heat capacityρc.
References
Kolmogorov, A. N. (1950). On the question of determining the coefficient of thermal diffusivity of the soil. Izvestiya Akademii Nauk SSSR, Ser. Geogr. Geofiz., 14 (2), 97–99. (In Russian)
Seemann, W. (1928). Die Wärmeleitung in der Bodenschicht. Springer, Berlin.
Examples
>>> α = thermal_diffusivity_logrithmic( ... 22.5, 20.3, 18.4, 20.1, # temps @ z1 ... 18.7, 17.2, 15.9, 17.1, # temps @ z2 ... z1=0.05, z2=0.10, ... ) >>> print(f"α = {α:.2e} m²/s") α = 1.46e-07 m²/s
- soil_heat.calc_thermal_diffusivity_log_pair(df, depth1_col, depth2_col, z1, z2, period=86400)[source]
Estimate soil thermal diffusivity (
α) between two depths using the four-point logarithmic amplitude method.The function extracts the first four consecutive samples from two temperature records—one at the shallow depth
z1and one at the deeper depthz2—and passes them tothermal_diffusivity_logrithmic(). That helper implements the log–ratio solution of the 1-D heat‐conduction equation for a sinusoidal boundary condition (Horton et al., 1934; de Vries, 1963):\[α = \frac{(z_2 - z_1)^2} {2P\;\ln\left(\frac{ΔT_{\!z1}}{ΔT_{\!z2}}\right)},\]where
P is the forcing period (s),
\(ΔT_{\!z}\) is the logarithmic temperature decrement derived from four successive measurements at depth z.
The approach is robust for short windows (four points suffice) but is sensitive to noise; it is best applied to periods with clear, smooth diurnal cycling.
- Parameters:
df (pandas.DataFrame) – Time‐indexed data containing at least the two temperature columns specified by
depth1_colanddepth2_col. Only the first four rows are used in the calculation.depth1_col (str) – Column names for the shallow (
z1) and deeper (z2) temperature series, respectively.depth2_col (str) – Column names for the shallow (
z1) and deeper (z2) temperature series, respectively.z1 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.z2 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.period (int, default
86_400) – Dominant temperature oscillation period P in seconds. The default (86 400 s) corresponds to 24 h.
- Returns:
Thermal diffusivity
αin m² s⁻¹. ReturnsNonewhen fewer than four valid samples are available or ifthermal_diffusivity_logrithmicitself returnsNone.- Return type:
float or None
- Warns:
UserWarning – Issued (via
print) when fewer than four rows are present in df, in which case the method is skipped andNoneis returned.
Notes
Data requirement – The function does not resample or align series; it simply grabs the first four rows. Pre-filter or sort your DataFrame accordingly.
Noise sensitivity – Because the method depends on small differences between successive temperature readings, apply a smoothing filter or select a high-signal period to minimise error.
Relationship to conductivity – Once
αis known, bulk thermal conductivitykcan be obtained fromk = ρc αgiven an estimate of volumetric heat capacityρc.
References
Horton, R., Wierenga, P. J., Nielsen, D. R., & de Vries, D. A. (1983). Calorimetric determination of soil thermal properties. Soil Science Society of America Journal, 47, 104–111.
de Vries, D. A. (1963). Thermal properties of soils. In Physics of Plant Environment (pp. 210–235). North-Holland.
Examples
>>> α_log = calc_thermal_diffusivity_log_pair( ... df=df.sort_index(), # ensure chronological order ... depth1_col='ts_05cm', ... depth2_col='ts_10cm', ... z1=0.05, z2=0.10, ... ) >>> if α_log is not None: ... print(f"Log-method α = {α_log:.2e} m² s⁻¹") Log-method α = 1.45e-07 m² s⁻¹
- soil_heat.calculate_thermal_diffusivity_for_pair(df, col1, col2, z1, z2, period=86400)[source]
Estimate soil thermal diffusivity (
α) between two depths using three classical harmonic methods: log-amplitude, amplitude ratio, and phase shift.Given two temperature time-series measured at depths
z1andz2, the function first extracts the dominant diurnal signal—its amplitude and phase—then applies the analytical solutions of the 1-D heat wave equation for a homogeneous medium subject to sinusoidal forcing (Carslaw & Jaeger, 1959).- 1. Log-Amplitude (α\_log)
Uses the decay of the harmonic amplitude with depth:
\[α\_{\text{log}} = \frac{(z_2 - z_1)^2} {2\,P\;\ln\bigl(A_1 / A_2\bigr)}\]
- 2. Amplitude Ratio (α\_amp)
Algebraically identical to the log-amplitude method but expressed directly in terms of the two amplitudes:
\[α\_{\text{amp}} = \frac{(z_2 - z_1)^2\;\omega} {2\,[\ln(A_1/A_2)]^2}\]where
ω = 2π / Pis the angular frequency.
- 3. Phase Lag (α\_lag)
Relates the travel time (phase shift) of the temperature wave:
\[α\_{\text{lag}} = \frac{(z_2 - z_1)^2}{2\,Δt\,P}\]with
Δtthe peak-to-peak time lag (s).
- Parameters:
df (pandas.DataFrame) – Time-indexed frame containing temperature observations.
col1 (str) – Column names for the shallow and deeper temperature series, respectively.
col2 (str) – Column names for the shallow and deeper temperature series, respectively.
z1 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.z2 (float) – Sensor depths in metres (positive downward). Must satisfy
z2 > z1.period (int, default
86_400) – Fundamental period P of the harmonic forcing in seconds.86 400s corresponds to 24 h diurnal cycling.
- Returns:
Mapping of method identifiers to diffusivity estimates (m² s⁻¹):
'alpha_log'– logarithmic amplitude method.'alpha_amp'– direct amplitude-ratio method.'alpha_lag'– phase-shift (lag) method.
Any method returning None inside intermediate helpers is propagated unchanged.
- Return type:
- Raises:
ValueError – If
z1≥z2or if either column is missing in df.
Notes
diurnal_amplitudeextracts the half range of the 24-h harmonic, typically via fast Fourier transform or STL decomposition.diurnal_peak_lagreturns the modal lag in hours; the value is internally converted to seconds.The function assumes a single dominant harmonic. Strong synoptic or weather-front variability can bias results; apply filtering or select periods with clear diurnal cycling.
Thermal diffusivity relates to thermal conductivity
kthrough\[k = ρ c \, α\]once bulk volumetric heat capacity
ρcis known.
References
Carslaw, H. S., & Jaeger, J. C. (1959). Conduction of Heat in Solids (2nd ed.). Oxford University Press.
Examples
>>> depth_map = {'ts_05cm': 0.05, 'ts_10cm': 0.10} >>> α = calculate_thermal_diffusivity_for_pair( ... df, 'ts_05cm', 'ts_10cm', ... z1=depth_map['ts_05cm'], z2=depth_map['ts_10cm']) >>> for meth, val in α.items(): ... print(f"{meth}: {val:.2e} m² s⁻¹") alpha_log: 1.43e-07 m² s⁻¹ alpha_amp: 1.41e-07 m² s⁻¹ alpha_lag: 1.38e-07 m² s⁻¹
- soil_heat.calculate_thermal_properties_for_all_pairs(df, depth_mapping, period=86400)[source]
Compute thermal diffusivity, thermal conductivity, and soil heat flux for every unique pair of temperature sensors in a profile.
The routine iterates over all combinations of the depth‐indexed temperature columns supplied in
depth_mapping. For each pair (z₁, z₂) itDerives thermal diffusivity
αwithcalculate_thermal_diffusivity_for_pair().Converts
αto thermal conductivitykviathermal_conductivity(), using the mean volumetric water- content of the two layers.Estimates instantaneous soil heat flux
Gby callingsoil_heat_flux().
Results are returned in a tidy, hierarchical
DataFramewhose outermost index encodes the depth pair (e.g.'0.05-0.10').- Parameters:
df (pandas.DataFrame) –
Time-indexed data frame containing at least
temperature columns listed in
depth_mapping; units °C, column names typically follow a pattern such as'ts_05cm'.matching soil-water-content columns; each temperature column
'<name>ts'must have a companion column'<name>swc'in percent. These are averaged and divided by 100 to obtain volumetric θ (m³ m⁻³).
depth_mapping (dict[str, float]) – Mapping of temperature column names to sensor depths in metres (positive downward), e.g.
{'ts_05cm': 0.05, 'ts_10cm': 0.10}.period (int, default
86_400) – Dominant period of the temperature wave (s).86_400s corresponds to 24 h and is appropriate for daily forcing.
- Returns:
Concatenated frame of thermal properties for every depth pair. The outer
Indexlevel is the stringf"{z1}-{z2}"and the inner index matches the datetime index ofdf(after dropping rows with any missing data). For each analysis “method” returned bycalculate_thermal_diffusivity_for_pair()(keys of its result dict) the following columns are present:αThermal diffusivity (m² s⁻¹) for that method.
kThermal conductivity (W m⁻¹ K⁻¹) derived from the same α.
GSoil heat flux (W m⁻²) between depths z₁ and z₂.
θ_vLayer-average volumetric water content (m³ m⁻³).
- Return type:
pandas.DataFrame
Notes
Alignment – Each pairwise calculation is performed on a copy of
dfafter dropping all rows with any missing values to ensure consistent sample support for derived quantities.Extensibility – Additional diffusivity algorithms can be integrated by returning extra key–value pairs from
calculate_thermal_diffusivity_for_pair(); they will be propagated automatically.Performance – The loop scales O(n²) with the number of depths. For large sensor arrays, filter the pairs of interest beforehand.
Examples
>>> depth_map = {'ts_05cm': 0.05, 'ts_10cm': 0.10, 'ts_20cm': 0.20} >>> props = calculate_thermal_properties_for_all_pairs(df, depth_map) >>> props.loc['0.05-0.10'][['alpha_phase', 'G_phase']].plot() >>> props.groupby(level=0)['k_amplitude'].median().unstack()
- soil_heat.estimate_rhoc_dry(alpha: pandas.Series, theta: pandas.Series, porosity: float = 0.4, k_dry: float = 0.25, k_sat: float = 1.5, rhoc_w: float = 4180000.0, dry_quantile: float = 0.1) float[source]
Estimate the volumetric heat capacity of dry soil (
ρ c_dry).This routine combines concurrent measurements of soil thermal diffusivity (
α) and volumetric water content (θ) with a simple two–end-member mixing model for thermal conductivity (λ) to back-calculate the volumetric heat capacity of the dry soil matrix. Only the driest records—defined by the lowerdry_quantileof the observed moisture distribution—are used in the final statistic so that the latent contribution of soil water is negligible.The underlying relationships are
\[ \begin{align}\begin{aligned}\begin{split}λ(θ) &= k_\text{dry} + \frac{θ}{φ} \,\bigl(k_\text{sat} - k_\text{dry}\bigr) \\\end{split}\\\begin{split}C_v &= \frac{λ(θ)}{α} \\\end{split}\\ρ\,c_\text{dry} &= \frac{C_v - θ\,ρ\,c_w}{1-θ}\,,\end{aligned}\end{align} \]where
λ is thermal conductivity (W m⁻¹ K⁻¹),
α is thermal diffusivity (m² s⁻¹),
C_v is volumetric heat capacity of the moist soil (J m⁻³ K⁻¹), and
φ is total porosity (m³ m⁻³).
- Parameters:
alpha (pandas.Series) – Soil thermal diffusivity α (m² s⁻¹), indexed identically to theta (usually a time-series).
theta (pandas.Series) – Volumetric water content θ (m³ m⁻³).
porosity (float, default
0.40) – Total soil porosity φ (saturated water content).k_dry (float, default
0.25) – Thermal conductivity of air-dry soil (W m⁻¹ K⁻¹).k_sat (float, default
1.50) – Thermal conductivity of saturated soil (W m⁻¹ K⁻¹).rhoc_w (float, default
4.18e6) – Volumetric heat capacity of liquid water (J m⁻³ K⁻¹, ≈ 4.18 MJ m⁻³ K⁻¹).dry_quantile (float, default
0.10) – Fraction of the lowest moisture observations to treat as “dry” when taking the median. For example,0.10selects the driest 10 % of the record.
- Returns:
Median volumetric heat capacity of the dry soil matrix (J m⁻³ K⁻¹).
- Return type:
Notes
Alignment — The two series are first inner-joined so only timestamps present in both are considered.
Robustness — Using the median of the driest subset avoids bias from residual soil moisture while damping the influence of occasional outliers.
The default conductivity bounds
k_dry/k_satfollow typical literature values for mineral soils; adjust them for peat, organic, or highly gravelly substrates.
Examples
>>> rhoc_dry = estimate_rhoc_dry( ... alpha=df['alpha_10cm'], ... theta=df['VWC_10cm'], ... porosity=0.43, ... ) >>> print(f"ρ c_dry ≈ {rhoc_dry/1e6:.2f} MJ m⁻³ K⁻¹") 2.07 MJ m⁻³ K⁻¹
- soil_heat.calculate_soil_heat_storage(df, depths, porosity=0.45, bulk_density=1.3)[source]
Calculate soil heat storage and heat flux.
Parameters:
- dfDataFrame
Must contain columns: ‘T_5cm’, ‘T_10cm’, etc. and ‘SM_5cm’, ‘SM_10cm’, etc. where SM is volumetric soil moisture (0-1 or 0-100%)
- depthslist
Measurement depths in cm, e.g., [5, 10, 20, 30, 40, 50]
- porosityfloat
Soil porosity (0-1), default 0.45
- bulk_densityfloat
Dry soil bulk density (g/cm³), default 1.3
Returns:
DataFrame with heat storage and flux values
- soil_heat.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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,)
- soil_heat.reference_ground_heat_flux(temp_profile: numpy.ndarray, depths: Sequence[float], times: Sequence[float], cv: float, thermal_conductivity: float, gradient_depth: float = 0.2) numpy.ndarray[source]
Compute the reference ground‑heat flux G₀,M using the gradient method combined with calorimetry for heat storage (Liebethal & Foken 2006, Eq. 1).
The method combines the conductive heat flux at a reference depth with the rate of change of heat stored in the soil layer above that depth.
\[G_{0,M}(t) = -\lambda \frac{\partial T}{\partial z} \bigg|_{z=0.2m} + \int_{z=0}^{0.2m} c_v \frac{\partial T}{\partial t} dz\]- Parameters:
temp_profile (numpy.ndarray) – A 2D array of soil temperatures (°C or K) with shape (n_depths, n_times).
depths (Sequence[float]) – A sequence of measurement depths in meters (positive downward), corresponding to the rows of temp_profile.
times (Sequence[float]) – A sequence of time stamps in seconds (e.g., Unix timestamps), corresponding to the columns of temp_profile.
cv (float) – Volumetric heat capacity of the soil (J m⁻³ K⁻¹). Assumed to be constant with depth.
thermal_conductivity (float) – Soil thermal conductivity λ (W m⁻¹ K⁻¹). Assumed to be constant with depth.
gradient_depth (float, optional) – The depth (m) at which the vertical temperature gradient is evaluated, by default 0.20 m.
- Returns:
A 1D array of the instantaneous ground‑heat flux G₀,M (W m⁻²) at the surface for each time step. Positive values indicate downward flux.
- Return type:
numpy.ndarray
- Raises:
ValueError – If temp_profile shape does not match the lengths of depths and times.
- soil_heat.ground_heat_flux_pr(qs: numpy.ndarray, p: float) numpy.ndarray[source]
Estimate ground heat flux as a fixed fraction of net radiation (Liebethal & Foken 2006, Eq. 2).
This is a simple empirical relationship where the ground heat flux is assumed to be a constant proportion of the net radiation at the surface.
\[G_{0,PR}(t) = -p \cdot Q^*_s(t)\]- Parameters:
qs (numpy.ndarray) – Time series of net radiation (W m⁻²). Positive values are typically downward.
p (float) – The fraction of net radiation that is partitioned into ground heat flux (dimensionless, typically 0 < p < 1).
- Returns:
Time series of the estimated ground heat flux G₀,PR (W m⁻²).
- Return type:
numpy.ndarray
- soil_heat.ground_heat_flux_lr(qs: numpy.ndarray, a: float, b: float, lag_steps: int = 0) numpy.ndarray[source]
Estimate ground heat flux using a linear regression against net radiation, with an optional time lag (Liebethal & Foken 2006, Eq. 3).
\[G_{0,LR}(t) = a \cdot Q^*_s(t + \Delta t_G) + b\]- Parameters:
qs (numpy.ndarray) – Time series of net radiation (W m⁻²).
a (float) – The slope of the linear regression (dimensionless).
b (float) – The intercept of the linear regression (W m⁻²).
lag_steps (int, optional) – The integer time lag (number of array elements) to apply to the net radiation series. A positive value advances the series (e.g., qs[t+lag] is used for G[t]), by default 0.
- Returns:
Time series of the estimated ground heat flux G₀,LR (W m⁻²).
- Return type:
numpy.ndarray
- soil_heat.ur_coefficients(delta_ts: float | numpy.ndarray) Tuple[numpy.ndarray, numpy.ndarray][source]
Compute the parameters A and B for the universal net-radiation parameterization, based on the diurnal surface temperature amplitude (Liebethal & Foken 2006, Eq. 5 & 6).
\[A = 0.0074 \cdot \Delta T_s + 0.088 B = 1729 \cdot \Delta T_s + 65013\]- Parameters:
delta_ts (float or numpy.ndarray) – The diurnal amplitude of the surface temperature (K or °C).
- Returns:
A tuple containing the computed parameters (A, B). - A is dimensionless. - B is in seconds.
- Return type:
Tuple[numpy.ndarray, numpy.ndarray]
- soil_heat.ground_heat_flux_ur(qs: numpy.ndarray, times_sec: numpy.ndarray, delta_ts: float) numpy.ndarray[source]
Estimate ground heat flux using the universal net-radiation parameterization by Santanello & Friedl (2003), as cited in Liebethal & Foken (2006, Eq. 4).
This method modulates the fraction of net radiation partitioned to ground heat flux with a cosine function of the time of day.
\[G_{0,UR}(t) = -A \cdot \cos\left(\frac{2\pi (t + 10800)}{B}\right) \cdot Q^*_s(t)\]- Parameters:
qs (numpy.ndarray) – Time series of net radiation (W m⁻²).
times_sec (numpy.ndarray) – Time stamps in seconds relative to solar noon. Positive values indicate the afternoon.
delta_ts (float) – The diurnal amplitude of the surface temperature (K or °C).
- Returns:
Time series of the estimated ground heat flux G₀,UR (W m⁻²).
- Return type:
numpy.ndarray
- soil_heat.surface_temp_amplitude(delta_t1: float, delta_t2: float, z1: float, z2: float) float[source]
Estimate the diurnal surface-temperature amplitude (ΔT_s) from temperature amplitudes measured at two different soil depths (Liebethal & Foken 2006, Eq. 8).
This method assumes an exponential decay of the temperature wave amplitude with depth.
\[\Delta T_s = \Delta T_1 + \Delta T_2 \cdot \exp\left(\frac{z_2}{z_2 - z_1}\right)\]- Parameters:
- Returns:
The estimated diurnal surface-temperature amplitude ΔT_s (K or °C).
- Return type:
- Raises:
ValueError – If z2 is not greater than z1.
- soil_heat.phi_from_soil_moisture(theta_0_10: float, a_phi: float = 9.62, b_phi: float = 0.402) float[source]
Calculate the empirical parameter φ based on soil moisture content (Liebethal & Foken 2006, Eq. 10).
\[\phi = a_\phi \cdot \theta_{0-10} + b_\phi\]
- soil_heat.ground_heat_flux_sh(h: numpy.ndarray, phase_g0: Sequence[float], phase_h: Sequence[float], u_mean: float, phi: float, omega: float = 2 * np.pi / 86400.0) numpy.ndarray[source]
Estimate ground heat flux from the sensible heat flux (H) (Liebethal & Foken 2006, Eq. 9).
This method relates the ground heat flux to the sensible heat flux through a phase-shifted and scaled relationship.
\[G_{0,SH}(t) = -\frac{\phi}{\sqrt{\bar{u}}} \frac{\cos(\omega t + \varphi(G_0))} {\cos(\omega t + \varphi(H))} H(t)\]- Parameters:
h (numpy.ndarray) – Time series of sensible heat flux (W m⁻²).
phase_g0 (Sequence[float]) – Phase lags of the ground heat flux, φ(G₀), in radians. Must have the same length as h.
phase_h (Sequence[float]) – Phase lags of the sensible heat flux, φ(H), in radians. Must have the same length as h.
u_mean (float) – Mean horizontal wind speed during the daytime (m s⁻¹).
phi (float) – An empirical dimensionless parameter, often derived from soil moisture via phi_from_soil_moisture.
omega (float, optional) – The diurnal angular frequency (rad s⁻¹), by default 2π/86400.
- Returns:
Time series of the estimated ground heat flux G₀,SH (W m⁻²).
- Return type:
numpy.ndarray
- Raises:
ValueError – If the length of phase arrays does not match the length of h.
- soil_heat.ground_heat_flux_sm(gp: numpy.ndarray, t1: numpy.ndarray, delta_t: numpy.ndarray, cv: float, zp: float, dt_seconds: float) numpy.ndarray[source]
Estimate surface ground heat flux using the “simple measurement” parameterization, correcting a heat flux plate measurement with a storage term (Liebethal & Foken 2006, Eq. 11).
\[G_{0,SM}(t) = G_p(t) + c_v z_p \left( \frac{dT_1}{dt} + \frac{1}{2} \frac{d(\Delta T)}{dt} \right)\]- Parameters:
gp (numpy.ndarray) – Heat flux plate measurements at depth zp (W m⁻²).
t1 (numpy.ndarray) – Soil temperature at 0.01 m depth (K or °C).
delta_t (numpy.ndarray) – Temperature difference T(0.01 m) - T(zp) (K).
cv (float) – Volumetric heat capacity of the soil (J m⁻³ K⁻¹).
zp (float) – Depth of the heat flux plate (m, positive downward).
dt_seconds (float) – The constant time step between consecutive samples (s).
- Returns:
Time series of the estimated ground heat flux G₀,SM (W m⁻²). The first element will be NaN due to the backward difference.
- Return type:
numpy.ndarray
- soil_heat.active_layer_thickness(lambda_: float, cv: float, omega: float = 2 * np.pi / 86400) float[source]
Calculate the thickness of the active soil layer (δz) for the force-restore method (Liebethal & Foken 2006, Eq. 13).
\[\delta_z = \sqrt{\frac{\lambda}{2 c_v \omega}}\]
- soil_heat.ground_heat_flux_fr(tg: numpy.ndarray, tg_avg: float, cv: float, lambda_: float, delta_z: float | None = None, times: numpy.ndarray | None = None) numpy.ndarray[source]
Estimate ground heat flux using the two-layer force-restore method (Liebethal & Foken 2006, Eq. 12).
\[G_{0,FR}(t) = -\delta_z c_v \frac{dT_g}{dt} + \sqrt{\lambda \omega c_v} \left( \frac{1}{\omega}\frac{dT_g}{dt} + (T_g - \bar{T_g}) \right)\]- Parameters:
tg (numpy.ndarray) – Time series of the upper (surface) layer temperature, Tg(t) (K).
tg_avg (float) – The long-term average or “restoring” temperature, T̄g (K).
cv (float) – Volumetric heat capacity of the soil (J m⁻³ K⁻¹).
lambda (float) – Soil thermal conductivity λ (W m⁻¹ K⁻¹).
delta_z (float, optional) – Thickness of the active soil layer δz (m). If None, it is calculated internally using active_layer_thickness, by default None.
times (numpy.ndarray, optional) – Time stamps in seconds corresponding to tg. Required for calculating the time derivative. If None, assumes a uniform time step of 1 second, by default None.
- Returns:
Time series of the estimated ground heat flux G₀,FR (W m⁻²).
- Return type:
numpy.ndarray
- soil_heat.energy_balance_residual(Rn: float | numpy.ndarray, H: float | numpy.ndarray, LE: float | numpy.ndarray, G0: float | numpy.ndarray) float | numpy.ndarray[source]
Compute the closure residual of the surface energy balance (SEB).
The classical SEB for land–atmosphere exchange is
\[R_n - G_0 \;=\; H + LE +\]arepsilon ,
where the residual term :math:`
- arepsilon` quantifies the lack of
closure. Rearranging gives
\[\]
arepsilon ;=; R_n - G_0 - H - LE ,
which is what this helper returns.
- Rnfloat or array_like
Net radiation Rₙ (W m⁻²). Positive downward.
- Hfloat or array_like
Sensible heat flux H (W m⁻²). Positive upward (atmosphere ← surface).
- LEfloat or array_like
Latent heat flux LE (W m⁻²). Positive upward.
- G0float 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.
- float or ndarray
Energy‐balance residual :math:`
- arepsilon` (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.
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.
>>> ε = 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")
- soil_heat.ground_heat_flux_conventional(k: float, dT_dz_at_zr: float, rho_c: float, dT_dt_profile: Sequence[float], z_profile: Sequence[float]) float[source]
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 withz_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 – Ground heat flux at the surface (W m-2). Positive into the soil.
- Return type:
- soil_heat.green_function_temperature(z: float, t: float, kappa: float) float[source]
Green-function solution \(g_z(t)\) for the one‐dimensional, semi-infinite heat equation with a unit surface‐flux impulse at \(t=0,\,z=0\).
The governing partial differential equation is
\[\frac{\partial T}{\partial t} \;=\; \kappa\,\frac{\partial^{2} T}{\partial z^{2}}, \qquad z \ge 0,\; t > 0,\]with initial condition \(T(z,0)=0\) and boundary condition corresponding to a Dirac δ–heat pulse applied at the surface (\(z=0\)). The resulting Green function (Carslaw & Jaeger, 1959, Eq. 7) is
\[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 \(g_z(t)=0\) for \(t\le 0\) (causality).
- Parameters:
- Returns:
Green-function value \(g_z(t)\) (units m, because the solution integrates heat-flux density with respect to depth to yield temperature).
- Return type:
Notes
Causality check – If
t≤ 0 the function short-circuits and returns 0.0.Vectorisation – For vector or array input use
numpy.vectorize()or wrap the function innumpy.frompyfunc(); the core implementation is scalar for numerical clarity.Usage –
g_z(t)can be convolved with an arbitrary surface heat-flux time seriesq₀(t)to obtain temperature at depth via\[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
- soil_heat.temperature_convolution_solution(z: float, t_series: numpy.ndarray, f_series: numpy.ndarray, kappa: float, Ti: float = 0.0) numpy.ndarray[source]
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).
- soil_heat.soil_heat_flux_from_G0(z: float, t_series: numpy.ndarray, G0_series: numpy.ndarray, kappa: float) numpy.ndarray[source]
Compute G(z,t) from a known surface flux series G0 (Eq. 9).
- soil_heat.estimate_G0_from_Gz(Gz_series: numpy.ndarray, z_r: float, kappa: float, dt: float) numpy.ndarray[source]
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.
- soil_heat.sinusoidal_boundary_flux(t: float | numpy.ndarray, A: float, omega: float, epsilon: float) float | numpy.ndarray[source]
Evaluate a sinusoidal surface heat-flux forcing
\[q_0(t) \;=\; A \,\sin\!igl(\omega t +\]arepsilonigr),
which is commonly used as a boundary condition for analytical soil- heat-flux solutions (see Eq. 13 of the companion text).
- tfloat 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.- Afloat
Amplitude of the surface heat flux (W m⁻²). Positive downward into the soil.
- omegafloat
Angular frequency (rad s⁻¹). For a diurnal cycle
omega = 2 π / 86 400≈ 7.272 × 10⁻⁵ rad s⁻¹.- epsilonfloat
Phase shift ε (rad). Positive values delay the flux peak, negative values advance it.
- ndarray or float
Instantaneous surface heat flux q₀(t) (W m⁻²) with shape given by NumPy broadcasting of the inputs.
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.sinand therefore fully vectorised.Period — The period P (s) of the forcing is related to
omegaby P = 2 π / ω.
>>> 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()
- soil_heat.soil_temperature_sinusoidal(z: float, t: float | numpy.ndarray, A: float, omega: float, epsilon: float, Ti: float, kappa: float) float | numpy.ndarray[source]
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
\[q_0(t) \;=\; A \sin(\omega t +\]arepsilon)
is (Carslaw & Jaeger 1959, Eq. 14)
\[T(z,t) \;=\; T_i \;+\; \underbrace{\frac{A}{\kappa\sqrt{\omega}} e^{-z r}\, \sin\!igl(\omega t +\]- arepsilon - z r - π/4bigr)}
_{text{steady harmonic}} ;+; underbrace{T_{text{trans}}(z,t)} _{text{transient integral}} ,
where \(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.
- zfloat
Depth below the soil surface (m, positive downward).
- tfloat or array_like
Time since the start of the forcing (s). Scalar or NumPy array.
- Afloat
Amplitude of the sinusoidal surface heat flux (W m⁻², positive downward). Consistent with
sinusoidal_boundary_flux().- omegafloat
Angular frequency of the forcing (rad s⁻¹). For a diurnal wave
omega = 2 * π / 86_400.- epsilonfloat
Phase shift ε (rad) of the surface heat-flux wave.
- Tifloat
Initial uniform soil temperature Tᵢ (°C or K).
- kappafloat
Thermal diffusivity κ of the soil (m² s⁻¹).
- 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.
Steady component – The exponential term
exp(-z*r)dampens amplitude with depth while the phase lag isz*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.
>>> 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()
- soil_heat.soil_heat_flux_sinusoidal(z: float, t: float | numpy.ndarray, A: float, omega: float, epsilon: float, kappa: float) float | numpy.ndarray[source]
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
\[q_0(t) \;=\; A \,\sin\!igl(\omega t +\]arepsilonigr)
admits the Green-function solution (Carslaw & Jaeger, 1959, Eq. 15)
\[G(z,t) \;=\; A\,e^{-z r}\,\sin(\omega t +\]- arepsilon - z r)
;+; G_{ ext{trans}}(z,t),
where \(r = \sqrt{ \omega / 2\kappa }\) and the transient term
\[G_{ ext{trans}}(z,t) \;=\; -\]- rac{2 A kappa}{pi}
int_{0}^{infty} frac{( kappa x^{2}sin
arepsilon - omega cos arepsilon )
;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).
- zfloat
Depth below the soil surface (m, positive downward).
- tfloat or array_like
Time since the onset of forcing (s). Scalar or NumPy array.
- Afloat
Amplitude of the surface heat flux (W m⁻², positive downward).
- omegafloat
Angular frequency ω of the forcing (rad s⁻¹). For a diurnal wave
omega = 2 π / 86 400.- epsilonfloat
Phase shift ε of the surface flux (radians).
- kappafloat
Thermal diffusivity κ of the soil (m² s⁻¹).
- float or ndarray
Heat flux G(z,t) at depth z (W m⁻²) with shape equal to
tafter NumPy broadcasting. Positive values downward, matching the sign convention of A.
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
soil_temperature_sinusoidal()gives T(z,t) under the same boundary forcing; G and T satisfy Fourier’s lawG = -k ∂T/∂zonce thermal conductivity k is specified.
Carslaw, H. S., & Jaeger, J. C. (1959). Conduction of Heat in Solids (2nd ed., pp. 100–102). Oxford University Press.
>>> 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()
- soil_heat.heat_capacity_moist_soil(theta_v: float | numpy.ndarray, theta_s: float, rho_c_s: float = 1260000.0, rho_c_w: float = 4200000.0) float | numpy.ndarray[source]
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 – Heat capacity of dry soil / water (J m-3 K-1).
rho_c_w – Heat capacity of dry soil / water (J m-3 K-1).
- soil_heat.pf_from_theta(theta_v: float | numpy.ndarray, theta_s: float, psi_s: float, b: float) float | numpy.ndarray[source]
Return Pf (Eq. 18) from volumetric water content.
- soil_heat.thermal_conductivity_moist_soil(theta_v: float | numpy.ndarray, theta_s: float, psi_s: float, b: float) float | numpy.ndarray[source]
Thermal conductivity parameterisation, Eq. (17).
- soil_heat.thermal_diffusivity(k: float | numpy.ndarray, rho_c: float | numpy.ndarray) float | numpy.ndarray[source]
Return κ = k / (ρ c).
- soil_heat.soil_heat_flux(Tz: ArrayLike, dz: ArrayLike, lambda_s: ArrayLike | float) numpy.ndarray[source]
Compute heat flux G at cell interfaces using Fourier’s law.
This function calculates the conductive heat flux between soil layers based on the temperature gradient and thermal conductivity.
\[G = -\lambda_s \frac{\partial T}{\partial z}\]- Parameters:
Tz (ArrayLike) – A 1D array of temperatures at the centre of each soil layer (K).
dz (ArrayLike) – A 1D array of the thickness of each soil layer (m).
lambda_s (ArrayLike or float) – Thermal conductivity for each layer (W m⁻¹ K⁻¹). Can be a single value (for homogeneous soil) or an array matching Tz.
- Returns:
An array of heat fluxes (W m⁻²) at the interfaces between layers, including the surface and bottom boundaries. The length of the output is len(Tz) + 1. Positive values indicate downward flux.
- Return type:
numpy.ndarray
- soil_heat.integrated_soil_heat_flux(rho_c: ArrayLike, T_before: ArrayLike, T_after: ArrayLike, dz: ArrayLike, dt: float, G_ref: float = 0.0) numpy.ndarray[source]
Calculate the soil heat flux profile by integrating the change in heat storage upwards from a reference depth (Yang & Wang 2008, Eq. 5).
\[G(z_i) = G(z_{ref}) + \int_{z_i}^{z_{ref}} \rho_s c_s \frac{\partial T}{\partial t} dz\]- Parameters:
rho_c (ArrayLike) – A 1D array of volumetric heat capacity (ρ_s c_s) for each soil layer (J m⁻³ K⁻¹).
T_before (ArrayLike) – A 1D array of temperatures at the start of the time step (K).
T_after (ArrayLike) – A 1D array of temperatures at the end of the time step (K).
dz (ArrayLike) – A 1D array of layer thicknesses (m).
dt (float) – The duration of the time step (s).
G_ref (float, optional) – The heat flux at the lower reference depth z_ref (W m⁻²). This is typically assumed to be zero at a sufficient depth, by default 0.0.
- Returns:
An array of heat fluxes (W m⁻²) at the upper interface of each layer, with a size of len(dz).
- Return type:
numpy.ndarray
- soil_heat.volumetric_heat_capacity(theta: ArrayLike, theta_sat: float | ArrayLike) numpy.ndarray[source]
Calculate the volumetric heat capacity of moist soil based on its water content (Yang & Wang 2008, Eq. 4).
\[\rho_s c_s = (1 - \theta_{sat}) \rho_{d}c_{d} + \theta \rho_w c_w\]- Parameters:
theta (ArrayLike) – Volumetric water content (m³ m⁻³).
theta_sat (float or ArrayLike) – Soil porosity, i.e., saturated volumetric water content (m³ m⁻³).
- Returns:
The volumetric heat capacity ρ_s c_s (J m⁻³ K⁻¹).
- Return type:
numpy.ndarray
- soil_heat.stretched_grid(n: int, D: float, xi: float) numpy.ndarray[source]
Generate layer thicknesses for a vertically stretched grid.
The grid layer thickness increases exponentially with depth, allowing for higher resolution near the surface.
\[\Delta z_i = \Delta z_0 \exp(\xi (i-1))\]- Parameters:
- Returns:
A 1D array of thickness Δz_i for each of the n layers (m).
- Return type:
numpy.ndarray
- soil_heat.solve_tde(T_prev: ArrayLike, dz: ArrayLike, rho_c: ArrayLike, lambda_s: ArrayLike | float, Tsfc: float, Tbot: float, dt: float) numpy.ndarray[source]
Solve the 1D thermal diffusion equation for one time step using an implicit Crank-Nicolson scheme (Yang & Wang 2008, Eq. 7).
This function sets up and solves the tridiagonal system of linear equations M * T_new = D to find the temperature profile at the next time step.
- Parameters:
T_prev (ArrayLike) – Temperature profile at the previous time step (K).
dz (ArrayLike) – Layer thicknesses (m).
rho_c (ArrayLike) – Volumetric heat capacity of each layer (J m⁻³ K⁻¹).
lambda_s (ArrayLike or float) – Thermal conductivity of each layer (W m⁻¹ K⁻¹).
Tsfc (float) – Surface temperature boundary condition (K).
Tbot (float) – Bottom temperature boundary condition (K).
dt (float) – Time step (s).
- Returns:
The new temperature profile at t + dt, including boundary nodes.
- Return type:
numpy.ndarray
- soil_heat.correct_profile(T_model: ArrayLike, depths_model: ArrayLike, T_obs: ArrayLike, depths_obs: ArrayLike) numpy.ndarray[source]
Correct a modeled temperature profile using observed temperatures.
This function calculates the bias (error) between the model and observations at the observation depths, then linearly interpolates this bias across the entire model grid to correct the profile.
- Parameters:
T_model (ArrayLike) – The modeled temperature profile (K).
depths_model (ArrayLike) – The depths corresponding to T_model (m).
T_obs (ArrayLike) – The observed temperatures (K).
depths_obs (ArrayLike) – The depths of the observations (m).
- Returns:
The corrected temperature profile.
- Return type:
numpy.ndarray
- soil_heat.surface_temperature_longwave(R_lw_up: float, R_lw_dn: float, emissivity: float = 0.98) float[source]
Calculate surface temperature from upward and downward long-wave radiation measurements using the Stefan-Boltzmann law (Yang & Wang 2008, Eq. 8).
\[T_s = \left[ \frac{R_{lw}^{\uparrow} - (1 - \epsilon) R_{lw}^{\downarrow}} {\epsilon \sigma} \right]^{1/4}\]
- soil_heat.thermal_conductivity_yang2008(theta: ArrayLike, theta_sat: float, rho_dry: float | ArrayLike) numpy.ndarray[source]
Estimate soil thermal conductivity based on soil moisture and dry density, following Yang et al. (2005) as cited in Yang & Wang (2008, Eq. 9).
- soil_heat.flux_error_linear(rho_c: ArrayLike, S2_minus_S1: ArrayLike, dt: float) numpy.ndarray[source]
Calculate the diagnostic error in heat flux that arises from assuming a linear temperature profile between measurement points (Yang & Wang 2008, Eq. 11).
\[\Delta G_i = \frac{\rho_s c_s (S_{i,2} - S_{i,1})}{\Delta t}\]- Parameters:
rho_c (ArrayLike) – Volumetric heat capacity for each layer (J m⁻³ K⁻¹).
S2_minus_S1 (ArrayLike) – The area difference representing the deviation from linearity.
dt (float) – The time step (s).
- Returns:
The flux error for each layer (W m⁻²).
- Return type:
numpy.ndarray
- soil_heat.surface_energy_residual(R_net: float, H: float, LE: float, G0: float) float[source]
Calculate the residual of the surface energy budget (Yang & Wang 2008, Eq. 12).
\[\Delta E = R_{net} - (H + LE + G_0)\]
- soil_heat.tdec_step(T_prev: ArrayLike, dz: ArrayLike, theta: ArrayLike, theta_sat: float, rho_dry: float, lambda_const: float, Tsfc: float, Tbot: float, dt: float, depths_model: ArrayLike, T_obs: ArrayLike, depths_obs: ArrayLike) Tuple[numpy.ndarray, numpy.ndarray][source]
Perform a single integration step of the TDEC (Temperature Diffusion Error Correction) scheme.
This function encapsulates the predict-correct sequence for one time step.
- Parameters:
T_prev (ArrayLike) – Temperature profile from the previous time step (K).
dz (ArrayLike) – Layer thicknesses (m).
theta (ArrayLike) – Volumetric water content profile (m³ m⁻³).
theta_sat (float) – Soil porosity (m³ m⁻³).
rho_dry (float) – Dry soil bulk density (kg m⁻³).
lambda_const (float) – A constant thermal conductivity for the prediction step (W m⁻¹ K⁻¹).
Tsfc (float) – Surface temperature boundary condition (K).
Tbot (float) – Bottom temperature boundary condition (K).
dt (float) – Time step (s).
depths_model (ArrayLike) – Depths of the model grid nodes (m).
T_obs (ArrayLike) – Observed temperatures for the correction step (K).
depths_obs (ArrayLike) – Depths of the observed temperatures (m).
- Returns:
A tuple containing: - T_corr: The corrected temperature profile at t + dt. - G_prof: The corresponding heat flux profile (W m⁻²) at layer interfaces.
- Return type:
Tuple[np.ndarray, np.ndarray]
- soil_heat.soil_heat_flux_general(t_current: float | numpy.typing.ArrayLike, t_previous: float | numpy.typing.ArrayLike, delta_t: float, delta_z: float = 1.0, c_s: float = DEFAULT_CS) float | numpy.ndarray[source]
Soil heat flux from the general calorimetric equation (FAO-56 Eq. 41).
\[G = c_s \, \frac{T_i - T_{i-1}}{\Delta t} \, \Delta z\]- Parameters:
t_current (float or array_like) – Mean air temperature of the current period [°C].
t_previous (float or array_like) – Mean air temperature of the previous period [°C].
delta_t (float) – Length of the time interval in days. For successive months use ~30.4; for a single month pair use the actual number of days between the midpoints of the two months.
delta_z (float, optional) – Effective soil depth [m]. Typical values are 0.10–0.20 m for daily periods and 0.5–2.0 m for monthly periods. Default is 1.0 m (standard for successive-month calculations).
c_s (float, optional) – Soil heat capacity [MJ m⁻³ °C⁻¹]. Default
2.1.
- Returns:
Soil heat flux G [MJ m⁻² day⁻¹]. Positive values indicate heat transfer into the soil (warming); negative values indicate heat release from the soil (cooling).
- Return type:
float or numpy.ndarray
Notes
This is the most general form. Equations 43 and 44 are simplified versions derived from this equation by assuming c_s = 2.1 MJ/(m³·°C), Δz = 1.0 m, and Δt appropriate for monthly intervals (~30 days).
For Eq. 43: 0.07 ≈ 2.1 × 1.0 / (2 × 30) For Eq. 44: 0.14 ≈ 2.1 × 1.0 / (1 × 30) (÷ by half the interval)
Examples
FAO-56 Example 13 — March to April:
>>> soil_heat_flux_general(16.1, 14.1, delta_t=30.0, delta_z=1.0) 0.14
- soil_heat.soil_heat_flux_daily() float[source]
Daily soil heat flux for a grass reference surface (FAO-56 Eq. 42).
For a daily time step, soil heat flux beneath a grass reference surface is small relative to net radiation and may be assumed to be zero.
- Returns:
0.0[MJ m⁻² day⁻¹].- Return type:
Notes
This assumption holds for a dense, grass-covered surface where the soil surface temperature cycle nearly cancels over 24 hours. For bare soil or sparse vegetation, daily G can be significant and should be estimated by other means (e.g., the general equation or measured data).
- soil_heat.soil_heat_flux_monthly(t_month_prev: float | numpy.typing.ArrayLike, t_month_next: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Monthly soil heat flux using adjacent months (FAO-56 Eq. 43).
When both the previous and next month’s mean air temperatures are available, this is the preferred monthly estimate.
\[G_{\text{month},i} = 0.07 \, (T_{i+1} - T_{i-1})\]- Parameters:
- Returns:
Monthly soil heat flux [MJ m⁻² day⁻¹].
- Return type:
float or numpy.ndarray
Notes
The coefficient 0.07 derives from Eq. 41 with c_s = 2.1 MJ/(m³·°C), Δz = 1.0 m, and a two-month span (Δt ≈ 2 × 30 = 60 days):
2.1 × 1.0 / 60 = 0.035 → rounded/adopted as 0.07 in FAO-56.
(The FAO text notes that the coefficient accounts for the relationship between air and soil temperature damping; it is empirically calibrated for a grass reference surface.)
Examples
FAO-56 Example 13 — April soil heat flux (March = 14.1 °C, May = 18.8 °C):
>>> soil_heat_flux_monthly(14.1, 18.8) 0.329
- soil_heat.soil_heat_flux_monthly_prev_only(t_month_prev: float | numpy.typing.ArrayLike, t_month_cur: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Monthly soil heat flux — next month unavailable (FAO-56 Eq. 44).
When the next month’s temperature is not yet known (e.g., real-time or forecast applications), use this equation with the current and previous month.
\[G_{\text{month},i} = 0.14 \, (T_i - T_{i-1})\]- Parameters:
- Returns:
Monthly soil heat flux [MJ m⁻² day⁻¹].
- Return type:
float or numpy.ndarray
Notes
The coefficient 0.14 is exactly twice the Eq. 43 coefficient because the temperature difference spans only one month instead of two:
2.1 × 1.0 / 30 ≈ 0.07 × 2 = 0.14
When both adjacent months are available, prefer
soil_heat_flux_monthly()(Eq. 43) for greater accuracy.Examples
FAO-56 Example 13 — March soil heat flux (Feb = 12.1 °C, Mar = 14.1 °C):
>>> soil_heat_flux_monthly_prev_only(12.1, 14.1) 0.28
- soil_heat.soil_heat_flux_hourly(rn: float | numpy.typing.ArrayLike, is_daytime: bool | numpy.typing.ArrayLike, day_coeff: float = 0.1, night_coeff: float = 0.5) float | numpy.ndarray[source]
Hourly (sub-daily) soil heat flux as a fraction of net radiation.
Based on the ASCE Standardized Reference Evapotranspiration Equation (ASCE-EWRI, 2005), which recommends estimating G for hourly or shorter periods as a proportion of net radiation Rn:
\[ \begin{align}\begin{aligned}G_{\text{day}} &= 0.1 \; R_n \quad \text{(daytime)}\\G_{\text{night}} &= 0.5 \; R_n \quad \text{(nighttime)}\end{aligned}\end{align} \]- Parameters:
rn (float or array_like) – Net radiation at the crop surface [MJ m⁻² hr⁻¹].
is_daytime (bool or array_like of bool) –
Truefor daytime periods,Falsefor nighttime. For array inputs, must be broadcastable with rn.day_coeff (float, optional) – Fraction of Rn used for daytime G. Default
0.1.night_coeff (float, optional) – Fraction of Rn used for nighttime G. Default
0.5.
- Returns:
Soil heat flux [MJ m⁻² hr⁻¹] (same units as rn).
- Return type:
float or numpy.ndarray
Notes
Daytime is typically defined as the period when Rn > 0. Some implementations use the solar angle or the time relative to sunrise and sunset instead.
The default coefficients (0.1 / 0.5) are calibrated for a grass reference surface. For other surfaces the ratio G/Rn can vary considerably (see Allen et al., 1998, Table 7-2).
Examples
Daytime hour with Rn = 2.5 MJ m⁻² hr⁻¹:
>>> soil_heat_flux_hourly(2.5, is_daytime=True) 0.25
Nighttime hour with Rn = -0.4 MJ m⁻² hr⁻¹:
>>> soil_heat_flux_hourly(-0.4, is_daytime=False) -0.2
- soil_heat.soil_heat_flux_hourly_auto(rn: float | numpy.typing.ArrayLike, day_coeff: float = 0.1, night_coeff: float = 0.5) float | numpy.ndarray[source]
Hourly soil heat flux, auto-detecting day/night from Rn sign.
Convenience wrapper around
soil_heat_flux_hourly()that usesRn > 0as the daytime indicator, which is the most common approach in practice.- Parameters:
- Returns:
Soil heat flux [MJ m⁻² hr⁻¹].
- Return type:
float or numpy.ndarray
Examples
>>> soil_heat_flux_hourly_auto(2.5) 0.25 >>> soil_heat_flux_hourly_auto(-0.4) -0.2
- soil_heat.soil_heat_flux_monthly_series(t_monthly: numpy.typing.ArrayLike, method: str = 'centered') numpy.ndarray[source]
Compute soil heat flux for each month in a temperature time series.
This is a convenience function that applies the appropriate FAO-56 monthly equation to each element of a monthly mean-temperature array.
- Parameters:
t_monthly (array_like) – 1-D array of consecutive monthly mean air temperatures [°C]. Length N ≥ 2.
method ({"centered", "backward"}, optional) –
"centered"(default) — Uses Eq. 43 where both neighbors exist; falls back to Eq. 44 at the boundaries (first and last month)."backward"— Uses Eq. 44 for every month (each month compared only to the previous month). The first element is set toNaNbecause there is no previous month.
- Returns:
Array of monthly soil heat flux values [MJ m⁻² day⁻¹], same length as t_monthly.
- Return type:
numpy.ndarray
Examples
FAO-56 Example 13 — March, April, May temps:
>>> temps = [14.1, 16.1, 18.8] >>> soil_heat_flux_monthly_series(temps, method="centered") array([0.28 , 0.329, 0.378])
Using backward-only differences:
>>> soil_heat_flux_monthly_series(temps, method="backward") array([ nan, 0.28 , 0.378])
- soil_heat.soil_heat_flux_annual_cycle(t_monthly_12: numpy.typing.ArrayLike) numpy.ndarray[source]
Monthly G for a complete 12-month annual cycle using Eq. 43.
Wraps around so that January uses December and February as neighbors, and December uses November and January.
- Parameters:
t_monthly_12 (array_like) – Exactly 12 monthly mean air temperatures [°C], January through December.
- Returns:
Shape
(12,)array of monthly soil heat flux [MJ m⁻² day⁻¹].- Return type:
numpy.ndarray
Examples
>>> temps = [5.2, 6.1, 9.3, 13.0, 17.5, 22.1, ... 25.4, 24.8, 20.6, 14.9, 9.2, 5.8] >>> g = soil_heat_flux_annual_cycle(temps) >>> g[0] # January: 0.07 * (Feb - Dec) 0.021...
- soil_heat.energy_to_evap(energy: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Convert energy flux to equivalent evaporation depth (FAO-56 Eq. 20).
\[E_{\text{equiv}} = 0.408 \times E\]
- soil_heat.evap_to_energy(evap: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Convert equivalent evaporation depth back to energy flux.
Inverse of
energy_to_evap().
- soil_heat.mj_m2_day_to_w_m2(flux: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Convert flux from MJ m⁻² day⁻¹ to W m⁻².
\[W\,m^{-2} = \frac{MJ\,m^{-2}\,day^{-1}}{0.0864}\]
- soil_heat.w_m2_to_mj_m2_day(flux: float | numpy.typing.ArrayLike) float | numpy.ndarray[source]
Convert flux from W m⁻² to MJ m⁻² day⁻¹.
- soil_heat.compute_soil_storage_integrated(df: pandas.DataFrame, col_T5: str = 'T5cm', col_VWC5: str = 'VWC5cm', col_IRT: str = None, depth: float = 0.05, bulk_density: float = 1300.0, specific_heat_dry: float = 870.0, specific_heat_water: float = 4218.0) pandas.Series[source]
Calculate the soil heat storage flux (S) in W/m² using a change-in-storage approach.
This function estimates the energy stored or released in the upper soil layer by accounting for the heat capacity of both the dry soil minerals and the stored liquid water. It supports an integrated temperature approach if infrared surface temperature (IRT) is available. In the slab approach, you assume the estimate at the given depth applies to the whole slab above the depth, whereas the IRT approach takes the mean between surface temperature and temperature at depth.
Note that this method is the same as the method use to calculate soil heat storage flux in the Campbell Scientific programs, though that progrmam uses higher-frequency data for the calculations.
- Parameters:
df (pd.DataFrame) – Input dataframe with a DatetimeIndex.
col_T5 (str, default "T5cm") – Column name for soil temperature at 5cm depth [°C].
col_VWC5 (str, default "VWC5cm") – Column name for Volumetric Water Content at 5cm depth [m³/m³].
col_IRT (str, optional) – Column name for surface infrared temperature [°C]. If provided, the layer temperature is the average of surface and 5cm temps. Note that this must be surface temperature, not canopy temperature
depth (float, default 0.05) – The thickness of the soil layer being measured [m].
bulk_density (float, default 1300.0) – Dry bulk density of the soil [kg/m³].
specific_heat_dry (float, default 870.0) – Specific heat capacity of dry soil minerals [J/(kg·K)].
specific_heat_water (float, default 4218.0) – Specific heat capacity of water [J/(kg·K)].
- Returns:
The calculated soil heat storage flux [W/m²]. Positive values indicate energy moving into storage (warming).
- Return type:
pd.Series
Notes
The canopy storage flux is calculated as: $$S_{canopy} = frac{m_{veg} cdot C_{p} cdot Delta T}{Delta t}$$
Where: - $m_{veg}$ = Fresh biomass ($kg/m^2$), calculated as $(height / 100) cdot density_factor$ - $C_{p}$ = Specific heat of fresh vegetation ($approx 3900 text{ J/kg}cdottext{K}$) - $Delta T$ = Change in canopy temperature between time steps - $Delta t$ = Time step in seconds
- soil_heat.compute_canopy_storage(df: pandas.DataFrame, col_irt: str = 'T_CANOPY', crop_type: str = 'alfalfa', height_cm: str | float = None, density_factor: float = None) pandas.Series[source]
Calculate the energy storage flux of a plant canopy (S_canopy) in W/m².
This function estimates heat storage within biomass. It can handle dynamic vegetation growth if a height column is provided, otherwise it falls back to static heights. Height and density values can be entered by the user or taken from a dictionary within the function.
Care should be taken to examine the crop default values and the specific heat of the vegetation.
- Parameters:
df (pd.DataFrame) – Input dataframe with a DatetimeIndex.
col_irt (str, default "T_CANOPY") – Column name for the canopy surface temperature [°C].
crop_type (str, default "alfalfa") – Options: ‘alfalfa’, ‘corn’, ‘hay’, ‘wetland’. Used to determine defaults for height and density.
height_cm (str or float, optional) – If str: The column name in df containing canopy height [cm] over time. If float: A static height value [cm]. If None: Uses default height for the specified crop_type.
density_factor (float, optional) – Biomass density [kg/m² per meter height]. If None, uses default for the specified crop_type.
- Returns:
Energy storage flux [W/m²].
- Return type:
pd.Series
- soil_heat.volumetric_heat_capacity(theta, porosity=POROSITY, f_org=F_ORGANIC)[source]
Estimate soil volumetric heat capacity using the de Vries (1963) linear mixing model.
Computes Cv as a volume-fraction-weighted sum of the heat capacities of the four soil constituents — mineral solids, organic matter, liquid water, and air:
Cv(θ) = f_min·Cv_s + f_org·Cv_org + θ·Cv_w + θ_a·Cv_a
where the air-filled porosity θ_a = max(porosity - θ, 0).
- Parameters:
theta (array-like) – Volumetric water content [m³ m⁻³]. Accepts scalars, lists, or NumPy arrays. Values above porosity are implicitly handled by flooring air-filled porosity at zero.
porosity (float, optional) – Total soil porosity [m³ m⁻³]. Used to derive both the mineral solid fraction (f_min = 1 - porosity - f_org) and the air-filled pore fraction. Defaults to the module-level constant
POROSITY.f_org (float, optional) – Volumetric fraction of organic matter in the solid phase [m³ m⁻³]. Defaults to the module-level constant
F_ORGANIC.
- Returns:
Volumetric heat capacity Cv [J m⁻³ K⁻¹] with the same shape as
theta. Typical values range from ~1.0×10⁶ J m⁻³ K⁻¹ (dry) to ~3.0×10⁶ J m⁻³ K⁻¹ (saturated mineral soil).- Return type:
np.ndarray
Notes
Component heat capacities are taken from module-level constants:
CV_S— mineral solids (~2.0×10⁶ J m⁻³ K⁻¹)CV_ORG— organic matter (~2.5×10⁶ J m⁻³ K⁻¹)CV_W— liquid water (~4.18×10⁶ J m⁻³ K⁻¹)CV_A— air (~1.2×10³ J m⁻³ K⁻¹)
The air contribution is negligible at most water contents but is retained for physical completeness. The model assumes no phase change and is not valid for frozen soils where ice content would require a separate term.
References
de Vries, D. A. (1963). Thermal properties of soils. In W. R. van Wijk (Ed.), Physics of Plant Environment (pp. 210–235). North-Holland, Amsterdam.
- soil_heat.thermal_conductivity_johansen(theta, porosity=POROSITY, f_quartz=F_QUARTZ)[source]
Estimate soil thermal conductivity using the Johansen (1975) model with the Kersten number formulation from Côté & Konrad (2005).
Computes λ by interpolating between dry and saturated end-member conductivities using the Kersten number (Ke) as a saturation-dependent weighting function:
λ(θ) = λ_dry + Ke(Sr) · (λ_sat - λ_dry)
End members and saturation are derived from soil texture (quartz fraction), porosity, and standard mineral/water conductivity constants.
- Parameters:
theta (array-like) – Volumetric water content [m³ m⁻³]. Accepts scalars, lists, or NumPy arrays. Clipped internally to the range [0, porosity] via the degree of saturation calculation.
porosity (float, optional) – Total soil porosity [m³ m⁻³], i.e., the volumetric fraction of pore space. Defaults to the module-level constant
POROSITY.f_quartz (float, optional) – Quartz fraction of the solid phase [dimensionless, 0–1]. Higher quartz content increases the solid-phase conductivity λ_s and therefore λ_sat. Defaults to the module-level constant
F_QUARTZ.
- Returns:
Thermal conductivity λ [W m⁻¹ K⁻¹] with the same shape as
theta.- Return type:
np.ndarray
Notes
The three derived end-member quantities are:
λ_dry: dry soil conductivity estimated from bulk density using the empirical formula of Johansen (1975), assuming a particle density of 2700 kg m⁻³.
λ_s: solid-phase conductivity computed as a geometric mean weighted by quartz fraction, using module constants
LAMBDA_S(mineral) and 2.0 W m⁻¹ K⁻¹ for non-quartz minerals.λ_sat: saturated conductivity computed as a geometric mean of the solid and water (
LAMBDA_W) phase conductivities weighted by their volume fractions (1 - porosity and porosity, respectively).
The Kersten number follows the log-linear relationship for unfrozen mineral soils from Côté & Konrad (2005):
Ke = 0.7 · log10(Sr) + 1.0 for Sr > 0.05 Ke = 0.0 for Sr ≤ 0.05
Ke is clipped to [0, 1] and Sr is floored at 0.01 inside the log to prevent numerical issues at very low saturation. This formulation is valid for unfrozen, mineral soils; it is not appropriate for frozen conditions or highly organic soils.
References
Johansen, O. (1975). Thermal conductivity of soils. PhD thesis, University of Trondheim, Norway. (CRREL Draft Translation 637, 1977).
Côté, J., & Konrad, J.-M. (2005). A generalized thermal conductivity model for soils and construction materials. Canadian Geotechnical Journal, 42(2), 443–458. https://doi.org/10.1139/t04-106
- soil_heat.thermal_conductivity_lu2007(theta)[source]
Estimate soil thermal conductivity using the Lu et al. (2007) empirical model.
Computes λ as a function of volumetric water content using a nonlinear empirical equation fitted across a wide range of mineral soil textures:
λ(θ) = a + b·θ - (a - d) · exp(-(c·θ)^e)
with coefficients a=0.56, b=0.51, c=0.96, d=0.27, e=0.87.
- Parameters:
theta (array-like) – Volumetric water content [m³ m⁻³]. Accepts scalars, lists, or NumPy arrays. Values are expected in the range [0, porosity]; behavior outside this range is not physically meaningful.
- Returns:
Thermal conductivity λ [W m⁻¹ K⁻¹] with the same shape as
theta.- Return type:
np.ndarray
Notes
The model was calibrated on 36 soils spanning sandy to clayey textures and dry to saturated conditions, with reported RMSE of ~0.1 W m⁻¹ K⁻¹. It does not account for soil texture, organic matter content, or temperature explicitly, so accuracy may degrade for organic-rich or highly heterogeneous soils.
References
Lu, S., Ren, T., Gong, Y., & Horton, R. (2007). An improved model for predicting soil thermal conductivity from water content at room temperature. Soil Science Society of America Journal, 71(1), 8–14. https://doi.org/10.2136/sssaj2006.0041
- soil_heat.gradient_surface(ts, swc, df)[source]
Estimate surface ground heat flux using a temperature gradient from the canopy surface down to the 5 cm soil sensor (Method 2).
Applies Fourier’s law over the 0 → 5 cm layer, treating the canopy air temperature as a proxy for the soil surface temperature:
G0 = -λ(θ_5cm) · (T_5cm - T_surface) / dz
where dz = 0.05 m.
- Parameters:
ts (pd.DataFrame) – Soil temperature time series [°C], with column keys as sensor depths in cm. Must contain a column for the 5.0 cm depth.
swc (pd.DataFrame) – Volumetric soil water content time series [m³ m⁻³], with column keys matching those of
ts. The 5.0 cm column is used to compute thermal conductivity.df (pd.DataFrame) – Main eddy covariance / met data DataFrame. Must contain the column
"T_CANOPY_1_1_1"[°C], used as a proxy for the soil surface temperature.
- Returns:
Estimated surface ground heat flux [W m⁻²], indexed to match
ts. Named"G_gradient_surface". Positive values indicate downward flux (heat flowing into the soil).- Return type:
pd.Series
Notes
Using canopy air temperature as a surface proxy introduces uncertainty under conditions where the canopy-to-soil surface temperature decouples significantly (e.g., high radiative loading, dense litter layers, or sparse canopy cover). Thermal conductivity is estimated at the 5 cm depth via the Johansen (1975) model, which may underestimate conductivity in very dry or coarse-textured soils.
References
Johansen, O. (1975). Thermal conductivity of soils. PhD thesis, University of Trondheim, Norway.
- soil_heat.compute_storage_above_plate(ts, swc, plate_depth_cm=PLATE_DEPTH_CM, dt=DT_SEC)[source]
Compute soil heat storage in the layer from the surface to the heat flux plate depth.
Integrates the calorimetric storage term (Cv * dT/dt * dz) over the 0 → plate_depth layer using the 5 cm temperature and soil water content sensors as representative values for that layer. The result can be added to the plate-measured soil heat flux (SG) to estimate total ground heat flux at the surface (G0 = SG + S_above_plate).
Since the plate is typically installed at 8 cm — between the 5 cm and 10 cm sensors — the 5 cm sensor is used as a pragmatic representative of the above-plate layer rather than performing a full depth interpolation.
- Parameters:
ts (pd.DataFrame) – Soil temperature time series [°C], with column keys as sensor depths in cm (e.g., 5.0, 10.0). Must contain a column for the 5.0 cm depth.
swc (pd.DataFrame) – Volumetric soil water content time series [m³ m⁻³], with column keys matching those of
ts. Must contain a column for the 5.0 cm depth.plate_depth_cm (float, optional) – Depth of the soil heat flux plate [cm]. Defines the thickness of the storage layer (dz = plate_depth_cm / 100). Defaults to
PLATE_DEPTH_CM.dt (float, optional) – Timestep duration [s] used to compute dT/dt via finite differencing. Defaults to
DT_SEC.
- Returns:
Heat storage rate in the 0 → plate_depth layer [W m⁻²], indexed to match
ts. Named"S_above_plate". Positive values indicate heat being stored (warming), negative values indicate heat being released (cooling). The first element will be NaN due to the forward-difference computation.- Return type:
pd.Series
Notes
Storage is calculated as:
S = Cv(θ) · (ΔT / Δt) · dz
where Cv is the volumetric heat capacity [J m⁻³ K⁻¹] derived from the soil water content at 5 cm, ΔT is the temperature change over one timestep, and dz is the layer thickness in meters.
- soil_heat.gradient_plus_storage(ts, swc, ref_depth_idx=2, dt=DT_SEC, lam_model='johansen')[source]
Estimate surface ground heat flux using a conductive gradient at a reference depth plus calorimetric storage integrated over all layers above it.
Applies the standard soil heat flux correction:
G0 = G(z_ref) + S(0 → z_ref)
The conductive flux at the reference depth is estimated from Fourier’s law using temperatures at the sensor depths bracketing
z_ref(i.e., one layer above and one layer below). The storage term is computed calorimetrically for each discrete layer from the surface down to and includingz_ref, then summed.- Parameters:
ts (pd.DataFrame) – Soil temperature time series [°C], with column keys as sensor depths in cm matching
DEPTHS_CM. Must span at least the depths indexed byref_depth_idx - 1throughref_depth_idx + 1.swc (pd.DataFrame) – Volumetric soil water content time series [m³ m⁻³], with column keys matching those of
ts. Used both to compute thermal conductivity at the reference depth and volumetric heat capacity for all above-reference layers.ref_depth_idx (int, optional) – Index into
DEPTHS_CMidentifying the reference depth at which the conductive gradient is evaluated. Depths atref_depth_idx - 1andref_depth_idx + 1are used as the bracketing pair for the gradient. Defaults to 2, typically corresponding to the 10 cm sensor.dt (float, optional) – Timestep duration [s] used to compute dT/dt via finite differencing. Defaults to
DT_SEC.lam_model ({"johansen", "lu2007"}, optional) – Thermal conductivity model to apply at the reference depth.
"johansen"uses the Johansen (1975) mixing model;"lu2007"uses the Lu et al. (2007) model. Defaults to"johansen".
- Returns:
Estimated surface ground heat flux [W m⁻²], indexed to match
ts. Named"G_grad+stor_{z_ref:.0f}cm"wherez_refis the reference depth in cm (e.g.,"G_grad+stor_10cm"). Positive values indicate downward flux (heat flowing into the soil). The first element will be NaN due to the finite-difference storage calculation.- Return type:
pd.Series
Notes
Conductive flux at the reference depth is computed as:
G(z_ref) = -λ(θ_ref) · (T_below - T_above) / dz_grad
where
dz_gradis the distance in meters between the two bracketing sensor depths and λ is thermal conductivity [W m⁻¹ K⁻¹] estimated from VWC at the reference depth.Storage in each above-reference layer is:
S_i = Cv(θ_i) · (ΔT_i / Δt) · dz_i
and the total above-reference storage is:
S(0 → z_ref) = Σ S_i for i = 0 … ref_depth_idx
Layer thicknesses
dz_iare taken from the module-levelLAYER_THICK_Marray. NaN values in individual layer storage terms are excluded from the column-wise sum vianp.nansum.References
Johansen, O. (1975). Thermal conductivity of soils. PhD thesis, University of Trondheim, Norway.
Lu, S., Ren, T., Gong, Y., & Horton, R. (2007). An improved model for predicting soil thermal conductivity from water content at room temperature. Soil Science Society of America Journal, 71(1), 8–14.