# Boundary and initial-state audit

The domain is connected, but its edges and deep layer contain strong physical assumptions. The continuing accumulation cannot be diagnosed from solver closure alone. This audit does not change those assumptions or run an alternative model.

All quantities below were measured by `audit.py` in this folder. Reproduce using the command in `receipt.json`; it binds the native graph, model text files and their external static binary arrays. The source is the original `combined_feedback2_dec206_precise` G2 generation, not a newer daily-coupling generation.

## Outer boundary and connections

- **One connected native graph:** 33,378 active nodes, comprising 16,689 cells in each layer, with a vertical connection at every active column. No separate disconnected aquifer block explains the accumulation. This is graph connectivity, not proof that all connections carry significant water when cells are dry.
- **No CHD, GHB or HFB package** is present. There are 1,500 missing lateral neighbour faces in each layer, touching 1,025 upper cells. These include active/inactive interfaces, not only the exterior outline. Groundwater cannot move laterally into the excluded cells through those faces. There is no deeper modeled layer below layer 2, nor a separate basal-flow boundary. Interpretation: side and basal faces without explicit exchange impose no flow, rather than a measured regional head or underflow relation. IDOMAIN zero excludes a cell from the solution ([MODFLOW DIS documentation](https://modflow6.readthedocs.io/en/latest/_mf6io/gwf-dis.html)).
- The water system is **not closed**: recharge, evaporation, wells, rivers and drains add/remove water. The distinction is between exchange through these packages and regional groundwater exchange across the truncated sides/base.
- The northern drain is **17 interior upper cells**, all away from the active/inactive edge. It is a finite one-way outlet at terrain minus 2.5 m, with conductance 500 m²/day, not a fixed-head open northern boundary. Source lineage: `docs/records/model-domain-review-2026-09-23/lowland_mf6_candidate/run.py:289` constructs an interior transect after removing its river-terminal endpoint cells.
- The prescribed mountain source has 134 upper cells, 65 at the edge; the pumping package has 3,279 lower cells. Antecedent sums are +322,675.01 and −643,051.58 m³/day respectively. These are prescribed source/sink terms, not head-dependent lateral boundaries. Read: the source `build_receipt.json` identifies the mountain term as a distributed 35% component of a mixed native deep-aquifer pool. Its paired bookkeeping does not independently establish the recipient flux's physical magnitude or placement.

`boundary_packages.csv` gives all current drain/river cell counts, edge intersections, elevations and conductances. The surface/drain catalogue covers the active upper footprint, but a drain only discharges when its activation elevation is exceeded. Complete coverage is therefore not evidence of effective outlets at a deep water table.

## Aquifer and storage representation

Measured from NPF/STO: the upper layer has ICELLTYPE=1 and ICONVERT=1 everywhere; the lower layer has both flags zero everywhere. The upper layer can lose saturated thickness as head falls. The lower layer always retains full hydraulic thickness and confined storage in the equations. Setting specific yield to 0.05 in the input for both layers does not make it an active lower-layer drainable-storage parameter. [MODFLOW NPF](https://modflow6.readthedocs.io/en/latest/_mf6io/gwf-npf.html), [STO](https://modflow6.readthedocs.io/en/latest/_mf6io/gwf-sto.html).

Measured: upper thickness ranges 50–250 m; lower thickness 80–250 m. Upper specific storage is 1e-6 per metre and specific yield 0.05; lower specific storage ranges 1e-6–6.25e-6 per metre. Lower thickness times specific storage ranges 0.00015–0.0005. Conductivities and all discrete property values are in `receipt.json`. Their presence and numerical positivity do not qualify their spatial assignment against geology.

Read: `lowland_mf6_candidate/run.py:229–242` assigns default properties to the expanded footprint and replaces them with zonal properties only within `data["old"]`; lines 270–272 apply the uniform layer-type flags. Thus whole-domain consistency of seepage does not remove inherited property-footprint assumptions, nor prove that a wholly confined lower layer is appropriate everywhere.

## What the initial state really is

Measured: the IC arrays equal terrain minus 10 m in both layers to roundoff. That is a numerical starting guess. It is **not** the state used directly for the dated historical run. The first period is steady state and subsequent periods transient. Read: `production_adoption/mf6/combined_candidate.py:256` and source `build_receipt.json.antecedent` define that first solve from full-calendar mean recharge, residual evaporation, pumping and prescribed mountain input, plus the mean capped river stage. The receipt explicitly says this is not proof of equilibrium. A steady solve under averaged stresses is also not a periodically equilibrated seasonal state.

Measured from the saved first solution: 3,043 upper heads are at/below the upper bottom; 670 lower heads are at/below the lower bottom. Before 2017 those counts remain 2,871 and 659. The assumed lower-layer full saturation remains in the equations even where head is below its bottom. This warrants a physical consistency review of aquifer extent, bottoms and layer type; it is not evidence of a physical negative-depth water table. `initial_state.csv` records these counts and the above-terrain counts separately.

## Discriminating interpretation with existing data

1. In the accumulating cells, separate native storage gain from changing below-bottom heads. A large numerical head movement in a dry upper layer is not proportional to stored water.
2. Compare annual recharge/mountain/river supply, lateral inflow and actual outlet discharge by accumulating cluster. If accumulation is mainly in cells well below outlet elevations, changing a surface drain's resistance may have little immediate effect there. The parent investigation's retained cell budgets can test this without a new physical assumption.
3. Compare original-control and R1000 repeated-year trajectories from their common starting state. Shared drift indicates a baseline forcing/boundary/storage response; the difference isolates the seepage change more closely. Neither result establishes that regional exchange is physically missing.
4. Use `outer_edge_cells.csv` to check whether accumulating clusters intersect truncated sides. A sidewall exchange sensitivity is defensible only after identifying the relevant external aquifer and independent head/underflow evidence. Do not add a general-head boundary merely to dispose of a surplus.
5. Check the lower confined assumption and layer-bottom geometry specifically where computed lower head is below its bottom. Check property changes against the old footprint separately from country. These are more discriminating than changing specific yield to accelerate settling.

No alternative physics was implemented. A closed model budget demonstrates correct accounting for the boundaries supplied; it cannot demonstrate that those boundaries represent the actual aquifer.

## 2017 is a different forcing regime from the initialization

Measured by `audit.py --forcing-comparison`, using actual native CBC rates at total solver time 1 (step 1, period 1), compared with the exact retained original-model 2017 control budget. The original dated times are 5116–5480, corresponding to 1 January–31 December 2017; the replay uses local times 1–365. Rates below are daily means, not cumulative volumes.

| Term | Steady antecedent, million m³/day | Original 2017 mean, million m³/day | Ratio |
| --- | ---: | ---: | ---: |
| Prescribed recharge RCHA | 4.900092 | 12.457172 | 2.54223 |
| Prescribed mountain input | 0.322675 | 0.460086 | 1.42585 |

The combined prescribed recharge-plus-mountain input increases by 7.694491 million m³/day. The original 2017 mean storage gain is 4.530710 million m³/day; drain discharge, groundwater evaporation, pumping and net river export also increase relative to the antecedent. These response terms are in `antecedent_vs_2017_budget.csv`. They must not be described as independently prescribed changes simply because their budget values differ.

**Interpretation:** repeating 2017 applies substantially greater prescribed recharge than the full-calendar mean used for initialization. Storage gain under those repeated stresses is therefore not, by itself, evidence of a missing outlet or incorrect storage coefficient. The candidate may need to adjust to a different forcing regime even with unchanged physics. The comparison establishes the inputs and resulting accounts, not whether 2017 recharge is physically correct, whether it is representative for conditioning, or whether the present boundary conditions are valid. Continue comparing repeated original-control and candidate trajectories; keep historical forcing differences distinct from seepage-treatment effects.

The first steady solution is a numerical solution under averaged stresses, not an observed state and not a seasonal equilibrium for 2017. `antecedent_comparison_receipt.json` binds extracted native arrays, original replay budget, exact time identity and units. Only the antecedent native records were read from the large CBC; no full-history budget scan or new simulation was needed.

### Check against a stale antecedent recharge field

Measured by `audit.py --recharge-calendar`: **the stale-antecedent-recharge hypothesis is ruled out for this source generation.** All 6,210 actual G2 recharge binary files were reread: the steady field plus 6,209 dated daily fields covering 2003–2019. The actual steady domain rate and the recomputed, day-weighted transient-calendar mean both equal 4,900,092.440508737 m³/day. The maximum difference at an individual cell is 1.47e-10 m³/day, numerical roundoff. Steady and 2017 input totals also reproduce the respective native budget totals exactly. This does not rely on the builder's label.

The 2017 actual mean is 12,457,172.33182567 m³/day, or 2.542232 times the actual transient-calendar mean. It is the largest annual prescribed recharge in this generation; 2018 is second. These are modeled recharge rankings, **not a finding that those years were meteorologically wet or that the recharge inputs are physically correct**. Repeated-2017 conditioning is consequently a deliberate repetition of the highest-recharge year in the available modeled calendar. Interpret that stress test separately from equilibrium under the historical sequence.

`native_recharge_calendar.csv` records each date, actual input total, filename and checksum. `native_recharge_annual.csv` records annual totals, exact day counts including leap years, and ratios. `native_recharge_calendar_receipt.json` binds the comparisons and maximum cell discrepancy. No source files or model parameters were changed.