Jointly fit Van Genuchten SWRC and hydraulic conductivity parameters

Fits the Van Genuchten (1980) soil water retention curve and Mualem-Van Genuchten hydraulic conductivity function simultaneously to observed (h, θ) and (h, K) data. This joint fit shares the common structural parameters (θ_r, θ_s, α, n) between both models, typically improving identifiability compared to fitting them separately.

Usage

fit_swrc_hcc(
  data,
  theta_col,
  K_col,
  h_col,
  lower = NULL,
  upper = NULL,
  fixed = NULL,
  wrc_weight = 1,
  hcc_weight = 1,
  workers = 1L
)

Arguments

data: A data frame or tibble, optionally grouped with dplyr::group_by().
theta_col: Bare column name of observed volumetric water content (m³/m³). Rows where this is NA are excluded from the WRC component.
K_col: Bare column name of observed hydraulic conductivity. Rows where this is NA are excluded from the HCC component. Must be positive; the log₁₀ transformation is applied internally.
h_col: Bare column name of matric potential / pressure head. Absolute values are used internally.
lower: Named numeric vector of lower parameter bounds. Parameter names: "theta_r", "theta_s", "alpha", "n", "Ks", "tau".
upper: Named numeric vector of upper parameter bounds.
fixed: Named numeric vector of parameters to hold fixed. E.g. fixed = c(tau = 0.5) (the default behaviour if tau is not supplied).
wrc_weight: Positive scalar weighting the WRC residual component in the joint objective. Defaults to 1.
hcc_weight: Positive scalar weighting the HCC (log K) residual component. Defaults to 1.
workers: Number of parallel workers. Defaults to 1.

Value

An object of class fit_swrc_hcc (a tibble subclass) with one row per group containing:

Group keys (if grouped)
theta_r, theta_s, alpha, n, Ks, tau: fitted values
std_error_* for each free parameter (approximate, from numerical Hessian)
convergence: TRUE if optim() returned code 0

Details

The objective function minimises the combined weighted residual: $$RSS = w_\theta \sum ({\theta_{obs}} - \theta_{pred})^2 + w_K \sum (\log_{10}K_{obs} - \log_{10}K_{pred})^2$$

K residuals are on the log₁₀ scale because K spans many orders of magnitude; this prevents high-K observations from dominating the fit. The wrc_weight and hcc_weight scalars allow the user to balance the two components.

Missing values (NA) in theta_col or K_col are silently dropped from the respective component, so WRC and HCC measurements need not be at the same h values and can be in the same data frame with NAs where data is absent.

Parameters estimated

Parameter	Meaning
`theta_r`	Residual water content (m³/m³)
`theta_s`	Saturated water content (m³/m³)
`alpha`	Van Genuchten shape (1/cm or 1/m)
`n`	Van Genuchten shape (dimensionless; > 1)
`Ks`	Saturated hydraulic conductivity (units of K_col)
`tau`	Mualem tortuosity (dimensionless; fixed at 0.5 by default)

References

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

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

Examples

library(tibble)

# Synthetic joint dataset — WRC and K at the same h values
joint <- tibble(
  h     = c(0, 10, 100, 500, 1000, 5000, 15000),
  theta = c(0.44, 0.40, 0.32, 0.25, 0.20, 0.12, 0.08),
  K     = c(25, 8.5, 0.9, 0.08, 0.012, 4e-4, 3e-5)
)

fit_swrc_hcc(joint, theta_col = theta, K_col = K, h_col = h)
#> # A tibble: 1 × 13
#>   theta_r theta_s   alpha     n    Ks   tau std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>             <dbl>             <dbl>
#> 1       0   0.401 0.00748  1.32  25.0 -1.13            0.0448            0.0252
#> # ℹ 5 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   std_error_Ks <dbl>, std_error_tau <dbl>, convergence <lgl>

# Fix tau at Mualem value; estimate the remaining 5 parameters
fit_swrc_hcc(joint, theta_col = theta, K_col = K, h_col = h,
             fixed = c(tau = 0.5))
#> # A tibble: 1 × 13
#>   theta_r theta_s   alpha     n    Ks   tau std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>             <dbl>             <dbl>
#> 1       0   0.347 0.00433  1.19  29.7   0.5             0.105            0.0385
#> # ℹ 5 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   std_error_Ks <dbl>, std_error_tau <dbl>, convergence <lgl>

# Mixed h values: WRC and HCC measured at different points
mixed <- tibble(
  h     = c(0, 10, 100, 500, 1000, 5000,  10, 100, 1000),
  theta = c(0.44, 0.40, 0.32, 0.25, 0.20, 0.12, NA, NA, NA),
  K     = c(NA, NA, NA, NA, NA, NA, 8.5, 0.9, 0.012)
)

fit_swrc_hcc(mixed, theta_col = theta, K_col = K, h_col = h)
#> # A tibble: 1 × 13
#>   theta_r theta_s  alpha     n    Ks   tau std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl>             <dbl>             <dbl>
#> 1  0.0276   0.415 0.0111  1.34  28.8 -1.73            0.0531            0.0107
#> # ℹ 5 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   std_error_Ks <dbl>, std_error_tau <dbl>, convergence <lgl>