diff --git a/README.md b/README.md
index 0e9db74..7639ad1 100644
--- a/README.md
+++ b/README.md
@@ -50,6 +50,12 @@ spatial split, censoring/detection settings, season exclusions):
 uv run snakemake --cores 4 --rerun-triggers=mtime -- model_all
 ```
 
+Each rule depends on a stamp of the config section its run_id hashes
+(`outputs/config_stamps/`, written by the Snakefile at parse time), so editing
+the grid section rebuilds the grid and everything after it, editing a split or
+train setting reruns only those stages, and a comment-only edit or a
+re-pinned snapshot triggers exactly the stages it affects.
+
 This builds the full-ice-sheet covariate grid, the spatially-blocked
 train/test split, trains the configured models with cross-validation, and
 writes per-model outputs to `outputs/model/<model>/`:
diff --git a/Snakefile b/Snakefile
index a086ccb..6a5b761 100644
--- a/Snakefile
+++ b/Snakefile
@@ -13,20 +13,56 @@
 # the filenames, so every rule can target a static path directly.
 
 STORE = config.get("store", "ase")
-CONFIG_PATH = config.get("config_path", f"config/{STORE}.yaml")
 OUT_DIR = config.get("out_dir", "outputs")
 
 # Cross-store modeling pipeline (grid -> split -> train -> benchmark), driven by
 # a single config that names its input stores and models. Invoked without store=:
 #   uv run snakemake --cores 4 model_all
 MODEL_CONFIG_PATH = config.get("model_config_path", "config/model.yaml")
-import yaml as _yaml
-with open(MODEL_CONFIG_PATH) as _f:
-    _model_cfg = _yaml.safe_load(_f) or {}
-MODEL_STORES = sorted({s for stores in _model_cfg.get(
-    "inputs", {"antarctic": ["ase", "utig"], "greenland": ["greenland"]}).values() for s in stores})
-MODELS = [m["name"] if isinstance(m, dict) else m
-          for m in _model_cfg.get("train", {}).get("models", [{"name": "linear"}])]
+import os as _os
+from pathlib import Path as _Path
+from radar_postproc.config import config_hash, load_config, load_model_config
+
+_model_cfg = load_model_config(MODEL_CONFIG_PATH)
+MODEL_STORES = sorted({s for stores in _model_cfg["inputs"].values() for s in stores})
+MODELS = [m["name"] for m in _model_cfg["train"]["models"]]
+
+
+def config_stamp(name: str, section: dict) -> str:
+    """Path of a stamp file that changes exactly when `section` changes.
+
+    Each stage hashes only its own config section into its run_id, so each
+    rule depends on a stamp of that same section rather than on the whole
+    config file: editing train settings does not rebuild the grid, and a
+    comment-only edit rebuilds nothing. Written at parse time, only when the
+    hash differs (so mtime moves only on a real change); a stamp created for
+    the first time gets an epoch mtime, so migrating existing outputs does
+    not trigger a rebuild.
+    """
+    path = _Path(OUT_DIR) / "config_stamps" / f"{name}.hash"
+    digest = config_hash(section)
+    if not path.exists():
+        path.parent.mkdir(parents=True, exist_ok=True)
+        path.write_text(digest)
+        _os.utime(path, (0, 0))
+    elif path.read_text() != digest:
+        path.write_text(digest)
+    return str(path)
+
+
+def augment_stamp(store: str) -> str:
+    """Stamp of a store's full augment config (what its run_id hashes)."""
+    path = config.get("config_path", f"config/{store}.yaml")
+    return config_stamp(f"augment_{store}", load_config(path))
+
+
+# Section hashes below mirror grid.py / split.py / train.py exactly.
+_tcfg = _model_cfg["train"]
+STAMP_GRID = config_stamp("grid", {"inputs": _model_cfg["inputs"], "grid": _model_cfg["grid"]})
+STAMP_SPLIT = config_stamp("split", {"inputs": _model_cfg["inputs"], "split": _model_cfg["split"]})
+STAMP_TRAIN = {m["name"]: config_stamp(f"train_{m['name']}",
+                                       {"train": {**_tcfg, "models": [m]}, "model": m["name"]})
+               for m in _tcfg["models"]}
 
 # Which trained model the mission design tool ships. Lives in its own config
 # section so it never enters a stage's section hash (and so never perturbs a
@@ -44,6 +80,11 @@ rule all:
 
 
 rule run:
+    # Depends on a stamp of the store's config, so re-pinning
+    # icechunk.snapshot_id (or any other augment change) re-runs the
+    # extraction without --forcerun.
+    input:
+        stamp=lambda w: augment_stamp(w.store),
     output:
         parquet=f"{OUT_DIR}/{{store}}/{{store}}.parquet",
         manifest=f"{OUT_DIR}/{{store}}/{{store}}.manifest.json",
@@ -52,10 +93,8 @@ rule run:
 
         # Derive the config from the wildcard (not the module-level CONFIG_PATH)
         # so cross-store DAGs (e.g. model_all) build the right store.
-        run_pipeline(
-            config.get("config_path", f"config/{wildcards.store}.yaml"),
-            out_dir=OUT_DIR,
-        )
+        run_pipeline(config.get("config_path", f"config/{wildcards.store}.yaml"),
+                     out_dir=OUT_DIR)
 
 
 rule csv:
@@ -94,7 +133,9 @@ rule model_all:
 
 rule grid:
     # Full-ice-sheet covariate grid (the prediction domain). Network + slow;
-    # only re-runs when config/model.yaml's grid section changes.
+    # only re-runs when config/model.yaml's inputs/grid sections change.
+    input:
+        STAMP_GRID,
     output:
         parquet=f"{OUT_DIR}/model/grid.parquet",
         manifest=f"{OUT_DIR}/model/grid.manifest.json",
@@ -108,6 +149,7 @@ rule split:
     # Attach the radar target to grid points, assign blocking cells + folds,
     # emit helper maps for hand-picking test cells.
     input:
+        STAMP_SPLIT,
         grid=f"{OUT_DIR}/model/grid.parquet",
         stores=expand(f"{OUT_DIR}/{{s}}/{{s}}.parquet", s=MODEL_STORES),
     output:
@@ -122,6 +164,7 @@ rule split:
 
 rule train:
     input:
+        lambda w: STAMP_TRAIN[w.model],
         f"{OUT_DIR}/model/split.parquet",
     output:
         metrics=f"{OUT_DIR}/model/{{model}}/metrics.json",
diff --git a/agent_notes/20260901-calibration-qc-results.md b/agent_notes/20260901-calibration-qc-results.md
new file mode 100644
index 0000000..848e1f0
--- /dev/null
+++ b/agent_notes/20260901-calibration-qc-results.md
@@ -0,0 +1,124 @@
+# Calibration QC filter — results (2026-09-01)
+
+Branch `qc-calibration-filter`. Plan: `agent_plans/20260901-calibration-qc-filter.md`.
+All numbers reproducible with `uv run python scripts/qc_filter_comparison.py`
+(needs `outputs/baseline_20260807/`, a copy of the pre-QC `outputs/model` +
+augment parquets) and `scripts/season_crossover_matrix.py [--calibration-qc]
+--out-dir outputs/qc_filter/crossover`. Full tables: `outputs/qc_filter/comparison.md`.
+
+## Data
+
+Pinned to the 2026-09-03 15:51 UTC method-0.4.1 snapshots (antarctica
+087JBD7NTAE8BTBTEYSG, greenland ACA8WY61ZF9W6VSBA3HG, ase F6TXWSEQ9RQCD1H7MSMG).
+Upstream verified no science values changed, so baseline vs new differs only
+by the filter. The 0.4.1 refresh (from the 2026-09-02 0.4.0 snapshots first
+used here) changed exactly one thing that reaches the model: 72 traces of
+2018_Antarctica_DC8 lost their at-ceiling flag when its rising img2 envelope
+stopped counting as a ceiling — all 72 were already rejected by the img2
+rule, so the split parquet is identical apart from run_ids and the model
+retrains to the same posterior (fixed seed).
+
+Share of traces entering the split stage (after the two season exclusions)
+rejected, exclusive attribution:
+
+| rule | Antarctica | Greenland |
+|---|---|---|
+| seam \|Δ\| ≥ 3 dB | 9.1% | 8.9% |
+| surface from img2+ | 3.1% | 3.0% |
+| margin < 2 dB | 2.2% | 0.7% |
+| any | 14.4% | 12.7% |
+
+Worst seasons: 2016_Greenland_P3 51% (49% seam), 2017_Antarctica_Basler 46%
+(already excluded), 2013_Antarctica_Basler 33%, 2017_Antarctica_P3 31%,
+2019_Antarctica_GV 24% (22% at its −30 dB ceiling), 2016_Antarctica_DC8 22%.
+Unmeasured seams (check could not run): 92% of 2014_Greenland_P3, 61% of
+2012_Antarctica_DC8 → `drop_unmeasured_seam: true` would reject 35% / 50% of
+all traces, i.e. whole seasons; kept false. Threshold sensitivity: seam 2 dB →
+27% / 18% rejected, 5 dB → 12% / 7%; margin 1–3 dB moves totals by ±2 pts;
+dropping the img2 rule saves 2.6 / 2.5 pts. img2-sourced traces have median
+RSSNR 8–17 dB below img1 traces of the same DC8 season (2012: 46.7 vs 38.4,
+2014: 59.5 vs 42.1, 2016: 46.7 vs 35.9, 2018: 65.3 vs 50.2 dB) — the bias the
+changelog warns of is visible in the target, so the rule stays.
+
+Training grid points 19,272 → 17,447 (−9.5%); 714 kept points now match a
+different, QC-passing trace. Median training target rises 61.9 → 64.6 dB
+(Antarctica) and 73.6 → 74.0 dB (Greenland): the filter preferentially
+removes low-RSSNR traces (saturated surface ⇒ RSSNR biased low).
+
+## Model fit
+
+| | atten_refl baseline | atten_refl QC | linear baseline | linear QC |
+|---|---|---|---|---|
+| CV RMSE [dB] (fold range) | 13.02 (11.9–14.4) | 12.87 (11.8–14.3) | 13.92 | 13.80 |
+| CV coverage 1σ | 0.680 | 0.679 | 0.704 | 0.702 |
+| CV log score [dB] | −4.019 | −4.009 | −4.103 | −4.096 |
+| test RMSE (own test set) | 12.67 | 12.72 | 13.95 | 13.98 |
+| σ residual [dB] | 12.91 | 12.71 | 14.33 | 14.15 |
+| θ / τ [dB] | −1.22 / 3.01 | −1.08 / 2.79 | | |
+| divergences / R̂max | 0 / 1.0025 | 0 / 1.0049 | | |
+
+Same points, both posteriors (the fair comparison — the CV populations differ):
+on the QC test set (1,330 uncensored) atten_refl baseline 12.75 dB vs QC 12.72
+dB; on the baseline test set (1,433) 12.67 vs 12.64 dB. Bias +1.9 → +2.0 dB.
+Differences of 0.03–0.05 dB are noise.
+
+Posteriors (`posterior_comparison.png`): every parameter within ~2 sd of its
+baseline. Largest moves: β_a[greenland] −2.1 → −3.6 dB/km, −β_r[greenland]
+−16.0 → −17.2 dB, −α_r +6.05 → +5.66 dB, −β_r[T_air] +0.25 → +0.28 dB/K,
+σ 12.9 → 12.7 dB, τ 3.0 → 2.8 dB. The sheet-level attenuation and
+reflectivity offsets move against each other, so predictions barely change.
+
+Predictions (`prediction_difference.png`): new − baseline posterior mean over
+the full grid, median +0.01 dB (Antarctica) / +0.07 dB (Greenland), 5th–95th
+percentile −0.2…+0.6 / −0.3…+0.6 dB; predictive std ratio 0.985. Largest
+positive shifts on the Siple Coast / Ross ice streams and NW Greenland
+(≈ +0.6–1 dB); slightly negative in interior SW Greenland.
+
+## Model-free check: season crossovers
+
+`outputs/qc_filter/crossover/crossover_summary.md` (off-diagonal cells present
+both before and after):
+
+| | mean \|median Δ\| pre → QC | mean sd pre → QC | diagonal sd pre → QC |
+|---|---|---|---|
+| Antarctica (18 cells) | 11.65 → 11.90 dB | 9.84 → 9.69 dB | 8.61 → 7.71 dB |
+| Greenland (15 cells) | 7.13 → 6.53 dB | 9.03 → 7.55 dB | 7.06 → 6.92 dB |
+
+- **2016_Greenland_P3** (49% seam rejects) is the clear success: vs
+  2014/2017/2018/2019 it goes from −4.3/+2.1/+2.5/−1.2 dB with sd 11–13 dB to
+  −0.2/−1.2/−0.6/−3.3 dB with sd 7.7–8.8 dB — i.e. it now looks like every
+  other P3 season, and its pairs with 2013 drop from −18.8 to −14.9 dB.
+- Within-season repeatability tightens where seams were flagged:
+  2013_Antarctica_Basler 10.4 → 8.4 dB, 2017_Antarctica_P3 8.2 → 5.6 dB,
+  2014_Antarctica_DC8 9.4 → 8.4 dB, 2017_Antarctica_Basler 17.9 → 13.7 dB.
+- **Season-level offsets do not move**: 2012_Antarctica_DC8 stays ~12 dB below
+  the 2014–2018 DC8/P3 seasons; 2017_Antarctica_Basler stays ~30 dB below its
+  neighbours (and 9–13 dB below the 2022/23 BaslerMKB); 2013_Antarctica_P3 vs
+  2019_GV stays +15 dB; 2013_Greenland_P3 stays 16–19 dB low. The QC catches
+  2% of 2013_Greenland_P3 and 46% of 2017_Antarctica_Basler, so
+  `exclude_collections` must stay as it is.
+
+Out-of-fold residuals by season (`scripts/residual_audit.py`, now in
+`docs/figures/residuals_by_season.png`): per-season medians move by <= 1.4 dB
+(2012 DC8 +7.7 -> +7.9, 2013 P3 -5.3 -> -6.7, 2019 GV +11.2 -> +12.2,
+2017 P3 -9.2 -> -9.6); the same seasons stay offset.
+
+## Assessment
+
+The suggested filter is sound hygiene and cheap (≈13% of traces, no
+retuning), and it demonstrably repairs the one season whose problem *is*
+trace-level (2016 Greenland seams). But the required-SNR model was already
+insensitive to it: at the 5 km grid / 1 km nearest-trace matching scale the
+rejected traces were a minority whose biases partly average out, and the
+dominant residual scatter (σ ≈ 12.7 dB) and the season-level calibration
+offsets are untouched. Adopt it as the default, keep the two exclusions, and
+treat season-level (crossover-derived) calibration as the next lever — not
+QC threshold tuning, which the sensitivity table shows only trades data
+volume for no measurable fit change.
+
+Open follow-ups: 2012_Antarctica_DC8's img1 ceiling is now fit_ok at −29.7 dB
+(4% of traces at the ceiling) but its ~12 dB offset to later DC8 seasons is
+a whole-season effect; 2019_Antarctica_GV loses 22% to its −30 dB ceiling
+(pileup 0.11, span 0.35 — a credible fit, worth a look upstream);
+2023_Antarctica_BaslerMKB's ceiling (pileup 0.01) is a weak fit that removes
+0.7% — harmless but not evidence of saturation.
diff --git a/agent_plans/20260901-calibration-qc-filter.md b/agent_plans/20260901-calibration-qc-filter.md
new file mode 100644
index 0000000..1eb673f
--- /dev/null
+++ b/agent_plans/20260901-calibration-qc-filter.md
@@ -0,0 +1,50 @@
+# Radiometric calibration QC filter (2026-09-01)
+
+Branch: `qc-calibration-filter`.
+
+## Context
+
+Upstream `radar-return-statistics` shipped calibration method 0.4.0 (refreshed to 0.4.1 on 2026-09-03: 2018_Antarctica_DC8's rising img2 envelope no longer counts as a ceiling, 49 Antarctic frames backfilled, variable attrs added)
+(`docs/dataset_changelog.md`, 2026-09): four per-trace diagnostics —
+`img_comb_offset_dB`, `img_comb_pair`, `surface_source_image_index`,
+`surface_ceiling_margin_dB` — plus a season-keyed `saturation` root attr. No
+pre-existing values changed (verified byte-identical upstream), so a model
+retrained on the re-pinned snapshots differs from the 2026-08-07 baseline only
+through the QC filter.
+
+Stores landed 2026-09-02 00:36 UTC (0.4.0: antarctica SYAKG11X8AFFY0H9H65G,
+greenland JWDABR34HPM816P70FD0, ase WQCXS05H226PX1KGRZ2G) and were refreshed
+2026-09-03 15:51 UTC (0.4.1, the pinned ones: antarctica 087JBD7NTAE8BTBTEYSG,
+greenland ACA8WY61ZF9W6VSBA3HG, ase F6TXWSEQ9RQCD1H7MSMG). utig / crosssystem not updated (not model inputs).
+
+## Plan
+
+1. Carry the four calibration columns through augment (`extract.carry_columns`
+   default; warn-skipped on stores without them).
+2. Implement the suggested filter as `split.calibration_qc` (split stage, next to
+   `exclude_collections`) so thresholds can be varied without re-running augment,
+   and so it applies to observations AND non-detections before matching:
+   - seam: `|img_comb_offset_dB| >= 3 dB` -> drop; NaN passes unless
+     `drop_unmeasured_seam` (dropping NaN would remove 93% of 2014_Greenland_P3
+     — `insufficient_overlap` is a geometry limitation, not evidence of a step).
+   - img2: `surface_source_image_index >= 2` -> drop (-1 unknown passes). Kept
+     as suggested: img2-sourced traces have median RSSNR 8–17 dB lower than
+     img1 traces of the same DC8 season (surface-power low bias), and the rule
+     costs only 3% of Antarctic / 4% of Greenland traces.
+   - saturated: `surface_ceiling_margin_dB < 2 dB` -> drop; NaN (no credible
+     season ceiling) passes.
+3. Back up the pre-QC outputs (`outputs/baseline_20260807/`: model/ + augment
+   parquets), re-pin snapshots, `snakemake model_all`.
+4. Compare with `scripts/qc_filter_comparison.py` (per-season removal, threshold
+   sensitivity, benchmark side by side, both posteriors on the SAME test points,
+   posterior overlays, prediction-difference maps) and
+   `scripts/season_crossover_matrix.py --calibration-qc` (model-free check).
+5. Keep `exclude_collections` unchanged for the main comparison (one change at a
+   time); report how much of each excluded season the QC alone would catch.
+
+## Status
+
+- [x] 1–2 implemented + unit tests (`tests/unit/test_calibration_qc.py`).
+- [x] 3 pipeline run (2026-09-02 on 0.4.0; re-run 2026-09-03 16:13–16:27 UTC on 0.4.1 — identical posteriors; run_ids split 4ff998437f3b, atten_refl af5514ae3b72, linear c0b385764662)
+- [x] 4 comparison + assessment (`agent_notes/20260901-calibration-qc-results.md`)
+- [x] docs (`docs/2_input_data.md`); Snakefile: store + model configs are now rule inputs
diff --git a/config/antarctica.yaml b/config/antarctica.yaml
index 6e559f7..309474f 100644
--- a/config/antarctica.yaml
+++ b/config/antarctica.yaml
@@ -8,7 +8,7 @@ store:
 
 icechunk:
   branch: "main"
-  snapshot_id: "VE13DY3F546ZE1J9KC60"  # pinned 2026-08-03: reprocessed with missing bed picks + pick-free noise stats
+  snapshot_id: "087JBD7NTAE8BTBTEYSG"  # pinned 2026-09-03: calibration method 0.4.1 refresh, no science values changed
 
 extract:
   qc_only: true
diff --git a/config/ase.yaml b/config/ase.yaml
index 023b6a8..5f7a8a2 100644
--- a/config/ase.yaml
+++ b/config/ase.yaml
@@ -7,7 +7,7 @@ store:
 
 icechunk:
   branch: "main"
-  snapshot_id: "CS9HEDFGE8FXDHVC9KM0"  # pinned 2026-07-30: tail noise window shifted off the end-of-record taper
+  snapshot_id: "F6TXWSEQ9RQCD1H7MSMG"  # pinned 2026-09-03: calibration method 0.4.1 refresh, no science values changed
 
 extract:
   qc_only: true
diff --git a/config/greenland.yaml b/config/greenland.yaml
index 2189039..26b3b19 100644
--- a/config/greenland.yaml
+++ b/config/greenland.yaml
@@ -7,7 +7,7 @@ store:
 
 icechunk:
   branch: "main"
-  snapshot_id: "GEAMAHQ7BRVPG9SQPK20"  # pinned 2026-07-31: reprocessed with missing bed picks + new noise metrics
+  snapshot_id: "ACA8WY61ZF9W6VSBA3HG"  # pinned 2026-09-03: calibration method 0.4.1 refresh, no science values changed
 
 extract:
   qc_only: true
diff --git a/config/model.yaml b/config/model.yaml
index 1ae5cf2..ed0d6e0 100644
--- a/config/model.yaml
+++ b/config/model.yaml
@@ -41,6 +41,15 @@ split:
   seed: 42
   test_cells: ["ant:-5:1", "ant:-3:-2", "ant:1:-3", "ant:2:-2", "ant:2:0", "grl:0:-6", "grl:1:-4", "grl:0:-2"] # pick from outputs/model/cell_maps/, e.g. ["ant:-3:1", "grl:0:-5"]
   exclude_collections: [2017_Antarctica_Basler, 2013_Greenland_P3]  # crossover-identified calibration outliers
+  # Per-trace radiometric QC from the upstream calibration fields (dataset
+  # changelog 2026-09, method 0.4.1). Rules pass where their input is
+  # unmeasured (NaN / unknown source image) — see docs/2_input_data.md.
+  calibration_qc:
+    enabled: true
+    max_seam_offset_dB: 3.0      # image-combine seam step |img_comb_offset_dB| >= this -> drop
+    drop_unmeasured_seam: false  # NaN offsets (check not run) are kept
+    require_img1_surface: true   # surface sampled from a higher-gain image (index >= 2) -> drop
+    min_ceiling_margin_dB: 2.0   # surface within this of its season's clip level -> drop
 
 train:
   seed: 42
diff --git a/docs/2_input_data.md b/docs/2_input_data.md
index 14385f2..ec98ecf 100644
--- a/docs/2_input_data.md
+++ b/docs/2_input_data.md
@@ -24,6 +24,28 @@ Traces with radar-derived thickness under 100 m are dropped. (This is configured
 
 **Non-detections** are deliberately included. These are defined as locations where there is no bed pick available but there are bed picks both before and after in the same radar segment. This is used as a heuristic to filter out anything where bed picking simply hasn't been done.
 
+### Radiometric calibration QC
+
+Two instrument effects can bias RSSNR at the trace level, and since 2026-09 the upstream stores ship per-trace diagnostics for both (calibration method 0.4.1; see the [dataset changelog](https://github.com/englacial/radar-return-statistics/blob/main/docs/dataset_changelog.md)):
+
+* **Image-combine seam steps.** The MCoRDS products stitch a low-gain image (surface) onto higher-gain images (deep ice/bed). A miscalibrated stitch puts a power step at the seam, which biases bed power relative to surface power by that step. `img_comb_offset_dB` is the residual step measured in each frame's individual images.
+* **Surface saturation.** Where the surface return clips the receiver, surface power is underestimated and RSSNR is biased low. `surface_ceiling_margin_dB` is the distance below the season's fitted clip level, and `surface_source_image_index` records whether the surface sample came from the low-gain image (1) or a higher-gain image (≥2, likely saturated with a season-dependent low bias).
+
+The split stage applies the upstream-suggested filter (`split.calibration_qc` in `config/model.yaml`) to observations and non-detections alike, before grid matching:
+
+| rule | drops traces where | share dropped, Antarctica / Greenland |
+|---|---|---|
+| seam step | \|`img_comb_offset_dB`\| ≥ 3 dB | 9.1% / 8.9% |
+| higher-gain surface | `surface_source_image_index` ≥ 2 | 3.1% / 3.0% |
+| at the ceiling | `surface_ceiling_margin_dB` < 2 dB | 2.2% / 0.7% |
+| any | | 14.4% / 12.7% |
+
+(Shares are of the traces entering the split stage, i.e. after the season exclusions above; a trace failing several rules is counted once, in the first row it fails.)
+
+Every rule **passes where its input is unmeasured** (NaN offset because the seam check could not run — no published images, insufficient overlap; NaN margin because no credible ceiling was fitted; unknown source image). That is the changelog's "no evidence either way" reading, and it matters: the seam check could not run on 92% of 2014_Greenland_P3 and 61% of 2012_Antarctica_DC8 at their flight geometry, so dropping unmeasured traces would remove whole seasons (35% of Antarctic and 50% of Greenland traces) rather than bad traces. The img2 rule is kept as suggested because img2-sourced traces have median RSSNR 8–17 dB below img1 traces of the same DC8 season, consistent with the documented surface-power low bias.
+
+The filter is a trace-level cleanup, not a season-level recalibration. Retraining with it (2026-09-01, `outputs/qc_filter/`, `uv run python scripts/qc_filter_comparison.py`) lowers the atten_refl spatial-CV RMSE from 13.02 to 12.87 dB, leaves held-out test RMSE within 0.05 dB when both posteriors are scored on the same points, and shifts full-grid predictions by less than 1 dB (5th–95th percentile −0.3 to +0.6 dB). The crossover matrices above are almost unchanged by it (`scripts/season_crossover_matrix.py --calibration-qc`): within-season scatter tightens for the seam-affected seasons, but the season-to-season offsets (e.g. 2012_Antarctica_DC8 ~12 dB below the later DC8 seasons) remain, and the QC alone catches only 2% of 2013_Greenland_P3 and 46% of 2017_Antarctica_Basler — so the two season exclusions stay.
+
 
 ## Do the surveys agree with each other?
 
@@ -68,7 +90,8 @@ cp outputs/model/analysis/residuals_by_season.png docs/figures/
 ```
 
 Note that this figure reflects the current configuration, in which the two
