Next steps — the modelling roadmap

Distinct from data_gaps.md, which lists data we cannot fetch. This lists work to do, in the order it has to be done, and says what each step unblocks. Written 2026-08-11.

The ordering is not a preference. Each step consumes the output of the one before it, so doing them out of order means fitting a model to inputs that are about to change.

1. Dynamic water modelling — rainfall and river level coupling

Do this first. Everything downstream needs a driver, and rainfall is the only true forcing in the system. Until the catchment responds to rain in time, every water balance we have is an annual average wearing a daily timestamp.

What exists already

• build_rainfall.py builds per-catchment daily rainfall series for the nine gauged catchments plus the ungauged remainder
• daily flow series for the gauges; calibrate_runoff.py already fits a soil-moisture bucket against them
• analyse_recession.py for the falling limb

What is missing — the coupling itself: rainfall → soil moisture → runoff → river level, at a daily step, calibrated per catchment rather than regionally. The bucket exists; what it lacks is a routing and level relation so that “river level on day n” is a modelled quantity and not an observation we happen to hold.

Done looks like: a daily river level at each intake candidate, driven by rainfall, reproducing the observed hydrograph within a stated error over the calibration period — and, critically, over a validation period it was not fitted to.

Why it must be first: step 2 needs a driver and step 3 needs river level to decide what is divertible and what a returned flow does.

It is now a HARD prerequisite, not an ordering preference

Step 2 was built and it hit a wall that only step 1 can move. The coupled model has no quick-runoff pathway - soil-bucket overflow all becomes recharge - so it identifies recharge / Sy and neither term alone. A three-arm sensitivity test confirmed the degeneracy: moving the AET prior by +/-90 mm/yr moved recharge by 15% and fitted Sy by 13%, conserving the ratio to within 2.5%.

Step 1 partitions overflow between runoff and recharge. That is the same split the river model needs, so it closes the degeneracy as a BY-PRODUCT rather than as extra work.

Rainfall is already areal, and the routing already has a scar

• rain_catchment_daily.csv is area-weighted per catchment through rain_thiessen_weights.csv, over the catchments delineated by build_catchments.py. Nothing point-based needs fixing.
• Water must be counted exactly once. compare_04c_nesting.py exists because run 04 summed nine gauges whose catchments NEST, so shared water was double-counted; 04c is the eight-gauge correction. That is the precedent - the routing test is not “does rain reach a river” but “does each parcel reach exactly one accounting point”.

Flood / high-flow data is the calibration lever

Mean flows are where the quick and slow pathways are LEAST distinguishable, and that is what we have been calibrating on. Flood peaks are where they separate maximally: a sharp peak with a muted groundwater response is direct evidence of the split; a muted peak with a long recession is evidence of the other.

fetch_river_levels.py pulls stage from the hydrology API (the flood-monitoring service keeps only ~4 weeks). One dataset does two jobs - it also turns h_base in step 2 from a constant into a series.

Carry this caution into calibration: flood peaks are often SUB-DAILY. A daily model under-fits them systematically, and tuning the runoff fraction to compensate absorbs a timestep error into a physical parameter. Check period per measure first.

Also worth having: EA Recorded Flood Outlines. They constrain WHERE water goes once channels are exceeded, which matters here in a way it would not for a pure water-resources scheme - this project’s premise is that the floodplain is storage.

Sub-catchments at key points — worth testing, with one rule

The idea: subdivide strategically where flood data is good, so that river level at points we care about becomes a modelled output rather than an interpolation. This is semi-distributed modelling and it is sound. One rule keeps it honest:

Subdividing where an OBSERVATION exists adds information. Subdividing where a DECISION needs an answer only relocates the output. Both are legitimate and they are not the same thing. A model cut into twenty sub-catchments with six gauges has six constraints and fourteen extrapolations, and will look more precise while being no better informed - the classic way a semi-distributed model becomes confidently wrong.

So: split at gauges and good flood records FIRST and check the fit improves; only then push sub-catchment outlets out to candidate intake sites, and label those outputs as extrapolated. The nesting scar above is the same lesson from the other direction.

2. Dynamic water table

What exists

