Design pass: the soil cascade needs no change

[dfalster]

theta_sat is already a barrier in the exact vector field, and error
control already rejects the overshooting steps. The excursions to
theta = 1.7 are explicit overshoot triggered by the step size a leg
inherits, not a water-balance defect.

Two faults instead: odelia accepts a NaN error estimate (odelia#52),
and a throw inside a stage evaluation preempts rejection entirely.
Capacity limiting is measured to cost stiffness for no benefit.

Design doc plus standalone scripts; no package code changes.
Refs #522, #599, #565.

🤖 Generated with Claude Code

[dfalster]

What this was convened to do, and why it concluded the opposite

Scoping #522 (rainfall pulses) produced a plan whose first stage was a redistribution fix in compute_rates: the cascade sets water_input[i] = K(θ_{i-1}) and never checks whether layer i can accept it, so a saturated layer appeared able to overfill the one below — and pulses looked unsafe until that was fixed. The gate before implementing was to design the scheme properly. It does not survive measurement.

Everything below is from notes/scripts/soil_cascade/ — a standalone R replica of compute_rates/soil_K_from_soil_theta plus odelia's exact OdeControl with plant's defaults. No plant build needed; Rscript exp1.R, exp2.R, exp4.R.

1. θ_sat is already a barrier

Inflow K(θ_{i-1}) and outflow K(θ_i) saturate at the same K_sat, so the rate is exactly zero at θ_sat, and the [0, θ_sat] clamp inside K stops it going positive above:

θ₁ (of θ_sat)	inflow	outflow	rate₁ (yr⁻¹)
0.050 (0.117)	163.041	0.000	+543.47
0.385 (0.900)	163.041	29.769	+444.24
0.424 (0.990)	163.041	138.627	+81.38
0.428 (1.000)	163.041	163.041	0.00

So θ cannot exceed saturation in exact arithmetic. One Euler step from θ₁ = 0.05: h = 1e-3 → 0.593; h = 0.08 → 43.5, reproducing the recorded θ = [-51.4, 52.0, …]. Pure explicit overshoot.

2. Error control already rejects those steps

h = 1e-3, 1e-2 and 0.08 all reject. From a cold start the worst case integrates in 15 steps with max θ = 0.428 exactly; ode_step_size_max barely matters (15 vs 21 at 0.05).

3. The trigger is the inherited step size

step_size_last persists across set_state_from_system() (stochastic_patch_runner.h:146-159), so a leg beginning after a discrete change starts at whatever the previous quiet leg grew to. Varying only that:

h at leg start	leaf throws	leaf returns NaN
1e-6 (fresh)	completes, 16 steps	completes, 16 steps
≥ 1e-3	dies at step 0	dies, NaN accepted

Two independent faults:

• odelia accepts a NaN error estimate — std::max(NaN, rmax) is NaN and NaN > 1.1 is false, so the accept branch runs. Filed as adjust_step_size accepts a NaN error estimate instead of rejecting the step odelia#52. Fixing it makes that path recoverable (14–15 steps, 9–12 rejections).
• A throw inside a stage evaluation preempts rejection — set_ode_state → check_finite_ode_state() runs at every RK stage, and phylloptim's collar root-find throws. The run dies having taken zero steps, and fixing the NaN does not help it at all.

4. Capacity limiting: measured, and rejected

scheme	max |upper-tri|	lower-tri?	max |λ|
baseline (donor-only)	0	yes	1.546e4
receiver rejects excess	0	yes	1.546e4
donor throttles	4.67e3	no	1.518e4

Donor-throttling — the literal reading of "a layer can't overfill the one below" — breaks lower-triangularity, which vindicates the original donor-only choice. Receiver-side rejection keeps it but costs 1.20× stiffness at 0.9 θ_sat and 1.45× at 0.99, and at 0.99 turns a +81 yr⁻¹ rate into −420 yr⁻¹ by discarding inflow the baseline retains — a real hydrological change needing a scientific_version bump and a scenario re-bless, to enforce a bound that already holds.

5. Relation to #565

Same criterion, independently: #565 uses donor control to triangularise the cross-patch Jacobian under the "design out stiffness rather than reach for implicit solvers" philosophy. This pass supplies measured backing for two things it asserts — the stiff diagonal (≈1.5e4 yr⁻¹, structure exactly triangular to machine precision) and the need for a spill term on large events. It also flags that its N = 15 bucket makes the §3 trigger worse: measured max |λ| at 0.99 θ_sat goes 1.77e4 (N=5, dz=0.3) → 5.30e4 (N=15, dz=0.1).

6. Consequence for #522

Pulses are no longer gated on a physics change. They need the event queue, the pulse action, and the R generator. The one physics-adjacent piece that is required is capping the jump at layer 0's free capacity, because a jump is applied outside the integrator where nothing can reject it — and a ~13 mm event already exceeds capacity from θ₀ = 0.40.

Caveat

These experiments isolate the soil column with a stub for physiology, so they establish the mechanism, not the 43%-of-seeds figure. The 40-seed harness should be re-run against real fixes before #599 closes.

Retractions from my own earlier notes on this

• I claimed on [TF24 hydraulics] [Env drivers] Enable event-driven hydraulics #522 that a general events queue would need care around "the fixed-step explicit solver". There is no fixed-step solver; that issue's body is now corrected.
• I claimed RODAS was unreachable without a multi-week scalar-templating project. True of RodasStep as gated, but odelia#42's BlockFdJacobian needs no rebind — see the addendum on Dense TF24 stochastic runs throw on an unattainable psi_soil (43% of seeds on develop) #599.