# Error check — THEMIS hot-spring project (2026-09-07)

Fresh-eyes audit of the analysis code, the numbers in REPORT.md / literature_review.md /
the PDF, all 22 figures, and notes/reconstruction_notes.md. Four independent audit passes
(physics scripts; statistics and counts; candidate/day-night chain; figures, captions and
unsupported claims) followed by targeted reruns and repairs.

Severity: **critical** = a stated conclusion was unsupported; **major** = a reported number,
method description or figure was wrong but the conclusion survives; **minor** = small numeric
drift or a missing qualifier; **cosmetic** = wording / plotting.

Published deliverables are preserved untouched: `THEMIS_hot_springs_report_v1_2026-07-06.pdf`
(original), `THEMIS_hot_springs_report_v2_2026-09-07_precheck.pdf` and
`REPORT_v2_2026-09-07_precheck.md` (as reconstructed, before this check),
`literature_review_v1_2026-06-18.md`. Corrected versions carry the base names; corrected
figures carry a `_v2` suffix beside the originals in `data/analysis/` and `figures/`.

## Headline: did any conclusion change?

**The top-line result stands.** No compact, localized, hot feature consistent with an
active hot spring exists in the nighttime archive; the Bayesian limit (< ~2 vigorous vents,
95%) is unchanged; the depth-to-liquid-water distribution is unchanged and, if anything,
conservative.