-outlier seasons above are already excluded.
+outlier seasons above are already excluded and the radiometric calibration QC
+is applied.
 
 There is certainly room for improvement here (or perhaps just further calibration), but there is no obvious pattern of dramatic outlier seasons or instruments.
 
diff --git a/docs/3_model.md b/docs/3_model.md
index fb600ce..114fd71 100644
--- a/docs/3_model.md
+++ b/docs/3_model.md
@@ -127,21 +127,21 @@ uv run python scripts/posterior_physical.py
 cp outputs/model/analysis/posterior_physical.png docs/figures/
 ```
 
-The headline accuracy and calibration numbers:
+The headline accuracy and calibration numbers (trained 2026-09-03 on the method-0.4.1 calibration snapshots with the radiometric calibration QC filter described in [Input data](2_input_data.md#radiometric-calibration-qc); 17,322 training grid points, 1,171 of them censored, plus 187 non-detections):
 
 | quantity | value |
 |---|---|
-| CV RMSE (5-fold, spatially blocked) | 13.02 dB (fold range 11.91–14.38) |
+| CV RMSE (5-fold, spatially blocked) | 12.87 dB (fold range 11.82–14.32) |
 | CV 1σ coverage | 0.68 |
-| Held-out test RMSE | 12.67 dB (n = 1,433 + 90 censored) |
-| Held-out test 1σ coverage | 0.72 |
-| Fully-linear baseline (same layers) | CV 13.92 dB / test 13.95 dB |
-| Sampler diagnostics | 0 divergences, R̂ ≤ 1.007 |
+| Held-out test RMSE | 12.72 dB (n = 1,330 + 88 censored) |
+| Held-out test 1σ coverage | 0.71 |
+| Fully-linear baseline (same layers) | CV 13.80 dB / test 13.98 dB |
+| Sampler diagnostics | 0 divergences, R̂ ≤ 1.005 |
 
 Posterior distributions of all 21 learned parameters, converted to physical units (the z-score normalization is an invertible affine transform, and the normalizer constants are stored in `posterior.nc`, so this conversion is exact). Attenuation-side parameters become two-way dB/km via σ_target/σ_thickness; reflectivity-side parameters become dB contributions to RSSNR (sign-flipped for the −refl convention); covariate effects are fully per-unit (e.g. dB/km/K, dB/km/(mW/m²)); θ, τ, and σ are natively in dB. Intercept-like values are referenced to the mean covariate conditions of the training set.
 
 ![Posterior distributions in physical units](figures/posterior_physical.png)
-*Posteriors in physical units, with the headline CV and held-out test RMSE. The 11.6 dB/km one-way (23.3 two-way) depth-averaged attenuation rate at mean conditions falls in the physically expected range. Posterior widths are small because n ≈ 21k; the meaningful uncertainty is the 13 dB residual σ. The 0/1 indicators are reported as the step between their states; the interaction panels are Greenland's *offset* from the Antarctic slope, not an absolute slope.*
+*Posteriors in physical units, with the headline CV and held-out test RMSE. The 11.6 dB/km one-way (23.2 two-way) depth-averaged attenuation rate at mean conditions falls in the physically expected range. Posterior widths are small because n ≈ 17k; the meaningful uncertainty is the 12.7 dB residual σ. The 0/1 indicators are reported as the step between their states; the interaction panels are Greenland's *offset* from the Antarctic slope, not an absolute slope.*
 
 The distribution of observed and posterior predicted RSSNR values are shown below by ice sheet:
 
@@ -166,3 +166,17 @@ Maps of the mean and 80th percentile predictions for both ice sheets are shown b
 
 ![80th percentile map](figures/map_q80.png)
 *The 80th percentile of the posterior predictive (recommended for instrument design). Note the different colorscale from the prior set of plots.*
+
+### Effect of the radiometric calibration QC
+
+The model above is the first trained after the upstream stores shipped per-trace seam-step and surface-saturation diagnostics, with the suggested filter applied (`split.calibration_qc`). The filter drops about 14% of Antarctic and 13% of Greenland traces but leaves the fit essentially unchanged: scored on the same held-out points, the pre- and post-QC posteriors differ by 0.03 dB in RMSE, every parameter stays within about two posterior standard deviations of its previous value, and the predicted maps move by less than 1 dB almost everywhere.
+
+![Prediction change from the calibration QC](figures/calibration_qc_prediction_difference.png)
+*Change in posterior-mean required surface SNR from adopting the calibration QC filter (new − previous model). Median +0.01 dB (Antarctica) and +0.07 dB (Greenland); 5th–95th percentile within ±0.6 dB.*
+
+The season-level offsets visible in the crossover matrices and in the out-of-fold residuals by season are not touched by a trace-level filter and remain the largest known calibration issue. To reproduce the comparison (requires a copy of the previous `outputs/model` and augment parquets at `outputs/baseline_20260807/`):
+
+```
+uv run python scripts/qc_filter_comparison.py
+cp outputs/qc_filter/prediction_difference.png docs/figures/calibration_qc_prediction_difference.png
+```
diff --git a/docs/figures/calibration_qc_prediction_difference.png b/docs/figures/calibration_qc_prediction_difference.png
new file mode 100644
index 0000000..9559939
Binary files /dev/null and b/docs/figures/calibration_qc_prediction_difference.png differ
diff --git a/docs/figures/cells_antarctic.png b/docs/figures/cells_antarctic.png
index 9a64b8f..05f78cd 100644
Binary files a/docs/figures/cells_antarctic.png and b/docs/figures/cells_antarctic.png differ
diff --git a/docs/figures/cells_greenland.png b/docs/figures/cells_greenland.png
index b1bfede..4c8228a 100644
Binary files a/docs/figures/cells_greenland.png and b/docs/figures/cells_greenland.png differ
diff --git a/docs/figures/hist_obs_vs_ppc_sheets.png b/docs/figures/hist_obs_vs_ppc_sheets.png
index 79cba41..77e1314 100644
Binary files a/docs/figures/hist_obs_vs_ppc_sheets.png and b/docs/figures/hist_obs_vs_ppc_sheets.png differ
diff --git a/docs/figures/map_pred_mean.png b/docs/figures/map_pred_mean.png
index 55667b8..1888c72 100644
Binary files a/docs/figures/map_pred_mean.png and b/docs/figures/map_pred_mean.png differ
diff --git a/docs/figures/map_q80.png b/docs/figures/map_q80.png
index 04e526c..c17f695 100644
Binary files a/docs/figures/map_q80.png and b/docs/figures/map_q80.png differ
diff --git a/docs/figures/posterior_physical.png b/docs/figures/posterior_physical.png
index bcc7e03..497f3c9 100644
Binary files a/docs/figures/posterior_physical.png and b/docs/figures/posterior_physical.png differ
diff --git a/docs/figures/region_histograms_ppc.png b/docs/figures/region_histograms_ppc.png
index bb9e736..e42f905 100644
Binary files a/docs/figures/region_histograms_ppc.png and b/docs/figures/region_histograms_ppc.png differ
diff --git a/docs/figures/residuals_by_season.png b/docs/figures/residuals_by_season.png
index 2620e50..20ae64f 100644
Binary files a/docs/figures/residuals_by_season.png and b/docs/figures/residuals_by_season.png differ
diff --git a/mission_design_tool/data/antarctic.bin.gz b/mission_design_tool/data/antarctic.bin.gz
index b75d1cd..ef49016 100644
Binary files a/mission_design_tool/data/antarctic.bin.gz and b/mission_design_tool/data/antarctic.bin.gz differ
diff --git a/mission_design_tool/data/coast.json.gz b/mission_design_tool/data/coast.json.gz
index e90c371..7c57447 100644
Binary files a/mission_design_tool/data/coast.json.gz and b/mission_design_tool/data/coast.json.gz differ
diff --git a/mission_design_tool/data/greenland.bin.gz b/mission_design_tool/data/greenland.bin.gz
index 698b90d..c52e7d0 100644
Binary files a/mission_design_tool/data/greenland.bin.gz and b/mission_design_tool/data/greenland.bin.gz differ
diff --git a/mission_design_tool/data/meta.json b/mission_design_tool/data/meta.json
index d2bcd79..f11f0e2 100644
--- a/mission_design_tool/data/meta.json
+++ b/mission_design_tool/data/meta.json
@@ -1,16 +1,16 @@
 {
   "model": "atten_refl",
-  "run_id": "85b6fe9d1c61",
-  "created_at": "2026-08-07T19:02:39.926018+00:00",
+  "run_id": "af5514ae3b72",
+  "created_at": "2026-09-03T16:27:20.771905+00:00",
   "stride": 1,
   "resolution_m": 5000,
   "mu_scale": 10.0,
   "sd_scale": 4.0,
   "detection": {
-    "theta_dB": -1.2173389585592402,
-    "tau_dB": 3.014524074451636
+    "theta_dB": -1.0778932306687852,
+    "tau_dB": 2.793973023761131
   },
-  "cv_rmse_dB": 13.017503319956557,
+  "cv_rmse_dB": 12.86812078837849,
   "features": [
     "bedmachine_thickness_m",
     "era5_t2m_mean_K",
@@ -39,7 +39,7 @@
       "crs": "EPSG:3031",
       "n": 527845,
       "bytes": 5278450,
-      "gz_bytes": 2266063
+      "gz_bytes": 2267374
     },
     "greenland": {
       "shape": [
@@ -53,7 +53,7 @@
       "crs": "EPSG:3413",
       "n": 68069,
       "bytes": 680690,
-      "gz_bytes": 292361
+      "gz_bytes": 293025
     }
   },
   "outlines": {
diff --git a/scripts/posterior_physical.py b/scripts/posterior_physical.py
index 39bca9c..75b0cd0 100644
--- a/scripts/posterior_physical.py
+++ b/scripts/posterior_physical.py
@@ -26,9 +26,9 @@
 COV_UNITS = {"era5_t2m_mean_K": "K", "surface_v_m_yr": "m/yr", "ghf_mW_m2": "mW/m²"}
 
 
-def main():
-    idata = az.from_netcdf(f"outputs/model/{MODEL}/posterior.nc")
-    metrics = json.load(open(f"outputs/model/{MODEL}/metrics.json"))
+def physical_panels(idata):
+    """(label, draws, unit) per learned parameter of an atten_refl posterior, in
+    physical units (see module docstring for the conversions)."""
     norm = json.loads(idata.attrs["normalizer"])
     sy = norm["required_surface_snr_dB"]["std"]
     st = norm["bedmachine_thickness_m"]["std"]
@@ -86,6 +86,13 @@ def step_note(c):
     panels.append(("σ — residual scatter", post["sigma"].values * sy, "dB"))
     panels.append(("θ — detection threshold", post["theta"].values, "dB"))
     panels.append(("τ — picker softness", post["tau"].values, "dB"))
+    return panels
+
+
+def main():
+    idata = az.from_netcdf(f"outputs/model/{MODEL}/posterior.nc")
+    metrics = json.load(open(f"outputs/model/{MODEL}/metrics.json"))
+    panels = physical_panels(idata)
 
     ncol = 5
     nrow = -(-(len(panels) + 1) // ncol)  # +1 leaves room for the accuracy panel
diff --git a/scripts/qc_filter_comparison.py b/scripts/qc_filter_comparison.py
new file mode 100644
index 0000000..6c1ffe0
--- /dev/null
+++ b/scripts/qc_filter_comparison.py
@@ -0,0 +1,385 @@
+"""Old-vs-new model comparison for the radiometric calibration QC filter.
+
+Baseline = a copy of outputs/ taken before the store re-pin + split.calibration_qc
+(model/ plus the antarctica/greenland augment parquets). New = outputs/ after
+re-running `snakemake model_all`. Writes to outputs/qc_filter/:
+  qc_removal_by_season.{csv,md,png}  per-season share of traces each rule rejects
+  qc_sensitivity.md                  trace share rejected under threshold variants
+  benchmark_comparison.md            metrics.json side by side, both models
+  same_points_test.md                old and new posteriors scored on the SAME
+                                     held-out test points (old and new test sets)
+  cv_folds.png                       per-fold CV RMSE, old vs new
+  posterior_comparison.png           atten_refl parameters in physical units
+  prediction_difference.png          new - old posterior-mean RSSNR maps + hist
+  target_hist.png                    training-target distribution before/after
+
+Usage: uv run python scripts/qc_filter_comparison.py [--baseline outputs/baseline_20260807]
+"""
+
+import argparse
+import json
+from pathlib import Path
+
+import arviz as az
+import matplotlib
+import numpy as np
+import pandas as pd
+import pyarrow.parquet as pq
+import xarray as xr
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt  # noqa: E402
+
+from radar_postproc.config import load_model_config  # noqa: E402
+from radar_postproc.models import get_model  # noqa: E402
+from radar_postproc.models.normalize import (  # noqa: E402
+    apply_normalizer, invert_normalizer, invert_scale)
+from radar_postproc.split import calibration_qc_flags  # noqa: E402
+from radar_postproc.train import (  # noqa: E402
+    _censored_mask, _design_matrix, _metrics, add_indicator_columns)
+
+from plot_style import C_ANT, C_GRL, INK, LS_OBS, SHEET_COLOR, style_axis  # noqa: E402
+from posterior_physical import physical_panels  # noqa: E402
+
+SHEETS = {"antarctic": "antarctica", "greenland": "greenland"}  # sheet -> store
+MODELS = ["linear", "atten_refl"]
+C_OLD, C_NEW = "0.55", "tab:orange"
+QC_COLS = ["img_comb_offset_dB", "surface_source_image_index", "surface_ceiling_margin_dB"]
+
+
+# --- trace-level QC accounting ------------------------------------------------
+
+def load_traces(out_dir: Path) -> pd.DataFrame:
+    frames = []
+    for sheet, store in SHEETS.items():
+        path = out_dir / store / f"{store}.parquet"
+        have = set(pq.read_schema(path).names)
+        df = pd.read_parquet(path, columns=[c for c in ["collection", *QC_COLS] if c in have])
+        for c in QC_COLS:  # pre-calibration parquets: nothing measured
+            if c not in df:
+                df[c] = np.nan
+        df["sheet"] = sheet
+        frames.append(df)
+    return pd.concat(frames, ignore_index=True)
+
+
+def removal_by_season(traces: pd.DataFrame, qc: dict, excluded: list[str]) -> pd.DataFrame:
+    flags = calibration_qc_flags(traces, qc)
+    d = traces[["sheet", "collection"]].copy()
+    d["seam"] = flags["seam"]
+    d["img2"] = flags["img2"] & ~flags["seam"]              # exclusive attribution
+    d["saturated"] = flags["saturated"] & ~flags["seam"] & ~flags["img2"]
+    d["any"] = flags.any(axis=1)
+    d["seam_unmeasured"] = traces["img_comb_offset_dB"].isna()
+    d["margin_unmeasured"] = traces["surface_ceiling_margin_dB"].isna()
+    g = d.groupby(["sheet", "collection"]).agg(
+        n=("any", "size"), **{c: (c, "mean") for c in
+                              ["seam", "img2", "saturated", "any",
+                               "seam_unmeasured", "margin_unmeasured"]}).reset_index()
+    g["excluded"] = g["collection"].isin(excluded)
+    return g.sort_values(["sheet", "collection"]).reset_index(drop=True)
+
+
+def plot_removal(tab: pd.DataFrame, qc: dict, out: Path):
+    tab = tab.iloc[::-1]
+    fig, ax = plt.subplots(figsize=(9, 0.42 * len(tab) + 1.8))
+    y = np.arange(len(tab))
+    for sheet_rows, color in ((tab["sheet"] == "antarctic", C_ANT),
+                              (tab["sheet"] == "greenland", C_GRL)):
+        r = tab[sheet_rows]
+        yy = y[sheet_rows.to_numpy()]
+        left = np.zeros(len(r))
+        for rule, alpha in (("seam", 1.0), ("img2", 0.6), ("saturated", 0.3)):
+            ax.barh(yy, 100 * r[rule], left=left, color=color, alpha=alpha,
+                    edgecolor="white", linewidth=0.5)
+            left += 100 * r[rule].to_numpy()
+    labels = [f"{c}{'  (excluded)' if e else ''}" for c, e in zip(tab["collection"], tab["excluded"])]
+    ax.set_yticks(y, labels, fontsize=8)
+    for yi, (_, row) in zip(y, tab.iterrows()):
+        ax.text(100 * row["any"] + 0.5, yi,
+                f"{100 * row['any']:.1f}%  (n={row['n']:,}; seam unmeasured "
+                f"{100 * row['seam_unmeasured']:.0f}%)", va="center", fontsize=7, color=INK)
+    ax.set_xlabel("traces rejected [%]", color=INK)
+    ax.set_xlim(0, max(60, 100 * tab["any"].max() + 28))
+    handles = [plt.Rectangle((0, 0), 1, 1, color="0.4", alpha=a) for a in (1.0, 0.6, 0.3)]
+    ax.legend(handles, [f"seam step |Δ| ≥ {qc['max_seam_offset_dB']:g} dB",
+                        "surface from higher-gain image (img2+)",
+                        f"surface within {qc['min_ceiling_margin_dB']:g} dB of season ceiling"],
+              loc="lower right", fontsize=8, frameon=False)
+    ax.set_title("Calibration QC: share of traces rejected per season (exclusive attribution)\n"
+                 "blue = Antarctica, green = Greenland; unmeasured seams/ceilings pass",
+                 color=INK, fontsize=10)
+    style_axis(ax)
+    fig.tight_layout()
+    fig.savefig(out, dpi=140, bbox_inches="tight")
+    plt.close(fig)
+
+
+def sensitivity_table(traces: pd.DataFrame, qc: dict) -> pd.DataFrame:
+    variants = [("adopted", {})]
+    variants += [(f"max_seam_offset_dB={v}", {"max_seam_offset_dB": v}) for v in (2.0, 5.0)]
+    variants += [(f"min_ceiling_margin_dB={v}", {"min_ceiling_margin_dB": v}) for v in (1.0, 3.0)]
+    variants += [("require_img1_surface=false", {"require_img1_surface": False}),
+                 ("drop_unmeasured_seam=true", {"drop_unmeasured_seam": True})]
+    rows = []
+    for name, over in variants:
+        flags = calibration_qc_flags(traces, {**qc, **over})
+        rej = flags.any(axis=1)
+        row = {"variant": name}
+        for sheet in SHEETS:
+            m = traces["sheet"] == sheet
+            row[f"{sheet}_rejected_pct"] = round(100 * rej[m].mean(), 1)
+        rows.append(row)
+    return pd.DataFrame(rows)
+
+
+# --- model-level comparison ---------------------------------------------------
+
+def load_run(model_dir: Path):
+    metrics = json.loads((model_dir / "metrics.json").read_text())
+    idata = az.from_netcdf(model_dir / "posterior.nc")
+    model = get_model(metrics["model"])
+    model.feature_names = json.loads(idata.attrs["features"])
+    return metrics, idata, model
+
+
+def benchmark_rows(runs: dict) -> pd.DataFrame:
+    rows = []
+    for (version, name), (m, idata, _) in runs.items():
+        cv, det = m["pooled_cv"], m["detection"]
+        post = idata.posterior
+        sy = json.loads(idata.attrs["normalizer"])["required_surface_snr_dB"]["std"]
+        rows.append({
+            "model": name, "version": version,
+            "n_train": cv["n_points"] + cv["n_censored"], "n_censored": cv["n_censored"],
+            "n_nondetect": det.get("n_nondetect_used"),
+            "cv_rmse_dB": f"{cv['rmse_dB']['mean']:.2f} [{cv['rmse_dB']['min']:.2f}–{cv['rmse_dB']['max']:.2f}]",
+            "cv_mae_dB": round(cv["mae_dB"]["mean"], 2),
+            "cv_cov1σ": round(cv["coverage_1sigma"]["mean"], 3),
+            "cv_logscore": round(cv["logscore_dB"]["mean"], 3),
+            "cv_nd_logscore": round(det["cv_nd_logscore"], 3) if "cv_nd_logscore" in det else None,
+            "test_rmse_dB": round(m["test"]["rmse_dB"], 2),
+            "test_logscore": round(m["test"]["logscore_dB"], 3),
+            "σ_dB": round(float(post["sigma"].values.mean() * sy), 2),
+            "θ_dB": round(det["theta_mean_dB"], 2), "τ_dB": round(det["tau_mean_dB"], 2),
+            "rhat_max": round(m["diagnostics"]["rhat_max"], 4),
+            "run_id": m["run_id"],
+        })
+    return pd.DataFrame(rows)
+
+
+def score_points(run, points: pd.DataFrame, target: str, censoring: dict) -> dict:
+    """Score a fitted posterior on given grid points (its own normalizer)."""
+    metrics, idata, model = run
+    norm = json.loads(idata.attrs["normalizer"])
+    df = points.copy()
+    add_indicator_columns(df, metrics["indicators"])
+    X = _design_matrix(df, metrics["features"], norm, tuple(metrics["indicators"]),
+                       tuple(metrics["interactions"]))
+    y = df[target].to_numpy()
+    mu, _, pstd = model.predict(idata, X)
+    mu_dB = invert_normalizer(mu, norm[target])
+    std_dB = invert_scale(pstd, norm[target])
+    cens = _censored_mask(df, censoring)
+    out = _metrics(y[~cens], mu_dB[~cens], std_dB[~cens])
+    out["bias_dB"] = float(np.mean(mu_dB[~cens] - y[~cens]))
+    out["logscore_dB"] = (model.logscore(idata, X, apply_normalizer(y, norm[target]), cens)
+                          - float(np.log(norm[target]["std"])))
+    out["n_censored"] = int(cens.sum())
+    return out
+
+
+def usable_points(split: pd.DataFrame, cfg: dict, mask) -> pd.DataFrame:
+    tcfg = cfg["train"]
+    ok = split[cfg["split"]["target"]].notna() & split[tcfg["features"]].notna().all(axis=1)
+    if tcfg["min_thickness_m"] is not None:
+        ok &= ~(split["bedmachine_thickness_m"] < tcfg["min_thickness_m"])
+    return split[ok & mask]
+
+
+def plot_cv_folds(runs: dict, out: Path):
+    fig, axes = plt.subplots(1, len(MODELS), figsize=(5 * len(MODELS), 3.6), sharey=True)
+    for ax, name in zip(np.atleast_1d(axes), MODELS):
+        for i, (version, color) in enumerate((("baseline", C_OLD), ("QC", C_NEW))):
+            folds = runs[(version, name)][0]["folds"]
+            k = np.array([f["fold"] for f in folds])
+            ax.bar(k + (i - 0.5) * 0.38, [f["rmse_dB"] for f in folds], width=0.38,
+                   color=color, label=f"{version} (pooled {runs[(version, name)][0]['pooled_cv']['rmse_dB']['mean']:.2f} dB)")
+        ax.set_title(name, color=INK)
+        ax.set_xlabel("spatial CV fold", color=INK)
+        ax.legend(fontsize=8, frameon=False)
+        style_axis(ax)
+    np.atleast_1d(axes)[0].set_ylabel("out-of-fold RMSE [dB]", color=INK)
+    fig.suptitle("Spatially-blocked CV RMSE per fold — before (grey) vs after (orange) calibration QC",
+                 color=INK, fontsize=10)
+    fig.tight_layout()
+    fig.savefig(out, dpi=140, bbox_inches="tight")
+    plt.close(fig)
+
+
+def plot_posteriors(old, new, out: Path):
+    po, pn = physical_panels(old), physical_panels(new)
+    ncol = 5
+    nrow = -(-len(po) // ncol)
+    fig, axes = plt.subplots(nrow, ncol, figsize=(3.3 * ncol, 2.9 * nrow))
+    for ax, (label, d_old, unit), (_, d_new, _) in zip(axes.ravel(), po, pn):
+        lo, hi = np.percentile(np.concatenate([d_old, d_new]), [0.1, 99.9])
+        bins = np.linspace(lo, hi, 60)
+        ax.hist(d_old, bins=bins, density=True, color=C_OLD, alpha=0.35)
+        ax.hist(d_old, bins=bins, density=True, histtype="step", color=C_OLD, linewidth=1.3)
+        ax.hist(d_new, bins=bins, density=True, color=C_NEW, alpha=0.35)
+        ax.hist(d_new, bins=bins, density=True, histtype="step", color=C_NEW, linewidth=1.5)
+        ax.set_title(f"{label}\n{d_old.mean():+.3g} ± {d_old.std():.2g} → "
+                     f"{d_new.mean():+.3g} ± {d_new.std():.2g} {unit}", color=INK, fontsize=8.5)
+        style_axis(ax)
+        ax.set_yticks([])
+    for ax in axes.ravel()[len(po):]:
+        ax.axis("off")
+    axes.ravel()[0].legend([plt.Rectangle((0, 0), 1, 1, color=C_OLD, alpha=0.5),
+                            plt.Rectangle((0, 0), 1, 1, color=C_NEW, alpha=0.5)],
+                           ["baseline", "with calibration QC"], fontsize=8, frameon=False)
+    fig.suptitle("atten_refl posteriors in physical units — baseline (grey) vs calibration QC (orange)\n"
+                 "titles: baseline → QC, mean ± sd", color=INK, fontsize=11)
+    fig.tight_layout()
+    fig.savefig(out, dpi=140, bbox_inches="tight")
+    plt.close(fig)
+
+
+def plot_prediction_difference(old_zarr: Path, new_zarr: Path, out: Path) -> pd.DataFrame:
+    diffs, stats = {}, []
+    for sheet in SHEETS:
+        a = xr.open_zarr(old_zarr, group=sheet)
+        b = xr.open_zarr(new_zarr, group=sheet)
+        d = (b["pred_mean"] - a["pred_mean"]).load()
+        diffs[sheet] = (d.values, d["x"].values, d["y"].values)
+        v = d.values[np.isfinite(d.values)]
+        s = (b["pred_std"] / a["pred_std"]).values
+        stats.append({"sheet": sheet, "n_grid": int(v.size),
+                      "median_shift_dB": round(float(np.median(v)), 2),
+                      "p5_dB": round(float(np.percentile(v, 5)), 2),
+                      "p95_dB": round(float(np.percentile(v, 95)), 2),
+                      "median_pred_std_ratio": round(float(np.nanmedian(s)), 3)})
+    lim = max(np.nanpercentile(np.abs(v[0]), 99) for v in diffs.values())
+    fig = plt.figure(figsize=(16, 7))
+    gs = fig.add_gridspec(1, 3, width_ratios=[1.1, 0.75, 0.8])
+    axes = [fig.add_subplot(gs[0, i]) for i in range(3)]
+    for ax, sheet in zip(axes[:2], SHEETS):
+        arr, x, y = diffs[sheet]
+        im = ax.imshow(arr, cmap="PRGn", vmin=-lim, vmax=lim,
+                       extent=[x[0], x[-1], y[-1], y[0]])
+        ax.set_aspect("equal")
+        ax.set_axis_off()
+        ax.set_title(sheet, color=INK)
+    fig.colorbar(im, ax=axes[:2], shrink=0.7, label="new − baseline posterior-mean RSSNR [dB]")
+    ax = axes[2]
+    bins = np.linspace(-lim, lim, 80)
+    for sheet in SHEETS:
+        v = diffs[sheet][0]
+        v = v[np.isfinite(v)]
+        ax.hist(v, bins=bins, density=True, histtype="step", linewidth=1.5,
+                color=SHEET_COLOR[sheet], linestyle=LS_OBS,
+                label=f"{sheet}: median {np.median(v):+.2f} dB")
+    ax.axvline(0, color="0.4", linewidth=0.8)
+    ax.set_xlabel("new − baseline predicted RSSNR [dB]", color=INK)
+    ax.set_yticks([])
+    ax.legend(fontsize=8, frameon=False)
+    style_axis(ax)
+    fig.suptitle("Change in predicted required surface SNR from the calibration QC filter (atten_refl)",
+                 color=INK, fontsize=11)
+    fig.savefig(out, dpi=140, bbox_inches="tight")
+    plt.close(fig)
+    return pd.DataFrame(stats)
+
+
+def plot_target_hist(old_split, new_split, target: str, out: Path):
+    fig, ax = plt.subplots(figsize=(7.5, 4))
+    bins = np.linspace(-20, 120, 71)
+    for sheet in SHEETS:
+        for split, lw, alpha, tag in ((old_split, 1.0, 0.5, "baseline"), (new_split, 1.8, 1.0, "QC")):
+            v = split.loc[(split["ice_sheet"] == sheet) & (split["fold"] >= 0), target].dropna()
+            ax.hist(v, bins=bins, histtype="step", color=SHEET_COLOR[sheet], linewidth=lw,
+                    alpha=alpha, linestyle=LS_OBS,
+                    label=f"{sheet} {tag}: n={len(v):,}, median {v.median():.1f} dB")
+    ax.set_xlabel("training target: required surface SNR [dB]", color=INK)
+    ax.set_ylabel("grid points", color=INK)
+    ax.legend(fontsize=8, frameon=False)
+    ax.set_title("Training targets before (thin) and after (thick) calibration QC", color=INK, fontsize=10)
+    style_axis(ax)
+    fig.tight_layout()
+    fig.savefig(out, dpi=140, bbox_inches="tight")
+    plt.close(fig)
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--baseline", default="outputs/baseline_20260807")
+    ap.add_argument("--new", default="outputs")
+    ap.add_argument("--out", default="outputs/qc_filter")
+    args = ap.parse_args()
+    base, new, out = Path(args.baseline), Path(args.new), Path(args.out)
+    out.mkdir(parents=True, exist_ok=True)
+    cfg = load_model_config("config/model.yaml")
+    qc, target = cfg["split"]["calibration_qc"], cfg["split"]["target"]
+    md = ["# Calibration QC filter: baseline vs new\n",
+          f"QC config: `{json.dumps(qc)}`; excluded collections: "
+          f"{cfg['split']['exclude_collections']}\n"]
+
+    # 1. trace-level accounting (new parquets carry the calibration columns)
+    traces = load_traces(new)
+    tab = removal_by_season(traces, qc, cfg["split"]["exclude_collections"])
+    tab.to_csv(out / "qc_removal_by_season.csv", index=False)
+    show = tab.copy()
+    for c in ["seam", "img2", "saturated", "any", "seam_unmeasured", "margin_unmeasured"]:
+        show[c] = (100 * show[c]).round(1)
+    (out / "qc_removal_by_season.md").write_text(show.to_markdown(index=False) + "\n")
+    plot_removal(tab, qc, out / "qc_removal_by_season.png")
+    sens = sensitivity_table(traces, qc)
+    (out / "qc_sensitivity.md").write_text(sens.to_markdown(index=False) + "\n")
+    md += ["## Traces rejected per season (%)\n", show.to_markdown(index=False), "",
+           "## Threshold sensitivity (% of traces rejected)\n", sens.to_markdown(index=False), ""]
+
+    # 2. metrics side by side
+    runs = {(v, m): load_run(d / "model" / m)
+            for v, d in (("baseline", base), ("QC", new)) for m in MODELS}
+    bench = benchmark_rows(runs)
+    (out / "benchmark_comparison.md").write_text(bench.to_markdown(index=False) + "\n")
+    md += ["## Benchmark (each run scored on its own CV folds / test cells)\n",
+           bench.to_markdown(index=False), ""]
+    plot_cv_folds(runs, out / "cv_folds.png")
+
+    # 3. same held-out points, both posteriors
+    old_split = pd.read_parquet(base / "model" / "split.parquet")
+    new_split = pd.read_parquet(new / "model" / "split.parquet")
+    assert len(old_split) == len(new_split), "grid changed between runs"
+    rows = []
+    for set_name, split in (("baseline test set", old_split), ("QC test set", new_split)):
+        pts = usable_points(split, cfg, split["is_test"])
+        for (version, name), run in runs.items():
+            s = score_points(run, pts, target, cfg["train"]["censoring"])
+            rows.append({"test points": set_name, "n": len(pts), "model": name,
+                         "posterior": version, **{k: round(v, 3) for k, v in s.items()}})
+    same = pd.DataFrame(rows)
+    (out / "same_points_test.md").write_text(same.to_markdown(index=False) + "\n")
+    md += ["## Same test points, both posteriors (uncensored RMSE/MAE/coverage; "
+           "log score over all points)\n", same.to_markdown(index=False), ""]
+
+    # 4. posteriors, predictions, targets
+    plot_posteriors(runs[("baseline", "atten_refl")][1], runs[("QC", "atten_refl")][1],
+                    out / "posterior_comparison.png")
+    shift = plot_prediction_difference(base / "model/atten_refl/predictions.zarr",
+                                       new / "model/atten_refl/predictions.zarr",
+                                       out / "prediction_difference.png")
+    md += ["## Prediction shift (atten_refl, full grid)\n", shift.to_markdown(index=False), ""]
+    plot_target_hist(old_split, new_split, target, out / "target_hist.png")
+    n_old = int((old_split[target].notna() & (old_split["fold"] >= 0)).sum())
+    n_new = int((new_split[target].notna() & (new_split["fold"] >= 0)).sum())
+    changed = int(((old_split[target] != new_split[target]) & old_split[target].notna()
+                   & new_split[target].notna()).sum())
+    md += [f"Training grid points: {n_old} → {n_new}; {changed} kept points now match a "
+           "different (QC-passing) trace.\n"]
+    (out / "comparison.md").write_text("\n".join(md))
+    print("\n".join(md))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/residual_audit.py b/scripts/residual_audit.py
index 30307f6..4c99fc2 100644
--- a/scripts/residual_audit.py
+++ b/scripts/residual_audit.py
@@ -102,9 +102,11 @@ def main():
                for c in C_INST.values()]
     ax.legend(handles, C_INST.keys(), frameon=False, fontsize=9)
     excl = config["split"]["exclude_collections"]
+    qc_on = config["split"]["calibration_qc"]["enabled"]
     ax.set_title(f"Out-of-fold residuals by season (uncensored points; {MODEL}, "
                  f"indicators={list(indicators)}"
-                 + (f"; excluded: {', '.join(excl)}" if excl else "") + ")",
+                 + (f"; excluded: {', '.join(excl)}" if excl else "")
+                 + ("; calibration QC on" if qc_on else "") + ")",
                  color=INK, fontsize=11)
     fig.tight_layout()
     out = "outputs/model/analysis/residuals_by_season.png"
diff --git a/scripts/season_crossover_matrix.py b/scripts/season_crossover_matrix.py
index 9bc1d90..637c433 100644
--- a/scripts/season_crossover_matrix.py
+++ b/scripts/season_crossover_matrix.py
@@ -9,9 +9,13 @@
 Rendered as one annotated matrix per sheet: cell colour = median (PRGn
 diverging, centred at 0), gray = fewer than MIN_PAIRS pairs.
 
-Usage: uv run python scripts/season_crossover_matrix.py
+Usage: uv run python scripts/season_crossover_matrix.py [--calibration-qc]
+  --calibration-qc applies config/model.yaml's split.calibration_qc rules to the
+  traces first (output suffix _qc), for a model-free before/after check of
+  inter-season agreement.
 """
 
+import argparse
 from pathlib import Path
 
 import matplotlib
@@ -23,6 +27,9 @@
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt  # noqa: E402
 
+from radar_postproc.config import load_model_config  # noqa: E402
+from radar_postproc.split import apply_calibration_qc  # noqa: E402
+
 from plot_style import INK  # noqa: E402
 
 RADIUS_M = 500.0
@@ -33,8 +40,11 @@
           "greenland": ("outputs/greenland/greenland.parquet", "EPSG:3413")}
 
 
-def load_sheet(path: str, crs: str) -> pd.DataFrame:
+def load_sheet(path: str, crs: str, qc: dict | None = None) -> pd.DataFrame:
     df = pd.read_parquet(path)
+    if qc is not None:
+        df, counts = apply_calibration_qc(df, qc)
+        print(f"{path}: calibration QC dropped {counts['any']}/{counts['n_before']} traces")
     noise = df.get("post_bed_noise_interp_dB")
     if noise is None or noise.isna().all():
         noise = df["post_bed_noise_dB"]
@@ -71,9 +81,18 @@ def pair_deltas(a: pd.DataFrame, b: pd.DataFrame, same_season: bool) -> np.ndarr
 
 
 def main():
-    out_dir = Path("outputs/model/analysis")
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--calibration-qc", action="store_true")
+    ap.add_argument("--out-dir", default="outputs/model/analysis")
+    args = ap.parse_args()
+    qc = None
+    if args.calibration_qc:
+        qc = {**load_model_config("config/model.yaml")["split"]["calibration_qc"], "enabled": True}
+    suffix = "_qc" if qc else ""
+    out_dir = Path(args.out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
     for sheet, (path, crs) in SHEETS.items():
-        df = load_sheet(path, crs)
+        df = load_sheet(path, crs, qc)
         seasons = sorted(df["collection"].unique())
         n = len(seasons)
         med = np.full((n, n), np.nan)
@@ -107,17 +126,23 @@ def main():
         short = [s.split("_", 1)[0] + " " + s.split("_")[-1] for s in seasons]
         ax.set_xticks(range(n), short, rotation=45, ha="right", fontsize=9)
         ax.set_yticks(range(n), short, fontsize=9)
-        ax.set_title(f"{sheet}: RSSNR at season crossovers — median(row − col) ± sd [dB]\n"
+        ax.set_title(f"{sheet}: RSSNR at season crossovers{' (calibration QC applied)' if qc else ''}"
+                     " — median(row − col) ± sd [dB]\n"
                      f"pairs within {RADIUS_M:.0f} m, both margins > {MIN_MARGIN_DB:.0f} dB, "
                      f"min {MIN_PAIRS} pairs; diagonal = within-season (cross-frame)",
                      color=INK, fontsize=11)
         fig.colorbar(im, ax=ax, shrink=0.8, label="median Δ RSSNR [dB]")
         fig.tight_layout()
-        out = out_dir / f"season_crossover_matrix_{sheet}.png"
+        out = out_dir / f"season_crossover_matrix_{sheet}{suffix}.png"
         fig.savefig(out, dpi=140, bbox_inches="tight")
         plt.close(fig)
         total_pairs = int(cnt[np.triu_indices(n)].sum())
         print(f"{sheet}: {n} seasons, {total_pairs} pairs total -> {out}")
+        # Long-format table of the upper triangle (incl. diagonal) for numeric summaries.
+        pd.DataFrame([{"row": seasons[i], "col": seasons[j], "median_dB": med[i, j],
+                       "sd_dB": sd[i, j], "n_pairs": int(cnt[i, j])}
+                      for i in range(n) for j in range(i, n)]
+                     ).to_csv(out_dir / f"season_crossover_pairs_{sheet}{suffix}.csv", index=False)
 
 
 if __name__ == "__main__":
diff --git a/src/radar_postproc/config.py b/src/radar_postproc/config.py
index 36e5b89..aaf91f9 100644
--- a/src/radar_postproc/config.py
+++ b/src/radar_postproc/config.py
@@ -67,6 +67,12 @@ def load_config(config_path: str | Path) -> dict:
             "bed_pick_attempted",
             "qc_surface_pass",
             "qc_pass",
+            # Radiometric calibration diagnostics (upstream method 0.4.1; see
+            # docs/dataset_changelog.md there). Consumed by split.calibration_qc.
+            "img_comb_offset_dB",
+            "img_comb_pair",
+            "surface_source_image_index",
+            "surface_ceiling_margin_dB",
         ],
     )
 
@@ -118,6 +124,18 @@ def load_model_config(config_path: str | Path) -> dict:
     split.setdefault("test_cells", [])  # e.g. ["ant:-3:1"]; empty -> warning, no test set
     # Collections dropped entirely before matching (crossover-identified outliers).
     split.setdefault("exclude_collections", [])
+    # Per-trace radiometric QC from the upstream calibration fields (image-combine
+    # seam steps + surface saturation). Applied to observations AND non-detections
+    # before matching. Every rule passes where its input is NaN / unknown: the
+    # check did not run there, which is not evidence either way (dropping
+    # unmeasured seams would remove whole seasons, e.g. 93% of 2014_Greenland_P3).
+    split.setdefault("calibration_qc", {})
+    qc = split["calibration_qc"]
+    qc.setdefault("enabled", False)
+    qc.setdefault("max_seam_offset_dB", 3.0)     # drop |img_comb_offset_dB| >= this
+    qc.setdefault("drop_unmeasured_seam", False)  # also drop NaN offsets
+    qc.setdefault("require_img1_surface", True)   # drop surface_source_image_index >= 2
+    qc.setdefault("min_ceiling_margin_dB", 2.0)   # drop surface_ceiling_margin_dB < this
 
     # Bayesian model training / prediction.
     config.setdefault("train", {})
diff --git a/src/radar_postproc/split.py b/src/radar_postproc/split.py
index 398537c..b05ab7b 100644
--- a/src/radar_postproc/split.py
+++ b/src/radar_postproc/split.py
@@ -126,7 +126,8 @@ def compute_ceiling(surface_power_dB, surface_twtt, noise_dB, thickness_m,
 _OBS_COLS = ["latitude", "longitude", "collection", "bed_power_dB", "post_bed_noise_dB",
              "post_bed_noise_interp_dB", "post_bed_peak_interp_dB",
              "surface_power_dB", "surface_twtt",
-             "bed_pick_available", "bed_pick_attempted"]
+             "bed_pick_available", "bed_pick_attempted",
+             "img_comb_offset_dB", "surface_source_image_index", "surface_ceiling_margin_dB"]
 
 _UTIG_SUFFIXES = ("_BaslerJKB", "_BaslerMKB")
 
@@ -138,8 +139,45 @@ def institution_of(collection) -> str | None:
     return "UTIG" if collection.endswith(_UTIG_SUFFIXES) else "CReSIS"
 
 
+def calibration_qc_flags(df: pd.DataFrame, qc: dict) -> pd.DataFrame:
+    """Per-rule rejection flags from the upstream radiometric calibration fields.
+
+    Columns (True = reject): seam (image-combine step |img_comb_offset_dB| >=
+    max_seam_offset_dB; NaN rejected only with drop_unmeasured_seam), img2
+    (surface sampled from a higher-gain image, index >= 2 — likely saturated,
+    season-dependent low bias; -1/NaN unknown passes), saturated
+    (surface_ceiling_margin_dB < min_ceiling_margin_dB; NaN = no credible season
+    ceiling, passes). Missing columns (older stores) reject nothing.
+    """
+    n = len(df)
+    nan = pd.Series(np.nan, index=df.index)
+    seam = df.get("img_comb_offset_dB", nan).astype("float64")
+    src = df.get("surface_source_image_index", nan).astype("float64")
+    margin = df.get("surface_ceiling_margin_dB", nan).astype("float64")
+    flags = pd.DataFrame(index=df.index)
+    flags["seam"] = (seam.abs() >= qc["max_seam_offset_dB"]).to_numpy()
+    if qc["drop_unmeasured_seam"]:
+        flags["seam"] |= seam.isna().to_numpy()
+    flags["img2"] = (src >= 2).to_numpy() if qc["require_img1_surface"] else np.zeros(n, bool)
+    flags["saturated"] = (margin < qc["min_ceiling_margin_dB"]).to_numpy()
+    return flags
+
+
+def apply_calibration_qc(obs: pd.DataFrame, qc: dict) -> tuple[pd.DataFrame, dict]:
+    """Drop traces failing any calibration_qc rule; return (kept, counts by rule)."""
+    if not qc["enabled"]:
+        return obs, {}
+    flags = calibration_qc_flags(obs, qc)
+    reject = flags.any(axis=1).to_numpy()
+    counts = {k: int(v) for k, v in flags.sum().items()}
+    counts["any"] = int(reject.sum())
+    counts["n_before"] = len(obs)
+    return obs[~reject].reset_index(drop=True), counts
+
+
 def _load_observations(sheet: str, stores: list[str], target: str, out_dir: Path,
-                       exclude_collections: tuple = ()) -> pd.DataFrame:
+                       exclude_collections: tuple = (),
+                       calibration_qc: dict | None = None) -> tuple[pd.DataFrame, dict]:
     """Pooled attempted radar traces for a sheet, in its native projected CRS.
 
     Rows are picked observations OR non-detections (attempted, no bed pick —
@@ -170,6 +208,14 @@ def _load_observations(sheet: str, stores: list[str], target: str, out_dir: Path
         obs = obs[~obs["collection"].isin(list(exclude_collections))].reset_index(drop=True)
         logger.info("Split %s: excluded %d traces from collections %s",
                     sheet, n0 - len(obs), list(exclude_collections))
+    qc_counts = {}
+    if calibration_qc is not None:
+        obs, qc_counts = apply_calibration_qc(obs, calibration_qc)
+        if qc_counts:
+            logger.info("Split %s: calibration QC dropped %d/%d traces (seam %d, img2 "
+                        "surface %d, saturated %d)", sheet, qc_counts["any"],
+                        qc_counts["n_before"], qc_counts["seam"], qc_counts["img2"],
+                        qc_counts["saturated"])
 
     # Old stores carry no pick flags: rows with a target are picked observations.
     avail = obs["bed_pick_available"].astype("float64")
@@ -185,7 +231,7 @@ def _load_observations(sheet: str, stores: list[str], target: str, out_dir: Path
     obs["delta_dB"] = obs["post_bed_peak_interp_dB"] - obs["post_bed_noise_interp_dB"]
     tx = Transformer.from_crs("EPSG:4326", REGION_CRS[sheet], always_xy=True)
     obs["x"], obs["y"] = tx.transform(obs["longitude"].to_numpy(), obs["latitude"].to_numpy())
-    return obs
+    return obs, qc_counts
 
 
 def _cells_summary(df: pd.DataFrame, target: str, cell_size_m: float) -> pd.DataFrame:
@@ -269,11 +315,14 @@ def run_split(config_path: str, out_dir: str | None = None, repo_dir: str = ".")
     ceiling = np.full(len(df), np.nan)
     collection = np.full(len(df), None, dtype=object)
     institution = np.full(len(df), None, dtype=object)
+    qc_summary: dict[str, dict] = {}
     for sheet, stores in config["inputs"].items():
         for store in stores:
             augment_run_ids[store] = read_run_id(out_dir / store / f"{store}.parquet")
-        obs = _load_observations(sheet, stores, target, out_dir,
-                                 exclude_collections=tuple(split_cfg["exclude_collections"]))
+        obs, qc_summary[sheet] = _load_observations(
+            sheet, stores, target, out_dir,
+            exclude_collections=tuple(split_cfg["exclude_collections"]),
+            calibration_qc=split_cfg["calibration_qc"])
         rows = (df["ice_sheet"] == sheet).to_numpy()
         # Nearest *attempted* trace: a grid point takes whichever attempted trace
         # is closest — a picked observation or a non-detection.
@@ -378,6 +427,7 @@ def run_split(config_path: str, out_dir: str | None = None, repo_dir: str = ".")
         folds={str(k): int(v) for k, v in fold_sizes.items()},
         n_train_points=n_train,
         n_test_points=n_test,
+        calibration_qc=qc_summary,
     )
     paths = write_stage_output(df, manifest, model_dir / "split.parquet")
     paths["cells"] = str(model_dir / "cells.csv")
diff --git a/tests/unit/test_calibration_qc.py b/tests/unit/test_calibration_qc.py
new file mode 100644
index 0000000..f04f70b
--- /dev/null
+++ b/tests/unit/test_calibration_qc.py
@@ -0,0 +1,63 @@
+"""Unit tests for the split-stage radiometric calibration QC rules."""
+
+import numpy as np
+import pandas as pd
+
+from radar_postproc.config import load_model_config
+from radar_postproc.split import apply_calibration_qc, calibration_qc_flags
+
+
+def _qc(**overrides):
+    qc = {"enabled": True, "max_seam_offset_dB": 3.0, "drop_unmeasured_seam": False,
+          "require_img1_surface": True, "min_ceiling_margin_dB": 2.0}
+    qc.update(overrides)
+    return qc
+
+
+def _obs():
+    return pd.DataFrame({
+        #                        clean  seam+  seam-  seamNaN  img2  img-1  sat   marginNaN
+        "img_comb_offset_dB":      [0.5, 3.0,  -4.0,  np.nan,  1.0,  1.0,  1.0,  1.0],
+        "surface_source_image_index": [1, 1,    1,     1,       2,    -1,   1,    1],
+        "surface_ceiling_margin_dB": [5.0, 5.0, 5.0,   5.0,     5.0,  5.0,  1.9,  np.nan],
+    })
+
+
+def test_rules_flag_expected_rows():
+    flags = calibration_qc_flags(_obs(), _qc())
+    assert flags["seam"].tolist() == [False, True, True, False, False, False, False, False]
+    assert flags["img2"].tolist() == [False, False, False, False, True, False, False, False]
+    assert flags["saturated"].tolist() == [False] * 6 + [True, False]
+
+
+def test_unmeasured_passes_by_default_and_can_be_dropped():
+    kept, counts = apply_calibration_qc(_obs(), _qc())
+    assert len(kept) == 4 and counts == {"seam": 2, "img2": 1, "saturated": 1,
+                                         "any": 4, "n_before": 8}
+    kept, counts = apply_calibration_qc(_obs(), _qc(drop_unmeasured_seam=True))
+    assert len(kept) == 3 and counts["seam"] == 3
+
+
+def test_img1_rule_optional():
+    kept, _ = apply_calibration_qc(_obs(), _qc(require_img1_surface=False))
+    assert len(kept) == 5
+
+
+def test_missing_columns_reject_nothing():
+    obs = pd.DataFrame({"bed_power_dB": [1.0, 2.0]})
+    kept, counts = apply_calibration_qc(obs, _qc())
+    assert len(kept) == 2 and counts["any"] == 0
+
+
+def test_disabled_is_passthrough():
+    obs = _obs()
+    kept, counts = apply_calibration_qc(obs, _qc(enabled=False))
+    assert kept is obs and counts == {}
+
+
+def test_config_defaults(tmp_path):
+    path = tmp_path / "m.yaml"
+    path.write_text("split: {calibration_qc: {enabled: true}}\n")
+    qc = load_model_config(path)["split"]["calibration_qc"]
+    assert qc == {"enabled": True, "max_seam_offset_dB": 3.0, "drop_unmeasured_seam": False,
+                  "require_img1_surface": True, "min_ceiling_margin_dB": 2.0}