Why lumped capacity matters in unsaturated lattice transport

OOFEM's van Genuchten implementation works with both the consistent and the lumped capacity matrix. On a wetting-front problem only one of them avoids the textbook overshoot above the initial suction.

The previous post showed the transport lattice reproducing the 1D diffusion analytical for a linear constant-capacity problem. Real porous media are rarely that easy to model: in concrete, soil, and almost any partially-saturated material the capacity c(p) = dθ/dp is a sharp peak near the air-entry value, and the relative permeability k_r(p) varies by orders of magnitude with saturation. This is the Richards-equation regime, and it’s where the choice between a consistent and a lumped capacity matrix is important.

Two ways to assemble the capacity

For a transport conduit between nodal pressures p₁, p₂, capacity at the integration point c, length L (along the Voronoi edge), and cross-section area A (the perpendicular Delaunay edge times the out-of-plane thickness), both forms preserve the total element mass c · L · A but distribute it differently. The consistent form keeps the off-diagonal couplings:

M_consistent = (c L A / 6) · [[2, 1], [1, 2]]

The lumped form moves the off-diagonal mass onto the diagonal:

M_lumped     = (c L A / 2) · [[1, 0], [0, 1]]

For a homogeneous linear-c problem the two give nearly identical results — the difference washes out over a step. Under van Genuchten, where c(p) spikes at the wetting front, the off-diagonal coupling amplifies high-frequency transients exactly where the nonlinearity is the strongest. This was first analysed in detail by Celia, Bouloutas & Zarba (1990), who showed that consistent-mass Richards solvers lose mass and oscillate near sharp fronts while their lumped counterparts stay monotone.

A small Celia-style benchmark

A 100×100 mm prism, van Genuchten with mild parameters (m = 0.5, a = 1000 Pa, θ_s = 1, θ_r = 0), initially partially unsaturated (u_IC = 2000 Pa, S ≈ 0.45). The bottom face is saturated (u_BC = 0) from t = 0 — a step BC, no ramp. Sign convention here is tension-positive: u = 0 is fully saturated and larger u is drier (opposite to Celia’s geotechnical convention, where the same physical state has negative pressure head). Adaptive time-stepping via deltatfunction (very small dt for the first 50 steps to absorb the initial Jacobian stiffness, then growing) lets the solver converge for both capacity choices.

Both runs use the same mesh and the same modified Newton solver — the only difference is the #@lumpedcapacity directive in control.in.

The overshoot in profile

The blue solid line is the lumped envelope (maximum suction in each horizontal band); the red dashed line is the consistent envelope. The grey dashed reference at the top is the initial-condition bound, p = 2000 Pa.

Around step 50, where the wetting front is roughly 4 mm into the specimen, the consistent capacity overshoots the IC bound by about 17% — peaks of p ≈ 2350 Pa appear in a thin band right at the front. Lumped stays strictly below p = 2000 everywhere. The overshoot dissipates by step 200 as the front advances and the high-frequency content gets damped.

The pressure in 2D

The 2D pressure fields of the two settings look very similar.

Lumped on the left, consistent on the right.

Why is there overshoot but not undershoot?

Two things combine. First, van Genuchten’s k_r(p) is highly asymmetric: on the wet side (low p, near BC) the relative permeability is near unity, so high-frequency oscillations get damped almost instantly by the strong diffusion. On the dry side (high p, near IC) k_r ≈ 0.7% — the same oscillations linger. Second, the BC nodes are clamped exactly at p = 0, so the undershoot mode would have to develop one cell in from the boundary and is immediately erased by the wet-side damping. The visible artefact is what’s left: a small wave overshooting above the IC.

For a linear constant-c problem with sharp BC step you would see both overshoot and undershoot symmetrically. The asymmetry here is purely numerical, driven by the lattice + the vG nonlinearity.

What to take away

For Richards-type problems with sharp wetting fronts, the lumped capacity is essentially required. The consistent matrix isn’t wrong in the sense of accuracy at converged equilibrium — both eventually agree — but the transient artefact at the front breaks the physical bound and makes the result harder to trust. For concrete-like parameters with much smaller a and paev, the consistent variant sometimes fails to converge at all and the choice isn’t even available.

Reproduce

git clone https://github.com/githubgrasp/oofem-examples.git
cd oofem-examples/lattice-flow-2d-lumped
docker run --rm -v "$PWD":/work ghcr.io/githubgrasp/oofem-public:lattice-flow-2d-lumped bash run-all.sh

The :lattice-flow-2d-lumped image tag is immutable — use :latest instead to track the current OOFEM build.

The example folder is at github.com/githubgrasp/oofem-examples/lattice-flow-2d-lumped; issues and questions go on the issue tracker.

Reference

Celia, M. A., Bouloutas, E. T. & Zarba, R. L. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 26(7), 1483–1496. doi:10.1029/WR026i007p01483

Built with OOFEM.

← All posts