Let soil layers have different widths

[dfalster]

Soil evaporation needs a very thin layer at the top of the profile, which a layer
count cannot express. set_soil_layer_widths() takes the WIDTH of each layer, so
nothing is left ambiguous between a top, a midpoint and a bottom; boundaries and
midpoints are derived from the widths. depth becomes read-only, because writing
it rebuilt no geometry and silently left the old column in place.

Needs phylloptim 0.9.0 (traitecoevo/phylloptim#134), whose per-layer root
resistance is a fix and not only an enabler: one column average inflated total
vertical root resistance 3.68x on a 2 cm surface layer.

No default output moves — run_scm on the default environment is bit-identical.
Soil evaporation itself is out of scope.

Closes #626.

[dfalster]

Layers or nodes? The framing question on the issue

It was put that plant "doesn't really have layers, but nodes at depth". Half true — and the half that is false is the structural half.

TF24_Environment was already a finite-volume layer discretisation carrying all three geometric descriptions, and each is read by exactly the part of the model that needs it:

member	meaning	read by
dz[i]	layer width — already a vector	the water balance; the pulse capacity cap
z[i]	cumulative depth to the bottom of layer i; z.back() == depth	root-mass distribution, Q(z) differenced across a layer
z_mid[i]	layer midpoint	the leaf's gravitational head

So the values were not "none of top/midpoint/bottom": z[i] is unambiguously the bottom, and both packages already treated it that way. The root-mass distribution is layer-integrated (Q(z[a-1]) − Q(z[a]) with prev_q seeded at 1) precisely because those are boundaries. And the state is θ, intensive, so it needs no width at all.

What was nodal was the vocabulary: the setter said set_soil_number_of_depths, R saw a vector called soil_depth, delta_z was commented "distance between layers" though used as a width, and R could not ask for the geometry at all. So the fix is definitional rather than structural — width is the only geometry a caller states, and nothing is left to be a top, a midpoint or a bottom.

The inter-layer cascade needs no change, and this is the note asking you not to "fix" it later

water_flux[i] = K(theta_i), fed to layer i+1, is free drainage: Darcy's q = −K(∂ψ/∂z + 1) with the matric-gradient term dropped. Two consequences:

• it contains no internode distance, so there is nothing to rescale when the layers stop being equal, and it is evaluated at the donor cell, which is correct upwinding;
• the units work at any width — q is a flux density (m yr⁻¹), rate = (q_in − q_out − uptake)/dz[i] is yr⁻¹, and sum(dz[i]·dθ_i/dt) = infiltration − drainage − uptake is exact for any distribution.

notes/plan-tf24-soil-redistribution.md reached the same conclusion from the redistribution side (#608), listing "any change to the inter-layer cascade" under not needed.

⚠️ But if anyone adds a capillary or diffusive flux — and soil evaporation will make them want to, since a top layer drying from above needs upward flow — that term's distance is z_mid[i+1] − z_mid[i], NOT dz[i], and the interfacial conductivity must be a distance-weighted harmonic or geometric mean of the two cells. The two distances coincide only for equal layers, so a uniform-profile budget test passes on the mistake. The graded-profile budget test added here does not.

Two floating-point constraints, both load-bearing

Measured over 180 (depth, n) pairs:

route to a uniform width vector	disagrees with today's values
pass the environment's own dz	0 / 180
difference the profile, diff(c(0, z))	122 / 180
rebuild z as cumsum(dz)	103 / 180

Hence two rules, both stated at their call sites:

1. TF24_Strategy passes environment.dz, never a differenced z. phylloptim::layer_thicknesses() is for a caller holding only a profile; from here it is actively wrong, and resistance scales with dz[i]². agents.md §3.2 updated — its old instruction ("do not open-code soil_depths_.back()/n") is obsolete and said the opposite of what is now required.
2. set_soil_number_of_depths keeps z[i] = (i+1)·delta_z and is deliberately NOT unified with the new setter's running sum. z feeds Q() and the gravitational head, so unifying them re-baselines every TF24 number for nothing. The geometry test compares the two paths with expect_equal, not expect_identical, and says why.

Bit-exactness, and the plant-side proof the upstream fix landed

Built both ways — per-layer, and with the scalar temporarily restored — in the same tree:

run_scm on the default uniform environment:

build	offspring_production
scalar dz	11.143684537309916
per-layer dz[i]	11.143684537309916

Bit-identical, which is why no default output moves and scientific_version is unchanged. On a graded column the same run gives 9.2796204566000107 (scalar) against 12.981712222339016 (per-layer).

At the individual level (θ = 0.25, height 5 m, 1.5 m column):

profile	opt_root_psi (MPa) scalar → per-layer	assimilation scalar → per-layer
uniform, 5 layers	1.661740 → 1.661740	15.464934 → 15.464934
2 cm surface layer	1.747280 → 0.598266	15.338550 → 17.251842
thick over thin	1.654759 → 4.211291	15.478450 → 11.179142

Read the middle row: a scalar thickness charges the 2 cm layer for 0.3 m of root segment — 225× over-resistance in that layer, 3.7× in total — so the plant appears to need more suction when it needs far less. Reversing the profile reverses the error. So test-tf24-root-pars.R diffs the sign, not a tolerance: a scalar keeps every graded profile inside 6% of the uniform value and on the wrong side of it, which no tolerance test could see.

The stiffness and runoff limits, which this PR does not fix

Shipped permissive by choice: widths are validated for finiteness and positivity, nothing more. The numbers, so they are on the record:

min dz	diagonal K'(θ)/dz at 0.99·θ_sat	stable step 2/λ
0.30 m (default)	1.8e4 yr⁻¹	1.1e-4 yr
0.10 m	5.3e4	3.8e-5
0.02 m	2.6e5	7.6e-6 yr ≈ 4 min
0.01 m	5.3e5	3.8e-6

The analytic scaling reproduces the 1.77e4 / 5.30e4 measured in plan-tf24-soil-redistribution.md. It only bites when the layer is wet — at dz = 0.02 the diagonal is 2.6e5 at 0.99·θ_sat but 8.5 at 0.5·θ_sat — so the hazard is the post-rain transient, which is exactly when an evaporation layer matters. And ode_step_size_min is 1e-6 yr, which the stable step reaches at dz ≈ 2.6 mm: a hard floor on layer thickness in this model.

Separately, runoff is recalibrated by a thin layer 0. Infiltration keys off that layer alone through (θ₀/θ_sat)^b_infil with b_infil = 8, so a thin layer fills far faster and holds the switch shut longer. In the R replica (notes/scripts/soil_cascade/), 20 events × 50 mm: annual runoff fraction 0.134 → 0.204 from thinning layer 0 alone, no parameter changed. Inert under constant rain — the steady state contains no dz — so the whole effect is in the event transient. a_infil/b_infil were calibrated against a 30 cm surface layer. This is the follow-up item most likely to be missed, and it is a consequence of the new capability rather than a regression.

Both are in soil_widths_graded()'s help, where a user deliberately asks for a thin layer, and in notes/plan-626-soil-layer-widths.md.

⚠️ The scenario gateway is red, and it was red before this branch

Every scenario's offspring_production differs from tests/testthat/test_data/scenario_baseline.rds:

	S01	S02	S03	S04	S05	S06	S07	S08
rel. change	−0.68	−0.88	+29.3	+26.8	+40.9	+19.1	−0.43	−0.37

19–41× on S03–S06 is the range #590's density-coordinate change is documented to produce, so the recorded baseline predates it. The bit-identical run_scm above is what says #626 did not move it — every scenario runs a uniform column.

Deliberately not re-blessed here. Filed separately so it stops costing reviewers an A/B; folding an unrelated re-baselining into this PR would hide what actually moved.

Verification

• make clean && make RcppR6 && make full_compile; full test_dir green, no failures.
• test-model-version.R (the scientific-surface hash) green under NOT_CRAN=true.
• Water budget closes to < 1e-12 relative on four graded profiles, including a single 2 cm layer.
• test-patch.R's pinned n+5 ODE layout and first-step rates unmoved, as expected — the default environment is untouched.

Release order

⚠️ phylloptim#134 must merge first, and its merge sha must replace the branch ref currently in Remotes:. The pin is == and plant compiles against the installed headers, so install 0.9.0 locally before bumping the pin — doing it the other way round fails in a way that reads like a broken build (it did here).

⚠️ Stacked on #632, which will squash-merge. Afterwards: git rebase --onto develop feature/events-522 feature/soil-layer-widths-626, and check the changed-file count is ~19 and not ~55 before merging.

[dfalster]

⚠️ No CI runs on this PR, and "no checks reported" is not a pass

R-CMD-check.yaml's pull_request trigger is restricted to base branches main, master and develop (.github/workflows/R-CMD-check.yaml:12-16). This PR's base is feature/events-522, so no workflow fires at all — GitHub shows no checks rather than a failure, which reads the same as green if you are skimming.

So the plant half is verified locally only, until it is rebased onto develop after #632 lands. What was run here:

• make clean && make RcppR6 && make full_compile, then the full test_dir — no failures.
• test-model-version.R (the scientific-surface hash) green under NOT_CRAN=true.
• The bit-identical run_scm A/B in the comment above, both arms built in the same tree.

The rebase is needed anyway (#632 squash-merges), so the sequence is: merge phylloptim#134 → swap Remotes: from the branch ref to the merge sha → merge #632 → git rebase --onto develop feature/events-522 → CI runs for the first time here.

phylloptim#134 is green, and the cross-platform figures did not move

All six jobs pass. The golden comparison's summary line is identical on gcc and clang and identical to the figures recorded in that repo's own guide:

3626 of 5184 values differ. Worst by class:
  profit (the maximum)     rel 1.86e-07  (tol 1e-05)
                           abs 1.35e-09
  others (from the argmax) rel 9.02e-06  (tol 0.005)

Same worst rows, same counts. So the per-layer change moved neither the macOS bit-exact baseline nor the cross-platform spread — which is what "bit-identical for every profile reachable before this change" has to mean in practice.

[dfalster]

Upstream merged; pin swapped to the commit

phylloptim#134 merged as 232a81ee3f9e3c7c6be620cef8611b5be09a8409, and Remotes: now points at that sha rather than the branch ref it carried during review. Rebuilt from clean against the merged headers and re-verified: run_scm on the default environment gives 11.143684537309916, identical to the value measured against the branch build, and the full suite is green.

This PR passes #640's stricter gateway

#640 re-blessed the baseline and added a numeric offspring_production diff at 1e-3 relative — which is a materially harder test than the two thresholded booleans this branch was measured against. Checked ahead of the rebase by borrowing #640's test-scenario-gateway.R and scenario_baseline.rds onto this branch without committing them:

[ FAIL 0 | WARN 0 | SKIP 0 | PASS 2 ]

Worth noting why it passes: the numbers #640 blessed on develop are, to every printed digit, the ones this branch produces —

	S01	S02	S03	S07
blessed on develop	4.797763201e-01	1.153314819e-02	1.665416219e-10	8.107050184e+02

So neither #626 nor #632 moves a scenario, which is the expected result for #626 (every scenario runs a uniform column, and the default column is bit-identical) and a useful free datapoint for #632.

Those two files were reverted rather than committed — they belong to develop and arrive via the rebase.

Correction to #639

I attributed the drift there to the density-coordinate change (#590), flagged as unconfirmed and with bisection as its first action item. #640 did the bisection and the answer is #633 and #619, not #590. The magnitude range fitted #590's documented 2.4–47x, but fitting a range is not identifying a cause.

Still blocked on #632

R-CMD-check.yaml only fires for base main/master/develop, so this PR still has no CI. It cannot be rebased onto develop early either — it genuinely depends on #632: test-events.R's width fixes are edits to code #632 adds, and aux_num = 5 is #632's. Sequence unchanged from above: merge #632 → git rebase --onto develop feature/events-522 → file count should read ~19, not ~55 → CI runs here for the first time.

[dfalster]

Rebased onto develop; #632 is in, so this is now a plain PR rather than a stacked one

git rebase --onto origin/develop a5b221c1 — the two commits replay onto develop's tip (32151e87). 19 files, 980 insertions, the same count as before the rebase, which is the check that matters: a squash-merged parent otherwise makes the child re-land the whole parent diff, and that would have read ~55 files here.

One conflict, in NEWS.md, where #636's breaking-change entry and mine both sit at the top of the section. Both kept, #636's first. Nothing else conflicted — tests/testthat/test-events.R and agents.md auto-merged.

Everything #632 added now arrives from develop rather than from this branch (aux_num = 5, add_water_pulse_to_layer), and test-scenario-gateway.R / scenario_baseline.rds are untouched by this PR — they come from #640.

Re-verified after the rebase, not carried over

Rebuilt from clean against merged phylloptim 0.9.0 and re-run on the rebased tree, so none of this is inherited from the pre-rebase measurements:

• Full test_dir: no failures.
• test-model-version.R (scientific-surface hash): 23 pass.
• test-scenario-gateway.R with Re-bless the scenario gateway and diff its numbers #640's numeric 1e-3 diff: 2 pass. This is the first time it has been run for real on this branch rather than by borrowing the files.
• run_scm on the default environment: 11.143684537309916 — the same value measured before the rebase, and before the phylloptim bump. Picking up Add a general discrete-event queue #632, Refuse an unknown trait name #637 and Re-bless the scenario gateway and diff its numbers #640 did not move it either.

A note for whoever picks this up next

The local branch label feature/soil-layer-widths-626 in the shared working copy points at 99902446 ("Name a trait each strategy actually carries"), which is unrelated to #626 — a second session reset it. The PR's commits only ever existed on origin, and this rebase was done from the remote ref in a separate worktree so the shared checkout was not disturbed. ⚠️ If that local label gets pushed, it will clobber this PR. Either delete it or reset it to origin/feature/soil-layer-widths-626 before working here.

Also worth knowing: develop still pins phylloptim == 0.8.0, so the site R library holds 0.8.0 and any sibling checkout reinstalling it will break a build of this branch. This branch was built against 0.9.0 in a private library rather than by reinstalling over the shared one, so other sessions on develop were left working.

[itowers1]

Wondering whether there would be an issue with a more narrow layer being beneath a wider layer. Small loss of water from above could overwhelm the water capacity of the layer below. Also, can't necessarily see a case where that would be necessary anyway. Maybe we should apply a block, enforcing small -> large?

[dfalster]

@itowers1 — measured it. The answer splits in two: the water-capacity worry does not materialise, but you are right that something breaks, and it is not what either of us expected.

Notation: a profile like 0.74 / 0.74 / 0.02 is a three-layer soil column, top down — two layers 74 cm thick over one 2 cm thick, 1.5 m in total. θ is volumetric soil moisture (m³ water per m³ soil) and θ_sat is its saturated value, so θ/θ_sat = 1 means the layer's pore space is completely full.

The water side is fine, and that part is clean

A thin layer's storage is not a bottleneck, because the water arriving from above does not depend on how thick the layer receiving it is. Drainage between layers is q = K(θ), a rate of water arrival per unit ground area (m yr⁻¹) — it has no thickness in it. The receiving layer settles where its own drainage matches what arrives, q_in = K(θ), and that balance point contains no thickness either. Since the layer above can never deliver more than its own saturated conductivity, the layer below settles at or under saturation whatever its thickness. Thickness changes only how fast a layer fills, not how full it ends up.

Measured, starting the whole column at 99% of saturation under 2 m yr⁻¹ rain — a saturated thick layer draining hard into a thin one, which is the case you describe:

profile	thinnest layer	how fast that layer responds	fullest any layer got	water-balance error
uniform 0.3 ×5	0.300 m	1.8e4 yr⁻¹	0.726 of saturation	−1.4e-15
thin TOP (2 cm)	0.020 m	2.6e5 yr⁻¹	0.726	8.4e-16
thin BOTTOM (2 cm)	0.020 m	2.6e5 yr⁻¹	0.726	1.8e-16
thin MIDDLE (2 cm)	0.020 m	2.6e5 yr⁻¹	0.726	5.2e-16
thin top and bottom	0.020 m	2.6e5 yr⁻¹	0.726	1.8e-15

No layer gets past 73% of saturation in any of them, and water is conserved to rounding. The middle column is the rate at which the thinnest layer's moisture responds to a change — it is what sets how short a step the solver has to take, and it comes out identical whether the thin layer sits at the top, the middle or the bottom. It depends only on that layer's own thickness, so which end of the column it sits at cannot enter it.

But something does break, and it is not the soil

0.74 / 0.74 / 0.02 does fail. Not intermittently — every single run, at every rainfall rate I tried from 1.5 to 6 m yr⁻¹, and at the same time each run (t = 4.37 yr). So it is a real property of that configuration rather than an occasional numerical flake. But the thing that fails is the plant population, not the water:

Non-finite cohort density in the SCM size-density (characteristic) equations:
species 1 has a node with density=4.03e52 (log_density=121.1, height=0.859) at time=4.37

That is the guard #552 added for #550, whose own comment says "the divergence is in the equations, not the stepper" — the number of plants of a given size runs away to infinity, while the plants themselves stay perfectly ordinary (that node is a healthy 0.86 m tall).

And the soil is genuinely fine underneath it. Running the same configuration to t = 4, which is 0.37 yr before it dies:

	fullest any layer got	driest any layer got	water-balance error
0.74/0.74/0.02	0.990 of saturation — exactly where it started, never exceeded	0.323	3.1e-15
uniform 0.3 ×5	0.990	0.299	0.0

So no layer is ever overwhelmed, nothing oversaturates, and water is conserved. The water balance is not what fails.

Where that leaves your proposal — including the part that favours it

Two findings pull opposite ways, so here are both.

For the block. Every failure I can find has a thin layer sitting under a thicker one, and the mirror image is clean. Scanning the bottom layer's thickness with two equal layers above it, total held at 1.5 m — then the same scan with the thin layer moved to the top instead:

bottom layer	thin at the bottom	same thickness, thin at the top
30 cm → 6 cm	works	works
5 cm	FAILS	works
4 cm	works	works
3, 2.5, 2, 1.5 cm	FAILS	works
1 cm, 5 mm	works	works

None of 18 profiles that get thicker with depth failed, including every soil_widths_graded() ladder down to a 5 mm top layer. So your rule would have blocked every failure I found, and it costs the motivating use case nothing — a thin surface layer over thicker ones already satisfies it.

Against it. Three things:

1. The boundary is ragged, so "decreasing" is not really the rule. Look down that first column: 5 cm fails, 4 cm works, 3 through 1.5 cm fail, 1 cm and 5 mm work again. A failure set interleaved like that is a knife-edge in a sensitive equation, not a property of the geometry. So thinning-with-depth is something the failures I happened to sample have in common — not an explanation of them. I have not found the mechanism and do not want to imply otherwise.
2. It would refuse configurations that work. 1.20/0.28/0.02, 1.45/0.045/0.005 (a 5 mm bottom layer), 0.70/0.70/0.10, 0.72/0.72/0.06 and 0.50/0.40/0.30/0.28/0.02 all run fine and would all be blocked.
3. It fixes a symptom in a different part of the model. The fragility is [TF24 hydraulics] SCM cohort-density blow-up under extreme seasonal drought #550's, in the population equations. A geometry rule in the soil setter would close one route into it and leave the others open. And the failure is already loud — fix: graceful failure on SCM cohort-density blow-up (#550) #552's guard stops the run with a clear message, so nobody gets quietly wrong numbers, only a refused run.

Your second point is the one I cannot settle from the code: whether a thin layer under a thick one is ever actually wanted. A thin bottom layer at a water table or a bedrock contact is a real land-surface-model pattern, and land models often thin at both ends — but I have no plant use case in hand, and if there is none, a block costs nothing and removes a foot-gun.

My recommendation: file the reproducer against #550 rather than restricting the geometry, because the fragility lives there, it fails loudly, and the ragged boundary means a geometry rule would give false reassurance rather than real safety. But this is a modelling call between the two of you, and the evidence for your instinct is stronger than I expected — so if you would rather close the door, say so and I will add it. It is a one-line check in set_soil_layer_widths, and either way soil_widths_graded() and every test here are unaffected.

Everything above is measured on the rebased branch (now based on develop, CI green).

[dfalster]

Done — set_soil_layer_widths() now refuses a profile in which any layer is thinner than the one above it, with equal allowed, so rep(w, n) and every uniform profile stay legal. Commit 3602a58.

@itowers1 you asked whether there is a case for thin-under-thick. I went looking rather than guessing, and the answer supports you — with one caveat that is worth having in writing.

The literature is one-directional

Refining towards the surface is the standard, and I found nothing advocating the reverse:

• Downer & Ogden 2004 put it most directly: cells "of the order of 1 cm are required near the soil surface, but not throughout the soil column".
• SiBUC (Tinumbang et al. 2024) uses a 2 cm first layer over a second "several meters thick" — our motivating shape, at a more extreme ratio than anything here.
• FATES-HYDRO (Fang et al. 2022) is the closest analogue to plant — plant hydrodynamics inside a vegetation-dynamics model. Coarsening the near-surface grid beyond 10 cm "significantly overestimates aboveground biomass" in most temperate sites, and they warn specifically against a coarse surface grid "when soil evaporation is a dominant fraction of evapotranspiration". That is exactly Allow soil depths to have variable depths #626's motivation, from the other end.
• CLM3, whose conductivity function TF24 actually uses (Zeng & Decker 2009 is the source of our soil_K_from_soil_theta), grades thin-to-thick over 10 layers.

And the two cases one might invoke for refining at depth are served by other mechanisms, not by geometry. A contrast in soil properties at a horizon boundary is per-layer parameters — which we already have from #558 — and is what the layered-soil literature (Berninger et al. 2015, Aryeni et al. 2022) is actually about. A water table is a bottom boundary condition (He et al. 2022, He et al. 2026), and TF24 does not have one. So the geometry was never the right place to express either.

⚠️ The caveat, which is now in the code and NEWS

A thickening profile is not automatically well resolved, and this rule must not be read as saying it is. Mackay et al. 2022 show the usual Richards discretisation carries an interpolation error that "vanish[es] on uniform grids" but "progressively worsens with increasing layer thickness differences" — i.e. with the contrast between adjacent layers, whichever direction it runs. They explicitly flag models where "users can specify soil layer thicknesses" freely, which is what #626 makes us. And Regenass et al. 2021 add that "a fast increase in layer thickness with depth is equivalent to a lower vertical resolution".

So grading steeply has its own cost at the bottom of the column: soil_widths_graded(1.5, 5, top = 0.005) buys a 5 mm surface layer at the price of a 1.11 m basal one. That is now in that function's help, with the advice to prefer more layers over a steeper ladder.

What the block costs, measured

You asked implicitly and I should answer it explicitly: it refuses a lot that works. Sampling 48 profiles that thin with depth (3–6 layers, 1.5 m total, widths ≥ 5 mm):

	count
run fine — refused for nothing	46 (96%)
actually fail	2 (4%)

Examples now refused that ran clean: 0.707/0.462/0.331, 0.677/0.670/0.153, 0.633/0.444/0.423, 0.981/0.324/0.195.

I have gone with the block anyway, on your reasoning rather than against it: the 96% is only an argument for permitting these if someone wants them, and on the evidence above nobody should. But it means the rule is a foot-gun removed, not a stability proof — the failures it blocks are ragged in the geometry (5 cm fails, 4 cm passes, 3–1.5 cm fail, 1 cm passes again), so thinning-with-depth is a property the observed failures share rather than a characterisation of them. Both that and the #550 pointer are in the comment at the guard, so the next person does not mistake it for a guarantee.

Verification

• Full suite: 3196 pass, 0 fail.
• run_scm on the default environment: 11.143684537309916 — unchanged by the guard, as it must be.
• Scenario gateway (Re-bless the scenario gateway and diff its numbers #640's numeric 1e-3 diff): passes.
• New coverage in test-environment-TF24.R: the refusal, that it names both offending layers and their widths, that non-monotone counts too, that equal widths and every thickening profile stay legal, and that a refusal leaves the previous geometry intact.
• Two of my own tests used thick-over-thin profiles and had to change. The one that mattered was the sign test in test-tf24-root-pars.R: it asserted the per-layer fix in both directions, and the thick-over-thin half is now unreachable. The thin-top half still discriminates on its own — a scalar thickness puts opt_root_psi at 1.747 MPa where per-layer gives 0.598, against 1.662 for uniform, so the sign of the difference still flips — and I have added a steeply-graded profile as a second increasing case. The retired numbers are kept in the comment.

One thing for someone else

make roxygen regenerated man/rainfall_pulse.Rd and man/run_scm.Rd, which are stale on develop: rainfall_pulse() gained an env argument and run_scm() reordered its arguments in the R source without a doc regen. Not included here — not my change, and it would muddy this diff — but they want regenerating on develop.