Compute soil water diffusivity (Van Genuchten / Mualem)

Calculates soil water diffusivity D(h) at given matric potential values: $$D(h) = \frac{K(h)}{C(h)}$$

Usage

soil_water_diffusivity(data, ks, theta_r, theta_s, alpha, n, m, h, tau = 0.5)

Arguments

data: A data frame or tibble. Passed as the first argument so the function is compatible with |> and %>% pipes.
ks: Saturated hydraulic conductivity (e.g. cm/day). Either a bare column name from data or a scalar numeric value. Must be strictly positive.
theta_r: Residual water content (m³/m³). Either a bare column name or a scalar numeric value.
theta_s: Saturated water content (m³/m³). Either a bare column name or a scalar numeric value.
alpha: Van Genuchten shape parameter (1/cm or 1/m, matching units of h). Either a bare column name or a scalar numeric value. Must be strictly positive.
n: Van Genuchten shape parameter (dimensionless). Must be > 1 for the Mualem model. Either a bare column name or a scalar numeric value.
m: Van Genuchten shape parameter (dimensionless). Typically 1 - 1/n. Must satisfy 0 < m < 1. Either a bare column name or a scalar numeric value. Used only in the K(h) formula; C(h) derives m from n internally.
h: Matric potential / pressure head (cm or m, must match units of alpha). Either a bare column name or a scalar numeric value. Sign convention: absolute value is used internally.
tau: Tortuosity/connectivity parameter (dimensionless). Defaults to 0.5 (Mualem 1976). Must satisfy tau > -2. Either a bare column name or a scalar numeric value.

Value

The input data as a tibble with an additional column .D containing soil water diffusivity in units of length²/time (e.g. cm²/day when ks is in cm/day and h in cm).

Details

where K(h) is the Mualem-Van Genuchten hydraulic conductivity and C(h) is the Van Genuchten soil water capacity (−dθ/dh). D(h) has units of length² per time (e.g. cm²/day when ks is in cm/day and h in cm).

D(h) appears in the diffusive form of the Richards equation: $$\frac{\partial \theta}{\partial t} = \nabla \cdot \left(D(\theta) \nabla \theta\right) - \frac{\partial K}{\partial z}$$

and is used in analytical solutions for infiltration and redistribution.

Note

D(h) is undefined at saturation (h = 0) because C(0) = 0. The function returns Inf at h = 0 unless K(0) is also zero, in which case R returns NaN via 0/0. Handle h = 0 rows in advance if needed.

References

Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3), 513–522. https://doi.org/10.1029/WR012i003p00513

Van Genuchten, M. Th. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x

Examples

library(tibble)

# Exclude h = 0 (D is undefined at saturation)
soil <- tibble(h = c(1, 10, 100, 1000, 15000))

soil_water_diffusivity(
  soil,
  ks      = 10,
  theta_r = 0.05,
  theta_s = 0.45,
  alpha   = 0.02,
  n       = 1.5,
  m       = 1 - 1/1.5,
  h       = h
)
#> # A tibble: 5 × 2
#>       h         .D
#>   <dbl>      <dbl>
#> 1     1 13078.    
#> 2    10  1975.    
#> 3   100    78.0   
#> 4  1000     1.47  
#> 5 15000     0.0128