Skip to content

Development history and validation record

This is the single record of how the current ecosys engine was built, how it was checked against the original Excel model, what was corrected along the way, and what is still different. The scientific decisions themselves live in scientific-decisions.md; the equations and the map from equations to code are in the LaTeX model reference. The step-by-step roadmaps, task logs and investigation notes this document was written from were removed from the working tree once their content was absorbed here. They remain in git history (last present at commit 477a718, and the legacy branch keeps that history).

1. Where this started

The first implementation, ecosys_legacy, read the original Excel parameter workbooks (ch1_orig.xls, ch2_orig.xls, ch3_orig.xls) directly and reproduced the saved TestCs137.xls results. Against that workbook it agreed to roughly 2.7e-6 on product tables, 2.7e-5 on cumulative ingestion, 6.2e-5 on ingestion categories and 7.9e-4 on annual dose rows. It is preserved on the legacy branch, with its tests, usage examples, algorithm audit and the original discrepancy study.

It was replaced because it was structured around the workbook, not the model: one mutable state object, global caches, routing by German workbook labels and VBA integer indexes, calculation mixed with reporting, duplicated logic, ambiguous time semantics, and unit-stripped floats. The engine was rewritten so that:

  • canonical SQLite data is the only source of parameters (production code never reads Excel, German labels or VBA indexes);
  • physical inputs carry astropy units at the API boundary;
  • kernels are pure, vectorized functions over batch axes, with validity masks separate from values;
  • the time support (the grid on which integrals are evaluated) is an explicit part of a request.

The Excel workbooks stay as ground-truth sources: ecosys.legacy_workbook loads them and export_excel_data writes the canonical database. Every source cell that was migrated is audited (dose coefficients, inhalation classes and skin-organ tables are checked cell by cell against the original sheets in tests_export_excel_data). See the reports in export_excel_data/.

2. Building the engine (22–24 September 2026)

The implementation ran as an ordered sequence of small, individually tested tasks, each gated on the full suite passing. Baseline before the work: 109 legacy tests, 12 exporter tests, 4 engine tests, no failures.

Phase What was built Gate evidence
Guardrails Import-boundary test (ecosys may not import legacy code or xlrd); scientific decision register; first decisions on physical air activity, year/calendar and cohort aging tests/test_architecture_boundaries.py, scientific-decisions.md
Typed data Immutable typed records for every canonical table; exact-ID ParameterCatalog; dense compiled arrays for plants, soil, processing graphs, animals, doses tests/data/
Domain inputs Unit-checked DepositionEvent, Landscape, PopulationCohorts, SimulationRequest; result axis metadata; batch-shape rules tests/domain/
Time and kernels Gregorian elapsed TimeAxis, seasonal schedules, decay, interpolation, exact exponential interval integrals tests/kernels/
Deposition, soil Vectorized dry/wet deposition, interception, LAI; soil migration, fixation, desorption, resuspension analytical tests
Plants and hay Grass, five plant categories, root uptake, translocation, harvest stocks, hay preparation analytical/property tests
Food chain Generic processing graph with storage; animal products with rations and retention; seasonal consumption and ingestion dose integration tests
Exposure Cloud/resuspension inhalation and ground external kernels analytical tests
Engine One public EcosysEngine.run(): multiple events, superposition, population dose tests/integration/, tests/validation/test_end_to_end_matrix.py
Scale Vectorization, chunk equivalence, masks, GIS-ready contract (arrays plus masks; no GIS dependency) tests/performance/, performance report

A first review found the engine returned deposition only, so a second pass (architecture remediation) closed these gaps before the release gate: an explicit absolute output origin for multi-event alignment; canonical consumption/activity profile resolution; complete typed result tree; compilation of soil-plant transfers, element-specific grass parameters, translocation curves and storage durations; one-time deposition at the event date propagated through soil, plants, feed, animals, foods, ingestion and exposure; per-event and collective population results; and end-to-end validation against independent hand-derived references (the earlier matrix only checked shapes). Suite at that point: 593 passed.

Time-support grid

