Methods and caveats

For the reader who wants the chain of reasoning, the numbers, and the places where the numbers are soft. The full report with all sections is in the PDF; the section numbers below follow it.

1. Objective and literature
2. Data
3. Sub-pixel detectability
4. Regional search
5. What stays invisible
6. Global warmest-pixel screen
7. Bayesian bound
8. Depth to liquid water
9. Day–night at the detected pixels
10. Thermal-model residual
11. Global mosaic
12. Caveats

1. Objective and literature

Christensen et al. (2004, Space Sci. Rev. 110, 85) list the search for "pre-dawn thermal anomalies associated with active sub-surface hydrothermal systems" as one of THEMIS's science objectives, to be pursued with global nighttime multispectral temperature maps at 100 m/pixel and a noise-equivalent ΔT of 1 K, initially focused on young volcanic sites. The published outcome was null (Christensen et al. 2003 and later citations). Diffuse crustal heat flow is ~20 mW/m² today and cannot raise the surface temperature by 1 K; even the radiogenic Eridania region only reached ~65 mW/m² in the Noachian (Ojha et al. 2021, Nat. Commun. 12, 1754). Only focused heat is detectable. The literature review has the citations.

2. Data

• THEMIS IR cumulative index (1,127,080 observations) and product indices from the ASU PDS archive.
• Regional: 230 nighttime band-9 (12.57 µm) map-projected radiance cubes (~16 GB) over Cerberus Fossae, Elysium, Olympus, Arsia, Pavonis and Ascraeus Montes plus a Margaritifer Terra control. Radiance is inverted to brightness temperature through the Planck function; the reader handles 3-band and 10-band cubes, ISIS-3 and PDS-3 formats, and the per-line rectification offsets of the older products.
• Global: every nighttime band-9 brightness-temperature record (BTR) product, 136,567 distinct images (137,907 index rows), plus 198,858 daytime BTR products, all cached locally (372 GB). Nighttime is incidence angle > 90°; daytime < 85°.
• Ancillary: Odyssey GRS thorium and potassium (5°), MOLA MEGDR topography (4 ppd), MGS TES albedo (8 ppd) and nighttime thermal inertia (20 ppd).

Two data-quality traps matter. Some BTR products are uncalibrated byte fills (scaling factor 1, offset 0, every pixel the same value) whose index "maximum temperature" of 255.0 or 250.0 is a fill value; and some index rows carry 32767 sentinels in place of coordinates. Both are now rejected by the readers and the screen.

3. Sub-pixel detectability

The band-9 radiance of a pixel containing a fraction f of vent at temperature T_v over regolith at T_bg is L = (1−f)B(T_bg) + fB(T_v), with B the band-integrated Planck function. Requiring the retrieved brightness temperature to rise by ΔT gives the minimum vent area. For T_bg = 190 K and a 100 m pixel:

The numbers scale with the background: at 170 K the boiling-water threshold is 10 m², at 210 K it is 25 m². The 1 K specification is defined near 245 K; at 190 K the radiometric floor is somewhat worse, which the empirical 1.3 K and the 5 K column cover.

4. Regional search

Detector (anomaly_search.py): remove a ~3 km block-median background (robust to point sources), measure the residual noise as 1.4826 × MAD, flag warm residuals above both 5 K and 5σ, reject clusters touching the data edge, larger than 2,000 px or smaller than 2 px, and classify morphology by principal-component elongation and bounding-box fill into compact, linear/artifact and diffuse. Over 230 regional images: median residual noise 1.29 K (an upper bound on NEΔT since it includes scene texture), 33,365 warm clusters, 53% linear, 29% diffuse, 19% compact, excess peaking at 6–10 K with a 99th percentile of 19 K.

Vetting cascade: 33,365 → 6,248 compact → 5,344 in clean images (σ ≤ 1.5 K) → 2,606 small (≤ 12 px) → 1,369 strong (≥ 8 K) → 69 recurring in the same 0.05° cell across two overpasses. The 69 (Cerberus 43, Margaritifer 15, Elysium 10, Olympus 1) have peak excess ≤ 14.6 K and absolute temperatures of 170–200 K; 45% coincide within 0.1° with cells of the independent global catalogue. They are persistent rocky surfaces. The strongest one-off detections (Ascraeus Mons +17.6 K and +26.3 K over 150 px; Cerberus +15.6 K) are compact bright spots at ~155–175 K absolute, consistent with outcrop or pit skylights, not water.

Note that the original version of this step georeferenced the map-projected cubes with raw-image corner coordinates, reversing latitude along each strip (errors of hundreds of km); the cluster statistics were unaffected but the "3 recurring cells at 162°E" it reported were an artefact. Positions were recomputed from each cube's own projection on 2026-09-07 (validated to ~0.5 km against overlapping cubes).

6. What stays invisible

A steady heat flux q reaching the surface re-radiates as a brightness-temperature offset ΔT = q/G with G = 4σT³, the radiative stiffness: 1.6 W m⁻² K⁻¹ at 190 K, 2.1 at 210 K, 2.8 at 230 K. Hence:

Integrated, up to ~20 GW over a 100 × 100 km region or ~3 × 10⁵ GW globally could be present as diffuse heat and be unseen. A buried water body at T_w under overburden of conductivity k leaks q = k(T_w − T̄_surf)/d and is invisible when deeper than d = k(T_w − T̄)/(G ΔT_det): for 273 K water, 0.9 m under air-fall dust (k = 0.03), 3 m under loose regolith, 15 m under duricrust, 66 m under ice-cemented rock (k = 2.2); for 373 K water, 2.3 to 171 m. A 10 m slab losing the ceiling flux stays warm at least ~40 years (100 m: ~400 years), two to three times longer with the latent heat of freezing, so persistently warm hidden water must be replenished by circulation, which is what makes vents.

7. Global warmest-pixel screen

The index lists each image's warmest pixel. Per 2° latitude × 30° L_s bin, the robust median and MAD of the warmest pixel define the natural bedrock ceiling; images with z > 5 are flagged: 833 rows (832 observations) of 137,952. Of the 828 with a product, 194 are byte fills and 4 have sentinel geometry. The 630 genuine images: 36 calibration artifacts (whole-image median > 260 K), 594 physical night scenes, 282 with a localized cluster (106 compact, 92 diffuse, 84 linear). The 106 compact features are all ≤ 256.5 K: high-latitude defrosted spots in ~170 K frost scenes, and the six hottest are perihelion equatorial bedrock with 10–15 K excess over ~1 km outcrops. Nothing compact exceeds the ~260 K bedrock ceiling.

Independently, all 136,564 analysable images were streamed through the detector: 7.88 million warm clusters (about 7.84 million distinct), median per-image residual noise 1.58 K globally and 1.32 K in the equatorial band. Persistent focused cells (≥ 2 overpasses, 0.05° grid): 21,385, all in the equatorial rocky belt.

8. Bayesian bound

Coverage. Rasterizing all footprints on a 0.5° grid: 99.8% of Mars imaged at night at least once, 94.9% at least five times, 67% at least ten; median 10 looks per cell, mean 11; 2.77 million cell-looks. (Removing the 1,340 duplicate index rows changes the ≥5 and ≥10 fractions by under one point.)

Efficiency. Synthetic vents (band-integrated Planck mixing; T_v ∈ {273, 300, 350, 400, 500, 700} K, area 2–20,000 m²) were injected into 310 real BTR images stratified by latitude and recovered through the same screen and cluster logic. A sub-pixel vent is caught only if it lifts its pixel over the bedrock ceiling (~239 K), because single pixels are rejected by the cluster stage; so the 50%-recovery area per look is 7,900 / 5,000 / 3,200 m² at 273 / 300 / 350 K. Ten looks bring completeness to ~1 above those areas.

Posterior. With zero detections, the 95% upper limit on the Poisson mean is 1.92 (Jeffreys prior), 3.00 (flat), 4.61 (flat, 99%). The effective searched area A_eff = Σ area × [1 − (1 − η)^N] over cells; N₉₅ = μ₉₅ × A_Mars / A_eff. For a vent of ~8,000 m² at ≥ 300 K (completeness 0.998): N < 1.92 (Jeffreys), 3.00 (flat), 4.61 (99%). Marginalized over T_v uniform in 273–373 K: N < 1.98, areal density < 1.36 × 10⁻⁸ km⁻², one vent per 7.3 × 10⁷ km². At 794 m² and 350 K: N < 4.4; at 79 m²: N < 86.

8b–9. Depth to liquid water

Below the seasonal skin T(z) = T̄_surf(lat) + (q/k) z, so liquid exists below d = (T_melt − T̄_surf) k/q. Six million area-weighted draws with T̄_surf from annual-mean insolation (220 K equator, 155 K pole, ±6 K), q lognormal with median 21 mW/m² (σ_ln 0.30, clipped 8–60), k ~ N(2.3, 0.5) W m⁻¹ K⁻¹ clipped to 1.2–3.6, and T_melt a mixture of 45% pure water (273 K), 35% perchlorate brine (252 K), 15% chloride brine (230 K), 5% extreme eutectic (210 K):

The 2-D atlas replaces the heat-flow prior with a GRS thorium-and-potassium heat-production map (mean 20.0 mW/m², range 16–23.5, high over Acidalia and Arabia, low over Hellas and Elysium) and adds a −1.5 K/km lapse on MOLA elevation: median 6.3 km, 4–6 km at the equator, ~14 km at the poles; P(< 1 m) = 2.9 × 10⁻², P(< 1 km) = 5.5 × 10⁻², P(< 5 km) = 0.39. The equatorial annual-mean temperature in the model is ~7 K warm (the equilibrium temperature of the mean flux exceeds the mean temperature); correcting it deepens the medians by ~10% and halves the brine tail, so the published numbers are conservative.

10. Day–night at the detected pixels