**One supporting claim was unsupported as published and has been redone.** The
per-candidate day/night and thermal-model results of §10–§11 ("every one of 21,385 persistent
candidates individually resolved as rock", "95th-percentile endogenic warming < 3 K") were
computed from samples taken a median of 3 pixels away from the detected warm pixels, i.e.
from near-candidate background. See finding C1 for the corrected analysis and its numbers.

**Two other analysis bugs were fixed and rerun:** the regional (§4) cluster positions were
flipped north–south (C2), and the global mosaic's TES albedo layer was misregistered by 180°
of longitude (C3). Neither changes the null result.

**Bookkeeping was wrong in several places:** 194 of the 828 "pixel-analysed" warmest-pixel
outliers were empty byte-fill products (B1); five different image counts were used
interchangeably (B4); the global noise figure was stale (B5).

## A. Physics scripts (detection_limit, detection_depth, depth_to_water, depth_atlas, thermal_model)

| id | severity | finding | fix | conclusion changed? |
|---|---|---|---|---|
| A1 | cosmetic | REPORT §3 table, 1200 K row, 5 K column said 9 m²; script gives 6.5 m² (8.7 is the 1000 K value) | table corrected | no |
| A2 | minor | §6 "~10⁵ GW globally"; script gives 3.0×10⁵ GW | text → ~3×10⁵ GW | no |
| A3 | minor | §6/§8b cooling time "40–400 yr", "decades": 2.1 W m⁻² is the *ceiling* of invisible flux (so a lower bound on lifetime) and latent heat of freezing (3.3×10⁸ J m⁻³) exceeds the sensible heat, lengthening it ~2–3× | text reworded ("at least", "centuries") | no |
| A4 | cosmetic | "×80 / ×100 / ×50–100" background multipliers used different heat-flow bases (25 vs 21 vs 20 mW m⁻²) | harmonized to ~×100 / ~40–80× with "~20–25 mW m⁻²" | no |
| A5 | minor | `depth_to_water.py` docstring claims a ~213 K equatorial anchor; the code yields 220.4 K (equilibrium T of the mean flux exceeds mean T by ~5–8 K, Jensen bias). Rerun with −7 K: pure-water median 7.2 → 8.0 km, marginalized 5.4 → 6.2 km, P(marg < 1 m) 2.6e-2 → 1.3e-2, P(pure < 1 km) stays 0. Bias makes the published depths *shallower* and the tail *fatter*, i.e. conservative for the report's claims | docstring corrected; caveat noted in §8b | no |
| A6 | minor | §9 prose "heat flow high over Tharsis": the atlas has Tharsis at the mean (19.9 vs 20.0 mW m⁻²); highs are Acidalia 22.9 / Arabia 21.1, lows Hellas 17.9 / Elysium 18.8 | text corrected | no |
| A7 | cosmetic | `depth_atlas.py` K heat-production constant 16% low (2.92e-9 vs 3.48e-9 W/kg); absorbed by the 20 mW m⁻² calibration (0.1 mW m⁻² effect) | documented only | no |
| A8 | **major (model limitation)** | `thermal_model.py` has no atmospheric back-radiation and no CO₂-frost floor: equatorial nights 10–40 K too cold at low inertia (128 K at I=50), amplitudes 25–50% too large. §11 differences each candidate against its own in-image background, which cancels the absolute bias; sensitivity rerun with F_down = 15 W m⁻² changes the §11 percentiles by < 0.1 K (95th 3.00 → 3.00 K, 99th 5.19 → 5.14 K) | optional `F_down` argument added (default 0 preserves reproducibility); limitation stated in §11 | no |
| A9 | cosmetic | `build_tref.py` final diagnostic print crashed (`nstep` undefined, LT index 48 out of range) after the lookup was already saved | fixed | no |
| A10 | cosmetic | `detection_depth.py` printed "~0 GW" for the 10×10 km case and "1e+14" for 1.45e14 m² (format widths) | fixed | no |

Verified correct: Planck band integration and mixed-pixel inversion (independent recomputation agrees to all digits); G = 4σT³ values; the invisible-flux and buried-water tables; the 210 K conduction sink vs 190/230 K radiative stiffness (different quantities, defensible); every §8b and §9 number against fresh seeded runs (byte-identical outputs); the eccentricity-weighted annual insolation; area weighting; GRS heat-production formula; MOLA lapse sign; the thermal model's orbit, declination, energy-balance signs, skin depth, I→k conversion, and explicit-scheme stability (κΔt/Δz² = 0.25).

## B. Statistics, bounds and counts

| id | severity | finding | fix | conclusion changed? |
|---|---|---|---|---|
| B1 | **major** | Of the 828 warmest-pixel outliers that were "fully pixel-analysed", **194 are uncalibrated byte-fill products** (SCALING_FACTOR 1, OFFSET 0, every DN identical; index MAX_BT = 255.0 or 250.0 is a fill value) and 4 more carry 32767 sentinel geometry. `btrlib` nulled only DN = 0, the screen only rejected MAX_BT outside 100–330 K, and `reprocess_outliers.py` filed them as "physical scenes, no cluster". Corrected split: 833 rows → 832 distinct → 828 with a BTR product → **630 genuine**: 36 calibration artifacts (median > 260 K), 594 physical scenes, 282 with a cluster (106 compact, 92 diffuse, 84 linear). 106 compact ≤ 256.5 K unchanged | `btrlib._reject_fill` (constant-DN / unit-scaling products → all-NaN); sentinel + fill rows dropped in `global_maxbt_screen.py`; new `global_figure_v2.py` → `global_quadrant_v2.png`; §7 rewritten | no |
| B2 | minor | The same fill products are counted among the 136,564 "analysed" global images: 908 stat rows have σ = 0 (constant DN), 50 more have < 5,000 valid px. Images with real pixel variance ≈ 135,600 (99.3%) | footnoted in §7 | no |
| B3 | minor | "829" double-counted one observation (I02021001 appears twice in the index); `global_figure.py` hard-coded "829 images flagged from 137,952" | v2 figure computes its counts from the data | no |
| B4 | **major (bookkeeping)** | Five image counts were used interchangeably: 137,952 (index rows with valid MAX_BT), 137,907 (worklist rows), 137,505 (hard-coded atlas title = done-list lines incl. 941 double-processed images), 136,567 (distinct), 136,564 (analysed). The headline value 136,564 appeared in no section body; §7 still said "we are streaming … growing to 100%" | §7 rewritten with the definitive accounting; atlas title → 136,564 (`global_anomaly_atlas_v2.png`); coverage title → distinct count and 99.8% (`coverage_map_v2.png`) | no |
| B5 | minor | Global "median NEΔT 1.3–1.4 K" came from small early test shards; the full sample (136,514 images) gives median σ = **1.58 K** (1.32 K at 0–30°N, 2–4 K at the poles). σ is a residual MAD after background removal and includes scene texture, so it is an upper bound on NEΔT | §7 corrected | no |
| B6 | minor | 1,340 duplicate worklist rows (two calibration versions per observation) were rasterized in `coverage_map.py` and processed twice by `btr_stream.py` (941 obsids): ≥5-look fraction 94.9 → 94.7%, ≥10 67.1 → 66.4%, mean looks 11.1 → 10.9; median 10 and 99.8% unchanged; ~7.84 M distinct clusters of the 7.88 M rows | noted in §7/§8 | no |
| B7 | minor | Injection-recovery: below 1.5 pixels of vent area only the MAX_BT screen can detect (single pixels are rejected by `MIN_CLUSTER_PIX = 2`), so the "screen ∨ cluster" description is literally true but the cluster branch never applies to sub-pixel vents. Among the 312 genuine no-cluster scenes the warmest pixel is 251.6 K, so nothing hot hides there | documented in §8 | no |
| B8 | cosmetic | Mars area 1.449e8 km² is 4πR² for a sphere; the oblate value is 1.4437e8 (0.4%, cancels in N₉₅). Marginalization over T_v by averaging A_eff over {273, 300, 350} K is the correct expectation for independently drawn vent temperatures; a uniform-in-T interpolation gives N < 1.94 instead of 1.98 | documented | no |
| B9 | cosmetic | "~1,366 non-recurring" → 1,363 (six detections form the three recurring cells); "tails to ~25 K" → 99th percentile 19 K, max 38 K; "37 uniform artifacts (no localized source)" — 8 of the 36 carry a cluster | text softened | no |

Verified correct: Poisson zero-count limits (Jeffreys Γ(½,1) 95% = 1.921; uniform 2.996 / 4.605); A_eff redundancy logic; every §8 number from a fresh run (N < 1.92 / 3.00 / 4.61; marginalized 1.98; 1.36e-8 km⁻²; 4.42 at 794 m²; 86.5 at 79 m²); coverage rerun byte-identical (99.8%, median 10, 2.77e6 cell-looks); injection ceiling 238.98 K and A50 ≈ 3,000–8,000 m²; §4 cascade 33,365 → 6,248 → 5,344 → 2,606 → 1,369 → 3 exact (positions aside, see C2); morphology 53/29/19%; regional σ median 1.29 K; 7,880,590 cluster rows; 21,385 persistent cells.

## C. Candidate / day–night chain

| id | severity | finding | fix | conclusion changed? |
|---|---|---|---|---|
| C1 | **critical** | **§10–§11 sampled the wrong pixels.** Each persistent candidate is a 0.05° cell whose catalogue position is the *mean* of its detections; `sample_night_bt.py` and `daytime_discriminate.py` took a 3×3 mean at the pixel nearest that mean (on a 4× downsampled grid). Detections within a cell are separated by a median 0.8 km, so the samplers landed a median 3 px (up to 12 px) from the warm pixels. Evidence: re-sampled night excess median −0.1 K vs catalogue 10 K (r = −0.06); on 38 cells with ≥ 15 K catalogue excess, the recorded pixel really is 10.4 K warm but the sampled position reads 2.1 K. Hence "84.8% cool by day", "r = +0.01", "ΔT_geo median 0.01 K / 95th 3 K", and "every candidate individually fit and found thermophysical" described background, not the anomalies | **Redone** at the detection pixels with window sampling on the day side, control-field bias calibration and per-image Ls/LT in the thermal model (`resample_candidates_v2.py`, `fit_residual_v2.py`, outputs in `data/analysis/v2/`). Results: see §C1-results below | **§10–§11 statements replaced**; overall null result unchanged (it never rested on §10–§11) |
| C2 | **major** | Regional SNU-cube georeferencing (§4) interpolated *raw-image* corner coordinates onto the *map-projected* array; raw line 1 is the south end, so latitude was reversed along the strip: position errors 22 km (289-row cube) to 545 km (10,692 rows). Global BTR work uses raw products and is unaffected | `anomaly_search.py` now uses the cube's own sinusoidal projection; regional search rerun → `candidates_v2.csv`, `focused_candidates_v2.csv`; see §C2-results | §4 recurring-cell positions replaced; regional null unchanged |
| C3 | **major** | `combine_grids.py` looked up TES albedo at longitude L−180 for grid longitude L (grid 0–360, tif −180..180): Tharsis got Syrtis's albedo and vice versa. With the fix the §12 night-excess floor and temperate exceedances change: see §C3-results | fixed in place; `global_daynight_mosaic_v2.png`, `daynight_grids_v2.npz` | no (screen result unchanged) |
| C4 | major (description) | `fit_residual.py` seasonal binning was dead code: the candidate file has no Ls column, so every candidate was modelled at Ls = 250; day and night images were not season-matched | v2 fit uses each image's own Ls and local time | subsumed by C1 |
| C5 | major (limitation) | Morphology gate: `MIN_CLUSTER_PIX = 2` plus the elongation test send every ≤ 2-px-wide cluster to "linear/artifact"; the smallest "compact" cluster is 3×3 px, so the focused catalogue is blind below ~300 m. The §8 bound rests on the MAX_BT screen and is unaffected | stated in §4/§8 | no |
| C6 | major (wording) | "Sunlit ceiling ~325–330 K" is the sub-solar *noon* perihelion equilibrium (326–335 K for A 0.15–0.05); THEMIS day passes are at ~15 h local time, where the project's own model tops out near 300–312 K, so the 306.8 K margin was not apples-to-apples | text corrected; the model fit (v2) is the actual test | no |
| C7 | minor | 4,214 daytime worklist rows have lt = ls = inc = 0 (bad metadata); `grid_accumulate.py` destriped them as midnight, biasing ~2% of day images (cells including one by a few K). Re-accumulating the day grid takes ~11 h and was not redone | excluded in the v2 candidate sampling; caveat in §12 | no |
| C8 | minor | §12 coverage "95.8 / 97.1%" is cnt > 0; at the MINCNT = 20 actually used it is 95.7 / 97.1% | text | no |
| C9 | minor | "top-100 nearest-neighbour 3.8° ≈ 225 km" used flat degrees; great-circle 3.52° = 208 km, and the top-100 are slightly *more* clustered than random (4.9°) | superseded by C1 v2 analysis | no |
| C10 | minor | "temperate blob aspect 4–6": the eight largest are 2.7–11.9 | text ("~3–12") | no |
| C11 | cosmetic | `sample_context.py` row index ≤ 1 px off; unused helpers in `combine_grids.py`/`anomaly_search.py`; day local time measured then dropped in `daytime_discriminate.py` | documented | no |
| C12 | minor | Literature review: the objectives "quote" splices the abstract and Section 1; "~2–5 a.m. orbit" is not in Christensen et al. 2004 (archive night local times are ~3–6 h); Eridania heat flow was cited to Tarnas et al. 2021 but is Ojha et al. 2021 and is a *Noachian* value | corrected in `literature_review.md` and REPORT §1 | no |
| C13 | minor | "Cerberus Fossae, where Mars's marsquakes originate" overstated (dominant located source; S1222a and dichotomy-boundary events elsewhere) | text | no |

Verified correct: the anomaly detector core (background, MAD σ, thresholds, morphology PCA, edge handling; synthetic tests recover a 3×3 +12 K spot and produce zero false clusters on a gradient strip); corner-based georeferencing of *raw* BTR products (r = +0.99 against a projected cube); persistent-cell definition (21,385 rows, cell ids consistent); TES albedo/TI orientation in `sample_context.py`; §10/§11 numbers reproduce from the CSVs exactly as written (the CSVs themselves are what was wrong).

### C1-results (day/night per-candidate redo)

`resample_candidates_v2.py` rebuilt the 21,385 cells from the shards (identical to the
catalogue), re-read each cell's strongest detecting night image from the cache and
recomputed the background at the detected pixel; the recomputed night excess matches the
catalogue to 0.1 K for 98.8% of cells (median 9.9 K). The detection coordinates were
georeferenced into 9,385 daytime images (exact footprint test; 4,214 bad-metadata images
excluded) and sampled in a 17×17-px window (cross-image georeferencing ≈ 1 km), with the
window-maximum bias calibrated against 20 random control sites per image. 21,223 cells have
night + day; 21,151 have a converged thermal fit. Runtime 12.5 min on 40 processes.

Two further flaws surfaced: `daytime_match.py` used a corner bounding-box test, so about
two thirds of the original day "samples" were edge pixels of images that do not contain the
point; and the original ΔT_geo differential implicitly assumed the warm pixel shares its
surroundings' albedo, so any dark rock registers as "endogenic".

| quantity | as published (§10–§11) | at the detected pixels |
|---|---|---|
| cool/neutral by day (< 2 K) | 84.8% (median day residual 0.13 K) | window-median: 88.3% (+0.05 K); centre pixel: 75.4% (+0.15 K); window-max (charitable, calibrated): 28.9% (+4.4 K, matched by −3.9 K window-min = terrain texture) |
| r(night excess, day excess) | +0.01 | −0.04 (window-max), −0.01 (window-median) |
| strongest 205 night features (> 20 K): median day excess | — | 3.3 K (window-max), 0.0 K (window-median): no proportional warm-both branch |
| hottest day pixel | 306.8 K (claimed below a 325–330 K noon ceiling) | 338 K, in a scene-wide 336 K calibration-outlier image (+1.1 K over its background); 4 cells > 310 K |
| ΔT_geo, original same-albedo definition: 50th / 95th / > 5 K | 0.01 K / 3.0 K / 1.1% | **8.4 K / 14.9 K / 93%** — driven by sunlit-slope day pixels (r = 0.86 with day window-max, 0.01 with night excess) |
| ΔT_geo, albedo-free ceiling: 95th, F_down = 0 / 15 W m⁻² | — | 5.9 K / 1.0 K (window-max day value); 1.0 K (centre pixel); 0.0 K (window median) |
| cells inside the flat-surface manifold with zero internal heat | "100%" | 94–98% |

Consequence: §10–§11 rewritten. The supportable statements are that night and day excesses
are uncorrelated with no co-located proportional day excess (the population is dark
high-inertia rock), and that a 1-D flat-surface fit bounds the endogenic warming of ≥95% of
candidates below ~1–6 K depending on model assumptions, with the tail dominated by sunlit
slopes. The "< 3 K for 95%, median 0.01 K" claim, the "100% thermophysical" headline row and
the abstract sentence are withdrawn and replaced. The §8 hot-spring bound never depended on
§10–§11 and is unchanged.

### C2-results (regional georeferencing rerun)

`anomaly_search.py` now georeferences from each cube's own sinusoidal projection (ISIS-3 v2
cubes via the rasterio affine + CRS; PDS-3 v1 cubes via the label projection, undoing the
per-line `RECTIFY_LEFTEDGE` shift and anchoring to the THMIDX corners because the v1 label
geometry is itself 1.4–8 km off, as shown by cross-correlation against overlapping v2 cubes).
Validation: north-up in 230/230 cubes (old path 0/229); footprint centres within 2 km of the
index centre for 219/230; overlapping-cube registration ≤ 4 px. Full evidence in
`data/analysis/regional_v2_README.md`.