The remediation exposed that integrated answers depend on the grid they are evaluated on (as they do in the VBA model). Instead of pretending otherwise (SD-18), the request now carries an explicit support:

  • SimulationRequest.on_default_grid(...) builds the deterministic VBA-style support (daily, two-day, five-day, monthly in year 3, then sparse long-term nodes; about 290 nodes), while a custom grid must start at the event origin;
  • animal intake uses left-held ration forcing on the support; delayed animal sources use a named piecewise-linear reconstruction, zero outside the horizon;
  • reporting (project_result) only selects from a completed run and never changes the integration support;
  • EcosysEngine.iter_chunks streams very large rasters (100,000 cells over 289 nodes in 40 chunks, about 3.8 GB traced peak) because the complete in-memory path is not intended for that size.

A later revision replaced the single-animal mast model by a continuously replenished herd (SD-10), which moved meat concentrations from far below to close to the workbook. The suite passed 641 tests at the release gate.

3. Validating against the Excel model

Two independent reference sets are used.

TestCs137.xls (legacy branch): a single scenario (Cs-137, 29 April 1986, 300 Bq/m³ for one day, 16,000 Bq/m² wet deposition, 4.5 mm rain, adult, arable), with tables of deposition, primary/feed/food concentrations, annual doses and cumulative ingestion by category. It was the reference during the first comparison campaign.

VM exports (tests/artifacts_ecosys_for_excel/run1–3.XLS): three independent runs of the Excel-for-ECOSYS-87 1.4D workbook (29 April, 2 May and 14 August 1986) with date-indexed chart series for grass, hay, milk, meat, root vegetables and oats. They are compared by scripts/compare_vm_exports.py and locked in tests/validation/test_vm_foodchain_exports.py.

The air-activity unit trap

All three exports display 300 Bq h/m³ but the workbook stores 25,920,000 Bq s/m³ (300 × 86,400): its unit-index cell is 0, and the VBA multiplies by 86,400 for every index other than the seconds and hours options. The dry ground deposition in the exports (12,960 Bq/m² at 0.0005 m/s) confirms it. Physically this is 300 Bq/m³ sustained for one day, equivalent to 7,200 Bq h/m³, 24 times the literal reading. The engine accepts a physical integral (SD-01) and no compensating multiplier is introduced; comparisons feed it the workbook's actual forcing. The close grass, root-vegetable and oat matches under that forcing corroborate the interpretation.

The same evidence disposed of a suspected "0.81 shortcut" in the old code: 0.81 = 12,960 / 16,000 is just dry over wet ground deposition for this case, not a VBA constant.

Evidence about the intended model

The Excel-model manual (chapter 3, Das Rechenmodell ECOSYS für Excel, 2000) documents a later model than the 1993 paper. Its prose (equations are embedded pictures and cannot be extracted) established:

  • plant types 1–5 and their between-harvest rules: hay is a mean over two preparation intervals; types 2 and 4 use the end-of-harvest concentration stored and decaying only; type 3 (leafy vegetables) is harvested year-round; type 5 uses the mean across the harvest period. The scenario workbook itself sets fruiting vegetables to type 4 and orchard fruit to type 5, opposite to the generic examples; the canonical data keeps the workbook's explicit type;
  • grass is harvested continuously in season and its direct contamination is neglected after a configured horizon (730 days);
  • hay preparation has a first-interval weight (multiplier 2 means 2/3 and 1/3; the paper uses 70/30) and a drying factor 5 for water loss;
  • fixation and desorption change plant-available soil activity; resuspended soil on plants carries a nuclide enrichment factor while grazed soil does not;
  • animals inhale as well as ingest, handled like feed intake in product transfer;
  • ingestion doses are potential doses assuming exclusively local food; the 10%-per-mm wet-deposition fallback is labelled arbitrary in the manual.

No numerical support rule or long-term animal equilibrium formula follows from the text.

4. Corrections and what they did

4.1 First alignment attempts