For each of the 21,385 persistent cells the strongest detection is re-read from the detecting night image at its recorded pixel (the recomputed excess matches the catalogue to 0.1 K in 98.8% of cells; median 9.9 K). The detection's coordinates are georeferenced into the covering daytime images (9,385 images, exact footprint test; 4,214 images with missing local-time metadata excluded). Because cross-image georeferencing is good to ~1 km, the day side is sampled in a 17 × 17 pixel window: the window maximum of the background-removed temperature (most charitable to a warm-both source) calibrated against 20 random control windows in the same image, plus the window median and centre pixel. 21,223 cells have both.

Result: candidate sites are rugged (calibrated window maximum +4.4 K, window minimum −3.9 K, symmetric texture); window median +0.05 K with 88% within 2 K; correlation between night and day excess −0.04 (window max) and −0.01 (median); the 205 strongest night features (> 20 K) have a median day maximum of 3.3 K and a window median of 0.0 K. There is no warm-both branch. The hottest day pixel in any window, 338 K, lies in one of 31 daytime images whose whole scene reads above 330 K, a calibration outlier; the candidate's excess there is +1.1 K.

11. Thermal-model residual

A vectorized 1-D diurnal model (40 layers, explicit diffusion, Newton surface balance, Kepler orbit, optional constant sky flux) is run on an inertia × albedo grid (I 30–1500, A 0.05–0.40) at each image's own season and local time. The diurnal amplitude fixes the inertia; a geothermal source adds a DC offset. Required warming ΔT_geo = observed level − maximum thermophysical level at that amplitude.

The same-albedo differential is not a test of internal heat when applied to a real anomaly: it assumes the warm pixel shares its surroundings' albedo, so any dark rock in lighter dust registers as endogenic (a 0.1 albedo drop alone raises the diurnal mean by ~8 K). It correlates with the daytime window maximum (r = 0.86), not with the night excess (r = 0.01). Allowing the candidate its own inertia and albedo, 94–98% of cells need no internal heat, and the tail is southern mid-latitude, low-sun scenes with a sunlit slope pixel in the window. The supportable per-candidate ceiling is ~1–6 K depending on sky flux and day statistic. The earlier claim of "< 3 K for 95%, median 0.01 K" was computed from samples that missed the detected pixels and is withdrawn.

12. Global mosaic

All 335,425 day and night images were binned to 0.05° after subtracting a reference-surface temperature T_ref(lat, L_s, local time) from the same model (I = 200, A = 0.15), which removes orbital striping. Coverage at ≥ 20 samples per cell: night 95.7%, day 97.1%. Thermal inertia and albedo trace a manifold night_anom = f(day_anom, TES albedo); a geothermal offset lifts a cell above it. The 95th percentile of the night excess is 13.7 K (the floor, from sub-grid heterogeneity and 7 km albedo); 1.64% of temperate cells exceed 15 K and 0.34% exceed 20 K. The largest connected exceedances are polar (74–83°); the ten largest temperate ones, 2,000–7,000 cells each, are all cooler than their surroundings by day (day anomaly −6 to −14 K against night anomalies of +13 to +32 K), and 94% of all temperate above-threshold cells have negative day anomalies: thermal inertia, not internal heat. (The original version of this screen misregistered the albedo layer by 180° of longitude; the corrected floor is 14 K rather than 17 K.)

Caveats

• Small, cool seeps are unconstrained. Below the bedrock confusion floor (~100 m² of near-freezing water per pixel) the bound relaxes to tens of vents. Only the conduction argument of section 6 speaks to that regime.
• Cluster morphology is blind below 3 × 3 pixels. The compact-cluster catalogue cannot see sources narrower than ~300 m; the sub-pixel regime is covered by the warmest-pixel screen, whose efficiency was measured directly.
• The thermal model is 1-D, flat and has no atmosphere beyond an optional constant sky flux. Nights at low inertia come out 10–40 K too cold. Differencing against in-image backgrounds cancels the absolute bias (a 15 W/m² sky flux changes the section 11 percentiles by < 0.1 K), but slopes are not represented, which sets the few-kelvin per-candidate ceiling.
• Georeferencing between day and night images is ~1 km, hence the window sampling and its texture-driven ±4 K spread.
• The mosaic's 14 K floor makes section 12 a broad-anomaly check only; about 2% of daytime images with missing local-time metadata were destriped as midnight and bias the cells they touch by a few kelvin.
• The depth atlas uses a uniform 50 km heat-producing column; crustal-thickness modulation would widen the heat-flow range somewhat but cannot move the latitude-set kilometre-scale floor. The K heat-production constant in the code is 16% low, absorbed by the 20 mW/m² calibration.
• Byte-fill products (about 0.7% of images, mostly polar) are counted among the images "analysed"; the fraction of the planet with real pixel data is ~99.3% rather than 99.998%.
• Prior sensitivity of the bound is modest (1.9 to 3.0 between Jeffreys and flat priors) and is reported rather than hidden.

Everything above was re-derived from the on-disk data on 2026-09-07; the reruns of the deterministic scripts reproduced the original outputs byte for byte, and the three analysis bugs found were fixed and rerun. See the error check.