| quantity | as published | corrected |
|---|---|---|
| cluster positions | flipped; median error 332 km, max 1,086 km | product frame, ~0.3–0.5 km |
| non-positional columns (33,365 rows), σ, morphology, cascade to 1,369 | — | byte-identical |
| recurring 0.05° cells across ≥2 overpasses | 3 (all "162°E") | **69** (Cerberus 43, Margaritifer 15, Elysium 10, Olympus 1; peak ≤ 14.6 K) |
| recurring cells within 0.1° of a global persistent cell | not checkable | 31/69 (45%) |
| focused detections within 0.1° of a global cell | 7.6% (chance 10–15%) | 29.1% |
| strongest one-offs | positions unusable | Ascraeus 11.60°N 254.67°E (+17.6 K; +26.3 K/150 px nearby, no global counterpart); Cerberus 3.69°N 159.96°E (+15.6 K; global cells 2–4 km away) |

The original "3 recurring locations, all at 162°E" statement and the §4 day/night
co-location check were artefacts of the flip and are withdrawn. Regional conclusion
(thermophysical, no vent) unchanged. Side effect: `focused_candidates.csv` was regenerated
with identical content during verification.

### C3-results (mosaic albedo fix)

Rerun of `combine_grids.py` with the longitude fix (`global_daynight_mosaic_v2.png`,
`daynight_grids_v2.npz`):