A trace of the 1,129-row comparison against TestCs137.xls located the early ingestion divergences: oats (harvest-day coordinate cliff; workbook dates are one-based, so first oat harvest is elapsed day 108, not 109), root vegetables (pre-growth foliar translocation must not occur for type-4/5 crops), early grazing products (ration forcing), early hay (a running primary series is not available feed), and food with zero reference but nonzero output. Fixing the first two moved raw oats at six months from 3.58 to 95.856 Bq/kg (workbook 95.861) and root-vegetable food from 34.4 to 1.595 Bq/kg (workbook 1.595); feed and food are now calculated in separate stages. One-year cumulative ingestion barely moved (2.075 → 2.061 mSv against 2.262), showing that downstream dose agreement is not a criterion for rejecting an upstream fix.

A parallel experiment aligned the model to the 1993 paper (accumulated grass growth dilution, additive leafy contamination, harvest-anchored carry-forward, no grass cutoff). It increased the disagreement (one-year ingestion 2.71 mSv) and was superseded once the manual showed the Excel model differs from the paper. The manual-informed configuration is the one in use (one-year 2.075 mSv at that point). Both runs were reproducible with the old benchmark script; the numbers are historical.

4.2 Food-chain corrections validated on the VM exports

Four corrections were implemented in order, each gated by focused tests and the three-export comparison (suite 665 → 679 tests). Values are clean Bq/kg against the VM value on the same exported support:

Correction Effect (example)
A. Recurring harvested stocks for types 2, 4, 5 (carry the harvest stock by physical decay until the next harvest) Run 3 primary oats day 715: 3.153 → 3.417 (VM 3.417)
B. Running hay preparation, then feed-domain storage Run 1 feed hay day 182: 10,866 → 13,839 (VM 13,889)
C. Linear seasonal ration interpolation by material identity Run 1 sheep milk day 7: 4.6 → 4,176 (VM 4,406)
D. Left-node intake during the transient, current-node equilibrium Run 2 cow milk day 7: 3,285 → 3,485 (VM 3,485)

Adult cumulative ingestion moved from 2.061 to 2.238 mSv at one year (VM 2.262) and from 2.975 to 3.124 mSv at 70 years (VM 3.108). Median absolute error over exported concentrations above 1 Bq/kg fell from 3.05% to 0.02% (feed), 8.9% to 0.20% (primary animal products) and 6.8% to 0.73% (food).

4.3 Divergence inventory (E1–E18)

Every divergence located in the three report-table comparisons is explained. Those that are model errors were fixed; the rest are conventions of the VBA/Excel implementation that were deliberately not reproduced.

