Note: All code in this vignette requires TensorFlow and Keras to be installed. See
?soilFluxfor installation instructions. The outputs shown are representative examples from a training run on theswrc_exampledataset.
The soilFlux package implements a physics-informed one-dimensional convolutional neural network (CNN1D-PINN) for estimating the complete soil water retention curve (SWRC) as a continuous function of matric potential. The architecture and physics constraints are adapted from Norouzi et al. (2025).
The soil water retention curve describes the relationship between volumetric water content (\(\theta\), m³/m³) and matric potential (expressed as pF = log₁₀(h [cm])). It is a key input to most hydrological and land-surface models.
| Feature | Description |
|---|---|
| Monotone by construction | θ decreases with pF by architecture (cumulative integral of positive slopes) |
| Physics-informed | Four residual constraints enforce known physical behaviour |
| Continuous output | Predicts θ at any pF value, not just standard pressure heads |
| Multi-covariate | Uses texture (sand/silt/clay fractions), bulk density, SOC, and depth |
| Ready to save/load | Weights persisted via Keras HDF5 for spatial prediction |
library(soilFlux)
# 1. Load and prepare data --------------------------------------------------
data("swrc_example") # example dataset bundled with the package
df <- prepare_swrc_data(
swrc_example,
depth_col = "depth",
fix_bd = TRUE,
fix_theta = TRUE
)
# Train / validation / test split (by PEDON_ID)
ids <- unique(df$PEDON_ID)
set.seed(42)
tr_ids <- sample(ids, floor(0.7 * length(ids)))
val_ids <- sample(setdiff(ids, tr_ids), floor(0.15 * length(ids)))
te_ids <- setdiff(ids, c(tr_ids, val_ids))
train_df <- df[df$PEDON_ID %in% tr_ids, ]
val_df <- df[df$PEDON_ID %in% val_ids, ]
test_df <- df[df$PEDON_ID %in% te_ids, ]
# 2. Fit model --------------------------------------------------------------
fit <- fit_swrc(
train_df = train_df,
x_inputs = c("clay", "silt", "bd_gcm3", "soc", "Depth_num"),
val_df = val_df,
epochs = 60,
batch_size = 256,
lambdas = norouzi_lambdas("norouzi"),
verbose = TRUE
)
# 3. Evaluate ---------------------------------------------------------------
m <- evaluate_swrc(fit, test_df)
print(m)
# 4. Dense SWRC curves ------------------------------------------------------
dense <- predict_swrc_dense(fit, newdata = test_df, n_points = 500)
# 5. Plot -------------------------------------------------------------------
plot_swrc(dense, obs_points = test_df,
obs_col = "theta_n",
facet_row = "Depth_label",
facet_col = "Texture")prepare_swrc_data() expects (at minimum):
| Column | Description | Units |
|---|---|---|
PEDON_ID |
Unique profile identifier | — |
depth |
Depth interval (e.g. "0-5") |
cm |
matric_head |
Matric head | cm |
water_content |
Volumetric water content | m³/m³ or % |
| Covariates | e.g. clay, silt, bd,
soc |
%, g/cm³ |
The package automatically detects and corrects common unit issues:
# Bulk density: if median > 10, assumes kg/m3 and converts to g/cm3
fix_bd_units(c(1200, 1450, 1300)) # kg/m3 -> g/cm3
# Theta: if max > 1.5, assumes % and divides by 100
theta_unit_factor(c(30, 40, 25)) # returns 100The model output is:
\[\hat{\theta}(\text{pF}) = \theta_s - \int_0^{\text{pF}} \text{softplus}(s(t)) \, dt\]
where \(s(t)\) is a 1-channel Conv1D output. Because softplus > 0 always, the integral is strictly increasing, so \(\hat{\theta}\) is strictly decreasing — monotonicity is guaranteed by construction, not by post-processing.
Each observation is encoded as a 3-D sequence: * Shape:
[N, K, p+1] * K = number of knot points
(default 64) * p = number of covariates * Last channel =
knot position in [0, 1]
This design broadcasts each covariate across all knot positions, allowing the Conv1D layers to learn shape-from-texture mappings.
Four residual constraints from Norouzi et al. (2025) are embedded in the loss function:
| Set | Constraint | pF range | Default weight |
|---|---|---|---|
| S1 | Linearity at dry end: ∣∂²θ/∂pF²∣ → 0 | [5.0, 7.6] | λ₃ = 1 |
| S2 | Non-negativity: θ(pF = 6.2) ≥ 0 | fixed | λ₄ = 1000 |
| S3 | Non-positivity: θ(pF = 7.6) ≤ 0 | fixed | λ₅ = 1000 |
| S4 | Flat plateau: ∣∂θ/∂pF∣ → 0 | [−2.0, −0.3] | λ₆ = 1 |
Plus two data-loss weights:
| Term | Description | Default |
|---|---|---|
| λ₁ | Wet-end data loss (pF ≤ 4.2) | 1 |
| λ₂ | Dry-end data loss (pF > 4.2) | 10 |
# Exact replication of Norouzi et al. (2025) Table 1
norouzi_lambdas("norouzi")
# Smoother dry-end (lambda3 = 10)
norouzi_lambdas("smooth")# Save
save_swrc_model(fit, dir = "./models", name = "model_5")
# Check it exists
swrc_model_exists("./models", "model_5")
# Reload (in a new session, after library(soilFlux))
fit2 <- load_swrc_model("./models", "model_5")
pred <- predict_swrc(fit2, newdata = test_df)classify_texture(
sand = c(70, 20, 10),
silt = c(15, 50, 30),
clay = c(15, 30, 60)
)
# Add as a column to a data frame
add_texture(df, sand_col = "sand_total")
# Texture triangle (requires ggtern)
texture_triangle(df, color_col = "Texture")Norouzi, A. M., Feyereisen, G. W., Papanicolaou, A. N., & Wilson, C. G. (2025). Physics-Informed Neural Networks for Estimating a Continuous Form of the Soil Water Retention Curve. Journal of Hydrology.