🌱 RootWaterSim

Interactive simulator of 1D soil water flow and root water uptake Built for teaching soil physics, hydrology, and precision agriculture.

---

Overview

RootWaterSim solves the Richards equation (mixed form) in one dimension, coupled with the Feddes root water uptake model. It provides:

• A Streamlit web app for interactive exploration of soil moisture dynamics.
• Jupyter notebooks that teach the underlying theory step-by-step.
• A clean, modular Python library that can be reused in research or teaching.

---

Governing Equations

1. Richards Equation (Mixed Form)

The one-dimensional Richards equation describes variably saturated water flow in soil:

∂θ/∂t = ∂/∂z [ K(h) ( ∂h/∂z + 1 ) ] - S(z, h)

where:

• θ = volumetric water content (cm³/cm³)
• t = time (days)
• z = depth (cm, positive downward)
• h = pressure head (cm; negative in the unsaturated zone)
• K(h) = unsaturated hydraulic conductivity (cm/day)
• S(z, h) = root water uptake sink term (cm/day)

The term (∂h/∂z + 1) accounts for both the hydraulic gradient and the gravitational component.

2. Van Genuchten–Mualem Hydraulic Properties

Soil hydraulic properties are computed using the van Genuchten (1980) model for water retention and the Mualem (1976) model for unsaturated conductivity.

Effective saturation:

Sₑ(h) = 1 / (1 + (α |h|)ⁿ)ᵐ

Water content:

θ(h) = θᵣ + (θₛ - θᵣ) · Sₑ(h)

Unsaturated hydraulic conductivity:

K(h) = Kₛ · Sₑ^(1/2) · [1 - (1 - Sₑ^(1/m))ᵐ]²

Capillary capacity:

C(h) = dθ/dh = m·n·αⁿ·|h|^(n-1)·(θₛ - θᵣ) · Sₑ^(1/m + 1)    for h < 0
C(h) = 0                                                        for h ≥ 0

where:

Symbol	Meaning	Unit
θᵣ	Residual water content	cm³/cm³
θₛ	Saturated water content	cm³/cm³
α	Inverse air-entry value	1/cm
n	Pore-size distribution index (>1)	—
m	1 - 1/n	—
Kₛ	Saturated hydraulic conductivity	cm/day
h	Pressure head	cm

The inverse relation h_from_theta(θ) is computed analytically from the van Genuchten equation.

3. Feddes Root Water Uptake Model

The actual root water uptake rate is:

S(z, h) = α(h) · Sₘₐₓ(z) · Tₚ

where:

• α(h) is the Feddes stress reduction factor (dimensionless, 0 to 1)
• Sₘₐₓ(z) is the normalized root length density distribution (dimensionless)
• Tₚ is the potential transpiration rate (cm/day)

Feddes stress function α(h):

α(h) = 0                    if h ≥ h₁  (too wet, anaerobic)
α(h) = (h - h₁)/(h₂ - h₁)  if h₂ ≤ h < h₁  (increasing uptake)
α(h) = 1                    if h₃ ≤ h < h₂  (optimal uptake)
α(h) = (h - h₄)/(h₃ - h₄)  if h₄ ≤ h < h₃  (decreasing uptake)
α(h) = 0                    if h < h₄  (too dry, permanent wilting)

Typical threshold values for many crops:

• h₁ = -10 cm (upper threshold — anaerobic zone)
• h₂ = -25 cm (lower threshold — full recovery)
• h₃ = -200 cm (upper stress threshold)
• h₄ = -8000 cm (permanent wilting point)

Root distribution Sₘₐₓ(z):

Two options are available:

• Linear: Sₘₐₓ(z) ∝ max(z_root - z, 0) — root density decreases linearly with depth
• Uniform: Sₘₐₓ(z) = constant for z ≤ z_root, zero otherwise

The distribution is normalized such that Σ Sₘₐₓ(z) · Δz = 1.

---

Numerical Method

The Richards equation is solved using a fully implicit finite-difference scheme with Picard iteration for nonlinearity.

• Spatial discretization: cell-centered grid with uniform spacing Δz
• Internode conductivity: arithmetic mean of adjacent nodes
• Time integration: backward Euler (fully implicit)
• Nonlinear solver: Picard iteration with relaxation factor ω = 0.8
• Convergence criterion: max(|h_new - h_old| / (|h_new| + ε)) < 1e-5
• Maximum Picard iterations: 50

Boundary conditions:

• Top: specified flux q_top (cm/day), positive into soil (irrigation, rainfall minus evaporation)
• Bottom: free drainage (unit gradient, q = K(h)) or fixed pressure head (Dirichlet)

---

Code Structure

RootWaterSim/
├── README.md                  # This file
├── LICENSE                    # MIT License
├── requirements.txt           # Python dependencies
├── .gitignore                 # Ignored files
├── .vscode/
│   └── settings.json          # VS Code settings
│
├── src/                       # Core library package
│   ├── __init__.py            # Package init (empty)
│   ├── soil.py                # VanGenuchtenSoil — soil hydraulic properties
│   ├── root_uptake.py         # feddes_stress, root_distribution — uptake model
│   └── solver.py              # RichardsSolver1D — PDE solver
│
├── app/                       # Streamlit web application
│   └── streamlit_app.py       # Interactive UI with Plotly visualisations
│
├── notebooks/                 # Jupyter notebooks for teaching
│   ├── 01_soil_hydraulic_properties.ipynb
│   └── 02_richards_equation_and_simulation.ipynb
│
├── tests/                     # Unit tests
│   ├── conftest.py            # Pytest configuration
│   ├── test_soil.py           # Tests for VanGenuchtenSoil
│   └── test_root_uptake.py    # Tests for feddes_stress and root_distribution
│
└── .pytest_cache/             # Pytest cache (auto-generated)