# Divergence Cause Outcome
E1 Ground external dose −49% The engine used the ground deposition only; the VBA (and the paper's Eq. 6/22) use ground plus dry deposition on lawn Fixed (SD-12); ground dose within 0.04% from 5 to 70 y
E2, E3 Early ground dose low by t/(t+1); run 1 vegetable ingestion −0.024 mSv VBA clock: dose at row t covers the whole calendar day t ("deposition at 0:00, dose at 24:00"); every seasonal switch applies one day earlier Not reproduced: one physical clock (node t is t days after deposition)
E4 Early animal products offset, nonzero at day 0 VBA spreads animal inhalation uniformly over the deposition day; the engine used an instantaneous impulse Fixed (SD-17): first month within 0.0034%
E5 Grass drifting to +107% at 70 y Root-zone loss rates were taken from the request landscape instead of each plant's own soil type Fixed (SD-12): grass within 0.1%
E6, E7, E10, E11, E13, E14 Type-2/4/5 crop, leafy, beet-leaf, hay and growth-month offsets of 1–30% at specific times VBA anchors and averages over nodes of its sparse grid (last node in the harvest window, unweighted node means, freezes, right-node month rates) Not reproduced; each is predicted within 0.1% by reproducing the VBA rule
E8, E17 Arable crops +4…7% and hay +8.5% from 5 y After the third year the VBA drops stocks and uses current uptake or 5 × current grass; the engine keeps the last stock decaying physically Not reproduced (maintainer decision 2026-09-28: keep the more physical regime)
E9 Stored cereal food +15…18% late VBA selects the source node nearest to t − storage on its sparse late grid Not reproduced
E12 Brewers' grain +66% Feed graph chained a stored copy onto the stored copy Fixed
E15, E16 Winter grass marker; +0.0033% early VBA "not available" marker; VBA decays the whole state from deposition No effect / not reproduced
E18 Late sheep milk +13.8% / −8.2% VBA uses a constant first-year mean ration after year 3 Fixed (SD-10 revision: interval-mean equilibrium ration)

The ingestion quadrature was then changed deliberately away from the VBA's endpoint rule to a logarithmic interval mean with an onset step (SD-18). The default support is converged (70 y within 0.4% of daily/monthly support), and comparison with the VM shifts accordingly.

4.4 Skin and cloud

The source audit showed that cloud coefficients are Sv/s per Bq/m³ and the skin-organ table is Sv/s per Bq from 1 Bq on skin/clothing (the skin-tissue table is Sv/s per Bq/cm²). Cloud submersion is now an opt-in pathway using physical Bq s/m³ integrated air (reference scenario: 0.000209952 mSv, identical to the workbook). Skin-tissue dose is a separate result excluded from effective-dose totals; the workbook's Haut/Kleidung cells are effective-dose values and are not skin-tissue references (SD-01, SD-11).

5. Current state

Default support, adult, clean engine against the VM exports (cumulative ingestion):

Run 1 year 70 years Main remaining contributors
1 (29 Apr) −1.1% −0.4% E3 clock convention; late grass/crop conventions
2 (2 May) −0.7% −0.6% meat −2%; late conventions
3 (14 Aug) +1.0% +1.0% cereals +8.8% (E9); vegetables +1.8%

Against TestCs137.xls the earlier one-year comparison (1.985 vs 2.262 mSv ingestion) predates the herd model and the corrections above and is only historical; the VM comparison is the current one. Supported pathways include ingestion, cloud and resuspension inhalation, ground and (on request) cloud external, plus skin tissue separately. The maintained suites have 688 tests (2026-09-29, after the fixes for the code audit).

Remaining differences and open questions

  • Support dependence of the default long-term ingestion is much reduced but not eliminated; quadrature between nodes remains an open part of SD-18.
  • VM animal products are zero at t = 0 while the engine spreads the deposition-day inhalation (SD-17); some delayed feed/food and stored-source node selections still differ from the VM; feed potatoes, rye and acid whey are not requested by the feed comparator.
  • Decision register items still open: SD-13 (combined nuclides and decay chains: no ingrowth is computed), SD-14 (definition of "annual" dose; the engine returns interval and cumulative dose), SD-15 (general pathway applicability policy).
  • The VBA seasonal consumption factors apply only while t ≤ 1096 days; the engine applies them for the whole horizon (not yet measured).
  • Comparison scripts that depended on the legacy time-grid helper (compare_vm_report_tables.py, benchmark_full_validation.py) were removed with the legacy code, so the report-table comparison behind E1–E18 cannot currently be re-run from this branch; rank_vm_report_table_divergences.py needs its CSV output. Note when porting: the legacy helper hard-coded "1 month" as day 30; the correct anniversary is day 31 for 2 May and 14 August. Chart-series comparison (compare_vm_exports.py) and support convergence (scripts/foodchain_support_convergence.py) still work.

6. Reproducing the evidence

python -m pytest -q                                    # full suite
python scripts/compare_vm_exports.py                   # chart series vs run1-3.XLS
python scripts/foodchain_support_convergence.py        # default vs refined support
python benchmarks/benchmark_engine.py                  # full engine timing/memory
python benchmarks/benchmark_kernels.py                 # kernel timing
python -m export_excel_data --help                     # workbook -> SQLite export

Timing is evidence, not a test. The parameter data itself is documented in data-model/.

7. Packaging decisions

The engine is the public ecosys package (the earlier ecosys_new name was dropped); it ships its canonical SQLite data. The legacy implementation was distributed temporarily alongside it and now lives only on the legacy branch. The workbook importer is kept as ecosys.legacy_workbook with the three source workbooks, the only part of ecosys allowed to depend on xlrd. No GIS framework or other new runtime dependency was introduced.