• build_water_table.py — a STATIC surface, wt = 0.05 + 0.768 × ground, n=100, R² 0.83, residual SD 21.1 m
• build_wt_dynamics.py — the observed dynamics characterised but not modelled: Chalk swings a median 11.9 m/yr, most stations peak in Jan–Feb, Chalk headroom 52.8 m at the seasonal high

What is missing — a water table that MOVES, driven by step 1. The current surface has no time dimension, so saturated thickness, and therefore every volume in aquifer_volumes.md, is a single annual figure for a quantity observed to swing by 12 m.

Method, when we get to it: cross-correlate each station against its catchment rainfall series and its nearest gauge; fit a lag and a gain per aquifer; drive the fitted surface with step 1’s output. The lag is the number that matters — the project’s surplus-to-deficit lag is 229 days, and an aquifer whose head peaks with the rivers offers no seasonal shift at all.

Reuse the rainfall→river pattern — with one change

The obvious move, and the right one: the water table is another STORE with an input, a loss and a lag, which is exactly what calibrate_runoff.py already fits for soil moisture. Same daily step, same calibrate-then-validate-on-held-out-years discipline, same shape of object. Do not invent a second framework for it.

Two things that do NOT carry over, and both matter:

1. The time constant is an order of magnitude longer. Runoff responds in days; the water table in months. Consequences: the calibration window must contain several full recharge cycles (the ~11-year record from 2010 is enough, but a 3-year station is not — see the coverage note above), and a fit that looks good on a single wet winter is meaningless.

2. THE COUPLING IS TWO-WAY, AND RAINFALL→RIVER IS NOT. Rain drives rivers and rivers do not drive rain. But the aquifer feeds the river as baseflow and a high river stage can recharge the aquifer through the bed. Modelling it one-way — river level as a driver of the water table — would silently double-count: the same water arrives as baseflow in the river AND is credited as storage in the aquifer.

This is the single place where copying the rainfall pattern wholesale would introduce an error rather than save work. The fix is to make exchange a function of the HEAD DIFFERENCE between aquifer and river, signed, so the direction falls out of the state rather than being assumed — which also makes step 3’s closed balance possible instead of retro-fitted.

Two things to carry forward, both already established

• Headroom, not storage, is the binding constraint on injection. Surplus water exists in winter; winter is when the aquifer is fullest. A dynamic water table is what turns that from an observation into a constraint the model enforces.
• Coverage is the weak point, and it is specific: only 47 of 115 stations have 3+ years, and the Otter Sandstone — the actual target — has one. Fix that before trusting a fitted response for HEY.

3. Injection and baseflow accountancy for leaky aquifers

Last, because it needs both of the above. This is the step that closes the loop: what happens to water we put in, and what comes back out where we did not ask for it.

What exists

• leakage per day, residence time and half-life tiers in aquifer_registry.yml; derive_leakage.py and estimate_aquifer_unknowns.py behind them
• core/model/aquifer.py treats the aquifer as a storage account with a dead fraction — a balance, not a groundwater simulation

What is missing — the accountancy proper:

• injection — how much can go in, at what rate, before head rise eats the headroom. Amplitude bounds it: an aquifer that naturally swings 2 m will not take a 20 m mound without pushing back. Nothing in the record contains an injection, so this cannot be validated against observation and must be presented as bounded, not predicted.
• baseflow return — leakage out of a leaky aquifer is not a loss if it surfaces in a river we are counting. Run 04 already credits 80% of leakage back as baseflow; that number is an assumption and should become a routed quantity once step 1 gives us river level.
• double-counting — the same water must not appear as aquifer storage AND as river baseflow in the same month. This is the specific failure mode this step exists to prevent.

Done looks like: a closed daily balance where injected volume, storage change, leakage and baseflow credit sum to zero within tolerance, with the recovery constrained by the compartment rule rather than by polygon area.

Deliberately parked

Not forgotten, and not being chased — the flows above are worth more than these:

Why parking these is defensible: they widen the range on volumes that are already labelled guess, and no decision currently turns on their central values. Why it stops being defensible: the moment a site is chosen and one of these bodies is near it. --site E,N on inventory_aquifers.py is what tells you which.

Also parked: the per-aquifer water-table fit currently applies to the Chalk only (n=48). Inferior Oolite (7), Great Oolite (6) and Carboniferous Limestone (6) all sit just under the n≥8 threshold — a handful more long records each would promote them.

← All notes · More from Aquifer Storage