Module Reference

src/soil.py — VanGenuchtenSoil

Method	Description
theta(h)	Compute water content θ for pressure head h
K(h)	Compute unsaturated hydraulic conductivity K(h)
h_from_theta(theta)	Inverse of the water retention curve
capillary_capacity(h)	Compute specific moisture capacity C(h) = dθ/dh

Constructor parameters: theta_r, theta_s, alpha, n, Ks

src/root_uptake.py

Function	Description
feddes_stress(h, h1, h2, h3, h4)	Feddes water stress reduction factor α(h)
root_distribution(z, z_root, dz, method)	Normalized root length density distribution

src/solver.py — RichardsSolver1D

Method	Description
_picard_step(h_old, q_top, h_bottom)	One Picard iteration step
simulate(h_init, n_steps, q_top_series, h_bottom)	Run the full time-stepping simulation

Constructor parameters: soil, dz, z_nodes, dt, root_params

The root_params dictionary accepts:

• z_root: maximum rooting depth (cm)
• Tp: potential transpiration rate (cm/day)
• method: root distribution type ("linear" or "uniform")

---

Installation

1. Clone the repository:

   git clone https://github.com/GreenSmart-DSS/RootWaterSim.git
   cd RootWaterSim

2. (Optional) Create and activate a virtual environment:

   python -m venv .venv
   source venv/bin/activate      # Linux / macOS
   venv\Scripts\activate         # Windows

3. Install the dependencies:

   pip install -r requirements.txt

Dependencies

Package	Minimum Version	Purpose
numpy	≥1.21	Numerical arrays and operations
scipy	≥1.7	Sparse linear algebra (spsolve)
matplotlib	≥3.5	Plotting (notebooks)
plotly	≥5.0	Interactive visualisation (app)
streamlit	≥1.20	Web application framework
pytest	≥7.0	Unit testing
jupyter	≥1.0	Notebook environment
ipympl	≥0.9	Interactive matplotlib backend

---

Usage

Streamlit Web App

The easiest way to explore the model is the interactive web app:

streamlit run app/streamlit_app.py

Then open the URL shown in your terminal (usually http://localhost:8501).

Adjust soil type, rooting depth, irrigation schedule, and watch the moisture profile evolve. The app provides three tabs:

• Water Content Profile: θ(z) at any time step, with field capacity and wilting point markers
• Pressure Head: h(z) at any time step
• Plant Stress: stress factor α(h) and actual vs. potential transpiration over time

Jupyter Notebooks

jupyter notebook notebooks/

Start with 01_soil_hydraulic_properties.ipynb to learn about water retention curves, then move to 02_richards_equation_and_simulation.ipynb to see the full simulation pipeline.

Python API

import numpy as np
from src.soil import VanGenuchtenSoil
from src.solver import RichardsSolver1D
from src.root_uptake import feddes_stress, root_distribution

# Define soil
soil = VanGenuchtenSoil(
    theta_r=0.078, theta_s=0.43, alpha=0.036, n=1.56, Ks=24.96
)

# Grid and time settings
dz = 2.0          # cell size (cm)
z_max = 80.0      # soil depth (cm)
z_nodes = np.arange(0, z_max + dz, dz)
dt = 0.01         # time step (days)

# Root and plant parameters
root_params = {"z_root": 40.0, "Tp": 0.5, "method": "linear"}

# Create solver
solver = RichardsSolver1D(soil, dz, z_nodes, dt, root_params)

# Initial condition (hydrostatic-like)
h_init = np.linspace(-100, -50, len(z_nodes))

# Irrigation schedule: 1 cm/day for 5 days starting at day 1
n_steps = int(15 / dt)
q_top_series = np.zeros(n_steps)
irr_start = int(1 / dt)
irr_end = int(6 / dt)
q_top_series[irr_start:irr_end] = 1.0

# Run simulation
h_hist, theta_hist = solver.simulate(h_init, n_steps, q_top_series, h_bottom=None)

---

Testing

Run the test suite with pytest:

pytest tests/ -v

The test suite covers:

• test_soil.py (8 tests): water retention, conductivity, inversion, capillary capacity
• test_root_uptake.py (5 tests): Feddes stress function, root distribution (uniform and linear)

All 13 tests should pass.

---

Contributing

Contributions are welcome! Please follow these guidelines:

1. Fork the repository and create a feature branch.
2. Write unit tests for any new functionality.
3. Ensure all tests pass before submitting a pull request.
4. Keep the code modular and well-documented.

---

License

This project is licensed under the MIT License. See the LICENSE file for details.

Copyright (c) 2025 GreenSmart-DSS (Morteza Khoshsimaie Chenar)

---

References

1. 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.
2. Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3), 513–522.
3. Feddes, R. A., et al. (1978). A mathematical model of the water uptake by plants. Simulation monographs, 18.
4. Richards, L. A. (1931). Capillary conduction of liquids through porous mediums. Physics, 1(5), 318–333.