Fit Van Genuchten parameters to observed retention data

Fits the Van Genuchten (1980) soil water retention model to observed (h, θ) data. Returns a tidy tibble of estimated parameters. If the input data is grouped (via dplyr::group_by()), the fit is performed independently for each group.

Usage

fit_swrc(
  data,
  theta_col,
  h_col,
  workers = 1L,
  lower = NULL,
  upper = NULL,
  fixed = NULL,
  weights = NULL
)

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³).
h_col: Bare column name of observed matric potential / pressure head. Absolute values are used internally.
workers: Number of parallel workers. Defaults to 1 (sequential). Values > 1 use parallel::mclapply() on Unix. On Windows, silently reduced to 1.
lower: Named numeric vector of lower bounds for parameters. Names should be a subset of c("theta_r", "theta_s", "alpha", "n"). Unspecified parameters receive physical defaults (0, 1e-6, 1e-6, 1.001 respectively).
upper: Named numeric vector of upper bounds. Same naming convention as lower. Unspecified parameters default to (1, 1, Inf, Inf).
fixed: Named numeric vector of parameters to hold at fixed values rather than estimate. E.g. fixed = c(theta_r = 0.05). Fixed parameters are excluded from the optimisation and reported with NA standard errors.
weights: Bare column name or numeric scalar; per-observation weights for the least-squares objective. Larger weights increase influence of that observation. Defaults to uniform weights.

Value

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

Group keys (if input was grouped)
theta_r, theta_s, alpha, n: fitted parameter values (or the supplied value for fixed parameters)
std_error_theta_r, std_error_theta_s, std_error_alpha, std_error_n: standard errors (NA for fixed parameters)
convergence: TRUE if optimisation converged

The result stores the original data and fitting options as attributes, enabling confint.fit_swrc() to compute profile-likelihood confidence intervals.

Details

Fitting paths

The function selects an optimisation back-end based on the arguments supplied:

Arguments	Back-end
Any combination without `fixed`	`nls(algorithm = "port")` — box-constrained NLS with physical default bounds; SE from NLS Jacobian
`fixed` specified (+ optional bounds + `weights`)	`stats::optim(method = "L-BFGS-B")` — reduced-dimension optimisation; SE from numerical Hessian (approximate)

Starting values

Auto-initialised from the data:

theta_r: 90% of minimum observed θ
theta_s: 105% of maximum observed θ
alpha: 0.1
n: 2.0

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

Examples

library(tibble)

observed <- 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)
)

# Default unconstrained fit
fit_swrc(observed, theta_col = theta, h_col = h)
#> # A tibble: 1 × 9
#>   theta_r theta_s  alpha     n std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>  <dbl> <dbl>             <dbl>             <dbl>
#> 1       0   0.425 0.0182  1.27             0.124            0.0164
#> # ℹ 3 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   convergence <lgl>

# Box-constrained fit (n must be between 1.1 and 4)
fit_swrc(observed, theta_col = theta, h_col = h,
         lower = c(n = 1.1), upper = c(n = 4))
#> # A tibble: 1 × 9
#>   theta_r theta_s  alpha     n std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>  <dbl> <dbl>             <dbl>             <dbl>
#> 1       0   0.425 0.0182  1.27             0.124            0.0164
#> # ℹ 3 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   convergence <lgl>

# Fix theta_r at a known value; estimate the other three parameters
fit_swrc(observed, theta_col = theta, h_col = h,
         fixed = c(theta_r = 0.05))
#> # A tibble: 1 × 9
#>   theta_r theta_s  alpha     n std_error_theta_r std_error_theta_s
#>     <dbl>   <dbl>  <dbl> <dbl>             <dbl>             <dbl>
#> 1    0.05   0.422 0.0144  1.37                NA            0.0124
#> # ℹ 3 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   convergence <lgl>

# Grouped fit
library(dplyr)
tibble(
  horizon = rep(c("A", "B"), each = 7),
  h       = rep(c(0, 10, 100, 500, 1000, 5000, 15000), 2),
  theta   = c(0.44, 0.40, 0.32, 0.25, 0.20, 0.12, 0.08,
              0.38, 0.33, 0.27, 0.21, 0.16, 0.10, 0.07)
) |>
  group_by(horizon) |>
  fit_swrc(theta_col = theta, h_col = h)
#> # A tibble: 2 × 10
#>   horizon theta_r theta_s  alpha     n std_error_theta_r std_error_theta_s
#>   <chr>     <dbl>   <dbl>  <dbl> <dbl>             <dbl>             <dbl>
#> 1 A             0   0.425 0.0182  1.27             0.124            0.0164
#> 2 B             0   0.361 0.0215  1.26             0.126            0.0171
#> # ℹ 3 more variables: std_error_alpha <dbl>, std_error_n <dbl>,
#> #   convergence <lgl>