| quantity | as published | corrected |
|---|---|---|
| night-excess 95th / 99th percentile | 16.9 / 27.3 K | **13.7 / 22.5 K** |
| temperate cells > 15 K / > 20 K | 2.77% / 0.86% | **1.64% / 0.34%** |
| coverage night / day (≥20 samples) | 95.8 / 97.1% | 95.7 / 97.1% |
| largest blobs overall | polar | polar (lat 74.8, 73.9, −82.7, 76.2) |
| ten largest temperate blobs, aspect | "4–6, never compact" | 2.1–46.7 (four are < 3) |
| ten largest temperate blobs, mean day anomaly | not checked | **−6 to −14 K (all cool by day)**; night anomaly +13 to +32 K |
| temperate > 15 K cells with day anomaly < 0 | not checked | 94% |

The original "elongated, never compact" argument was wrong for four of the ten largest
blobs; the day/night sense replaces it and is decisive (thermal inertia lifts night and
lowers day; a geothermal offset lifts both). Conclusion unchanged.

## D. Figures, captions, PDF, notes

| id | severity | finding | fix | conclusion changed? |
|---|---|---|---|---|
| D1 | **major** | `global_quadrant.png` drew the hot-spring boundary at 250 K and wrote "EMPTY" while six compact points at 250.5–256.5 K sat inside; the text acknowledged them, the figure and caption did not | `global_quadrant_v2.png` (boundary at the 260 K bedrock ceiling, the six points labelled, fill products removed) | no |
| D2 | **major** | `daynight_discrimination.png` caption said "none in the endogenic regime" while the shaded region held 15% of candidates; and the whole figure is superseded by C1 | caption rewritten; v2 figure from the redo | see C1 |
| D3 | **major** | `global_daynight_mosaic.png` caption/§12 said "no coherent broad anomaly" while the middle panel showed province-scale ±5–10 K structure (manifold-fit residual over dust provinces); three thresholds (15/17/20 K) never reconciled; and the albedo layer was misregistered (C3) | v2 figure and caption qualified "above the screen floor" | no |
| D4 | major | Atlas title hard-coded "137,505 images"; coverage title said "100%" (rounding of 99.82) | `_v2` figures | no |
| D5 | major | The two strongest one-off regional candidates (`zoom_I01474006_7458_1102.png`, `zoom_I54939008_6465_1013.png`; +15.6 K and +17.6 K, and a +26.3 K/150-px cluster in the same Ascraeus image) were never discussed; §4 dismissed all 1,363 non-recurring detections in one sentence | discussed in §4 with corrected positions from C2 | no (absolute T 155–175 K, far below any hot-spring signature) |
| D6 | major | `depth_to_water_hist.png`: components each normalized to 1 under a "fraction of Mars surface" label, and a 2.5% spike at 0.1 m (depth-0 draws clipped to the log floor) under the title "note the empty shallow tail" | `depth_to_water_hist_v2.png` (mixture-weighted components; depth-0 spike drawn and labelled as ambient-temperature eutectic brine) | no |
| D7 | minor | Headline table said "84.8% of the 21,385" (denominator was 21,308); "4,043 day images cover all 21,385" (21,308 sampled); "≈60 m" vs table 66 m; 1.9 vs 1.92; "95% ≥ 5 looks" vs 94.85% | harmonized | no |
| D8 | minor | The `daynight_I05456012_I44697047.png` pair overlaps only in a ~0.3° corner; caption now says so | caption | no |
| D9 | cosmetic | Clipped title in `img_I05456012.png`; unexplained "?" region bar in `fig_by_region.png`; colorbars capped (persistence 8, ΔT_geo 8); §12 figures land after the Figures/Reproduce subsections; PDF Title metadata "themis_report.html"; `notes/reconstruction_notes.md` said 46 scripts (41) | notes corrected; others documented | no |
| D10 | minor | Unsupported literals checked: Yellowstone ~2 W m⁻² (OK, Morgan et al. 1977), Mars background 20–25 mW m⁻² (OK), Athabasca lava < 20 Ma (OK, conservative), TES albedo 7 km (OK, 8 ppd = 7.4 km) | — | no |

## E. Reproducibility

Re-running `detection_limit`, `detection_depth`, `coverage_map`, `global_maxbt_screen`,
`bayes_bound`, `vet_candidates`, `depth_to_water`, `depth_atlas` and `build_tref` from the
on-disk inputs produced byte-identical outputs (md5) before any fix was applied. All fixes
are in new `_v2` scripts/outputs or in-place code changes documented at the top of each
changed file.